An algorithm for the least-squares solution of rank-deficient linear - PowerPoint PPT Presentation

An algorithm for the least-squares solution of rank-deficient linear systems G. Rodriguez Dipartimento di Matematica e Informatica, Universit` a di Cagliari viale Merello 92, 09123 Cagliari, Italy in collab. with Antonio Aric` o (Dip. di Matematica e Informatica, Cagliari) Applied Inverse Problems 2009 Computational Aspects and Emerging Applications Vienna, July 20–24 2009 G. Rodriguez Least-squares solution of rank-deficient linear systems

Displacement structure The idea was originally conceived for Toeplitz matrix ∆ ξ,η ( T ) := Z ξ T − TZ η = G T H ∗ T T ∈ C m × n , G ∈ C m × δ , H ∈ C n × δ , and   0 0 · · · 0 φ   1 0 · · · 0 0     0 1 · · · 0 0 Z φ = .     . . . . . . . .   . . . . 0 0 · · · 1 0 [Kailath et al. 1979] G. Rodriguez Least-squares solution of rank-deficient linear systems

Displacement structure The idea was originally conceived for Toeplitz matrix ∆ ξ,η ( T ) := Z ξ T − TZ η = G T H ∗ T Toeplitz-like matrices form an algebra any Schur complement inherits the displacement structure it is possible to perform an LU factorization operating only on the generators with complexity O ( n 2 ) ( generalized Schur algorithm ) G. Rodriguez Least-squares solution of rank-deficient linear systems

Cauchy-like matrices & pivoting It was later noted that Cauchy-like structure is pivoting-invariant various classical structures (Toeplitz, Hankel, Vandermonde) can be efficiently converted into Cauchy-like [Heinig 1995], [Gohberg, Kailath, Olshevsky 1995] C ij = φ ∗ i · ψ j D t C − CD s = G C H ∗ , C t i − s j � � G ∗ D t = diag( t 1 , . . . , t m ) , C = · · · φ 1 φ m This leads to a fast and more stable algorithm, known as GKO, with complexity O ( n 2 ) and storage O ( n 2 ). G. Rodriguez Least-squares solution of rank-deficient linear systems

Schur complements of augmented matrices Various matrix problems can be reduced to the computation of Schur complements of suitable augmented matrices [Kailath, Chun 1994]. Tikhonov regularization J ( λ, x ) = � A x − b � 2 + λ 2 � H x � 2 , λ > 0   I m 0 A 0   0 I p λ H 0 S m + n + p   → ( A ∗ A + λ 2 H ∗ H ) − 1 A ∗ b = x λ − − − −   A ∗ λ H ∗ A ∗ b 0 0 0 I n 0 G. Rodriguez Least-squares solution of rank-deficient linear systems

Schur complements of augmented matrices Generalized Cross Validation (GCV) 1 m � b − A x λ � 2 A ( λ ) = A ( A ∗ A + λ 2 H ∗ H ) − 1 A ∗ V ( λ ) = 2 , [ 1 m trace( I − A ( λ )) ]   I m 0 A 0   0 I p λ H 0 S m + n + p   → A ( A ∗ A + λ 2 H ∗ H ) − 1 A ∗ − − − −  A ∗ λ H ∗ A ∗  0 0 0 A 0 G. Rodriguez Least-squares solution of rank-deficient linear systems

Schur complements of augmented matrices Multi-parameter regularization � � A x − b � 2 + λ 2 � H x � 2 + µ 2 � K x � 2 � J ( λ , x ) = 2 [Brezinski, Redivo-Zaglia, R, Seatzu Numer. Math. ]   I m 0 0 A 0   0 I p 0 λ H 0   S   → ( A ∗ A + λ 2 H ∗ H + µ 2 K ∗ K ) − 1 A ∗ b 0 0 I q µ K 0 −     A ∗ λ H ∗ µ K ∗ A ∗ b 0 0 0 0 I n 0 G. Rodriguez Least-squares solution of rank-deficient linear systems

