Rounding Error Analysis in Solving M-Matrix Linear Systems of Block Hessenberg Form

The problem of solving large M-matrix linear systems with sparse coefficient matrix in block Hessenberg form is here addressed. In previous work of the authors a divide-and-conquer strategy was proposed and a backward error analysis of the resulting algorithm was presented showing its effectiveness for the solution of computational problems of queueing theory and Markov chains. In particular, it was shown that for block Hessenberg M-matrices the algorithm is weakly backward stable in the sense that the computed solution is the exact solution of a nearby linear system, where the norm of the perturbation is proportional to the condition number of the coefficient matrix. In this note a better error estimate is given by showing that for block Hessenberg M-matrices the algorithm is even backward stable.     

A note on the solution of not balanced banded Toeplitz systems

A direct algorithm is presented for the solution of linear systems having banded Toeplitz coefficient matrix with unbalanced bandwidths. It is derived from the cyclic reduction algorithm, it makes use of techniques based on the displacement rank and it relies on the Morrison–Sherman–Woodbury formula. The algorithm always equals and sometimes outperforms the already known direct ones in terms of asymptotic computational cost. The case where the coefficient matrix is a block banded block Toeplitz matrix in block Hessenberg form is analyzed as well. The algorithm is numerically stable if applied to M‐matrices that are point diagonally dominant by columns. Copyright © 2007 John Wiley & Sons, Ltd.     

Rounding Errors in Solving Block Hessenberg Systems

A rounding error analysis is presented for a divide-and-conquer algorithm to solve linear systems with block Hessenberg matrices. Conditions are derived under which the algorithm computes a stable solution. The algorithm is shown to be stable for block diagonally dominant matrices and for M-matrices.

     

A new implementation of the CMRH method for solving dense linear systems

The CMRH method [H. Sadok, Méthodes de projections pour les systèmes linéaires et non linéaires, Habilitation thesis, University of Lille1, Lille, France, 1994; H. Sadok, CMRH: A new method for solving nonsymmetric linear systems based on the Hessenberg reduction algorithm, Numer. Algorithms 20 (1999) 303–321] is an algorithm for solving nonsymmetric linear systems in which the Arnoldi component of GMRES is replaced by the Hessenberg process, which generates Krylov basis vectors which are orthogonal to standard unit basis vectors rather than mutually orthogonal. The iterate is formed from these vectors by solving a small least squares problem involving a Hessenberg matrix. Like GMRES, this method requires one matrix–vector product per iteration. However, it can be implemented to require half as much arithmetic work and less storage. Moreover, numerical experiments show that this method performs accurately and reduces the residual about as fast as GMRES. With this new implementation, we show that the CMRH method is the only method with long-term recurrence which requires not storing at the same time the entire Krylov vectors basis and the original matrix as in the GMRES algorithm. A comparison with Gaussian elimination is provided.     

Efficient eigenvalue computation for quasiseparable Hermitian matrices under low rank perturbations

In this paper we address the problem of efficiently computing all the eigenvalues of a large N×N Hermitian matrix modified by a possibly non Hermitian perturbation of low rank. Previously proposed fast adaptations of the QR algorithm are considerably simplified by performing a preliminary transformation of the matrix by similarity into an upper Hessenberg form. The transformed matrix can be specified by a small set of parameters which are easily updated during the QR process. The resulting structured QR iteration can be carried out in linear time using linear memory storage. Moreover, it is proved to be backward stable. Numerical experiments show that the novel algorithm outperforms available implementations of the Hessenberg QR algorithm already for small values of N.      

On the inertia of the block H‐matrices

The problem of determining matrix inertia is very important in many applications, for example, in stability analysis of dynamical systems. In the (point‐wise) H‐matrix case, it was proven that the diagonal entries solely reveal this information. This paper generalises these results to the block H‐matrix cases for 1, 2, and matrix norms. The usefulness of the block approach is illustrated on 3 relevant numerical examples, arising in engineering.     

A new algorithm for generalized Sylvester‐observer equation and its application to state and velocity estimations in vibrating systems

We propose a new algorithm for block‐wise solution of the generalized Sylvester‐observer equation XA?FXE = GC, where the matrices A, E, and C are given, the matrices X, F, and G need to be computed, and matrix E may be singular. The algorithm is based on an orthogonal decomposition of the triplet (A, E, C) into the observer‐Hessenberg‐triangular form. It is a natural generalization of the widely known observer‐Hessenberg algorithm for the Sylvester‐observer equation: XA?FX = GC, which arises in state estimation of a standard first‐order state‐space control system. An application of the proposed algorithm is made to state and velocity estimations of second‐order control systems modeling a wide variety of vibrating structures. For dense un‐structured data, the proposed algorithm is more efficient than the recently proposed SVD‐based algorithm of the authors. Copyright © 2010 John Wiley & Sons, Ltd.     

