Self-learning projective quantum Monte Carlo simulations guided by restricted Boltzmann machines

TL;DR

This paper introduces a self-learning projective quantum Monte Carlo method guided by restricted Boltzmann machines, which adaptively optimizes the wave function during simulation, eliminating the need for separate variational steps.

Contribution

The novel approach integrates machine learning into PQMC, enabling adaptive wave function optimization directly during simulations, improving accuracy and efficiency.

Findings

Accurate ground state simulations of the quantum Ising chain

High agreement with Jordan-Wigner and loop quantum Monte Carlo results

Eliminates separate variational optimization step

Abstract

The projective quantum Monte Carlo (PQMC) algorithms are among the most powerful computational techniques to simulate the ground state properties of quantum many-body systems. However, they are efficient only if a sufficiently accurate trial wave function is used to guide the simulation. In the standard approach, this guiding wave function is obtained in a separate simulation that performs a variational minimization. Here we show how to perform PQMC simulations guided by an adaptive wave function based on a restricted Boltzmann machine. This adaptive wave function is optimized along the PQMC simulation via unsupervised machine learning, avoiding the need of a separate variational optimization. As a byproduct, this technique provides an accurate ansatz for the ground state wave function, which is obtained by minimizing the Kullback-Leibler divergence with respect to the PQMC samples, rather than by minimizing the energy expectation value as in standard variational optimizations. The high accuracy of this self-learning PQMC technique is demonstrated for a paradigmatic sign-problem-free model, namely, the ferromagnetic quantum Ising chain, showing very precise agreement with the predictions of the Jordan-Wigner theory and of loop quantum Monte Carlo simulations performed in the low-temperature limit.

Equations

$$\hat{H}=-J\sum _{j=1}^{N-1}{\sigma }_j^z{\sigma }_{j+1}^z-\Gamma \sum _{j=1}^N{\sigma }_j^x-h\sum _{j=1}^N{\sigma }_j^z,$$

$$f(\mathbf{x},\tau +\Delta \tau )=\sum _{{\mathbf{x}}^{\prime }}\tilde{G}(\mathbf{x},{\mathbf{x}}^{\prime },\Delta \tau )f({\mathbf{x}}^{\prime },\tau ).$$

$$E_{\tau \to \infty }\equiv \frac{{\sum }_{\mathbf{x}}f\left(\mathbf{x},\tau \to \infty \right)E_{\mathrm{loc}}(\mathbf{x})}{{\sum }_{\mathbf{x}}f\left(\mathbf{x},\tau \to \infty \right)}.$$

$$H_{\mathrm{RBM}}\left(\mathbf{x},\mathbf{h}\right)=-\sum _{i,j}{J}_{ij}{h}_i{x}_j-\sum _j{a}_j{x}_j-\sum _i{b}_i{h}_i.$$

$$P(\mathbf{x})\equiv \sum _{\mathbf{h}}P(\mathbf{x},\mathbf{h})=\frac{1}{Z}\sum _{\mathbf{h}}\exp \left[-H_{\mathrm{RBM}}(\mathbf{x},\mathbf{h})\right].$$

$$\mathrm{KL}\left(q\left|\right|p\right)=\sum _{\mathbf{x}}q(\mathbf{x})\ln \left(q(\mathbf{x})/p(\mathbf{x})\right).$$

$$\frac{\partial L(\mathbb{W})}{\partial J_{ij}}\propto {⟨x_j{h}_i⟩}_{\mathrm{data}}-{⟨x_j{h}_i⟩}_{\mathrm{model}}.$$

$$E_{{\psi }_{G_s}}\equiv \frac{⟨{\psi }_{G_s}|\hat{H}|{\psi }_{G_s}⟩}{⟨{\psi }_{G_s}|{\psi }_{G_s}⟩}.$$

Full text

I Introduction

The similarity between the imaginary-time Schrödinger equation and the diffusion equation allows one to simulate quantum many-body systems by stochastically evolving a (typically large) population of random walkers Anderson (1975). This is the basis of so-called projective quantum Monte Carlo algorithms (PQMC). These algorithms are particularly efficient when they simulate the ground state of sign-problem-free Hamiltonians, since the corresponding wave function is real and nonnegative in a suitable basis. Using fixed-node and released node methods Reynolds et al. (1982); Ceperley and Alder (1984, 1986); Booth et al. (2009), PQMC simulations provide accurate predictions even when the sign problem occurs, in particular for many-fermion systems, albeit in general at a larger computational cost. PQMC algorithms have been used to simulate fundamental quantum systems, including the electron gas Ceperley and Alder (1980), molecular systems Reynolds et al. (1982), liquid and solid Helium Boronat (2002), electrons in solids Foulkes et al. (2001); Booth et al. (2013), ultracold gases Giorgini et al. (1999), quantum spin and lattice models Trivedi and Ceperley (1990); Buonaura and Sorella (1998), and nuclear matter Carlson et al. (2015). However, these algorithms are efficient only if they are provided with a sufficiently accurate trial wave function, which is used to guide the random walkers towards the most relevant regions of the configuration space. If this guiding wave function is not included, or if it is not sufficiently accurate, the computational cost of a PQMC simulation scales exponentially with the system size Nemec (2010); Boninsegni and Moroni (2012); Inack et al. (2018a). For various relevant systems, the general form of a suitable guiding wave function is provided by some physical theory, and the details can be further tuned via a variational minimization of the energy expectation value. However, an appropriate theory is not always available, and the variational optimization might turn out to be an extremely challenging computational task. Even more, an inaccurate ansatz for the guiding wave function might lead to biased predictions, in particular when different phases closely compete, as in many strongly-correlated and/or disordered many-body systems. This problem is especially relevant in the field of adiabatic quantum computing, where quantum Monte Carlo (QMC) algorithms are being used to simulate how quantum annealers solve complex optimization problems Santoro et al. (2002); Boixo et al. (2014); Inack and Pilati (2015); Boixo et al. (2016); Heim et al. (2015). In fact, the Hamiltonians corresponding to typical instances of hard optimization problems can be written in the form of random Ising models. These models are characterized by glassy ground states, for which an accurate ansatz is hard to guess Stella and Santoro (2007).