Square linear systems Let C x = b , with C ∈ C n × n and b ∈ C n . Then, � C � b S n ( A C , b ) = C − 1 b . A C , b = ⇒ − I O In [Aric` o, R 2009] this method has been implemented in Matlab/C-MEX with O ( n 2 ) complexity and O ( n ) storage. This suite of programs allows to convert various structures into Cauchy-like ; works naturally with complex matrices; heavily based on the BLAS library; various pivoting strategies (partial, total, S&B’s, Gu’s); Toeplitz-like: no assumptions; Cauchy-like: reconstructability, can deal with multiple nodes. G. Rodriguez Least-squares solution of rank-deficient linear systems

Execution time 3 10 2 10 1 10 0 10 −1 10 −2 10 −3 10 partial Gu toms729 −4 10 LU α n 2 −5 10 128 256 512 1024 2048 4096 8192 16384 Random Toeplitz matrix - n = 2 k , k = 7 , . . . , 14 toms729 only works with real matrices drsolve always uses complex arrays & arithmetics G. Rodriguez Least-squares solution of rank-deficient linear systems

Errors w.r.to a prescribed solution −8 10 −9 10 −10 10 −11 10 −12 10 −13 10 no piv. partial total −14 10 Gu toms729 LU −15 10 128 256 512 1024 2048 4096 8192 16384 Random Toeplitz matrix - n = 2 k , k = 7 , . . . , 14 toms729 only works with real matrices drsolve always uses complex arrays & arithmetics G. Rodriguez Least-squares solution of rank-deficient linear systems

Least squares by Gaussian elimination A ∈ R m × n , A x = b , rank ( A ) = n ≤ m Peters-Wilkinson method [1970] � L 1 � L 1 , U ∈ R n × n PAQ = U , L 2 L T L y = L T ( P b ) , U ( Q T x ) = y [Bj¨ orck et al. 1988] contains an application to sparse matrices. G. Rodriguez Least-squares solution of rank-deficient linear systems

Least squares by Gaussian elimination A ∈ R m × n , A x = b , rank ( A ) = n ≤ m Peters-Wilkinson method [1970] Augmented matrix method [Siegel 1965, Bartel et al. 1970, Bj¨ orck 1992] � I � � r � � b � A = , r = b − A x A ∗ 0 x 0 Successfully applied in [Arioli et al. 1989] and [Duff et al. 1990] to sparse matrices, with Bunch-Kaufman pivoting. G. Rodriguez Least-squares solution of rank-deficient linear systems

Least squares via Schur complement This method was originally proposed for Toeplitz matrices in [Kailath et al. 1994]. It was further studied (and implemented) in [R 2006] for Toeplitz- and Cauchy-like matrices.   I m A 0  A ∗ A ∗ b  0 A x = b − → M A = 0 I n 0 S m + n ( M A ) = ( A ∗ A ) − 1 A ∗ b = x LS Since M A has displacement structure, it is possible to apply the generalized Schur algorithm. Moreover, the pseudoinverse matrix can be constructed S m + n ( M A ) = A † . b = I m ⇒ G. Rodriguez Least-squares solution of rank-deficient linear systems

Computation of Schur complement of M A Computation scheme the structure is converted from Toeplitz to generalized Cauchy, via fast transforms the generalized Schur algorithm is applied G. Rodriguez Least-squares solution of rank-deficient linear systems

Computation of Schur complement of M A Computation scheme the structure is converted from Toeplitz to generalized Cauchy, via fast transforms the generalized Schur algorithm is applied Advantages M A , has displacement structure, but it’s not Toeplitz-like, while M C is Cauchy-like Cauchy-like structure allows to apply partial pivoting fast algorithm O (27 r ( m + n ) 2 ) flops ( r = rank ∆ ( M A )) small storage required O (3 r ( m + n )) easily applicable to other kinds of displacement structure G. Rodriguez Least-squares solution of rank-deficient linear systems

Pivoting restriction While computing the Schur complement S r ( M ), the action of partial (or total) pivoting must be restricted to the first r rows (and columns) of M . In fact, if � M 11 � � P r � � Q r � M 12 0 0 M = , P = e Q = M 21 M 22 0 I m − r 0 I n − r then � P r M 11 Q r � P r M 12 PMQ = M 21 Q r M 22 and S r ( PMQ ) = S r ( M ) . G. Rodriguez Least-squares solution of rank-deficient linear systems

Partially reconstructable blocks C ij = φ ∗ i · ψ j t i − s j D t C − CD s = G C H ∗ C G. Rodriguez Least-squares solution of rank-deficient linear systems

Partially reconstructable blocks C ij = φ ∗ i · ψ j t i − s j D t C − CD s = G C H ∗ D L M C − M C D R = G M C H ∗ ⇒ C M C   I m C 0 D L = D t ⊕ D s ⊕ D s ,  ,  C ∗ C ∗ M C = 0 D R = D t ⊕ D s ⊕ D t , 0 I n 0 G. Rodriguez Least-squares solution of rank-deficient linear systems

Partially reconstructable blocks C ij = φ ∗ i · ψ j t i − s j D t C − CD s = G C H ∗ D L M C − M C D R = G M C H ∗ ⇒ C M C   I m C 0 D L = D t ⊕ D s ⊕ D s ,  ,  C ∗ C ∗ M C = 0 D R = D t ⊕ D s ⊕ D t , 0 I n 0 Some blocks are partially reconstructable: block (3 , 2) − → D s I n − I n D s = 0 It is necessary to store separately non recostructable entries and keep pivoting into account. G. Rodriguez Least-squares solution of rank-deficient linear systems

Optimal values of ξ and η Z ξ, m A − AZ η, n = G A H ∗ A , ξ, η ∈ C \ { 0 }   0 0 · · · 0 φ   1 0 · · · 0 0     0 1 · · · 0 0 ∈ C k × k . Z φ, k =    . . . .  . . . .   . . . . 0 0 · · · 1 0 G. Rodriguez Least-squares solution of rank-deficient linear systems

Optimal values of ξ and η Z ξ, m A − AZ η, n = G A H ∗ A , ξ, η ∈ C \ { 0 } D t C − CD s = G C H ∗ C C ij = φ ∗ i · ψ j t m = ξ, s n , j = η i t i − s j G. Rodriguez Least-squares solution of rank-deficient linear systems

Optimal values of ξ and η Z ξ, m A − AZ η, n = G A H ∗ A , ξ, η ∈ C \ { 0 } C ij = φ ∗ i · ψ j t m = ξ, s n , j = η i t i − s j Setting ξ = 1 e η = e i πϕ , let ϕ ∗ = arg max ϕ min i , j | t i − s j | . Then, i , j | t i − s j | = 2 sin πϕ ∗ ϕ ∗ = gcd( m , n ) , min 2 n . m G. Rodriguez Least-squares solution of rank-deficient linear systems