Stochastic extension of the Lanczos method for nuclear shell-model calculations with variational Monte Carlo method

We propose a new variational Monte Carlo (VMC) approach based on the Krylov subspace for large-scale shell-model calculations. A random walker in the VMC is formulated with the $M$-scheme representation, and samples a small number of configurations from a whole Hilbert space stochastically. This VMC framework is demonstrated in the shell-model calculations of $^{48}$Cr and $^{60}$Zn, and we discuss its relation to a small number of Lanczos iterations. By utilizing the wave function obtained by the conventional particle-hole-excitation truncation as an initial state, this VMC approach provides us with a sequence of systematically improved results.


I. INTRODUCTION
Large-scale shell model calculations in the nuclear structure study are performed by solving an eigenvalue problem of a large-dimension sparse matrix.The Krylovsubspace method is one of the best tools, and practically the only solution to solve the problem efficiently [1].The Krylov subspace is spanned by a starting vector v and the product of the first p powers of the Hamiltonian matrix H, namely {v, Hv, H 2 v, • • •, H p v}.An eigenvalue of the Hamiltonian matrix in Krylov subspace is called the Ritz value, and it is known to converge in a small number of iterations to the largest (or smallest) eigenvalue, since H p v with large p is dominated by the eigenvectors which have the large absolute eigenvalue.In the Krylovsubspace method, the converged Ritz value becomes a good approximation to the exact eigenvalue.Since the construction of the Krylov subspace requires only the matrix-vector product, the Krylov-subspace method has been extensively used and developed to solve an eigenvalue problem of a huge sparse matrix.
In the case of nuclear shell-model calculations, the Hamiltonian matrix in M -scheme basis is very sparse since the Hamiltonian consists of one-body and two-body interactions.The needed dimension of the Hamiltonian matrix is often quite huge, therefore, the Krylov-subspace iteration algorithm is quite efficient.The Lanczos algorithm, one of the most famous Krylov subspace algorithms, was introduced in 70's, [2][3][4] and has been widely used in shell-model calculations.Nowadays, it is implemented to take advantage of massively parallel computations [5].
Nevertheless, the application is still hampered by the exponential growth of the dimension of the Hilbert space.The size of a state vector surpasses the capacity of memory and the truncation of the model space is required.The most naive truncation is to assume a Fermi level for single-particle occupation and to restrict the number of particle-hole excitation across the Fermi level up to t particles.It is called the t-particle t-hole truncation and frequently used in practical calculations.As t increases, the eigenstate in the truncated subspace approaches the true eigenstate rather gradually.However, the significance of the large t component remains and is difficult to estimate (e.g.[6]) due to the limitation of computational resources.
On the other hand, much effort has been paid to circumvent this problem by introducing the Markov Chain Monte Carlo (MCMC) in the context of quantum Monte Carlo [7,8].However, the notorious "sign problem" prohibits to use realistic shell-model interaction practically.
In the present work, to estimate the omitted contribution of t-particle t-hole truncation for the exact shellmodel energy in full Hilbert space we introduce the variational Monte Carlo (VMC) approach into the shell-model calculations, and show that this VMC can overcome the limitation of the truncation scheme without treating the basis vectors of the full Hilbert space explicitly.