Wave functions based on a certain type of generative neural network, specifically, restricted Boltzmann machines (RBM) [see Fig. 1, upper panel], have been recently proposed as generic variational ansatzes for quantum spin models Carleo and Troyer (2017). Their representational power, their entanglement properties, and more elaborate versions are currently the subject of very intense research activity Deng et al. (2017); Gao and Duan (2017); Saito (2017); Chen et al. (2018); Glasser et al. (2018); Freitas et al. (2018); Cai and Liu (2018); Saito and Kato (2018); Carleo et al. (2018); Luo and Clark (2018); Kochkov and Clark (2018); Ruggeri et al. (2018); Beach et al. (2019); McBrian et al. (2019); Collura et al. (2019); Kessler et al. (2019); Nagy and Savona (2019). They appear to be particularly useful when no other physics-inspired ansatz is available. In this paper, we show how to employ them as guiding wave functions in PQMC simulations. In a previous article Inack et al. (2018b), a neural network state that mimics a type of an unrestricted Boltzmann machine (see Fig. 1, lower panel) has been employed as guiding wave function. In the case of Ref. Inack et al. (2018b), the neural network parameters had to be obtained via a separate variational optimization, and the PQMC algorithm had to be overloaded with the sampling of additional hidden variables. Here we show that, using unsupervised machine learning algorithms, the RBM can be trained along the PQMC simulations directly from random-walker configurations. One obtains an adaptive guiding wave function that improves during the PQMC simulation, allowing one to completely eliminate the bias due to the finite random-walker population Inack et al. (2018a) and to drastically reduce the statistical fluctuations. As a testbed, we consider a paradigmatic sign-problem free model, namely, the ferromagnetic quantum Ising chain, making comparison against the exact predictions of the Jordan-Wigner theory and of loop quantum Monte Carlo (QMC) simulations performed with the ALPS library Bauer et al. (2011) in the low-temperature limit. Excellent agreement is obtained for ground-state energies and for the average magnetization, above, below, and at the ferromagnetic quantum critical point. A (possibly useful) byproduct of the self-learning procedure described here is an optimized RBM ansatz. This ansatz is developed by minimizing a Kullback-Leibler divergence, rather than the common minimization of the variational energy Torlai et al. (2018). In the model considered here, the self-learned RBM wave function turns out to be comparatively as accurate as the one obtained via variational minimization using a sophisticated optimization method such as the stochastic reconfiguration algorithm Sorella (1998) as implemented in the NetKet library Carleo et al. (2019).

The paper is organized as follows. Section II introduces the model Hamiltonian we address and describes the PQMC algorithm, the RBM variational wave-function, as well as the adaptive unsupervised learning procedure. Numerical results and comparison with previous theories are presented in Section III. Section IV summarizes our findings and discusses possible future extensions.

II Method

As a test bed for the self-learning PQMC simulations we consider the one-dimensional ferromagnetic quantum Ising Hamiltonian, defined as:

[TABLE]

where $\sigma^{x}_{j}$ and $\sigma^{z}_{j}$ indicate conventional Pauli matrices acting on the spin at the lattice site $j=1,\dots,N$. $N$ is the total number of spins, and we consider open boundary conditions. The parameter $J>0$ fixes the strength of the ferromagnetic interactions between nearest-neighbor spins. In the following, we use $J=1$ as unit of the energy scale. $\Gamma$ and $h$ are the intensities of the transverse and longitudinal magnetic field, respectively. Below, the eigenstates of the Pauli matrix ${\sigma^{z}_{j}}$ are denoted as $\left|x_{j}\right>$. The eigenvalue is $x_{j}=1$ when $\left|x_{j}\right>=\left|\uparrow\right>$ and $x_{j}=-1$ when $\left|x_{j}\right>=\left|\downarrow\right>$. The quantum states of $N$ spins $\left|\mathit{\boldsymbol{x}}\right>=\left|x_{1}x_{2}...x_{N}\right>$, with $\mathit{\boldsymbol{x}}=(x_{1},\dots,x_{N})$, form the computational basis considered in the present work. We denote as $\left|\psi\right>$ the quantum state corresponding to the wave-function $\left<x\right|\left.\psi\right>=\psi(\mathit{\boldsymbol{x}})$.

II.1 Projective Quantum Monte Carlo simulations

The PQMC algorithms allow one to simulate the ground-state of a generic Hamiltonian by stochastically evolving the imaginary-time Schrödinger equation. The algorithm’s accuracy and efficiency greatly improve if one introduces a suitable ansatz for the ground-state wave-function — usually called guiding (or trial) wave function and indicated as $\psi_{G}({\mathit{\boldsymbol{x}}})$ — and let the product $f({\mathit{\boldsymbol{x}}},\tau)=\psi({\mathit{\boldsymbol{x}}},\tau)\psi_{G}({\mathit{\boldsymbol{x}}})$ evolve according to the modified imaginary-time Schrödinger equation. This equation reads:

[TABLE]

Here, $\tilde{G}({\mathit{\boldsymbol{x}}},{\mathit{\boldsymbol{x}}}^{\prime},\Delta\tau)=G(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{x}}^{\prime},\Delta\tau)\frac{\psi_{G}({\mathit{\boldsymbol{x}}})}{\psi_{G}({\mathit{\boldsymbol{x}}^{\prime}})}$, where $G({\mathit{\boldsymbol{x}}},{\mathit{\boldsymbol{x}}}^{\prime},\Delta\tau)=\left<\mathit{\boldsymbol{x}}\left|\exp\left[{-\Delta\tau(\hat{H}-E_{\mathrm{ref}})}\right]\right|\mathit{\boldsymbol{x}}^{\prime}\right>$ is the imaginary-time Green’s function for a (short) time step $\Delta\tau$ (we set $\hbar=1$ throughout this paper). $E_{\mathrm{ref}}$ is a reference energy introduced to stabilize the numerics, as explained below.

