An accelerated linear method for optimizing non-linear wavefunctions in variational Monte Carlo

TL;DR

This paper introduces an accelerated linear optimization method for non-linear wavefunctions in variational Monte Carlo, combining iterative eigenvalue solving with adaptive optimization strategies to enhance efficiency and scalability.

Contribution

It presents a novel approach using the Jacobi-Davidson algorithm and adaptive switching between AMSGrad and the linear method for improved optimization of complex wavefunctions.

Findings

Reduced memory usage and improved scaling with parameters.

Effective optimization across diverse atomic, molecular, and model systems.

Demonstrated efficiency gains over traditional methods.

Abstract

Although the linear method is one of the most robust algorithms for optimizing non-linearly parametrized wavefunctions in variational Monte Carlo, it suffers from a memory bottleneck due to the fact at each optimization step a generalized eigenvalue problem is solved in which the Hamiltonian and overlap matrices are stored in memory. Here we demonstrate that by applying the Jacobi-Davidson algorithm, one can solve the generalized eigenvalue problem iteratively without having to build and store the matrices in question. The resulting direct linear method greatly lowers the cost and improves the scaling of the algorithm with respect to the number of parameters. To further improve the efficiency of optimization for wavefunctions with a large number of parameters, we use the first order method AMSGrad far from the minimum as it is very inexpensive, and only switch to the direct linear method near the end of the optimization where methods such as AMSGrad have long convergence tails. We apply this improved optimizer to various wavefunctions with both real and orbital space Jastrow factors for atomic systems such as Beryllium and Neon, molecular systems such as the Carbon dimer and Iron(II) Porphyrin, and model systems such as the Hubbard model and Hydrogen chains.

Tables

Table 1. Table 1: Real space Beryllium atom optimization using a J- K ​ S z 𝐾 subscript 𝑆 𝑧 KS_{z} GHF ansatz, in a Slater TZP all electron basis Van Lenthe and Baerends ( 2003 ) . Energies are in units of Hartrees, and time is in seconds. The same number of samples were used for each macroiteration to ensure an error of ∼ 0.5 similar-to absent 0.5 \sim 0.5 mHa was reached for all calculations. Aggressive AMSGrad parameters were used consisting of α = 0.01 , β 1 = 0.1 , β 2 = 0.01 formulae-sequence 𝛼 0.01 formulae-sequence subscript 𝛽 1 0.1 subscript 𝛽 2 0.01 \alpha=0.01,\beta_{1}=0.1,\beta_{2}=0.01 . As a comparison, the energy from exact diagonalization in an aug-cc-pVQZ basis is − 14.640 14.640 -14.640 Ha.

Optimizer	Energy	Iteration	Wall Time	Time/Iteration
LM	$-14.667$	$21$	$123.84$	$5.90$
aLM	$-14.667$	$8$	$23.52$	$2.94$
AMSGrad	$-14.666$	$181$	$128.61$	$0.71$

Table 2. Table 2: Optimized real space Neon atom energies with a J- K ​ S z 𝐾 subscript 𝑆 𝑧 KS_{z} GHF ansatz, in a Slater DZ all electron basis Van Lenthe and Baerends ( 2003 ) . Energies are in units of Hartrees, and time is in seconds. dTol is shorthand for tolerance of the davidson solver, and likewise for cgTol with respect to the conjugate gradient solver. The same number of samples were used for each macroiteration to ensure an error of ∼ 0.5 similar-to absent 0.5 \sim 0.5 mHa was reached for all calculations. As a comparison, the energy from exact diagonalization in an aug-cc-pV5Z basis is − 128.900 128.900 -128.900 Ha.

Optimizer	Energy	dTol	cgTol	Time/Iteration
aLM	$-128.907$	${10}^{-3}$	${10}^{-2}$	$555$
aLM	$-128.908$	${10}^{-3}$	${10}^{-3}$	$664$
aLM	$-128.909$	${10}^{-4}$	${10}^{-3}$	$761$
LM	$-128.910$	-	-	$2489$

Table 3. Table 3: Optimized orbital space C 2 subscript 𝐶 2 C_{2} molecule energies with a J- K ​ S z 𝐾 subscript 𝑆 𝑧 KS_{z} GHF ansatz, in a 6-31G all electron basis. Lowdin’s orbital localization scheme was used to define the sampled Hilbert space Löwdin ( 1950 ) . Energies are in units of Hartrees, and time is in seconds. dTol is shorthand for tolerance of the davidson solver, and likewise for cgTol with respect to the conjugate gradient solver. Aggressive parameters were used for SR and AMSGrad, with a stepsize of 0.1 and diagonal shift of 0.001 for SR, and α = 0.01 , β 1 = 0.1 , β 2 = 0.01 formulae-sequence 𝛼 0.01 formulae-sequence subscript 𝛽 1 0.1 subscript 𝛽 2 0.01 \alpha=0.01,\beta_{1}=0.1,\beta_{2}=0.01 for AMSGrad. The initial energy was − 75.470 75.470 -75.470 Ha, and the same number of samples were used for each macroiteration to ensure an error of ∼ 0.5 similar-to absent 0.5 \sim 0.5 mHa was reached for all calculations. As a comparison, the energy from exact diagonalization in the same basis is − 75.643 75.643 -75.643 Ha.

Optimizer	Energy	dTol	cgTol	Iteration	Time/Iteration
aLM	$-75.630$	${10}^{-3}$	${10}^{-2}$	$24$	$76.2$
aLM	$-75.630$	${10}^{-3}$	${10}^{-3}$	$16$	$84.9$
aLM	$-75.630$	${10}^{-4}$	${10}^{-3}$	$17$	$113.3$
LM	$-75.630$	-	-	$23$	$1485.7$
AMSGrad	$-75.630$	-	-	$397$	$2.5$
SR	$-75.630$	-	-	$319$	$3.8$

Table 4. Table 4: H 50 subscript 𝐻 50 H_{50} chain optimization using the orbital space J- K ​ S z 𝐾 subscript 𝑆 𝑧 KS_{z} GHF ansatz, in a STO-6G basis. Lowdin’s orbital localization scheme was used to obtain the localized orbitals for the Jastrow factor Löwdin ( 1950 ) . Energies are in units of Hartrees, and time is in seconds. The same number of samples were used for each macroiteration to ensure an error of ∼ 0.5 similar-to absent 0.5 \sim 0.5 mHa was reached for all calculations. Aggressive parameters were used for SR and AMSGrad, with a stepsize of 0.1 and diagonal shift of 0.001 for SR, and α = 0.01 , β 1 = 0.1 , β 2 = 0.01 formulae-sequence 𝛼 0.01 formulae-sequence subscript 𝛽 1 0.1 subscript 𝛽 2 0.01 \alpha=0.01,\beta_{1}=0.1,\beta_{2}=0.01 for AMSGrad. For aLM the loosest tolerances corresponding to dTol = 10 − 3 dTol superscript 10 3 \text{dTol}=10^{-3} , cgTol = 10 − 2 cgTol superscript 10 2 \text{cgTol}=10^{-2} were used. 25 AMSGrad iterations were done to start off the aLM run, using the less aggressive parameters (this is what results in the differing Avg. Time per Iteration values for the aLM), the energy at which LM iterations start is tabulated with an asterisk.

