LSMR: An iterative algorithm for sparse least-squares problems - - PowerPoint PPT Presentation

▶
lsmr an iterative algorithm for sparse least squares
SMART_READER_LITE
LIVE PREVIEW

LSMR: An iterative algorithm for sparse least-squares problems - - PowerPoint PPT Presentation

Dec 25, 2023 367 likes •1.09k views

LSMR: An iterative algorithm for sparse least-squares problems David Fong Michael Saunders Institute for Computational and Mathematical Engineering (iCME) Stanford University 2nd IMA Conference on Numerical Linear Algebra and Optimisation

slide-1
SLIDE 1

LSMR: An iterative algorithm for sparse least-squares problems

David Fong Michael Saunders

Institute for Computational and Mathematical Engineering (iCME) Stanford University

2nd IMA Conference on Numerical Linear Algebra and Optimisation University of Birmingham, UK September 13–15, 2010

David Fong, Michael Saunders LSMR Algorithm 1/38

slide-2
SLIDE 2

LSMR in one slide

solve Ax = b min Ax − b2

David Fong, Michael Saunders LSMR Algorithm 2/38

slide-3
SLIDE 3

LSMR in one slide

solve Ax = b min

  • A

λI

  • x −
  • b
  • 2

min Ax − b2

David Fong, Michael Saunders LSMR Algorithm 2/38

slide-4
SLIDE 4

LSMR in one slide

solve Ax = b min

  • A

λI

  • x −
  • b
  • 2

min Ax − b2

LSQR ≡ CG on the normal equation LSMR ≡ MINRES on the normal equation

David Fong, Michael Saunders LSMR Algorithm 2/38

slide-5
SLIDE 5

LSMR in one slide

solve Ax = b min

  • A

λI

  • x −
  • b
  • 2

min Ax − b2

LSQR ≡ CG on the normal equation LSMR ≡ MINRES on the normal equation

  • Almost same complexity as LSQR
  • Better convergence properties for

inexact solves

David Fong, Michael Saunders LSMR Algorithm 2/38

slide-6
SLIDE 6

LSQR

Iterative algorithm for

min

  • A

λI

  • x −

b

  • 2

David Fong, Michael Saunders LSMR Algorithm 3/38

slide-7
SLIDE 7

LSQR

Iterative algorithm for

min

  • A

λI

  • x −

b

  • 2

Properties

  • A is rectangular (m × n) and often sparse
  • A can be an operator
  • CG on the normal equation

(ATA + λ2I)x = ATb

  • Av, ATu plus O(m + n) operations per iteration

David Fong, Michael Saunders LSMR Algorithm 3/38

slide-8
SLIDE 8

Monotone convergence of residual

Measure of Convergence

  • rk = b − Axk
  • rk → ˆ

r, ATrk → 0

David Fong, Michael Saunders LSMR Algorithm 4/38

slide-9
SLIDE 9

Monotone convergence of residual

Measure of Convergence

  • rk = b − Axk
  • rk → ˆ

r, ATrk → 0

LSQR rk

50 100 150 200 250 300 350 400 450 500 1.69 1.695 1.7 1.705 1.71 1.715 1.72 1.725 1.73 1.735

Name:lp fit1p, Dim:1677x627, nnz:9868, id=80

iteration count ||r||

LSQR log ATrk

50 100 150 200 250 300 350 400 450 500 −5 −4 −3 −2 −1 1 2

Name:lp fit1p, Dim:1677x627, nnz:9868, id=80

iteration count log||ATr||

David Fong, Michael Saunders LSMR Algorithm 4/38

slide-10
SLIDE 10

Monotone convergence of residual

Measure of Convergence

  • rk = b − Axk
  • rk → ˆ

r, ATrk → 0 — LSQR — LSMR

rk

50 100 150 200 250 300 350 400 450 500 1.69 1.695 1.7 1.705 1.71 1.715 1.72 1.725 1.73 1.735 iteration count ||r||

