A thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems
NNoname manuscript No. (will be inserted by the editor)
A thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems Ken-Ichi Ishikawa · Tomohiro Sogabe
Received: date / Accepted: date
Abstract
A thick-restart Lanczos type algorithm is proposed for Hermitian J -sym-metric matrices. Since Hermitian J -symmetric matrices possess doubly degeneratespectra or doubly multiple eigenvalues with a simple relation between the degenerateeigenvectors, we can improve the convergence of the Lanczos algorithm by restrictingthe search space of the Krylov subspace to that spanned by one of each pair of thedegenerate eigenvector pairs. We show that the Lanczos iteration is compatible with the J -symmetry, so that the subspace can be split into two subspaces that are orthogonalto each other. The proposed algorithm searches for eigenvectors in one of the twosubspaces without the multiplicity. The other eigenvectors paired to them can be easilyreconstructed with the simple relation from the J -symmetry. We test our algorithm onrandomly generated small dense matrices and a sparse large matrix originating from aquantum field theory. Keywords
Eigensolver · Lanczos method · Hermitian matrix · J -symmetric matrix Mathematics Subject Classification (2010) · · In theoretical physics, symmetry is an important guiding principle to search for lawsof nature. When numerically simulating a physics system that has a symmetry, it ispreferable to retain the symmetry property in the numerical simulation algorithm.
K.-I. IshikawaCore of Research for the Energetic Universe, Graduate School of Science, Hiroshima University, Higashi-Hiroshima 739-8526, JapanORCID: 0000-0001-6425-1894Tel.: +81-82-424-7363, Fax: +81-82-424-0717E-mail: [email protected]. SogabeDepartment of Applied Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, JapanE-mail: [email protected] a r X i v : . [ m a t h . NA ] S e p Ken-Ichi Ishikawa, Tomohiro Sogabe
Hermitian matrices often appear in analyzing physics systems, and spectra analysisis important to understand the nature. In this case the Hermitian symmetry mustbe considered in the spectra analysis and the eigensolver algorithm. Developingeigensolver algorithms for spectra analysis is one of the major research areas inapplied mathematics.In physics, matrices having both Hermiticity symmetry and J -symmetry exist. Let A be a Hermitian and J -symmetric matrix in C n × n . The J -symmetric property of A isdefined by JAJ − = A T , J T = − J , J T = J − , (1)where J is a skew-symmetric matrix in R n × n . We take this definition from [15] for the J -symmetry . Throughout this paper, we employ the following notations: T for matrixtransposition, H for Hermitian conjugation, and ∗ for complex conjugation. We use I and O to denote the identity and null matrices with appropriate sizes, respectively. Theeigenvalues of A are doubly degenerated, or doubly multiple, with the J -symmetry.Let x be an eigenvector of A associated to the eigenvalue λ , satisfying Ax = λ x . Thevector defined by y ≡ Jx ∗ , (2)is also the eigenvector of A associated to λ because of the symmetry (1).Any Hermitian J -symmetric matrix A has the following block structure. Withoutloss of generality, J is J = (cid:20) O − II O (cid:21) , (3)where the size of each block is ( n / ) × ( n / ) and n is an even number. Based on this,the symmetry (1) requires A = (cid:20) A A A H A T (cid:21) , (4)where A and A are ( n / ) × ( n / ) matrices satisfying A H = A and A T = − A .We can also explicitly construct eigenvectors x and y = Jx ∗ associated with an eigen-value λ in the block structure. Then A can be diagonalized as A = U (cid:20) Λ OO Λ (cid:21) U H , U = (cid:20) X − X ∗ X X ∗ (cid:21) , (5)where Λ = diag ( λ , λ , . . . , λ n / ) , and X and X are ( n / ) × ( n / ) matrices satisfying X H X + X H X = I and X T X − X T X = O .To investigate the spectrum of a sparse matrix in a large dimension, iterativeeigensolver algorithms are generally employed. Iterative eigensolver algorithms inthe literature, such as the Lanczos algorithm for Hermitian matrices, do not considerthe J -symmetry. Therefore, even after convergence on one of a pair of the degenerate The literature shows another definition, ( JA ) = ( JA ) T , for the J -symmetry; however, this is oppositefrom ours because JAJ − = − A T . thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems 3 eigenvector pairs, the algorithm continues to search for the other eigenvector asso-ciated to the converged eigenvalue even though the other eigenvector can be easilyreconstructed from the converged one. In this article, we restrict ourselves to improv-ing the Lanczos algorithm for large Hermitian J -symmetric matrices by incorporatingthe J -symmetry property.We could improve the convergence of the Lanczos algorithm by restricting thesearch space of the Krylov subspace of A to that spanned by one of each pair ofthe doubly degenerate eigenvector pairs by imposing complete orthogonality to thesubspace spanned by the other eigenvectors paired to them. However, this strategycannot be directly realized before knowing the invariant subspace of A . We find thatthe Krylov subspace K k ( A , v ) = span { v , Av , . . . , A k − v } generated with the Lanczosalgorithm with a starting vector v is orthogonal to K k ( A , Jv ∗ ) = span { Jv ∗ , AJv ∗ , . . . , A k − Jv ∗ } with the same Lanczos algorithm with the starting vector Jv ∗ . On the basisof this property, we can obtain only one half of the degenerate eigenvectors in K k ( A , v ) and can reconstruct the other half via (2). This can achieve the strategy just statedabove.Early attempts have been made to incorporate the symmetry property or structureof matrices into eigensolver algorithms in [7,15] in which dense matrix algorithmswere developed for the Hermitian J -symmetric matrices that appeared in a quantummechanical system with time-reversal and inversion symmetry. The dense matrixalgorithms retain or respect the matrix structure (4) during the computation [7,15].To focus on the algorithm for large sparse matrices, we do not discuss algorithmsfor dense matrices in this paper. The one similar to our algorithm that incorporatesthe symmetry properties of matrices was studied in [3,4,5,13], in which the Lanczostype eigensolvers for complex J -skew-symmetric matrices or Hamiltonian matriceswere investigated. Because the symmetry treated in [3,4,5,13] is different fromthe Hermitian J -symmetry, their algorithm cannot be directly applied to Hermitian J -symmetric matrices, even though it is based on the Lanczos algorithm. See also [12]for numerical algorithms for structured eigenvalue problems.This paper is organized as follows. In the next section, we prove the orthogonal-ity between the Krylov subspace K k ( A , v ) generated with the Lanczos iteration to K k ( A , Jv ∗ ) with the same Lanczos iteration. Having observed the orthogonality, wedescribe the restarting method based on the Krylov–Schur transformation method [19]and the thick-restart method [20,21]. Then, we propose the thick-restart Lanczos algo-rithm for Hermitian J -symmetric matrices in Section 3. We estimate the computationalcost in terms of the matrix-vector multiplication. In Section 4, we test the proposedalgorithm for two types of the Hermitian J -symmetric matrix. The one type is anartificial randomly generated matrix satisfying the structure of (4) and (5), and theother matrix originates from quantum field theory. The convergence behavior andthe computational cost are compared with those of the standard thick-restart Lanczosalgorithm. We summarize this paper in the last section. The definition of the J -symmetry employed in [3,4,5,13] is opposite to that of [7,15] and ours;therefore, their matrices are J -skew-symmetric compared to our definition. Ken-Ichi Ishikawa, Tomohiro Sogabe J -symmetry The Lanczos iteration transforms A to a tridiagonal form by generating the orthonormalbasis vectors. The Lanczos decomposition after m -step iteration starting with a unitvector v is given by AV m = V m T m + β m v m + e Tm , (6)where V m = [ v , v , . . . , v m ] and v j ∈ C n , and e m is the m -dimensional unit vector inthe m -th direction. The basis vectors are orthonormal, V Hm + V m + = I . T m denotes the m × m tridiagonal matrix that is given by T m = α β β α β β α β . . . . . . . . . β m − α m − β m − β m − α m . (7)The approximate eigenpairs are obtained from the eigenpairs of T m . Because of theHermiticity property of A and the recurrence structure of the Lanczos iteration, all α j and β j can be taken to be real. Thus, T m becomes a real symmetric matrix. In thestandard Lanczos iteration, the Hermitian symmetry of A is respected in the form T m .We also investigate the J -symmetry property of decomposition (6). After takingthe complex conjugate of (6) followed by multiplying it by J from the left-hand side,we have JA ∗ V ∗ m = JV ∗ m T m + β m Jv ∗ m + e Tm . (8)Using the J -symmetry and Hermiticity, JA ∗ = A H J = AJ , and by defining w j ≡ Jv ∗ j ,we obtain AW m = W m T m + β m w m + e Tm , (9) W m = [ w , w , . . . , w m ] . (10)The columns of W m + are also orthonormal, W Hm + W m + = I . Consequently, the vectors W m + have the same Lanczos decomposition as that of V m + when it starts from w = Jv ∗ . Furthermore, we can prove the orthogonality between W m + and V m + .To implement the thick-restart method [20,21], we show the J -symmetry propertyfor a generalized decomposition similar to the Krylov–Schur decomposition [19] in-stead of the Lanczos decomposition in the following part of the paper. After preparingtwo lemmas related to the J -symmetry, we prove the main theorem for the orthogonal-ity property between W m + and V m + on the generalized decomposition. Subsequently,the orthogonality properties for the Lanczos decomposition and the thick-restartmethod are proved as corollaries.We first show orthogonal properties between v and w = Jv ∗ . thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems 5 Lemma 1
Let v be an arbitrary vector in C n , and J and A be matrices satisfying (1) .Then, the vector w defined by w ≡ Jv ∗ satisfiesw H v = , (11) w H Av = . (12) Proof
From the definition of w , it follows that w H v = ( Jv ∗ ) H v = v T J H v = v T J T v = − v T Jv , (13)where J T = − J is used. The identity v T Jv = v T J T v and (13) yield w H v = − w H v . (14)Thus, w H v =
0. Similarly, w H Av = ( Jv ∗ ) H Av = v T J H Av = v T J T Av = − v T JAv = − v T A T Jv , (15)where J T = − J and JA = A T J are used. The identity v T A T Jv = v T J T Av and (15) yield w H Av = − w H Av . (16)Thus, w H Av = (cid:117)(cid:116) We have the following lemma that is a generalization of the relation between (6)and (9).
Lemma 2
Let v , v , . . . , v k , v k + be the vectors having the following relation:AV k = V k S k + v k + b T , (17) where V k = [ v , . . . , v k ] , S k is a matrix in R k × k , and b is a vector in R k for a HermitianJ-symmetric matrix A satisfying (1) . Then, the vectors w j = Jv ∗ j ( j = , . . . , k + ) satisfy the following decomposition:AW k = W k S k + w k + b T , (18) where W k = JV ∗ k = [ w , . . . , w k ] .Proof By taking the complex conjugate of (17) followed by multiplying it by J fromthe left-hand side, we obtain (18) using JA T = AJ and A ∗ = A T . (cid:117)(cid:116) Using Lemmas 1 and 2, we can show the orthogonality properties between W k + and V k + as follows. Theorem 1
Let V k + = [ v , . . . , v k + ] be a matrix satisfying (17) , and W k + = JV ∗ k + . Ifthe matrices V k and W k are orthogonal and A-orthogonal to each other: V Hk W k = O andV Hk AW k = O, then the matrices V k + and W k + also satisfy the following orthogonalityrelations: V Hk + W k + = O , (19) V Hk + AW k + = O . (20) Ken-Ichi Ishikawa, Tomohiro Sogabe
Proof
The decomposition (18) follows from Lemma 2. Multiplying W Hk and V Hk to (17) and (18), respectively, yields W Hk AV k = W Hk V k S k + W Hk v k + b T , (21) V Hk AW k = V Hk W k S k + V Hk w k + b T . (22)From the premise that W Hk AV k = V Hk AW k = O and W Hk V k = V Hk W k = O , it follows that W Hk v k + = O , (23) V Hk w k + = O . (24)Together with Lemma 1 and the premise, we find V Hk + W k + = O . (25)Multiplying w Hk + and v Hk + to (17) and (18), respectively, yields w Hk + AV k = w Hk + V k S k + w Hk + v k + b T , (26) v Hk + AW k = v Hk + W k S k + v Hk + w k + b T . (27)Because of (23) and (24) as well as Lemma 1, the right-hand sides of (26) and (27)vanish. Thus, w Hk + AV k = O , (28) v Hk + AW k = O . (29)Together with Lemma 1 and the premise, we find V Hk + AW k + = O . (30)Therefore, V k + and W k + are orthogonal and A -orthogonal to each other. (cid:117)(cid:116) Using Theorem 1 as well as Lemmas 1 and 2, we can show that the Lanczosvectors V m + generated with (6) are orthogonal and A -orthogonal to W m + = JV ∗ m + . Corollary 1
The m-step Lanczos vectors V m + with m ≥ generated with a unitvector v for a Hermitian J-symmetric matrix A satisfyV Hm + W m + = O , V Hm + AW m + = O , (31) with W m + = JV ∗ m + , (32) when no breakdown occurs.Proof Because the Lanczos decomposition (6) is a particular form of (17) with k → m ,a real matrix S m → T m and a real vector b → β m e m , Lemma 2 can be applied toobtain (9). Because w H v = w H Av = m =
1. Moreover, the Lanczos decomposition retainsits form applicable to Theorem 1 for any m >
1. Therefore, the corollary follows fromTheorem 1 and Lemmas 1 and 2 by induction. (cid:117)(cid:116) thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems 7 We cannot simultaneously find degenerate pairs of eigenvectors with the standardsingle-vector Lanczos process. This is true for computation with exact arithmetic.However, with finite precision arithmetic, the single-vector Lanczos process wouldgenerate a small overlap to W k via round-off errors, so that even after the convergenceof an eigenvector, a late convergence to the other paired eigenvector would be possible.Because this behavior is rather accidental, a block-type Lanczos algorithm has to beapplied to accelerate the convergence for degenerate eigenvectors [1,2,8,18,22].We further investigate the structure of the Lanczos decomposition. According toCorollary 1 and the Lanczos decompositions (6) and (9), the block-type decompositioncan be constructed as: AV [ m ] = V [ m ] T [ m ] + V m + B m E T [ m ] , (33)where we define V [ m ] ≡ (cid:2) V V . . . V m − V m (cid:3) , T [ m ] ≡ A B B A B . . . . . . . . . B m − A m − B m − B m − A m , V j ≡ (cid:2) v j , w j (cid:3) , A j ≡ diag ( α j , α j ) , B j ≡ diag ( β j , β j ) , E T [ m ] ≡ (cid:2) O O . . .
O I (cid:3) , (34)where the size of E T [ m ] is 2 × m . When m = n /
2, the decomposition should termi-nate, because V [ n / ] completely block-tridiagonalize A and the Krylov subspaces K n / ( A , v ) and K n / ( A , Jv ∗ ) span the entire eigenspace of A . Because K n / ( A , v ) and K n / ( A , Jv ∗ ) are orthonormal and (6) and (9) are independent iterations, the Lanczositeration for (6) terminates at m = n / K n / ( A , v ) without the eigenvalue multiplicity associ-ated with the J -symmetry. In other words, the standard Lanczos iteration with exactprecision arithmetic is enough to find all the eigenvectors without multiplicity. How-ever, this is impractical because it requires exact precision arithmetic. With the finiteprecision, the orthogonality to K n / ( A , Jv ∗ ) is not maintained because of round-offerrors and, eventually, eigenvectors, including multiplicity, could be extracted fromthe single-vector Lanczos iteration, as stated previously.By using Corollary 1 and with the above analysis, we can construct a Lanczostype algorithm in which the orthogonality to W k is enforced to search for eigenvectorswithout multiplicity of eigenvalues associated to the J -symmetry. Additionally, theother vectors paired to them can be easily reconstructed. However, for a practicalnumerical algorithm of the Lanczos type iteration, the iteration should terminateat a finite step, and a restarting mechanism is required [6,20,21]. The most usefuland simplest but effective restarting method is the so-called thick-restart method [20,21] that is a specialization of the Krylov–Schur transformation [19] to Hermitianmatrices. To involve the thick-restart method to the Lanczos algorithm with the J -symmetry, we have to prove the orthogonality between V k and W k after the Krylov–Schur transformation and restarting. To achieve this, we have the following corollary. Ken-Ichi Ishikawa, Tomohiro Sogabe
Corollary 2
Let V m + and W m + be the orthonormal matrices containing basis vec-tors generated with the m-step Lanczos process, (6) and (9) , and Z k be an orthonormalmatrix in R k × k . The Krylov–Schur transformation with Z k on the Lanczos decomposi-tion is defined by AU m = U m ( Z − m T m Z m ) + v m + b Tm , (35) AQ m = Q m ( Z − m T m Z m ) + w m + b Tm , (36) U m ≡ V m Z m , (37) Q m ≡ W m Z m , (38) b Tm ≡ β m e Tm Z m . (39) Then, the matrices U m + = [ U m , v m + ] and Q m + = [ Q m , w m + ] satisfy the orthogonalrelations: U Hm + Q m + = O and U Hm + AQ m + = O.Proof
Because U m and Q m satisfy U Hm Q m = O and Q Hm AU m = O and the decompo-sitions (35) and (36) are particular forms of the decomposition in Lemma 2, theorthogonality and A -orthogonality between U m + and Q m + simply follow from Theo-rem 1. (cid:117)(cid:116) For the thick-restart method, Z m is chosen to diagonalize T m , and the dimensionof the decomposition is reduced from m to k < m with a selection criterion forvectors V k ← V m . In the reduction, the last vectors are kept hold as v k + ← v m + and w k + ← w m + to retain the decomposition form properly. The orthogonality propertiesof the new basis ( U m + , Q m + ) and the reduced basis ( U k + , Q k + ) still hold accordingto Corollary 2. After restarting, the Lanczos iteration continues to keep the decomposedform applicable to Theorem 1. Because the structures of the Lanczos and Krylov–Schur decompositions for W k and Q k are the same as those for V k and U k , respectively,we do not need to explicitly iterate the Lanczos algorithm for W k and Q k . Consequently,we can continue the Lanczos thick-restart cycle only in the subspace K k ( A , v ) that isorthogonal and A -orthogonal to K k ( A , Jv ∗ ) . We note that according to Corollaries 1and 2, all the eigenpairs without multiplicity can be obtained with the thick-restartLanczos algorithm using exact precision arithmetic. This is impractical and we mustenforce the orthogonality between K k ( A , v ) and K k ( A , Jv ∗ ) for a practical algorithm. J -symmetry Based on Theorem 1 as well as Corollaries 1 and 2, we construct a thick-restartLanczos algorithm for Hermitian J -symmetric matrices (TRLAN–JSYM) which effi-ciently searches for eigenvectors without the multiplicity of eigenvalues in K k ( A , v ) .Algorithm 1 shows the TRLAN–JSYM algorithm. The Lanczos iteration with the J -symmetry is described in Algorithm 2. We include the invert mode for the smalleigenvalues.The main difference from the standard thick-restart Lanczos algorithm (TRLAN)is in Algorithm 2, where we simultaneously construct w j using w j = Jv ∗ j and enforcethe orthogonality of V j + to W j = JV ∗ j to avoid the contamination of the search space K k ( A , v ) from K k ( A , Jv ∗ ) . thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems 9 We will compare the efficiency between the TRLAN–JSYM and the standardTRLAN algorithms in Section 4. For the comparison, we estimate the computationalcost of the TRLAN–JSYM and the standard TRLAN algorithms as follows. We countthe total number of matrix-vector multiplication Av j or A − v j contained in the Lanczosstep Algorithm 2. For the invert mode with a large sparse matrix, it could requirean iterative linear solver and a computational cost to obtain A − v j . To focus on thecomputational cost comparison between the TRLAN–JSYM and TRLAN algorithms,assuming the computation cost of a single inversion is identical between the twoalgorithms, we do not count the cost involved in the inversion and regard the singleinversion operation A − v j as one matrix-vector multiplication for the invert mode. Weneglect the cost that explicitly computes the true residual at line 28 in Algorithm 1. Forthe first outer iteration, the count is m , and after restarting, it is m − k . The thickness k for restarting is defined by k = min ( icnv + mwin , m − ) , (40)as shown in line 48 of Algorithm 1, where icnv is the number of converged eigenvectorsand mwin is the initial thickness for restarting. When all desired eigenvectors areobtained at an outer iteration N conv , the upper and lower bounds of the total number ofmatrix-vector multiplication N MV is estimated as m + ( m − mwin − nev )( N conv − ) < N MV < m + ( m − mwin )( N conv − ) , (41)where nev is the number of desired eigenvectors without the multiplicity of J -symmetry.The inequality (41) follows from the fact that icnv increases monotonically from zeroto nev toward the convergence.The same cost estimate can be derived for the standard TRLAN algorithm. TheTRLAN algorithm can be obtained by removing W m from Algorithms 1 and 2. There-fore, the cost bound for the TRLAN algorithm is identical to (41). However, to findall eigenvectors paired with the J -symmetry using the TRLAN algorithm, nev for theTRLAN must be double that of the TRLAN–JSYM. Thus it is natural to double allparameters for the TRLAN algorithm than those of the TRLAN–JSYM. Therefore,the upper and lower bounds of the total number of matrix-vector multiplication N MV of the TRLAN algorithm is2 ( m + ( m − mwin − nev )( N (cid:48) conv − )) < N MV < ( m + ( m − mwin )( N (cid:48) conv − )) , (42)where the parameters ( nev , mwin , m ) are those of the TRLAN–JSYM, and we intro-duced N (cid:48) conv for the number of outer iterations because it could be different from thatof the TRLAN–JSYM.Although the computational cost of the Gram–Schmidt orthonormalization of theLanczos step is minor compared to that of the matrix-vector multiplication, we brieflydiscuss the cost here. As shown in Algorithm 2 for the TRLAN–JSYM, the number ofvectors to orthonormalize is 2 m and the cost scales with O (( m ) ) . On the other hand,it scales with O ( m ) to orthonormalize V m with the standard Lanczos algorithm. Asdescribed above, it is natural to double the parameters for the TRLAN. The cost of the Lanczos part to orthonormalize V m then becomes O (( m ) ) . Therefore, the scalingof the cost is the same for both algorithms.Our naive estimates on the computational cost are (41) and (42), where we as-sume that the parameters of the TRLAN is twice as large as those of the TRLAN–JSYM. Although N conv and N (cid:48) conv depend on ( nev , mwin , m ) and the algorithm itself,if N conv (cid:39) N (cid:48) conv holds, the TRLAN–JSYM algorithm has a better performance thanthe TRLAN algorithm. We will see whether the condition N conv (cid:39) N (cid:48) conv holds or notin the numerical tests on the two types of the matrices in the next section.So far, we have described the single-vector Lanczos iteration type algorithm tointroduce the TRLAN–JSYM. If the matrix A has a dense cluster of eigenvalues ormultiple eigenvalues other than those with the J -symmetry, we need to incorporatethe block type Lanczos iteration in the algorithm for efficiency. We can extend theproposed algorithm to the blocked version in the same manner as it was conducted forthe standard thick-restart Lanczos algorithm [18,22]. The study on the block versionwill be addressed in future studies and we have only shown the single vector versionto demonstrate the idea for simplicity. In this section, we show two numerical tests to explore the efficiency of the TRLAN–JSYM algorithm compared to the TRLAN algorithm for Hermitian J -symmetricmatrices. The first test is conducted for randomly generated Hermitian J -symmetricmatrices satisfying the structure (4). The second test is applied to a matrix in quantumfield theory. We refer to these two test cases as Case A and Case B, respectively.We implement both the algorithms, TRLAN–JSYM and TRLAN, with Fortran2003. The numerical tests were performed on a single node of the subsystem A of theITO supercomputer system of Kyushu university [16]. The code is parallelized usingOpenMP and the Intel MKL library and executed with 36 threads.4.1 Case A We generated ten matrices with a size of 2000 × A with the structure (4), we employ (5). Theeigenvalues Λ are generated from uniformly distributed random real numbers in { x ∈ R : 0 < x < } , and the elements of matrices X and X are constructed from uniformlydistributed random complex numbers in { z ∈ C : − < Re ( z ) < , − < Im ( z ) < } .The constraints X H X + X H X = I and X T X − X T X = O are then imposed by aGram–Schmidt algorithm similar to that used in Algorithm 2. Fig. 1 shows largeeigenvalues for the ten random matrices. We solve the largest several eigenvalues withthe TRLAN–JSYM and TRLAN algorithms and compare the convergence behavior. thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems 11 S a m p l enu m be r Eigenvalue
Fig. 1: Large eigenvalue distribution of the ten random matricesTable 1: Algorithm parameters and statistics for convergence (Case A) nev mwin m N conv N MV TRLAN–JSYM 5 10 50 [ 6, 9.1, 14] [ 248, 369.6, 559]100 [ 3, 4.1, 6] [ 280, 377.5, 546]150 [ 2, 2.8, 4] [ 290, 400.8, 567]200 [ 2, 2.1, 3] [ 389, 408.7, 578]10 20 50 [ 10, 12.6, 16] [ 306, 378.2, 484]100 [ 4, 4.9, 6] [ 332, 404.6, 497]150 [ 3, 3.2, 4] [ 401, 432.0, 537]200 [ 2, 2.3, 3] [ 379, 432.1, 557]TRLAN 10 20 100 [ 6, 8.3, 12] [ 480, 655.6, 935]200 [ 3, 3.8, 5] [ 555, 693.7, 904]300 [ 2, 2.6, 3] [ 572, 740.3, 856]400 [ 2, 2.1, 3] [ 773, 813.1,1150]20 40 100 [ 10, 12.6, 15] [ 569, 699.2, 844]200 [ 4, 4.6, 6] [ 652, 746.7, 958]300 [ 3, 3.1, 4] [ 794, 828.4,1051]400 [ 2, 2.1, 3] [ 748, 787.7,1110]
Table 1 shows the algorithmic parameters used in this test. For the stopping conditionof the algorithms, we employ tol = − for tolerance. We also tabulate the numberof outer iteration counts N conv and matrix-vector multiplications N MV for convergencein the table. The minimal, average, and maximal values from the ten samples areshown in square brackets, respectively.Figs. 2 and 3 are the convergence histories of the eigenvalues and correspondingresiduals for the 1st random matrix, respectively. The left panels show the resultwith the TRLAN–JSYM algorithm with ( nev , mwin , m ) = ( , , ) , and the rightpanels show the result with the TRLAN algorithm with ( , , ) . We employ thisdoubled parameter for the TRLAN, as discussed in the previous section. The behaviorof the TRLAN–JSYM is smooth, while it reorders several times for the TRLAN,even though the TRLAN algorithm successfully captures all the eigenvalues with E i ngen V a l ue s Outer Iterations (a) TRLAN–JSYM E i ngen V a l ue s Outer Iterations (b) TRLAN
Fig. 2: Convergence behavior of the large eigenvalues from the random matrix − − − − − − − − || A x j − λ j x j || Outer Iterations (a) TRLAN–JSYM − − − − − − − − || A x j − λ j x j || Outer Iterations (b) TRLAN
Fig. 3: Residual history for the large eigenvalues from the random matrix J -symmetry cause the reordering of the eigenvalues for theTRLAN. As mentioned in Section 2, the TRLAN tends to evaluate one eigenvector ofa pair of the doubly degenerate eigenvalues in the early stage of the iterations. Becauseof the finite precision arithmetic, it loses complete orthogonality to the other half of thedegenerate eigenspace during the iterations, so that the other eigenvalue paired to theconverged eigenvalue emerges in the later stage. Similar behaviors are also observedfor other random matrices at the same algorithmic parameter. For the cases with alarger m , we do not observe the eigenvalue reordering with the TRLAN algorithm,because they quickly converge. Increasing m , N conv rapidly decreases, as shown inTable 1. However, N MV is almost constant or slightly increasing.We compare the computational cost of the algorithms according to the discussiondone in the previous section. The natural parameter choice for the TRLAN algorithm isto employ numbers twice as large as those of the TRLAN–JSYM algorithm. The N conv are comparable among algorithms paired with the doubled parameters, revealing that thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems 13 N conv (cid:39) N (cid:48) conv holds for (41) and (42). The number of matrix-vector multiplicationsof the TRLAN–JSYM algorithm is smaller by roughly a factor of two than that ofthe TRLAN algorithm, as seen in the table. Even with the same maximum Krylovdimension m for both algorithms, e.g. the TRLAN–JSYM with ( nev , mwin , m ) =( , , ) and the TRLAN with ( , , ) , the TRLAN–JSYM algorithm stillbeats the TRLAN algorithm because m does not drastically change the number ofmatrix-vector multiplications.4.2 Case B We evaluate the proposed algorithm 1 for a matrix that appears in a quantum fieldtheory called the twisted Eguchi-Kawai (TEK) model with adjoint fermions [9,10,11].The equation of motion for the adjoint fermions follows from a matrix D called theWilson–Dirac operator in the physics literature. The matrix D = ( D i , j ) is defined by D i , j = δ α , β δ a , b − κ ∑ µ = (cid:2)(cid:0) δ α , β − ( γ µ ) α , β (cid:1) ( V µ ) a , b + (cid:0) δ α , β + ( γ µ ) α , β (cid:1) ( V µ ) b , a (cid:3) , (43)where i and j are collective indices of i = ( a , α ) and j = ( b , β ) , respectively. V µ are ( N − ) × ( N − ) matrices satisfying V T µ = V H µ , i.e. real orthonormal matrices in theadjoint representation of the SU( N ) group, and a , b denote the group indices runningin 1 , . . . , N − γ µ denotes 4 × { γ µ , γ ν } = δ µ , ν I , and α , β denote spin indices running from 1 to 4. Anexplicit form of γ µ can be seen in [17]. The parameter κ implicitly determines themass of the fermion. For more details of D , we refer to [11,14].The matrix D satisfies the following properties: γ D γ = D H , (44) CDC T = D T , (45)where γ = γ γ γ γ and C = γ γ . The matrices γ and C act only on the spin indicesin this notation. We employ the definition for γ µ from [17] and give the explicit formin Appendix A. In this case, γ is real and symmetric and C is real and skew symmetric C T = − C . Monte Carlo methods have been used for simulating quantum field theories.For the system considered herein, the quantum field V µ corresponds to the stochasticvariable in a Monte Carlo algorithm. The spectrum of D becomes stochastic becauseit depends on V µ .We test the proposed algorithm for the matrix A defined by A ≡ ( D γ ) = DD H . (46)The matrix A is Hermitian and J -symmetric with J = C γ . The distribution of thesmall eigenvalues of A is physically important because it carries the information aboutthe dynamics of the theory. The details of the algebraic property of D and A are givenin Appendix A. Table 2: Algorithm parameters and statistics for convergence (Case B)
Large eigenvalues Small eigenvaluesnev mwin m N conv N MV Time[sec] N conv N MV Time[sec]TRLAN–JSYM 4 8 24 54 817 103 . .
248 17 673 133 . .
48 16 24 711 1576 334 . .
148 25 770 173 . .
316 32 48 145 1128 352 . .
796 16 986 403 . . . .
696 15 1170 228 . .
316 32 48 628 2087 512 . .
596 24 1408 316 . .
332 64 96 337 2332 866 . . . . We set N =
289 of SU( N ) for the test. The dimension of A is 4 × ( − ) = V µ is generated with a Monte Carlo algorithm at a parameter set ofthe TEK model. We employ a single Monte Carlo sample of V µ for the test.We compare the convergence behavior of the eigenvalues between the proposedalgorithm (TRLAN–JSYM) and the standard (single vector) thick-restart Lanczosalgorithm (TRLAN). We use the normal and invert modes for solving large and smalleigenvalues, respectively. The conjugate–gradient (CG) algorithm is used in the invertmode. The algorithmic parameters, the number of desired eigenvalues nev, the restartwindow size mwin, and the maximum size of the search dimension m are shownin Table 2. We also tabulate the results of the outer iteration count, the number ofmatrix-vector multiplications, and the computational time for the convergence. Thetimings are shown as reference values, showing how the cost of the matrix-vectormultiplication dominates the computational time in actual applications. We note thatthe convergence behavior of the CG in the invert mode was almost identical betweenthe two algorithms as has been assumed in Section 3, justifying the cost comparisonin terms of the number of matrix-vector multiplications of A − v j in the invert mode.Compared with the TRLAN–JSYM, we double the parameters of the TRLAN tofind all doubly degenerate eigenvalues. We set the tolerance to be 10 − for theeigensolvers.Figs. 4 and 5 show the convergence behavior and residual history of the largeeigenvalues, respectively. The algorithmic parameters are ( nev , mwin , m ) = ( , , ) for the TRLAN–JSYM and ( , , ) for the TRLAN, respectively. We observesimilar convergence behavior as in Case A, where several reorderings occur amongapproximate eigenvalues during the iterations in the TRLAN algorithm. The sameconvergence behavior is seen in Figs. 6 and 7 for the small eigenvalues. thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems 15 E i ngen V a l ue s Outer Iterations (a) TRLAN–JSYM E i ngen V a l ue s Outer Iterations (b) TRLAN
Fig. 4: Convergence behavior of the large eigenvalues (Case B) − − − − − − − − || A x j − λ j x j || Outer Iterations (a) TRLAN–JSYM − − − − − − − − || A x j − λ j x j || Outer Iterations (b) TRLAN
Fig. 5: Residual history for the large eigenvalues (Case B)The computational costs are compared in Table 2. According to the discussiondone in Section 3, most cases satisfy N conv (cid:39) N (cid:48) conv for (41) and (42) among al-gorithms paired with doubled parameters, and N MV for the TRLAN–JSYM algo-rithm is approximately twice as small as that for the TRLAN algorithm. Three cases, ( nev , mwin , m ) = ( , , ) for both modes, and ( nev , mwin , m ) = ( , , ) for theinvert mode, are the exceptions. In these cases, one or two eigenvalues, which are thelargest for the invert mode or the smallest for the normal mode among nev eigenvalues,show slow convergence. Even for these cases, however, the TRLAN–JSYM showssmoother convergence behavior than that of the TRLAN. All the cases we have inves-tigated show that the TRLAN–JSYM algorithm has better computational cost thanthat of TRLAN regarding N MV . The computational timings are roughly proportionalto N MV , indicating that the matrix-vector multiplication dominates the timings. Withthe invert mode for small eigenvalue problems, the timings are well proportional to thenumber of matrix-vector multiplications, because the CG algorithm is used for A − v and the cost of the Lanczos and the true residual computing parts become negligible. E i gen V a l ue s Outer Iterations (a) TRLAN–JSYM E i gen V a l ue s Outer Iterations (b) TRLAN
Fig. 6: Convergence behavior of the small eigenvalues (Case B) − − − − − − − − || A x j − λ j x j || Outer Iterations (a) TRLAN–JSYM − − − − − − − − || A x j − λ j x j || Outer Iterations (b) TRLAN
Fig. 7: Residual history for the small eigenvalues (Case B)
In this study, we have shown the orthogonality and A -orthogonality between the twoKrylov subspaces, K k ( A , v ) and K k ( A , Jv ∗ ) that are generated with the Lanczos algo-rithm for Hermitian J -symmetric matrices A . By employing this property, we proposedthe thick-restarted Lanczos algorithm for Hermitian J -symmetric matrices (TRLAN–JSYM) using which we could efficiently search for one half of the doubly degenerateeigenvectors in K k ( A , v ) without the need to explicitly construct K k ( A , Jv ∗ ) . Theother half of the degenerate eigenvectors are simply constructed from the convergedeigenvectors by utilizing the J -symmetry property.We demonstrated the proposed algorithm TRLAN–JSYM for two test cases, therandom matrices, and the fermion matrix from the quantum field theory called the TEKmodel. The convergence observed for the TRLAN–JSYM algorithm was smootherthan that for the TRLAN algorithm, as expected. The TRLAN–JSYM algorithm per-formed better than the TRLAN algorithm regarding the matrix-vector multiplication.The TRLAN algorithm shows the reordering of eigenvalues among the degenerated thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems 17 eigenvalues caused by the loss of orthogonality between the eigenvectors paired with J -symmetry during the standard Lanczos iteration with the finite precision arithmetic.We did not discuss the mathematical background on the loss of orthogonality in theLanczos algorithm. If exact arithmetic was employed for both of the algorithms, theTRLAN algorithm becomes identical to the TRLAN–JSYM algorithm, according toCorollary 1. However, our algorithm enforces the orthogonality to achieve the smoothconvergence behavior at finite precision arithmetic, resulting in better performance. Acknowledgements
Numerical computations are performed on the ITO supercomputer system of Kyushuuniversity. This work is partly supported by Priority Issue 9 to be tackled by using Post K Computer. Wethank Antonio Gonz´alez-Arroyo and Masanori Okawa for comments on the draft and the details of theWilson-Dirac operator in the adjoint representation. We are very grateful to the anonymous reviewer forhis/her comments that enhanced the quality of the manuscript.
A Properties of the matrices D (43) and A (46) In this appendix we show the algebraic properties of the matrices D (43) and A (46), including the J -symmetry. We employ the following explicit form for γ µ : γ = − i − i i i , γ = −
10 0 1 00 1 0 0 − , γ = − i
00 0 0 ii − i , γ = − − . (47)In addition to these, we also have γ = γ γ γ γ as γ = . (48)These γ µ matrices have the following properties: { γ µ , γ ν } = δ µ , ν I , for µ , ν = ,..., , (49) γ H µ = γ µ , for µ = ,..., , (50) γ T = − γ , γ T = γ , γ T = − γ , γ T = γ , γ T = γ . (51)The matrix C = γ γ has the following properties: C − = γ γ = − γ γ = − C , (52) C − = γ γ = γ T γ T = ( γ γ ) T = C T , (53) C γ µ C T = γ γ γ γ γ = − γ γ γ = γ = − γ T ( µ = ) γ γ γ γ γ = γ γ γ = − γ = − γ T ( µ = ) γ γ γ γ γ = − γ γ γ = γ = − γ T ( µ = ) γ γ γ γ γ = − γ γ γ = − γ = − γ T ( µ = ) = − γ T µ . (54)We show the properties of (44) and (45) in detail. To simplify the proof, we suppress the matrix indicesof (43) and write it as D = I − κ ∑ µ = (cid:104) ( − γ µ ) V µ + ( + γ µ ) V T µ (cid:105) , (55)8 Ken-Ichi Ishikawa, Tomohiro Sogabewhere the direct product of the spinor index and the color index is implicit.Equation (44) is shown as: γ D γ = γ (cid:32) I − κ ∑ µ = (cid:104) ( − γ µ ) V µ + ( + γ µ ) V T µ (cid:105)(cid:33) γ = I − κ ∑ µ = (cid:104) γ ( − γ µ ) γ V µ + γ ( + γ µ ) γ V T µ (cid:105) = I − κ ∑ µ = (cid:104) ( + γ µ ) V µ + ( − γ µ ) V T µ (cid:105) , where { γ , γ µ } = V T µ = V H µ for the matrices in the adjoint representation of the SU( N )group, the last line is identical to = (cid:32) I − κ ∑ µ = (cid:104) ( − γ µ ) V µ + ( + γ µ ) V T µ (cid:105)(cid:33) H = D H . (56)Next, we show (45) in: CDC T = C (cid:32) I − κ ∑ µ = (cid:104) ( − γ µ ) V µ + ( + γ µ ) V T µ (cid:105)(cid:33) C T = I − κ ∑ µ = (cid:104) C ( − γ µ ) C T V µ + C ( + γ µ ) C T V T µ (cid:105) = I − κ ∑ µ = (cid:104) ( + γ T µ ) V µ + ( − γ T µ ) V T µ (cid:105) , where we used (54). The last line is identical to = (cid:32) I − κ ∑ µ = (cid:104) ( − γ µ ) V µ + ( + γ µ ) V T µ (cid:105)(cid:33) T = D T . (57)We finally show the J -symmetry of A (46). The Hermiticity of A is apparent from (44) and (46). Theproperties of J ≡ C γ = γ γ γ are: J − = γ γ γ = − γ γ γ = γ γ γ = − γ γ γ = − J , (58) J − = γ γ γ = γ T γ T γ T = ( γ γ γ ) T = J T . (59)Because γ , γ , and γ are real matrices, J is real. Therefore, the properties of J ≡ C γ follows those in (1).The J -symmetry of A (46) is shown as: JAJ − = C γ D γ D γ γ C T = C γ D γ DC T = γ CD γ DC T = γ CDC T C γ DC T = γ CDC T γ CDC T = γ D T γ D T = ( D γ ) T ( D γ ) T = [( D γ )( D γ )] T = A T , (60)where we used C γ = γ C , C T C = I , γ T = γ , ( γ ) = I , and (45). Therefore, A of (46) is J -symmetric. thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems 19 References
1. Baglama, J., Calvetti, D., Reichel, L.: IRBL: An Implicitly Restarted Block-Lanczos Methodfor Large-Scale Hermitian Eigenproblems. SIAM J. Sci. Comput. (5), 1650–1677 (2003).doi:10.1137/S1064827501397949. URL https://doi.org/10.1137/S1064827501397949
2. Bai, Z., Demmel, J., Dongarra, J., Ruhe, A., van der Vorst, H.: Templates for the Solu-tion of Algebraic Eigenvalue Problems. Society for Industrial and Applied Mathematics(2000). doi:10.1137/1.9780898719581. URL https://epubs.siam.org/doi/abs/10.1137/1.9780898719581
3. Benner, P., Faßbender, H.: An implicitly restarted symplectic Lanczos method for the Hamiltonianeigenvalue problem. Linear Algebra Appl. , 75 – 111 (1997). doi:https://doi.org/10.1016/S0024-3795(96)00524-1. URL
4. Benner, P., Faßbender, H., Stoll, M.: A Hamiltonian Krylov–Schur-type method basedon the symplectic Lanczos process. Linear Algebra Appl. (3), 578 – 600 (2011).doi:https://doi.org/10.1016/j.laa.2010.04.048. URL
5. Benner, P., Faßbender, H., Yang, C.: Some remarks on the complex J -symmetric eigenproblem.Linear Algebra Appl. , 407 – 442 (2018). doi:https://doi.org/10.1016/j.laa.2018.01.014. URL
6. Calvetti, D., Reichel, L., Sorensen, D.C.: An implicitly restarted Lanczos method for large symmetriceigenvalue problems. Electron. Trans. Numer. Anal. , 1–21 (1994)7. Dongarra, J., Gabriel, J., Koelling, D., Wilkinson, J.: The eigenvalue problem for Hermi-tian matrices with time reversal symmetry. Linear Algebra Appl. , 27 – 42 (1984).doi:https://doi.org/10.1016/0024-3795(84)90068-5. URL
8. Golub, G., Underwood, R.: The Block Lanczos Method for Computing Eigenvalues. In: J.R. Rice(ed.) Mathematical Software, pp. 361 – 377. Academic Press (1977). doi:https://doi.org/10.1016/B978-0-12-587260-7.50018-2. URL
9. Gonz´alez-Arroyo, A., Okawa, M.: A twisted model for large N lattice gauge theory. Phys. Lett. ,174–178 (1983). doi:10.1016/0370-2693(83)90647-010. Gonz´alez-Arroyo, A., Okawa, M.: Twisted-Eguchi-Kawai model: A reduced model for large- N latticegauge theory. Phys. Rev. D27 , 2397 (1983). doi:10.1103/PhysRevD.27.239711. Gonz´alez-Arroyo, A., Okawa, M.: Twisted space-time reduced model of large N QCD with two adjointWilson fermions. Phys. Rev.
D88 , 014514 (2013). doi:10.1103/PhysRevD.88.01451412. Kressner, D.: Numerical Methods for General and Structured Eigenvalue Problems,
Lecture Notes inComputational Science and Engineering , vol. 46, 1st edn. Springer-Verlag Berlin Heidelberg (2005).doi:10.1007/3-540-28502-4. URL
13. Mehrmann, Volker. and Watkins, David.: Structure-Preserving Methods for Computing Eigenpairs ofLarge Sparse Skew-Hamiltonian/Hamiltonian Pencils. SIAM J. Sci. Comput. (6), 1905–1925 (2001).doi:10.1137/S1064827500366434. URL https://doi.org/10.1137/S1064827500366434
14. Montvay, I.: SUPERSYMMETRIC YANG–MILLS THEORY ON THE LATTICE. Int. J. Mod. Phys.
A17 , 2377–2412 (2002). doi:10.1142/S0217751X0201090X15. Petkov, M.G., Ivanov, I.G.: Solution of Symmetric and Hermitian J -symmetric Eigenvalue Problem.Mathematica Balkanica. New series , 337 – 349 (1994). URL
16. Research Institute for Information Technology, Kyushu University, “Introduction of ITO,” 2018, URL
17. Rothe, H.J.: Lattice Gauge Theories, 4th edn. WORLD SCIENTIFIC (2012). doi:10.1142/8229. URL
18. Shimizu, N., Mizusaki, T., Utsuno, Y., Tsunoda, Y.: Thick-restart block Lanczos methodfor large-scale shell-model calculations. Comput. Phys. Commun. , 372 – 384 (2019).doi:https://doi.org/10.1016/j.cpc.2019.06.011. URL
19. Stewart, G.W.: A Krylov–Schur Algorithm for Large Eigenproblems. SIAM J. Matrix Anal. Appl. (3), 601–614 (2001). doi:10.1137/S0895479800371529. URL http://dx.doi.org/10.1137/S0895479800371529 (1), 156 – 173 (1999). doi:https://doi.org/10.1006/jcph.1999.6306.URL
21. Wu, K., Simon, H.: Thick-Restart Lanczos Method for Large Symmetric Eigenvalue Problems. SIAMJ. Matrix Anal. Appl. (2), 602–616 (2000). doi:10.1137/S0895479898334605. URL https://doi.org/10.1137/S0895479898334605
22. Zhou, Y., Saad, Y.: Block Krylov–Schur method for large symmetric eigenvalue problems. Numer.Algorithms (4), 341–359 (2008). doi:10.1007/s11075-008-9192-9. URL https://doi.org/10.1007/s11075-008-9192-9 thick-restart Lanczos type method for Hermitian J -symmetric eigenvalue problems 21 Algorithm 1
The thick-restart Lanczos algorithm for a Hermitian J -symmetric matrix A (TRLAN–JSYM). Require:
Maximum Krylov subspace dimension size m , restart window size mwin, and number of desiredeigenpairs nev. Ensure:
Eigenpairs ( x i , λ i ) and residual norms || r i || = || Ax i − λ i x i || for i = ,..., nev in V ( : , ) , ev ( ) , res ( ) .1: k = v = v = v / || v || (cid:46) Initial unit vector3: w = Jv ∗ (cid:46) Initial dual vector4: loop
5: LANCZOS JSYM( m , k , V m + , W m + , ¯ T m + ) (cid:46) m − k -step Lanczos with J -symmetry6: T m = Z m Λ m Z Tm (cid:46) Compute eigenpairs of T m
7: Move desired eigenpairs in the top dimensions of Z m and Λ m by sorting.8: V m : = V m Z m (cid:46) Compute approximate eigenvectors9: T m = (cid:46) Compute the Krylov–Schur transformation for ¯ T m + for i = ,..., m do t i , i = λ i t m + , i = t m + , m z m , i end for if Normal Mode then (cid:46)
Compute estimated residuals15: for i = ,..., nev do
16: ev ( i ) = λ i
17: res est ( i ) = | t m + , i | end for else if Invert Mode then (cid:46)
Compute estimated residuals20: c = || Av m + || for i = ,..., nev do
22: ev ( i ) = / λ i
23: res est ( i ) = c | t m + , i ev ( i ) | end for end if for i = ,..., nev do if res est ( i ) < tol then (cid:46) Check true residuals28: res ( i ) = || Av i − v i ev ( i ) || if res ( i ) < tol then
30: is convd ( i ) = . TRUE . else
32: is convd ( i ) = . FALSE . end if end if end for
36: Move converged eigenpairs in the top dimensions of Z , T , V , res , res est , ev with the key is convdby sorting.37: icnv = for i = ,..., nev do if is convd ( i ) == . TRUE . then t m + , i = (cid:46) Decouple converged subspace41: icnv = icnv + (cid:46) Count number of converged eigenpairs42: end if end for if icnv == nev then Exit Loop end if (cid:46)
Shrink Krylov–Schur decomposition to k + k = MIN ( icnv + mwin , m − ) for i = ,..., k do t k + , i = t m + , i end for v k + = v m + W k + = JV ∗ k + end loop Algorithm 2 m − k -step Lanczos iteration for a Hermitian J -symmetric matrix A . procedure LANCZOS JSYM( n , m , k , V m + , W m + , ¯ T m ) Require: V k + , W k + , ¯ T k Ensure: V m + , W m + , ¯ T m γ = √ (cid:46) Reorthogonalization threshold parameter3: for j = k + ,..., m do if Normal Mode then v j + = Av j else if Invert Mode then v j + = A − v j end if t j , j = v Hj v j + v j + : = v j + − v j t j , j b = || v j + || loop for i = ,..., j do c = w Hi v j + v j + : = v j + − w i c (cid:46) Orthogonalization to W j c = v Hi v j + v j + : = v j + − v i c end for b = || v j + || if b γ > b then Exit Loop end if b = b end loop v j + : = v j + / b t j + , j = b w j + = Jv ∗ j + (cid:46) Construct W j + end for29: