Gaussian Process States: A Data-Driven Representation of Quantum Many-Body Physics

We present a novel, nonparametric form for compactly representing entangled many-body quantum states, which we call a “Gaussian process state.” In contrast to other approaches, we define this state explicitly in terms of a configurational data set, with the probability amplitudes statistically inferred from this data according to Bayesian statistics. In this way, the nonlocal physical correlated features of the state can be analytically resummed, allowing for exponential complexity to underpin the ansatz, but efficiently represented in a small data set. The state is found to be highly compact, systematically improvable, and efficient to sample, representing a large number of known variational states within its span. It is also proven to be a “universal approximator” for quantum states, able to capture any entangled many-body state with increasing data-set size. We develop two numerical approaches which can learn this form directly—a fragmentation approach and direct variational optimization—and apply these schemes to the fermionic Hubbard model. We find competitive or superior descriptions of correlated quantum problems compared to existing state-of-the-art variational ansatzes, as well as other numerical methods.

Popular Summary

Simulating many interacting quantum particles is a problem at the heart of a broad range of key scientific challenges, from the promise of chemistry by computational design to manipulating the emergent properties of novel materials. The central quantum variable, the wave function, is, however, an exponentially complex object, with amplitudes associated with every possible classical arrangement of the particles in the system. Therefore, the history of advances in interacting quantum systems is often punctuated by the discovery of accurate and compact representations of these quantum states. We introduce a new representation of quantum states via a probabilistic, machine-learned model over the distribution of possible classical arrangements of particles.

The central paradigm of our work stems from asking the following question: If we knew the amplitudes on a small number of these classical configurations, what is the optimal statistical model to infer the amplitudes on all other configurations, in order to define the true quantum state? We show that this “Gaussian process state” explicitly encompasses many other established wave-function forms in a highly compact manner, and we establish numerical algorithms to allow optimization of these chosen classical configurations in order to compactly find the most accurate state for a general interacting quantum system.

This new framework has the potential to impact a diverse range of fields, from quantum-state tomography, which reconstructs quantum states from limited experimental measurements, to the computational modeling of quantum chemistry and condensed matter physics.

Figure 1

Exemplification of the kernel function feature space. The kernel counts correlation features common to the test configuration (bottom) and training configuration (top), including, e.g., the single-site Gutzwiller feature (green), two-site features explicit in the Jastrow wave function (blue), and nearest-neighbor plaquette features used in EPS ansatzes (red).

Figure 2

Left panel: Mean-squared GPS amplitude error for the ground state of the 10-site, half-filled 1D Hubbard chain at $U=8t$, as the kernel power hyperparameter is varied ($p$, with $\theta =0$). Also shown is the size of the training set of configurations selected by the RVM from the exact data to construct the GPS with this underlying kernel complexity. Right panel: Real part of the spin structure factor $\mathrm{Re}\mathbf{(}S(k)\mathbf{)}$ corresponding to the exact wave function (black dotted line), the GPS with $p=0$ and $p=1$ (red solid line), and the single Slater determinant representation obtained with the Hartree-Fock method (green dashed line) for the same system. The inset shows the absolute error of $\mathrm{Re}\mathbf{(}S(k=\pi )\mathbf{)}$ for the GPS with respect to the kernel hyperparameter $p$, as well as the error obtained for a single Slater determinant representation.

Figure 3

Two-dimensional t-SNE projection of the Hilbert space of a half-filled 1D Hubbard chain of eight sites. The four clusters are specified by the occupancy of the reference site of the configurations (clockwise from the top left: spin-up electron, holon, doublon, and spin-down electron). The colored points are the 30 configurations selected by the RVM as the most relevant, with red and blue colors indicating positive and negative weights.

Figure 4

