- Home
- Documents
*The State-of-the-Art of Preconditioners for Sparse Linear ... iterative methods. In Section 4, we...*

prev

next

out of 35

View

0Download

0

Embed Size (px)

36

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

NICHOLAS GOULD and JENNIFER SCOTT, Rutherford Appleton Laboratory

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.

Categories and Subject Descriptors: G.1.3 [Numerical Linear Algebra]: Linear Systems (Direct and Iter- ative Methods); G.1.6 [Optimization]: Least-squares Methods

General Terms: Algorithms, Performance

Additional Key Words and Phrases: Least-squares problems, normal equations, augmented system, sparse matrices, direct solvers, iterative solvers, preconditioning

ACM Reference Format: Nicholas Gould and Jennifer Scott. 2017. The state-of-the-art of preconditioners for sparse linear least- squares problems. ACM Trans. Math. Softw. 43, 4, Article 36 (January 2017), 35 pages. DOI: http://dx.doi.org/10.1145/3014057

1. INTRODUCTION

The method of least-squares is a commonly used approach to find an approximate solution of overdetermined or inexactly specified systems of equations. Since the de- velopment of the principle of least-squares by Gauss in 1795 [Farebrother 1999], the solution of least-squares problems has been, and continues to be, a fundamental method in scientific data fitting. Least-squares solvers are used across a wide range of disci- plines, for everything from simple curve fitting through the estimation of satellite image sensor characteristics, data assimilation for weather forecasting and for climate modelling, to powering Internet mapping services, exploration seismology, NMR spec- troscopy, piezoelectric crystal identification (used in ultrasound for medical imaging), aerospace systems, and neural networks.

In this study, we are interested in the important special case of the linear least- squares problem,

min x

‖b − Ax‖2, (1)

This work was supported by EPSRC grants EP/I013067/1 and EP/M025179/1. Authors’ addresses: N. Gould and J. Scott, Scientific Computing Department, STFC Rutherford Appleton Laboratory, Harwell Campus, Oxfordshire, OX11 0QX, UK; emails: nick.gould@stfc.ac.uk, jennifer.scott@ stfc.ac.uk. Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies show this notice on the first page or initial screen of a display along with the full citation. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers, to redistribute to lists, or to use any component of this work in other works requires prior specific permission and/or a fee. Permissions may be requested from Publications Dept., ACM, Inc., 2 Penn Plaza, Suite 701, New York, NY 10121-0701 USA, fax +1 (212) 869-0481, or permissions@acm.org. c© 2017 ACM 0098-3500/2017/01-ART36 $15.00 DOI: http://dx.doi.org/10.1145/3014057

ACM Transactions on Mathematical Software, Vol. 43, No. 4, Article 36, Publication date: January 2017.

http://dx.doi.org/10.1145/3014057 http://dx.doi.org/10.1145/3014057

36:2 N. Gould and J. Scott

where A ∈ IRm×n with m ≥ n is large and sparse, and b ∈ IRm. Solving Equation (1) is mathematically equivalent to solving the n × n normal equations

Cx = AT b, C = AT A; (2) this, in turn, is equivalent to solving the (m+ n) × (m+ n) augmented system

Ky = c, K = [

Im A AT 0

] , y =

[ r(x)

x

] , c =

[ b 0

] , (3)

where r(x) = b − Ax is the residual vector and Im is the m × m identity matrix. In- creasingly, the sizes of the problems that scientists and engineers wish to solve are getting larger (problems in many millions of variables are becoming typical); they are also often ill-conditioned. In other applications, it is necessary to solve many thou- sands of problems of modest size; thus, efficiency in this case is essential. The normal equations are attractive in that they are always consistent and positive semidefinite (positive definite if A is of full column rank). However, a well-known drawback is that the condition number of C is the square of the condition number of Aso that the normal equations are often highly ill-conditioned [Björck 1996]. Furthermore, the density of C can be much greater than that of A (if A has a single dense row, C will be dense). The main disadvantages of working with the augmented system are that K is symmetric indefinite and is much larger than C (particularly if m � n).

Two main classes of methods may be used in attempting to solve these linear systems: direct methods and iterative methods. A direct method proceeds by computing an ex- plicit factorization, either a sparse LLT Cholesky factorization of the normal equations (2) (assuming that A is of full column rank so that C is positive definite) or a sparse LDLT factorization of the augmented system (3). Alternatively, a QR factorization of A may be used, that is, a “thin” QR factorization of the form

A = Q [

R 0

] ,

where Q is an m × m orthogonal matrix and R is an n × n upper triangular matrix. While direct solvers are generally highly reliable, iterative methods may be preferred because they often require significantly less storage and, in some applications, it may not be necessary to solve the system with the high accuracy offered by a direct solver. However, the successful application of an iterative method often requires a suitable preconditioner to achieve acceptable (and ideally, fast) convergence rates. Currently, there is much less knowledge of preconditioners for least-squares problems than there is for sparse symmetric linear systems and, as remarked in Bru et al. [2014], “the problem of robust and efficient iterative solution of least-squares problems is much harder than the iterative solution of systems of linear equations.” This is, at least in part, because A does not have the properties of differential problems that can make standard preconditioners effective.

In the past decade or so, a number of different techniques for preconditioning Krylov subspace methods for least-squares problems have been developed. A brief overview with a comprehensive list of references is included in the introduction to the recent paper by Bru et al. [2014]. However, in the literature, the reported experiments on the performance of different preconditioners are often limited to a small set of problems, generally arising from a specific application. Moreover, they may use prototype codes that are not available for others to test and they may only be run using MATLAB. Our aim is to perform a wider study in which we use a large set of test problems to evaluate the performance of a range of preconditioners for which software has been made available. The intention is to gain a clearer understanding of when particular

ACM Transactions on Mathematical Software, Vol. 43, No. 4, Article 36, Publication date: January 2017.

Preconditioners for Sparse Linear Least-Squares Problems 36:3

Table I. Test Machine Characteristics

Processor 8× Intel i7-4790 (3.6 GHz) Memory 15.6 GB Compiler gfortran version 4.7 with option -O BLAS open BLAS (serial) or Intel MKL (serial vs. parallel)

preconditioners perform well (or perform poorly); we will use this to influence our future work on linear least-squares. Our attention is limited to preconditioners for which software in Fortran or C is available. It is beyond the scope of this work to provide efficient and robust implementations for all the approaches that have appeared in the literature (although, even then, as we discuss in Section 8, we have found it necessary in some cases to modify and possibly reengineer some of the existing software to make it suitable for use in this study).

The rest of the article is organized as follows. In Section 2, we describe our test envi- ronment, including the set of problems used in this study. Direct solvers for solving the normal equations and/or the augmented system are briefly recalled in Section 3. One of these (HSL_MA97) is used for comparison with the performance of the preconditioned iterative methods. In Section 4, we report on experiments with two methods, LSQR and LSMR, that are mathematically equivalent to applying conjugate gradients and MINRES, respectively, to the normal equations but have favorable numerical proper- ties. On the basis of our findings, LSMR is used in the rest of our experiments. Precon- ditioning strategies are briefly described in Sections 5 to 7. The software used in our experiments is discussed in Section 8. We discuss the benefits of simple parallelism in Section 9. We present numerical results in Section 10 and our conclusions in Section 11.

2. TEST ENVIRONMENT

The characteristics of the machine used to perform our tests are given in Table I. In our experiments, the direct solvers (see Section 3) are run in parallel, using four processors. Our initial experiments on iterative methods (those in Sections 4 and 8) are run on a single processor, although, where Basic Linear Algebra Subprograms (BLASs) are used, these may take advantage of parallel