Optimizer	Energy	Iteration	Avg. Time per Iteration
SR	$-26.886$	$247$	$108.3$
AMSGrad	$-26.886$	$289$	$104.7$
AMSGrad	$-26.905$	$1146$	$104.5$
AMSGrad	$-26.908$	$1898$	$104.6$
aLM*	$-26.713$	$25$	$109.4$
aLM	$-26.886$	$36$	$335.4$
aLM	$-26.908$	$69$	$622.2$
aLM	$-26.910$	$89$	$796.9$

Table 5. Table 5: 98-site Hubbard model optimization using a J- K ​ S z 𝐾 subscript 𝑆 𝑧 KS_{z} GHF ansatz. Energies are in units of t 𝑡 t , and time is in seconds. Errors on energies are ∼ 0.001 similar-to absent 0.001 \sim 0.001 t 𝑡 t . The loosest tolerances corresponding to dTol = 10 − 3 dTol superscript 10 3 \text{dTol}=10^{-3} , cgTol = 10 − 2 cgTol superscript 10 2 \text{cgTol}=10^{-2} were used for aLM. 75 AMSGrad iterations were done to start off the aLM run, using the less aggressive parameters.

Iteration	Energy
$0$	$-113.032$
$75$	$-116.258$
$93$	$-116.824$

Table 6. Table 6: Fe(II) porphyrin energy results in Hartrees for the A 1 ​ g 5 superscript subscript A 1 𝑔 5 {}^{5}\text{A}_{1g} and B 1 ​ g 3 superscript subscript B 1 𝑔 3 {}^{3}\text{B}_{1g} states in a cc-pVDZ basis. RHF is restricted open-shell Hartree Fock. A HCI-SCF calculation was performed in the active space with a tight tolerance to get a highly accurate variational result. Semi-stochastic perturbation theory was performed on top of the HCI wavefunction to obtain a nearly exact, albeit non-variational, result in the active space. The active space orbitals from the HCI-SCF calculation were localized in an intrinsic bonding orbital scheme Knizia ( 2013 ) , and then used as the basis for the VMC and GFMC calculations with a J- K ​ S z 𝐾 subscript 𝑆 𝑧 KS_{z} GHF ansatz. Sufficient samples were taken to obtain an error of less than 0.5 mHa.

Method	${}_{}{}^5\text{A}_{1g}^{}$	${}_{}{}^3\text{B}_{1g}^{}$	$\Delta \left({}_{}{}^5\text{A}_{1g}^{}-{}_{}{}^3\text{B}_{1g}^{}\right)$
RHF	-2244.753	-2244.394	-0.359
VMC	-2245.025	-2244.998	-0.027
GFMC	-2245.029	-2245.002	-0.027
Exact(SHCI)	-2245.029	-2245.002	-0.027

Equations

$$|\Psi ⟩=\hat{J}\hat{P}|\Phi ⟩,$$

$$|\text{GHF}⟩=\prod _{k=1}^{n_{\text{elec}}}{\hat{a}}_k^{†}|0⟩,$$

$${\hat{a}}_k^{†}=\sum _{l=1}^{n_{\text{orb}}}\sum _{\sigma =\uparrow ,\downarrow }{\theta }_{l\sigma ,k}\cdot {\hat{c}}_{l\sigma }^{†}.$$

$$\hat{J}=\exp \sum _{p\ge q}^{n_{orb}}\sum _{\sigma ,\gamma =\uparrow ,\downarrow }{J}_{p\sigma ,q\gamma }{\hat{n}}_{p\sigma }{\hat{n}}_{q\gamma },$$

$$\hat{J}=\exp \sum _{i,I}f(r_{iI})+\sum _{i>j}g(r_{ij})+\sum _{i,j,I}h(r_{ij},r_{iI},r_{jI}),$$

