A thick-restart Lanczos type method for Hermitian $J$-symmetric eigenvalue problems

A thick-restart Lanczos type algorithm is proposed for Hermitian $J$-symmetric matrices. Since Hermitian $J$-symmetric matrices possess doubly degenerate spectra or doubly multiple eigenvalues with a simple relation between the degenerate eigenvectors, we can improve the convergence of the Lanczos algorithm by restricting the search space of the Krylov subspace to that spanned by one of each pair of the degenerate 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 orthogonal to each other. The proposed algorithm searches for eigenvectors in one of the two subspaces without the multiplicity. The other eigenvectors paired to them can be easily reconstructed with the simple relation from the $J$-symmetry. We test our algorithm on a matrix originating from a quantum field theory.


Introduction
In theoretical physics, symmetry is an important guiding principle to search for laws of nature. When numerically simulating a physics system that has a symmetry, it is preferable to retain the symmetry property in the numerical simulation algorithm. Hermitian matrices often appear in analyzing physics systems, and spectra analysis is important to understand the nature. In this case the Hermitian symmetry must be considered in the spectra analysis and the eigensolver algorithm. Developing eigensolver algorithms for spectra analysis is one of the major research areas in applied 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 is defined by where J is a skew-symmetric matrix in R n×n . Throughout this paper, we employ the following notations: T for matrix transposition, H for Hermitian conjugation, and * for complex conjugation. We use I and O to denote the identity and null matrices with appropriate sizes, respectively. The eigenvalues 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. The vector defined by is also the eigenvector of A associated to λ because of the symmetry (1).
To investigate the spectrum of a sparse matrix in a large dimension, iterative eigensolver algorithms are generally employed. Iterative eigensolver algorithms existing in the literature, such as the Lanczos algorithm for Hermitian matrices, do not consider the J-symmetry. Therefore, even after convergence on one of a pair of the degenerate eigenvector pairs, the algorithm continues to search for the other eigenvector associated to the converged eigenvalue even though the other eigenvector can be easily reconstructed from the converged one. In this article, we restrict ourselves to improving the Lanczos algorithm for large Hermitian J-symmetric matrices by incorporating the J-symmetry property.
We could improve the convergence of the Lanczos algorithm by restricting the search space of the Krylov subspace of A to that spanned by one of each pair of the doubly degenerate eigenvector pairs by imposing complete orthogonality to the subspace spanned by the other eigenvectors paired to them. However, this strategy cannot be directly realized before knowing the invariant subspace of A. We find that the Krylov subspace K k (A, v) = span{v, Av, . . . , A k−1 v} generated with the Lanczos algorithm with a starting vector v is orthogonal to K k (A, Jv * ) = span{Jv * , AJv * , . . . , A k−1 Jv * } with the same Lanczos algorithm with the starting vector Jv * . On the basis of 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 stated above.
Early attempts have been made to incorporate the symmetry property or structure of matrices into eigensolver algorithms in [7,15] in which dense matrix algorithms were developed for the Hermitian J-symmetric matrices that appeared in a quantum mechanical system with time-reversal and inversion symmetry. The one similar to our algorithm that incorporates the symmetry properties of matrices was studied in [3,4,5,13], in which the Lanczos type eigensolvers for complex J-symmetric matrices or Hamiltonian matrices were investigated. Our algorithm for Hermitian J-symmetric matrices is more restrictive or specialized than those of [3,4,5,13], because our algorithm adheres to two symmetries, i.e., Hermiticity and J-symmetry. See also [12] for numerical algorithms for structured eigenvalue problems. This paper is organized as follows. In the next section, we prove the orthogonality 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, we describe the restarting method based on the Krylov-Schur transformation method [18] and the thick-restart method [19,20]. Then, we propose the thick-restart Lanczos algorithm for Hermitian J-symmetric matrices in Section 3. In Section 4, we test the proposed algorithm for a Hermitian J-symmetric matrix that originates from a quantum field theory. Thus, the convergence behavior is compared with that of the standard thick-restart Lanczos algorithm. We summarize this paper in the last section.