Name:lp fit1p, Dim:1677x627, nnz:9868, id=80

lsqr lsmr

log ATrk

50 100 150 200 250 300 350 400 450 500 −5 −4 −3 −2 −1 1 2 iteration count log||ATr||

Name:lp fit1p, Dim:1677x627, nnz:9868, id=80

lsqr lsmr

David Fong, Michael Saunders LSMR Algorithm 4/38

slide-11
SLIDE 11

LSMR Algorithm

David Fong, Michael Saunders LSMR Algorithm 5/38

slide-12
SLIDE 12

Golub-Kahan bidiagonalization

Given A (m × n) and b (m × 1)

Direct bidiagonalization

U T b A

  • V = B

David Fong, Michael Saunders LSMR Algorithm 6/38

slide-13
SLIDE 13

Golub-Kahan bidiagonalization

Given A (m × n) and b (m × 1)

Direct bidiagonalization

U T b A

  • V = B

Iterative bidiagonalization

1 β1u1 = b, α1v1 = ATu1 2 for k = 1, 2, . . . , set

βk+1uk+1 = Avk − αkuk αk+1vk+1 = ATuk+1 − βk+1vk

David Fong, Michael Saunders LSMR Algorithm 6/38

slide-14
SLIDE 14

Golub-Kahan bidiagonalization (2)

The process can be summarized by b = Vk(β1e1) AVk = Uk+1Bk ATUk = VkBT

k

Ik

  • where

Bk =        α1 β2 α2 ... ... βk αk βk+1       

David Fong, Michael Saunders LSMR Algorithm 7/38

slide-15
SLIDE 15

Golub-Kahan bidiagonalization (3)

Vk spans the Krylov subspace: span{v1, . . . , vk} = span{ATb, (ATA)ATb, . . . , (ATA)k−1ATb}

David Fong, Michael Saunders LSMR Algorithm 8/38

slide-16
SLIDE 16

Golub-Kahan bidiagonalization (3)

Define xk = Vkyk

Subproblem to solve

min

yk rk = min yk β1e1 − Bkyk

(LSQR) min

yk AT rk = min yk

  • ¯

β1e1 −

  • BT

k Bk

¯ βk+1eT

k

  • yk
  • (LSMR)

where rk = b − Axk, ¯ βk = αkβk

David Fong, Michael Saunders LSMR Algorithm 8/38

slide-17
SLIDE 17

Least squares subproblem

min

yk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − BT

k Bk

¯ βk+1eT

k

! yk ‚ ‚ ‚ ‚ ‚

David Fong, Michael Saunders LSMR Algorithm 9/38

slide-18
SLIDE 18

Least squares subproblem

min

yk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − BT

k Bk

¯ βk+1eT

k

! yk ‚ ‚ ‚ ‚ ‚ = min

yk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − RT

k Rk

qT

k Rk

! yk ‚ ‚ ‚ ‚ ‚ Qk+1Bk = „Rk « , RT

k qk = ¯

βk+1ek

David Fong, Michael Saunders LSMR Algorithm 9/38

slide-19
SLIDE 19

Least squares subproblem

min

yk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − BT

k Bk

¯ βk+1eT

k

! yk ‚ ‚ ‚ ‚ ‚ = min

yk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − RT

k Rk

qT

k Rk

! yk ‚ ‚ ‚ ‚ ‚ Qk+1Bk = „Rk « , RT

k qk = ¯

βk+1ek = min

tk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − RT

k

ϕkeT

k

! tk ‚ ‚ ‚ ‚ ‚ tk = Rkyk, qk = (¯ βk+1/(Rk)k,k)ek = ϕkek

David Fong, Michael Saunders LSMR Algorithm 9/38

slide-20
SLIDE 20

Least squares subproblem

min

yk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − BT

k Bk

¯ βk+1eT

k

! yk ‚ ‚ ‚ ‚ ‚ = min

yk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − RT

k Rk

qT

k Rk