We employ the continuous-time algorithm of Refs. Becca and Sorella (2017); Sorella and Capriotti (2000), which allows one to stochastically simulate the exact (modified) imaginary-time Green’s function $\tilde{G}({\mathit{\boldsymbol{x}}},{\mathit{\boldsymbol{x}}}^{\prime},\Delta\tau)$, avoiding any finite time-step error which would occur with the use of Trotter approximations. This is achieved by using the formal zero-time limit of the Green’s function and by appropriately sampling the time interval $\delta\tau$ between consecutive single-spin flips. More flips might occur within a time-step $\Delta\tau$. Since the modified imaginary-time Green’s function does not define a standard Markov process, i.e., one driven by a column normalized transition matrix, the simulation has to be performed by stochastically evolving a (large) population of random walkers; these walkers are subjected to updates in configuration space and to a branching process in which they are annihilated or replicated. The spin-configuration updates $\mathit{\boldsymbol{x}}^{\prime}\rightarrow\mathit{\boldsymbol{x}}$ (with $\mathit{\boldsymbol{x}}^{\prime}\neq\mathit{\boldsymbol{x}}$) are randomly selected from a probability distribution proportional to the modified imaginary-time Green’s function. The guiding wave function favors updates towards relevant regions of the configuration space. In the branching process, the random walkers are annihilated or replicated according to the weight factor $w_{\mathit{\boldsymbol{x}}^{\prime}}=\exp\left[{-\delta\tau(E_{\mathrm{loc}}(\mathit{\boldsymbol{x}}^{\prime})-E_{\mathrm{ref}})}\right]$, where the local energy is $E_{\mathrm{loc}}(\mathit{\boldsymbol{x}}^{\prime})=\sum_{\mathit{\boldsymbol{x}}}H_{\mathit{\boldsymbol{x}}^{\prime},\mathit{\boldsymbol{x}}}\frac{\psi_{G}({\mathit{\boldsymbol{x}}})}{\psi_{G}({\mathit{\boldsymbol{x}}^{\prime}})}$ 111One has to accumulate the weight factors of the updates occurring within the time-step $\Delta\tau$. See Refs. Sorella and Capriotti (2000); Becca and Sorella (2017); Inack et al. (2018b) for more details.. This weight accounts for the column normalization of the transition matrix. By dynamically tuning $E_{\mathrm{ref}}$, the size of the random-walker population can be kept very close to a desired target value $N_{\mathrm{w}}$. To implement the branching process and for the tuning of $E_{\mathrm{ref}}$ we use the textbook recipe of Ref. Thijssen (2007).

The continuous-time PQMC algorithm sketched here is more exhaustively described in Refs. Becca and Sorella (2017); Sorella and Capriotti (2000); Inack et al. (2018b), and we refer the interested readers to those references for additional implementation details. The projection can reach long imaginary times $\tau=N_{\mathrm{MC}}\Delta\tau$ by iterating a large number $N_{\mathrm{MC}}$ of Monte Carlo steps, each step corresponding to a short time step $\Delta\tau$. In the long imaginary-time limit, the walkers sample spin configurations with a probability distribution proportional to $f(\mathit{\boldsymbol{x}},\tau\rightarrow\infty)=\psi_{0}({\mathit{\boldsymbol{x}}})\psi_{G}({\mathit{\boldsymbol{x}}})$ (if $N_{w}$ is large enough, as explained below), where $\psi_{0}({\mathit{\boldsymbol{x}}})$ is the ground-state wave function. Unbiased predictions of ground-state energies are obtained via Monte Carlo integration of the sum

[TABLE]

Unbiased prediction of ground-state expectation values of other observables that commute with the Hamiltonian are obtained with analogous formulas. However, for operators $\hat{O}$ that do not commute with the Hamiltonian, the analogous formulas would predict so-called mixed estimators, namely $\left<\psi_{0}\left|\hat{O}\right|\psi_{G}\right>/\left<\psi_{0}\right|\left.\psi_{G}\right>$. In general, these are affected by a bias due to the guiding wave function, unless the latter coincides with the ground-state wave function. Nevertheless, if the operator $\hat{O}$ is diagonal in the chosen computational basis, the pure estimator corresponding to $\left<\psi_{0}\left|\hat{O}\right|\psi_{0}\right>/\left<\psi_{0}\left|\right.\psi_{0}\right>$ can be determined via the standard forward-walking technique (see, e.g., Ref. Boronat (2002)). An accurate guiding wave function reduces the computational cost of the forward walking technique.

In the large random-walker population limit, $N_{\mathrm{w}}\rightarrow\infty$, the above Monte Carlo estimates are unbiased, being affected only by statistical fluctuations, which can be systematically reduced. For finite $N_{\mathrm{w}}$, a systematic bias might arise due to the spurious correlations among replicated walkers Hetherington (1984); Golinelli et al. (1992); Nemec (2010); Boninsegni and Moroni (2012); Inack et al. (2018a); Pollet et al. (2018). In Ref. Inack et al. (2018a), it was indeed shown that if one does not introduce a guiding wave function, which corresponds to setting $\psi_{G}({\mathit{\boldsymbol{x}}})=1$ in the equations above, the random-walker population required to keep this systematic bias below a chosen (small) threshold increases exponentially with the system size. This, in turn, implies a computational cost which scales exponentially. Instead, if the guiding wave function is exact, i.e. if $\psi_{G}({\mathit{\boldsymbol{x}}})=\psi_{0}({\mathit{\boldsymbol{x}}})$, the local energy $E_{\mathrm{loc}}(\mathit{\boldsymbol{x}})$ is a constant function. This freezes the annihilations and replications of the random walkers in the branching process, eliminating any potential bias. If $\psi_{G}({\mathit{\boldsymbol{x}}})$ is a reasonable approximation for $\psi_{0}({\mathit{\boldsymbol{x}}})$, the fluctuations of the random-walker number are anyway reduced compared to the case of the simple PQMC algorithm performed with $\psi_{G}({\mathit{\boldsymbol{x}}})=1$, giving a much faster convergence to the exact $N_{w}\rightarrow\infty$ limit.

II.2 Boltzmann machines for PQMC algorithms

