Solving large and sparse linear least-squares problems by conjugate gradient algorithms

Abstract

Let A ε ℛm × n(with m ⩾ n and rank (A) = n) and b ε ℛm × 1 be given. Assume that an approximation to x = A†b = (ATA)−1ATb is to be calculated. This problem is transformed into an equivalent problem with a symmetric and positive definite matrix C = D−1 (RT)−1ATAR−1D−1 where D is a diagonal matrix and R is an upper triangular matrix. The conjugate gradient (CG) algorithm is applied in the solution of the system of linear algebraic equations Cy = d (y = DRx, d = D−1(RT)−1ATb). If D = R = I, then the CG algorithm is in fact applied to the system of normal equations ATAx = ATb and the speed of convergence could be very slow. If D and R are obtained by some kind of orthogonalization DR = QTA with Q εℛm × n satisfying QTQ = I; D is often equal to I and Q is never used in the CG algorithm) and if the calculations are performed without rounding errors, then C = I and the CG algorithm converges in one iteration only. Even if the orthogonalization process is carried out with rounding errors, the matrix C is normally close to the identity matrix I and the CG algorithm is quickly convergent. However, for large m and n the orthogonal decomposition is an expensive process (both in regard to storage and in regard to computing time). Therefore it may be profitable to calculate an incomplete orthogonal decomposition. This is achieved by introducing a special parameter T, a drop-tolerance, such that all elements which in the course of the computations become smaller in absolute value than T are removed. Numerical examples are given to illustrate that the CG algorithm applied to system Cy = d and used with incomplete factors D and R is very efficient for some classes of problems. It should be emphasized that matrix C is never calculated explicitly; the whole work is carried out by the use of A, D and R only.