! yk ‚ ‚ ‚ ‚ ‚ Qk+1Bk = „Rk « , RT

k qk = ¯

βk+1ek = min

tk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − RT

k

ϕkeT

k

! tk ‚ ‚ ‚ ‚ ‚ tk = Rkyk, qk = (¯ βk+1/(Rk)k,k)ek = ϕkek = min

tk

‚ ‚ ‚ ‚ „ zk ¯ ζk+1 « − „ ¯ Rk « tk ‚ ‚ ‚ ‚ ¯ Qk+1 RT

k

¯ β1e1 ϕkeT

k

! = ¯ Rk zk ˜ ζk+1 !

David Fong, Michael Saunders LSMR Algorithm 9/38

slide-21
SLIDE 21

Least squares subproblem

min

yk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − BT

k Bk

¯ βk+1eT

k

! yk ‚ ‚ ‚ ‚ ‚ = min

yk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − RT

k Rk

qT

k Rk

! yk ‚ ‚ ‚ ‚ ‚ Qk+1Bk = „Rk « , RT

k qk = ¯

βk+1ek = min

tk

‚ ‚ ‚ ‚ ‚ ¯ β1e1 − RT

k

ϕkeT

k

! tk ‚ ‚ ‚ ‚ ‚ tk = Rkyk, qk = (¯ βk+1/(Rk)k,k)ek = ϕkek = min

tk

‚ ‚ ‚ ‚ „ zk ¯ ζk+1 « − „ ¯ Rk « tk ‚ ‚ ‚ ‚ ¯ Qk+1 RT

k

¯ β1e1 ϕkeT

k

! = ¯ Rk zk ˜ ζk+1 !

Things to note

xk = Vkyk, tk = Rkyk, zk = ¯ Rktk, two cheap QRs

David Fong, Michael Saunders LSMR Algorithm 9/38

slide-22
SLIDE 22

Least squares subproblem (2)

Remember xk = Vkyk, tk = Rkyk, zk = ¯ Rktk Rk and ¯ Rk both upper-bidiagonal

David Fong, Michael Saunders LSMR Algorithm 10/38

slide-23
SLIDE 23

Least squares subproblem (2)

Remember xk = Vkyk, tk = Rkyk, zk = ¯ Rktk Rk and ¯ Rk both upper-bidiagonal

Key steps to compute xk

xk = Vkyk

David Fong, Michael Saunders LSMR Algorithm 10/38

slide-24
SLIDE 24

Least squares subproblem (2)

Remember xk = Vkyk, tk = Rkyk, zk = ¯ Rktk Rk and ¯ Rk both upper-bidiagonal

Key steps to compute xk

xk = Vkyk = Wktk RT

k W T k = V T k

David Fong, Michael Saunders LSMR Algorithm 10/38

slide-25
SLIDE 25

Least squares subproblem (2)

Remember xk = Vkyk, tk = Rkyk, zk = ¯ Rktk Rk and ¯ Rk both upper-bidiagonal

Key steps to compute xk

xk = Vkyk = Wktk RT

k W T k = V T k

= ¯ Wkzk ¯ RT

k ¯

W T

k = W T k

David Fong, Michael Saunders LSMR Algorithm 10/38

slide-26
SLIDE 26

Least squares subproblem (2)

Remember xk = Vkyk, tk = Rkyk, zk = ¯ Rktk Rk and ¯ Rk both upper-bidiagonal

Key steps to compute xk

xk = Vkyk = Wktk RT

k W T k = V T k

= ¯ Wkzk ¯ RT

k ¯

W T

k = W T k

= xk−1 + ζk ¯ wk where zk =

  • ζ1

ζ2 · · · ζk T

David Fong, Michael Saunders LSMR Algorithm 10/38

slide-27
SLIDE 27

Flow chart of LSMR

A, b

Golub-Kahan bidiagonalization

QR decompositions Estimate backward error large

Solution x

small Update x

