Moments, Model Reduction and Nonlinearity in Solving Linear - PowerPoint PPT Presentation

Moments, Model Reduction and Nonlinearity in Solving Linear Algebraic Problems Zdenˇ ek Strakoš Charles University, Prague http://www.karlin.mff.cuni.cz/˜strakos 16th ILAS Meeting, Pisa, June 2010.

Thanks I wish to give my thanks to many authors and collaborators whose work and help has enabled the presented view. Related to moments, I feel more and more indebted to Gene Golub. Iveta Hnˇ etynková, Petr Tichý, Dianne O’Leary, Gerard Meurant. Z. Strakoš 2

Stieltjes (1893-1894) Given a sequence of numbers ξ k , k = 0 , 1 , . . . , a non-decreasing distribution function ω ( λ ) is sought such that � ∞ λ k dω ( λ ) = ξ k , k = 0 , 1 , . . . , 0 where ξ k represents the k -th (generalized) mechanical moment of the λ ≥ 0 . distribution of positive mass on the half line Key tool: continued fractions; cf. Stieltjes (1884, 1893-94), Chebyshev (1855, 1859, ... ). The story goes back at least to Gauss (1814) ..... Z. Strakoš 3

Matching moments formulation Consider a non-decreasing distribution function ω ( λ ) , λ ≥ 0 with the moments given by the Riemann-Stieltjes integral � ∞ λ k dω ( λ ) , ξ k = k = 0 , 1 , . . . . 0 λ ( n ) Find the distribution function ω ( n ) ( λ ) with n points of increase i i = 0 , 1 , . . . , which matches the first 2 n moments for the distribution function ω ( λ ) , � ∞ n λ k dω ( n ) ( λ ) ≡ ) k = ξ k , ω ( n ) ( λ ( n ) � k = 0 , 1 , . . . , 2 n − 1 . i i 0 i =1 Z. Strakoš 4

Outline 1. Matching moments model reduction 2. CG and Gauss-Christoffel quadrature 3. Non-Hermitian generalizations 4. Noise level in discrete ill-posed problems 5. Conclusions 5 Z. Strakoš

1 : Model reduction in dynamical systems Consider a linear dynamical system (here we assume A HPD) z ′ = A z + bu , b ∗ z . y = The transfer function b ∗ ( λI − A ) − 1 b given by the Laplace transform describes the impulse-response in the frequency domain. Z. Strakoš 6

ω ( λ ) 1 : Distribution function ω i = | ( s i , w 1 ) | 2 , w 1 = b/ � b � λ i , s i A , are the eigenpairs of ω N 1 . . . ω 4 ω 3 ω 2 ω 1 0 . . . . . . L λ 1 λ 2 λ 3 λ N U 7 Z. Strakoš

1 : Stieltjes recurrence Let p 0 ( λ ) ≡ 1 , p 1 ( λ ) , . . . , p n ( λ ) be the first n + 1 orthonormal ω ( λ ) . polynomials corresponding to the distribution function P n ( λ ) = [ p 0 ( λ ) , . . . , p n − 1 ( λ )] T , Then, writing λ P n ( λ ) = T n P n ( λ ) + δ n +1 p n ( λ ) e n represents the Stieltjes recurrence (1893-4), see Chebyshev (1855), Brouncker (1655), Wallis (1656), with the Jacobi matrix   γ 1 δ 2 ...   δ 2 γ 2   T n ≡ , δ l > 0 , ℓ = 2 , . . . , n .   ... ... δ n       δ n γ n Z. Strakoš 8

1 : Jacobi matrix and quadrature Consider the n th Gauss-Christoffel quadrature approximation ω ( n ) ( λ ) of ω ( λ ) . the Riemann-Stieltjes distribution function Its algebraic degree is 2 n − 1 , i.e., it matches the first 2 n moments � U n λ ℓ − 1 dω ( λ ) = ω ( n ) { λ ( k ) � j } ℓ − 1 , ξ ℓ − 1 = ℓ = 1 , . . . , 2 n , j L j =1 where , e 1 ) | 2 , ω ( n ) = | ( z ( n ) λ ( n ) = θ ( n ) , i i i i z ( n ) T n and is the normalized eigenvector of corresponding to the i θ ( n ) eigenvalue . The (orthogonal) polynomial p n ( λ ) has the roots i θ ( n ) , i = 1 , . . . , n . i Z. Strakoš 9

ω ( λ ) 1 : Continued fraction corresponding to 1 F N ( λ ) ≡ δ 2 2 λ − γ 1 − δ 2 3 λ − γ 2 − ... λ − γ 3 − . . . δ 2 N λ − γ N − 1 − λ − γ N The entries γ 1 , . . . , γ N and δ 2 , . . . , δ N represent coefficients of the Stieltjes recurrence. Z. Strakoš 10

1 : Partial fraction decomposition Consider (for simplicity of notation) � b � = 1 . Using the spectral decomposition, � U N dω ( µ ) ω j = R N ( λ ) � b ∗ ( λI − A ) − 1 b = λ − µ = P N ( λ ) , λ − λ j L j =1 R N ( λ ) P N ( λ ) ≡ F N ( λ ) The denominator P n ( λ ) corresponding to the n th convergent F n ( λ ) of F N ( λ ) , n = 1 , 2 , . . . is the n th orthogonal polynomial in the sequence determined by ω ( λ ) ; see Chebyshev (1855). Z. Strakoš 11

1 : Expansion at infinity Recall ( � b � = 1 ) � U N dω ( µ ) ω j � b ∗ ( λI − A ) − 1 b = λ − µ = ≡ F N ( λ ) , λ − λ j L j =1 and consider λ + µ 2 λ − µ = 1 1 1 − µ � − 1 = 1 � 1 + µ � � λ 2 + . . . , λ λ λ � − 1 � � λ + λ 2 � 1 = 1 1 − λ j = 1 1 + λ j j λ 2 + . . . . λ − λ j λ λ λ Z. Strakoš 12

1 : Minimal partial realization This gives the ninimal partial realization 2 n � � ξ ℓ − 1 1 � b ∗ ( λI − A ) − 1 b = F N ( λ ) = + O . λ ℓ λ 2 n +1 ℓ =1 Using the same expansion of the n th convergent F n ( λ ) of F N ( λ ) , 2 n � � ξ ℓ − 1 1 � F n ( λ ) = + O , λ ℓ λ 2 n +1 ℓ =1 Where the moments in the numerators are identical due to the Gauss quadrature property. Z. Strakoš 13

1 : Minimal partial realization 1859 - 94 The n th convergent F n ( λ ) of F N ( λ ) matches 2 n moments ξ 0 , . . . , ξ 2 n − 1 , F N ( λ ) and it approximates with the error proportional to λ − (2 n +1) . This represents the minimal partial realization; see Chebyshev (1859), Stieltjes (1894). The minimal partial realization was rediscovered in the engineering literature by Kalman (1979). The works of Krylov (1931), Hestenes and Stiefel (1952), Vorobyev (1958, 1965) (see Brezinski (1991, ... )), were not studied and recalled. The links with Chebyshev and Stieltjes were pointed out by Gragg (1974), Bultheel and Van Barel (1997). Z. Strakoš 14

Is it of any good to recall these line of thoughts in modern NLA? We wish to solve large systems of linear algebraic equations etc. by modern methods and algorithms and not to be bothered by some stories and the thoughts of Chebyshev or Stieltjes about some moments .... They can not be used in computations .... Are they of any use? Whatever we think, the moments are going to get us. Z. Strakoš 15

Outline 1. Matching moments model reduction 2. CG and Gauss-Christoffel quadrature 3. Non-Hermitian generalizations 4. Noise level in discrete ill-posed problems 5. Conclusions Z. Strakoš 16

λ = 0 2 : Specific problem for b ∗ ( λI − A ) − 1 b λ = 0 . Then the minimal partial Consider with realization representing the expansion at infinity is not applicable. We do either b ∗ A − 1 b = b ∗ ( A − 1 b ) ≈ b ∗ x n , x n A x = b , or with being an approximation to the solution of b ∗ A − 1 b ≈ b ∗ A − 1 n b , A n A . Mathematically, both with being an approximation to approaches are equivalent. Computationally, however, they can give different results. Z. Strakoš 17

2 : Caution A possible source of confusion with model reduction in dynamical systems. Approximation of the scalar value b ∗ A − 1 b will lead to matching the same moments as in the minimal partial realization given by the expansion of the function F N ( λ ) at infinity. A − 1 F N ( λ ) at zero that uses matching moments with The expansion of is in approximating the scalar b ∗ A − 1 b of no use. The error of the approximation will be expressed in an elegant way with no relation to the error estimates in dynamical systems model reduction (see below). Z. Strakoš 18

2 : Krylov subspace methods A x = b A n x n = b n x n approximates the solution x using the subspace of small dimension. S n ≡ K n ( A, r 0 ) ≡ span { r 0 , Ar 0 , · · · , A n − 1 r 0 } − → moments ! Z. Strakoš 19

2 : Conjugate gradients (CG): A HPD � x − x n � A = u ∈ x 0 + K n ( A,r 0 ) � x − u � A min w 1 = r 0 / � r 0 � , with the formulation via the Lanczos process, T n = W ∗ A W n = W n T n + δ n +1 w n +1 e T n , n ( A ) A W n ( A ) , and the CG approximation given by T n y n = � r 0 � e 1 , x n = x 0 + W n y n . A n = Q n A Q n = W n W ∗ n A W n W ∗ n = W n T n W ∗ n , Z. Strakoš 20

2 : Computational algorithm b ∗ A − 1 b we set x 0 = 0 ), r 0 = b − Ax 0 , Given x 0 (in approximating p 0 = r 0 For n = 1 , 2 , . . . γ n − 1 = ( r n − 1 , r n − 1 ) / ( p n − 1 , Ap n − 1 ) x n = x n − 1 + γ n − 1 p n − 1 r n − 1 − γ n − 1 Ap n − 1 r n = δ n = ( r n , r n ) / ( r n − 1 , r n − 1 ) p n = r n + δ n p n − 1 . Search directions are given by the modified residuals, γ n − 1 gives the line δ n ensures the local A -orthogonality of the direction search minimum, vectors. No moments are visible. If we wish to get an insight, we need them. Z. Strakoš 21

2 : Numerical PDE connection of CG Find u ≡ u ( ξ 1 , ξ 2 ) , where ξ 1 , ξ 2 denote the space variables, such that in a bounded domain Ω ⊂ R 2 , −∇ 2 u = f ∂u u = g D on ∂ Ω D , and ∂n = g N on ∂ Ω N , ∂ Ω D ∪ ∂ Ω N = ∂ Ω , ∂ Ω D ∩ ∂ Ω N = ∅ . where and For the Galerkin FEM approximation h ) � 2 = �∇ ( u − u h ) � 2 + � x − x n � 2 �∇ ( u − u ( n ) A . Z. Strakoš 22

An insight from viewing CG through moments? Z. Strakoš 23

CG ≡ matrix formulation of the Gauss q. 2 : � ξ ( λ ) − 1 dω ( λ ) ← → Ax = b , x 0 ζ ↑ ↑ n � − 1 � ω ( n ) θ ( n ) � T n y n = � r 0 � e 1 ← → i i i =1 x n = x 0 + W n y n ω ( n ) ( λ ) − → ω ( λ ) Z. Strakoš 24