Accelerating Eigenvalue Computation for
Nuclear Structure Calculations via Perturbative Corrections

Nuclear structure calculations require solving $A$-body Schrödinger equations where $A=Z+N$ is the number of nucleons consisting of $Z$ protons and $N$ neutrons. The Configuration Interaction (CI) method or no-core shell method (NCSM) [1], which represents the solution to the Schrödinger’s equation by a linear combination of $A$-body basis functions, reduces the problem to an algebraic eigenvalue problem:

where $H\in\mathbb{R}^{n\times n}$ is the matrix representation of the $A$-body nuclear Hamiltonian operator in a configuration space (spanned by a set of $n$ $A$-body basis functions), $E_{k}$ is the $k$th eigenvalue of $H$ representing an approximation to a discrete energy level, and $\Psi_{k}$ is the corresponding eigenvector that contains the coefficients of the $A$-body basis function in the expansion of the approximate eigenfunction in the $A$-body basis.

The dimension ($n$) of the matrix $H$ depends on the number of nucleons $A$ and the size of the CI model space (determined by a truncation parameter $N_{\max}$). Although $n$ can be quite large, $H$ is very sparse, and often only a few of its eigenpairs at the low end of the spectrum are of interest, making iterative methods suitable for solving (1).

The construction of the matrix $H$ by the CI method is typically done in a hierarchical fashion where the leading submatrix of $H$ corresponds to a matrix constructed from a smaller configuration space. Because the $A$-body basis functions that form the lower dimension CI space associated with a small $N_{\max}$ are typically more important than basis functions outside of such a configuration space, the eigenvectors of $H$ tend to be localized; i.e., the leading components of the $\Psi_{k}$ tend to be the larger in magnitude, while the tailing components are relatively small in magnitude. We illustrate these properties of a nuclear Hamiltonian and its wavefunction using nucleus ¹²C as an example. The Hamiltonian is constructed with a nucleon-nucleon interaction Daejeon16 [2] where hbar-omega value of $20$ describes the harmonic oscillator basis functions. Figure 1(a) shows that the leading submatrix of $H$ that is $63$ times smaller corresponds to a matrix constructed from a smaller configuration space. Figure 1(b) shows that the eigenvector corresponding to the lowest eigenvalue is localized in its first few components. Therefore, the vector formed by padding the eigenvector of the leading submatrix with zeros can serve as a good initial guess for many iterative methods used to solve large-scale eigenvalue problems (1). Previous works [3, 4] select such initial guesses for algorithms like the Lanczos algorithm [5], the block Lanczos algorithm [6], the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) algorithm [7], and the Residual Minimization Method with Direct Inversion of Iterative Subspace (RMM-DIIS) correction [8, 9, 10]. Another study [11] uses greedy algorithms to incrementally enlarge the submatrix and use the eigenvector of the enlarged submatrix as an improved starting guess.

In this paper, we propose solving the eigenvalue problem (1) with an iterative subspace projection method constructed from basis vectors obtained from eigenvalue and eigenvector perturbation theory. We call this method a Subspace Projection with Perturbative Corrections (SPPC) method. In this approach, the Hamiltonian matrix to be partially diagonalized is viewed as the sum of two matrices, i.e., $H=H_{0}+V$, where the eigenpairs of $H_{0}$ are relatively easy and inexpensive to compute, and $V$ is viewed as a perturbation to $H_{0}$. Perturbative corrections to eigenpairs of $H_{0}$ in successively higher order can be computed by solving a sequence of linear systems of equations. Together with the initial approximation to the desired eigenvector obtained from $H_{0}$, these correction vectors form a low dimensional subspace from which approximate eigenpairs of $H$ are extracted through the Rayleigh-Ritz procedure. When $H_{0}$ is chosen to be a block diagonal matrix diag($\hat{H}_{0}$,0), where $\hat{H}_{0}$ is the matrix representation of the $A$-body nuclear Hamiltonian in a smaller configuration space associated with a smaller $N_{\max}$ value, these linear systems can be solved efficiently in the small configuration space. The overall computational cost of the SPPC method grows linearly with respect to the highest order of perturbation included in the correction subspace. Adding each perturbative correction vector and performing the Rayleigh-Ritz calculation requires multiplying the sparse matrix $H$ with a single vector. We show numerically that the low dimensional subspace constructed in a low order SPPC method provides more accurate approximation to a few lowest eigenvalues of $H$ than a Krylov subspace of the same dimension constructed from the same starting guess. The method also appears to be more efficient than the locally optimal block preconditioned conjugate gradient (LOBPCG) method in early iterations even when a good preconditioner is available for the LOBPCG method. Although the SPPC method can stagnate as higher perturbative corrections are included, convergence stagnation can be mitigated by combining SPPC with other iterative algorithms for solving large scale eigenvalue problems using the SPPC’s eigenvector approximation as the starting guess for secondary algorithms.