David Fong, Michael Saunders LSMR Algorithm 11/38

slide-28
SLIDE 28

Flow chart of LSMR

Computational cost

Av, ATv, 3m + 3n O(1) 3n O(1)

A, b

Golub-Kahan bidiagonalization

QR decompositions Estimate backward error large

Solution x

small Update x

David Fong, Michael Saunders LSMR Algorithm 11/38

slide-29
SLIDE 29

Computational and storage requirement

Storage Work m n m n MINRES on ATAx = ATb Av1 x, v1, v2, w1, w2 8 LSQR Av, u x, v, w 3 5 LSMR Av, u x, v, h, ¯ h 3 6 where hk, ¯ hk are scalar multiples of wk, ¯ wk

David Fong, Michael Saunders LSMR Algorithm 12/38

slide-30
SLIDE 30

Numerical experiments

Test Data

  • University of Florida Sparse Matrix Collection
  • LPnetlib: Linear Programming Problems
  • A = (Problem.A)’

b = Problem.c (127 problems)

David Fong, Michael Saunders LSMR Algorithm 13/38

slide-31
SLIDE 31

Numerical experiments

Test Data

  • University of Florida Sparse Matrix Collection
  • LPnetlib: Linear Programming Problems
  • A = (Problem.A)’

b = Problem.c (127 problems)

Solve min Ax − b2

with LSQR and LSMR

  • Examples of rk
  • Backward error tests: nnz(A) ≤ 63220
  • Reorthogonalization: nnz(A) ≤ 15977

David Fong, Michael Saunders LSMR Algorithm 13/38

slide-32
SLIDE 32

rk for LSQR and LSMR – typical

50 100 150 200 250 300 600 700 800 900 1000 1100 1200 1300 1400 iteration count ||r||

Name:lp greenbeb, Dim:5598x2392, nnz:31070, id=100

lsqr minres lsmr

David Fong, Michael Saunders LSMR Algorithm 14/38

slide-33
SLIDE 33

rk for LSQR and LSMR – rare

10 20 30 40 50 60 70 80 90 1.975 1.98 1.985 1.99 1.995 2 iteration count ||r||

Name:lp woodw, Dim:8418x1098, nnz:37487, id=104

lsqr minres lsmr

David Fong, Michael Saunders LSMR Algorithm 15/38

slide-34
SLIDE 34

Backward error – estimates

(A + Ei)T (A + Ei)x = (A + Ei)T b

David Fong, Michael Saunders LSMR Algorithm 16/38

slide-35
SLIDE 35

Backward error – estimates

(A + Ei)T (A + Ei)x = (A + Ei)T b Two estimates given by Stewart (1975 and 1977) E1 = exT

x2

E1 = e

x

e = ˆ r − r E2 = −rrTA

r2

E2 = ATr

r

where ˆ r is the residual for the exact solution

David Fong, Michael Saunders LSMR Algorithm 16/38

slide-36
SLIDE 36

Backward error – estimates

(A + Ei)T (A + Ei)x = (A + Ei)T b Two estimates given by Stewart (1975 and 1977) E1 = exT

x2

E1 = e

x

e = ˆ r − r E2 = −rrTA

r2

E2 = ATr

r

where ˆ r is the residual for the exact solution

Note

E2 is computable

David Fong, Michael Saunders LSMR Algorithm 16/38

slide-37
SLIDE 37

log10 E2 for LSQR and LSMR – typical

200 400 600 800 1000 1200 1400 1600 1800 −7 −6 −5 −4 −3 −2 −1 iteration count log(E2)

Name:lp pilot ja, Dim:2267x940, nnz:14977, id=88

E2 LSQR E2 LSMR

David Fong, Michael Saunders LSMR Algorithm 17/38

slide-38
SLIDE 38

log10 E2 for LSQR and LSMR – rare

20 40 60 80 100 120 −8 −7 −6 −5 −4 −3 −2 −1 1 iteration count log(E2)