In Ref. Inack et al. (2018b), it has been shown that a guiding wave function based on an unrestricted Boltzmann machine is sufficiently accurate to drastically reduce the required random-walker population, leading to a polynomially scaling computational cost, at least for the ferromagnetic quantum Ising chain. Remarkably, this has been achieved by optimizing just three variational parameters. The small number of variational parameters to be optimized is the main benefit of the unrestricted architecture of Refs. Inack et al. (2018b); Collura et al. (2019). However, using the unrestricted Boltzmann machine requires to sample additional hidden spins, both during the variational optimization 222In one and in two dimensional models with short-range interactions, the unrestricted Boltzmann machine ansatzes can be mapped to constrained matrix product states and to constrained tensor network states, respectively Collura et al. (2019). In one dimension this mapping allows one to implement alternative variational minimization techniques, which avoid the sampling over hidden spins. and, chiefly, during the PQMC simulations. In fact, the PQMC algorithm has to be extended, combining visible-spin and hidden-spin sampling as explained in Refs. Vitiello and Whitlock (1991); Inack et al. (2018b). An inefficient sampling of the additional hidden spins might introduce further correlations among replicated walkers, leading to a larger finite-$N_{w}$ bias. This deleterious effect has indeed been observed in PQMC simulation performed close to the ferromagnetic quantum critical point Inack et al. (2018b), where the statistical correlations among subsequent hidden-spin configurations are larger. This effect did not lead to an exponentially scaling computational cost in the ferromagnetic quantum Ising chain, but it might be more detrimental for more challenging Hamiltonians.

Below, we show that a guiding wave function based on an RBM allows one to eliminate the finite-$N_{\mathrm{w}}$ bias. The additional benefit the RBM ansatz provides is that the optimization of the guiding function can be performed via unsupervised machine learning directly from the random-walker population, with no need for a separate variational optimization. Furthermore, as opposed to the unrestricted Boltzmann machines, the RBM ansatz does not require to sample additional classical hidden-spin variables. This allows one to avoid overloading the PQMC algorithm with the Monte Carlo updates for the hidden spins. We argue that these two benefits might be particularly relevant in two setups: in simulation of disordered models, where frustration effects might increase the correlations among successive hidden spin configurations, and in quantum annealing simulations, where the Hamiltonian varies in time and an automatically adaptive guiding wave function would allow one to perform more efficient simulations of the quantum annealers’ dynamics.

Boltzmann machines are generative stochastic neural networks often used to approximate the probability distribution corresponding to a given population of stochastic samples. Beyond the visible spin variables $\mathit{\boldsymbol{x}}$, they include $N_{h}$ additional hidden spin variables $h_{i}=\pm 1$, with $i=1,\dots,N_{h}$, collectively denoted in the following as $\mathit{\boldsymbol{h}}=(h_{1},\dots,h_{N_{h}})$. The probability distribution is written in the form of the Boltzmann weight $P(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{h}})$ corresponding to a classical Ising Hamiltonian. In the case of the RBM, this classical Hamiltonian function reads:

[TABLE]

Notice that no visible-visible interactions nor hidden-hidden interactions are allowed. Instead, all visible spins interact with all hidden spins via the coupling parameters $J_{ij}$. One obtains two layers, respectively referred to as visible and hidden layer, with the bipartite all-to-all connectivity illustrated in the upper panel of Fig. 1. The parameters $a_{j}$ and $b_{i}$, called biases, play the role of longitudinal local magnetic fields. The couplings $J_{ij}$, together with the bias terms $a_{j}$ and $b_{i}$, define the RBM. In the following, all these parameters will be collectively denoted as $\mathbb{W}=\{W_{m}\}=\{J_{ij},a_{j},b_{i}\}$, for $i=1,\dots,N_{h}$ and $j=1,\dots,N$. The index $m$ labels all the RBM parameters. In practical applications the required number of hidden units is typically $N_{h}\sim N$. Thus, the number of parameters scales, to leading order, as $N^{2}$. This quadratic scaling can be reduced to a linear scaling in translationally invariant models Sohn and Lee (2012); Carleo and Troyer (2017). However, in view of future studies on disordered systems, we avoid translational invariance adopting open boundary conditions. This allows us to test the unsupervised learning of RBM wave functions in the most generic setup. In the unrestricted Boltzmann machine considered in Ref. Inack et al. (2018b), intra-layer visible-visible and hidden-hidden direct interactions were included, reflecting the same (nearest-neighbor) connectivity of the quantum Hamiltonian in Eq. 1. Furthermore, $N_{h}$ was taken equal to $N$, and inter-layer interactions were allowed only between visible and hidden spins corresponding to the same index, i.e., for $i=j$. The resulting connectivity structure is shown in the lower panel of Fig. 1. This structure is analogous to that of the shadow wave-function used to describe the liquid and solid phases of Helium-4 Reatto and Masserini (1988); Vitiello et al. (1988).

In general, the probability to sample a visible-spin configuration $\mathit{\boldsymbol{x}}$ is the marginal distribution over all possible hidden-spin configurations $\mathit{\boldsymbol{h}}$:

[TABLE]

Notice that the fictitious temperature and the Boltzmann constant $k_{B}$ are here equal to unity. The normalization factor is the partition function $Z=\sum_{\mathit{\boldsymbol{x}},\mathit{\boldsymbol{h}}}\exp\left[-H_{{\rm RBM}}\left(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{h}}\right)\right]$. Due to the absence of intra-layer interactions in the RBM, the hidden-spin configurations can be analytically traced out, resulting in a marginal distribution $P(\mathit{\boldsymbol{x}})\propto\exp\left(\sum_{j}a_{j}x_{j}\right)\prod_{i}F_{i}(\mathit{\boldsymbol{x}})$, where $F_{i}(\mathit{\boldsymbol{x}})=2\cosh\left[b_{i}+\sum_{j}J_{ij}x_{j}\right]$. In Ref. Carleo and Troyer (2017), it has been proposed to use the function $P(\mathit{\boldsymbol{x}})$ to define an (unnormalized) ground-state wave-function. In order to describe both amplitude and phase, the RBM parameters $\bf{W}$ should, in general, be complex valued. However, we consider here models whose ground-state wave function can be assumed to be real and nonnegative in a suitable basis; therefore, the RBM parameters can be restricted to have real parameters. Extensions to complex-valued ground-states for, e.g., fermionic systems, have recently been addressed Luo and Clark (2018); Choo et al. (2019); Ferrari et al. (2019). In Ref. Carleo and Troyer (2017), the RBM parameters $\bf{W}$ were determined via variational minimization of the expectation value $\left<\psi_{{\rm RBM}}\left|\hat{H}\right|\psi_{{\rm RBM}}\right>/\left<\psi_{{\rm RBM}}\left|\right.\psi_{{\rm RBM}}\right>$, where $\left|\psi_{{\rm RBM}}\right>$ indicates the quantum state corresponding to the (unnormalized) wave-function $\left<\mathit{\boldsymbol{x}}|\psi_{{\rm RBM}}\right>\propto P(\mathit{\boldsymbol{x}})$. The optimization of the variational parameters was performed using the stochastic reconfiguration algorithm Sorella and Capriotti (2000). In the following, the corresponding minimal variational energy, obtained using the Netket Carleo et al. (2019) library, will be indicated as $E_{\mathrm{var.min.}}$. In Ref. Carleo and Troyer (2017), the variational-minimization approach has been referred to as reinforcement learning, due to the close similarity with the reinforcement-learning techniques used in the field of machine learning. However, in typical machine learning applications of RBMs the parameters $\bf{W}$ are determined via unsupervised machine-learning techniques Ackley et al. (1985), trying to learn the (unknown) distribution corresponding to a (typically large) dataset of stochastic samples. The unsupervised learning approach has already been used in Ref. Torlai et al. (2018) to perform quantum states tomography. Here, we adopt it to extract reasonably accurate variational ansatzes from the random-walker population generated by the PQMC algorithm. As discussed below, these ansatzes can then be used in an adaptive scheme as guiding wave functions to boost the efficiency of the PQMC simulation itself. On the one hand, this scheme allows one to eliminate the bias originating from the finite $N_{w}$. On the other hand, it allows one to eliminate the residual error in the RBM ansatz, providing unbiased predictions of ground-state properties.

