Computing the lowest eigenvalues of the Fermion matrix by subspace iterations

TL;DR

The paper presents a subspace iteration method for efficiently computing multiple lowest eigenvalues of large Hermitean matrices, showing that increasing the number of eigenvalues does not significantly raise computational costs.

Contribution

It introduces a subspace iteration approach that efficiently computes multiple eigenvalues simultaneously without substantial additional computational effort.

Findings

Computational cost remains nearly constant as the number of eigenvalues increases.

Multiple eigenvalues can be obtained at almost no extra cost compared to computing a single eigenvalue.

The method is effective for large, sparse Hermitean matrices.

Abstract

Subspace iterations are used to minimise a generalised Ritz functional of a large, sparse Hermitean matrix. In this way, the lowest $m$ eigenvalues are determined. Tests with $1 \leq m \leq 32$ demonstrate that the computational cost (no. of matrix multiplies) does not increase substantially with $m$. This implies that, as compared to the case of a $m=1$, the additional eigenvalues are obtained for free.

Tables

Table 1. Table 1: Number of A ​ v 𝐴 𝑣 Av operations to obtain m 𝑚 m lowest eigenvalues with precision δ = 10 − 6 𝛿 superscript 10 6 \delta=10^{-6} on three random SU(2) lattices

m	$4^4$	$6^4$	$8^4$
1	2992	28500	29966
2	2816	11911	40075
4	4011	24414	37521
8	4921	19322	50947
16	5268	21494	64801
32	6000	25089	66994

Equations

$$q(x)=\frac{x^{†}Ax}{x^{†}x},x\in C^N.$$

$$P(x)=x(x^{†}x{)}^{-1}{x}^{†}.$$

$$q(x)=Tr\{P(x)A\}.$$

$$g=(1-P(x))Ax(x^{†}x{)}^{-1}\in C^{N\times m}.$$

$$x^{\prime }=x+h\alpha$$

$$x^{†}h=0,$$

$$\displaystyle\left(\begin{array}[]{cc}x^{\dagger}Ax&x^{\dagger}Ah\\ h^{\dagger}Ax&h^{\dagger}Ah\end{array}\right)\left(\begin{array}[]{c}\beta_{1}\\ \beta_{2}\end{array}\right)\quad\quad\quad\quad$$

$$\displaystyle=\left(\begin{array}[]{cc}x^{\dagger}x&0\\ 0&h^{\dagger}h\end{array}\right)\left(\begin{array}[]{c}\beta_{1}\\ \beta_{2}\end{array}\right)\Lambda.$$

$$\left(\begin{array}[]{c}\beta_{1}\\ \beta_{2}\end{array}\right)\in C^{2m\times m}.$$

$$\alpha ={\beta }_2{\beta }_1^{-1}.$$

$$h^{\prime }=g^{\prime }+(1-P(x^{\prime }))h\gamma$$

$$\gamma =(g^{†}g{)}^{-1}\left\{g^{\prime †}{g}^{\prime }-g^{†}{g}^{\prime }\right\}$$

$$q(x^{(n)})\approx q(x^{(\infty )})+a{f}^n,$$

$$q(x^{(n)})-q(x^{(\infty )})\approx \frac{1}{1-f}[q(x^{(n)})-q(x^{(n+1)})]$$

$$\frac{1}{1-f}|{\lambda }_i^{(n)}-{\lambda }_i^{(n+1)}|<\delta |{\lambda }_i^{(n+1)}|.$$

Full text

1 Introduction

For a large, sparse Hermitean matrix $A\in C^{N\times N}$, minimisation of the Ritz functional

[TABLE]

provides the smallest eigenvalue, together with the corresponding eigenvector. The minimum of a functional can be found iteratively by the method of conjugate gradients (CG)[1], its application to the special case of the Ritz functional has been worked out by Geradin[2] and Fried[3].

If a small number of lowest eigenvalues is required instead, $q(x)$ may be minimised repeatedly, restricting $x$ to the space orthogonal to the previous eigenvectors. An efficient implementation of this idea is described in [4].

2 Subspace iterations

An alternative approach to compute the $m$ smallest eigenvalues (with eigenvectors) simultaneously considers the corresponding subspace as a whole. An $m$–dimensional subspace of $C^{N}$ is spanned by $m$ non–degenerate (column) vectors, which are combined into a rectangular matrix $x\in C^{N\times m}$. From now on, $N$ is considered ‘large’, $10^{5}$ say, but $m$ is ‘small’, e.g. $O(10)$.

The projector onto the subspace is

[TABLE]

Note that the inversion of $x^{\dagger}x\in C^{m\times m}$ is a small problem and can be solved by standard techniques. Minimisation is now applied to the generalised Ritz functional[5]

[TABLE]

At its minimum, $P(x)$ projects onto the space of the $m$ lowest eigenvalues.