Name:lp sc205, Dim:317x205, nnz:665, id=17

E2 LSQR E2 LSMR

David Fong, Michael Saunders LSMR Algorithm 18/38

slide-39
SLIDE 39

Backward error - optimal

µ(x) ≡ min

E E

st (A + E)T (A + E)x = (A + E)T b Exact µ(x)

(Waldén, Karlson, & Sun 1995, Higham 2002)

C ≡

  • A

r x

  • I − rrT

r2

  • µ(x) = σmin(C)

David Fong, Michael Saunders LSMR Algorithm 19/38

slide-40
SLIDE 40

Backward error - optimal

µ(x) ≡ min

E E

st (A + E)T (A + E)x = (A + E)T b Cheaper estimate ˜ µ(x)

(Grcar, Saunders, & Su 2007)

K =

  • A

r xI

  • v =

r

  • miny Ky − v

˜ µ(x) = Ky

x

David Fong, Michael Saunders LSMR Algorithm 19/38

slide-41
SLIDE 41

Backward error - optimal

µ(x) ≡ min

E E

st (A + E)T (A + E)x = (A + E)T b Cheaper estimate ˜ µ(x)

(Grcar, Saunders, & Su 2007)

K =

  • A

r xI

  • v =

r

  • miny Ky − v

˜ µ(x) = Ky

x

r = b - A*x; p = colamd(A); eta = norm(r)/norm(x); K = [A(:,p); eta*speye(n)]; v = [r; zeros(n,1)]; [c,R] = qr(K,v,0); mutilde = norm(c)/norm(x);

David Fong, Michael Saunders LSMR Algorithm 19/38

slide-42
SLIDE 42

Backward errors for LSQR – typical

200 400 600 800 1000 1200 1400 1600 1800 −8 −7 −6 −5 −4 −3 −2 −1 1 iteration count log(Backward Error)

Name:lp cre a, Dim:7248x3516, nnz:18168, id=93

E1 LSQR E2 LSQR Optimal LSQR

David Fong, Michael Saunders LSMR Algorithm 20/38

slide-43
SLIDE 43

Backward errors for LSQR – rare

100 200 300 400 500 600 −7 −6 −5 −4 −3 −2 −1 1 iteration count log(Backward Error)

Name:lp pilot, Dim:4860x1441, nnz:44375, id=107

E1 LSQR E2 LSQR Optimal LSQR

David Fong, Michael Saunders LSMR Algorithm 21/38

slide-44
SLIDE 44

Backward errors for LSMR – typical

200 400 600 800 1000 1200 1400 1600 1800 −8 −7 −6 −5 −4 −3 −2 −1 1 iteration count log(Backward Error)

Name:lp cre a, Dim:7248x3516, nnz:18168, id=93

E1 LSMR E2 LSMR Optimal LSMR

David Fong, Michael Saunders LSMR Algorithm 22/38

slide-45
SLIDE 45

Backward errors for LSMR – rare

100 200 300 400 500 600 −7 −6 −5 −4 −3 −2 −1 1 iteration count log(Backward Error)

Name:lp pilot, Dim:4860x1441, nnz:44375, id=107

E1 LSMR E2 LSMR Optimal LSMR

David Fong, Michael Saunders LSMR Algorithm 23/38

slide-46
SLIDE 46

For LSMR

E2 ≈ optimal BE almost always

Typical: E2 ≈ ˜

µ(x)

50 100 150 200 250 −8 −7 −6 −5 −4 −3 −2 −1 1 iteration count log(Backward Error)

Name:lp ken 11, Dim:21349x14694, nnz:49058, id=108

E1 LSMR E2 LSMR Optimal LSMR

Rare: E1 ≈ ˜

µ(x)

10 20 30 40 50 60 70 80 90 −9 −8 −7 −6 −5 −4 −3 −2 −1 1 iteration count log(Backward Error)

Name:lp ship12l, Dim:5533x1151, nnz:16276, id=91

