Approximating Runge-Kutta matrices by triangular matrices

The implementation of implicit Runge-Kutta methods requires the solution of large systems of non-linear equations. Normally these equations are solved by a modified Newton process, which can be very expensive for problems of high dimension. The recently proposed triangularly implicit iteration methods for ODE-IVP solvers [5] substitute the Runge-Kutta matrixA in the Newton process for a triangular matrixT that approximatesA, hereby making the method suitable for parallel implementation. The matrixT is constructed according to a simple procedure, such that the stiff error components in the numerical solution are strongly damped. In this paper we prove for a large class of Runge-Kutta methods that this procedure can be carried out and that the diagnoal entries ofT are positive. This means that the linear systems that are to be solved have a non-singular matrix. The research reported in this paper was supported by STW (Dutch Foundation for Technical Sciences).     

Generalization of Strang's Preconditioner with Applications to Toeplitz Least Squares Problems

In this paper, we propose a method to generalize Strang's circulant preconditioner for arbitrary n-by-n matrices A_n. The th column of our circulant preconditioner Sn is equal to the th column of the given matrix A_n. Thus if A_n is a square Toeplitz matrix, then Sn is just the Strang circulant preconditioner. When Sn is not Hermitian, our circulant preconditioner can be defined as . This construction is similar to the forward-backward projection method used in constructing preconditioners for tomographic inversion problems in medical imaging. We show that if the matrix A_n has decaying coefficients away from the main diagonal, then is a good preconditioner for An. Comparisons of our preconditioner with other circulant-based preconditioners are carried out for some 1-D Toeplitz least squares problems: min ∥ b - Ax∥₂. Preliminary numerical results show that our preconditioner performs quite well, in comparison to other circulant preconditioners. Promising test results are also reported for a 2-D deconvolution problem arising in ground-based atmospheric imaging.     

BTTB preconditioners for BTTB least squares problems

In this paper, we consider solving the least squares problem min_x‖b-Tx‖₂ by using preconditioned conjugate gradient (PCG) methods, where T is a large rectangular matrix which consists of several square block-Toeplitz-Toeplitz-block (BTTB) matrices and b is a column vector. We propose a BTTB preconditioner to speed up the PCG method and prove that the BTTB preconditioner is a good preconditioner. We then discuss the construction of the BTTB preconditioner. Numerical examples, including image restoration problems, are given to illustrate the efficiency of our BTTB preconditioner. Numerical results show that our BTTB preconditioner is more efficient than the well-known Level-1 and Level-2 circulant preconditioners.     

Applications of shortest path algorithms to matrix scalings

Summary A symmetric scaling of a nonnegative, square matrixA is a matrixXAX ^–1, whereX is a nonsingular, nonnegative diagonal matrix. By associating a family of weighted directed graphs with the matrixA we are able to adapt the shortest path algorithms to compute an optimal scaling ofA, where we call a symmetric scalingA ofA optimal if it minimizes the maximum of the ratio of non-zero elements.Dedicated to Professor F.L. Bauer on the occasion of his 60th birthdayThe first author was partially supported by the Deutsche Forschungsgemeinschaft under grant GO 270/3, the second author by the U.S. National Science Foundation under grand MCS-8026132     

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.     

Local refinement techniques for elliptic problems on cell-centered grids

Summary Algebraic multilevel analogues of the BEPS preconditioner designed for solving discrete elliptic problems on grids with local refinement are formulated, and bounds on their relative condition numbers, with respect to the composite-grid matrix, are derived. TheV-cycle and, more generally,v-foldV-cycle multilevel BEPS preconditioners are presented and studied. It is proved that for 2-D problems theV-cycle multilevel BEPS is almost optimal, whereas thev-foldV-cycle algebraic multilevel BEPS is optimal under a mild restriction on the composite cell-centered grid. For thev-fold multilevel BEPS, the variational relation between the finite difference matrix and the corresponding matrix on the next-coarser level is not necessarily required. Since they are purely algebraically derived, thev-fold (v>1) multilevel BEPS preconditioners perform without any restrictionson the shape of subregions, unless the refinement is too fast. For theV-cycle BEPS preconditioner (2-D problem), a variational relation between the matrices on two consecutive grids is required, but there is no restriction on the method of refinement on the shape, or on the size of the subdomains.     

A generalized successive overrelaxation method for least squares problems

