On short recurrences for generating orthogonal Krylov subspace bases Petr Tichý joint work with Vance Faber, Jörg Liesen, Zdeněk Strakoš Institute of Computer Science AS CR March 19, 2009, Podbanské, Slovakia Algoritmy 2009 1
Outline Introduction 1 Formulation of the problem 2 The Faber-Manteuffel theorem 3 4 Ideas of a new proof Barth-Manteuffel ( ℓ, m ) -recursion 5 6 Generating a B -orthogonal basis 2
Outline Introduction 1 Formulation of the problem 2 The Faber-Manteuffel theorem 3 4 Ideas of a new proof Barth-Manteuffel ( ℓ, m ) -recursion 5 6 Generating a B -orthogonal basis 3
Krylov subspace methods Basis Methods based on projection onto the Krylov subspaces K j ( A , v ) ≡ span( v, A v, . . . , A j − 1 v ) j = 1 , 2 , . . . . A ∈ R n × n , v ∈ R n . 4
Krylov subspace methods Basis Methods based on projection onto the Krylov subspaces K j ( A , v ) ≡ span( v, A v, . . . , A j − 1 v ) j = 1 , 2 , . . . . A ∈ R n × n , v ∈ R n . Each method must generate a basis of K j ( A , v ) . The trivial choice v, A v, . . . , A j − 1 v is computationally infeasible (recall the Power Method). For numerical stability: Well conditioned basis. For computational efficiency: Short recurrence. 4
Krylov subspace methods Basis Methods based on projection onto the Krylov subspaces K j ( A , v ) ≡ span( v, A v, . . . , A j − 1 v ) j = 1 , 2 , . . . . A ∈ R n × n , v ∈ R n . Each method must generate a basis of K j ( A , v ) . The trivial choice v, A v, . . . , A j − 1 v is computationally infeasible (recall the Power Method). For numerical stability: Well conditioned basis. For computational efficiency: Short recurrence. Best of both worlds: Orthogonal basis computed by short recurrence. 4
Optimal Krylov subspace methods with short recurrences CG (1952), MINRES, SYMMLQ (1975) based on three-term recurrences r j +1 = γ j A r j − α j r j − β j r j − 1 , 5
Optimal Krylov subspace methods with short recurrences CG (1952), MINRES, SYMMLQ (1975) based on three-term recurrences r j +1 = γ j A r j − α j r j − β j r j − 1 , generate orthogonal (or A -orthogonal) Krylov subspace basis, 5
Optimal Krylov subspace methods with short recurrences CG (1952), MINRES, SYMMLQ (1975) based on three-term recurrences r j +1 = γ j A r j − α j r j − β j r j − 1 , generate orthogonal (or A -orthogonal) Krylov subspace basis, optimal in the sense that they minimize some error norm: � x − x j � A in CG, � x − x j � A T A = � r j � in MINRES, � x − x j � in SYMMLQ -here x j ∈ x 0 + A K j ( A , r 0 ) . 5
Optimal Krylov subspace methods with short recurrences CG (1952), MINRES, SYMMLQ (1975) based on three-term recurrences r j +1 = γ j A r j − α j r j − β j r j − 1 , generate orthogonal (or A -orthogonal) Krylov subspace basis, optimal in the sense that they minimize some error norm: � x − x j � A in CG, � x − x j � A T A = � r j � in MINRES, � x − x j � in SYMMLQ -here x j ∈ x 0 + A K j ( A , r 0 ) . An important assumption on A : A is symmetric (MINRES, SYMMLQ) & pos. definite (CG). 5
Gene Golub By the end of the 1970s it was unknown if such methods existed also for general unsymmetric A . Gatlinburg VIII (now Householder VIII) held in Oxford from July 5 to 11, 1981. “A prize of $500 has been offered by Gene Golub for the construction of a 3-term conjugate gradient like descent method for non-symmetric real matrices or a proof that there G. H. Golub, 1932–2007 can be no such method”. 6
What kind of method Golub had in mind We want to solve A x = b using CG-like descent method: error is minimized in some given inner product norm, � · � B = �· , ·� 1 / 2 B . 7
What kind of method Golub had in mind We want to solve A x = b using CG-like descent method: error is minimized in some given inner product norm, � · � B = �· , ·� 1 / 2 B . Starting from x 0 , compute x j +1 = x j + α j p j , j = 0 , 1 , . . . , p j is a direction vector, α j is a scalar (to be determined), span { p 0 , . . . , p j } = K j +1 ( A , r 0 ) , r 0 = b − A x 0 . 7
What kind of method Golub had in mind We want to solve A x = b using CG-like descent method: error is minimized in some given inner product norm, � · � B = �· , ·� 1 / 2 B . Starting from x 0 , compute x j +1 = x j + α j p j , j = 0 , 1 , . . . , p j is a direction vector, α j is a scalar (to be determined), span { p 0 , . . . , p j } = K j +1 ( A , r 0 ) , r 0 = b − A x 0 . � x − x j +1 � B is minimal iff α j = � x − x j , p j � B � p j , p i � B = 0 . and � p j , p j � B 7
What kind of method Golub had in mind We want to solve A x = b using CG-like descent method: error is minimized in some given inner product norm, � · � B = �· , ·� 1 / 2 B . Starting from x 0 , compute x j +1 = x j + α j p j , j = 0 , 1 , . . . , p j is a direction vector, α j is a scalar (to be determined), span { p 0 , . . . , p j } = K j +1 ( A , r 0 ) , r 0 = b − A x 0 . � x − x j +1 � B is minimal iff α j = � x − x j , p j � B � p j , p i � B = 0 . and � p j , p j � B p 0 , . . . , p j has to be a B -orthogonal basis of K j +1 ( A , r 0 ) . 7
Faber and Manteuffel, 1984 Faber and Manteuffel gave the answer in 1984: For a general matrix A there exists no short recurrence for generating orthogonal Krylov subspace bases. What are the details of this statement ? 8
Outline Introduction 1 Formulation of the problem 2 The Faber-Manteuffel theorem 3 4 Ideas of a new proof Barth-Manteuffel ( ℓ, m ) -recursion 5 6 Generating a B -orthogonal basis 9
Formulation of the problem B -inner product, Input and Notation Without loss of generality, B = I . Otherwise change the basis: � x, y � B = � B 1 / 2 x, B 1 / 2 y � , A ≡ B 1 / 2 AB − 1 / 2 , ˆ v ≡ B 1 / 2 v . ˆ 10
Formulation of the problem B -inner product, Input and Notation Without loss of generality, B = I . Otherwise change the basis: � x, y � B = � B 1 / 2 x, B 1 / 2 y � , A ≡ B 1 / 2 AB − 1 / 2 , ˆ v ≡ B 1 / 2 v . ˆ Input data : A ∈ C n × n , a nonsingular matrix. v ∈ C n , an initial vector. 10
Formulation of the problem B -inner product, Input and Notation Without loss of generality, B = I . Otherwise change the basis: � x, y � B = � B 1 / 2 x, B 1 / 2 y � , A ≡ B 1 / 2 AB − 1 / 2 , ˆ v ≡ B 1 / 2 v . ˆ Input data : A ∈ C n × n , a nonsingular matrix. v ∈ C n , an initial vector. Notation: d min ( A ) . . . the degree of the minimal polynomial of A . d = d ( A , v ) . . . the grade of v with respect to A , the smallest d s.t. K d ( A , v ) is invariant under mult. with A . 10
Formulation of the problem Our Goal Generate a basis v 1 , . . . , v d of K d ( A , v ) s.t. 1. span { v 1 , . . . , v j } = K j ( A, v ) , for j = 1 , . . . , d , 2. � v i , v j � = 0 , for i � = j , i, j = 1 , . . . , d . 11
Formulation of the problem Our Goal Generate a basis v 1 , . . . , v d of K d ( A , v ) s.t. 1. span { v 1 , . . . , v j } = K j ( A, v ) , for j = 1 , . . . , d , 2. � v i , v j � = 0 , for i � = j , i, j = 1 , . . . , d . Arnoldi’s algorithm: Standard way for generating the orthogonal basis (no normalization for convenience): v 1 ≡ v , j � h i,j = � A v j , v i � v j +1 = A v j − h i,j v i , � v i , v i � , i =1 j = 0 , . . . , d − 1 . 11
Formulation of the problem Arnoldi’s algorithm - matrix formulation In matrix notation: v 1 = v , · · · h 1 , 1 h 1 ,d − 1 . ... . 1 . ... A [ v 1 , . . . , v d − 1 ] = [ v 1 , . . . , v d ] , h d − 1 ,d − 1 � �� � � �� � ≡ V d − 1 ≡ V d 1 � �� � ≡ H d,d − 1 V ∗ d V d is diagonal , d = dim K n ( A , v ) . 12
Formulation of the problem Arnoldi’s algorithm - matrix formulation In matrix notation: v 1 = v , · · · h 1 , 1 h 1 ,d − 1 . ... . 1 . ... A [ v 1 , . . . , v d − 1 ] = [ v 1 , . . . , v d ] , h d − 1 ,d − 1 � �� � � �� � ≡ V d − 1 ≡ V d 1 � �� � ≡ H d,d − 1 V ∗ d V d is diagonal , d = dim K n ( A , v ) . j � ( s + 2 )-term recurrence: v j +1 = A v j − h i,j v i . i = j − s 12
Formulation of the problem Optimal short recurrences (Definition - Liesen and Strakoš, 2008) A admits an optimal ( s + 2) -term recurrence, if for any v , H d,d − 1 is at most ( s + 2) -band Hessenberg, and for at least one v , H d,d − 1 is ( s + 2) -band Hessenberg. s + 1 � �� � • · · · • ... ... • ... ... • A V d − 1 = V d . ... ... . . ... • • � �� � d − 1 13
Recommend
More recommend