E1 LSMR E2 LSMR Optimal LSMR

David Fong, Michael Saunders LSMR Algorithm 24/38

slide-47
SLIDE 47

For LSMR, optimal BE ˜ µ(x) seems to be monotonic For LSQR, usually not Typical for LSQR and LSMR

1000 2000 3000 4000 5000 6000 7000 −7 −6 −5 −4 −3 −2 −1 iteration count log(||ATr||/||r||)

Name:lp maros, Dim:1966x846, nnz:10137, id=81

Optimal LSQR Optimal LSMR

Rare LSQR, typical LSMR

200 400 600 800 1000 1200 1400 1600 −8 −7 −6 −5 −4 −3 −2 −1 iteration count log(||ATr||/||r||)

Name:lp cre c, Dim:6411x3068, nnz:15977, id=90

Optimal LSQR Optimal LSMR

David Fong, Michael Saunders LSMR Algorithm 25/38

slide-48
SLIDE 48

Optimal backward errors ˜ µ(xLSMR) ≤ ˜ µ(xLSQR) almost always Typical

100 200 300 400 500 600 −7 −6 −5 −4 −3 −2 −1 iteration count log(||ATr||/||r||)

Name:lp pilot, Dim:4860x1441, nnz:44375, id=107

Optimal LSQR Optimal LSMR

Rare

20 40 60 80 100 120 140 −9 −8 −7 −6 −5 −4 −3 −2 −1 iteration count log(||ATr||/||r||)

Name:lp standgub, Dim:1383x361, nnz:3338, id=50

Optimal LSQR Optimal LSMR

David Fong, Michael Saunders LSMR Algorithm 26/38

slide-49
SLIDE 49

Errors

  • xLSQR − x∗ is monotonic
  • xLSMR − x∗ seems to be monotonic
  • xLSQR − x∗ ≤ xLSMR − x∗

David Fong, Michael Saunders LSMR Algorithm 27/38

slide-50
SLIDE 50

Errors

  • xLSQR − x∗ is monotonic
  • xLSMR − x∗ seems to be monotonic
  • xLSQR − x∗ ≤ xLSMR − x∗

LSQR error ≤ LSMR error

10 20 30 40 50 60 70 80 90 −3 −2 −1 1 2 3 4 5 6 iteration count log||xk − x*|

Name:lp ship12l, Dim:5533x1151, nnz:16276, id=91

lsqr lsmr

Both give min-length x

10 20 30 40 50 60 70 −2 −1 1 2 3 4 5 6 7 iteration count log||xk − x*|

Name:lp pds 02, Dim:7716x2953, nnz:16571, id=92

lsqr lsmr

David Fong, Michael Saunders LSMR Algorithm 27/38

slide-51
SLIDE 51

Space-time trade-offs

LSMR is well-suited for limited memory computations. What if we have

  • more memory
  • Av expensive

Can we speed things up?

David Fong, Michael Saunders LSMR Algorithm 28/38

slide-52
SLIDE 52

Space-time trade-offs

LSMR is well-suited for limited memory computations. What if we have

  • more memory
  • Av expensive

Can we speed things up? Some ideas:

  • Reorthogonalization
  • Restarting
  • Local reorthogonalization

David Fong, Michael Saunders LSMR Algorithm 28/38

slide-53
SLIDE 53

Reorthogonalization

Golub-Kahan process

Infinite precision Finite precision Uk, Vk orthonormal Lose orthogonality At most min(m, n) iterations Could take 10n or more

David Fong, Michael Saunders LSMR Algorithm 29/38

slide-54
SLIDE 54

Reorthogonalization

Golub-Kahan process

Infinite precision Finite precision Uk, Vk orthonormal Lose orthogonality At most min(m, n) iterations Could take 10n or more Apply modified Gram-Schmidt to uk+1 and/or vk+1: u ← u − (uT

j u)uj

j = k, k−1, k−2, . . . (similarly for v)

David Fong, Michael Saunders LSMR Algorithm 29/38