$$|{\bar{\Psi }}_i⟩=\{\begin{array}{ll}|\Psi ⟩,& \text{if}i=0\\ |{\Psi }_i⟩-|\Psi ⟩⟨\Psi ||{\Psi }_i⟩,& \text{if}i\ne 0\end{array}$$

$$\bar{\mathbf{H}}\cdot \delta \mathbf{p}=E\cdot \bar{\mathbf{S}}\cdot \delta \mathbf{p},$$

$$\delta p_i^{\prime }=\frac{\delta p_i}{1-{\sum }_{i=1}^{n_{var}}{N}_i\cdot \delta p_i},$$

$$N_i=⟨\Psi ||{\Psi }_i⟩,$$

$$N_i=-\frac{(1-\xi )\cdot {\sum }_j^{nonlin}{\bar{\mathbf{S}}}_{ij}\cdot \delta p_j}{(1-\xi )+\xi \cdot \sqrt{1+{\sum }_{ij}^{nonlin}\delta p_i\cdot {\bar{\mathbf{S}}}_{ij}\cdot \delta p_j}}.$$

$${\bar{\mathbf{H}}}_{ij}$$

$$=\sum _n\frac{⟨{\bar{\Psi }}_i||n⟩}{⟨\Psi ||n⟩}\cdot \frac{⟨n|\hat{H}|{\bar{\Psi }}_j⟩}{⟨n||\Psi ⟩}\cdot |⟨n||\Psi ⟩{|}^2$$

$${\bar{\mathbf{H}}}_{ij}=\text{expectationvalue}{\frac{⟨{\bar{\Psi }}_i||n⟩}{⟨\Psi ||n⟩}\cdot \frac{⟨n|\hat{H}|{\bar{\Psi }}_j⟩}{⟨n||\Psi ⟩}}_{{\rho }_n},$$

$${\bar{\mathbf{S}}}_{ij}=\text{expectationvalue}{\frac{⟨{\bar{\Psi }}_i||n⟩}{⟨\Psi ||n⟩}\cdot \frac{⟨n||{\bar{\Psi }}_j⟩}{⟨n||\Psi ⟩}}_{{\rho }_n}.$$

$${\mathbf{g}}_i[n]=\frac{⟨n||{\Psi }_i⟩}{⟨n||\Psi ⟩},$$

$${\mathbf{h}}_i[n]=\frac{⟨n|\hat{H}|{\Psi }_i⟩}{⟨n||\Psi ⟩}.$$

$$\left(\begin{array}{cc}{E}_0& {\mathbf{G}}_r^T\\ {\mathbf{G}}_c& \mathbf{H}\end{array}\right)\left(\begin{array}{c}1\\ \delta \mathbf{p}\end{array}\right)=E\left(\begin{array}{cc}1& 0^T\\ 0& \mathbf{S}\end{array}\right)\left(\begin{array}{c}1\\ \delta \mathbf{p}\end{array}\right),$$

$${{\mathbf{G}}_{\mathbf{r}}}_i=\text{expectationvalue}{\mathbf{h}}_i[n]-E_0\cdot \text{expectationvalue}{\mathbf{g}}_i[n]$$

$${{\mathbf{G}}_{\mathbf{c}}}_i=\text{expectationvalue}{E}_L[n]\cdot {\mathbf{g}}_i[n]-E_0\cdot \text{expectationvalue}{\mathbf{g}}_i[n].$$

$${\mathbf{S}}_{ij}=\text{expectationvalue}{\mathbf{g}}_i[n]\cdot {\mathbf{g}}_j[n]-\text{expectationvalue}{\mathbf{g}}_i[n]\cdot \text{expectationvalue}{\mathbf{g}}_j[n],$$

$${\mathbf{H}}_{ij}=\text{expectationvalue}{\mathbf{g}}_i[n]\cdot {\mathbf{h}}_j[n]-\text{expectationvalue}{\mathbf{g}}_i[n]\cdot E_L[n]\cdot \text{expectationvalue}{\mathbf{g}}_j[n]$$

$$-\text{expectationvalue}{\mathbf{g}}_i[n]\cdot \text{expectationvalue}{\mathbf{h}}_j[n]+\text{expectationvalue}{\mathbf{g}}_i[n]\cdot E_o\cdot \text{expectationvalue}{\mathbf{g}}_j[n].$$

$${\bar{\mathbf{g}}}_i[n]=\{\begin{array}{l}0,\text{ for }i=0\\ {\mathbf{g}}_i[n],\text{ for }i\ne 0\end{array},$$

$${\bar{\mathbf{h}}}_i[n]=\{\begin{array}{l}0,\text{ for }i=0\\ {\mathbf{h}}_i[n],\text{ for }i\ne 0\end{array},$$

$${\bar{\mathbf{G}}}_{\mathbf{r}}{}_i=\{\begin{array}{l}0,\text{ for }i=0\\ {{\mathbf{G}}_{\mathbf{r}}}_i,\text{ for }i\ne 0\end{array},$$

$${\bar{\mathbf{G}}}_{\mathbf{c}}{}_i=\{\begin{array}{l}0,\text{ for }i=0\\ {{\mathbf{G}}_{\mathbf{c}}}_i,\text{ for }i\ne 0\end{array},$$

$$\bar{\mathbf{S}}={\mathbf{e}}_0\cdot {{\mathbf{e}}_0}^T+\text{expectationvalue}\bar{\mathbf{g}}\cdot {\bar{\mathbf{g}}}^T-\text{expectationvalue}\bar{\mathbf{g}}\cdot \text{expectationvalue}{\bar{\mathbf{g}}}^T,$$

$$\bar{\mathbf{H}}=E_0\cdot {\mathbf{e}}_0\cdot {{\mathbf{e}}_0}^T+{\mathbf{e}}_0\cdot {{\bar{\mathbf{G}}}_{\mathbf{r}}}^T+{\bar{\mathbf{G}}}_{\mathbf{c}}\cdot {{\mathbf{e}}_0}^T$$

$$+\text{expectationvalue}\bar{\mathbf{g}}\cdot {\bar{\mathbf{h}}}^T-\text{expectationvalue}\bar{\mathbf{g}}\cdot E_L\cdot \text{expectationvalue}{\bar{\mathbf{g}}}^T$$

$$-\text{expectationvalue}\bar{\mathbf{g}}\cdot \text{expectationvalue}{\bar{\mathbf{h}}}^T+\text{expectationvalue}\bar{\mathbf{g}}\cdot E_0\cdot \text{expectationvalue}{\bar{\mathbf{g}}}^T.$$

$$\bar{\mathbf{S}}\cdot \mathbf{z}={\mathbf{e}}_0\cdot ({{\mathbf{e}}_0}^T\cdot \mathbf{z})+\text{expectationvalue}\bar{\mathbf{g}}\cdot ({\bar{\mathbf{g}}}^T\cdot \mathbf{z})-\text{expectationvalue}\bar{\mathbf{g}}\cdot (\text{expectationvalue}{\bar{\mathbf{g}}}^T\cdot \mathbf{z}),$$

$$\mathbf{A}\cdot \mathbf{v}=\lambda \cdot \mathbf{B}\cdot \mathbf{v}.$$

Full text

I Introduction

Many wavefunction based electronic structure methods rely upon the variational principle, which states that the expectation value of the Hamiltonian for a given wavefunction is an upper bound to the true ground state energy of the Hamiltonian. Thus, searching for an approximate ground state solution amounts to minimizing the energy of an appropriately chosen ansatz. Given such an ansatz $\ket{\Psi(\mathbf{p})}$, we seek the stationary point of the energy functional $E=\frac{\bra{\Psi(\mathbf{p})}\hat{H}\ket{\Psi(\mathbf{p})}}{\bra{\Psi(\mathbf{p})}\ket{\Psi(\mathbf{p})}}$, with respect to the parameters $\{\mathbf{p}_{i}\}$. If the wavefunction ansatz is a linear function of the parameters, the resulting optimization problem becomes relatively easy and many robust algorithms are available. However, in variational Monte Carlo (VMC)Becca and Sorella (2017); Nightingale and Umrigar (1999) the most commonly used wavefunctions depend non-linearly on the parameters. Although this non-linearity leads to great flexibility in the wavefunction, it also introduces strong coupling between the parameters and results in a very challenging non-linear optimization problem. In addition, in Monte Carlo algorithms, all quantities such as energies, gradients and Hessians are only known up to a stochastic error which renders several commonly used algorithms for deterministic optimization ineffective.

Only in the past decade have approaches been developed within VMC that are able to optimize the energy of a wavefunction with more than a few thousand parameters. These techniques include the stochastic reconfiguration (SR) algorithmSorella (2001, 2005); Beaudet et al. (2008) and its low-memory extension introduced by Neuscamman et al.Neuscamman et al. (2012), and adaptive stochastic gradient descent (SGD) algorithms developed in the machine learning community utilized by Booth and coworkersSchwarz et al. (2017) and our groupSabzevari and Sharma (2018). This has led the way to optimization of wavefunction ansatzes with up to $\sim 10^{5}$ parameters. The above algorithms are termed first order optimization methods, as they only rely upon information from first derivatives.

Another very promising approach that has worked extremely well in the face of these challenges is the linear method (LM), which distinguishes itself from the former algorithms by using some information from second derivativesNightingale and Melik-Alaverdian (2001); Umrigar et al. (2007); Toulouse and Umrigar (2007, 2008). This algorithm amounts to solving a generalized eigenvalue problem in the self plus tangent space of the current wavefunction. The LM utilizes a strong zero-variance property, and is the method of choice for fast and robust optimization for relatively small systemsNeuscamman (2016a); Goetz and Neuscamman (2017); Brown et al. (2007); Petruzielo et al. (2012). It has been a challenge to apply the LM to larger systems due to the substantial memory requirement of storing the matrices and the large expense of solving the generalized eigenvalue problem. Neuscamman and coworkers have recently introduced a variant of the LM called the Blocked Linear Method (BLM)Zhao and Neuscamman (2017) that has been applied to wavefunctions containing up to $\sim 25,000$ parameters. In this work we follow an alternative route to extending the LM for large systems. We do so by by employing Davidson’s method with the Jacobi orthogonal component correction to solve the generalized eigenvalue problem iteratively. This yields a solution that only requires repeated action of the Hamiltonian and overlap matrices onto a vector, which can be done at a low scaling cost as discussed by Neuscamman et al.Neuscamman et al. (2012). In addition, by sampling the matrices after a correlation length of moves have been made, one can reduce the space and time complexity of applying these matrices. Further, to tackle truly large systems, we also use AMSGrad, a significantly cheaper optimizer, far from the minimum where the first order optimization methods work very wellReddi et al. (2019) and only switch over to the direct linear method near convergence. These ideas allow us to apply this accelerated LM (aLM) to wavefunctions with up to $\sim 60,000$ parameters, with promising outlook for $10^{5}$ parameters.

The rest of this article is organized as follows. In section II we outline the theory and methods employed, beginning with the VMC framework our group has developed in real and orbital space, and then moving onto the theory of the aLM. Section III goes over the finer computational details of the overall algorithm, and in section IV we show numerical experiments for various systems from chemistry and physics. We end with our final remarks in section V, the conclusion, and VI, acknowledgments.

II Theory and Methods

II.1 Variational Monte Carlo Framework

In the past year, we have developed an improved orbital space variational Monte Carlo (VMC) framework in which the overall cost of the algorithm for obtaining a system size independent stochastic error was lowered from $O(N^{5})$ to $O(N^{4})$Sabzevari and Sharma (2018) (Neuscamman and Wei have achieved the same reduction in cost with a semi-stochastic algorithmWei and Neuscamman (2018)). Here $N$ is a general measure of system size, and is of the same order as the number of orbitals and the number of electrons. The key improvements include:

The use of an effective Hamiltonian that reduces the cost of the local energy calculation from $O(N^{4})$ to $O(N^{2})$; inspired by the Heat-bath Configuration Interaction (HCI) algorithmHolmes et al. (2016). 2. 2.

The use of continuous time Monte Carlo for sampling occupation number vectors that ensures ergodic and efficient movesBortz et al. (1975); Gillespie (1976). 3. 3.

The use of AMSGrad to optimize parametersReddi et al. (2019). This is a stochastic gradient descent algorithm developed in the machine learning community that yields comparable convergence to the popular stochastic reconfiguration (SR) algorithmNeuscamman et al. (2012); Sabzevari and Sharma (2018).

We have recently published an article expanding this framework with highly accurate wavefunction ansatzes consisting of a Jastrow factor and a projected symmetry broken mean-field referenceMahajan and Sharma (2019). These wavefunctions have the form

[TABLE]

where $\ket{\Phi}$ is a reference ansatz with greater variational freedom that breaks certain symmetries, $\hat{P}$ is a symmetry projector that projects out the parts of the reference that satisfy the symmetries of the Hamiltonian, and $\hat{J}$ is a Jastrow operator that adds explicit correlation between electrons. It has been demonstrated that a variational ansatz of the form $\hat{P}\ket{\Phi}$ has more flexibility than an ansatz $\ket{\Phi}$ that obeys the symmetries of the Hamiltonian by constructionBach et al. (1994); Lykos and Pratt (1963); Löwdin (1955). We follow prior work done in the community, and project the correct symmetries before varying the parametersScuseria et al. (2011); Rodríguez-Guzmán et al. (2012); Jiménez-Hoyos et al. (2012). For this paper, the symmerty broken reference will be taken to be a generalized Hartree Fock (GHF) ansatz. GHF wavefunctions are a product of fermionic creation operators acting on the vacuum state,

[TABLE]

with the creation operators defined as linear combinations of local atomic orbital creation operators,

[TABLE]

Note that the sum is over all spin projections, hence the GHF wavefunction breaks both $\hat{S}^{2}$ and $\hat{S}_{z}$ symmetries. By using complex variables for $\theta$, complex conjugation ($\hat{K}$) symmetry can also be broken. We then restore $\hat{K}$ and $\hat{S}_{z}$ symmetries with the aforementioned projector. The orbital space Jastrow factor has the form

[TABLE]

where $\hat{n}_{p\sigma}$ and $\hat{n}_{q\gamma}$ are number operators for the spin orbitals $p\sigma$ and $q\gamma$, respectively, and $J_{p\sigma,q\gamma}$ are variational parameters. A Jastrow factor such as this is part of a larger class of correlators which include correlated product statesMezzacapo et al. (2009); Neuscamman and Chan (2012), and Gutzwiller factorsGutzwiller (1963). This form ensures that the full ansatz is size consistent by using local number projectorsNeuscamman (2012, 2013a) and also imparts some dynamical correlationNeuscamman (2013b, 2016b). The inclusion of a linearized Jastrow-type factor with a spin symmetry projected GHF wavefunction has been used before by Scuseria and HendersonHenderson and Scuseria (2013). In the context of VMC, many studies have utilized an ansatz such as this for model HamiltoniansTahara and Imada (2008a, b); Kurita et al. (2015); Zhao et al. (2017); Darmawan et al. (2018); Misawa et al. (2019).

In addition to the above work, we have recently extended symmetry broken references to real space VMC. The GHF ansatz has the same expression, however the Jastrow factor has a different form in this basis,

[TABLE]

where $r_{ij}$ is the distance between electrons $i$ and $j$, and $r_{iI}$ is the distance between electron $i$ and nucleus $I$. The first summation represents the one body electron-nucleus correlations, the second summation represents the two body electron-electron correlations, and the third summation represents the three body electron-electron-nucleus correlations. The functions $f$, $g$, and $h$ are usually of some polynomial form, here we take them to be the polynomials introduced by Boys and HandyBoys and Handy (1969). This Jastrow factor is chosen to satisfy the electron-electron and electron-nuclear cusp conditions of the true many-body wavefunction, which account for much of the shortcomings of the approximate wavefunctions used throughout quantum chemistryKato (1957); Drummond et al. (2004); Trygve Helgaker (2013). Details of the real space algorithm and wavefunctions will be discussed in a forthcoming publication. For more information on the orbital space algorithm, please refer to our previous publicationsMahajan and Sharma (2019); Sabzevari and Sharma (2018).

II.2 The Linear Method

The variational parameters defining the wavefunction ansatz need to be optimized in accordance with the variational principle; we accomplish this with the accelerated Linear Method (aLM). Here we outline our approach, following previous work done in the communityToulouse and Umrigar (2007, 2008); Umrigar et al. (2007); Nightingale and Melik-Alaverdian (2001); Neuscamman et al. (2012). The LM begins by defining the following $n_{\text{var}}+1$ dimensional basis at any point in parameter space for a given normalized ansatz $\ket{\Psi}$.

[TABLE]

where $\ket{\Psi_{i}}=\frac{\partial}{\partial p_{i}}\ket{\Psi}$. The wavefunction can then be expanded as $\ket{\Psi_{\text{lin}}}=\sum_{i}\ket{\bar{\Psi}_{i}}\delta p_{i}$, with $\delta p_{0}=1$. By minimizing the energy of this linearized wavefunction with respect to the parameters $\{\delta p_{i}\}$, one obtains the following generalized eigenvalue problem:

[TABLE]

where we have defined $\mathbf{\bar{H}}_{ij}=\bra{\bar{\Psi}_{i}}\hat{H}\ket{\bar{\Psi}_{j}}$ and $\mathbf{\bar{S}}_{ij}=\bra{\bar{\Psi}_{i}}\ket{\bar{\Psi}_{j}}$. This eigenvalue problem can be solved to find the ideal parameter updates $\{\delta p_{i}\}$. However, there is an arbitrariness in this procedure, as the normalization constant of $\ket{\Psi}$ is a function of the parameters and how it is defined can lead to unwieldy parameter updates. To take into account this ambiguity, we follow Umrigar and Toulouse and perform the following uniform scaling on the parameter updatesToulouse and Umrigar (2008, 2007):

[TABLE]

where for linear parameters,

[TABLE]

and for non-linear parameters,

[TABLE]

In this work, $\xi=0.5$ for all calculations.

The matrices, $\mathbf{\bar{H}}$ and $\mathbf{\bar{S}}$ can be calculated via Monte Carlo sampling by inserting the resolution of identity.

[TABLE]

[TABLE]

where we have defined a probability distribution $\rho_{n}=\lvert\bra{n}\ket{\Psi}\rvert^{2}$. It is important to note the matrix $\mathbf{\bar{H}}$ will be symmetric when the complete sum is performed, however, this will not be the case when a Monte Carlo average is calculatedNightingale and Melik-Alaverdian (2001). A similar result is obtained with respect to the overlap matrix,

[TABLE]

Inserting explicit expressions for the basis in Eq.(6), it is evident that $\mathbf{\bar{S}}_{00}=1$ and $\mathbf{\bar{H}}_{00}=\expectationvalue{E_{L}[n]}=E_{0}$. Where we recognize a local energy, $E_{L}[n]=\frac{\bra{n}\hat{H}\ket{\Psi}}{\bra{n}\ket{\Psi}}$, which when averaged is the energy of the zeroth order wavefunction $\ket{\Psi}$. It is convenient to introduce the quantities $\mathbf{g}_{i}[n]$, and $\mathbf{h}_{i}[n]$ which are defined as:

[TABLE]

[TABLE]

We can then write the matrix equation as:

[TABLE]

where $\mathbf{G_{r}}$ and $\mathbf{G_{c}}$ are two different estimates of the energy gradient (neglecting a factor of $2$), with dimension $n_{\text{var}}$,

[TABLE]

and

[TABLE]

The newly introduced matrices $\mathbf{H}$ and $\mathbf{S}$ are $n_{\text{var}}\cdot n_{\text{var}}$ dimensional matrices with the form,

[TABLE]

[TABLE]

This rewriting of Eq.(7), together with the introduction of

[TABLE]

[TABLE]

and of

[TABLE]

[TABLE]

allow the matrices $\mathbf{\bar{H}}$ and $\mathbf{\bar{S}}$ to be compactly written as a sum of outer products of vectors,

[TABLE]

[TABLE]

Here $\mathbf{e_{0}}$ is the orthogonal unit vector of index [math] and dimension $n_{\text{var}}+1$. This form suggests that, instead of building these matrices at a cost of $O(n_{\text{s}}(n_{\text{var}})^{2})$ (where $n_{\text{s}}$ is the number of samples), the action of the matrices on a vector $\mathbf{z}$ can be calculated at a cost of $O(n_{\text{s}}n_{\text{var}})$:

[TABLE]

and likewise with the $\mathbf{\bar{H}}$ matrix. This amounts to a series of dot products, scalar multiplications, and addition of vectors. Although this improves the time complexity of the action of the matrices, it introduces a large memory cost by requiring the storage of the primitive values $E_{L}[n]$, $\mathbf{g}_{i}[n]$, and $\mathbf{h}_{i}[n]$ for the entire Monte Carlo simulation over the sampled states $\{n\}$. While this seems excessive, in practice this seemingly large memory cost is still usually smaller than that of the non-direct algorithm in which the matrices $\mathbf{\bar{H}}$ and $\mathbf{\bar{S}}$ are stored in memory. This is because as the VMC algorithm is parallelized across $n_{\text{proc}}$ processors, the memory requirement per processor for the direct algorithm reduces to $O(n_{\text{s}}n_{\text{var}}/n_{\text{proc}})$, while the memory requirement per processor for the non-direct linear method remains $O((n_{\text{var}})^{2})$. As a result, as the degree of parallelization increases, the memory requirement per process for the direct algorithm improves relative to the non direct algorithm. Furthermore, the matrices can be sampled after a correlation length of moves have been made. This leads to a $n_{\text{corr}}$ factor reduction in the space complexity of storing the primitives and the time complexity of applying the matrices onto a vector, yielding a $O(n_{\text{s}}n_{\text{var}}/n_{\text{corr}}/n_{\text{proc}})$ asymptotic CPU and memory cost per process.

II.3 Davidson’s method and the Jacobi orthogonal component correction

In this subsection, we will discuss how to solve the introduced generalized eigenvalue problem with only access to the action of the matrices onto a vector. In this formulation of the LM, we can make use of iterative solvers such as Davidson’s methodDavidson (1975), which allow one to calculate a single or a few close together eigenvalues and eigenvectors without explicit construction of the matricesDavidson (1975); Morgan and Scott (1986); Morgan (1992). We will only focus on calculating a single eigenpair $(\lambda,\mathbf{v})$ for the generalized eigenproblem,

[TABLE]

$\mathbf{A}$ is a general square matrix and $\mathbf{B}$ is the overlap matrix or metric, which is positive semi-definite. Assuming we have a reasonable guess $\mathbf{x}$ to the eigenvector $\mathbf{v}$, we seek a correction $\mathbf{\delta}$,

[TABLE]

where $\lambda$ is replaced by the Rayleigh quotient of our guess $\rho=\frac{\mathbf{x}^{\dagger}\mathbf{A}\mathbf{x}}{\mathbf{x}^{\dagger}\mathbf{B}\mathbf{x}}$. Rearranging the above equation yields,

[TABLE]

where $\mathbf{r}$ is the so called residual vector equal to $(\mathbf{A}-\rho\cdot\mathbf{B})\cdot\mathbf{x}$. This can conceivably be solved for the correction, but in reality gives a poor update due to the fact $\mathbf{x}$ and $\mathbf{\delta}$ can be linearly dependent. In practice, an approximate correction equation is solved,

[TABLE]

where $\mathbf{F}$ is a preconditioner, ie. some matrix similar to $(\mathbf{A}-\rho\cdot\mathbf{B})$. Davidson’s original paper was concerned with problems where $\mathbf{B}=\mathbf{I}$, and by assuming $\mathbf{A}$ is well approximated by it’s diagonal values, he came to the following expression for the correction,

[TABLE]

where $\mathbf{A_{diag}}$ is the matrix containing only the diagonal of $\mathbf{A}$. Eq.(31) can be recognized as a perturbative correction to $\mathbf{x}$ and is very easy to calculate. For diagonally dominant matrices like many of those in quantum chemistry, this correction is remarkably successfulDavidson (1975); Van Dam et al. (1996); Trygve Helgaker (2013). The $\mathbf{\bar{H}}$ matrix is not diagonally dominant and the $\mathbf{\bar{S}}$ matrix is not equal to the identity, due to the aforementioned strong coupling of the wavefunction parameters. This forces us to find other expressions for the matrix $\mathbf{F}$.

The Jacobi orthogonal component correction (JOCC) attempts to find a correction $\mathbf{\delta}$, that is orthogonal to the current best guess $\mathbf{x}$. It does this by solving a projected form of Eq.(29) with

[TABLE]

Where $\mathbf{P}=(\mathbf{I}-\mathbf{x}\cdot\mathbf{y}^{\dagger})$ and $\mathbf{y}=\mathbf{B}\cdot\mathbf{x}$Sleijpen et al. (1996). This $\mathbf{F}$ matrix happens to be positive semi-definite (with the addition of a diagonal shift to $\mathbf{A}$) and only requires a handful of iterations of conjugate gradient to generate a reasonable correction. Alternatively, one can also set a loose convergence threshold for $\mathbf{\delta}$ compared to the tolerance required of $\mathbf{x}$. Davidson’s method with the JOCC is termed the Jacobi-Davidson method, and generalizes Davidson’s method to any matrix, no matter the magnitude of the off-diagonal elementsSleijpen et al. (1996); Sleijpen and Van der Vorst (2000); Bai et al. (2000).

Here we’ll outline the algorithmic details of Davidson’s method. The goal is to solve the eigenvalue problem in a gradually growing subspace of orthonormal vectors with the hope that the solution exists in the span of this subspace. A general Davidson’s algorithm with input matrices $\mathbf{A}$ and $\mathbf{B}$, a target value $t$, and a reasonable guess vector $\mathbf{x_{0}}$, has the following series of steps.

Initialize the subspace matrices $\mathbf{V}=[\mathbf{x_{0}}]$, $\mathbf{V_{A}}=[\mathbf{A}\mathbf{x_{0}}]$, and $\mathbf{V_{B}}=[\mathbf{B}\mathbf{x_{0}}]$. 2. 2.

Solve the subspace problem,

[TABLE]

where $\mathbf{A^{\prime}}=\mathbf{V}^{\dagger}\mathbf{A}\mathbf{V}=\mathbf{V}^{\dagger}\mathbf{V_{A}}$ and $\mathbf{B^{\prime}}=\mathbf{V}^{\dagger}\mathbf{B}\mathbf{V}=\mathbf{V}^{\dagger}\mathbf{V_{B}}$. Select the eigenpair $(\theta,\mathbf{s})$ in which $\theta$ is closest to $t$. 3. 3.

Transform $\mathbf{s}$ into the original problem space, $\mathbf{u}=\mathbf{V}\cdot\mathbf{s}$, and calculate action of matrices, $\mathbf{u_{A}}=\mathbf{A}\cdot\mathbf{u}=\mathbf{V_{A}}\cdot\mathbf{s}$ and $\mathbf{u_{B}}=\mathbf{B}\cdot\mathbf{u}=\mathbf{V_{B}}\cdot\mathbf{s}$ 4. 4.

Calculate the residual vector $\mathbf{r}=(\mathbf{A}-\theta\cdot\mathbf{B})\cdot\mathbf{u}=\mathbf{u_{A}}-\theta\cdot\mathbf{u_{B}}$. If $\norm{\mathbf{r}}$ is less than some predefined threshold, return $(\theta,\mathbf{u})$. 5. 5.

Approximately solve the correction equation for $\mathbf{\delta}$,

[TABLE]

with a chosen preconditioner $\mathbf{F}$. In the Jacobi-Davidson scheme, $\mathbf{F}=(\mathbf{I}-\mathbf{u_{B}}\cdot\mathbf{u}^{\dagger})(\mathbf{A}-\theta\cdot\mathbf{B})(\mathbf{I}-\mathbf{u}\cdot\mathbf{u_{B}}^{\dagger})$ 6. 6.

Orthonormalize $\mathbf{\delta}$ with respect to the subspace, the columns of $\mathbf{V}$. This can be done with a Gram-Schmidt procedure:

1. (a)

for $\mathbf{v}\in\text{Cols}(\mathbf{V})$

$\delta:=\delta-\mathbf{v}\cdot(\mathbf{v}^{\dagger}\cdot\mathbf{B}\cdot\mathbf{\delta})$ 2. (b)

$\mathbf{x_{t}}=\mathbf{\delta}/\norm{\mathbf{\delta}}$, where $\norm{\mathbf{\delta}}=\sqrt{\mathbf{\delta}^{\dagger}\cdot\mathbf{B}\cdot\mathbf{\delta}}$ 7. 7.

Append $\mathbf{x_{t}}$ to the subspace. $\mathbf{V}:=[\mathbf{V},\mathbf{x_{t}}]$, $\mathbf{V_{A}}:=[\mathbf{V_{A}},\mathbf{A}\mathbf{x_{t}}]$, and $\mathbf{V_{B}}:=[\mathbf{V_{B}},\mathbf{B}\mathbf{x_{t}}]$ 8. 8.

Repeat steps 2-8 until a desired convergence for $\norm{\mathbf{r}}$ has been achieved.

To ensure the size of the subspace doesn’t grow too large, restarts should be taken advantage of. This can be accomplished very easily by taking a number of the eigenvectors in the subspace problem that are closest to the desired target value as the new subspace. In detail, if $|\text{Cols}(\mathbf{V})|$ is greater than or equal to some $\mathbf{V_{max}}$, define $\mathbf{N}$ as a matrix whose columns are a $\mathbf{V_{restart}}$ number of eigenvectors of the subspace problem (ie. the eigenproblem with matrices $\mathbf{A^{\prime}}$, $\mathbf{B^{\prime}}$) with eigenvalues closest to $t$. The new subspace matrices can then be calculated, $\mathbf{V}:=\mathbf{V}\cdot\mathbf{N}$, $\mathbf{V_{A}}:=\mathbf{V_{A}}\cdot\mathbf{N}$, and $\mathbf{V_{B}}:=\mathbf{V_{B}}\cdot\mathbf{N}$. We find $\mathbf{V_{max}}=25$ and $\mathbf{V_{restart}}=5$ work reasonably well for matrices of any dimension. It is also useful to update the target value $t$, when the residual norm $\norm{\mathbf{r}}$, is smaller than some threshold, and to continue updating it to $\theta$ if the corresponding eigenvector has a smaller residual norm.

III Computational Details

Given a wavefunction, a VMC run utilizing the LM optimizer can be performed efficiently if the primitives $E_{L}[n]$, $\mathbf{g}_{i}[n]$, and $\mathbf{h}_{i}[n]$ can be computed with a low computational cost.

For the orbital space algorithm, previous publications have demonstrated that the local energy and $\mathbf{g}_{i}[n]$ can be computed at a cost of $O(N^{2})$ if intermediate values are pre-calculatedSabzevari and Sharma (2018); Mahajan and Sharma (2019). The third term $\mathbf{h}_{i}[n]$, can be naively calculated at a cost of $O(N^{4})$, but by considering the gradient of the local energy with respect to the $i\text{th}$ parameter,

[TABLE]

one can calculate the last primitive using,

[TABLE]

Noting that the local energy can be calculated directly from the Jastrow and Hartree-Fock parameters at a cost of $O(N^{3})$, a reverse accumulation algorithmic differentiation scheme will calculate the gradient at the same asymptotic cost, albeit with a fairly large pre-factorLinnainmaa (1976); Griewank and Walther (2008). To do this we turn to the Stan Math library, an excellent C++ implementation of algorithmic differentiation, that was very easily extendable to linear algebra with complex numbersCarpenter et al. (2015).

For the real space algorithm, we comment on the computational costs of calculating each primitive, the details of which will be a part of a forthcoming publication. The local energy can be calculated at a cost of $O(N^{2})$, while $\mathbf{g}_{i}[n]$ and $\mathbf{h}_{i}[n]$ are calculated at a cost of $O(N^{3})$. The expressions for these values are very similar to those found in reference Filippi et al. (2016), with the added complication that the reference is composed of complex variables and is projected onto the correct symmetry before any derivatives are taken.

For both the real space and orbital space wavefunctions, the VMC algorithm for sampling the energy and the matrices in the direct linear method scales as $O(N^{4})$ when a system size independent stochastic error is required. This scaling comes about by using the fact that the primitives are serially correlated, with a correlation length that can be assumed to increase linearly with the size of the system. Thus by only sampling the expensive primitives every $O(N)$ Monte Carlo moves, and having to perform a total of $O(N^{2})$ moves to obtain a system-size independent error, yields an overall scaling of $O(N^{4})$.

To stabilize the linear method we add a diagonal shift to the $\mathbf{\bar{H}}$ matrix, which acts as a trust radius to the eigenvector solution. This diagonal shift decays by a multiplicative factor every iteration of the linear method until it reaches a value of $10^{-6}$. We find that an initial shift of $0.1$ and a decay factor of $0.65$ yield excellent results, and will be used for the majority of the calculations in this paper. Furthermore, to guarantee a large step in a poor update direction is never taken, we employ a short correlated sampling run with wavefunctions created from various step sizes multiplied by the LM update. The step size that generated the lowest relative energy is taken as the wavefunction to begin the next VMC iteration. In detail, we calculate new wavefunctions with the original parameters plus $0.01,0.05,0.1,0.5,1.0$ times the LM update, and perform a VMC run to calculate the energy using the median value step size, and then estimate the energies for the other wavefunctions with,

[TABLE]

where $\ket{\Psi}$ is the wavefunction that the random walk is performed with, and $\ket{\Psi^{\prime}}$ is a “slave” wavefunction that should not be too different from $\ket{\Psi}$. This expression minimizes the relative error between the averaged energies, and allows a short VMC run to be used to calculate energies that can then be compared between each other with high accuracy. We use $0.35$ times $n_{\text{s}}$ number of iterations to perform correlated sampling.

When optimizing a large system with the LM, the linear expansion of the wavefunction can be a poor representation of the working non-linear ansatz, particularly far from the global minumum. When this occurs, none of the strategies discussed above will result in decent convergence. To remedy this, we turn to the first order SGD method, AMSGrad, at the beginning of the optimization. This optimizer, developed in the machine learning community, has given us comparable convergence to SR in every practical example our group has attempted, in addition to requiring significantly less overhead. Although the algorithm has been discussed in a previous publication, here we briefly restate the AMSGrad method, since it is a critical part of the presented algorithm and remedies a known problem with the LMBecca and Sorella (2017). A general algorithm for adaptive SGD methods is the following,

[TABLE]

[TABLE]

[TABLE]

where $\alpha$ is the stepsize and $f$ and $g$ are some functions that take in all the past gradients and generate vectors $\mathbf{m}^{(i)}$ and $\mathbf{n}^{(i)}$. These are then used to calculate the updates $\{\delta p_{i}\}$. Adaptive SGD methods differ in the chosen functions, $f$ and $g$. In AMSGrad, $f$ and $g$ respectively calculate the exponentially decaying moving average of the first and second moment of the gradients $\{\mathbf{G}^{(i)}\}$,

[TABLE]

[TABLE]

with the caveat that the second moment at iteration $(i)$ is always greater than at iteration $(i-1)$, i.e. $\mathbf{n}^{(i)}_{j}\geq\mathbf{n}^{(i-1)}_{j}$. This ensures that the learning rate, $\alpha/\sqrt{\mathbf{n}^{(i)}_{j}}$, is a monotonically decreasing function of the number of iterations. This was shown to improve the convergence of AMSGrad for a synthetic problem over other adaptive SGD methodsReddi et al. (2019). When optimizing with AMSGrad, unless otherwise specified, we have used the aggressive parameters, $\alpha=0.01,\beta_{1}=0.1,\beta_{2}=0.01$. When starting off an aLM optimization with AMSGrad we use the more stable parameters, $\alpha=0.001,\beta_{1}=0.1,\beta_{2}=0.001$.

IV Results and Discussion

We begin by presenting results on a couple of atomic systems with real space Jastrow wavefunctions, and then discuss results for various realistic and model molecular systems with orbital space Jastrow wavefunctions. The focus of the analysis will be on comparing aLM to LM and the first order optimization methods, AMSGrad and SR. Integrals and initial GHF parameters for all calculations were generated by PyscfSun et al. (2018).

IV.1 Beryllium Atom - Real Space

Table 1 shows VMC results for the Beryllium atom using various optimizers. To investigate the direct LM implementation, no AMSGrad iterations were used in the aLM optimization. This is a small system consisting of 4 electrons in 17 orbitals yielding 364 total variables, yet it is already evident that the aLM outperforms the original LM, as well as the cheaper but less robust SGD method, AMSGrad.

IV.2 Neon Atom - Real Space

Figure 1 and Table 2 illustrate optimization of the Neon atom. This is a system of 10 electrons in 10 orbitals giving a wavefunction that consists of 492 total variables. We compare aLM curves with varying tolerances on the Davidson and conjugate gradient solvers, to the usual LM. To investigate the direct LM implementation, no AMSGrad iterations were used in the aLM optimization.

Inspecting Figure 1, we see very little difference in the optimization curves between the different implementations of the linear method, however, the difference in timings and energies in Table 2 is more pronounced. This suggests that a schedule of progressively tighter tolerances would ensure fast and robust optimization.

IV.3 Carbon Dimer - Orbital Space

Here we pivot to orbital space VMC, beginning with the carbon dimer at equilibrium geometry. The bond length is taken to be the experimentally determined value of 1.24244 angstromsXu and Dunning Jr (2013); Douay et al. (1988). This system consists of 12 electrons in 18 orbitals and the resulting wavefunction has a total of 1,530 parameters. Note again, no AMSGrad iterations were used in the aLM optimization to focus on the direct LM implementation.

Figure 2 shows the optimized energy of the Carbon dimer with the number of VMC iterations, while Figure 3 shows the same optimization but with the wall time plotted on the x-axis. It’s quite clear from inspecting the two plots that both first order methods (SR and AMSGrad) behave very similarly whether one considers the number of iterations or the total CPU time. The convergence of the various flavors of aLM with different tolerances and LM as a function of the number of iterations is very similar. However, the CPU time for the LM is significantly greater than that of the aLM method and even the first iteration of the LM method takes $1000+$ seconds and does not appear in the plot of Figure 3.

Table 3 shows that all tolerances optimized to the same energy. It should be noted that tighter tolerance does ensure convergence to the final answer in a smaller number of iterations, although this seems to be balanced out by the fact that each iteration of the looser tolerances costs less. A schedule of progressively tightening tolerances again seems like the most efficient course of action. In Figure 2, we also see that the first order methods optimize very well in the first $\sim 20$ iterations, but are hindered by long tails, needing $\sim 300$ more iterations to converge to the final value. This motivates the use of a first order method in the initial optimization of large systems, as the time per iteration of these methods is at least an order of magnitude cheaper than any of the variants of LM. We also remark on the fact that SR and AMSGrad give very comparable optimization curves for this system, again noting the significantly smaller memory overhead required by AMSGrad. AMSGrad amounts to a single loop over the parameters of the wavefunction, while SR requires storing the primitives for the overlap matrix, and then solving a system of linear equations. While this extra asymptotic cost for SR is comparable to the proposed scheme for aLM, aLM’s unmatched optimization make this cost worth the effort.

IV.4 $\mathbf{H}_{50}$ - Orbital Space

Table 4 illustrates the optimization of a chain of 50 Hydrogen atoms at a bond length of 2.0 Bohr (1.058354 angstroms) using various optimizers. This is a system of 50 electrons in 50 orbitals and results in an ansatz with 15,050 parameters.

The lowest energy for each optimizer is the converged value for that optimizer. SR and AMSGrad optimize to $-26.886$ Ha in a similar number of iterations, and have very comparable timings per iteration. However, AMSGrad continues to optimize to an energy more than 20 mHa lower than SR, but it’s quite clear there is a very long tail as it takes more than 700 more iterations to lower the energy by another 3 mHa. As expected, aLM optimizes significantly faster reaching the SR energy in 36 iterations (25 AMSGrad with the less aggressive parameters plus 16 iterations with aLM and the loose tolerances), then reaching the AMSGrad energy in 33 more iterations, and reaching an energy lower than either of the other optimizers in 89 total iterations. Remarkably, when comparing total time to reach a common converged energy, the aLM is cheaper than SR and AMSGrad. We expect this to be further improved if hand-coded analytic expressions for all the primitives can be implemented as the algorithmic differentiation scheme carries a large computational prefactor.

IV.5 Hydrogen chains - Orbital Space

In this section, we investigate practical scaling of the aLM algorithm, considering only the time spent to solve the generalized eigenvalue problem. Calculations were performed on linear hydrogen chains consisting of 10 to 60 atoms at a bond length of 2.0 Bohr in a STO-6G basis.

The computational experiments suggest that the solver has an asymptotic cost that is slightly higher than the application of the matrix on a vector, a very satisfactory result. Inserting $n_{\text{var}}=O(N^{2})$, we attain an $O(n_{\text{s}}N^{2.32}/n_{\text{proc}})$ cost as a function of system size.

IV.6 98-site Hubbard model - Orbital Space

Here we discuss the 98 site Hubbard model with $U/t=2$. A system of 98 electrons in 98 orbitals, this ansatz results in 57,722 parameters. We use this example to investigate any bias in the aLM optimizer. Due to the fact the matrix elements of $\mathbf{\bar{H}}$ and $\mathbf{\bar{S}}$ have stochastic noise, any eigenvector will inherit a non-linear bias, which is expected to increase with the size of the matrix. Checkpoints in the optimization can be seen in Table 5. The converged energy per electron using the aLM optimizer is -1.19208 $t$ and Green’s function Monte Carlo (GFMC) reference energy is -1.19622 $t$LeBlanc et al. (2015). Error on the energy per electron values are 0.00001 $t$.

With the converged wavefunction, we perform 500 iterations using plain SGD with a small step size of 0.0001. SGD is an unbiased optimizer and if the aLM results in any bias it should be evident in the distribution of energies. In Table 5, we see that the distribution of energies result in a Gaussian, and can be presumed to be from noise in the gradient as the wavefunction is bounced around the minimum. The converged aLM energy falls directly in the middle bin of the histogram.

IV.7 Iron(II) Porphyrin - Orbital Space

We apply the orbital space J-$KS_{z}$GHF wavefunction and the aLM optimizer to calculate the energy of the Fe(II) porphyrin (Fe(P)) cluster in the 32 electrons, 29 orbitals active space of Li ManniLi Manni et al. (2016). The main aim of the current calculation is to investigate whether VMC calculates the same energy gap between the ${}^{5}\text{A}_{1g}$ and ${}^{3}\text{B}_{1g}$ states as the near exact Semi-stochastic Heat-bath Configuration Interaction (SHCI) algorithmSharma et al. (2017); Holmes et al. (2016); Smith et al. (2017).

Table 6 presents results on the Fe(P) system, for various ab-initio methods. The predicted energies of the VMC calculation are quite good, the small discrepancies can be attributed to stochastic noise, and lack of dynamical correlation in the ansatz. GFMCRunge (1992); Trivedi and Ceperley (1990); Sorella (1998); Van Bemmel et al. (1994) on top of the optimized J-$KS_{z}$GHF wavefunction, give results equivalent to the SHCI energies. The gap predicted by either of the Monte Carlo calculations agrees with the SHCI calculations up to 1 mHa. It is quite remarkable that a wavefunction with merely 5,423 parameters gives results comparable to a CI wavefunction with more than 5 million parameters, a true testament to the ability of these non-linear ansatzes to capture electron correlation. We want to note however, the Jastrow factors of the ansatz are heavily dependent on the orbitals that they are defined in. We localized active space orbitals from an HCI-SCF calculation, and differing localization schemes would result in differing energies, and it is not clear how to choose a basis to achieve the correct answer irrespective of the system, the orbitals, and the localization method. Alternatively, Neuscamman et al. have demonstrated that one can optimize the local orbitals simultaneously with the other variational parameters, and further improve the effectiveness of the Jastrow factorNeuscamman (2016a). We also note that the correlation length for sampling the optimized wavefunctions is $\sim 1.5$; almost every move yields an uncorrelated sample. We believe this is due to the efficient moves generated by the continuous time Monte Carlo algorithm.

V Conclusion

In this article, we have detailed an efficient algorithm of the linear method to optimize non-linear wavefunction ansatzes that exhibit strong coupling between their parameters. We have demonstrated that the optimizer functions equally well for wavefunctions with real or orbital space Jastrow factors and results in an order of magnitude reduction in wall time compared to the naive implementation.

For wavefunctions with real space Jastrow factors, our work here indicates that the aLM is the optimizer of choice, as the first order methods struggle to optimize consistently. During the writing of this article we became aware of a related study by the group of Neuscamman et al. comparing various optimizersOtis and Neuscamman (2019). They suggest that the most efficient convergence is obtained when the bulk of the iterations are carried out using an adaptive stochastic gradient descent interspersed with the linear method (in their work they use the BLM). This conclusion is not in disagreement with ours here, as both studies argue that exclusively using a gradient descent algorithm often results in unmanageably slow and poor convergence.

For wavefunctions with orbital space Jastrow factors, the aLM optimizer brings the computer time required to use the linear method to the same order of magnitude as that of the first order methods: AMSGrad and SR. Our results do suggest that for challenging optimization problems the linear method will more consistently optimize to a lower energy. However, the low cost per iteration of the SGD methods, like AMSGrad, suggests that using AMSGrad away from the minimum and switching over to the aLM closer to the minimum will lead to convergence with the least computational cost. We do expect a further reduction in cost of the aLM if an analytic expression for the local energy gradient can be implemented, as the majority of the computational time in the orbital space algorithm is spent performing algorithmic differentiation, which scales well but comes with the largest prefactor in the algorithm. Future work will focuses on implementing hand-coded analytic expressions for the gradient of the local energy, to bring this prefactor down making the aLM algorithm more computationally favorable relative to the first order methods.

The good agreement between VMC results and the near exact SHCI calculations for the calculated triplet-quintet splitting of the Fe(P) molecule illustrates that the J-$KS_{z}$GHF ansatz is flexible enough to be effectively used as an active space solver for challenging strongly correlated chemical systems. The polynomial scaling of the VMC algorithm opens up the possibility of using the algorithm to tackle all valence active spaces that can’t be treated with methods such as full configuration interaction quantum Monte Carlo, selected configuration interaction, or density matrix renormalization group, all of which scale exponentially with the size of the active space.

VI Acknowledgements

The funding for this project was provided by the National Science Foundation through the grant CHE-1800584. IS was also partially supported by the GAANN program through the DOEd.