- Home
- Documents
*LSMR: AN ITERATIVE ALGORITHM FOR SPARSE LSMR: AN ITERATIVE ALGORITHM FOR SPARSE LEAST-SQUARES...*

prev

next

out of 23

View

4Download

0

Embed Size (px)

LSMR: AN ITERATIVE ALGORITHM FOR SPARSE LEAST-SQUARES PROBLEMS∗

DAVID CHIN-LUNG FONG† AND MICHAEL SAUNDERS‡

Abstract. An iterative method LSMR is presented for solving linear systems Ax = b and least- squares problem min ‖Ax− b‖2, with A being sparse or a fast linear operator. LSMR is based on the Golub-Kahan bidiagonalization process. It is analytically equivalent to the MINRES method applied to the normal equation ATAx = ATb, so that the quantities ‖ATrk‖ are monotonically decreasing (where rk = b − Axk is the residual for the current iterate xk). In practice we observe that ‖rk‖ also decreases monotonically. Compared to LSQR, for which only ‖rk‖ is monotonic, it is safer to terminate LSMR early. Improvements for the new iterative method in the presence of extra available memory are also explored.

Key words. least-squares problem, sparse matrix, LSQR, MINRES, Krylov subspace method, Golub-Kahan process, conjugate-gradient method, minimum-residual method, iterative method

AMS subject classifications. 15A06, 65F10, 65F20, 65F22, 65F25, 65F35, 65F50, 93E24

DOI. xxx/xxxxxxxxx

1. Introduction. We present a numerical method called LSMR for computing a solution x to the following problems:

Unsymmetric equations: solve Ax = b Linear least squares: minimize ‖Ax− b‖2 Regularized least squares: minimize

∥∥∥∥(AλI ) x−

( b 0

)∥∥∥∥ 2

where A ∈ Rm×n, b ∈ Rm, and λ ≥ 0. The matrix A is used as an operator for which products of the form Av and ATu can be computed for various v and u. Thus A is normally large and sparse and need not be explicitly stored.

LSMR is similar in style to the well known method LSQR [11, 12] in being based on the Golub-Kahan bidiagonalization of A [4]. LSQR is equivalent to the conjugate- gradient (CG) method applied to the normal equation (ATA+λ2I)x = ATb. It has the property of reducing ‖rk‖ monotonically, where rk = b − Axk is the residual for the approximate solution xk. (For simplicity, we are letting λ = 0.) In contrast, LSMR is equivalent to MINRES [10] applied to the normal equation, so that the quantities ‖ATrk‖ are monotonically decreasing. In practice we observe that ‖rk‖ also decreases monotonically, and is never very far behind the corresponding value for LSQR. Hence, although LSQR and LSMR ultimately converge to similar points, it is safer to use LSMR in situations where the solver must be terminated early.

Stopping conditions are typically based on backward error : the norm of some per- turbation to A for which the current iterate xk solves the perturbed problem exactly. Experiments on many sparse least-squares test problems show that for LSMR, a cer- tain cheaply computable backward error for each xk is close to the optimal (smallest possible) backward error. This is an unexpected but highly desirable advantage.

∗Version of May 28, 2010. Technical Report SOL 2010-2 http://www.siam.org/journals/sisc/ for Copper Mountain Special Issue 2010 †iCME, Stanford University (clfong@stanford.edu). Partially supported by a Stanford Graduate

Fellowship. ‡Systems Optimization Laboratory, Department of Management Science and Engineering, Stan-

ford University, CA 94305-4026 (saunders@stanford.edu). Partially supported by Office of Naval Research grant N00014-08-1-0191 and by the U.S. Army Research Laboratory, through the Army High Performance Computing Research Center, Cooperative Agreement W911NF-07-0027.

1

2 D. C.-L. FONG AND M. A. SAUNDERS

1.1. Notation. Matrices are denoted by A,B, . . . , vectors by v, w, . . . , and scalars by α, β, · · · . Two exceptions are c and s, which denote the significant compo- nents of a plane rotation matrix, with c2 + s2 = 1. For a vector v, ‖v‖ always denotes the 2-norm of v. For a matrix A, ‖A‖ usually denotes the Frobenius norm, and the condition number of a matrix A is defined by cond(A) = ‖A‖‖A+‖, where A+ denotes the pseudoinverse of A.

2. Derivation of LSMR.

2.1. Golub-Kahan process. We begin with the Golub-Kahan process [4]. 1. Set β1u1 = b (shorthand for β1 = ‖b‖, u1 = b/β1) and α1v1 = ATu1. 2. For k = 1, 2, . . . , set

βk+1uk+1 = Avk − αkuk and αk+1vk+1 = ATuk+1 − βk+1vk. (2.1)

After k steps, we have

AVk = Uk+1Bk and ATUk+1 = Vk+1LTk+1,

where we define

Lk =

α1 β2 α2

. . . . . . βk αk

, Bk = α1 β2 α2

. . . . . . βk αk

βk+1

= (

Lk βk+1e

T k

) .

Now consider

ATAVk = ATUk+1Bk = (Vk+1LTk+1)Bk

= Vk+1

( LTk βk+1ek 0 αk+1

)( Lk

βk+1e T k

) = Vk+1

( LTk Lk + β

2 k+1eke

T k

αk+1βk+1e T k

) = Vk+1

( BTkBk

αk+1βk+1e T k

) .

This is equivalent to what would be generated by the symmetric Lanczos process with matrix ATA and starting vector ATb.

