We obtain several results which may be useful in determining the convergence behavior of eigenvalue algorithms based upo n Arnoldi and nonsymmetric Lanczos recursions. We derive a relationship between nonsymmetric Lanczos eigenvalue procedures and Arnoldi eigenvalue procedures. We demonstrate that the Arnoldi recursions preserve a property which characterizes normal matrices, and that if we could determine the appropriate starting vectors, we could mimic the nonsymmetric Lanczos eigenvalue convergence on a general diagonalizable matrix by its convergence on related normal matrices. Using a unitary equivalence for each of these Krylov subspace methods, we define sets of test problems where we can easily vary certain spectral properties of the matrices. We use these and other test problems to examine the behavior of an Arnoldi and of a nonsymmetric Lanczos procedure. Mathematical Sciences Department, IBM Research Division, T.J. Watson Research Center, Yorktown Heights, NY 10598, USA, a...

Abstract. In this paper we analyze the bi-conjugate gradient algorithm in finite precision arithmetic, and suggest reasons for its often observed robustness. By using a tridiagonal structure, which is preserved by the finite precision bi-conjugate gradient iteration, we are able to bound its residual norm by a minimum polynomial of a perturbed matrix (i.e. the residual norm of the exact GMRES applied to a perturbed matrix) multiplied by an amplification factor. This shows that occurrence of near-breakdowns or loss of biorthogonality does not necessarily deter convergence of the residuals provided that the amplification factor remains bounded. Numerical examples are given to gain insights into these bounds. 1.

. 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 the n-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 GramSchmidt (IMGSA), and iterated classical Gram-Schmidt (ICGSA) implementations. Under the obvious assumption on the numerical nonsingularity of the system m...

The standard formulation of the conjugate gradient algorithm involves two inner product computations. The results of these two inner products are needed to update the search direction and the computed solution. Since these inner products are mutually interdependent, in a distributed memory parallel environment their computation and subsequent distribution requires two separate communication and synchronization phases. In this paper, we present three related mathematically equivalent rearrangements of the standard algorithm that reduce the number of communication phases. We present empirical evidence that two of these rearrangements are numerically stable. This claim is further substantiated by a proof that one of the empirically stable rearrangements arises naturally in the symmetric Lanczos method for linear systems, which is equivalent to the conjugate gradient method. 1 Introduction The conjugate gradient (CG) method is an effective iterative method for solving large s...

Sparsity is a fundamental concept of modern statistics, and often the only general principle available at the moment to address novel learning applications with many more variables than observations. While much progress has been made recently in the theoretical understanding and algorithmics of sparse point estimation, higher-order problems such as covariance estimation or optimal data acquisition are seldomly addressed for sparsity-favouring models, and there are virtually no algorithms for large scale applications of these. We provide novel approximate Bayesian inference algorithms for sparse generalized linear models, that can be used with hundred thousands of variables, and run orders of magnitude faster than previous algorithms in domains where either apply. By analyzing our methods and establishing some novel convexity results, we settle a long-standing open question about variational Bayesian inference for continuous variable models: the Gaussian lower bound relaxation, which has been used previously for a range of models, is proved to be a convex optimization problem, if and only if the posterior mode is found by convex programming. Our algorithms reduce to the same computational primitives than commonly used sparse estimation methods do, but require Gaussian marginal variance estimation as well. We show how the Lanczos algorithm from numerical mathematics can be employed to compute the latter. We are interested in Bayesian experimental design here (which is mainly driven by efficient approximate inference), a powerful framework for optimizing measurement architectures of complex signals, such as natural images. Designs

Many algorithms for solving linear systems, least squares problems or eigenproblems need to compute an orthonormal basis. The computation is commonly performed using a QR factorization computed using the classical or the modified Gram-Schmidt algorithm, the Householder algorithm, the Givens algorithm or the Gram-Schmidt algorithm with iterative reorthogonalization. For linear systems, it is well-known that the backward stability of the process depends crucially on the algorithm used for the QR factorization. For the eigenproblem, although textbooks warn users about the possible instability of eigensolvers due to loss of orthogonality, few theoretical results exist. In this paper we show that the loss of orthogonality of the computed basis can affect the reliability of the computed eigenpair. We also show that the classical residual norm kAx \Gamma xk and the specific computed one using a Krylov method can differ because of the loss of orthogonality of the computed basis of the Krylov s...

An explicit procedure is presented for computing both model and data resolution matrices within a PaigeSaunders LSQR algorithm for iterative inversion in seismic tomography. These methods are designed to avoid the need for an additional singular value decomposition of the ray-path matrix. The techniques discussed are completely general since they are based on the multiplicity of equivalent exact formulas that may be used to define the resolution matrices. Thus, resolution matrices may also be computed for a wide variety of iterative inversion algorithms using the same ideas. Keywords: seismic tomography, inversion, resolution matrices 1 INTRODUCTION Linear tomographic reconstruction schemes have well-defined resolution properties. 1 Yet, the resolution matrices summarizing these properties are not computed as often as they might be --- at least in part --- because of the common misconception that resolution matrices can only be found using singular value decomposition (SVD). Since ...

The Generalized minimal residual method (GMRES) is known as an efficient iterative method for solving large nonsymmetric systems of linear equations. In this thesis, we study numerical stability of the GMRES method. For the construction of the Arnoldi basis, we consider the Householder orthogonalization and the frequently used modified Gram-Schmidt process. While for the more expensive Householder implementation the orthogonality of the computed basis is preserved close to the machine precision level, for the modified Gram-Schmidt Arnoldi process the computed vectors gradually lose their orthogonality. Using the bound on the loss of orthogonality, it is proved that, under certain assumptions on the numerical nonsingularity of the system matrix, the GMRES implementation based on the Householder orthogonalization is backward stable. It produces an approximate solution with the residual which is of the same order as that one obtained from the direct solving of the system Ax = b by the Ho...

SVDPACKC comprises four numerical (iterative) methods for computing the singular value decomposition (SVD) of large sparse matrices using ANSI C. This software package implements Lanczos and subspace iteration-based methods for determining several of the largest singular triplets (singular values and corresponding left- and right-singular vectors) for large sparse matrices. The package has been ported to a variety of machines ranging from supercomputers to workstations: CRAY Y-MP, IBM RS/6000-550, DEC 5000100, HP 9000-750, SPARCstation 2, and Macintosh II/fx. This document (i) explains each algorithm in some detail, (ii) explains the input parameters for each program, (iii) explains how to compile/execute each program, and (iv) illustrates the performance of each method when we compute lower rank approximations to sparse term-document matrices from information retrieval applications. A user-friendly software interface to the package for UNIX-based systems and the Macintosh II/fx is als...

We present an e#cient and accurate variant of the conjugate gradient method for solving families of shifted systems. In particular we are interested in shifted systems that occur in Tikhonov regularization for inverse problems since these problems can be sensitive to roundo# errors. The success of our method in achieving accurate approximations is supported by theoretical arguments as well as several numerical experiments and we relate it to other implementations proposed in literature.