Randomized Kaczmarz solver for noisy linear systems

The Kaczmarz method is an iterative algorithm for solving systems of linear equations Ax=b. Theoretical convergence rates for this algorithm were largely unknown until recently when work was done on a randomized version of the algorithm. It was proved that for overdetermined systems, the randomized Kaczmarz method converges with expected exponential rate, independent of the number of equations in the system. Here we analyze the case where the system Ax=b is corrupted by noise, so we consider the system Ax≈b+r where r is an arbitrary error vector. We prove that in this noisy version, the randomized method reaches an error threshold dependent on the matrix A with the same rate as in the error-free case. We provide examples showing our results are sharp in the general context.     

Accelerated projection methods for computing pseudoinverse solutions of systems of linear equations

Iterative methods are developed for computing the Moore-Penrose pseudoinverse solution of a linear systemAx=b, whereA is anm ×n sparse matrix. The methods do not require the explicit formation ofA ^T A orAA ^T and therefore are advantageous to use when these matrices are much less sparse thanA itself. The methods are based on solving the two related systems (i)x=A ^T y,AA ^T y=b, and (ii)A ^T Ax=A ^T b. First it is shown how theSOR-andSSOR-methods for these two systems can be implemented efficiently. Further, the acceleration of theSSOR-method by Chebyshev semi-iteration and the conjugate gradient method is discussed. In particular it is shown that theSSOR-cg method for (i) and (ii) can be implemented in such a way that each step requires only two sweeps through successive rows and columns ofA respectively. In the general rank deficient and inconsistent case it is shown how the pseudoinverse solution can be computed by a two step procedure. Some possible applications are mentioned and numerical results are given for some problems from picture reconstruction.     

A Multilevel AINV Preconditioner

In this paper we describe an algebraic multilevel extension of the approximate inverse AINV preconditioner for solving symmetric positive definite linear systems Ax=b with the preconditioned conjugate gradient method. The smoother is the approximate inverse M and the coarse grids and the interpolation operator are constructed by looking at the entries of M. Numerical examples are given for problems arising from discretization of partial differential equations.     

How residual bounds for restarted GMRES describe the real behaviour

Given a system of linear equations Ax=b with a nonsingular matrix A, new bounds for the norm of the residual vector are presented and shown by example. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Convergence theorems for block splitting iterative methods for linear systems

In this paper, we will present the block splitting iterative methods with general weighting matrices for solving linear systems of algebraic equations Ax=b$\mathit{Ax}=b$ when the coefficient matrix A is symmetric positive definite of block form, and establish the convergence theories with respect to the general weighting matrices but special splittings. Finally, a numerical example shows the advantage of this method.     

On extensions of Cramer's rule for solutions of restricted linear systems 1

For the unique solution of a special consistent restricted linear system Ax=bx?M we derive two different determinantal forms, which both reduce to Cramer's classical rule if A is nonsingular. The representation for the minimum Euclidean-norm solution of Ax=b given recently by Ben-Israel [2] is reobtained in our first approach as special case. A determinantal formula for any {1, 2}-generalized inverse is also given.     

The continued fraction methods for the solution of systems of linear equations

A class of iterative methods is presented for the solution of systems of linear equationsAx=b, whereA is a generalm ×n matrix. The methods are based on a development as a continued fraction of the inner product (r, r), wherer=b-Ax is the residual. The methods as defined are quite general and include some wellknown methods such as the minimal residual conjugate gradient method with one step.     

Efficient preconditioning of linear systems arising from the discretization of hyperbolic conservation laws

In this paper, we describe a novel formulation of a preconditioned BiCGSTAB algorithm for the solution of ill-conditioned linear systems Ax=b. The developed extension enables the control of the residual r _m =b–Ax _m of the approximate solution x _m independent of the specific left, right or two-sided preconditioning technique considered. Thereby, the presented modification does not require any additional computational effort and can be introduced directly into existing computer codes. Furthermore, the proceeding is not restricted to the BiCGSTAB method, hence the strategy can serve as a guideline to extend similar Krylov sub-space methods in the same manner. Based on the presented algorithm, we study the behavior of different preconditioning techniques. We introduce a new physically motivated approach within an implicit finite volume scheme for the system of the Euler equations of gas dynamics which is a typical representative of hyperbolic conservation laws. Thereupon a great variety of realistic flow problems are considered in order to give reliable statements concerning the efficiency and performance of modern preconditioning techniques.     

On the solution of a nonlinear programming problem by decomposition

A decomposition method, used in least-weight plastic design, is extended to solve problems with nonlinearity arising from variable structure geometry. The problem considered is that of finding vectorsx ₁,x ₂, andq that minimize [l max{|x ₁|, |x ₂|}], subject toAx ₁=b ₁ andAx ₂=b ₂, where both the vectorl and the matrixA are nonlinear functions ofq.     

Incomplete factorizations of singular linear systems

Recent works have shown that, whenA is a Stieltjes matrix, its so-called modified incomplete factorizations provide effective preconditioning matrices for solvingAx=b by polynomially accelerated iterative methods. We extend here these results to the singular case with the conclusion that the latter techniques are able to solve singular systems at the same speed as regular systems.Research supported by the Fonds National de la Recherche Scientifique (Belgium) — Aspirant.     

Equivalence of Linear Programming and Basis Pursuit