Lanczos iteration and J-symmetry
The Lanczos iteration transforms A to a tridiagonal form by generating the orthonormal basis vectors. The Lanczos decomposition after m-step iteration starting with a unit vector v 1 is given by where and v j ∈ C n , and e m is the m-dimensional unit vector in the m-th direction. The basis vectors are orthonormal, V H m+1 V m+1 = I. T m denotes the m × m tridiagonal matrix that is given by The approximate eigenpairs are obtained from the eigenpairs of T m . Because of the Hermiticity 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 the standard Lanczos iteration, the Hermitian symmetry of A is respected in the form T m . We also investigate the J-symmetry property of decomposition (3). After taking the complex conjugate of (3) followed by multiplying it by J from the left-hand side, we have Using the J-symmetry and Hermiticity, JA * = A H J = AJ, and by defining w j ≡ Jv * j , we obtain The basis vectors W m+1 are also orthonormal, W H m+1 W m+1 = I. Consequently, the vectors W m+1 have the same Lanczos decomposition as that of V m+1 when it starts from w 1 = Jv * 1 . Furthermore, we can prove the orthogonality between W m+1 and V m+1 . To implement the thick-restart method [19,20], we show the J-symmetry property for a generalized decomposition similar to the Krylov-Schur decomposition [18] instead of the Lanczos decomposition in the following part of the paper. After preparing two lemmas related to the J-symmetry, we prove the main theorem for the orthogonality property between W m+1 and V m+1 on the generalized decomposition. Subsequently, the orthogonality properties for the Lanczos decomposition and the thickrestart method are proved as corollaries.
We first show orthogonal properties between v and w = Jv * .

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 * satisfies Proof From the definition of w, it follows that where (10) yield where J T = −J and JA = A T J are used. The identity v T A T Jv = v T J T Av and (12) yield Thus, w H Av = 0.
We have the following lemma that is a generalization of the relation between (3) and (6).
. . , v k , v k+1 be the vectors having the following relation: where V k = [v 1 , . . . , v k ], S k is a matrix in R k×k , and b is a vector in R k for a Hermitian J-symmetric matrix A satisfying (1). Then, the vectors w j = Jv * j ( j = 1, . . . , k + 1) satisfy the following decomposition: Proof By taking the complex conjugate of (14) followed by multiplying it by J from the left-hand side, we obtain (15) using JA T = AJ and A * = A T .
Using Lemmas 1 and 2, we can show the orthogonality properties between W k+1 and V k+1 as follows.
] be the vectors satisfying (14), and W k+1 = JV * k+1 . If the vectors V k and W k are orthogonal and A-orthogonal to each other: V H k W k = O and V H k AW k = O, then the vectors V k+1 and W k+1 also satisfy the following orthogonality relations: Proof The decomposition (15) (14) and (15), respectively, yields From the premise that Together with Lemma 1 and the premise, we find Multiplying w H k+1 and v H k+1 to (14) and (15), respectively, yields Because of (20) and (21) as well as Lemma 1, the right-hand sides of (23) and (24) vanish. Thus, Together with Lemma 1 and the premise, we find Therefore, V k+1 and W k+1 are orthogonal and A-orthogonal to each other.
Using Theorem 1 as well as Lemmas 1 and 2, we can show that the Lanczos vectors V m+1 generated with (3) are orthogonal and A-orthogonal to W m+1 = JV * m+1 .

Corollary 1
The m-step Lanczos vectors V m+1 with m ≥ 1 generated with a unit vector v 1 for a Hermitian J-symmetric matrix A satisfy with when no breakdown occurs.
Proof Because the Lanczos decomposition (3) is a particular form of (14) with k → m, a real matrix S m → T m and a real vector b → β m e m , Lemma 2 can be applied to obtain (6). Because w H 1 v 1 = 0 and w H 1 Av 1 = 0 hold from Lemma 1, we can apply Theorem 1 to (3) and (6) when m = 1. Moreover, the Lanczos decomposition retains its form applicable to Theorem 1 for any m > 1. Therefore, the corollary follows from Theorem 1 and Lemmas 1 and 2 by induction.
We cannot simultaneously find degenerate pairs of eigenvectors with the standard single-vector Lanczos process. This is true for computation with exact arithmetic. However, with finite precision arithmetic, the single-vector Lanczos process would generate a small overlap to W k via round-off errors, so that even after the convergence of 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 be applied to accelerate the convergence for degenerate eigenvectors [1,2,8,17,21].
By using Corollary 1, we can construct a Lanczos type algorithm in which the orthogonality to W k is enforced to search for eigenvectors without multiplicity of eigenvalues associated to the J-symmetry. Additionally, the other vectors paired to them can be easily reconstructed. However, for a practical numerical algorithm of the Lanczos type iteration, the iteration should terminate at a finite step, and a restarting mechanism is required [6,19,20]. The most useful and simplest but effective restarting method is the so-called thick-restart method [19,20] that is a specialization of the Krylov-Schur transformation [18] to Hermitian matrices. 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.
Corollary 2 Let V m+1 and W m+1 be the orthonormal basis vectors generated with the m-step Lanczos process, (3) and (6), and Z k be an orthonormal matrix in R k×k . The Krylov-Schur transformation with Z k on the Lanczos decomposition is defined by  (31) are particular forms of the decomposition in Lemma 2, the orthogonality and A-orthogonality between U m+1 and Q m+1 simply follow from Theorem 1.
For the thick-restart method, Z m is chosen to diagonalize T m , and the dimension of the decomposition is reduced from m to k < m with a selection criterion for vectors V k ← V m . In the reduction, the last vectors are kept hold as v k+1 ← v m+1 and w k+1 ← w m+1 to retain the decomposition form properly. The orthogonality properties of the new basis (U m+1 , Q m+1 ) and the reduced basis (U k+1 , Q k+1 ) still hold according to Corollary 2. After restarting, the Lanczos iteration continues to keep the decomposed form 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 is orthogonal and A-orthogonal to K k (A, Jv * ).

