Solving the Quantum Many-Body Problem with Artificial Neural Networks

Neural-Network Quantum States —

Consider a quantum system with $N$ discrete-valued degrees of freedom $\mathcal{S}=(\mathcal{S}_{1},\mathcal{S}_{2}\dots\mathcal{S}_{N})$, which may be spins, bosonic occupation numbers, or similar. The many-body wave function is a mapping of the $N-$dimensional set $\mathcal{S}$ to (exponentially many) complex numbers which fully specify the amplitude and the phase of the quantum state. The point of view we take here is to interpret the wave function as a computational black box which, given an input many-body configuration $\mathcal{S}$, returns a phase and an amplitude according to $\Psi(\mathcal{S})$. Our goal is to approximate this computational black box with a neural network, trained to best represent $\Psi(\mathcal{S})$. Different possible choices for the artificial neural-network architectures have been proposed to solve specific tasks, and the best architecture to describe a many-body quantum system may vary from one case to another. For the sake of concreteness, in the following we specialize our discussion to restricted Boltzmann machines (RBM) architectures, and apply them to describe spin $1/2$ quantum systems. In this case, RBM artificial networks are constituted by one visible layer of $N$ nodes, corresponding to the physical spin variables in a chosen basis (say for example $\mathcal{S}=\sigma_{1}^{z},\dots\sigma_{N}^{z}$) , and a single hidden layer of $M$ auxiliary spin variables ($h_{1}\dots h_{M}$) (see Fig. 1). This description corresponds to a variational expression for the quantum states which reads:

where $h_{i}=\{-1,1\}$ is a set of $M$ hidden spin variables, and the weights $\mathcal{W}=\{a_{i},b_{j},W_{ij}\}$ fully specify the response of the network to a given input state $\mathcal{S}$. Since this architecture features no intra-layer interactions, the hidden variables can be explicitly traced out, and the wave function reads $\Psi(\mathcal{S};\mathcal{W})=e^{\sum_{i}a_{i}\sigma_{i}^{z}}\times\Pi_{i=1}^{M}F_{i}(\mathcal{S}),$ where $F_{i}(\mathcal{S})=2\cosh\left[b_{i}+\sum_{j}W_{ij}\sigma_{j}^{z}\right]$. The network weights are, in general, to be taken complex-valued in order to provide a complete description of both the amplitude and the wave-function’s phase.

The mathematical foundations for the ability of NQS to describe intricate many-body wave functions are the numerously established representability theorems kolmogorov1961onthe ; hornik1991approximation ; leroux2008representational , which guarantee the existence of network approximates of high-dimensional functions, provided a sufficient level of smoothness and regularity is met in the function to be approximated. Since in most physically relevant situations the many-body wave function reasonably satisfies these requirements, we can expect the NQS form to be of broad applicability. One of the practical advantages of this representation is that its quality can, in principle, be systematically improved upon increasing the number of hidden variables. The number $M$ (or equivalently the density $\alpha=M/N$) then plays a role analogous to the bond dimension for the MPS. Notice however that the correlations induced by the hidden units are intrinsically non local in space and are therefore well suited to describe quantum systems in arbitrary dimension. Another convenient point of the NQS representation is that it can be formulated in a symmetry-conserving fashion. For example, lattice translation symmetry can be used to reduce the number of variational parameters of the NQS ansatz, in the same spirit of shift-invariant RBM’s sohn2012learning ; norouzi2009stacksof . Specifically, for integer hidden variable density $\alpha=1,2,\dots$, the weight matrix takes the form of feature filters $W_{j}^{(f)}$ , for $f\in[1,\alpha]$. These filters have a total of $\alpha N$ variational elements in lieu of the $\alpha N^{2}$ elements of the asymmetric case (see Supp. Mat. for further details).

Given a general expression for the quantum many-body state, we are now left with the task of solving the many-body problem upon machine learning of the network parameters $\mathcal{W}$. In the most interesting applications the exact many-body state is unknown, and it is typically found upon solution either of the static Schrödinger equation $\mathcal{H}\left|\Psi\right\rangle=E\left|\Psi\right\rangle$, either of the time-dependent one $i\mathcal{H}\left|\Psi\right(t)\rangle=\frac{d}{dt}\left|\Psi(t)\right\rangle$, for a given Hamiltonian $\mathcal{H}$. In the absence of samples drawn according to the exact wave function, supervised learning of $\Psi$ is therefore not a viable option. Instead, in the following we derive a consistent reinforcement learning approach, in which either the ground-state wave function or the time-dependent one are learned on the basis of feedback from variational principles.

Ground State —

