Sparse Linear Least-squares Problems Research Articles

On greedy randomized block Gauss–Seidel method with averaging for sparse linear least-squares problems

On greedy randomized block Gauss–Seidel method with averaging for sparse linear least-squares problems

A Robust Algebraic Domain Decomposition Preconditioner for Sparse Normal Equations

Solving the normal equations corresponding to large sparse linear least-squares problems is an important and challenging problem. For very large problems, an iterative solver is needed and, in general, a preconditioner is required to achieve good convergence. In recent years, a number of preconditioners have been proposed. These are largely serial, and reported results demonstrate that none of the commonly used preconditioners for the normal equations matrix is capable of solving all sparse least-squares problems. Our interest is thus in designing new preconditioners for the normal equations that are efficient and robust and can be implemented in parallel. Our proposed preconditioners can be constructed efficiently and algebraically without any knowledge of the problem and without any assumption on the least-squares matrix except that it is sparse. We exploit the structure of the symmetric positive definite normal equations matrix and use the concept of algebraic local symmetric positive semidefinite splittings to introduce two-level Schwarz preconditioners for least-squares problems. The condition number of the preconditioned normal equations matrix is shown to be theoretically bounded independently of the number of subdomains in the splitting. This upper bound can be adjusted using a single parameter $\tau$ that the user can specify. We discuss how the new preconditioners can be implemented on top of the PETSc library using only 150 lines of Fortran, C, or Python code. Problems arising from practical applications are used to compare the performance of the proposed new preconditioner with that of other preconditioners.

A Schur complement approach to preconditioning sparse linear least-squares problems with some dense rows

The effectiveness of sparse matrix techniques for directly solving large-scale linear least-squares problems is severely limited if the system matrix A has one or more nearly dense rows. In this paper, we partition the rows of A into sparse rows and dense rows (As and Ad) and apply the Schur complement approach. A potential difficulty is that the reduced normal matrix AsTAs is often rank-deficient, even if A is of full rank. To overcome this, we propose explicitly removing null columns of As and then employing a regularization parameter and using the resulting Cholesky factors as a preconditioner for an iterative solver applied to the symmetric indefinite reduced augmented system. We consider complete factorizations as well as incomplete Cholesky factorizations of the shifted reduced normal matrix. Numerical experiments are performed on a range of large least-squares problems arising from practical applications. These demonstrate the effectiveness of the proposed approach when combined with either a sparse parallel direct solver or a robust incomplete Cholesky factorization algorithm.

The State-of-the-Art of Preconditioners for Sparse Linear Least-Squares Problems

In recent years, a variety of preconditioners have been proposed for use in solving large sparse linear least-squares problems. These include simple diagonal preconditioning, preconditioners based on incomplete factorizations, and stationary inner iterations used with Krylov subspace methods. In this study, we briefly review preconditioners for which software has been made available, then present a numerical evaluation of them using performance profiles and a large set of problems arising from practical applications. Comparisons are made with state-of-the-art sparse direct methods.

On Using Cholesky-Based Factorizations and Regularization for Solving Rank-Deficient Sparse Linear Least-Squares Problems

By examining the performance of modern parallel sparse direct solvers and exploiting our knowledge of the algorithms behind them, we perform numerical experiments to study how they can be used to efficiently solve rank-deficient sparse linear least-squares problems arising from practical applications. The Cholesky factorization of the normal equations breaks down when the least-squares problem is rank-deficient, while applying a symmetric indefinite solver to the augmented system can give an unacceptable level of fill in the factors. To try to resolve these difficulties, we consider a regularization procedure that modifies the diagonal of the unregularized matrix. This leads to matrices that are easier to factorize. We consider both the regularized normal equations and the regularized augmented system. We employ the computed factors of the regularized systems as preconditioners with an iterative solver to obtain the solution of the original (unregularized) problem. Furthermore, we look at using limited-me...

Using Perturbed $QR$ Factorizations to Solve Linear Least-Squares Problems

We propose and analyze a new tool to help solve sparse linear least-squares problems $\min_{x}\|Ax-b\|_{2}$. Our method is based on a sparse $QR$ factorization of a low-rank perturbation $\hat{A}$ of $A$. More precisely, we show that the $R$ factor of $\hat{A}$ is an effective preconditioner for the least-squares problem $\min_{x}\|Ax-b\|_{2}$, when solved using LSQR. We propose applications for the new technique. When $A$ is rank deficient, we can add rows to ensure that the preconditioner is well conditioned without column pivoting. When $A$ is sparse except for a few dense rows, we can drop these dense rows from $A$ to obtain $\hat{A}$. Another application is solving an updated or downdated problem. If $R$ is a good preconditioner for the original problem $A$, it is a good preconditioner for the updated/downdated problem $\hat{A}$. We can also solve what-if scenarios, where we want to find the solution if a column of the original matrix is changed/removed. We present a spectral theory that analyzes the generalized spectrum of the pencil $(A^{*}A,R^{*}R)$ and analyze the applications.

