Numerical stability of GMRES

The Generalized Minimal Residual Method (GMRES) is one of the significant methods for solving linear algebraic systems with nonsymmetric matrices. It minimizes the norm of the residual on the linear variety determined by the initial residual and then-th Krylov residual subspace and is therefore optimal, with respect to the size of the residual, in the class of Krylov subspace methods. One possible way of computing the GMRES approximations is based on constructing the orthonormal basis of the Krylov subspaces (Arnoldi basis) and then solving the transformed least squares problem. This paper studies the numerical stability of such formulations of GMRES. Our approach is based on the Arnoldi recurrence for the actually, i.e., in finite precision arithmetic, computed quantities. We consider the Householder (HHA), iterated modified Gram-Schmidt (IMGSA), and iterated classical Gram-Schmidt (ICGSA) implementations. Under the obvious assumption on the numerical nonsingularity of the system matrix, the HHA implementation of GMRES is proved backward stable in the normwise sense. That is, the backward error for the approximation is proportional to machine precision . Additionally, it is shown that in most cases the norm of the residual computed from the transformed least squares problem (Arnoldi residual) gives a good estimate of the true residual norm, until the true residual norm has reached the level A x.This work was supported by NSF contract Int921824.     

A simpler GMRES

The generalized minimal residual (GMRES) method is widely used for solving very large, nonsymmetric linear systems, particularly those that arise through discretization of continuous mathematical models in science and engineering. By shifting the Arnoldi process to begin with Ar₀ instead of r₀, we obtain simpler Gram–Schmidt and Householder implementations of the GMRES method that do not require upper Hessenberg factorization. The Gram–Schmidt implementation also maintains the residual vector at each iteration, which allows cheaper restarts of GMRES(m) and may otherwise be useful.     

A Newton basis GMRES implementation

The GMRES method by Saad and Schultz is one of the most populariterative methods for the solution of large sparse non-symmetriclinear systems of equations. The implementation proposed bySaad and Schultz uses the Arnoldi process and the modified Gram-Schmidt(MGS) method to compute orthonormal bases of certain Krylovsubspaces. The MGS method requires many vector-vector operations,which can be difficult to implement efficiently on vector andparallel computers due to the low granularity of these operations.We present a new implementation of the GMRES method in which,for each Krylov subspace used, we first determine a Newton basis,and then orthogonalize it by computing a QR factorization ofthe matrix whose columns are the vectors of the Newton basis.In this way we replace the vector-vector operations of the MGSmethod by the task of computing a QR factorization of a densematrix. This makes the implementa tion more flexible, and providesa possibility to adapt the computations to the computer at handin order to achieve better performance. E-mail addresses: bai@ms.uky.edu dyhu@rice.edu reichel@mcs.kent.edu Research supported in part by the IBM Bergen Scientific Center,the Center for Research on Parallel Computation at Rice University,a National Research Council fellowship and NSF Grant DMS-9002884.     

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.     

Entire vectors and time reversible Cauchy problems

Our major aim in this paper is to show the algebraic and topological equivalence between the set of all entire vectors ε(A) for a linear operatorA and the solution space of the corresponding abstract Cauchy problem, when—A is the generator of a strongly continuous holomorphic semigroup.     

A new adaptive GMRES algorithm for achieving high accuracy

GMRES(k) is widely used for solving non-symmetric linear systems. However, it is inadequate either when it converges only for k close to the problem size or when numerical error in the modified Gram–Schmidt process used in the GMRES orthogonalization phase dramatically affects the algorithm performance. An adaptive version of GMRES(k) which tunes the restart value k based on criteria estimating the GMRES convergence rate for the given problem is proposed here. This adaptive GMRES(k) procedure outperforms standard GMRES(k), several other GMRES-like methods, and QMR on actual large scale sparse structural mechanics postbuckling and analog circuit simulation problems. There are some applications, such as homotopy methods for high Reynolds number viscous flows, solid mechanics postbuckling analysis, and analog circuit simulation, where very high accuracy in the linear system solutions is essential. In this context, the modified Gram–Schmidt process in GMRES, can fail causing the entire GMRES iteration to fail. It is shown that the adaptive GMRES(k) with the orthogonalization performed by Householder transformations succeeds whenever GMRES(k) with the orthogonalization performed by the modified Gram–Schmidt process fails, and the extra cost of computing Householder transformations is justified for these applications. © 1998 John Wiley & Sons, Ltd.     