The rest of this paper is organized as follows. In Section II we describe the basic Subspace Projection with Perturbative Corrections (SPPC) method for computing one eigenpair of $H$, and show how successively higher order perturbative corrections can be obtained from solutions of a set of linear systems of equations. We draw the connection between SPPC and previously developed eigenvector continuation methods. In Section III, we present a version of the SPPC algorithm that can be used to compute a few eigenpairs. A few practical implementation details of the SPPC method are discussed in Section IV. Numerical examples that demonstrate the efficiency of the SPPC method relative to other conventional large-scale eigenvalue computation methods are presented in Section V. We also show the effectiveness of combining SPPC with a conventional eigensolver for computing a few lowest eigenpairs of several nuclei.

We split the matrix $H$ into the sum of two matrices $H_{0}$ and $V$, i.e.,

where

with the matrix $\hat{H}_{0}\in\mathbb{R}^{n_{0}\times n_{0}}$ (where $n_{0}\ll n$) being the leading submatrix of $H$ that corresponds to the representation of the Hamiltonian within a smaller configuration space.

The eigenvectors of $H_{0}$, which can be obtained from the eigenvectors of $\hat{H}_{0}$, are computed with a much lower computational cost compared to those of the full matrix $H$. They can be used as good initial guesses for conventional algorithms, such as the Lanczos and the LOBPCG algorithms, for computing the eigenpairs of $H$.

The SPPC method uses the eigenvectors of $H_{0}$ to initiate a subspace construction procedure to produce a subspace from which improved approximations of the desired eigenpairs can be obtained.

Instead of using the Lanczos algorithm or the LOBPCG method to construct the subspace, we use perturbative corrections to the initial eigenvector approximation to construct the subspace. In this approach, we view $H$ as a matrix obtained from perturbing $H_{0}$ by $cV$ with $c=1$, i.e., we can write $H$ as

It follows from the Rayleigh-Schrödinger perturbation theory [12] that the eigenvalues and eigenvectors of $H(c)$, which we denote by $E_{k}(c)$ and $\Psi_{k}(c)$ respectively, can be written in terms of perturbative corrections to the eigenvalues and eigenvectors of $H_{0}$, which we denote by $E_{k}^{(0)}$ and $\Psi_{k}^{(0)}$, i.e.,

Here, ${(E_{k}^{(p)},\Psi_{k}^{(p)})}_{p\geq 1}$ represent the $p$th order perturbative corrections to the eigenpair $(E_{k}^{(0)},\Psi_{k}^{(0)})$, and are independent of the parameter $c$.

Substituting (5) into the equation

and matching coefficients of the same degree order yields the following set of equations

that allow us to compute $E_{k}^{(p)}$ and $\Psi_{k}^{(p)}$ in a recursive fashion.

The asymptotic expansion used in (5) assumes that $c$ is a small parameter. As a result, the expansion serves as a good approximation to the desired eigenpair only when $c$ is sufficiently small, i.e., when $c$ falls within the radius of convergence for (5), which is generally much smaller than 1. As a result, (5) cannot be used directly in general to approximate the $k$th eigenpair of $H$ [13, 14]. However, the perturbative vectors $\Psi_{k}^{(p)}$ can be used to construct a subspace

from which approximation to $E_{k}$ and $\Psi_{k}$ can be obtained.