II.3 Unsupervised learning of adaptive RBM guiding wave functions

In the unsupervised learning approach, the RBM parameters are determined by maximizing the log-likelihood of a (typically large) training set: $L(\mathbb{W})=\sum_{w}\ln P(\mathit{\boldsymbol{x}}_{w})$, where $w$ labels the instances in the training set. It can be shown that this corresponds to the minimization of the so-called Kullback-Leibler divergence (see, e.g., Fischer and Igel (2012)). In general, the Kullback-Leibler divergence between two distributions $p(\mathit{\boldsymbol{x}})$ and $q(\mathit{\boldsymbol{x}})$ is defined as:

[TABLE]

It represents a measure of the distance between two distributions (non-symmetric with respect to the exchange $q\leftrightarrow p$). In the case discussed here, $q(\mathit{\boldsymbol{x}})$ is identified with the distribution of the random walkers $f(\mathit{\boldsymbol{x}},\tau\rightarrow\infty)$ obtained via PQMC simulations, while $p(\mathit{\boldsymbol{x}})$ with the RBM marginal distribution $P(\mathit{\boldsymbol{x}})$. The optimization of the RBM parameters $\bf{W}$ can be performed using the gradient ascent algorithm. It consists in performing iterative updates starting from an initial (random) guess $\mathbb{W}^{0}$. At the step $n=0,1,\dots,N_{\mathrm{steps}}$, one applies the rule: $W_{m}^{n+1}=W_{m}^{n}+\eta\frac{\partial}{\partial W_{m}}L(\mathbb{W}^{n})$, where the coefficient $\eta$ is the learning rate. This plain vanilla rule can be improved in various ways, e.g., by including a momentum term proportional to the update performed at the $n-1$ step, by annealing the learning rate, by adopting the adaptive gradient algorithm (AdaGrad) Duchi et al. (2011) or the adaptive moment estimation algorithm (Adam) Kingma and Ba (2014). For all results reported in this manuscript, the plain vanilla rule is augmented only by adding a momentum term corresponding to $\nu\frac{\partial}{\partial W_{m}}L(\mathbb{W}^{n-1})$, with $\nu$ a rate tuned as discussed below. Also the learning-rate annealing is adopted, as explained in the following.

The main task in the implementation of the gradient ascent algorithm is the computation of the gradients of the log-likelihood. It is common to use stochastic estimates computed on mini-batches of $N_{b}$ instances (typically, $N_{b}\sim 10-100$) randomly sampled from the (much larger) training set. The formula for the gradients with respect to the couplings $J_{ij}$ is usually stated in the following from:

[TABLE]

The first term on the right-hand size indicates the average obtained when the visible units are clamped to the data in the mini-batch. To determine its value, one uses the probability distribution of the hidden spins $h_{i}$. Due to the bipartite connectivity, it depends only on the visible-spin configurations. One has $h_{i}=1$ with probability $p_{h_{i}=1}(\mathit{\boldsymbol{x}})=1/\left[1+\exp\left(-2\sum_{j}x_{j}J_{ij}-2b_{i}\right)\right]$, and $h_{i}=-1$ with probability $1-p_{h_{i}=1}(\mathit{\boldsymbol{x}})$.

Determining the second term on the right-hand side of eq. (7) is more challenging. It represents the average obtained when both visible and hidden variables are sampled according to the probability distribution $P(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{h}})$ defined by the RBM model. In machine learning jargon, this term is often referred to as the dreaming phase of the learning process, in contrast to the first term which would correspond to the awake phase. The average can be determined via Monte Carlo estimation, starting from the visible variables corresponding to the mini-batch instances, and then alternating Gibbs sampling of all hidden variables with fixed visible variables, followed by the sampling of all visible variables with hidden variables fixed at the previously sampled values. The probability to sample $x_{j}=1$ is analogous to the formula for the hidden spins: $p_{x_{j}=1}(\mathit{\boldsymbol{h}})=1/\left[1+\exp\left(-2\sum_{i}h_{j}J_{ij}-2a_{j}\right)\right]$.

In principle, this alternated Gibbs sampling should be iterated till the Markov chain equilibrates. This would provide an unbiased estimate of the average. In practice, repeating a sufficient number of iterations to guarantee that equilibrium has been achieved is often computationally overwhelming. The algorithm corresponding to a finite number $k$ of iterations is referred to as $k$-step contrastive divergence, since it can be derived by minimizing the difference of two Kullback-Leibler divergences Hinton (2002). In the $k\rightarrow\infty$ limit, this algorithm is provably unbiased, meaning that it corresponds to maximizing the log-likelihood of the dataset Bengio and Delalleau (2009). However, it often turns out to be very accurate also for small $k$. Formulas for the log-likelihood derivatives with respect to the other RBM parameters, namely the biases $a_{j}$ and $b_{i}$, can also be derived, analogously to eq. (7). Rather than reporting them here, we provide in Algorithm 1 the detailed procedure, which is adapted from Ref. Fischer and Igel (2012) to the case of interest to us with binary values $x_{j},h_{i}=\pm 1$, rather than the values $1$ and [math], more common in the machine-learning literature. The input of the algorithm is a mini-batch of $N_{b}$ randomly sampled instances, while the output provides the partial derivatives of the log-likelihood $L_{W_{m}}\equiv\frac{\partial L(\mathbb{W})}{\partial W_{m}}$. It is convenient to normalize the output, dividing by $N_{b}$. For a recent review on the topic, the interested reader is referred to Ref. Torlai and Melko (2019).