In this note, we show that linear programming and the prominent Basis Pursuit problem (i.e., minimizing the ℓ₁-norm of a vector x subject to an underdetermined linear equation system Ax = b) are theoretically equivalent, and briefly discuss possible ramifications regarding computational complexity and practical applicability. (© 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Simultaneous backward stability of Gauss and Gauss–Jordan elimination

It is well known that some pivoting strategies are backward stable for Gauss elimination but are not backward stable for Gauss–Jordan elimination (GJE) when these procedures are used to solve a linear system Ax=b. We analyse the simultaneous backward stability for Gauss and GJE of several pivoting strategies, including a pivoting strategy which we call double partial pivoting. Copyright © 2002 John Wiley & Sons, Ltd.     

Convergence analysis for a multiplicatively relaxed EM algorithm

The expectation maximization (EM) algorithm is an iterative procedure used to determine maximum likelihood estimators in situations of incomplete data. In the case of independent Poisson variables it converges to a solution of a problem of the form min ∑[〈aⁱ,x〉 ? b_i log 〈aⁱ, x〉] such that x ?0. Thus, it can be used to solve systems of the form Ax = b, x?0 (with A stochastic and b positive.) It converges to a feasible solution if it exists and to an approximate one otherwise (the one that minimizes d (b, Ax), where d is the Kullback–Leibler information divergence). We study the convergence of the multiplicatively relaxed version proposed by Tanaka for use in positron emission tomography. We prove global convergence in the underrelaxed and unrelaxed cases. In the overrelaxed case we present a local convergence theorem together with two partial global results: the sequence generated by the algorithm is bounded and, if it converges, its limit point is a solution of the aforementioned problem.     

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.     

A note on preconditioned GMRES for solving singular linear systems

For solving a singular linear system Ax=b by GMRES, it is shown in the literature that if A is range-symmetric, then GMRES converges safely to a solution. In this paper we consider preconditioned GMRES for solving a singular linear system, we construct preconditioners by so-called proper splittings, which can ensure that the coefficient matrix of the preconditioned system is range-symmetric.     

Solving linear programs with finite precision: I. Condition numbers and random programs

We define a condition number (A,b,c) for a linear program min x s.t. Ax=b,x0 and give two characterizations via distances to degeneracy and singularity. We also give bounds for the expected value, as well as for higher moments, of log (A,b,c) when the entries of A,b and c are i.i.d. random variables with normal distribution. This work has been substantially funded by a grant from the Research Grants Council of the Hong Kong SAR (project number CityU 1085/02P)     

Krylov sequences of maximal length and convergence of GMRES

In most practical cases, the convergence of the GMRES method applied to a linear algebraic systemAx=b is determined by the distribution of eigenvalues ofA. In theory, however, the information about the eigenvalues alone is not sufficient for determining the convergence. In this paper the previous work of Greenbaum et al. is extended in the following direction. It is given a complete parametrization of the set of all pairs {A, b} for which GMRES(A, b) generates the prescribed convergence curve while the matrixA has the prescribed eigenvalues. Moreover, a characterization of the right hand sidesb for which the GMRES(A, b) converges exactly inm steps, wherem is the degree of the minimal polynomial ofA, is given. This work was supported by AS CR Grant A2030706. Part of the work was performed while the third author visited Instituto di Analisi Numerica (IAN CNR).     

Trust-region and other regularisations of linear least-squares problems

We consider methods for regularising the least-squares solution of the linear system Ax=b. In particular, we propose iterative methods for solving large problems in which a trust-region bound ‖x‖≤Δ is imposed on the size of the solution, and in which the least value of linear combinations of ‖Ax−b‖₂ ^q and a regularisation term ‖x‖₂ ^p for various p and q=1,2 is sought. In each case, one or more “secular” equations are derived, and fast Newton-like solution procedures are suggested. The resulting algorithms are available as part of the ALAHAD optimization library. This work was partially supported by EPSRC grants EP/E053351/1 and EP/F005369/1.     

Treatment of near-breakdown in the CGS algorithm

Lanczos' method for solving the system of linear equationsAx=b consists in constructing a sequence of vectors (x _k ) such thatr _k =b–Ax _k =P _k (A)r ₀ wherer ₀=b–Ax ₀.P _k is an orthogonal polynomial which is computed recursively. The conjugate gradient squared algorithm (CGS) consists in takingr _k =P _k ² (A)r₀. In the recurrence relation forP _k , the coefficients are given as ratios of scalar products. When a scalar product in a denominator is zero, then a breakdown occurs in the algorithm. When such a scalar product is close to zero, then rounding errors can seriously affect the algorithm, a situation known as near-breakdown. In this paper it is shown how to avoid near-breakdown in the CGS algorithm in order to obtain a more stable method.     

Perturbation theory for pseudo-inverses

A perturbation theory for pseudo-inverses is developed. The theory is based on a useful decomposition (theorem 2.1) ofB ⁺ -A ⁺ whereB andA arem ×n matrices. Sharp estimates of B ⁺ -A ⁺ are derived for unitary invariant norms whenA andB are of the same rank and B -A is small. Under similar conditions the perturbation of a linear systemAx=b is studied. Realistic bounds on the perturbation ofx=A ⁺ b andr=b=Ax are given. Finally it is seen thatA ⁺ andB ⁺ can be compared if and only ifR(A) andR(B) as well asR(A ^H ) andR(B ^H ) are in the acute case. Some theorems valid only in the acute case are also proved.This work was sponsored in part by The Swedish Institute of Applied Mathematics.