Recent improvements in lattice QCD calculations have been achieved by making use of the eigenvalue spectrum of lattice operators, such as the lattice Dirac operator and the lattice Laplacian operator. The Lanczos iteration method is an eigensolver algorithm employed for that purpose [1]. It enables calculation of both eigenvalues and eigenvectors of a Hermitian matrix. It is most useful when a partial eigenvalue spectrum is required. Additional iterations provide more eigenmodes with higher accuracy.
There are several improvement methods for the Lanczos algorithm: (1) Polynomial acceleration applies a polynomial to the matrix of interest so that the Lanczos iteration converges faster with a transformed eigenvalue spectrum. Chebyshev polynomials are a popular choice [2]. (2) There are two kinds of restarted Lanczos algorithms that take advantage of information learned from the previous run. The Implicitly Restarted Lanczos (IRL) algorithm suppresses unwanted eigenvalues at the restart [3]. Alternatively, the Thick-Restart Lanczos (TRL) algorithm suppresses or explicitly deflates converged wanted eigenvalues [4]. (3) There are approaches that utilize blocking to improve the computational e ciency of the Lanczos algorithm. The Block Lanczos (BL) algorithm can enhance the compute-to-communication ratio with the Split Grid method [5].
Grid [6,7] and QUDA [8][9][10] are very popular and highly performing lattice QCD libraries, especially for modern parallel computing systems. For the Lanczos eigensolver, Grid has IRL and BL eigensolvers. QUDA has a TRL eigensolver and its block variant. All of them utilize Chebyshev polynomials to enhance the convergence.
We investigate these Lanczos eigensolver algorithms in the Grid and the QUDA libraries measuring and analyzing their performance for the Highly Improved Staggered Quark (HISQ) Dirac operator. In Section 2, we describe our simulation details. In Section 3, we describe some basics of the Lanczos iteration method. Section 4 describes the general strategy of utilizing Chebyshev polynomials and studies its optimization. In Sections 5 and 6, we illustrate ideas of IRL and TRL, respectively, and discuss their optimizations. In Section 7, we examine the performance of BL with the Split-Grid method. We conclude in Section 8.
We calculate eigenvalues and eigenvectors of the low modes of the massless HISQ Dirac operator ⇡. The Lanczos eigensolvers are run using a Hermitian variation ⇡ † ⇡, which is (semi)-positive definite that allows even-odd splitting. The eigensolver performance is measured in elapsed time taken by the eigensolver routine itself. We use an # 5 = 2 + 1 + 1 MILC HISQ gauge configuration, with lattice size 24 3 ⇥ 64, lattice spacing 0 ⇡ 0.12 fm, and 0< ; /0< B = 0.00507/0.0507 [11].
We run Grid on Intel CPUs and QUDA on NVIDIA GPUs, where they perform their best so that we can focus on algorithmic performance. All Grid runs are done on a single node of dual Intel Xeon E5-2650v2 system that has 2 ⇥ 8 = 16 CPU cores, where we use only half, i.e., eight, of them for performance stability. We apply OpenMP for their parallelization, except for the BL eigensolver in Section 7, where we apply MPI. QUDA runs are done on the same system but with two NVIDIA K80 GPUs.
Let us consider an < ⇥ < Hermitian matrix , and an <-dimensional vector 1. The Gram-Schmidt process on a set of vectors {1, 1, 2 1, • • • , = 1 1} for = <, which spans a Krylov subspace K = ( , 1), gives a three term recurrence relation
where @ 8 are orthonormal basis vectors of K = ( , 1), and C 8 9 = ⌦ @ 8 @ 9 ↵
. Combining recurrences of 8 = 1, • • • , =, we have
where
) is an = ⇥ = Hermitian tridiagonal matrix. The Lanczos algorithm is to construct ) = and & = by iterating the Lanczos recurrence relation Eq. ( 1) [1]. When = = <, & < is a unitary transformation, so ) < has the same eigenvalue spectrum as . Even for = < <, eigenvalues of ) = are approximate to some = eigenvalues of [12]. They converge to true eigenvalues of as = increases. The largest, the smallest, or the most isolated (the least dense) extreme eigenvalues appear first and converge first. When : (< <) eigenvalues are wanted, we need to perform = (: = <) iterations of the Lanczos recurrence (Eq. ( 1)) with the result that ) = has eigenvalues approximate to : eigenvalues of of interest within the target precision.
A QR iteration is usually employed to diagonalize a tridiagonal matrix ) = [13,14]. Then its eigenvalue decomposition gives
where ⇤ is a diagonal matrix with eigenvalues _ 8 of ) = , and + = is composed of their corresponding eigenvectors E 8 . With this, we calculate eigenvector estimates (called Ritz vectors) F 8 and eigenvalue estimates (called Ritz values) _8 by
The convergence of each eigenvector is measured by a residual
or dividing it by some normalization factors. Note that _8 can be di erent from _ 8 since we usually apply a polynomial to , as will be discussed in Section 4.
The convergence of a Lanczos iteration depends on the eigenvalue density. Less dense (including the largest and the smallest) eigenvalues appear early and converge fast. One can control the density of eigenvalues, and, thus, the convergence of Lanczos iteration, by applying a polynomial P to the matrix . It maps eigenvalues _ 8 of into P (_ 8 ) but preserves their eigenvectors.
Chebyshev polynomials ⇠ ? (G) of degree ? are bounded within [ 1,1] for |G| 1, while diverging rapidly for a high ? as |G| increases for |G| > 1 [15]. We can transform it to map unwanted eigenvalues into the bounded region and wanted eigenvalues into the diverging region. This Chebyshev acceleration can boost the convergence of Lanczos iteration. With suitable choices of the location of the bounded region and the polynomial degree, its benefit more than compensates for the cost of applying the polynomial.
Eigenvalues of ⇡ † ⇡ are non-negative, and our usual interests lie in its low modes. We transform the domain of ⇠ ? (G) so that [ 1, 1] becomes [U, V] that covers the unwanted eigenvalues. The value of U should be set to be greater than the largest desired eigenvalue, and V to be greater than the largest eigenvalue in the whole spectrum. A few power iterations can help to estimate the size of the largest eigenvalue. Choosing U needs some heuristics, but we may use the same value within a gauge ensemble.
In Fig. 1, we show the convergence of Lanczos iterations with di erent Chebyshev acceleration parameters. Here, Grid's ImplictlyRestartedLanczos eigensolver is used with a little modification not to restart. We run 1000 Lanczos iterations, and the convergences are checked by Eq. ( 4) in double precision. Figure 1a presents the number of converged eigenvalues for di erent U and V values. The same Chebyshev polynomial degree ? = 100 is used here, which means all the runs take the same time in principle. On this gauge configuration, the largest eigenvalue of ⇡ † ⇡ is around 21, and the 1000th eigenvalue is around 0.046. The result indicates that we get the best convergence when both U and V are the closest to the 1000th eigenvalue and the largest eigenvalue, respectively. The 800th eigenvalue is around 0.034. Hence, U = 0.04 might give a better convergence on this gauge configuration. However, a tight choice may not work for other gauge configurations in the same ensemble.
In Fig. 1b, we present the results for di erent polynomial degrees ?. Here, we used (U, V) = (0.05, 25), the best pair found in Fig. 1a. Although higher ? makes more eigenvalues converge, it demands more computations. In this example, ? = 100 gives the best e ciency. In general, the most e cient ? depends on the number of eigenvalues we want. We need heuristics to find an optimized ?.
Every improved Lanczos algorithm discussed in this paper utilizes Chebyshev acceleration. We use (U, V) = (0.05, 25), and ? = 50 or 100.
A basic Lanczos algorithm must carry on its iteration by increasing the Krylov subspace dimension = until eigenvalues of the desired number converge. As = increases, it takes more memory to store the eigenvectors as well as more computing time for the iteration1 and for the diagonalization of a larger ) = . Moreover, it is di cult to anticipate an optimized (or the smallest) = without monitoring the convergence, which includes diagonalizing ) 8 (8 < =) and computing eigenvectors and their residuals.
Suppose we have performed = Lanczos iterations for Hermitian matrix (= ⇠ ? (⇡ † ⇡)) with a starting vector 1. The Implicitly Restarted Lanczos (IRL) algorithm allows us to restart from the (= :)-th step of the Lanczos algorithm for with a new starting vector 1 = ( `:) • • • ( `1)1, where the {`8, 8 = 1, . . . :} are some eigenvalues [3]. The Krylov subspace K = ( , 1) for the restarted run does not contain eigenvectors corresponding to the eigenvalues `8. It not only excludes unwanted eigenvalues and eigenvalues near them, but it also improves the convergence of the remaining wanted eigenvalues far from the {`8 }.
When we are interested in only the C lowest eigenmodes, a typical implementation of IRL runs the Lanczos iteration to the =-th step, does the implicit restart that excludes the : largest eigenvalues found in the previous run, and repeats them until the C smallest eigenvalues converge. In this way, we can constrain the Krylov subspace dimension to any = > C, as much as the system memory allows. However, since : is usually set to make = : slightly larger than C, = determines how many restarts we need to do. Each restart must perform the QR decomposition : times and do the basis rotations using them. This cost can be comparable to that of many Lanczos iterations. Hence, we need to find an optimized = that maximizes the restarting e ciency.
The Grid library has an IRL eigensolver named ImplicitlyRestartedLanczos. In Fig. 2, we present its performance profiles for various Krylov subspace dimensions = with Chebyshev acceleration of ? = 50, 100. For non-restarted cases, a modified version mentioned in Section 4 is used. Note that in the Lanczos iteration, eigenvalues converge simultaneously, not one by one in order. Hence, it could be less e cient for calculating a small number of eigenvalues. In Fig. 2a, where we calculate 100 eigenvalues, the non-restarted Lanczos converges at = = 400 and = = 700 for ? = 50 and ? = 100, respectively. In a (relatively) slowly converging case like this, IRL can perform better than the basic non-restarted Lanczos. For both ? = 50 and 100, we obtain the best performance by restarting once.
On the other hand, in Fig. 2b, where we calculate 500 eigenvalues, the non-restarted Lanczos takes the least time. Here, = = 800, 1050 are only about twice 500, so the restarting costs for = = 600, 700, 800 are comparable to or more expensive than running Lanczos on a bigger Krylov subspace. However, the di erences between = = 700, 800, 1050 are not huge for ? = 50, so one may choose = = 700 to save memory. In addition, unlike the non-restarted Lanczos where the convergence is not guaranteed for a given =, IRL assures its convergence when the same = is used for other gauge configurations, even though it may need more restarts.
1In practice, Lanczos vectors @ 8 's orthogonality gets inexact as iteration continues due to the accumulation of numerical errors. One can correct it by re-orthogonalizing all @ 8 's per some period of iterations. This is expensive for a large =. One last remark is that we should optimize the Chebyshev acceleration parameters together with =. In Fig. 2a, the non-restarted Lanczos converges faster with ? = 100 than ? = 50, while the IRL's best performance comes with ? = 50.
In this section, we investigate another restarted Lanczos algorithm called Thick-Restart Lanczos (TRL) [4]. Unlike IRL, which suppresses unwanted eigenvectors to improve convergence, TRL suppresses (nearly) converged wanted eigenvectors from the Krylov subspace for the restarted run.
One can rewrite Eq. ( 2) as
where + : is composed of some : (< =) eigenvectors of ) = , and
) is composed of the corresponding : eigenvector estimates (Ritz vectors) of . ⇤ : is a diagonal matrix of : eigenvalues of ) = that are eigenvalue estimates (Ritz values) of . Note that F 8 's are basis vectors of the same Krylov subspace K = ( , 1). Hence, we can replace & = with , : , and ) = with ⇤ : in Eq. ( 2). The TRL algorithm restarts the Lanczos iteration upon , : and ⇤ : , but by appending the vector @ =+1 .2 One problem is that F 8 's are not orthogonal to @ =+1 , so h @ =+1 |F 8 i terms survive in the Gram-Schmidt process generating @ =+2 . This results in a arrow-like matrix ) =+1 . However, thanks to the symmetry of , the generation of @ =+2+8 (8 1) returns to the normal three-term Lanczos recurrence. Hence, we need to deal only with diagonalizing the arrow-like matrix, which can be done in several ways including the general matrix diagonalization. Now that @ =+2 is orthogonal to F 1 , • • • , F : , approximate eigenvectors of , the restarted run searches on the Krylov subspace orthogonal to them. It focuses on finding other eigenvectors rather than F 8 's, though it still slowly improves F 8 's convergences as well. With a proper choice of :, it can balance the convergence of eigenvalues e ciently.
The QUDA library has a TRL eigensolver named TRLM. This eigensolver determines : as the number of converged eigenvalues in the desired precision plus the half of the remaining eigenvalues that are still converging. It also implements the locking method, which explicitly deflates some converged eigenvectors from the restarted run. TRLM locks eigenvectors converged to the machine precision. It reduces memory usage and computing cost.
In Fig. 3, we present QUDA TRLM's performance for various Krylov subspace dimensions = with Chebyshev acceleration of ? = 50, 100. In both Fig. 3a and Fig. 3b, we find TRL can perform better than the non-restarted Lanczos with optimization. One benefit of TRL, compared with IRL, is that its restarting cost can be absorbed into the Ritz vector calculation. That might be why TRL takes the least time for = = 600 with six restarts in Fig. 3b.
Restarted Lanczos algorithms (both TRL and IRL) allow us to use less memory. It is especially useful for GPU systems, where memory is usually lacking. Figure 3b does not include the nonrestarted Lanczos result for ? = 50 because it does not converge within the given GPU memory size.
The Block Lanczos (BL) algorithm runs the Lanczos iteration with multiple (D) starting vectors {1 1 , 1 2 , • • • , 1 D }, where 1 8 's are orthogonal to each other [5]. At each 9-th iteration of BL, a block version of the Lanczos recurrence Eq. ( 1) constructs D basis vectors @
The convergence of a BL iteration is usually slower than that of the basic Lanczos for the same subspace dimension = = =D, because the Krylov subspace dimension = per starting vector 1 8 is smaller than =. However, on a parallel computing system, BL may outperform the basic Lanczos by parallelizing the multiplications @ 8 9 1 for 8 = 1, • • • , D. The Split Grid method is one way of doing it [5]. It splits a communication grid such as MPI and distributes parallel jobs into the split grids, so that each job runs with a lower surface-to-volume ratio.
The Grid library has a BL eigensolver implementing the Split Grid method, named Implic-itlyRestartedBlockLanczos.3 In Fig. 4, we compare its performance with the unblocked basic Lanczos and BL without the Split Grid method. The block size for BL is set to D = 8. In this measurement, we run 8 MPI processes connected by an intra-node network. We calculate 500 eigenvalues, for which the unblocked Lanczos converges at = = 800, while the BL converges at = = 904. Hence, without the Split-Grid method, the BL algorithm itself is slower than the basic Lanczos algorithm, as the results show.
However, the Split Grid method can enhance the BL's performance. In Fig. 4, BL's performance is improved by 12% with the Split Grid method, even on the intra-node network. Although it is still a little slower than the unblocked Lanczos, its advantage could become significant on a slow inter-node network.4
We have discussed improved Lanczos algorithms and their optimizations. A well-tuned Chebyshev polynomial improves the Lanczos iteration's convergence significantly. All other improved 3Although it provides both non-restarted and implicitly restarted versions of BL, the latter is not pursued. 4It is under investigation. A preliminary result shows that with 8 MPI inter-nodes, BL with Split Grid gives around 200% better performance than the unblocked Lanczos.
Lanczos algorithms utilize it. Restarted Lanczos algorithms, such as Implicitly Restarted Lanczos and Thick-Restart Lanczos, do not always perform better than the non-restarted Lanczos algorithm. Still, we can optimize them to perform better than or comparable to non-restarted Lanczos. It is advantageous for small memory systems such as GPUs. The performance of Block Lanczos utilizing the Split Grid method is comparable to that of the unblocked Lanczos on an intra-node network. It may perform better on inter-node networks.
2One cannot start the Lanczos algorithm with an eigenvector, because it generates a rank-1 Krylov subspace.