In this paper a new iterative method is given for solving large sparse least squares problems and computing the minimum norm solution to underdetermined consistent linear systems. The new scheme is called the generalized successive overrelaxation (GSOR) method and is shown to be convergent ifA is full column rank. The GSOR method involves a parameter ρ and an auxiliary matrixP. One can choose matrix P so that the GSOR method only involves matrix and vector operations; therefore the GSOR method is suitable for parallel computations. Besides, the GSOR method can be combined with preconditioning techniques, and therefore can be expected to be more effective. This author's work was supported by Natural Science Foundation of Liaoning Province, China.     

Matrix enlarging methods and their application

This paper explores several methods for matrix enlarging, where an enlarged matrixà is constructed from a given matrixA. The methods explored include matrix primitization, stretching and node splitting. Graph interpretations of these methods are provided. Solving linear problems using enlarged matrices yields the answer to the originalAx=b problem.à can exhibit several desirable properties. For example,à can be constructed so that the valence of any row and/or column is smaller than some desired number (≥4). This is beneficial for algorithms that depend on the square of the number of entries of a row or column. Most particularly, matrix enlarging can results in a reduction of the fill-in in theR matrix which occurs during orthogonal factorization as a result of dense rows. Numerical experiments support these conjectures.     

On the eigenvalue distribution of a class of preconditioning methods

Summary A class of preconditioning methods depending on a relaxation parameter is presented for the solution of large linear systems of equationAx=b, whereA is a symmetric positive definite matrix. The methods are based on an incomplete factorization of the matrixA and include both pointwise and blockwise factorization. We study the dependence of the rate of convergence of the preconditioned conjugate gradient method on the distribution of eigenvalues ofC ^–1 A, whereC is the preconditioning matrix. We also show graphic representations of the eigenvalues and present numerical tests of the methods.     

Incomplete block factorization preconditioning for indefinite elliptic problems

Summary. The application of the finite difference method to approximate the solution of an indefinite elliptic problem produces a linear system whose coefficient matrix is block tridiagonal and symmetric indefinite. Such a linear system can be solved efficiently by a conjugate residual method, particularly when combined with a good preconditioner. We show that specific incomplete block factorization exists for the indefinite matrix if the mesh size is reasonably small, and that this factorization can serve as an efficient preconditioner. Some efforts are made to estimate the eigenvalues of the preconditioned matrix. Numerical results are also given. Received November 21, 1995 / Revised version received February 2, 1998 / Published online July 28, 1999     

The spectral properties of the preconditioned matrix for nonsymmetric saddle point problems

In this paper, on the basis of matrix splitting, two preconditioners are proposed and analyzed, for nonsymmetric saddle point problems. The spectral property of the preconditioned matrix is studied in detail. When the iteration parameter becomes small enough, the eigenvalues of the preconditioned matrices will gather into two clusters—one is near (0,0) and the other is near (2,0)—for the PPSS preconditioner no matter whether A is Hermitian or non-Hermitian and for the PHSS preconditioner when A is a Hermitian or real normal matrix. Numerical experiments are given, to illustrate the performances of the two preconditioners.     

On the effective computation of the inertia of a non-hermitian matrix

Summary given a complex lower Hessenberg matrixA with unit codiagonal, a hermitian matrixH is constructed such that, ifH is non-singular InA= InH. IfA is real,H is real symmetric. Classical results of Fujiwara on the root-separation problem and of Schwarz on the eigenvalue-separation problem are included as special cases.The authors' research was conducted at the Universidade Estadual de Campinas and supported by the Fundação de Amparo à Pesquisa do Estado de São Paulo, Brasil, under grant n⁰ 78/0490.     

Optimal preconditioners of a given sparsity pattern

We consider the problem of finding the symmetric positive definite preconditionerM of a given form- e.g., having nonzero elements only in specified positions — which minimizes the ratio of the largest to smallest eigenvalue ofM ^–1 A, for a given symmetric positive definitive matrixA. We show how this problem can be expressed as one of minimizing a convex function and how an optimization code can be used to solve the problem numerically. Results are presented showing optimal preconditioners of various sparsity patterns and comparing these to preconditioners that have been proposed in the literature. Several conjectures are made, based on results from the optimization code.This work was supported by the Advanced Research Projects Agency of the Department of Defense under contract F49620-87-C-0065 and by the Applied Mathematical Sciences Program of the U.S. Department of Energy under contract DE-AC02-76ER03077.     