In the adaptive scheme presented in this paper, the above unsupervised learning algorithm is used to train the RBM to describe the (unnormalized) probability distribution $f(\mathit{\boldsymbol{x}},\tau\rightarrow\infty)$ corresponding to a large random-walker population produced by a PQMC simulation at equilibrium. It is in fact known that, in principle, for a sufficiently large $N_{h}$ an RBM can approximate any discrete distribution (in the sense of the Kullback-Leibler divergence), and many researchers have found that the training algorithms described above are capable of finding optimal, or close to optimal, RBM parameters.

The scheme we propose here involves many consecutive stints, labelled in the following by the index $s=0,1,2,\dots,N_{{\rm stints}}-1$, each including a PQMC simulation with a guiding function $\psi_{G_{s}}(\mathit{\boldsymbol{x}})$ for a long imaginary time $\tau_{s}$, followed by a learning stage which is used to construct the new guiding function $\psi_{G_{s+1}}(\mathit{\boldsymbol{x}})$. The imaginary time $\tau$ runs up to $\tau_{s}N_{{\rm stints}}$. In each learning stage, occurring at the imaginary times $(s+1)\tau_{s}$, the RBM is trained to describe the equilibrium random-walker distribution produced during the last stint. For the initial PQMC stint $s=0$, a quite crude guiding wave function $\psi_{G_{s=0}}(\mathit{\boldsymbol{x}})$ can be chosen, e.g., a constant function or a wave function based on (the square root of) an RBM with random parameters, $\psi_{G_{s=0}}(\mathit{\boldsymbol{x}})\propto\sqrt{P(\mathit{\boldsymbol{x}})}$. After the $s$ stint, the learned RBM distribution is $P_{s}(\mathit{\boldsymbol{x}})\propto\psi_{G_{s}}(\mathit{\boldsymbol{x}})\psi_{0}(\mathit{\boldsymbol{x}})$ (assuming $N_{w}$ and $\tau_{s}$ large enough). Notice that the learning process converges much faster if in each learning stage the RBM parameters are initialized at the values found in the previous learning stage.

The new guiding function for the stint $s+1$ is obtained by setting $\psi_{G_{s+1}}(\mathit{\boldsymbol{x}})=\sqrt{P_{s}(\mathit{\boldsymbol{x}})}$. With this choice the accuracy of the guiding function increases stint after stint. Indeed, under the idealized conditions of $N_{h}$ sufficiently large, successful unsupervised optimization, and $N_{w}$ and $\tau_{s}$ large enough, at the $s$ stint the random walkers sample the (unnormalized) distribution $f(\mathit{\boldsymbol{x}},\tau\rightarrow\infty)=\psi_{G_{s=0}}^{1/2^{s}}(\mathit{\boldsymbol{x}})\psi_{0}^{2-1/2^{s}}(\mathit{\boldsymbol{x}})$. This indicates an extremely rapid convergence $\psi_{G_{s\rightarrow\infty}}(\mathit{\boldsymbol{x}})\rightarrow\psi_{0}(\mathit{\boldsymbol{x}})$. As a consequence, in less idealized conditions, the (possible) bias due to the finite $N_{w}$ is expected to disappear after a relatively small number of stints/training stages. In the next section, we demonstrate that this is indeed the case. This is achieved with an affordable number of hidden spins $N_{h}$, using a simple implementation of the stochastic gradient ascent algorithm described above.

II.4 Simulation details

For all results reported in this article, the simulation details are the following. The PQMC simulations are performed with a target number of walkers $N_{w}=10^{4}$ and time step $\Delta\tau=0.04$. In each stint, the PQMC simulation runs for an imaginary time $\tau_{s}=20$, which turns out to be sufficiently large to represent the infinite imaginary-time limit $\tau\rightarrow\infty$. A small (randomly selected) fraction of the walker population, namely $N_{w}/20$, is stored at each PQMC step for the learning stage, excluding the initial time segment of each stint corresponding to $\tau\in\left[s\tau_{s},s\tau_{s}+8\right]$. This avoids correlations and non-equilibrium effects. This protocol provides $N_{\mathrm{train}}\simeq 15\times 10^{4}$ instances for each unsupervised learning stage. The number of stints ranges from $N_{{\rm stints}}=20$ to $N_{{\rm stints}}=50$. This appears to be sufficient to approach the $s\rightarrow\infty$ limit. Unsupervised learning is performed after each PQMC stint using the $k$-step contrastive divergence algorithm, performing a number $N_{\mathrm{steps}}=5\times 10^{4}$ to $10^{5}$ of stochastic gradient ascent steps, computed on mini-batches of size $N_{b}=10$ to $50$. The coefficient of the momentum term is $\nu=\eta/10$ (see description above). The learning rate $\eta$ is kept fixed within each learning stage, but is reduced stage after stage following the simple empirical annealing protocol $\eta(s)=\eta_{0}0.75^{s}$, where $\eta_{0}=0.01$ is the learning rate at the first learning stage $s=0$. Remarkably, $k=1$ is found to suffice. Tests with $k\approx 30$ provide comparable results. This strongly suggests that in this problem the $k$-step contrastive divergence algorithm is correctly minimizing the Kullback-Leibler divergence. The initial guiding function $\psi_{G_{s=0}}(\mathit{\boldsymbol{x}})$ is the square root of an RBM distribution with uniform random couplings $J_{ij}\in\left[-0.025:0.025\right]$. In the absence of longitudinal magnetic field, i. e. with $h=0$, the bias terms are initialized to zero, while for finite $h$ they are set to uniform random values $a_{j},b_{i}\in\left[0:0.05\right]$. It is worth mentioning that, while the setup described here turns out to be adequate for the problems addressed in the present work, it is plausible that even more efficient simulations could be implemented, e.g., by performing more frequent learning stages, storing larger training sets, or changing the optimization algorithm.

III Results