slide-55
SLIDE 55

Effects of reorthogonalization on various problems

10 20 30 40 50 60 70 80 90 −8 −7 −6 −5 −4 −3 −2 −1 1 iteration count log(E2)

Name:lp ship12l, Dim:5533x1151, nnz:16276, id=91

NoOrtho OrthoU OrthoV OrthoUV 100 200 300 400 500 600 700 800 900 −7 −6 −5 −4 −3 −2 −1 1 iteration count log(E2)

Name:lp scfxm2, Dim:1200x660, nnz:5469, id=62

NoOrtho OrthoU OrthoV OrthoUV 20 40 60 80 100 120 140 160 180 −6 −5 −4 −3 −2 −1 1 iteration count log(E2)

Name:lp agg2, Dim:758x516, nnz:4740, id=58

NoOrtho OrthoU OrthoV OrthoUV 2000 4000 6000 8000 10000 12000 −7 −6 −5 −4 −3 −2 −1 iteration count log(E2)

Name:lpi gran, Dim:2525x2658, nnz:20111, id=94

NoOrtho OrthoU OrthoV OrthoUV

David Fong, Michael Saunders LSMR Algorithm 30/38

slide-56
SLIDE 56

Orthogonality of Uk

100 200 300 400 500 600 5 6

Original

100 200 300 400 500 600 5 6 7 8 9 10 11 1213 14 15 16

MGS on U

100 200 300 400 500 600 5 6 7 8 9 10 11 12 13 14 15 16

MGS on V

100 200 300 400 500 600 5 6 7 8 9 10 11 12 13 14 15 16

MGS on U,V

David Fong, Michael Saunders LSMR Algorithm 31/38

slide-57
SLIDE 57

Orthogonality of Vk

100 200 300 400 500 600 5 6

Original

100 200 300 400 500 600 5 6 7 8 9 10 11 1213 14 15 16

MGS on U

100 200 300 400 500 600 5 6 7 8 9 10 11 12 13 14 15 16

MGS on V

100 200 300 400 500 600 5 6 7 8 9 10 11 12 13 14 15 16

MGS on U,V

David Fong, Michael Saunders LSMR Algorithm 32/38

slide-58
SLIDE 58

What we learnt so far

  • Reorthogonalizing Vk (only) is sufficient
  • Reorthogonalizing Uk (only) is nearly as good
  • xk converges the same for all options

What can be improved

  • May still use too much memory
  • Need more flexibility for space-time trade-off

David Fong, Michael Saunders LSMR Algorithm 33/38

slide-59
SLIDE 59

Reorthogonalization with Restarting

Restarting LSMR

rk = b − Axk min A∆x − rk

David Fong, Michael Saunders LSMR Algorithm 34/38

slide-60
SLIDE 60

Reorthogonalization with Restarting

Restarting LSMR

rk = b − Axk min A∆x − rk

Restarting leads to stagnation

500 1000 1500 2000 2500 3000 3500 4000 4500 −7 −6 −5 −4 −3 −2 −1 iteration count Backward Error

Name:lp maros, Dim:1966x846, nnz:10137, id=81

NoOrtho Restart5 Restart10 Restart50 NoRestart 2000 4000 6000 8000 10000 12000 14000 16000 −7 −6 −5 −4 −3 −2 −1 1 iteration count Backward Error

Name:lp cre c, Dim:6411x3068, nnz:15977, id=90

NoOrtho Restart5 Restart10 Restart50 NoRestart

David Fong, Michael Saunders LSMR Algorithm 34/38

slide-61
SLIDE 61

Local reorthogonalization

  • Reorthogonalize wrto only the last l vectors
  • Partial speed-up
  • Less memory
  • Depends on efficiency of Av and ATu

David Fong, Michael Saunders LSMR Algorithm 35/38

slide-62
SLIDE 62

Local reorthogonalization

  • Reorthogonalize wrto only the last l vectors
  • Partial speed-up
  • Less memory
  • Depends on efficiency of Av and ATu