On pointwise and analytic similarity of matrices

LetA(ε) andB(ε) be complex valued matrices analytic in ε at the origin.A(ε)≈ _p B(ε) ifA(ε) is similar toB(ε) for any |ε|<r,A(ε)≈_a B(ε) ifB(ε)=T(ε)A(ε)T ⁻¹(ε) andT(ε) is analytic and |T(ε)|≠0 for |ε|<r! In this paper we find a necessary and sufficient conditions onA(ε) andB(ε) such thatA(ε)≈ _a B(ε) provided thatA(ε)≈ _p B(ε). This problem arises in study of certain ordinary differential equations singular with respect to a parameter ε in the origin and was first stated by Wasow. Sponsored by the United States Army under Contract No. DAAG29-75-C-0024     

CMRH: A new method for solving nonsymmetric linear systems based on the Hessenberg reduction algorithm

The Generalized Minimal Residual (GMRES) method and the Quasi-Minimal Residual (QMR) method are two Krylov methods for solving linear systems. The main difference between these methods is the generation of the basis vectors for the Krylov subspace. The GMRES method uses the Arnoldi process while QMR uses the Lanczos algorithm for constructing a basis of the Krylov subspace. In this paper we give a new method similar to QMR but based on the Hessenberg process instead of the Lanczos process. We call the new method the CMRH method. The CMRH method is less expensive and requires slightly less storage than GMRES. Numerical experiments suggest that it has behaviour similar to GMRES. This revised version was published online in June 2006 with corrections to the Cover Date.     

On the loss of orthogonality in the second-order Arnoldi process

This note is devoted to the rounding error analysis of the second-order Arnoldi process for constructing an orthonormal basis of the second-order Krylov subspace. The effect of the rounding errors on the orthogonality of the derived vectors is given.     

A hybrid Arnoldi-Faber iterative method for nonsymmetric systems of linear equations

Summary We present here a new hybrid method for the iterative solution of large sparse nonsymmetric systems of linear equations, say of the formAx=b, whereA ^{N, N} , withA nonsingular, andb ^N are given. This hybrid method begins with a limited number of steps of the Arnoldi method to obtain some information on the location of the spectrum ofA, and then switches to a Richardson iterative method based on Faber polynomials. For a polygonal domain, the Faber polynomials can be constructed recursively from the parameters in the Schwarz-Christoffel mapping function. In four specific numerical examples of non-normal matrices, we show that this hybrid algorithm converges quite well and is approximately as fast or faster than the hybrid GMRES or restarted versions of the GMRES algorithm. It is, however, sensitive (as other hybrid methods also are) to the amount of information on the spectrum ofA acquired during the first (Arnoldi) phase of this procedure.     

Memory‐efficient Arnoldi algorithms for linearizations of matrix polynomials in Chebyshev basis

Novel memory‐efficient Arnoldi algorithms for solving matrix polynomial eigenvalue problems are presented. More specifically, we consider the case of matrix polynomials expressed in the Chebyshev basis, which is often numerically more appropriate than the standard monomial basis for a larger degree d. The standard way of solving polynomial eigenvalue problems proceeds by linearization, which increases the problem size by a factor d. Consequently, the memory requirements of Krylov subspace methods applied to the linearization grow by this factor. In this paper, we develop two variants of the Arnoldi method that build the Krylov subspace basis implicitly, in a way that only vectors of length equal to the size of the original problem need to be stored. The proposed variants are generalizations of the so‐called quadratic Arnoldi method and two‐level orthogonal Arnoldi procedure methods, which have been developed for the monomial case. We also show how the typical ingredients of a full implementation of the Arnoldi method, including shift‐and‐invert and restarting, can be incorporated. Numerical experiments are presented for matrix polynomials up to degree 30 arising from the interpolation of nonlinear eigenvalue problems, which stem from boundary element discretizations of PDE eigenvalue problems. Copyright © 2013 John Wiley & Sons, Ltd.     