An algorithm to compute a sparse basis of the null space

Summary LetA be a realm×n matrix with full row rankm. In many algorithms in engineering and science, such as the force method in structural analysis, the dual variable method for the Navier-Stokes equations or more generally null space methods in quadratic programming, it is necessary to compute a basis matrixB for the null space ofA. HereB isn×r, r=n–m, of rankr, withAB=0. In many instancesA is large and sparse and often banded. The purpose of this paper is to describe and test a variation of a method originally suggested by Topcu and called the turnback algorithm for computing a banded basis matrixB. Two implementations of the algorithm are given, one using Gaussian elimination and the other using orthogonal factorization by Givens rotations. The FORTRAN software was executed on an IBM 3081 computer with an FPS-164 attached array processor at the Triangle Universities Computing Center and on a CYBER 205 vector computer. Test results on a variety of structural analysis problems including two- and three-dimensional frames, plane stress, plate bending and mixed finite element problems are discussed. These results indicate that both implementations of the algorithm yielded a well-conditioned, banded, basis matrixB whenA is well-conditioned. However, the orthogonal implementation yielded a better conditionedB for large, illconditioned problems.The research by these authors was supported by the U.S. Air Force under grant No. AFOSR-83-0255 and by the National Science Foundation under grant No. MCS-82-19500The research by these authors was supported by the Applied Mathematical Sciences Program of the U.S. Department of Energy, under contract to Martin Marietta Energy Systems, Inc.     

Perturbation analysis of the orthogonal procrustes problem

Given two arbitrary real matricesA andB of the same size, the orthogonal Procrustes problem is to find an orthogonal matrixM such that the Frobenius norm MA – B is minimized. This paper treats the common case when the orthogonal matrixM is required to have a positive determinant. The stability of the problem is studied and supremum results for the perturbation bounds are derived.     

Circulant and skew-circulant preconditioners for skew-hermitian type Toeplitz systems

We study the solutions of Toeplitz systemsA _n x=b by the preconditioned conjugate gradient method. Then ×n matrixA _n is of the forma ₀ I+H _n wherea ₀ is a real number,I is the identity matrix andH _n is a skew-Hermitian Toeplitz matrix. Such matrices often appear in solving discretized hyperbolic differential equations. The preconditioners we considered here are the circulant matrixC _n and the skew-circulant matrixS _n whereA _n =1/2(C _n +S _n ). The convergence rate of the iterative method depends on the distribution of the singular values of the matricesC ^–1 _n A_n andS ^–1 _n A _n . For Toeplitz matricesA _n with entries which are Fourier coefficients of functions in the Wiener class, we show the invertibility ofC _n andS _n and prove that the singular values ofC ^–1 _n A _n andS ^–1 _n A _n are clustered around 1 for largen. Hence, if the conjugate gradient method is applied to solve the preconditioned systems, we expect fast convergence.     

Backward error analysis for linear systems associated with inverses ofH-matrices

In this paper, bounds on the growth factors resulting from Gaussian elimination applied to inverses ofH-matrices are developed and investigated. These bounds are then used in the error analysis for solving linear systemsAx =b whose coefficient matricesA are of this type. For each such system our results show that the Gaussian elimination without pivoting can proceed safely provided that the elements of the inverse of a certainM-matrix (associated with the coefficient matrixA) are not excessively large. We exhibit a particularly satisfactory situation for the special case whenA itself is an inverse of anM-matrix. Part of the first section of this paper is devoted to a discussion on some results of de Boor and Pinkus for the stability of triangular factorizations of systemsAx =b, whereA is a nonsingular totally nonnegative matrix, and to the explanation of why the analysis of de Boor and Pinkus is not applicable to the case when the coefficient matrixA is an inverse of anM-matrix.Research supported in part by NSF Grant MCS-8102114.Research supported in part by the U.S. Army Research Office under contract No. DAAG-29-81-K-0132 and in part by NSF Grant MCS-8219500.     

The symmetric procrustes problem

The following symmetric Procrustes problem arises in the determination of the strain matrix of an elastic structure: find the symmetric matrixX which minimises the Frobenius (or Euclidean) norm ofAX — B, whereA andB are given rectangular matrices. We use the singular value decomposition to analyse the problem and to derive a stable method for its solution. A perturbation result is derived and used to assess the stability of methods based on solving normal equations. Some comparisons with the standard, unconstrained least squares problem are given.