QMR: a quasi-minimal residual method for non-Hermitian linear systems

Summary The biconjugate gradient (BCG) method is the natural generalization of the classical conjugate gradient algorithm for Hermitian positive definite matrices to general non-Hermitian linear systems. Unfortunately, the original BCG algorithm is susceptible to possible breakdowns and numerical instabilities. In this paper, we present a novel BCG-like approach, the quasi-minimal residual (QMR) method, which overcomes the problems of BCG. An implementation of QMR based on a look-ahead version of the nonsymmetric Lanczos algorithm is proposed. It is shown how BCG iterates can be recovered stably from the QMR process. Some further properties of the QMR approach are given and an error bound is presented. Finally, numerical experiments are reported.This work was supported in part by DARPA via Cooperative Agreement NCC 2-387 between NASA and the Universities Space Research Association (USRA).     

Computation of matrix-valued formally orthogonal polynomials and applications

We present a computational procedure for generating formally orthogonal polynomials associated with a given bilinear Hankel form with rectangular matrix-valued moments. Our approach covers the most general case of moments of any size and is not restricted to square moments. Moreover, our algorithm has a built-in deflation procedure to handle linearly dependent or almost linearly dependent columns and rows of the block Hankel matrix associated with the bilinear form. Possible singular or close-to-singular leading principal submatrices of the deflated block Hankel matrix are avoided by means of look-ahead techniques. Applications of the computational procedure to eigenvalue computations, reduced-order modeling, the solution of multiple linear systems, and the fast solution of block Hankel systems are also briefly described.     

A look-ahead Bareiss algorithm for general Toeplitz matrices

Summary. The Bareiss algorithm is one of the classical fast solvers for systems of linear equations with Toeplitz coefficient matrices. The method takes advantage of the special structure, and it computes the solution of a Toeplitz system of order~ with only~ arithmetic operations, instead of~ operations. However, the original Bareiss algorithm requires that all leading principal submatrices be nonsingular, and the algorithm is numerically unstable if singular or ill-conditioned submatrices occur. In this paper, an extension of the Bareiss algorithm to general Toeplitz systems is presented. Using look-ahead techniques, the proposed algorithm can skip over arbitrary blocks of singular or ill-conditioned submatrices, and at the same time, it still fully exploits the Toeplitz structure. Implementation details and operations counts are given, and numerical experiments are reported. We also discuss special versions of the proposed look-ahead Bareiss algorithm for Hermitian indefinite Toeplitz systems and banded Toeplitz systems. Received August 27, 1993 / Revised version received March 1994     

On conjugate gradient type methods and polynomial preconditioners for a class of complex non-hermitian matrices

Summary We consider conjugate gradient type methods for the solution of large linear systemsA x=b with complex coefficient matrices of the typeA=T+iI whereT is Hermitian and a real scalar. Three different conjugate gradient type approaches with iterates defined by a minimal residual property, a Galerkin type condition, and an Euclidean error minimization, respectively, are investigated. In particular, we propose numerically stable implementations based on the ideas behind Paige and Saunder's SYMMLQ and MINRES for real symmetric matrices and derive error bounds for all three methods. It is shown how the special shift structure ofA can be preserved by using polynomial preconditioning, and results on the optimal choice of the polynomial preconditioner are given. Also, we report on some numerical experiments for matrices arising from finite difference approximations to the complex Helmholtz equation.This work was supported in part by Cooperatives Agreement NCC 2-387 between the National Aeronautics and Space Administration (NASA) and the Universities Space Research Association (USRA) and by National Science Foundation Grant DCR-8412314     

Look-ahead Levinson and Schur algorithms for non-Hermitian Toeplitz systems

Summary. We present generalizations of the nonsymmetric Levinson and Schur algorithms for non-Hermitian Toeplitz matrices with some singular or ill-conditioned leading principal submatrices. The underlying recurrences allow us to go from any pair of successive well-conditioned leading principal submatrices to any such pair of larger order. If the look-ahead step size between these pairs is bounded, our generalized Levinson and Schur recurrences require $ operations, and the Schur recurrences can be combined with recursive doubling so that an $ algorithm results. The overhead (in operations and storage) of look-ahead steps is very small. There are various options for applying these algorithms to solving linear systems with Toeplitz matrix. Received July 26, 1993     

Recurrence relation for biorthogonal polynomials

We use a geometric approach to obtain a recurrence relation for two families of biorthogonal polynomials associated to a nonsingular, strongly regular matrix M. We propose a “look-ahead procedure” for computing the biorthogonal polynomials when M has singular or ill-conditioned leading principal submatrices. These polynomials lead to two recursive triangular factorizations for the inverse of a nonsingular matrix M which is not necessarily strongly regular.     