A global harmonic Arnoldi method for large non-Hermitian eigenproblems with an application to multiple eigenvalue problems

The global Arnoldi method can be used to compute exterior eigenpairs of a large non-Hermitian matrix A, but it does not work well for interior eigenvalue problems. Based on the global Arnoldi process that generates an F-orthonormal basis of a matrix Krylov subspace, we propose a global harmonic Arnoldi method for computing certain harmonic F-Ritz pairs that are used to approximate some interior eigenpairs. We propose computing the F-Rayleigh quotients of the large non-Hermitian matrix with respect to harmonic F-Ritz vectors and taking them as new approximate eigenvalues. They are better and more reliable than the harmonic F-Ritz values. The global harmonic Arnoldi method inherits convergence properties of the harmonic Arnoldi method applied to a larger matrix whose distinct eigenvalues are the same as those of the original given matrix. Some properties of the harmonic F-Ritz vectors are presented. As an application, assuming that A is diagonalizable, we show that the global harmonic Arnoldi method is able to solve multiple eigenvalue problems both in theory and in practice. To be practical, we develop an implicitly restarted global harmonic Arnoldi algorithm with certain harmonic F-shifts suggested. In particular, this algorithm can be adaptively used to solve multiple eigenvalue problems. Numerical experiments show that the algorithm is efficient for the eigenproblem and is reliable for quite ill-conditioned multiple eigenproblems.     

A numerical procedure for computing the Moore-Penrose inverse

An algorithm for computing the Moore-Penrose inverse of an arbitraryn×m real matrixA is presented which uses a Gram-Schmidt like procedure to form anA-orthogonal set of vectors which span the subspace perpendicular to the kernel ofA. This one procedure will work for any value ofn andm, and for any value of rank (A).     

± sign pattern matrices that allow orthogonality

A sign pattern A is a ± sign pattern if A has no zero entries. A allows orthogonality if there exists a real orthogonal matrix B whose sign pattern equals A. Some sufficient conditions are given for a sign pattern matrix to allow orthogonality, and a complete characterization is given for ± sign patterns with n − 1 ⩽ N₋(A) ⩽ n + 1 to allow orthogonality.     

Accurate downdating of a modified Gram-Schmidt QR decomposition

A new algorithm for downdating a QR decomposition is presented. We show that, when the columns in the Q factor from the Modified Gram-Schmidt QR decomposition of a matrixX are exactly orthonormal, the Gram-Schmidt downdating algorithm for the QR decomposition ofX is equivalent to downdating the full Householder QR decomposition of the matrixX augmented by ann ×n zero matrix on top. Using this relation, we derive an algorithm that improves the Gram-Schmidt downdating algorithm when the columns in the Q factor are not orthonormal. Numerical test results show that the new algorithm produces far more accurate results than the Gram-Schmidt downdating algorithm for certain ill-conditioned problems.This work was partially supported in part by the National Science Foundation grants CCR-9209726 and CCR-9509085.     

On the modified Gram-Schmidt algorithm for weighted and constrained linear least squares problems

A framework and an algorithm for using modified Gram-Schmidt for constrained and weighted linear least squares problems is presented. It is shown that a direct implementation of a weighted modified Gram-Schmidt algorithm is unstable for heavily weighted problems. It is shown that, in most cases it is possible to get a stable algorithm by a simple modification free from any extra computational costs. In particular, it is not necessary to perform reorthogonalization.Solving the weighted and constrained linear least squares problem with the presented weighted modified Gram-Schmidt algorithm is seen to be numerically equivalent to an algorithm based on a weighted Householder-likeQR factorization applied to a slightly larger problem. This equivalence is used to explain the instability of the weighted modified Gram-Schmidt algorithm. If orthogonality, with respect to a weighted inner product, of the columns inQ is important then reorthogonalization can be used. One way of performing such reorthogonalization is described.Computational tests are given to show the main features of the algorithm.     