The idea of using perturbative corrections to construct an approximating subspace was proposed in [13, 15] in the context of an eigenvector continuation (EC) method [16, 17]. In an EC method, the eigenvectors of $H(c)$ for some choices of $c$’s are used to construct a subspace from which approximations to the eigenvectors of $H(c^{\prime})$ for $c^{\prime}\neq c$ are obtained from the projection of $H(c^{\prime})$ into such a subspace.

It was found in  [13, 15] that instead of using eigenvectors of $H(c)$ for several choices of $c\neq 1$ to construct a subspace from which approximate eigenvectors of $H(1)$ are extracted through the standard Rayleigh-Ritz procedure, more accurate approximations to the desired eigenpairs of $H(1)$ can be obtained from the subspace constructed from the eigenvector of $H_{0}$ as well as perturbative eigenvector corrections as discussed above.

In a Rayleigh-Ritz procedure, we compute an orthonormal basis matrix $Q_{k}^{(P)}$ of $\mathcal{M}_{k}^{(P)}$ and form the projected matrix

If $(\theta,v)$ is an eigenpair of $\tilde{H}$, then $(\theta,z)$ where $z=Q_{k}^{(P)}v$, yields an approximate eigenpair of $H$. We consider an approximate eigenpair to have converged if its relative residual norm

is smaller than a preset tolerance value.

To obtain an approximate $k$th eigenpair $(\tilde{E}_{k},\tilde{\Psi}_{k})$ of $H$, we choose $(\theta,v)$ as the lowest eigenpair of $\tilde{H}$ in the above mentioned Rayleigh-Ritz procedure. We present the SPPC method in Algorithm 1.

We should point out that the main distinction between the SPPC method proposed here and the EC method lies in the subspaces constructed in these methods and the cost of construction. Because EC projects $H$ onto a subspace constructed from the eigenvectors of $H(c)$ for several (nonzero) $c$’s, the cost of subspace construction may be just as expensive as the cost of solving the target eigenvalue problem with a particular choice of $c$. On the other hand, because the subspace constructed in the SPPC method uses perturbative corrections that can be obtained by solving much smaller linear systems, the cost of subspace construction is significantly lower.

Although the basic idea of SPPC was presented in  [13, 15], the computational cost of this method was not carefully analyzed and compared with those of the state-of-the-art large-scale eigensolvers. In [13], the eigenvectors of $H_{0}$ can be computed analytically for the one-dimensional quartic anharmonic oscillator. However, in general, identifying a $H_{0}$ that can be diagonalized analytically is not possible. In [15], SPPC is used to perform a $A$-body nuclear structure calculation. However, a different $H=H_{0}+V$ splitting scheme is used, and the eigenvectors of $H_{0}$ are not easier to compute than those of $H$.

Although we can use Algorithm 1 to compute each of the first $k_{ev}$ eigenpairs one by one (or in parallel), approximations to larger eigenvalues and the corresponding eigenvectors appear to converge slowly, and sometimes to wrong values, as we will show in section VI. A more effective way to obtain approximations to the first $k_{ev}$ eigenpairs is to combine the SPPC subspace $\mathcal{M}_{k}^{(P)}$ (8) constructed for each eigenpair to create a larger subspace

from which $k_{ev}$ approximate eigenpairs can be extracted simultaneously through the Rayleigh-Ritz procedure, by targeting the $k_{ev}$ lowest eigenpairs of the projected matrix. Algorithm 2 shows how this approach works.

In this section, we describe a few practical implementation details of the SPPC algorithm.

As we will show in the next section, the perturbative correction is typically effective when the order of perturbation $p$ is relatively low. The convergence of SPPC can stagnate when $p$ increases, i.e., adding higher order perturbative correction may not help because they may be linearly dependent with respect to the basis vectors included in $\mathcal{M}^{(p)}$ already. In this case, it is useful to combine SPPC with another algorithm that can take the eigenvector approximation produced by SPPC as the starting guess.

Algorithms that can be combined with the SPPC method in a hybrid algorithm include but are not limited to the following:

• •

  The Lanczos algorithm, which is a classical algorithm that generates an orthonormal basis of a Krylov subspace using the Gram-Schmidt procedure. It is initialized with an approximate eigenvector produced from SPPC or a linear combination of $k_{ev}$ approximate eigenvectors. It uses one SpMV per iteration to compute the next basis vector. The eigenpairs are approximated using Ritz pairs obtained from the Rayleigh-Ritz procedure with the basis vectors of the Krylov subspace.