The stable computation of formal orthogonal polynomials

For many applications — such as the look-ahead variants of the Lanczos algorithm — a sequence of formal (block-)orthogonal polynomials is required. Usually, one generates such a sequence by taking suitable polynomial combinations of a pair of basis polynomials. These basis polynomials are determined by a look-ahead generalization of the classical three term recurrence, where the polynomial coefficients are obtained by solving a small system of linear equations. In finite precision arithmetic, the numerical orthogonality of the polynomials depends on a good choice of the size of the small systems; this size is usually controlled by a heuristic argument such as the condition number of the small matrix of coefficients. However, quite often it happens that orthogonality gets lost.We present a new variant of the Cabay-Meleshko algorithm for numerically computing pairs of basis polynomials, where the numerical orthogonality is explicitly monitored with the help of stability parameters. A corresponding error analysis is given. Our stability parameter is shown to reflect the condition number of the underlying Hankel matrix of moments. This enables us to prove the weak and strong stability of our method, provided that the corresponding Hankel matrix is well-conditioned.This work was partially supported by the HCM project ROLLS, under contract CHRX-CT93-0416.     

Chebyshev rational interpolation

The Lanczos method and its variants can be used to solve efficiently the rational interpolation problem. In this paper we present a suitable fast modification of a general look-ahead version of the Lanczos process in order to deal with polynomials expressed in the Chebyshev orthogonal basis. The proposed approach is particularly suited for approximating analytic functions by means of rational interpolation at certain nodes located on the boundary of an elliptical region of the complex plane. In fact, in this case it overcomes some of the numerical difficulties which limited the applicability of the look-ahead Lanczos process for determining the coefficients both of the numerators and of the denominators with respect to the standard power basis. This revised version was published online in August 2006 with corrections to the Cover Date.     

A linear system solver based on a modified Krylov subspace method for breakdown recovery

Despite its usefulness in solving eigenvalue problems and linear systems of equations, the nonsymmetric Lanczos method is known to suffer from a potential breakdown problem. Previous and recent approaches for handling the Lanczos exact and near-breakdowns include, for example, the look-ahead schemes by Parlett-Taylor-Liu [23], Freund-Gutknecht-Nachtigal [9], and Brezinski-Redivo Zaglia-Sadok [4]; the combined look-ahead and restart scheme by Joubert [18]; and the low-rank modified Lanczos scheme by Huckle [17]. In this paper, we present yet another scheme based on a modified Krylov subspace approach for the solution of nonsymmetric linear systems. When a breakdown occurs, our approach seeks a modified dual Krylov subspace, which is the sum of the original subspace and a new Krylov subspaceK _m (w _j ,A ^T ) wherew _j is a newstart vector (this approach has been studied by Ye [26] for eigenvalue computations). Based on this strategy, we have developed a practical algorithm for linear systems called the MLAN/QM algorithm, which also incorporates the residual quasi-minimization as proposed in [12]. We present a few convergence bounds for the method as well as numerical results to show its effectiveness.Research supported by Natural Sciences and Engineering Research Council of Canada.     

Immanant positivity for Catalan-Stieltjes matrices

We give some sufficient conditions for the nonnegativity of immanants of square submatrices of Catalan-Stieltjes matrices and their corresponding Hankel matrices. To obtain these sufficient conditions, we construct new planar networks with a recursive nature for Catalan-Stieltjes matrices. As applications, we provide a unified way to produce inequalities for many combinatorial polynomials, such as the Eulerian polynomials, Schröder polynomials, and Narayana polynomials.     

Look-ahead methods for block Hankel systems

In this paper, the updating formulas used by three look-ahead methods for solving Hankel systems are generalized to the square block case. Each of the original methods was described in the literature using the terminology of different but strongly related fields: formal orthogonal polynomials, Padé approximants, structured matrices. This paper gives several of these connections generalized to the block case and shows that each viewpoint has its own merits.     

Fast algorithms for Toeplitz and Hankel matrices

The paper gives a self-contained survey of fast algorithms for solving linear systems of equations with Toeplitz or Hankel coefficient matrices. It is written in the style of a textbook. Algorithms of Levinson-type and Schur-type are discussed. Their connections with triangular factorizations, Padè recursions and Lanczos methods are demonstrated. In the case in which the matrices possess additional symmetry properties, split algorithms are designed and their relations to butterfly factorizations are developed.     

Lanczos-Based Exponential Filtering for Discrete Ill-Posed Problems

We describe regularizing iterative methods for the solution of large ill-conditioned linear systems of equations that arise from the discretization of linear ill-posed problems. The regularization is specified by a filter function of Gaussian type. A parameter determines the amount of regularization applied. The iterative methods are based on a truncated Lanczos decomposition and the filter function is approximated by a linear combination of Lanczos polynomials. A suitable value of the regularization parameter is determined by an L-curve criterion. Computed examples that illustrate the performance of the methods are presented.     