II. FORMULATION
We describe the form of a trial wave function of this VMC approach.At the beginning, we calculate the lowest eigenstate in the t-particle t-hole truncated subspace, |ψ t .Then, we introduced the VMC approach combined with the Krylov-subspace method with |ψ t being a starting vector to improve the approximation systematically.
The Ritz vector, or the approximated eigenvector on the Krylov subspace is taken as a trial wave function such as where H is the Hamiltonian and c = {c 0 , c The energy expectation value of the trial wave function is written by inserting the complete set of the M -scheme subspace, 1 = m∈{M π } |m m|, where |m has a fixed z-component of angular momentum M and parity π, such as Since the number of M -scheme states rapidly increase as a function of the particle number and the size of the single-particle space, the practical summation of |m in Eq.( 2) often becomes difficult to perform.Therefore, we apply the Monte Carlo technique to this summation and calculate the energy expectation value E(c) stochastically.The Monte Carlo random walker of m state is formulated in the M -scheme basis obeying the probability proportional to | m|Ψ(c) | 2 .The complete set of |m is represented by the N samples of the random walkers as with the local energy and the sampling density ( The m∈M.C. denotes the summation of N samples, m, which are generated by the M -scheme random walker with the probability density ρ c (m).
In order to compute the E(c) in Eq.( 3) stochastically, we briefly describe the MCMC process to generate a random walker in M -scheme basis states.It was first introduced in Ref. [9].
The transition of the random walker is controlled by the Metropolis algorithm.The candidate of the transition is created by two-particle two-hole operator like being the creation operator of a single-particle state i.The indices, i, j, k, and l, are restricted so that |m and |m ′ having the same zcomponent of angular momentum and parity.Whether or not a random walker |m moves to |m ′ depends on the ratio p c (m ′ ) = ρ(m ′ )/ρ(m) as If p c (m ′ ) ≥ 1, the walker |m always moves to |m ′ .If p c (m ′ ) < 1, according to the p c (m ′ ), we determine whether or not the walker |m moves to |m ′ .This procedure satisfies the detailed balance condition and ergodicity.
The M -scheme walker is automatically restricted to good M , parity, z-component of isospin subspace because the initial wave function |ψ t is already an eigenstate of these quantum numbers, and the sampling density not having suitable values of M and parity becomes zero.In addition, because the |ψ t is an eigenstate of the total angular momentum J 2 and the Hamiltonian commutes with the J 2 , the expectation value of J 2 of the resultant state in Eq.( 1) is the same as that of the |ψ t without statistical error.Therefore, unlike the formulation of the VMC in Ref. [9], we do not need angular-momentum projection which is represented by three-dimensional integral over the Euler angles and causes a time-consuming numerical computation.
In practical calculations, we implemented a cache algorithm for efficient computation.Namely, we store any calculated matrix elements m|H q |ψ t with 0 ≤ q ≤ p in memory, and reuse them if needed.This prescription shortens the computational time drastically, however it requires a large amount of memory in compensation.This is a trade off between the computational time and the memory usage, and the introduction of an efficient cache algorithm would ease this problem.
Here, we discuss the relation with the power method [10,11], which has a simple form of the trial wave function and has been used widely.In the power method, the wave function is approximated by with a constant value Λ.This wave function is equal to Eq.( 1) if c is taken as a set of binomial coefficients.Since the VMC has a larger variational space than the power method with the same p, the VMC provides a better approximation.Equation( 7) is extended to p q=1 (Λ q − H)|ψ t , which is equivalent to Eq.( 1).The optimization of the variational parameters requires negligible computational cost compared to that of calculation without variation by utilizing the reweighting method discussed below.
For the energy minimization with respect to the variational parameters c, the Nelder-Mead downhill simplex method is utilized [12].In this method, the energy functional is minimized by morphing the polytope consisting of p + 1 vertexes in the p-dimension parameter space.No other information of the energy functional is required, e.g.energy gradient.Thus we have to obtain the energy expectation values with many sets of samples with various c's.We describe the reweighting method [13] to determine the best variational parameters with a single set of samples for a certain set of c.
We suppose that a set of Monte Carlo samples and E(c 0 ) has already been obtained with appropriate initial parameters c 0 .Then, we calculate the energy expectation value E(c ′ ) where c ′ is close enough to c 0 .The E(c ′ ) is written as By using the random walker of ρ c0 (m), it is evaluated stochastically using with the reweighting factor Note that m∈M.C. in Eq.( 8) denotes the Monte Carlo summation of the M -scheme random walkers, of which the sampling density is ρ c0 (m), not ρ c ′ (m).In practice, when we evaluate the E(c 0 ) we generate a set of samples with the density distribution ρ c ′ (m) and keep all the computed matrix elements, m|H q |ψ t .By reusing this random walker and the matrix elements, we do not need additional computations to compute the E(c ′ ).

III. NUMERICAL RESULTS
To present the feasibility of the present method, we demonstrate the shell-model calculation of 48 Cr in pf shell as an example.Figure 1 shows the energy expectation values of the ground-state energy of the 48 Cr with GXPF1A interaction [14].Its M -scheme dimension is 1,963,461.We generate 48 random walkers utilizing the Metropolis algorithm.For each walker, after we run 1000 steps as burn-in process, a random walker moves 10000 steps in the M -scheme space.The variational parameters c in Eq.( 1) were determined to minimize the energy by the Nelder-Mead method and the reweighting technique discussed as before.
Figure 1 also shows the Ritz value of the Lanczos method with p-step iterations and |ψ t being the initial vector.The p and t in Fig. 1(a) denote the power p and the number of particle-hole truncation of the wave function |ψ t in Eq.( 1).The exact shell-model energy is −99.578MeV, which is almost reached at p = 3, 4 and t ≥ 3.If the variational parameters of the VMC were determined without statistical error, this trial wave function would equal the wave function obtained by the Lanczos method with p-iterations.More generally, this is a Monte Carlo formulation of Krylov subspace technique because the trial wave function is a linear combination of H q |ψ terms with 0 ≤ q ≤ N .The VMC results (error bars) agree quite well with those of Lanczos method (lines), which means that the Nelder-Mead optimization with reweighting works quite well.While the energy of the initial state |ψ t , or p = 0 error bars in Fig. 1(a), does not converge well as a function of t, the VMC values of p = 3, 4 converges well even for the small t.
Since the Hamiltonian only contains two-body interactions, H p |ψ t is considered to be a state in (t + 2p)particle (t + 2p)-hole truncated subspace.Figure 1(b) shows the energy as a function of t + 2p.The left ends of these lines on Fig. 1(b) correspond to the exact solution with t-particle t-hole truncation, and therefore the variational lower bound of the truncated subspace.With increasing p, the energy converges well and close to the Lanczos value.Physical observables other than the energy, e.g.moments and transition probabilities, are also obtained in a similar manner to Ref. [9]. Figure 2 shows the B(E2) transition probability obtained by the VMC.The effective charge is (e p , e n ) = (1.2,0.8)e.The convergence behaves similar to the case of the energy in Fig. 1(a) and the VMC dramatically improve the E2 values in the region of small t.
For the demonstration of larger-scale calculations, we discuss the case of 60 Zn in pf shell with the GXPF1A interaction [14].The M -scheme dimension is huge, 2.2 × 10 9 , though it is somehow tractable by the recent shellmodel code [15].Figure 3 shows the energy expectation values of the ground state obtained by the VMC and the corresponding Lanczos method.Though the dimension is about 10 3 times larger than that of 48 Cr, the energy convergence is similar to Fig. 1. Figure 4 shows the excitation energies obtained by the VMC approach and corresponding p-iterated Lanczos method.The excitation energies converged far faster than the energy itself, and the VMC significantly improves the convergence.

IV. SUMMARY
In summary we proposed a new VMC formalism to improve the conventional particle-hole truncation shellmodel calculations systematically, and demonstrated its feasibility at 60 Zn in pf In this formalism, we do not have the sign problem because it is implemented by the variational Monte Carlo, namely the density probability is defined by the wave function squared.However, the power of the Hamiltonian in the Monte Carlo implementation of the Lanczos method is restricted to be rather small, such as p = 3 or 4. Nevertheless, it provides us with good convergence of the energy and is expected to be useful to go beyond the conventional diagonalization method.The energy expectation value of our p-power VMC scheme agrees with that of p-th step Lanczos method in full space within statistical error.It means that the MCMC procedure and reweighting method to determine the variational parameters work well and stable.The present VMC approach provides us with the facile estimation of the exact energy eigenvalue by using the t-particle t-hole truncated wave functions.The energy variance can also be calculated stochastically with the same formulation, which is expected to be helpful for the energy-variance extrapolation technique [16], and it remains for future study.
The cache algorithm to store m|H q |φ t with a random walker m drastically reduces the computation time At present, since we store all matrix elements on memory, it requires a large amount of memory usage.This difficulty can be eased by sophisticated cache algorithm, the implementation of which is important for further applications.

FIG. 1 :
FIG. 1: (Color online) Energy expectation values of the J π = 0 + trial wave functions of 48 Cr in pf shell.(a) the blue, green, purple, and red error bars denote the VMC results with p = 0, 1, 2, 3, and 4, respectively.(b) the blue, green, purple, and red error bars denote the VMC results with t = 0, 1, 2, 3, and 4 as functions of t + 2p.The lines denote the results of the corresponding p-iteration Lanczos method.

FIG. 3 :
FIG. 3: (Color online) Energy expectation values of the J π = 0 + trial wave functions of 60 Zn in pf shell.See the caption of Fig.1.