LSQR: An Algorithm For Sparse Linear Equations And Sparse .

2y ago
9 Views
2 Downloads
1.43 MB
29 Pages
Last View : 2m ago
Last Download : 2m ago
Upload by : Aiyana Dorn
Transcription

LSQR: An Algorithm for Sparse LinearEquations and Sparse Least SquaresCHRISTOPHER C. PAIGEMcGill University, CanadaandMICHAEL A. SAUNDERSStanford UniversityAn iterative method is given for solving Ax ffib and minUAx - b 112,where the matrix A is large andsparse. The method is based on the bidiagonalization procedure of Golub and Kahan. It is analyticallyequivalent to the standard method of conjugate gradients, but possesses more favorable numericalproperties.Reliable stopping criteria are derived, along with estimates of standard errors for x and thecondition number of A. These are used in the FORTRAN implementation of the method, subroutineLSQR. Numerical tests are described comparing I QR with several other conjugate-gradient algorithms, indicating that I QR is the most reliable algorithm when A is ill-conditioned.Categories and Subject Descriptors: G.1.2 [Numerical Analysis]: ApprorJmation--least squaresapproximation; G.1.3 [Numerical Analysis]: Numerical Linear Algebra--linear systems (direct andtterative methods); sparse and very large systemsGeneral Terms: AlgorithmsAdditional Key Words and Phrases: analysis of varianceThe Algorithm: LSQR: Sparse Linear Equations and Least Square Problems. ACM Trans. Math.Softw. 8, 2 (June 1982).1. INTRODUCTIONA n u m e r i c a l m e t h o d is p r e s e n t e d h e r e for c o m p u t i n g a s o l u t i o n x t o e i t h e r of t h efollowing p r o b l e m s :Unsymmetric equations:solve A x ffi bLinear least squares:m i n i m i z e ][A x - b 112This work was supported in part by the Natural Sciences and Engineering Research Council ofCanada Grant A8652, the New Zealand Department of Scientific and Industrial Research, the U.S.Office of Naval Research under Contract N00014-75-C-0267, National Science Foundation GrantsMCS 76-20019 and ENG 77-06761, and the Department of Energy under Contract DE.AC0376SF00326, PA No. DE-AT03-76ER72018.Authors' addresses: C.C. Paige, School of Computer Science, McGill University, Montreal, P.Q.,Canada H3C 3G1; M.A. Saunders, Department of Operations Research, Stanford University,Stanford,CA 94305.Permission to copy without fee all or part of this material is granted provided that the copies are notmade or distributed for direct commercial advantage, the ACM copyright notice and the title of thepublication and its date appear, and notice is given that copying is by permission of the Associationfor Computing Machinery. To copy otherwise, or to republish, requires a fee and/or specificpermission. 1982 ACM 0098-3500/82/0300-0043 00.75ACMTransactionson MathematwalSoftware, Vol 8, No. 1, March 1982,Pages 43-71.

44 C . C . Paige and M. A. Saunderswhere A is a real matrix with m rows and n columns and b is a real vector. It willusually be true that m n and rank(A) n, but these conditions are notessential. The method, to be called algorithm LSQR, is similar in style to thewell-known method of conjugate gradients (CG) as applied to the least-squaresproblem [10]. The matrix A is used only to compute products of the form A v andATu for various vectors v and u. Hence A will normally be large and sparse or willbe expressible as a product of matrices that are sparse or have special structure.A typical application is to the large least-squares problems arising from analysisof variance (6.g., [12]).CG-like methods are iterative in nature. They are characterized by t h e i r needfor only a few vectors of working storage and by their theoretical convergencewithin at most n iterations (if exact arithmetic could be performed). In practicesuch methods may require far fewer or far more than n iterations to reach anacceptable approximation to x. The methods are most useful when A is wellconditioned and has many nearly equal singular values. These properties occurnaturally in many applications. In other cases it is often possible to divide thesolution procedure into a direct and an iterative part, such that the iterative parthas a better conditioned matrix for which CG-like methods will converge morequickly. Some such transformation methods are considered in [21].Algorithm LSQR is based on the bidiagonalization procedure of Golub andKahan [9]. It generates a sequence of approximations {xk } such that the residualnorm II rk [[2 decreases monotonically, where r k b - A x k . Analytically, thesequence (xh} is identical to the sequence generated by the standard CG algorithm and by several other published algorithms. However, LSQR is shown byexample to be numerically more reliable in various circumstances than the othermethods considered.The FORTRAN implementation of LSQR [22] is designed for practical application. It incorporates reliable stopping criteria and provides the user withcomputed estimates of the following quantities: x , r b - A x , A Tr, II r 112, It A IIF,standard errors for x , and the condition number of A.1.1NotationMatrices are denoted by A, B . . . . . vectors by v, w , . . . , and scalars by , f l , . . . .Two exceptions are c and s, which denote the significant components of anelementary orthogonal matrix, such that c 2 s 2 1. For a vector v, Hv [[ alwaysdenotes the Euclidean norm Hv [[2 (vTv) 1/2. For a matrix A, [[A [[ will usuallymean the Frobenius norm, HA HF ( a 2 ) 1/2, and the condition number for anunsymmetric matrix A is defined by cond(A) ]]A ]] IIA ]], where A denotes thepseudoinverse of A. The relative precision of floating-point arithmetic is e, thesmallest machine-representable number such that 1 e 1.2. MOTIVATION VIA THE LANCZOS PROCESSIn this section we review the symmetric Lanczos process [13] and its use insolving symmetric linear equations B x b. Algorithm LSQR is then derived byapplying the Lanczos process to a particular symmetric system. Although a moredirect development is given in Section 4, the present derivation may remainuseful for a future error analysis of LSQR, since many of the rounding errorproperties of the Lanczos process are already known [18].ACM Transactions on Mathematmal Software, Vol 8, No 1, March 1982.

LSQR: An Algorithm for Sparse Linear Equattons and Sparse Least Squares 45Given a symmetric matrix B and a starting vector b, the Lanczos process is amethod for generating a sequence of vectors {v,) and scalars {a, ), (fli} such thatB is reduced to tridiagonal form. A reliable computational form of the method isthe following.T h e L a n c z o s p r o c e s s (reduction to tridiagonal form):fll vl -- b,w, Bvt - fl, v,- ]a, t l Vt lI'vTw'Wt i 1, 2 , . . . ,(2.1)OltVtwhere vo - 0 and each fl, 0 is chosen so that II v, II 1 (i 0). The situation afterk steps is summarized byB V k Vk Tk flk lvk le T(2.2)where Tk - tridiag(fl, a,, fl, l) and Vk [Vl, v2. . . . . vk]. If there were no roundingerror we would have V T Vk I, and the process would therefore terminate withfl, 0 for some k n. Some other stopping criterion is needed in practice, sinceilk , is unlikely to be negligible for any k. In any event, eq. (2.2) holds to withinmachine precision.Now suppose we wish to solve the symmetric system B x b. Multiplying (2.2)by an arbitrary k-vector yk, whose last element is /h, gives B V k y k - Vk Tkyk flk Vk k. Since Vk (fl el) b by definition, it follows that ifyk and xk are definedby the equationsThyk file1,(2.3)xk Vkyk,then we shall have B x k b lkflk lVk to working accuracy. Hence xk may betaken as the exact solution to a perturbed system and will solve the originalsystem whenever 7/kflk l is negligibly small.The above arguments are not complete, but they provide at least somemotivation for defining the sequence of vectors {Xk} according to (2.3). It is nowpossible to derive several iterative algorithms for solving B x b, each characterized by the manner in which yk is eliminated from (2.3) (since it is usually notpractical to compute each Yk explicitly). In particular, the method of conjugategradients is known to be equivalent to using the Cholesky factorization TkLkDk L and is reliable when B (and hence Tk) is positive definite, while algorithmSYMMLQ employs the orthogonal factorization Tk --/:k Q to retain stability forarbitrary symmetric B. (See [20] for further details of these methods.)2.1 The Least-Squares SystemWe now turn to a particular system that arises from the "damped least-squares"problemmioll[:l]X-[:]ll:,2,,ACM Transactmns on Mathematmal Software, Vol. 8, No. 1, March 1982

46C.C. Paige and M. A. Saunderswhere A and b are given data and is an arbitrary real scalar. The solution of(2.4) satisfies the symmetric (but indefinite) systemArwhere r is the residual vector b - Ax. When the Lanczos process is applied to thissystem, a certain structure appears in relations (2.1)-(2.3). In particular, (2.1)reduces to the procedure defined as Bidiag 1 in Section 3, while (2.3) can bepermuted after 2k 1 iterations to the form['1P ,l Irk] ffi [ U I 0 lrt. ,l,(2.6)V,j L Y' Jx,where Bk is (k 1) x k and lower bidiagonal. We see that y is the solution ofanother damped least-squares problem,minl][B;]yk-[fl l][[2 ,(2.7)which can be solved reliably using orthogonal transformations. This observationforms the basis of algorithm LSQR.2.2 A Different Starting VectorFor completeness we note that a second least-squares algorithm can be derivedin an analogous way. Defining s -Ax, we can write (2.5) asAsand apply the Lanczos process again, using the same matrix as before but withthe new starting vector shown. This time, (2.1) reduces to Bidiag 2 as defined inSection 3, while (2.3) can be permuted to the form-x IJLy Jx,-O,e, '(2.9)V, JLyk j 'after 2k iterations, where Rk is k k and upper bidiagonal. (The odd-numberedsystems are not useful because T2k- is singular when 0.) We see that yksatisfies the system(RTRk 2I)yk ffi 01el,(2.10)which could easily be solved. However, (2.10) proves to be one form of the normalequations associated with (2.7), and the algorithm it suggests for computing xACM Transactlov.s on Mathematical Software, VoL 8, No. 1, March 1982.

LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares 47therefore has unsatisfactory numerical properties. We clarify this m a t t e r inSection 7.4.2.3 The Role ofThe quantities generated by the Lanczos process from (2.5) and (2.8) are Bk, Uk ,Vk and R k , P k , Vk, respectively. These are all i n d e p e n d e n t o f k , which meanst h a t they are the same as those generated when k ffi 0. We shall therefore assumefrom now on t h a t k ffi 0. A given k can be accommodated when solving (2.7), asshown in [22]. Methods for actually choosing ?, are beyond the scope of thispaper.3. THE BIDIAGONALIZATION PROCEDURESThe preceding use of the Lanczos process results in two forms of a bidiagonalization procedure due to Golub and K a h a n [9]. We state these forms as proceduresBidiag 1 and 2, and then give some unexpected relationships between them.B i d i a g I (starting vector b; reduction to lower bidiagonal form):f l l u l b,a l v l A Tul.fl, lu, A v , - a,u,0 1V l ATu I ,i ffi 1, 2 . . . . .flt lVz )(3.1)The scalars a, 0 and fl, 0 are chosen so t h a t [[u,[I -- [[ viii ffi 1. With thedefinitionsU k - [ U l , U 2 , . ,Uk],Vk [Vl, V2, ' ' , V, ],B2Bk --012[33 "".' . Otk'( w h e r e Bk is the rectangular matrix introduced in Section 2), the recurrencerelations (3.1) m a y be rewritten as(3.2)Uk (fl e ) ffi b,A Vk - Uk lBk,(3.3)TATUk I - V k B T Olk lVk lek l.(3.4)If exact arithmetic were used, then we would also have UT Uk ffi I and VkT V k /, but, in any event, the previous equations hold to within machine precision.B i d i a g 2 (starting vector A T b ; reduction to upper bidiagonal form):Oavl ffi A T b ,p pl ffi A v l .O, lV, ffi A T p , -- p,v,p, lp, l Av, l -- 0, 1p, ,i ffi 1, 2 . . . . .(3.5)JACM Transactions on Mathematmal Software, VoL 8, No 1, March 1982.

48 C.C. Paige and M. A. SaundersAgain, p, 0 and O, 0 are chosen so that IIP, II ]1vi II 1. In this case, ifplp2P k -- [ p l , p2, . . . ,pk],Yk --- I v , , v2, . . . ,Rk 03QQmvk],"'opk-101,we may rewrite (3.5) as(3.6)V k ( O l e l ) ffi A Tb,(3.7)A Vk ffi P k R k ,TA T p k ffi V,k R kT 0 k lVk lek,(3.8)and with exact arithmetic we would also have P T P k V kTVk LBidiag 2 is the procedure originally given by Golub and Kahan (with theparticular starting vector A Tb ). Either procedure may be derived from the otherby choosing.the appropriate starting vector and interchanging A and A T.3.1 Relationship Between the BidiagonalizationsThe principal connection between the two bidiagonalization procedures is thatthe matrices Vk are the same for each, and that the identityBTBkffi(3.9)RTRkholds. This follows from the fact that v is the same in both eases, and Vk ismathematically the result of applying k steps of the Lanezos process (2.2) withB A TA. The rather surprising conclusion is that Rk must be identical to thematrix that would be obtained from the conventional QR faetorization of Bk.Thuswhere Qk is orthogonal. In the presence of rounding errors, these identities will,of course, cease to hold. However, they throw light on the advantages of algorithmLSQR over two earlier methods, LSCG and L S L Q , as discussed in Section 7.4.The relationship between the orthonormal matrices Uh and Pk can be shownto be(3.11)for some vector rk. We also have the identities0/12 '4" 2 .p2,0 1 1 0l,(3.12)for i 1.ACM TransacUons on Mathematlcal Software, Vol 8, No. 1, March 1982

LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares 494. ALGORITHM LSQRThe quantities generated from A and b by Bidiag 1 will now be used to solve theleast-squares problem, min II b - A x II.Let the quantitiesxk V k y k ,(4.1)rk b - A x k ,(4.2)tk l ffi file1 - B k y k ,(4.3)be defined in terms of some vector yk. It readily follows from (3.2) and (3.3) t h a tthe equation(4.4)rk Uk ltk lholds to working accuracy. Since we want Hrk [[ to be small, and since Uk l isbounded and theoretically orthonormal, this immediately suggests choosing yk tominimize Utk H. Hence we are led naturally to the least-squares problemminII/ , el-- B k y k II(4.5)which forms the basis for LSQR.Computationally, it is advantageous to solve (4.5) using the standard QRfactorization of Bk [8], t h a t is, the same factorization (3.10) t h a t links the twobidiagonalizations. This takes the formpl 82p2 203t,Qo (4.6)pk-1 84Pk k- where Q k Qk, k l . . . Q2,3QI,2 is a product of plane rotations (e.g., [25]) designedto eliminate the subdiagonals f12, fla. . . . of Bk. The vectors yk and tk l could thenbe found fromRkyk h,tk Q [[o] h "(4.7)(4.8)However, yk in (4.7) will normally have r/o elements in common with yk-1. Insteadwe note t h a t [Rk f ] is the same as [Rk-1 fk-1] with a new row and columnadded. Hence, one way of combining (4.1) and (4.7) efficiently is according toxk V k R l f k -- D k f k ,where the columns of Dk -- [all(4.9)d2 . - . dk] can be found successively from theA C M T r a n s a c t i o n s on M a t h e m a t m a l Software, Vol. 8, N o 1, M a r c h 1982

50C . C . Paige and M. A. Saunders systemR kTD k T V T by forward substitution. With do x0 0, this givesdk ! ( v , - Okdk-1),(4.10)xk Xk-I kdk,(4.11)pkand only the most recent iterates need be saved. T h e broad outline of algorithmL S Q R is now complete.4.1 Recurrence RelationsThe QR factorization (4.6) is determined by constructing the k t h plane rotationQk,k l to operate on rows k and k 1 of the transformed [Bk / le ] to annihilatefl l. This gives the following simple recurrence relation:s,-c, JLP, lo , 1t'] o,, , ,,k] , , -k , '(4.12)where #51 --- al, 1 -- ill, and the scalars ca and sk are the nontrivial elements ofQk,k J. T h e quantities #54, k are intermediate scalars t h a t are subsequentlyreplaced by pk, d .T h e rotations Qk.k I are discarded as soon as t h e y have been used in (4.12),since Q, itself is not required. We see t h a t negligible work is involved in computingthe QR factorization to obtain Rk, fk, and k l.Some of the work in (4.10) can be eliminated by using vectors Wk -- pkdk inplace of dk. T h e main steps of L S Q R can now be summarized as follows. (Asusual the scalars a, 0 and fl, 0 are chosen to normalize the correspondingvectors; for example, al vl ffi A T e 1 implies the computations 61 ATul, a] II 1 1II,vl ( 1 / a l ) 1 . )Algortthm L S Q R(1)(Initialize.)fl u b,a v - - - - A Tul,Wl Vl,Xo O,(2) For i -- 1, 2, 3. . . . repeat steps 3-6.(3) (Continue the bidiagonalization.)(a)/ , 1u, 1 ffi A t , ,-a,u,(b) a, lv, l ATu, I --/ , v,.(4) (Construct and apply next orthogonal transformation.)(a)p, (g ( b ) c, F,/p,(c) s, #, ,tp,P,%,)'/ (d) O, l s,a, l(e) E -c,., (f) , c, ,(g) , , s, ,.ACM Transactions on Mathematical Software, Vol 8, No 1, March 1982

LSQR: An Algorithm for Sparse Linear Equat,ons and Sparse Least Squares 51(5) (Update x, w.)(a) x, x,-I ( ,lp,)w,(b) w, , v, - (8, l/p,)w,.(6) (Test for convergence.)Exit if some stopping criteria (yet to be discussed) have been met.4.2 Historical NoteThe previous algorithm was derived by the authors in 1973. An independentderivation is included in the work of Van Heijst et al. [24].5. ESTIMATION OF NORMSHere we show that estimates of the quantities II rk II, IIATrkll, Ilxkll, IIA II, andcond(A) can be obtained at minimal cost from items already required by LSQR.All five quantities will be used later to formulate stopping rules.Knowledge of [[A H and perhaps cond(A) can also provide useful debugginginformation. For example, a user must define his matrix A by providing twosubroutines to compute products of the form A v and A T u . These subroutines willtypically use data derived from earlier computations, and may employ rathercomplex data structures in order to take advantage of the sparsity of A. If theestimates of IIA I] and/or cond(A) prove to be unexpectedly high or low, then atleast one of the subroutines is likely to be incorrect. As a rule of thumb, werecommend that all columns of A be scaled to have unit length (11Aej]1 ffi 1, j 1 , . . . , n), since this usually removes some unnecessary ill-conditioning from theproblem. Under these circumstances, a programming error should be suspectedif the estimate of [[A Hdiffers by a significant factor from n 1/2 (since the particularnorm estimated will be IIA NF).For the purposes of estimating norms, we shall often assume that the orthogonality relations U k T U k I and v k T V k I hold, and that I[ Uk 112-- HVk 112--- 1. Inactual computations these are rarely true, but the resulting estimates have provedto be remarkably reliable.5.1 Estimates ofHrk I] and HA Trk IIFrom (4.4) and (4.8) we have--T(5.1)rk qJk l U k l Q k ek l(which explains the use of rk in (3.11)); hence, by assuming U T I U k I ffi L weobtain the estimateIIr, II Sk , ffi ,SkSk-1" ' 'Sl,(5.2)where the form of k fOllOWSfrom (4.12). LSQR is unusual in not having theresidual vector rk explicitly present, but we see that IIrk II is available essentiallyfree. Clearly the product of sines in (5.2) decreases monotonically. It shouldconverge to zero if the system A x ffi b is compatible. Otherwise it will converge toa positive finite limit.For least-squares problems a more important quantity is ATrk, which would bezero at the final iteration if exact arithmetic were performed. From (5.1), (3.4),A C M T r a n s a c t t o n s on M a t h e m a t i c a l Software, VoL 8, No. 1, M a r c h 1982

52 C.C. Paige and M. A. Saundersand (4.6) we haveATr4 4 l(V4B T ak lV4 l e T4 l ) Q kTe 4 l [?k ,V4[R T 0]e4 , 4 lak l(eW lQWe4 l)V4 l.The first term vanishes, and it is easily seen that the (k 1)st diagonal of Q4 is-c4. Hence we haveA Wr4 -- ( 4 l O14 l Ck )O4 l(5.3)IIATrkll k l k ,l ckl(5.4)andto working accuracy. No orthogonality assumptions are needed here.5.2 An Estimate of Ilxk IIThe upper bidiagonal matrix Rk may be reduced to lower bidiagonal form by theorthogonal factorizationR 4( kw 4,(5.5)where Q4 is a suitable product of plane rotations. Defining 5k by the system/ 4 4 -- fk,(5.6)it follows that xk --- ( V4R[1 )f4 ( Vk0 T) 5k -------1Y 454. Hence, under the assumptionthat v T v k ----I, we can obtain the estimateII x, II II e, II.(5.7)Note that the leading parts of L4, 4, l r4, and Ek do not change after iteration k.Hence we find that estimating II xk II via (5.5)-(5.7) costs only 13 multiplicationsper iteration, which is negligible for large n.5.3 Estimation of II A lit and cond(A)It is clear from (3.1) that all the v, lie in the range of A T and are thereforeorthogonal to the null space of A and A TA. With appropriate orthogonalityassumptions we have from (3.3) thatBTBk VTATAVk,and so from the Courant-Fischer minimax theorem the eigenvalues of BTBk areinterlaced by those of A T A and are bounded above and below by the largest andsmallest nonzero eigenvalues at A TA. The same can therefore be said of thesingular values of Bk compared with those of A. It follows that for the 2- and Fnorms,II Bk II - II A II,(5.8)where equality will be obtained in the 2-norm for some k rank(A) if b is notorthogonal to the left-hand singular vector of A corresponding to its largestsingular value. Equality will only be obtained for the F-norm if b containscomponents of all left-hand singular vectors of A corresponding to nonzeroACM Transactions on Mathematmal Software, Vol 8, No 1, March 1982.

LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares 53singular values. Nevertheless, we will use HBk [IF as a monotonically increasingestimate of the size of A.T h e foregoing also implies t h a t B T B k R kTRk is nonsingular and for the 2- andF-norms[IR ; [[ I[S [[ [[A [[.(5.9)T h e remarks on equality are the same, except now "largest singular value" isreplaced by "smallest nonzero singular value."Combining these results with the definition Dk V k R 1 in (4.9) now gives1 IIB II IID II- IIAII IIA ll-- cond(A)(5.10)for the 2- and F-norms. Hence we take I[Bk IIFIIDk Hr as a monotonically increasingestimate of cond(A), which starts at the optimistic estimate HB1 NFI[D HF -- 1.Use of Frobenius norms allows the estimates to be accumulated cheaply, since[[B , [[ [[B,-1 [[ a fl and [[Dk I[ IID,-1 ][ [[dk lie. T h e individualterms in the sum [[ dk [[2 y,,,.1 6 k can be used further for estimating standarderrors, as we show next.5.4 Standard ErrorsIn regression problems with m n rank(A), the standard error in the i t hc o m p o n e n t of the true solution x is taken to be & wheres'z 11b - A x [I2 o,,(5.11)m--nand o, - e T ( A T A ) - e , is the ith diagonal element of (ATA) - . Now from (3.3) and(3.10) we have V T A T A V k R T R k , which with (4.9) givesD T A TADk LAssuming t h a t p r e m a t u r e termination does not occur, it follows t h a t with exact T r k r kT e ,.arithmetic D n D T (ATA) -l, and we can approximate the a , b y o} - e,Since D k D [ Dk DT I d k d T, we haveo , o , -' , ,o , - O,and the e are monotonically increasing estimates of the o,,.In the implementation of L S Q R we accumulate o!,k for each i, and upon termination at iteration k we set 1 max(m - n, 1) and o u t p u t the squareroots ofs!k 1[b - Axkl[ 2 o k )1as estimates of the s, in (5.11). T h e a c c u r a c y of these estimates cannot beguaranteed, especially if termination occurs early for some reason. However, wehave obtained one reassuring comparison with the statistical package G L I M [16].On a m o d e r a t e l y ill-conditioned problem of dimensions 171 b y 38 (cond(A) 103, relative machine precision 10-11), an accurate solution xk was obtainedafter 69 iterations, and at this stage all s! k agreed to at least one digit with the s,o u t p u t by GLIM, and m a n y components agreed more closely.A C M T r a n s a c t m n s on M a t h e m a t m a l Software, Vol. 8, No. 1, M a r c h 1982

54 C . C . Paige and M. A. SaundersA further comparison was obtained from the 1033 by 434 gravity-meter problemdiscussed in [21]. For this problem a sparse QR factorization was constructed,Q A [0s], and the quantities a, were computed accurately using R T v i ei, ai, ffi[[ v, [[2. Again, the estimates of s}k) from LSQR proved to be accurate to at leastone significant figure, and the larger values were accurate to three or more digits.Note that s,2 estimates the variance of the ith component of x, and that s k) approximates this variance estimate. In an analogous manner, we could approximate certain covariance estimates by accumulating o,for any specific pairs (i,j), and then computing][ b - A x k [I 2 (, 1 qon termination. This facility has not been implemented in LSQR and we havenot investigated the accuracy of such approximations. Clearly only a limitednumber of pairs (i, j ) could be dealt with efficiently on large problems.6. STOPPING CRITERIAAn iterative algorithm must include rules for deciding whether the current iteratex is an acceptable approximation to the true solution x. Here we formulatestopping rules in terms of three dimensionless quantities ATOL, BTOL, andCONLIM, which the user will be required to specify. The first two rules apply tocompatible and incompatible systems, respectively. The third rule applies toboth. They are 1:Stop if I[r [[ BTOL[[ b [[ ATOI [ A [[ [I Xk U" 2:] 4 Trk [[ ATOL.Stop ifHA[[ []rk[[- 3:Stop if cond(A) CONLIM.We can implement these rules efficiently using the estimates of ][rk [[, [[A []F, andso forth, already described.The criteria 1 and 2 are based on allowable perturbations in the data. Theuser may therefore set ATOL and BTOL according to the accuracy of the data.For example, if (A, b) are the given data and (Z , ) represents the (unknown)true values, thenATOL [[A - [[IIA IIshould be used if an estimate of this is available; similarly for BTOL.Criterion 3 represents an attempt to regularize ill-conditioned systems.ACM Transactions on Mathematmal Software, Vol 8, No, 1, March 1982

LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares 556.1 Compatible SystemsTo justify S1, let r b - Axh as usual, and define the quantities64 - BTOL[I b II ATOLII AII [Ixk U,gk--BTOLII b IIrk - A T O M d II IIIIrkThen rk gk hk, and( A h Xx k ) Xk b - gkSO that xk is the exact solution for a system with both A and b perturbed. It canbe seen that these perturbations are within their allowable bounds when theinequality of S 1 holds. Hence, criterion S 1 is consistent with the ideas of backwardrounding error analysis and with knowledge of data accuracy. Since this argumentdoes not depend on orthogonality, 1 can be used in any method for solvingcompatible linear systems.6.2 Incompatible SystemsStewart [23] has observed that ifrk ffi b -- Axkandfk b - (A E k ) X kwhererkrWAiir, 'then (A E k ) T rk 0, SO that Xk and fk are the exact solution and residual for asystem with A perturbed. Since HEk 112 ]1A Trk I]/]] rk I], the perturbation to A willbe negligible if the test in 2 is satisfied.In our particular method it happens that Ekxk ffi 0, since (5.3) shows that x ffiVkyk is theoretically orthogonal to ATrk. Hence fk rk, SO both xk and rk areexact for the perturbed system. This strengthens the case for using rule 2.In practice we find that [IA Trk [I/I] rk I] can vary rather dramatically with k, butit does tend to stabilize for large k, and the stability is more apparent for LSQRthan for the standard method of conjugate gradients (see [[ATrk I] in Figures 3 and4, Section 8). Criterion 2 is sufficient to ensure that xk is an acceptable solutionto the least-squares problem, but the existence of an easily computable test thatis both sufficient and necessary remains an open question [23].6.3 Ill-Conditioned SystemsStopping rule S3 is a heuristic based on the following arguments. Suppose that Ahas singular values a o 2 - - .-. -- a . 0. It has been observed in some problemsA C M Transactions on M a t h e m a t i c a l Software, Vol. 8, No. 1, M a r c h 1982.

56C.C. Pmge and M A. Saunderst h a t as k increases the estimate ]]Bk lit HDk [Iv - cond(A) in (5.10) temporarilylevels off near some of the values of the ordered sequence ol/ol, ol/o2, . . . ,a l / o , , with varying numbers of iterations near each level. This tends to happenwhen the smaller o, are very close together, and therefore suggests criterion 3as a means of regularizing such problems when t h e y are very ill-conditioned, asin the discretization of ill-posed problems (e.g., [15]).For example, if the singular values of A were known to be of order 1, 0.9, 10 -3,10 -6, 10 -7, the effect of the two smallest singular values could probably besuppressed by setting C O N L I M 104.A more direct interpretation of rule 3 can be obtained from the fact t h a t xk D k f h . First, suppose t h a t the singular value decomposition of A is A U Z V Twhere u T u V T v W T I , - diag(al, o2, . . . , On), and letA(r) u ( r ) v Tbe defined by setting or . . . . .On ----0. A common m e t h o d for regularizing theleast-squares problem is to compute x (r) -- V ( r ) U T b for some r n, since it canreadily be shown t h a t the size of x (r) is bounded according toIIA [12[[x (r) 1[ cond(A(r)).II b IIorIn the case of LSQR, we have II xk II -- II Dk I1 II b II, and so if rule 3 has not yetcaused termination, we know t h a t II Bk IIF II x II / lib II --- II Bk IIF II D II CONLIM.Since UBk IIF usually increases to order UA [Iv quite early, we effectively haveIlAIIr Ilxkll CONLIM,II b IIwhich is exactly analogous to the bound above.6.4 Singular SystemsIt is sometimes the case t h a t rank(A) n. Known dependencies can often beeliminated in advance, but others m a y remain if only through errors in the dataor in formulation of the problem.With

LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares CHRISTOPHER C. PAIGE McGill University, Canada and MICHAEL A. SAUNDERS Stanford University An iterative method is given for solving Ax ffi b and minU Ax - b 112, where the matrix A is large and sparse.

Related Documents:

Bruksanvisning för bilstereo . Bruksanvisning for bilstereo . Instrukcja obsługi samochodowego odtwarzacza stereo . Operating Instructions for Car Stereo . 610-104 . SV . Bruksanvisning i original

approach based on a sparse mixture of sparse Gaussian graphical models (GGMs). Unlike existing fused- and group-lasso-based approaches, each task is represented by a sparse mixture of sparse GGMs, and can handle multi-modalities. We develop a variational inference algorithm combined with a novel sparse mixture weight selection algorithm. To handle

10 tips och tricks för att lyckas med ert sap-projekt 20 SAPSANYTT 2/2015 De flesta projektledare känner säkert till Cobb’s paradox. Martin Cobb verkade som CIO för sekretariatet för Treasury Board of Canada 1995 då han ställde frågan

service i Norge och Finland drivs inom ramen för ett enskilt företag (NRK. 1 och Yleisradio), fin ns det i Sverige tre: Ett för tv (Sveriges Television , SVT ), ett för radio (Sveriges Radio , SR ) och ett för utbildnings program (Sveriges Utbildningsradio, UR, vilket till följd av sin begränsade storlek inte återfinns bland de 25 största

Hotell För hotell anges de tre klasserna A/B, C och D. Det betyder att den "normala" standarden C är acceptabel men att motiven för en högre standard är starka. Ljudklass C motsvarar de tidigare normkraven för hotell, ljudklass A/B motsvarar kraven för moderna hotell med hög standard och ljudklass D kan användas vid

LÄS NOGGRANT FÖLJANDE VILLKOR FÖR APPLE DEVELOPER PROGRAM LICENCE . Apple Developer Program License Agreement Syfte Du vill använda Apple-mjukvara (enligt definitionen nedan) för att utveckla en eller flera Applikationer (enligt definitionen nedan) för Apple-märkta produkter. . Applikationer som utvecklas för iOS-produkter, Apple .

a key cost, and thereby a system performance bottleneck in many large SpMV computations. C. TAMU Sparse Matrix Collection The TAMU Sparse Matrix Suite Collection [5], is the largest, and the most diverse representation suite of sparse matrices available. It is an actively growing set of sparse matrices that arise in real applications.

2 The Adventures of Tom Sawyer. already through with his part of the work (picking up chips), for he was a quiet boy, and had no adventurous, troublesome ways. While Tom was eating his supper, and stealing sugar as opportunity offered, Aunt Polly asked him questions that were full of guile, and very deep for she wanted to trap him into damaging revealments. Like many other simple-hearted souls .