To demonstrate the accuracy of the NQS in the description of complex many-body quantum states, we first focus on the goal of finding the best neural-network representation of the unknown ground state of a given Hamiltonian $\mathcal{H}$. In this context, reinforcement learning is realized through minimization of the expectation value of the energy $E(\mathcal{W})=\langle\Psi_{M}|\mathcal{H}|\Psi_{M}\rangle/\langle\Psi_{M}|\Psi_{M}\rangle$ with respect to the network weights $\mathcal{W}$. In the stochastic setting, this is achieved with an iterative scheme. At each iteration $k$, a Monte Carlo sampling of $\left|\Psi_{M}(S;\mathcal{W}_{k})\right|^{2}$ is realized, for a given set of parameters $\mathcal{W}_{k}$. At the same time, stochastic estimates of the energy gradient are obtained. These are then used to propose a next set of weights $\mathcal{W}_{k+1}$ with an improved gradient-descent optimization sorella2007weakbinding . The overall computational cost of this approach is comparable to that of standard ground-state Quantum Monte Carlo simulations (see Supp. Material).

To validate our scheme, we consider the problem of finding the ground state of two prototypical spin models, the transverse-field Ising (TFI) model and the anti-ferromagnetic Heisenberg (AFH) model. Their Hamiltonians are

respectively, where $\sigma^{x},\sigma^{y},\sigma^{z}$ are Pauli matrices.

In the following, we consider the case of both one and two dimensional lattices with periodic boundary conditions (PBC). In Fig. 2 we show the optimal network structure of the ground states of the two spin models for a hidden variables density $\alpha=4$ and with imposed translational symmetries. We find that each filter $f=[1,\dots\alpha]$ learns specific correlation features emerging in the ground state wave function. For example, in the 2D case it can be seen (Fig. 2, rightmost panels) how the neural network learns patterns corresponding to anti-ferromagnetic correlations. The general behavior of the NQS is completely analogous to what observed in convolutional neural networks, where different layers learn specific structures of the input data.

In Fig. 3 we show the accuracy of the NQS states, quantified by the relative error on the ground-state energy $\epsilon_{\textrm{rel}}=\left(E_{\textrm{NQS}}(\alpha)-E_{\textrm{exact}}\right)/\left|E_{\textrm{exact }}\right|$, for several values of $\alpha$ and model parameters. In the left panel, we compare the variational NQS energies with the exact result obtained by fermionization of the TFI model, on a one-dimensional chain with PBC. The most striking result is that NQS achieve a controllable and arbitrary accuracy which is compatible with a power-law behavior in $\alpha$. The hardest to learn ground-state is at the quantum critical point $h=1$, where nonetheless a remarkable accuracy of one part per million can be easily achieved with a relatively modest density of hidden units. The same remarkable accuracy is obtained for the more complex one-dimensional AFH model (central panel). In this case we observe as well a systematic drop in the ground-state energy error, which for a small $\alpha=4$ attains the same very high precision obtained for the TFI model at the critical point. Our results are compared with the accuracy obtained with the spin-Jastrow ansatz (dashed line in the central panel), which we improve by several orders of magnitude. It is also interesting to compare the value of $\alpha$ with the MPS bond dimension $M$, needed to reach the same level of accuracy. For example, on the AFH model with PBC, we find that with a standard DMRG implementation dolfi2014matrixproduct we need $M\sim 160$ to reach the accuracy we have at $\alpha=4$. This points towards a more compact representation of the many-body state in the NQS case, which features about $3$ orders of magnitude less variational parameters than the corresponding MPS ansatz.

We next study the AFH model on a two-dimensional square lattice, comparing in the right panel of Fig. 3 to QMC results sandvik1997finitesize . As expected from entanglement considerations, the 2D case proves harder for the NQS. Nonetheless, we always find a systematic improvement of the variational energy upon increasing $\alpha$, qualitatively similar to the 1D case. The increased difficulty of the problem is reflected in a slower convergence. We still obtain results at the level of existing state-of-the-art methods or better. In particular, with a relatively small hidden unit density $(\alpha\sim 4)$ we already obtain results at the same level than the best known variational ansatz to-date for finite clusters (the EPS of Ref. mezzacapo2009groundstate and the PEPS states of Ref. lubasch2014algorithms ). Further increasing $\alpha$ then leads to a sizable improvement and consequently yields the best variational results so-far-reported for this 2D model on finite lattices.

Unitary Dynamics —

NQS are not limited to ground-state problems but can be extended to the time-dependent Schrödinger equation. For this purpose we define complex-valued and time-dependent network weights $\mathcal{W}(t)$ which at each time $t$ are trained to best reproduce the quantum dynamics, in the sense of the Dirac-Frenkel time-dependent variational principle dirac1930noteon ; frenkel1934wavemechanics . In this context, the variational residuals