An iterative CG minimisation of $q(x)$ starts with the computation of the gradient at $x$:

[TABLE]

By construction, $x^{\dagger}g=0$. A steepest descent method would search a new minimum along $g$, but CG algorithms use an ansatz with a more general shift ‘direction’ $h\in C^{N\times m}$ (to be discussed later):

[TABLE]

with a (small) coefficient matrix $\alpha\in C^{m\times m}$. It is helpful to assume

[TABLE]

i.e. all shift vectors (columns of $h$) are orthogonal to the current subspace (spanned by the columns of $x$). The matrix $\alpha$ is to be determined to minimise $q(x+h\alpha)$. To bring this to a manageable form, consider the projection of $A$ into the space spanned by the columns of $x$ and $h$. This defines a $(2m)$–dimensional eigenvalue problem (in a non–orthogonal basis), which is to be solved for the $m$ lowest eigenvalues. The Jacobi method[1] can deal with such a ‘small’ problem. In a matrix notation using $m\times m$ blocks, the result can be expressed as

[TABLE]

The eigenvalues are in $\Lambda=\mbox{diag}(\lambda_{1}\ldots\lambda_{m})$ and the eigenvectors are the columns of

[TABLE]

The new minimising subspace is spanned by $x\beta_{1}+h\beta_{2}$. Comparision with eq.(5) reveals that the minimum of $q(x+h\alpha)$ is found at

[TABLE]

3 Update of the search space

Following the general CG scheme, one starts with $h=g$, but subsequently $h$ is updated to a superposition of the actual gradient and the previous search direction. An appropriate ansatz is

[TABLE]

($g^{\prime}$ is the gradient evaluated at $x^{\prime}$). The projection in front of $h$ insures the orthogonality eq.(6) for $h^{\prime}$. $\gamma\in C^{m\times m}$ is another (small) matrix. Its choice is somewhat arbitrary: in the linear case (minimisation of a quadratic form), it is determined by the requirement that the search vectors should be conjugate, but this argument does not carry over to the nonlinear situation. Experience has shown that the following form

[TABLE]

is reasonably efficient and robust. It is a matrix version of the Polak–Ribiere prescription[1]. Note that matrices like $g^{\dagger}g\in C^{m\times m}$ are small.

4 Convergence and stopping

By construction, $q(x)$ will decrease monotonically with the iteration number $n$, but the rate of convergence and the behavior of individual eigenvalues is fairly irregular and errors have to be estimated with care. As an example, Fig. 1 shows the relative shifts of the four lowest eigenvalues of the fermion matrix $M_{f}^{\dagger}M_{f}$ for Wilson fermions on an $8^{4}$ random $SU(2)$ gauge configuration.

The following procedure is proposed: assuming geometric convergence

[TABLE]

$f$ is computed from three widely separated $n$’s. Then the relation

[TABLE]

shows how ‘local’ shifts are to be rescaled to provide error estimates: an eigenvalue $\lambda_{i}$ is considered converged within relative error $\delta$ if

[TABLE]

The corresponding eigenvector is frozen and only used in further iterations to keep the remaining subspace orthogonal to it.

For the example given above, Fig. 2 demonstrates that this procedure leads to correct stopping as the eigenvalues approach the precise values within a prescribed relative error of $\delta=10^{-6}$.

5 Computational cost

Subspace iterations require to store the large arrays $x$, $h$, $Ax$, and one more work matrix, this amounts to $4mN$ variables in total.

As with many sparse matrix algorithms, the number of multiplications $Av$ (with a single column vector $v\in C^{N}$) is considered a fair measure for the cost in cpu time. To estimate the performance of the subspace iterations, tests with realistic matrices have been performed: Wilson fermions in four dimensions with random $SU(2)$ gauge fields and (critical) hopping parameter $\kappa=0.25$. In this case, $A\sim M_{f}^{\dagger}M_{f}$. For lattice sizes $4^{4},6^{4},8^{4}$, $m=1,2,4,8,16,32$ eigenvalues and a relative error bound $\delta=10^{-6}$, the numbers of operations $Av$ needed for all eigenvalues to converge are shown in Table 1.

The absolute number of iterations needed to achieve a given precision grows with the matrix size (density of the eigenvalues), but the behavior is difficult to understand in detail. Since diagonalisation within the subspace is done non–iteratively, one might speculate that the level spacings around the $m$–th eigenvalue are crucial. This requires further studies.

In any event, the counts in Table 1 show one encouraging feature: the number of $Av$ operations does not grow substantially with $m$. The $Av$ part is due to one computation of $Ah$ per iteration and costs $\sim Nm$. It will be dominated for larger $m$ by the computation of small matrices like $x^{\dagger}x$ ($\sim Nm^{2}$ per iteration), but these scalar products are much cheaper and can be neglected for the moderate values of $m$ considered.