Block Bidiagonal Decomposition and Least Squares Problems ˚ Bjorck ¨ Ake Department of Mathematics ¨ Linkoping University Perspectives in Numerical Analysis, Helsinki, May 27–29, 2008 Outline Bidiagonal Decomposition and LSQR Block Bidiagonalization Methods Multiple Right-Hand Sides in LS and TLS Rank Deficient Blocks and Deflation Summary The Bidiagonal Decomposition In 1965 Golub and Kahan gave two algorithms for computing a bidiagonal decomposition (BD) L V T ∈ Rm×n , A=U 0 where B is lower bidiagonal, α1 β2 α2 . β3 . . B= .. . ∈ R(n+1)×n , αn βn+1 and U and V are square orthogonal matrices. Given u1 = Q1 e1 , the decomposition is uniquely determined if αk βk +1 6= 0, k = 1 : n, The Bidiagonal Decomposition The Householder Algorithm For given Householder matrix Q1 , set A0 = Q1 A, and apply Householder transformations alternately from right and left: Ak = Qk +1 (Ak −1 Pk ), k = 1, 2, . . . . Pk zeros n − k elements in kth row; Qk +1 zeros m − k − 1 elements in the k th column. If Q1 is chosen so that Q1 b = β1 e1 , then U T b AV = β1 e1 B , 0 0 where B is lower bidiagonal and U = Q1 Q2 · · · Qn+1 , V = P1 P2 · · · Pn−1 . The Bidiagonal Decomposition Lanczos algorithm In the Lanczos process, successive columns of U = (u1 , u2 , . . . , un+1 ), V = (v1 , v2 , . . . , vn ), with β1 u1 = b, β1 = kbk2 , are generated recursively. Equating columns in AT U = VB T and AV = UB, gives the coupled two-term recurrence relations (v0 ≡ 0): αk vk = AT uk − βk vk −1 , βk +1 uk+1 = Avk − αk uk , k = 1 : n. Here αk and βk +1 are determined by the condition kvk k2 = kuk k2 = 1. The Bidiagonal Decomposition Since the BD is unique, both algorithms generate the same decomposition in exact arithmetic. From the Lanczos recurrences it follows that Uk = (u1 , . . . , uk ) and Vk = (v1 , . . . , vk ), form orthonormal bases for the left and right Krylov subspaces h i Kk (AAT , u1 ) = span b, AAT b, . . . , (AAT )k−1 b , h i Kk (AT A, v1 ) = span AT b, . . . , (AT A)k −1 AT b . If either αk = 0 or βk+1 = 0, the process terminates. Then the maximal dimensioned Krylov space has been reached. Least Squares and LSQR The LSQR algorithm (Paige and Saunders 1982) is a Krylov subspace method for the LS problem min kb − Axk2 . x After k steps, one seeks an approximate solution xk = Vk yk ∈ Kk (ATA, AT b). From the recurrence formulas it follows that AVk = Uk +1 Bk . Thus b − AVk yk = Uk +1 (β1 e1 − Bk yk ), where Bk is lower bidiagonal. By the orthogonality of Uk +1 the optimal yk is obtained by solving a bidiagonal subproblem min kβ1 e1 − Bk y k2 , y Least Squares and LSQR This is solved using Givens rotations to transform Bk into upper bidiagonal form Rk fk G(β1 e1 | Bk ) = . φ¯k +1 0 giving xk = VRk−1 fk , kb − Axk k2 = |φ¯k +1 |. Then xk , k = 1, 2, 3, . . . gives approximations on a nested sequence of Krylov subspaces. Often superior to truncated SVD. LSQR uses the Lanczos recursion and interleaves the solution of the LS subproblems. A Householder implementation could be used for dense A. Least Squares and LSQR Bidiagonalization can also be used for total least squares (TLS) problems. Here the error function to be minimized is kb − Axk k22 kβ1 e1 − Bk yk k22 = . kxk k22 + 1 kyk k22 + 1 The solution of the lower bidiagonal TLS subproblems are obtained from the SVD of matrix β1 e1 Bk , k = 1, 2, 3 . . . . Partial least squares (PLS) is widely used in data analysis and statistics. PLS is mathematically equivalent to LSQR but differently implemented. The predictive variables tend to be many and are often highly correlated. Outline Bidiagonal Decomposition and LSQR Block Bidiagonalization Methods Multiple Right-Hand Sides in LS and TLS Rank Deficient Blocks and Deflation Summary Block Bidiagonalization Methods Let Q1 be a given m × m orthogonal matrix and set I U1 = Q1 p ∈ Rm×p . 0 Next form Q1T A and compute the LQ factorization of its first p rows (Ip 0)(Q1T A)P1 = (L1 0), where L1 is lower triangular. Proceed by alternately performing QR factorizations of blocks of p columns and LQ factorizations of blocks of p rows. Block Bidiagonalization Methods After k steps we have computed a block bidiagonal matrix L1 R2 L2 .. .. Tk = ∈ R(k+1)p×kp , . . Rk Lk Rk +1 and two block matrices with orthogonal columns I(k+1)p (U1 , U2 , . . . , Uk +1 ) = Qk +1 · · · Q2 Q1 , 0 (V1 , V2 , . . . , Vk ) = Pk · · · P2 P1 Ikp . 0 The matrix Tk is a banded lower triangular matrix with (p + 1) nonzero diagonals. Block Bidiagonalization Methods The Householder algorithm gives a constructive proof of the existence of such a block bidiagonal decomposition. If the LQ and QR factorizations have full rank then the decomposition is uniquely determined by U1 . To derive a block Lanczos algorithm we use the identities A(V1 V2 · · · Vn ) = (U1 U2 · · · Un+1 )Tn T A (U1 U2 · · · Un+1 ) = (V1 V2 · · · Vn )TnT . Start by forming Z1 = AT U1 and compute its thin QR factorization Z1 = V1 LT1 . Block Bidiagonalization Methods For k = 1 : n, we then do: Compute the residuals Wk = AVk − Uk Lk ∈ Rm×p , and compute the QR factorization Wk = Uk +1 Rk +1 . Compute the residuals Zk +1 = AT Uk +1 − Vk RkT+1 ∈ Rn×p , and compute the QR factorization Zk +1 = Vk +1 LTk+1 ; Householder QR factorizations can be used to guarantee orthogonality within each block Uk and Vk . Block Bidiagonalization Methods The block Lanczos bidiagonalization was given by Golub, Luk, and Overton in 1981. From the Lanczos recurrence relations it follows by induction that for k = 1, 2, 3, . . ., the block algorithm generates orthogonal bases for the Krylov spaces span(U1 , . . . , Uk ) = Kk (AAT , U1 ), span(V1 , . . . , Vk ) = Kk (AT A, AT U1 ), The process can be continued as long as the residual matrices have full column rank. Rank deficiency will be considered later. Outline Bidiagonal Decomposition and LSQR Block Bidiagonalization Methods Multiple Right-Hand Sides in LS and TLS Rank Deficient Blocks and Deflation Summary Multiple Right-Hand Sides in LS and TLS In many applications one needs to solve least squares problems with multiple right-hand sides B = (b1 , . . . , bp ). B = (b1 , . . . , bd ), d ≥ 2. There are two possible approaches; • Select a seed right-hand side and use the Krylov subspace generated to start up the solution of the second, etc. • Use a block Krylov solver where all right-hand sides are treated simultaneously. Block Krylov methods: • use a much larger search space from the start • introduce matrix–matrix multiplies into the algorithm. Multiple Right-Hand Sides in LS and TLS A natural generalization of LSQR is to compute a sequence of approximate solutions of the form Xk = Vk Yk , k = 1, 2, 3, . . . , are determined. That is, Xk is restricted to lie in a Krylov subspace Kk (ATA, AT B). It follows that B − AXk = B − AVk Yk = Uk +1 (R1 E1 − Tk Yk ). Using the orthogonality of the columns of Uk +1 gives kB − AXk kF = kR1 E1 − Tk Yk kF . Block LSQR Hence kB − AXk kF is minimized for Xk ∈ R(Vk ) by taking Xk = Vk Yk , where Yk solves min kR1 E1 − Tk Yk kF . Yk The approximations Xk are the optimal solutions on the nested sequence of block Krylov subspaces. Kk (AT A, AT B), k = 1, 2, 3, . . . .. The block bidiagonal LS subproblems are solved by using orthogonal transformations to bring the lower triangular banded matrix Tk into upper triangular banded form. Multiple Right-Hand Sides in LS and TLS As in LSQR the solution can be interleaved with the block bidiagonalization. For example, (p = 2) × × × ∗ ∗ ∗ × ⊗ × ∗ ∗ ∗ ∗ ⊗ ⊗ × ∗ ∗ ⊗ ⊗ × ∗ ∗ ∗ ∗ ∗ ⊗ ⊗ × ∗ ∗ ⊗ × × × × × = C1 RX C2 0 ∗ ∗ × × × × × × × × × × × × × S1 . S2 At this step the approximate solution to min kAX − BkF , is X = Vk (RX−1 C1 ), kB − AX kF = kC2 kf . Multiple Right-Hand Sides in LS and TLS The block algorithm can be used also for TLS problems with multiple right-hand sides which is min k F E kF , (A + E)X = B + F , E, F This cannot, as the LS problem, be reduced to p separate problems. Using the orthogonal invariance, the error function to be minimized is kR1 E1 − TY k2F kB − AX k2F = . kX k2F + 1 kY k2F + 1 Multiple Right-Hand Sides in LS and TLS The TLS solution can be expressed in terms of the SVD T Σ1 V1 B A = U1 U2 , Σ2 V2T where Σ1 = diag(σ1 , . . . , σn ). The solution of the lower triangular block bidiagonal TLS problem can be constructed from the SVD of the block bidiagonal matrix. Outline Bidiagonal Decomposition and LSQR Block Bidiagonalization Methods Multiple Right-Hand Sides in LS and TLS Rank Deficient Blocks and Deflation Summary Rank Deficient Blocks and Deflation When rank deficient blocks occur the bidiagonalization must be modified. For example (p = 1), at a particular step β1 α1 β2 α2 β3 Q3 Q2 Q1 (b, AP1 P2 ) = × × ⊗ ⊗ ⊗ ⊗ × × × × ⊗ × × × × , If α3 = 0 or β4 = 0, the reduced matrix has into block diagonal form. Then the LS problem decomposes and the bidiagonalization can be terminated. Rank Deficient Blocks and Deflation If αk = 0, then UkT (b, AVk ) = β1 e1 Tk 0 0 0 Ak and the problem is reduced to min kβ1 e1 − Tk y k2 . y Then the right Krylov vectors AT b, (AT A)AT b, . . . , (AT A)k −1 AT b are linearly dependent. Rank Deficient Blocks and Deflation If βk +1 = 0, then UkT+1 (b, AVk ) = β1 e1 T˜k 0 0 0 Ak where T˜k = (Tk αk ek ) is square. Then the original system is consistent and the solution satisfies T˜k y = β1 e1 . Then the left Krylov vectors b, (AAT )b, . . . , (AAT )k b are linearly dependent. Rank Deficient Blocks and Deflation From well known properties of tridiagonal matrices it follows: • The matrix Tk has full column rank and its singular values are simple. • The right-hand side βe1 has nonzero components along each left singular vector of Tk . Paige and Strakoˇs call this a core subproblem and shows that it is minimally dimensioned. If we scale b := γb, then only β1 changes. Thus, we have a core subproblem for any weighted TLS problem. min k(γr , E)kF s.t. (A + E)x = b + r . Rank Deficient Blocks and Deflation How do we proceed in the block algorithm if a rank deficient block occurs in an LQ or QR factorization? The triangular factor then typically has the form (p = 4 and r = rank(Rk ) = 3): × × × × 0 × × ∈ Rr ×p . Rk = 0 × 0 (Recall that the factorizations are performed without pivoting). All following elements in this diagonal of Tk will be zero. The bandwidth and block size are reduced by one. Rank Deficient Blocks and Deflation Example p = 2: Assume that the first element in the diagonal of the block R2 becomes zero. × × × × × × × × × ⊗ × × ⊗ × ⊗ ⊗ × × × × × × × × × × × × × × × × × × × × . If a diagonal element (say in L3 ) becomes zero the problem decomposes. In general, the reduction can be terminated when the bandwidth has been reduced p times. Rank Deficient Blocks and Deflation Deflation is related to linear dependencies in the associated Krylov subspaces: If (AAT )k bj is linearly dependent on previous vectors, then all left and right Krylov vectors of the form (AAT )p bj , (AT A)p AT bj , p ≥ k, should be removed. If (AT A)k AT bj is linearly dependent on previous vectors, then all vectors of the form (AT A)p AT bj , (AAT )p+1 bj , p ≥ k, should be removed. When the reduction terminates when both left and right Krylov subspaces have reached their maximal dimensions. Rank Deficient Blocks and Deflation The block Lanczos algorithm is modified similarly when a rank deficient block occurs. Recall the block Lanczos recurrence: For k = 1 : n, Vk LTk = AT Uk − Vk −1 RkT , Uk +1 Rk +1 = AVk − Uk Lk . Suppose the first rank deficient block is Lk ∈ Rp×r , r < p. Then Vk is n × r , i.e. only r < p vectors are determined. In the next step AVk − Uk Lk ∈ Rm×r . Recurrence still works! Block size has been reduced to r . Outline Bidiagonal Decomposition and LSQR Block Bidiagonalization Methods Multiple Right-Hand Sides in LS and TLS Rank Deficient Blocks and Deflation Summary • We have emphasized the similarity and (mathematical) equivalence of the Householder and Lanczos algorithms for block bidiagonalization. For dense problems in data analysis and statistics the Householder algorithm should be used because of its backward stability. Today quite large dense problems can be handled in seconds on a desktop computer! The relationship of block LSQR to multivariate PLS needs further investigation. G. H. Golub and W. Kahan. Calculating the singular values and pseudo-inverse of a matrix. SIAM J. Numer. Anal. Ser. B, 2:205–224, 1965. G. H. Golub, F. T. Luk, and M. L. Overton. A block Lanczos method for computing the singular values and corresponding singular vectors of a matrix. ACM Trans. Math. Software, 7:149–169, 1981. S. Karimi and F. Toutounian. The block least squares method for solving nonsymmetric linear systems with multiple right-hand sides. Appl. Math. Comput., 177:852–862, 2006. C. C. Paige and M. A. Saunders. LSQR. An algorithm for sparse linear equations and sparse least squares. ACM Trans. Math. Software, 8:43–71, 1982.

© Copyright 2019 ExploreDoc