are the objective functions to be minimized as a function of the time derivatives of the weights $\dot{\mathcal{W}}(t)$ (see Supp. Mat.) In the stochastic framework, this is achieved by a time-dependent VMC method carleo2012localization ; carleo2014lightcone , which samples $\left|\Psi_{M}(S;\mathcal{W}(t))\right|^{2}$ at each time and provides the best stochastic estimate of the $\dot{\mathcal{W}}(t)$ that minimize $R^{2}(t)$, with a computational cost $\mathcal{O}(\alpha N^{2})$. Once the time derivatives determined, these can be conveniently used to obtain the full time evolution after time-integration.

To demonstrate the effectiveness of the NQS in the dynamical context, we consider the unitary dynamics induced by quantum quenches in the coupling constants of our spin models. In the TFI model we induce a non-trivial quantum dynamics by means of an instantaneous change in the transverse field: the system is initially prepared in the ground-state of the TFI model for some transverse field, $h_{i}$, and then let evolve under the action of the TFI Hamiltonian with a transverse field $h_{f}\neq h_{i}$. We compare our results with the analytical solution obtained from fermionization of the TFI model for a one-dimensional chain with PBC. In the left panel of Fig. 4 the exact results for the time-dependent transverse spin polarization are compared to NQS with $\alpha=4$. In the AFH model, we study instead quantum quenches in the longitudinal coupling $J_{z}$ and monitor the time evolution of the nearest-neighbors correlations. Our results for the time evolution (and with $\alpha=4$ ) are compared with the numerically-exact MPS dynamics white2004realtime ; vidal2004efficient ; daley2004timedependent for a system with open boundaries (see Fig. 4, right panel).

The high accuracy obtained also for the unitary dynamics further confirms that neural network-based approaches can be fruitfully used to solve the quantum many-body problem not only for ground-state properties but also to model the evolution induced by a complex set of excited quantum states. It is all in all remarkable that a purely stochastic approach can solve with arbitrary degree of accuracy a class of problems which have been traditionally inaccessible to QMC methods for the past $50$ years. The flexibility of the NQS representation indeed allows for an effective solution of the infamous phase problem plaguing the totality of existing exact stochastic schemes based on Feynman’s path integrals.

Outlook —

Variational quantum states based on artificial neural networks can be used to efficiently capture the complexity of entangled many-body systems both in one a two dimensions. Despite the simplicity of the restricted Boltzmann machines used here, very accurate results for both ground-state and dynamical properties of prototypical spin models can be readily obtained. Potentially many novel research lines can be envisaged in the near future. For example, the inclusion of the most recent advances in machine learning, like deep network architectures, might be further beneficial to increase the expressive power of the NQS. Furthermore, the extension of our approach to treat quantum systems other than interacting spins is, in principle, straightforward. In this respect, applications to answer the most challenging questions concerning interacting fermions in two-dimensions can already be anticipated. Finally, at variance with Tensor Network States, the NQS feature intrinsically non-local correlations which can lead to substantially more compact representations of many-body quantum states. A formal analysis of the NQS entanglement properties might therefore bring about substantially new concepts in quantum information theory.

References

Appendix A Stochastic Optimization For The Ground State

In the first part of our Paper we have considered the goal of finding the best representation of the ground state of a given quantum Hamiltonian $\mathcal{H}$. The expectation value over our variational states $E(\mathcal{W})=\langle\Psi_{M}|\mathcal{H}|\Psi_{M}\rangle/\langle\Psi_{M}|\Psi_{M}\rangle$ is a functional of the network weights $\mathcal{W}$. In order to obtain an optimal solution for which $\nabla E(\mathcal{W}^{\star})=0$, several optimization approaches can be used. Here, we have found convenient to adopt the Stochastic Reconfiguration (SR) method of Sorella et al. sorella2007weakbinding , which can be interpreted as an effective imaginary-time evolution in the variational subspace. Introducing the variational derivatives with respect to the $k$-th network parameter,

the SR updates at the $p-$th iteration are of the form

where we have introduced the (positive-definite) covariance matrix

and a scaling parameter $\gamma(p)$. Since the covariance matrix can be non-invertible, $S^{-1}$ denotes its Moore-Penrose pseudo-inverse. Alternatively, an explicit regularization can be applied, of the form $S_{k,k^{\prime}}^{\text{reg}}=S_{k,k^{\prime}}+\lambda(p)\delta_{k,k^{\prime}}S_{k,k}$ . In our work we have preferred the latter regularization, with a decaying parameter $\lambda(p)=\max(\lambda_{0}b^{p},\lambda_{\text{min}})$ and typically take $\lambda_{0}=100$, $b=0.9$ and $\lambda_{\text{min}}=10^{-4}$.

Initially the network weights $\mathcal{W}$ are set to some small random numbers and then optimized with the procedure outlined above. In Fig. 5 we show the typical behavior of the optimization algorithm, which systematically approaches the exact energy upon increasing the hidden units density $\alpha$.

Appendix B Time-Dependent Variational Monte Carlo