Relative energy errors for the half-filled Hubbard chain (compared to DMRG) at $U=8t$, $L=32$, as a function of the size of the chosen set of configurations ($N_b$), as the kernel complexity $p$ is increased from 0 to 6 ($\theta =0,L^{\prime }=6$). For other methods, $N_{\text{par}}$ denotes the number of variational parameters optimized in the state, as a measure of compactness and variational freedom of the state. The GPS is trained on the features of the ground state of a small 12-site lattice using configurations with all particle and spin fluctuations from an $L^{\prime }=6$ lattice. Subsequent variational optimization of the training weights of the GPS gives rise to the “GPS (opt)” values. For comparison, variationally optimized and translationally symmetric Gutzwiller, Jastrow, and RBM wave-function results (with the same fixed free-fermion state) are included, with the RBM ratio of hidden to visible nodes ($\alpha$) chosen as 1, 2, 4, 8, and 16 [23], taking the best result of 10 independent seeds. Also included are results from iterative bootstrapping of the GPS, as introduced in Sec. 3b. The bootstrapped results are obtained by starting from a random initial set of data configurations. For these data points, the $p\to \infty$ limit of the kernel is used where the scaled distance weights $f(d)/\theta$ in the kernel function are optimized variationally.

Figure 5

Relative energy error (compared to DMRG) against $U/t$ for GPS, Gutzwiller, Jastrow, and RBM states, for a half-filled 32-site Hubbard chain. The GPS ansatz is obtained from a fit on a 12-site chain solution with training configurations coming from an $L^{\prime }=6$ lattice with all possible particle and spin fluctuations using kernel hyperparameters ($p=4$, $\theta =0$). All ansatze are built on the same fixed free-fermion orbitals. Monte Carlo optimizations for the Gutzwiller, Jastrow, and RBM (with $\alpha =5$) states are performed for 10 seeds, with the variationally lowest result shown. The red curve shows the GPS result, with the weights solely obtained from the fit on the fragment amplitudes; the blue curve represents additional variational optimization of the weights, with the $L^{\prime }=6$ data set fixed.

Figure 6

Illustration of the bootstrapping algorithm for the 2D target function shown in the top-left corner. (1) An initial set of configurations is randomly selected (blue points), with $w_t=0$, ensuring no features are described. (2) The weights of each configuration are optimized variationally, reasonably describing the desired function. (3) The RVM is used to drastically sparsify the data set without changing the learned function. (4) Configurations with high local energy variance are added to the sparsified set (green crosses). Note that the final three panels model the function to essentially the same accuracy, allowing systematic improvement in refining the data set.

Figure 7

Energy errors obtained for the bootstrapped GPS for the half-filled $6\times 6$ 2D square Hubbard lattice at $U=8t$ (compared to benchmark AFQMC values [44]), plotted against the number of variational parameters (excluding those in the Pfaffian state). Results are systematically improved via use of a more complex kernel or a decrease of the ${\tilde{\sigma }}^2$ regularization parameter, which increases the chosen data-set size. The GPS with blue background shading was obtained with the $k^{(1)}$ kernel; values with brown shading correspond to the $p\to \infty$ kernel, where the distance weighting was fixed at $f(d)=1/|d|$ and $\theta$ was variationally optimized; for the pink highlighted results, all scaled distance hyperparameters $f(d)/\theta$ were optimized variationally (giving $N_b+35$ variational parameters). Also displayed for comparison are the energy errors obtained with the Gutzwiller, Jastrow, and RBM ansatzes (with $\alpha$ chosen as 1, 2, 4, 8, and 16). Each of these displayed comparison values corresponds to the lowest energy obtained for repeated calculations with five different random seeds.

Figure 8

Convergence of the ground-state energy for the $8\times 8$ Hubbard model, at $U=8t$ and with 0.875 electrons per site. The bootstrapped GPS representation is initialized at the optimized Jastrow wave function, resulting in rapid convergence. Vertical grey lines indicate the iterations at which the configurational data set of the GPS is updated. The convergence is shown for Gutzwiller (1 parameter), Jastrow (34 parameters), the RBM with $\alpha =16$ (2064 parameters), and the bootstrapped GPS with optimized distance weights ($N_b+63$ parameters). The number of data configurations specifying the GPS after each application of the RVM in the algorithm range from about 1000 to 2000 configurations.