In this paper we design a fast bisection algorithm for computing a prescribed eigenvalue of a symmetric quasiseparable matrix represented by its Cholesky-like factorization A0=L0Ω0L0T, where L0 is a lower triangular quasiseparable matrix and Ω0 is a signature matrix. The workhorse of the algorithm is an efficient linear-time scheme for updating the factorization under shifting. Numerical experiments illustrate the effectiveness and the robustness of the proposed method.

This paper presents a Jacobi-type iteration for computing a given specified eigenpair of a symmetric matrix. For a certain class of diagonally dominant matrices, the procedure is shown to converge at a linear rate depending on how the matrix is significantly dominated. The cost per iteration is generally quadratic. Therefore, the proposed procedure can compute an approximation of the desired eigenpair in quadratic time.

In this paper we develop fast numerical algorithms for solving shifted linear systems with semidefinite quasiseparable matrices. A combination of Givens and hyperbolic plane rotations is used to update the Cholesky-type factorization of the input quasiseparable matrix by determining a factorization of its shifted version of the form LDL^T, where L is lower triangular and D is a signature matrix. If the shifted matrix is also definite then the Cholesky factorization of the shifted matrix is computed in a stable way by using orthogonal transformations. Since quasiseparability is maintained under diagonal shifting, a fast variant of the updating procedure using computations with generators is also devised. Numerical experiments show the effectiveness and robustness of the proposed algorithm.

This paper deals with efficient numerical methods for computing the action of the matrix generating function of Bernoulli polynomials, say q(τ ,A) , on a vector when A is a large and sparse matrix. This problem occurs when solving some non-local boundary value problems. Methods based on the Fourier expansion of q(τ ,w) have already been addressed in the scientific literature. The contribution of this paper is twofold. First, we place these methods in the classical framework of Krylov-Lanczos (polynomial-rational) techniques for accelerating Fourier series. This allows us to apply the convergence results developed in this context to our function. Second, we design a new acceleration scheme. Some numerical results are presented to show the effectiveness of the proposed algorithms.

This paper aims to develop efficient numerical methods for computing the inverse of matrix φ-functions, ψ_ℓ(A) := (φ_ℓ(A))^-1, for ℓ =1,2,…, when A is a large and sparse matrix with eigenvalues in the open left half-plane. While φ-functions play a crucial role in the analysis and implementation of exponential integrators, their inverses arise in solving certain direct and inverse differential problems with non-local boundary conditions. We propose an adaptation of the standard scaling-and-squaring technique for computing ψ_ℓ(A), based on the Newton-Schulz iteration for matrix inversion. The convergence of this method is analyzed both theoretically and numerically. In addition, we derive and analyze Padé approximants for approximating ψ_1(A/2^s), where s is a suitably chosen integer, necessary at the root of the squaring process. Numerical experiments demonstrate the effectiveness of the proposed approach.

We provide a new approach to obtain solutions of linear differential problems set in a Banach space and equipped with nonlocal boundary conditions. From this approach we derive a family of numerical schemes for the approximation of the solutions. We show by numerical tests that these schemes are numerically robust and computationally efficient.

The paper is concerned with efficient numerical methods for solving a linear system ϕ(A)x=b, where ϕ(z) is a ϕ-function and A∈RN×N. In particular in this work we are interested in the computation of ϕ(A)−1b for the case where ϕ(z)=ϕ1(z)=ez−1z and ϕ(z)=ϕ2(z)=ez−1−zz2. Under suitable conditions on the spectrum of A we design fast algorithms for computing both ϕℓ(A)−1 and ϕℓ(A)−1b based on Newton's iteration and Krylov-type methods, respectively. Adaptations of these schemes for structured matrices are considered. In particular the cases of banded and more generally quasiseparable matrices are investigated. Numerical results are presented to show the effectiveness of our proposed algorithms.

We present some accelerated variants of fixed point iterations for computing the minimal non-negative solution of the unilateral matrix equation associated with an M/G/1-type Markov chain. These variants derive from certain staircase regular splittings of the block Hessenberg M-matrix associated with the Markov chain. By exploiting the staircase profile, we introduce a two-step fixed point iteration. The iteration can be further accelerated by computing a weighted average between the approximations obtained at two consecutive steps. The convergence of the basic two-step fixed point iteration and of its relaxed modification is proved. Our theoretical analysis, along with several numerical experiments, shows that the proposed variants generally outperform the classical iterations.

This paper presents the results of a preliminary experimental investigation of the performance of a stationary iterative method based on a block staircase splitting for solving singular systems of linear equations arising in Markov chain modelling. From the experiments presented, we can deduce that the method is well suited for solving block banded or more generally localized systems in a parallel computing environment. The parallel implementation has been benchmarked using several Markovian models.

We present some accelerated variants of fixed point iterations for computing the minimal nonnegative solution of the unilateral matrix equation associated with an M/G/1-type Markov chain. These schemes derive from certain staircase regular splittings of the block Hessenberg M-matrix associated with the Markov chain. By exploiting the staircase profile we introduce a two-step fixed point iteration. The iteration can be further accelerated by computing a weighted average between the approximations obtained in two consecutive steps. The convergence of the basic two-step fixed point iteration and of its relaxed modification is proved. Our theoretical analysis along with several numerical experiments show that the proposed variants generally outperform the classical iterations.