In the second part of our Paper we have considered the problem of solving the many-body Schrödinger equation with a variational ansatz of the NQS form. This task can be efficiently accomplished by means of the Time-Dependent Variational Monte Carlo (t-VMC) method of Carleo et al.

are a functional of the variational parameters derivatives, $\dot{\mathcal{W}}(t)$, and can be interpreted as the quantum distance between the exactly-evolved state and the variationally evolved one. Since in general we work with unnormalized quantum states, the correct Hilbert-space distance is given by the Fubini-Study metrics, given by

where the correlation matrix and the forces are defined analogously to the previous section. In this case the diagonal regularization, in general, cannot be applied, and $S^{-1}(t)$ strictly denotes the Moore-Penrose pseudo-inverse.

The outlined procedure is globally stable as also already proven for other wave functions in past works using the t-VMC approach. In Fig. 6 we show the typical behavior of the time-evolved physical properties of interest, which systematically approach the exact results when increasing $\alpha$.

Appendix C Efficient Stochastic Sampling

We complete the supplementary information giving an explicit expression for the variational derivatives previously introduced and of the overall computational cost of the stochastic sampling. We start rewriting the NQS in the form

In our stochastic procedure, we generate a Markov chain of many-body configurations $\mathcal{S}^{(1)}\rightarrow\mathcal{S}^{(2)}\rightarrow\dots\mathcal{S}^{(P)}$ sampling the square modulus of the wave function $\left|\Psi_{M}(\mathcal{S})\right|^{2}$ for a given set of variational parameters. This task can be achieved through a simple Metropolis-Hastings algorithm metropolis1953equation , in which at each step of the Markov chain a random spin $s$ is flipped and the new configuration accepted according to the probability

In order to efficiently compute these acceptances, as well as the variational derivatives, it is useful to keep in memory look-up tables for the effective angles $\theta_{j}(\mathcal{S}^{(k)})$ and update them when a new configuration is accepted. These are updated according to

when the spin $s$ has been flipped. The overall cost of a Monte Carlo sweep (i.e. of $\mathcal{O}(N)$ single-spin flip moves) is therefore $\mathcal{O}(N\times M)=\mathcal{O}(\alpha N^{2}).$ Notice that the computation of the variational derivatives comes at the same computational cost as well as the computation of the local energies after a Monte Carlo sweep.

Appendix D Iterative Solver

The most time-consuming part of both the SR optimization and of the t-VMC method is the solution of the linear systems (6 and 11) in the presence of a large number of variational parameters $N_{\textrm{var}}$. Explicitly forming the correlation matrix $S$, via stochastic sampling, has a dominant quadratic cost in the number of variational parameters, $\mathcal{O}(N_{\textrm{var}}^{2}\times N_{\textrm{MC}})$, where $N_{\textrm{MC}}$ denotes the number of Monte Carlo sweeps. However, this cost can be significantly reduced by means of iterative solvers which never form the covariance matrix explicitly. In particular, we adopt the MINRES-QLP method of Choi and Saunders choi2014algorithm , which implements a modified conjugate-gradient iteration based on Lanczos tridiagonalization. This method iteratively computes the pseudo-inverse $S^{-1}$ within numerical precision. The backbone of iterative solvers is, in general, the application of the matrix to be inverted to a given (test) vector. This can be efficiently implemented due to the product structure of the covariance matrix, and determines a dominant complexity of $\mathcal{O}(N_{\textrm{var}}\times N_{\textrm{MC}})$ operations for the sparse solver. For example, in the most challenging case when translational symmetry is absent, we have $N_{\textrm{var}}=\alpha N^{2}$, and the dominant computational cost for solving (6 and 11) is in line with the complexity of the previously described Monte Carlo sampling.

Appendix E Implementing Symmetries

Very often, physical Hamiltonians exhibit intrinsic symmetries which must be satisfied also by their ground- and dynamically-evolved quantum states. These symmetries can be conveniently used to reduce the number of variational parameters in the NQS.

where the network weights have now a different dimension with respect to the standard NQS. In particular, $a^{(f)}$ and $b^{(f)}$ are vectors in the feature space with $f=1,\dots\alpha_{s}$ and the connectivity matrix $W_{j}^{(f)}$ contains $\alpha_{s}\times N$ elements. Notice that this expression corresponds effectively to a standard NQS with $M=S\times\alpha_{s}$ hidden variables. Tracing out explicitly the hidden variables, we obtain

In the specific case of site translation invariance, we have that the symmetry group has an orbit of $S=N$ elements. For a given feature $f$, the matrix $W_{j}^{(f)}$ can be seen as a filter acting on the $N$ translated copies of a given spin configuration. In other words, each feature has a pool of $N$ associated hidden variables that act with the same filter on the symmetry-transformed images of the spins.