2.2. Using Golub-Kahan to solve the normal equation. Krylov subspace methods for solving linear equations form solution estimates xk = Vkyk for some yk, where (as above) the columns of Vk are an expanding set of theoretically orthonormal vectors.

For the equation ATAx = ATb, any solution x has the property of minimizing ‖r‖, where r = b − Ax is the corresponding residual vector. Thus, in the development of LSQR it was natural to choose yk to minimize ‖rk‖ at each stage. Since

rk = b−AVkyk = β1u1 − Uk+1Bkyk = Uk+1(β1e1 −Bkyk),

where Uk+1 is theoretically orthonormal, the subproblem minyk ‖β1e1 −Bkyk‖ easily arose. In contrast, for LSMR we wish to minimize ‖ATrk‖. Let β̄k = αkβk for all k.

LSMR: AN ITERATIVE ALGORITHM FOR LEAST-SQUARES 3

Since

ATrk = ATb−ATAxk = β1ATu1 −ATAxk = β1α1v1 −ATAVkyk

= β̄1v1 − Vk+1 (

BTkBk αk+1βk+1e

T k

) yk

= Vk+1

( β̄1e1 −

( BTkBk β̄k+1e

T k

) yk

) ,

we are led to the subproblem

min yk ‖AT rk‖ = min

yk

∥∥∥∥β̄1e1 − ( BTkBkβ̄k+1eTk ) yk

∥∥∥∥ . (2.2) Efficient solution of this least-squares subproblem is the heart of algorithm LSMR.

2.3. Two QR factorizations. As in LSQR, we form the QR factorization

Qk+1Bk = ( Rk 0

) , Rk =

ρ1 θ2

ρ2 . . . . . . θk

ρk

. (2.3)

If we define tk = Rkyk and solve RTkqk = β̄k+1ek, we have qk = (β̄k+1/ρk)ek = ϕkek with ρk = (Rk)kk and ϕk = β̄k+1/ρk. Then we perform a second QR factorization

Q̄k+1

( RTk β̄1e1 ϕke

T k 0

) = ( R̄k zk 0 ζ̄k+1

) , R̄k =

ρ̄1 θ̄2

ρ̄2 . . . . . . θ̄k

ρ̄k

. (2.4)

Combining what we have gives

min yk

∥∥∥∥β̄1e1 − ( BTkBkβ̄k+1eTk ) yk

∥∥∥∥ = minyk ∥∥∥∥β̄1e1 − (RTkRkqTk Rk

) yk

∥∥∥∥ = min

tk

∥∥∥∥β̄1e1 − ( RTkϕkeTk ) tk

∥∥∥∥ = min

tk

∥∥∥∥( zkζ̄k+1 ) − ( R̄k 0

) tk

∥∥∥∥ . (2.5) The subproblem is solved by choosing tk from R̄ktk = zk.

2.4. Recurrence for xk. Let Wk and W̄k be computed by forward substitution from RTkW

T k = V

T k and R̄

T k W̄

T k = W

T k . Then from xk = Vkyk, Rkyk = tk, and

R̄ktk = zk, we have

xk = WkRkyk = Wktk = W̄kR̄ktk = W̄kzk = xk−1 + ζkw̄k.

4 D. C.-L. FONG AND M. A. SAUNDERS

2.5. Recurrence for Wk and W̄k. If we write

Vk = ( v1 v2 · · · vk

) , Wk =

( w1 w2 · · · wk

) ,

W̄k = ( w̄1 w̄2 · · · w̄k

) zk =

( ζ1 ζ2 · · · ζk

)T ,

an important fact is that when k increases to k + 1, all quantities remain the same except for one additional term.

The first QR factorization proceeds as follows. At iteration k, we write

Ql,l+1 =

Il−1

cl sl −sl cl

Ik−l−1

. Now if Qk+1 = Qk,k+1 . . . Q3,2Q1,2, we have

Qk+1Bk+1 = Qk

( Bk αk+1ek+1

βk+2

) =

Rk θk+1ek0 ᾱk+1 βk+2

Qk+2Bk+1 = Qk+1,k+2

Rk θk+1ek0 ᾱk+1 βk+2

= Rk θk+1ek0 ρk+1

0 0

and we see that θk+1 = skαk+1 = (βk+1/ρk)αk+1 = β̄k+1/ρk = ϕk. Therefore we can now write θk+1 instead of ϕk.

For the second QR factorization, if Q̄k+1 = Q̄k,k+1 . . . Q̄3,2Q̄1,2, we know that

Q̄k+1

( RTk

θk+1e T k

) = ( R̄k 0

) .

Therefore we would have

Q̄k+2

( RTk+1

θk+2e T k+1

) = Q̄k+1,k+2

R̄k θ̄k+1ekc̄kρk+1 θk+2

= R̄k θ̄k+1ekρ̄k+1

0

. By considering the last row of the matrix equation RTk+1W

T k+1 = V

T k+1 we obtain

θk+1w T k + ρk+1w

T k+1 = v

T k+1,

and from the last row of the matrix equation R̄Tk+1W̄ T k+1 = W

T k+1 we obtain

θ̄k+1w̄ T k + ρ̄k+1w̄

T k+1 = w

T k+1.

These equations serve to define wk+1 and w̄k+1.

2.6. The two rotations. To summarize, the rotations Qk,k+1 and Q̄k,k+1 have the following effects on our computation:(

ck sk −sk ck

)( ᾱk βk+1 αk+1

) = ( ρk θk+1 0