We present a class of fast subspace algorithms based on orthogonal iterations for structured matrices/pencils that can be expressed as small rank perturbations of unitary matrices. The representation of the matrix by means of a new data-sparse factorization—named LFR factorization—using orthogonal Hessenberg matrices is at the core of these algorithms. The factorization can be computed at the cost of O(n k^2) arithmetic operations, where n and k are the sizes of the matrix and the small rank perturbation, respectively. At the same cost from the LFR format we can easily obtain suitable QR and RQ factorizations where the orthogonal factor Q is a product of orthogonal Hessenberg matrices and the upper triangular factor R is again given into the LFR format. The orthogonal iteration reduces to a hopping game where Givens plane rotations are moved from one side to the other side of these two factors. The resulting new algorithms approximate an invariant subspace of size s associated with a set of s leading or trailing eigenvalues using only O ( nks ) operations per iteration. The number of iterations required to reach an invariant subspace depends linearly on the ratio |λ _s+1|/|λ _s| . Numerical experiments confirm the effectiveness of our adaptations.

The algebraic characterization of dual univariate interpolating subdivision schemes is investigated. Specifically, we provide a constructive approach for finding dual univariate interpolating subdivision schemes based on the solutions of certain associated polynomial equations. The proposed approach also makes it possible to identify conditions for the existence of the sought schemes.

Some variants of the (block) Gauss–Seidel iteration for the solution of linear systems with M -matrices in (block) Hessenberg form are discussed. Comparison results for the asymptotic convergence rate of some regular splittings are derived: in particular, we prove that for a lower-Hessenberg M-matrix ρ (P_GS)≥ρ (P_S)≥ρ (P_AGS) , where P_GS, P_S, P_AGS are the iteration matrices of the Gauss–Seidel, staircase, and anti-Gauss–Seidel method. This is a result that does not seem to follow from classical comparison results, as these splittings are not directly comparable. It is shown that the concept of stair partitioning provides a powerful tool for the design of new variants that are suited for parallel computation.

In this paper we introduce a family of rational approximations of the reciprocal of a $ϕ$-function involved in the explicit solutions of certain linear differential equations, as well as in integration schemes evolving on manifolds. The derivation and properties of this family of approximations applied to scalar and matrix arguments are presented. Moreover, we show that the matrix functions computed by these approximations exhibit decaying properties comparable to the best existing theoretical bounds. Numerical examples highlight the benefits of the proposed rational approximations w.r.t.~the classical Taylor polynomials and other rational functions.

We present a class of fast subspace tracking algorithms based on orthogonal iterations for structured matrices/pencils that can be represented as small rank perturbations of unitary matrices. The algorithms rely upon an updated data sparse factorization -- named LFR factorization -- using orthogonal Hessenberg matrices. These new subspace trackers reach a complexity of only $O(nk^2)$ operations per time update, where $n$ and $k$ are the size of the matrix and of the small rank perturbation, respectively.

It is shown that the problem of balancing a nonnegative matrix by positive diagonal matrices can be recast as a nonlinear eigenvalue problem with eigenvector nonlinearity. Based on this equivalent formulation some adaptations of the power method and Arnoldi process are proposed for computing the dominant eigenvector which defines the structure of the diagonal transformations. Numerical results illustrate that our novel methods accelerate significantly the convergence of the customary Sinkhorn–Knopp iteration for matrix balancing in the case of clustered dominant eigenvalues.

We present fast numerical methods for computing the Hessenberg reduction of a unitary plus low-rank matrix $A=G+U V^H$, where $G\in \mathbb C^{n\times n}$ is a unitary matrix represented in some compressed format using $O(nk)$ parameters and $U$ and $V$ are $n\times k$ matrices with $k< n$. At the core of these methods is a certain structured decomposition, referred to as a LFR decomposition, of $A$ as product of three possibly perturbed unitary $k$ Hessenberg matrices of size $n$. It is shown that in most interesting cases an initial LFR decomposition of $A$ can be computed very cheaply. Then we prove structural properties of LFR decompositions by giving conditions under which the LFR decomposition of $A$ implies its Hessenberg shape. Finally, we describe a bulge chasing scheme for converting the initial LFR decomposition of $A$ into the LFR decomposition of a Hessenberg matrix by means of unitary transformations. The reduction can be performed at the overall computational cost of $O(n^2 k)$ arithmetic operations using $O(nk)$ storage. The computed LFR decomposition of the Hessenberg reduction of $A$ can be processed by the fast QR algorithm presented in [8] in order to compute the eigenvalues of $A$ within the same costs.

Some fast algorithms for computing the eigenvalues of a (block) companion matrix have recently appeared in the literature. In this paper we generalize the approach to encompass unitary plus low rank matrices of the form $$A=U + XY^H$$ where U is a general unitary matrix. Three important cases for applications are U unitary diagonal, U unitary block Hessenberg and U unitary in block CMV form. Our extension exploits the properties of a larger matrix $$\hat{A}$$ obtained by a certain embedding of the Hessenberg reduction of A suitable to maintain its structural properties. We show that $$\hat{A}$$ can be factored as product of lower and upper unitary Hessenberg matrices possibly perturbed in the first k rows, and, moreover, such a data-sparse representation is well suited for the design of fast eigensolvers based on the QR iteration. The resulting eigenvalue algorithm is fast and backward stable.

In this paper, we focus on the solution of shifted quasiseparable systems and of more general parameter-dependent matrix equations with quasiseparable representations. We propose an efficient algorithm exploiting the invariance of the quasiseparable structure under diagonal shifting and inversion. This algorithm is applied to compute various functions of matrices. Numerical experiments show that this approach is fast and numerically robust.