Equivalence between modified symplectic Gram-Schmidt and Householder SR algorithms

The SR factorization for a given matrix A is a QR-like factorization A=SR, where the matrix S is symplectic and R is J-upper triangular. This factorization is fundamental for some important structure-preserving methods in linear algebra and is usually implemented via the symplectic Gram-Schmidt algorithm (SGS). There exist two versions of SGS, the classical (CSGS) and the modified (MSGS). Both are equivalent in exact arithmetic, but have very different numerical behaviors. The MSGS is more stable. Recently, the symplectic Householder SR algorithm has been introduced, for computing efficiently the SR factorization. In this paper, we show two new and important results. The first is that the SR factorization of a matrix A via the MSGS is mathematically equivalent to the SR factorization via Householder SR algorithm of an embedded matrix. The later is obtained from A by adding two blocks of zeros in the top of the first half and in the top of the second half of the matrix A. The second result is that MSGS is also numerically equivalent to Householder SR algorithm applied to the mentioned embedded matrix.     

The rate of convergence of GMRES on a tridiagonal Toeplitz linear system

The Generalized Minimal Residual method (GMRES) is often used to solve a nonsymmetric linear system Ax = b. But its convergence analysis is a rather difficult task in general. A commonly used approach is to diagonalize A = XΛ X ⁻¹ and then separate the study of GMRES convergence behavior into optimizing the condition number of X and a polynomial minimization problem over A’s spectrum. This artificial separation could greatly overestimate GMRES residuals and likely yields error bounds that are too far from the actual ones. On the other hand, considering the effects of both A’s spectrum and the conditioning of X at the same time poses a difficult challenge, perhaps impossible to deal with in general but only possible for certain particular linear systems. This paper will do so for a (nonsymmetric) tridiagonal Toeplitz system. Sharp error bounds on and sometimes exact expressions for residuals are obtained. These expressions and/or bounds are in terms of the three parameters that define A and Chebyshev polynomials of the first kind.     

Iterative refinement of linear least squares solutions I

An iterative procedure is developed for reducing the rounding errors in the computed least squares solution to an overdetermined system of equationsAx =b, whereA is anm ×n matrix (m n) of rankn. The method relies on computing accurate residuals to a certain augmented system of linear equations, by using double precision accumulation of inner products. To determine the corrections, two methods are given, based on a matrix decomposition ofA obtained either by orthogonal Householder transformations or by a modified Gram-Schmidt orthogonalization. It is shown that the rate of convergence in the iteration is independent of the right hand side,b, and depends linearly on the condition number, 2135;(A), of the rectangular matrixA. The limiting accuracy achieved will be approximately the same as that obtained by a double precision factorization.In a second part of this paper the case whenx is subject to linear constraints and/orA has rank less thann is covered. Here also ALGOL-programs embodying the derived algorithms will be given.This work was sponsored by the Swedish Natural Science Research Council.     

Asymptotic analysis of a parabolic semilinear problem with nonlinear boundary multiphase interactions in a perforated domain

We consider a parabolic semilinear problem with rapidly oscillating coefficients in a domain Ω_ε that is ε-periodically perforated by small holes of size O\mathcal {O}(ε). The holes are divided into two ε-periodical sets depending on the boundary interaction at their surfaces, and two different nonlinear Robin boundary conditions σ_ε(u _ε) + εκ _m (u _ε) = εg ^(m) _ε, m = 1, 2, are imposed on the boundaries of holes. We study the asymptotics as ε → 0 and establish a convergence theorem without using extension operators. An asymptotic approximation of the solution and the corresponding error estimate are also obtained. Bibliography: 60 titles. Illustrations: 1 figure.