Speed up with local reorthogonalization

50 100 150 200 250 300 350 400 450 −7 −6.5 −6 −5.5 −5 −4.5 −4 −3.5 −3 −2.5 −2 iteration count Backward Error

Name:lp fit1p, Dim:1677x627, nnz:9868, id=80

NoOrtho Local5 Local10 Local50 NoLocal 200 400 600 800 1000 1200 1400 −6 −5 −4 −3 −2 −1 1 2 iteration count Backward Error

Name:lp bnl2, Dim:4486x2324, nnz:14996, id=89

NoOrtho Local5 Local10 Local50 NoLocal

David Fong, Michael Saunders LSMR Algorithm 35/38

slide-63
SLIDE 63

Conclusions

LSMR has the good properties of LSQR and more

  • rk seems monotonic (nearly as small as for LSQR)
  • xk − x∗ also
  • ATrk is monotonic

David Fong, Michael Saunders LSMR Algorithm 36/38

slide-64
SLIDE 64

Conclusions

LSMR has the good properties of LSQR and more

  • rk seems monotonic (nearly as small as for LSQR)
  • xk − x∗ also
  • ATrk is monotonic

Stewart backward errors E1 = ˆ

r−rk xk

E2 = ATrk

rk

  • E1 monotonic if rk monotonic (theorem)

David Fong, Michael Saunders LSMR Algorithm 36/38

slide-65
SLIDE 65

Conclusions

LSMR has the good properties of LSQR and more

  • rk seems monotonic (nearly as small as for LSQR)
  • xk − x∗ also
  • ATrk is monotonic

Stewart backward errors E1 = ˆ

r−rk xk

E2 = ATrk

rk

  • E1 monotonic if rk monotonic (theorem)
  • E2 usually monotonic (1 exception in 127 cases)

David Fong, Michael Saunders LSMR Algorithm 36/38

slide-66
SLIDE 66

Conclusions

LSMR has the good properties of LSQR and more

  • rk seems monotonic (nearly as small as for LSQR)
  • xk − x∗ also
  • ATrk is monotonic

Stewart backward errors E1 = ˆ

r−rk xk

E2 = ATrk

rk

  • E1 monotonic if rk monotonic (theorem)
  • E2 usually monotonic (1 exception in 127 cases)
  • E2 ≈ optimal backward error (seems monotonic)

David Fong, Michael Saunders LSMR Algorithm 36/38

slide-67
SLIDE 67

Conclusions

LSMR has the good properties of LSQR and more

  • rk seems monotonic (nearly as small as for LSQR)
  • xk − x∗ also
  • ATrk is monotonic

Stewart backward errors E1 = ˆ

r−rk xk

E2 = ATrk

rk

  • E1 monotonic if rk monotonic (theorem)
  • E2 usually monotonic (1 exception in 127 cases)
  • E2 ≈ optimal backward error (seems monotonic)
  • ⇒ reliable rule for stopping early

David Fong, Michael Saunders LSMR Algorithm 36/38

slide-68
SLIDE 68

Acknowledgements

  • David Fong
  • Chris Paige
  • Stanford Graduate Fellowship
  • ONR and AHPCRC

David Fong, Michael Saunders LSMR Algorithm 37/38

slide-69
SLIDE 69

Acknowledgements

  • David Fong
  • Chris Paige
  • Stanford Graduate Fellowship
  • ONR and AHPCRC
  • First presented at 2010 Copper Mountain Conference
  • Jon Claerbout (Stanford Geophysics)

Nov 2009: Computational success ≡ 1 − AT rk/AT b

David Fong, Michael Saunders LSMR Algorithm 37/38

slide-70
SLIDE 70

Paper and Implementations

http://www.stanford.edu/group/SOL/software.html Report SOL 2010-2 submitted to SISC Matlab and F90 code

David Fong, Michael Saunders LSMR Algorithm 38/38