This section focuses on verifying the accuracy of the self-learning PQMC scheme described in the previous Section. The ferromagnetic quantum Ising chain (1) is used as a testbed. The PQMC predictions for the ground-state energy are compared against the exact value computed via the Jordan-Wigner transformation Cabrera and Jullien (1987). Notice that the comparison is made with the Jordan-Wigner results corresponding to the same system size and the same (open) boundary condition employed in the PQMC simulations. Let’s first focus on a chain with $N=80$ spins, in the absence of a longitudinal field (i.e., with $h=0$), with a transverse field of intensity $\Gamma=0.85$. This value is in the ferromagnetic phase, close to the quantum critical point, which corresponds to $\Gamma=1$ . This is the regime where the systematic bias due to the finite random-walker population is more sizable Inack et al. (2018a, b). [For moderately large $N_{w}$, the bias is larger around $\Gamma\approx 0.85$. However, in the $N_{w}\rightarrow\infty$ limit, the slowest vanishing of this bias occurs at the quantum critical point $\Gamma=1$, in the thermodynamic limit.]

Indeed, as shown in Fig. 2, a simple PQMC simulation performed without guiding wave function, i.e., with $\psi_{G}(\mathit{\boldsymbol{x}})=1$, and also one guided by a fixed RBM wave function $\psi_{G_{s=0}}(\mathit{\boldsymbol{x}})$ with $N_{h}=20$ hidden units and random parameters (see previous section for more details), provide energies with a noticeable upward bias of $\sim 1-2\%$. This bias occurs despite the relatively large population of random walkers $N_{w}=10^{4}$. Instead, the self-learning PQMC simulations guided by the adaptive wave function $\psi_{G_{s}}(\mathit{\boldsymbol{x}})$ display a radically different behavior. In the first stint ($s=0$), the guiding wave function $\psi_{G_{s=0}}(\mathit{\boldsymbol{x}})$ is the RBM with random parameters described in the previous section; the systematic bias is, again, noticeable. However, this bias completely disappears after a few stints/learning-stages $s\gtrsim 3$. Not only the bias disappears, but also the statistical fluctuations are drastically suppressed, leading to an extremely accurate and precise agreement with the Jordan-Wigner prediction. This indicates that the RBM is capable of accurately learning the random-walker distribution, leading to an adaptive guiding wave function $\psi_{G_{s\rightarrow\infty}}(\mathit{\boldsymbol{x}})$ sufficiently similar to the ground-state $\psi_{0}(\mathit{\boldsymbol{x}})$ to eliminate the systematic bias in the PQMC simulation. To quantify the accuracy of the learned adaptive wave functions $\psi_{G_{s}}(\mathit{\boldsymbol{x}})$, we compute the corresponding energy expectation values

[TABLE]

These computations are performed via Monte Carlo integration using a simple Metropolis sampling algorithm. The corresponding results are displayed in Fig. 2, placed at the final imaginary-time of the corresponding PQMC stint $\tau=\tau_{s}(s+1)$. For small $s$, due to the memory of the initial guiding wave function $\psi_{G_{s=0}}(\mathit{\boldsymbol{x}})$ and, to a lesser extent, to the bias affecting the PQMC algorithm, this expectation value largely deviates from the exact ground-state energy. This deviation might also originate from limitations in the learning process due to the finite number of learning updates $N_{\mathrm{steps}}$ and/or of the instances in the training set $N_{b}$. After a few stints/learning-stages the deviation from the ground-state energy significantly reduces, indicating that the adaptive RBM wave function better approximates the ground state. Remarkably, after a number $s\gg 1$ of stints and learning stages, $E_{\psi_{G_{s}}}$ approaches the optimal variational estimate $E_{\mathrm{var.min.}}$ (defined in the previous section) for an RBM wave function obtained using the NetKet library Carleo et al. (2019) for the same value of $N_{h}$. This library implements various sophisticated optimization algorithms. The results reported here are obtained using the stochastic reconfiguration method, which appears to be the most effective algorithm for the model addressed here. The agreement with the NetKet variational prediction indicates that the unsupervised learning algorithm is capable of identifying an optimal, or close to optimal, ground-state wave function from the random walkers sampled by the PQMC algorithm.

The adaptive RBM wave function $\psi_{G_{s}}(\mathit{\boldsymbol{x}})$ with $N_{h}=20$ is sufficiently accurate to eliminate the bias and to boost the efficiency of the PQMC simulation. However, the corresponding energy expectation value $E_{\psi_{G_{s}}}$ does not precisely agree with the exact ground-state energy, even after a number of stints $s\rightarrow\infty$.

As shown in Fig. 3, this residual deviation can be systematically reduced by increasing the number of hidden neurons $N_{h}$ in the RBM. This analysis is shown for the quantum critical point $\Gamma=1$, where it is more challenging to accurately approximate the ground-state with neural network ansatzes Carleo and Troyer (2017); Inack et al. (2018b); Beach et al. (2019); Collura et al. (2019). The system size is, again, $N=80$. One also observes that here the NetKet variational estimates $E_{\mathrm{var.min.}}$ for small $N_{h}$ are slightly lower than $E_{\psi_{G_{s\rightarrow\infty}}}$. This might indicate that the stochastic reconfiguration algorithm is able to better optimize the RBM parameters. However, one should also consider that when a variational ansatz is not very accurate, and the exact ground state is not contained in the manifold defined by the variational parameters, minimizing the energy expectation value does not necessarily provide a better overall description of the ground state. One cannot exclude that the adaptive RBM wave function $\psi_{G_{s\rightarrow\infty}}(\mathit{\boldsymbol{x}})$ produced by the unsupervised learning protocol, which aims at minimizing the Kullback-Leibler divergence with respect to the PQMC samples rather than at minimizing the energy expectation value, would provide more accurate predictions of physical properties other than the energy. As a reference, Fig. 3 displays also the minimal variational energy corresponding to the unrestricted Boltzmann machine ansatz of Ref. Inack et al. (2018b). This is obtained by minimizing the energy expectation value with an implementation of the stochastic gradient descent algorithm which samples the additional hidden spin variables. This variational estimate is comparable to the one corresponding to an RBM wave-function with $N_{h}=80-100$ hidden spins. More importantly, if the number of hidden units is too small, the adaptive RBM wave function might not be sufficiently accurate to eliminate the systematic bias of the PQMC simulation, even after many stints/learning-stages. This effect is visualized in Fig. 4, again for $\Gamma=1$. Indeed, one notices that, when the adaptive guiding function has $N_{h}=5-10$ hidden spins, the PQMC results are affected by a small bias $\sim 0.1\%$. This bias completely disappears for $N_{h}\geqslant 20$.