Preconditioned GMRES methods with incomplete Givens orthogonalization method for large sparse least-squares problems

We propose to precondition the GMRES method by using the incomplete Givens orthogonalization (IGO) method for the solution of large sparse linear least-squares problems. Theoretical analysis shows that the preconditioner satisfies the sufficient condition that can guarantee that the preconditioned GMRES method will never break down and always give the least-squares solution of the original problem. Numerical experiments further confirm that the new preconditioner is efficient. We also find that the IGO preconditioned BA-GMRES method is superior to the corresponding CGLS method for ill-conditioned and singular least-squares problems.

Solving sparse linear least-squares problems on some supercomputers by using large dense blocks

Efficient subroutines for dense matrix computations have recently been developed and are available on many high-speed computers. On some computers the speed of many dense matrix operations is near to the peak-performance. For sparse matrices storage and operations can be saved by operating only and storing only nonzero elements. However, the price is a great degradation of the speed of computations on supercomputers (due to the use of indirect addresses, to the need to insert new nonzeros in the sparse storage scheme, to the lack of data locality, etc.).

Sparse QR factorization in MATLAB

In the recently presented sparse matrix extension of MATLAB, there is no routine for sparse QR factorization. Sparse linear least-squares problems are instead solved by the augmented system method. The accuracy in computed solutions is strongly dependent on a scaling parameter δ. Its optimal value is expensive to compute, and it must therefore be approximated by a simple heuristic. We describe a multifrontal method for sparse QR factorization and its implementation in MATLAB. It is well known that the multifrontal approach is suitable for vector machines. We show that it is also attractive in MATLAB. In both cases, scalar operations are expensive, and the reformulation of the sparse problem into dense subproblems is advantageous. Using the new routine, we implement two methods for the solution of sparse linear least-squares problems and compare these with the built-in MATLAB function. We show that the QR-based methods normally are much faster and more accurate than the MATLAB implementation of the augmented system method. A better choice of the parameter δ or iterative refinement must be used to make the augmented system method as accurate as the methods based on QR factorization.

SIRT- and CG-type methods for the iterative solution of sparse linear least-squares problems

SIRT and CG-type methods have been successfully employed for the approximate solution of least-squares problems arising in tomography. In this paper we study and compare their convergence and regularization properties. It is pointed out that SIRT methods apply an uncontrollable implicit rescaling which affects the statistical characteristics of the system, whereas CG-type methods do not. For a large class of model problems it is shown that virtually the same solutions as obtained by SIRT methods can be obtained by applying a CG-type method to a properly rescaled system, but with an amount of work proportional to the square root of the amount of work with SIRT.

On the augmented system approach to sparse least-squares problems

We study the augmented system approach for the solution of sparse linear least-squares problems. It is well known that this method has better numerical properties than the method based on the normal equations. We use recent work by Arioli et al. (1988) to introduce error bounds and estimates for the components of the solution of the augmented system. In particular, we find that, using iterative refinement, we obtain a very robust algorithm and our estimates of the error are accurate and cheap to compute. The final error and all our error estimates are much better than the classical or Skeel's error analysis (1979) indicates. Moreover, we prove that our error estimates are independent of the row scaling of the augmented system and we analyze the influence of the Bjorck scaling (1967) on these estimates. We illustrate this with runs both on large-scale practical problems and contrived examples, comparing the numerical behaviour of the augmented systems approach with a code using the normal equations. These experiments show that while the augmented system approach with iterative refinement can sometimes be less efficient than the normal equations approach, it is comparable or better when the least-squares matrix has a full row, and is, in any case, much more stable and robust.

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

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.

A Comparison of Some Methods for Solving Sparse Linear Least-Squares Problems

The method of normal equations, the Peters–Wilkinson algorithm and an algorithm based on Givens rotations for solving large sparse linear least squares problems are discussed and compared. Numerical experiments show that the method of normal equations should be considered when the observation matrix is sparse and well conditioned. For ill-conditioned problems, the algorithm based on Givens rotations is preferable.