Updating the QR decomposition of block tridiagonal and block Hessenberg matrices

We present an efficient block-wise update scheme for the QR decomposition of block tridiagonal and block Hessenberg matrices. For example, such matrices come up in generalizations of the Krylov space solvers MinRes, SymmLQ, GMRes, and QMR to block methods for linear systems of equations with multiple right-hand sides. In the non-block case it is very efficient (and, in fact, standard) to use Givens rotations for these QR decompositions. Normally, the same approach is also used with column-wise updates in the block case. However, we show that, even for small block sizes, block-wise updates using (in general, complex) Householder reflections instead of Givens rotations are far more efficient in this case, in particular if the unitary transformations that incorporate the reflections determined by a whole block are computed explicitly. Naturally, the bigger the block size the bigger the savings. We discuss the somewhat complicated algorithmic details of this block-wise update, and present numerical experiments on accuracy and timing for the various options (Givens vs. Householder, block-wise vs. column-wise update, explicit vs. implicit computation of unitary transformations). Our treatment allows variable block sizes and can be adapted to block Hessenberg matrices that do not have the special structure encountered in the above mentioned block Krylov space solvers.     

Queueing models,block hessenberg matrices and the method of neuts

We introduce a numerical method to compute the stationary probability vector of queueing models whose infinitesimal generator is of block Hessenberg form. It is shown that the stationary probability vector is equal to the first column of the inverse of the coefficient matrix. Furthermore, it is shown that the first column of the inverse of an upper (or lower) Hessenberg matrix may be obtained in a relatively small number of operations. Together, these results allow us to define a powerful algorithm for solving certain queueing models. The efficiency of this algorithm is discussed and a comparison with the method of Neuts is undertaken. A relationship with the method of Gaussian elimination is established and used to develop some stability results.This work was supported in part by NSF Grant MCS-83-00438.     

On Using Hölder Norms in the Quasi-Minimal Residual Approach

Given a process to span a basis for the underlying Krylov subspaces, the quasi-minimal residual (QMR-)approach is often used to derive iterative methods for the solution of linear systems. The QMR-approach is only reasonable if the resulting methods are based on short recurrences. The key ingredient of the QMR-approach is the efficient solution of a sequence of least-squares problems by computing the QR decomposition of an upper Hessenberg matrix by means of Givens rotations. Since (Hölder) p-norms are not unitarily invariant, a generalization of the minimization problem from the Euclidean norm to general p-norms while still leading to methods based on short recurrences appeared infeasible. Here, it is shown that this generalization is possible if the upper Hessenberg matrix reduces to a lower bidiagonal matrix.This revised version was published online in October 2005 with corrections to the Cover Date.     

Preordering saddle‐point systems for sparse LDLT factorization without pivoting

This paper focuses on efficiently solving large sparse symmetric indefinite systems of linear equations in saddle‐point form using a fill‐reducing ordering technique with a direct solver. Row and column permutations partition the saddle‐point matrix into a block structure constituting a priori pivots of order 1 and 2. The partitioned matrix is compressed by treating each nonzero block as a single entry, and a fill‐reducing ordering is applied to the corresponding compressed graph. It is shown that, provided the saddle‐point matrix satisfies certain criteria, a block LDL^T factorization can be computed using the resulting pivot sequence without modification. Numerical results for a range of problems from practical applications using a modern sparse direct solver are presented to illustrate the effectiveness of the approach.     

Block s‐step Krylov iterative methods

Block (including s‐step) iterative methods for (non)symmetric linear systems have been studied and implemented in the past. In this article we present a (combined) block s‐step Krylov iterative method for nonsymmetric linear systems. We then consider the problem of applying any block iterative method to solve a linear system with one right‐hand side using many linearly independent initial residual vectors. We present a new algorithm which combines the many solutions obtained (by any block iterative method) into a single solution to the linear system. This approach of using block methods in order to increase the parallelism of Krylov methods is very useful in parallel systems. We implemented the new method on a parallel computer and we ran tests to validate the accuracy and the performance of the proposed methods. It is expected that the block s‐step methods performance will scale well on other parallel systems because of their efficient use of memory hierarchies and their reduction of the number of global communication operations over the standard methods. Copyright © 2009 John Wiley & Sons, Ltd.     

A new recursive algorithm for inverting Hessenberg matrices