By exploiting the adaptive RBM wave function, unbiased PQMC results can be obtained also for larger system sizes. Figure 5 shows data for different system sizes (up to $N=160$). With a number of hidden spins $N_{h}=N/2$, one sees from the inset that the relative error with respect to the exact result is as small as $\sim 10^{-3}\%$ even for the largest system size.

The PQMC ground-state energy predictions are shown in Fig. 6 as a function of the transverse field $\Gamma$. This plot explores a moderately broad region around the critical point. The system size is $N=80$. Perfect agreement with the Jordan-Wigner prediction is found over the whole region. Results obtained in the presence of a longitudinal field of intensity $h=0.2$ are also shown. In this case, the comparison is made against the predictions of loop QMC simulations Todo and Kato (2001) performed with the ALPS library Bauer et al. (2011). Notice that the loop QMC algorithm provides finite-temperature results, but we consider here a sufficiently low temperature $T=0.05$ (we set $k_{B}=1$) such that thermal effects are negligible. In the PQMC simulations, the guiding function is, for both $h=0$ and $h=0.2$, an optimized adaptive RBM wave function $\psi_{G_{s\rightarrow\infty}}(\mathit{\boldsymbol{x}})$ with $N_{h}=40$ hidden spins. For $h=0$, the bias terms of the RBM are set to zero; in fact, we find that optimizing them does not lead to a more accurate guiding wave function. This is expected, since nonzero bias terms would break the $Z_{2}$ symmetry of the (finite system) ground state. Instead, for $h=0.2$ the unsupervised learning algorithm leads to sizable bias terms.

The PQMC algorithm provides predictions also for observables other than the energy. In Fig. 7, the estimates of the ground-state magnetization $M\equiv\left<\psi_{0}\right|\sum_{j}\sigma_{j}^{z}\left|\psi_{0}\right>$, obtained using the forward-walking technique, are shown as a function of the transverse field. A finite longitudinal field $h=0.2$ is included, so that the $Z_{2}$ symmetry, which would lead to a vanishing average magnetization in finite systems, is broken. The adaptive RBM wave function has $N_{h}=N/2=40$ hidden units. The agreement with the loop QMC predictions is extremely precise. This can be appreciated in the inset of Fig. 7, where the loop QMC results are plotted as a function of the temperature $T$, showing that thermal effects become indeed irrelevant at the lowest temperatures we consider. The inset displays also the prediction provided by the PQMC mixed estimator (see the definition in the previous section). This is affected by a remarkably small bias $\sim 0.1\%$, indicating that the optimized adaptive RBM guiding function is a good approximation for the ground state. As expected, the pure estimator removes even this small bias.

IV Conclusions

Machine learning techniques have been employed in the recent past to implement smart updates in classical and in path-integral Monte Carlo simulations Huang and Wang (2017); Liu et al. (2017a, b); Shen et al. (2018); Li et al. (2019), leading to a substantial efficiency improvement. In this paper, we used them to boost the performance of QMC simulations based on projective algorithms. Specifically, we have introduced a self-learning PQMC algorithm in which the guiding wave function is adaptively developed along the PQMC simulation, avoiding the need of a separate variational optimization. This adaptive guiding wave function is based on an RBM, and it is optimized via unsupervised machine learning using the $k$-step contrastive divergence algorithm. If the number of hidden units is sufficiently large, this adaptive guiding wave function allows one to eliminate any systematic bias in the PQMC simulation and to drastically boost its performance by suppressing the statistical fluctuations.

The self-learning PQMC algorithm appears to be particularly suitable when it is otherwise challenging to guess an accurate ansatz for the ground-state wave function. This is the case, e.g., of the disordered Ising Hamiltonians that describe the combinatorial optimization problems commonly addressed in the field of adiabatic quantum optimization. As a side product of the self-learning PQMC simulations, one obtains an accurate approximation for the ground-state wave function written in the form of an RBM. By construction, this function minimizes the Kullback-Leibler divergence with respect to the random-walker distribution. We argue that this unsupervised learning approach is complementary to other common techniques for wave function optimization based on the variational minimization of energy expectation values. Following Ref. Carleo and Troyer (2017), one might refer to these latter techniques as reinforcement learning. It is indeed plausible that, when the available ansatz is not sufficiently flexible to exactly describe the ground state, the minimization of the Kullback-Leibler divergence might lead to a fairer description of ground state properties other than the energy. Our unsupervised learning approach is also complementary to the optimization techniques for neural network states based on supervised machine learning introduced very recently in Ref. Kochkov and Clark (2018).

The self-learning PQMC algorithm has been tested in the ferromagnetic quantum Ising chain, obtaining excellent agreement with the Jordan-Wigner theory and with the loop QMC simulations performed in the low temperature limit. The model considered here is not affected by a sign problem. As a perspective, one can envision the extension of the self-learning PQMC technique to models where the sign problem occurs. Examples are many-fermion models and other non-stoquastic Hamiltonians relevant in adiabatic quantum computing. Fermionic PQMC simulations can be implemented using random walkers that carry a sign. They have been performed in the past adopting fixed-node and released node approaches, eventually performing annihilation of walkers with opposite sign Booth et al. (2009) to alleviate the suppression of the signal-to-noise ratio associated to the sign problem. In such fermionic simulations, the adaptive RBM guiding function would remain nonnegative, since its aim is to learn the probability distribution of the random walkers. It might allow one to further alleviate the sign problem by learning the nodal regions, and so reducing the random-walker crossings that lead to the suppression of the signal-to-noise ratio. We leave these endeavors to future investigations.

Acknowledgements

We acknowledge useful discussions with Ferran Mazzanti, Giacomo Torlai, and Giuseppe Santoro. E. M. I. thanks Dustin Lang for his help in performing the simulations on the Perimeter Institute HPC architectures.

S. P. and P. P. acknowledge financial support from the FAR2018 project of the University of Camerino and from the Italian MIUR under the project PRIN2017 CEnTraL 20172H2SC4. S. P. also acknowledges the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Economic Development, Job Creation and Trade. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.