• •

  The block Lanczos algorithm, which is a variation of the Lanczos algorithm that operates in blocks. It is initialized with a block of vectors approximating several eigenvectors of $H$, and builds a Krylov subspace in blocks, with each iteration performing $k_{ev}$ SpMVs. The main advantage of the block Lanczos algorithm over the Lanczos algorithm is that it can make use of approximations to several eigenvectors more effectively and most dense linear algebra operations can take advantage of level 3 BLAS.

• •

  The LOBPCG algorithm, which is an iterative method that solves the equivalent trace minimization formulation of the eigenvalue problem. Similar to the block Lanczos algorithm, it can be initialized with approximations to several eigenvectors and each iteration performs $k_{ev}$ SpMVs. One advantage of the LOBPCG algorithm is that it can utilize a preconditioner if it is available. The use of a preconditioner can accelerate convergence. A common choice for the preconditioner is a block diagonal part of the Hamiltonian matrix or a shifted matrix of the Hamiltonian for some appropriately chosen shift [7, 4]. For our numerical experiments, we use a shifted preconditioner that involves a specific block diagonal part of the Hamiltonian matrix and a constant shift that approximates the lowest eigenvalue.

• •

  The residual minimization method (RMM) combined with the Direct Inversion of Iterative Subspace (DIIS) refinement algorithm (RMM-DIIS), which is a quasi-Newton algorithm for improving a specific eigenpair without computing other eigenpairs, provided that the initial approximation to the desired eigenpairs is sufficiently accurate. Each RMM-DIIS iteration performs one SpMV. The RMM-DIIS can also incorporate a preconditioner to accelerate convergence. For our numerical experiments, we choose a preconditioner that is a specific diagonal part of the Hamiltonian matrix.

In this section, we demonstrate the effectiveness of SPPC algorithm and compare it with other existing algorithms for computing the ground and a few low excited states of several light nuclei. We also show that SPPC can be effectively combined with RMM-DIIS to yield an efficient and accurate hybrid eigensolver for nuclear configuration interaction calculations. We call this hybrid eigensolver the SPPC+RMM-DIIS method.

For all algorithms tested in this section, we use the relative residual norm (10) as the stopping criteria and set the convergence tolerance to be $10^{-6}$. To ensure a fair comparison with the SPPC, we use the eigenvectors of the zero-order part $H_{0}$ as the initial guesses for the algorithms. All experiments were conducted using MATLAB.

In conclusion, the Subspace Projection with Perturbative Corrections (SPPC) method, combined with the Residual Minimization Method with Direct Inversion of Iterative Subspace (RMM-DIIS), presents a significant advancement in the efficient computation of eigenpairs for large Hamiltonian matrices in nuclear structure calculations. The SPPC method leverages perturbative correction vectors to enhance the accuracy of eigenpair approximations in the initial iterations, substantially reducing the number of sparse matrix-vector multiplications (SpMVs) required for convergence. Although the SPPC may experience stagnation in subsequent iterations, this challenge is effectively mitigated by integrating it with the RMM-DIIS algorithm, which provides robust refinement of eigenvector approximations. Our numerical experiments across several nuclear Hamiltonians demonstrate that the SPPC+RMM-DIIS hybrid approach outperforms traditional methods in terms of SpMVs. This hybrid method offers a promising solution for large-scale nuclear structure calculations, providing a reliable and efficient approach to solving the $A$-body Schrödinger equation.

While the preliminary results of the SPPC are promising, we have not provided a theoretical background explaining why it works. We plan to address this in our future work, discussing the convergence behavior in detail. From a practical standpoint, w e aim to develop a method to automatically detect stagnation, eliminating the need for manual decisions on when to switch to the RMM-DIIS. Additionally, we are interested in implementing the SPPC in a hybrid MPI/OpenMPI code, such as the software MFDn (Many-Fermion Dynamics for nuclear structure) [18, 19, 20], to be run at high-performance computing centers. We also want to explore further optimizations and applications of the SPPC to other large-scale eigenvalue problems in nuclear physics and beyond.