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
LSMR in one slide solve Ax = b min � Ax − b � 2 David Fong, Michael Saunders LSMR Algorithm 2/38
LSMR in one slide � � � � �� solve Ax = b A b � � min x − � � min � Ax − b � 2 λI 0 � � 2 David Fong, Michael Saunders LSMR Algorithm 2/38
LSMR in one slide � � � � �� solve Ax = b A b � � min x − � � min � Ax − b � 2 λI 0 � � 2 LSQR ≡ CG on the normal equation LSMR ≡ MINRES on the normal equation David Fong, Michael Saunders LSMR Algorithm 2/38
LSMR in one slide � � � � �� solve Ax = b A b � � min x − � � min � Ax − b � 2 λI 0 � � 2 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
LSQR Iterative algorithm for � A � � � b �� � � min x − � � λI 0 � � 2 David Fong, Michael Saunders LSMR Algorithm 3/38
LSQR Iterative algorithm for � A � � � b �� � � min x − � � λI 0 � � 2 Properties • A is rectangular ( m × n ) and often sparse • A can be an operator ( A T A + λ 2 I ) x = A T b • CG on the normal equation • Av , A T u plus O ( m + n ) operations per iteration David Fong, Michael Saunders LSMR Algorithm 3/38
Monotone convergence of residual Measure of Convergence • r k = b − Ax k r � , � A T r k � → 0 • � r k � → � ˆ David Fong, Michael Saunders LSMR Algorithm 4/38
Monotone convergence of residual Measure of Convergence • r k = b − Ax k r � , � A T r k � → 0 • � r k � → � ˆ log � A T r k � LSQR LSQR � r k � Name:lp fit1p, Dim:1677x627, nnz:9868, id=80 Name:lp fit1p, Dim:1677x627, nnz:9868, id=80 2 1.735 1.73 1 1.725 0 1.72 −1 1.715 log||A T r|| ||r|| 1.71 −2 1.705 −3 1.7 −4 1.695 1.69 −5 0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 iteration count iteration count David Fong, Michael Saunders LSMR Algorithm 4/38
Monotone convergence of residual Measure of Convergence — LSQR • r k = b − Ax k — LSMR r � , � A T r k � → 0 • � r k � → � ˆ log � A T r k � � r k � Name:lp fit1p, Dim:1677x627, nnz:9868, id=80 Name:lp fit1p, Dim:1677x627, nnz:9868, id=80 2 1.735 lsqr lsqr lsmr lsmr 1.73 1 1.725 0 1.72 −1 1.715 log||A T r|| ||r|| 1.71 −2 1.705 −3 1.7 −4 1.695 1.69 −5 0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 iteration count iteration count David Fong, Michael Saunders LSMR Algorithm 4/38
LSMR Algorithm David Fong, Michael Saunders LSMR Algorithm 5/38
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
Golub-Kahan bidiagonalization Given A ( m × n ) and b ( m × 1 ) Direct bidiagonalization U T � � b A V = B Iterative bidiagonalization 1 β 1 u 1 = b , α 1 v 1 = A T u 1 2 for k = 1 , 2 , . . . , set β k +1 u k +1 = Av k − α k u k α k +1 v k +1 = A T u k +1 − β k +1 v k David Fong, Michael Saunders LSMR Algorithm 6/38
Golub-Kahan bidiagonalization (2) The process can be summarized by b = V k ( β 1 e 1 ) AV k = U k +1 B k � I k � A T U k = V k B T k 0 where α 1 β 2 α 2 B k = ... ... β k α k β k +1 David Fong, Michael Saunders LSMR Algorithm 7/38
Golub-Kahan bidiagonalization (3) V k spans the Krylov subspace: span { v 1 , . . . , v k } = span { A T b, ( A T A ) A T b, . . . , ( A T A ) k − 1 A T b } David Fong, Michael Saunders LSMR Algorithm 8/38
Golub-Kahan bidiagonalization (3) Define x k = V k y k Subproblem to solve min y k � r k � = min y k � β 1 e 1 − B k y k � (LSQR) � � � � B T k B k y k � A T r k � = min � � ¯ min β 1 e 1 − y k (LSMR) � � ¯ β k +1 e T � � y k k � � where r k = b − Ax k , ¯ β k = α k β k David Fong, Michael Saunders LSMR Algorithm 8/38
Least squares subproblem ‚ ! ‚ B T k B k ‚ ‚ ¯ min β 1 e 1 − y k ‚ ‚ ¯ β k +1 e T ‚ ‚ y k ‚ k ‚ David Fong, Michael Saunders LSMR Algorithm 9/38
Least squares subproblem ‚ ! ‚ B T k B k ‚ ‚ ¯ min β 1 e 1 − y k ‚ ‚ ¯ β k +1 e T ‚ ‚ y k ‚ k ‚ ‚ ! ‚ R T k R k „ R k « ‚ ‚ ¯ R T k q k = ¯ = min β 1 e 1 − y k Q k +1 B k = , β k +1 e k ‚ ‚ 0 q T ‚ k R k ‚ y k ‚ ‚ David Fong, Michael Saunders LSMR Algorithm 9/38
Least squares subproblem ‚ ! ‚ B T k B k ‚ ‚ ¯ min β 1 e 1 − y k ‚ ‚ ¯ β k +1 e T ‚ ‚ y k ‚ k ‚ ‚ ! ‚ R T k R k „ R k « ‚ ‚ ¯ R T k q k = ¯ = min β 1 e 1 − y k Q k +1 B k = , β k +1 e k ‚ ‚ 0 q T ‚ k R k ‚ y k ‚ ‚ ‚ ! ‚ R T ‚ ‚ k ¯ q k = (¯ = min β 1 e 1 − t k t k = R k y k , β k +1 / ( R k ) k,k ) e k = ϕ k e k ‚ ‚ ϕ k e T ‚ ‚ t k ‚ k ‚ David Fong, Michael Saunders LSMR Algorithm 9/38
Least squares subproblem ‚ ! ‚ B T k B k ‚ ‚ ¯ min β 1 e 1 − y k ‚ ‚ ¯ β k +1 e T ‚ ‚ y k ‚ k ‚ ‚ ! ‚ R T k R k „ R k « ‚ ‚ ¯ R T k q k = ¯ = min β 1 e 1 − y k Q k +1 B k = , β k +1 e k ‚ ‚ 0 q T ‚ k R k ‚ y k ‚ ‚ ‚ ! ‚ R T ‚ ‚ k ¯ q k = (¯ = min β 1 e 1 − t k t k = R k y k , β k +1 / ( R k ) k,k ) e k = ϕ k e k ‚ ‚ ϕ k e T ‚ ‚ t k ‚ k ‚ „ z k „ ¯ ¯ ¯ ! ! R T R k z k ‚ ‚ β 1 e 1 « R k « k ¯ ‚ ‚ = min t k Q k +1 = ¯ ‚ − ‚ ζ k +1 0 ˜ ϕ k e T 0 0 ζ k +1 t k ‚ ‚ k David Fong, Michael Saunders LSMR Algorithm 9/38
Least squares subproblem ‚ ! ‚ B T k B k ‚ ‚ ¯ min β 1 e 1 − y k ‚ ‚ ¯ β k +1 e T ‚ ‚ y k ‚ k ‚ ‚ ! ‚ R T k R k „ R k « ‚ ‚ ¯ R T k q k = ¯ = min β 1 e 1 − y k Q k +1 B k = , β k +1 e k ‚ ‚ 0 q T ‚ k R k ‚ y k ‚ ‚ ‚ ! ‚ R T ‚ ‚ k ¯ q k = (¯ = min β 1 e 1 − t k t k = R k y k , β k +1 / ( R k ) k,k ) e k = ϕ k e k ‚ ‚ ϕ k e T ‚ ‚ t k ‚ k ‚ „ z k „ ¯ ¯ ¯ ! ! R T R k z k ‚ ‚ β 1 e 1 « R k « k ¯ ‚ ‚ = min t k Q k +1 = ¯ ‚ − ‚ ζ k +1 0 ˜ ϕ k e T 0 0 ζ k +1 t k ‚ ‚ k Things to note z k = ¯ x k = V k y k , t k = R k y k , R k t k , two cheap QRs David Fong, Michael Saunders LSMR Algorithm 9/38
Least squares subproblem (2) z k = ¯ Remember x k = V k y k , t k = R k y k , R k t k R k and ¯ R k both upper-bidiagonal David Fong, Michael Saunders LSMR Algorithm 10/38
Least squares subproblem (2) z k = ¯ Remember x k = V k y k , t k = R k y k , R k t k R k and ¯ R k both upper-bidiagonal Key steps to compute x k x k = V k y k David Fong, Michael Saunders LSMR Algorithm 10/38
Least squares subproblem (2) z k = ¯ Remember x k = V k y k , t k = R k y k , R k t k R k and ¯ R k both upper-bidiagonal Key steps to compute x k x k = V k y k R T k W T k = V T = W k t k k David Fong, Michael Saunders LSMR Algorithm 10/38
Least squares subproblem (2) z k = ¯ Remember x k = V k y k , t k = R k y k , R k t k R k and ¯ R k both upper-bidiagonal Key steps to compute x k x k = V k y k R T k W T k = V T = W k t k k = ¯ R T ¯ k ¯ W T k = W T W k z k k David Fong, Michael Saunders LSMR Algorithm 10/38
Least squares subproblem (2) z k = ¯ Remember x k = V k y k , t k = R k y k , R k t k R k and ¯ R k both upper-bidiagonal Key steps to compute x k x k = V k y k R T k W T k = V T = W k t k k = ¯ R T ¯ k ¯ W T k = W T W k z k k = x k − 1 + ζ k ¯ w k � T � where z k = ζ 1 ζ 2 · · · ζ k David Fong, Michael Saunders LSMR Algorithm 10/38
Flow chart of LSMR A, b Golub-Kahan bidiagonalization QR decompositions Update x Estimate large backward error small Solution x David Fong, Michael Saunders LSMR Algorithm 11/38
Flow chart of LSMR A, b Computational cost Av, A T v, 3 m + 3 n Golub-Kahan bidiagonalization O (1) QR decompositions Update x 3 n Estimate large backward error O (1) small Solution x David Fong, Michael Saunders LSMR Algorithm 11/38
Computational and storage requirement Storage Work m n m n MINRES on A T Ax = A T b Av 1 x, v 1 , v 2 , w 1 , w 2 8 LSQR Av, u x, v, w 3 5 x, v, h, ¯ Av, u h LSMR 3 6 where h k , ¯ h k are scalar multiples of w k , ¯ w k David Fong, Michael Saunders LSMR Algorithm 12/38
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
Recommend
More recommend