In the current article, the authors present a new recursive symbolic computational Hessenberg matrix algorithm, for inverting general Hessenberg matrices.     

Non‐recursive solution of sparse block Hessenberg systems

A double‐phase algorithm, based on the block recursive LU decomposition, has been recently proposed to solve block Hessenberg systems with sparsity properties. In the first phase the information related to the factorization of A and required to solve the system, is computed and stored. The solution of the system is then computed in the second phase. In the present paper the algorithm is modified: the two phases are merged into a one‐phase version having the same computational cost and allowing a saving of storage. Moreover, the corresponding non‐recursive version of the new algorithm is presented, which is especially suitable to solve infinite systems where the coefficient matrix dimension is not a priori fixed and a subsequent size enlargement technique is used. A special version of the algorithm, well suited to deal with block Hessenberg matrices having also a block band structure, is presented. Theoretical asymptotic bounds on the computational costs are proved. They indicate that, under suitable sparsity conditions, the proposed algorithms outperform Gaussian elimination. Numerical experiments have been carried out, showing the effectiveness of the algorithms when the size of the system is of practical interest. Copyright © 2004 John Wiley & Sons, Ltd.     

Comparison theorems for splittings of M-matrices in (block) Hessenberg form

BIT Numerical Mathematics - 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...     

A method of congruent type for linear systems with conjugate-normal coefficient matrices

Minimal residual methods, such as MINRES and GMRES, are well-known iterative versions of direct procedures for reducing a matrix to special condensed forms. The method of reduction used in these procedures is a sequence of unitary similarity transformations, while the condensed form is a tridiagonal matrix (MINRES) or a Hessenberg matrix (GMRES). The algorithm CSYM proposed in the 1990s for solving systems with complex symmetric matrices was based on the tridiagonal reduction performed via unitary congruences rather than similarities. In this paper, we construct an extension of this algorithm to the entire class of conjugate-normal matrices. (Complex symmetric matrices are a part of this class.) Numerical results are presented. They show that, on many occasions, the proposed algorithm has a superior convergence rate compared to GMRES.     

求解一类分块二阶线性方程组的QHSS迭代方法

本文将QHSS迭代方法运用于求解一类分块二阶线性方程组.通过适当地放宽QHSS迭代方法的收敛性条件,我们给出了用QHSS迭代方法求解一类分块二阶线性方程组的具体迭代格式,并证明了当系数矩阵中的(1,1)块对称半正定时该QHSS迭代方法的收敛性.我们还用数值实验验证了QHSS迭代方法的可行性和有效性.     

Spectral methods for orthogonal rational functions

We present an operator theoretic approach to orthogonal rational functions based on the identification of a suitable matrix representation of the multiplication operator associated with the corresponding orthogonality measure. Two alternatives are discussed, leading to representations which are linear fractional transformations with matrix coefficients acting on infinite Hessenberg or five-diagonal unitary matrices. This approach permits us to recover the orthogonality measure throughout the spectral analysis of an infinite matrix depending uniquely on the poles and the parameters of the recurrence relation for the orthogonal rational functions. Besides, the zeros of the orthogonal and para-orthogonal rational functions are identified as the eigenvalues of matrix linear fractional transformations of finite Hessenberg or five-diagonal matrices. As an application we use operator perturbation theory results to obtain new relations between the support of the orthogonality measure and the location of the poles and parameters of the recurrence relation for the orthogonal rational functions.     

An inequality for nonnegative matrices and the inverse eigenvalue problem

We present two versions of the same inequality, relating the maximal diagonal entry of a nonnegative matrix to its eigenvalues. We demonstrate a matrix factorization of a companion matrix, which leads to a solution of the nonnegative inverse eigenvalue problem (denoted the nniep) for 4×4 matrices of trace zero, and we give some sufficient conditions for a solution to the nniep for 5×5 matrices of trace zero. We also give a necessary condition on the eigenvalues of a 5×5 trace zero nonnegative matrix in lower Hessenberg form. Finally, we give a brief discussion of the nniep in restricted cases.     

Necessary and sufficient conditions for GMRES complete and partial stagnation

In this paper we give necessary and sufficient conditions for the complete or partial stagnation of the GMRES iterative method for solving real linear systems. Our results rely on a paper by Arioli, Pták and Strakoš (1998), characterizing the matrices having a prescribed convergence curve for the residual norms. We show that we have complete stagnation if and only if the matrix A is orthonormally similar to an upper or lower Hessenberg matrix having a particular first row or column or a particular last row or column. Partial stagnation is characterized by a particular pattern of the matrix Q in the QR factorization of the upper Hessenberg matrix generated by the Arnoldi process.