Gallai–Edmonds structure theorem for weighted matching polynomial

In this paper, we prove the Gallai–Edmonds structure theorem for weighted matching polynomials. Our result implies the Parter–Wiener theorem and its recent generalization about the existence of principal submatrices of a Hermitian matrix whose graph is a tree.     

Modified moments for indefinite weight functions

Summary The problem of generating the recurrence coefficients of orthogonal polynomials from the moments or from modified moments of the weight function is well understood for positive weight distributions. Here we extend this theory and the basic algorithms to the case of an indefinite weight function. While the generic indefinite case is formally not much different from the positive definite case, there exist nongeneric degenerate situations, and these require a different more complicated treatment. The understanding of these degenerate situations makes it possible to construct a stable approximate solution of an ill-conditioned problem.The application to adaptive iterative methods for linear systems of equations is anticipated.Dedicated to R.S. Varga on the occasion of his sixtieth birthdayG.H. Golub (USA, Canada), M.H. Gutknecht (other countries)     

Two Kinds of Numbers and Their Applications

C. Radoux （J. Comput. Appl. Math., 115 （2000） 471-477） obtained a computational formula of Hankel determinants on some classical combinatorial sequences such as Catalan numbers and polynomials, Bell polynomials, Hermite polynomials, Derangement polynomials etc. From a pair of matrices this paper introduces two kinds of numbers. Using the first kind of numbers we give a unified treatment of Hankel determinants on those sequences, i.e., to consider a general representation of Hankel matrices on the first kind of numbers. It is interesting that the Hankel determinant of the first kind of numbers has a close relation that of the second kind of numbers.     

Backward stability of a pivoting strategy for sign-regular linear systems

A matrixA issign-regular if, for each orderk, allk×k submatrices ofA have determinant with the same sign. In this paper, a pivoting strategy ofO(n) operations for the Gaussian elimination of linear systems whose coefficient matrices are sign-regular is proposed. Backward error analysis of this pivoting strategy is performed and small error bounds are obtained. Our results can also be applied to linear systems whose coefficient matrices have sign-regular inverses.     

GPBi‐CGstab(L): A Lanczos‐type product method unifying Bi‐CGstab(L) and GPBi‐CG

Lanczos‐type product methods (LTPMs), in which the residuals are defined by the product of stabilizing polynomials and the Bi‐CG residuals, are effective iterative solvers for large sparse nonsymmetric linear systems. Bi‐CGstab(L) and GPBi‐CG are popular LTPMs and can be viewed as two different generalizations of other typical methods, such as CGS, Bi‐CGSTAB, and Bi‐CGStab2. Bi‐CGstab(L) uses stabilizing polynomials of degree L, while GPBi‐CG uses polynomials given by a three‐term recurrence (or equivalently, a coupled two‐term recurrence) modeled after the Lanczos residual polynomials. Therefore, Bi‐CGstab(L) and GPBi‐CG have different aspects of generalization as a framework of LTPMs. In the present paper, we propose novel stabilizing polynomials, which combine the above two types of polynomials. The resulting method is referred to as GPBi‐CGstab(L). Numerical experiments demonstrate that our presented method is more effective than conventional LTPMs.     

Generalized biorthogonal bases and tridiagonalisation of matrices

Summary. In connection with the breakdown problem of the Lanczos algorithm a theory of generalized biorthogonal bases is developed. The connection between the generalized biorthogonal bases of Krylov chains and look-ahead Lanczos recursions is worked out in detail. It is shown how generalized biorthogonal bases with “antidiagonal blocks” can be constructed with small computational effort. Finally a special look-ahead Lanczos algorithm is derived which requires minimal computational effort and storage. Received September 21, 1995 / Revised version received August 12, 1996     

A fast solver for linear systems with displacement structure

We describe a fast solver for linear systems with reconstructible Cauchy-like structure, which requires O(rn ²) floating point operations and O(rn) memory locations, where n is the size of the matrix and r its displacement rank. The solver is based on the application of the generalized Schur algorithm to a suitable augmented matrix, under some assumptions on the knots of the Cauchy-like matrix. It includes various pivoting strategies, already discussed in the literature, and a new algorithm, which only requires reconstructibility. We have developed a software package, written in Matlab and C-MEX, which provides a robust implementation of the above method. Our package also includes solvers for Toeplitz(+Hankel)-like and Vandermonde-like linear systems, as these structures can be reduced to Cauchy-like by fast and stable transforms. Numerical experiments demonstrate the effectiveness of the software.