Thick-restart Lanczos algorithm with J-symmetry
Based on Theorem 1 as well as Corollaries 1 and 2, we construct a thick-restart Lanczos algorithm for Hermitian J-symmetric matrices (TRLAN-JSYM) which efficiently 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 small eigenvalues.
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 enforce the orthogonality of V j+1 to W j = JV * j to avoid the contamination of the search space So far, we have described the single-vector Lanczos iteration type algorithm to introduce the TRLAN-JSYM. If the matrix A has a dense cluster of eigenvalues or multiple eigenvalues other than those with the J-symmetry, we need to incorporate the block type Lanczos iteration in the algorithm for efficiency. We can extend the proposed algorithm to the blocked version in the same manner as it was conducted for the standard thick-restart Lanczos algorithm [17,21]. The study on the block version will be addressed in future studies and we have only shown the single vector version to demonstrate the idea for simplicity.

Definition of the Test Matrix
We test the proposed algorithm 1 for a matrix that appears in a quantum field theory 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 the Wilson-Dirac operator in the physics literature. The matrix D = (D i, j ) is defined by where i and j are collective indices of i = (a, α) and j = (b, β ), respectively. V µ are (N 2 − 1) × (N 2 − 1) matrices in the adjoint representation of the SU(N) group, and a, b denote the group indices running in 1, . . . , N 2 − 1. γ µ denotes 4 × 4 Hermitian matrices satisfying the anti-commuting relation {γ µ , γ ν } = 2δ µ,ν I, and α, β denote spin indices running from 1 to 4. An explicit form of γ µ can be seen in [16]. The parameter κ implicitly determines the mass of the fermion. For more details of D, we refer to [11,14].
The matrix D satisfies the following properties: where γ 5 = γ 4 γ 1 γ 2 γ 3 and C = γ 4 γ 2 . The matrices γ 5 and C act only on the spin indices in this notation. We consider the explicit form of γ µ from [16]. In this case, γ 5 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 stochastic variable in a Monte Carlo algorithm. The spectrum of D becomes stochastic because it depends on V µ . We test the proposed algorithm for the matrix A defined by The matrix A is Hermitian and J-symmetric with J = Cγ 5 . The distribution of the small eigenvalues of A is physically important because it carries the information about the dynamics of the theory.

Numerical Results
We set N = 289 of SU(N) for the test. The dimension of A is 4×(289 2 −1) = 334080. The ensemble for V µ is generated with a Monte Carlo algorithm at a parameter set of the TEK model. We employ a single Monte Carlo sample of V µ for the test. Now, we compare the convergence behavior of the eigenvalues between the proposed algorithm (TRLAN-JSYM) and the standard (single vector) thick-restart Lanczos algorithm (TRLAN). We use the normal and invert modes for solving large and  iterations, while several reorderings occur for the TRLAN. The reordering of the eigenvalues for the TRLAN is caused by the sorting algorithm for the eigenvalues and the property of the Lanczos iteration to the J-symmetry. As we mentioned in Section 2, the TRLAN tends to evaluate one eigenvecotr of a pair of the doubly degenerate eigenvalues in the early stage of the iterations. It loses complete orthogonality to the other half of the degenerate eigenspace during the iterations, so that the other eigenvalue paired to the converged eigenvalue emerges in the later stage. Even though both the algorithms yield the desired number of eigenvalues within a similar iteration count, the improvement on the computational cost of the TRLAN-JSYM is far more evident because the subspace dimension of the TRLAN-JSYM is one half of that of the TRLAN. We observe the same convergence behavior as for the small eigenvalues in Figs. 3 and 4.

Summary
In this study, we have shown the orthogonality and A-orthogonality between the two Krylov subspaces, K k (A, v) and K k (A, Jv * ) that are generated with the Lanczos algorithm for Hermitian J-symmetric matrices A. By employing this property, we proposed the thick-restarted Lanczos algorithm for Hermitian J-symmetric matrices (TRLAN-JSYM) using which we could efficiently search for one half of the doubly degenerate eigenvectors in K k (A, v) without the need to explicitly construct K k (A, Jv * ). The other half of the degenerate eigenvectors are simply constructed from the converged eigenvectors by utilizing the J-symmetry property. We demonstrated the proposed algorithm TRLAN-JSYM for the fermion matrix from the quantum field theory called the TEK model. The convergence observed for the TRLAN-JSYM algorithm was smother than that for the TRLAN algorithm as expected.