Integrating Neural Networks with a Quantum Simulator for State Reconstruction

TL;DR

This paper presents a method combining neural networks with a quantum simulator to reconstruct quantum states from experimental data, effectively handling measurement errors and capturing complex observables.

Contribution

It introduces a neural network approach, specifically using RBMs, integrated with quantum hardware for accurate state reconstruction considering experimental errors.

Findings

Successfully reconstructed quantum states from Rydberg simulator data.

Captured complex observables like Rénnyi mutual information.

Mitigated measurement errors with a novel regularization technique.

Abstract

We demonstrate quantum many-body state reconstruction from experimental data generated by a programmable quantum simulator, by means of a neural network model incorporating known experimental errors. Specifically, we extract restricted Boltzmann machine (RBM) wavefunctions from data produced by a Rydberg quantum simulator with eight and nine atoms in a single measurement basis, and apply a novel regularization technique to mitigate the effects of measurement errors in the training data. Reconstructions of modest complexity are able to capture one- and two-body observables not accessible to experimentalists, as well as more sophisticated observables such as the R\'enyi mutual information. Our results open the door to integration of machine learning architectures with intermediate-scale quantum hardware.

Equations

$$\hat{H}(\Omega ,\Delta )=-\Delta \sum _{i=1}^N{\hat{n}}_i-\frac{\Omega }{2}\sum _{i=1}^N{\hat{\sigma }}_i^x+\sum _{i<j}\frac{V_{nn}}{|i-j{|}^6}{\hat{n}}_i{\hat{n}}_j,$$

$$|\psi ⟩=\frac{1}{\sqrt{2}}|rgrggrgr⟩+\frac{1}{2}|rgrgrggr⟩+\frac{1}{2}|rggrgrgr⟩.$$

$$|\psi ⟩=\sum _{\mathbf{\sigma }}\sqrt{p(\mathbf{\sigma })}|\mathbf{\sigma }⟩.$$

$$p_{\mathbf{\lambda }}(\mathbf{\sigma })$$

$$\hat{H}(\Omega ,\Delta )=-\Delta \sum _{i=1}^N{\hat{n}}_i-\frac{\Omega }{2}\sum _{i=1}^N{\hat{\sigma }}_i^x+\sum _{i<j}\frac{V_{nn}}{|i-j{|}^6}{\hat{n}}_i{\hat{n}}_j$$

$$|e_1⟩$$

$$|e_2⟩$$

$$|e_3⟩$$

$$⟨e_1|H_{\text{eff}}|e_2⟩=⟨e_1|H_{\text{eff}}|e_3⟩$$

$$p(\mathbf{\tau }|\mathbf{\sigma })={\delta }_{\mathbf{\tau },\mathbf{\sigma }}$$

$${\mathcal{L}}_{\mathbf{\lambda }}=-\frac{1}{|\mathcal{D}|}\sum _{\mathbf{\tau }\in \mathcal{D}}\log p_{\mathbf{\lambda }}(\mathbf{\tau })$$

$${\nabla }_{\mathbf{\lambda }}{\mathcal{L}}_{\mathbf{\lambda }}$$

$${\mathcal{E}}_{\text{eff}}(\mathbf{\sigma })=\mathbf{b}\cdot \mathbf{\sigma }+\sum _j\log \left(1+e^{W_{ji}{\sigma }_i+c_j}\right)$$

$$\frac{1}{|\mathcal{D}|}\sum _{\mathbf{\tau }\in D}{\nabla }_{\mathbf{\lambda }}{\mathcal{E}}_{\text{eff}}(\mathbf{\tau })\approx \frac{1}{M}\sum _{i=1}^M{\nabla }_{\mathbf{\lambda }}{\mathcal{E}}_{\text{eff}}({\mathbf{\tau }}_i)$$

$$⟨{\nabla }_{\mathbf{\lambda }}{\mathcal{E}}_{\text{eff}}(\mathbf{\sigma }){⟩}_{p_{\mathbf{\lambda }}(\mathbf{\sigma })}\approx \frac{1}{M}\sum _{i=1}^M{\nabla }_{\mathbf{\lambda }}{\mathcal{E}}_{\text{eff}}({\mathbf{\sigma }}_i)$$

$${\mathcal{L}}_{\mathbf{\lambda }}=-\frac{1}{|\mathcal{D}|}\sum _{\mathbf{\tau }\in \mathcal{D}}\log {\tilde{p}}_{\mathbf{\lambda }}(\mathbf{\tau })$$

$${\nabla }_{\mathbf{\lambda }}{\mathcal{L}}_{\mathbf{\lambda }}$$

$${\tilde{p}}_{\mathbf{\lambda }}(\mathbf{\sigma }|\mathbf{\tau })=\frac{p(\mathbf{\tau }|\mathbf{\sigma })p_{\mathbf{\lambda }}(\mathbf{\sigma })}{{\tilde{p}}_{\mathbf{\lambda }}(\mathbf{\tau })}$$

$$p(\mathbf{\tau }|\mathbf{\sigma })=\prod _{i}p(\mathbf{\tau }|{\sigma }_i)$$

$${\tilde{p}}_{\mathbf{\lambda }}(\mathbf{\sigma }|\mathbf{\tau },\mathbf{h})$$

$$=\prod _i{\tilde{p}}_{\mathbf{\lambda }}({\sigma }_i|\mathbf{\tau },\mathbf{h}),$$

$$p(\mathbf{\tau }|\mathbf{\sigma })$$

$$\tilde{W}$$

$${\tilde{b}}_{\sigma ,i}$$

$${\tilde{b}}_{\tau ,i}$$

$$⟨\hat{\mathcal{O}}{⟩}_{{\psi }_{\mathbf{\lambda }}}$$

$$=\sum _{\mathbf{\sigma }}|{\psi }_{\mathbf{\lambda }}(\mathbf{\sigma }){|}^2\sum _{{\mathbf{\sigma }}^{\prime }}⟨\mathbf{\sigma }|\hat{\mathcal{O}}|{\mathbf{\sigma }}^{\prime }⟩\frac{{\psi }_{\mathbf{\lambda }}({\mathbf{\sigma }}^{\prime })}{{\psi }_{\mathbf{\lambda }}(\mathbf{\sigma })}$$

$$:=⟨{\mathcal{O}}_L(\mathbf{\sigma }){⟩}_{p_{\mathbf{\lambda }}(\mathbf{\sigma })}$$

$$\simeq n_{\mathrm{mc}}^{-1}\sum _{k=1}^{n_{\mathrm{mc}}}{\mathcal{O}}_L({\mathbf{\sigma }}_k)$$

$$P_{\text{FD}}(\mathbf{\tau })=\frac{1}{N_s}\sum _{{\mathbf{\tau }}_i\in \mathcal{D}}{\delta }_{\mathbf{\tau },{\mathbf{\tau }}_i}$$

$$|\Psi ⟩=\sum _{\mathbf{\tau }}\sqrt{P_{\text{FD}}(\mathbf{\tau })}|\mathbf{\tau }⟩$$

$$F\left(P_{\text{FD}},P_{\text{GT}}\right)\le \sqrt{N_s}{e}^{-H_2/4}$$

$$\begin{array}{rl}\frac{d\hat{\rho }}{dt}& =-i[\hat{H}(\Omega (t),\Delta (t))+{\hat{H}}_{dis},\hat{\rho }]\\ & +\sum _{i=1}^N\sum _{t=rg,gg}{\gamma }_t\left({\tilde{\sigma }}_i^t\hat{\rho }{\tilde{\sigma }}_i^{t†}-\frac{1}{2}\left\{{\tilde{\sigma }}_i^{t†}{\tilde{\sigma }}_i^t,\hat{\rho }\right\}\right)\end{array}$$

$$F(\hat{\rho },\hat{\sigma })=\text{Tr}\left[\sqrt{\sqrt{\hat{\rho }}\hat{\sigma }\sqrt{\hat{\rho }}}\right]$$

$${\hat{\rho }}_i^{(s)}$$

Full text

Quantum state tomography Banaszek et al. (2013) is an important tool for reconstructing generic quantum states, but traditional techniques require a number of measurements scaling exponentially in the system size Häffner et al. (2005). In certain cases, methods that exploit particular entanglement or symmetry properties Cramer et al. (2010); Lee and Landon-Cardinal (2015); Riofrío et al. (2017); Lanyon et al. (2017); Tóth et al. (2010) allow for more efficient tomography of states prepared in experiments. However, such approaches still involve explicit reconstruction of local density operators Ježek et al. (2003); Cramer et al. (2010), incurring a significant computational overhead in the estimation of nontrivial observables from experimental data – especially in the presence of measurement errors introduced by realistic experimental hardware. In order to facilitate the characterization of near-term quantum hardware Preskill (2018), a state reconstruction method which can efficiently extract physical quantities of interest directly from noisy experimental datasets is highly desirable.

Neural network-based machine learning has recently emerged as a powerful technique for learning compact representations of high-dimensional data Hinton (2006); Graves et al. (2013); LeCun et al. (2015). In experimental quantum science, these tools have already been applied profitably to the classification of experimental snapshots Rem et al. (2019); Bohrdt et al. (2019) and qubit readout Seif et al. (2018). The same data-driven approach can be applied to tomographic tasks. Recent theoretical work has demonstrated that a generative model called a restricted Boltzmann machine (RBM) is capable of accurate reconstruction of quantum states and observables directly from synthetic datasets generated by numerical algorithms Torlai and Melko (2019).

In this Letter, we present a proof-of-principle demonstration of neural network quantum state reconstruction from experimental data. Our experimental system consists of a one-dimensional array of strongly interacting Rydberg atoms Bernien et al. (2017); Endres et al. (2016). Leveraging the high purity and approximate positivity of the experimental state, we train RBMs using single measurement basis data consisting of bit-strings obtained via repeated, simultaneous single-shot readout of the ground and Rydberg populations of all atoms. The RBMs learn a higher-fidelity and more efficient representation of the underlying bit-string probability distributions than standard inference from the limited size training dataset. This approach also enables us to implement an efficient procedure for denoising the full probability distribution from bit-flip-type measurement errors, by incorporating a dedicated “noise layer” in the network architecture. We test the validity of our approach by comparing predictions of the trained RBMs with numerical results for observables that are off-diagonal in the measurement basis, including the quantum mutual information. These results demonstrate the utility of RBMs in reconstructing approximately pure, positive states from experimental data, and pave the way to further integration of neural network models with quantum hardware.

Experimental system.

Our experimental approach Bernien et al. (2017); Endres et al. (2016) involves a programmable Rydberg atom quantum simulator, a flexible neutral-atom system for realizing Ising-type quantum spin models Schauß et al. (2015); Labuhn et al. (2016); Zeiher et al. (2016, 2017); Bernien et al. (2017); Guardado-Sanchez et al. (2018); Barredo et al. (2018). In the present experiments (Fig. 1a), a one-dimensional array of $N$ trapped Rubidium atoms is prepared; $N=8$ atoms are used below, but we have also applied our protocol to arrays of $N=9$ atoms SI . Each atom can occupy a ground state $|g\rangle$ or an excited (Rydberg) state $|r\rangle$, and two atoms excited to the Rydberg state at a distance $r$ interact with a van der Waals-type potential $V(r)\propto r^{-6}$. When subjected to a uniform laser drive, the effective Hamiltonian of the many-body system can be written as Bernien et al. (2017); Pohl et al. (2010); van Bijnen et al. (2011); Schauß et al. (2015)

[TABLE]

where $V_{nn}$ is the interaction strength between Rydberg atoms at adjacent sites, $\hat{\sigma}_{i}^{\alpha}$, with $\alpha=x,y,z$, are the Pauli pseudo-spin operators at site $i$ (defined as $\hat{\sigma}_{i}^{z}=|r_{i}\rangle\langle r_{i}|-|g_{i}\rangle\langle g_{i}|$, $\hat{\sigma}_{i}^{x}=|r_{i}\rangle\langle g_{i}|+h.c.$, etc), and $\hat{n}_{i}=\frac{1}{2}\left(1+\hat{\sigma}_{i}^{z}\right)$ projects onto the Rydberg state at site $i$. The parameters $\Omega,\Delta$ denote the effective Rabi frequency and detuning, respectively, which characterize the laser drive, and can be varied in time as $\Omega(t),\Delta(t)$ to drive the system into nontrivial ordered phases Pohl et al. (2010); Bernien et al. (2017); Weimer and Büchler (2010); Sela et al. (2011).

We focus on the transition into the $\mathbb{Z}_{2}$ phase Bernien et al. (2017), where a high density of Rydberg excitations is energetically favorable, subject to the constraint that no two adjacent atoms are excited. The atoms are initially pumped into the fiducial state $|g\,g\,g\,g\,g\dots\rangle$, coinciding with the ground state of Hamiltonian (5) at $t=0$. They then evolve adiabatically under a “sweep” of the laser parameters $\Omega(t),\Delta(t)$ for a time $T_{\textnormal{ev}}$, with $\hat{H}(\Omega(T_{\textnormal{ev}}),\Delta(T_{\textnormal{ev}}))$ lying deep in the $\mathbb{Z}_{2}$ phase (Fig. 1c). For our eight-atom system, the final $\mathbb{Z}_{2}$-ordered state at $t=T_{\textnormal{ev}}$ is well approximated by the ground state of the Rydberg Hamiltonian with a small transverse field and short-range interactions only SI :

[TABLE]

This state exhibits quantum fluctuations on two pairs of adjacent atoms, as indicated in Fig. 1b.

Pure state ansatz.

The ground state of the Hamiltonian (5) has real, positive amplitudes in the occupation number basis $|\bm{\sigma}\rangle=|\sigma_{1},\dots,\sigma_{N}\rangle$ – defined as the simultaneous eigenstates of $\hat{n}_{1},\dots,\hat{n}_{N}$ – as long as $\Omega>0$ Bravyi et al. (2008), which can always be arranged by applying a suitable global unitary 111The exact phases required to render the Rydberg Hamiltonian in the form (5) will vary between different experimental realizations, but as the final measurements are always taken in the occupation number basis this has no effect on observables, provided no variation in the laser phase occurs during evolution.. Therefore, if the quantum state of the simulator evolves perfectly adiabatically and with negligible loss of purity, it is uniquely characterized by its probability distribution $p(\bm{\sigma})$ over projective measurements in the $|\bm{\sigma}\rangle$ basis, and at any time may be written as the pure state

[TABLE]

Of course, some loss of purity is inevitable – in our experiments, due primarily to single-atom decay and dephasing processes Levine et al. (2018) – and the true state is described by a mixed density operator $\hat{\rho}$. Although this pure state approximation cannot capture all of the physics of the experimental state, it can in principle accurately describe local subsystems, to the extent that the corresponding reduced density operators of the true and reconstructed states agree SI . We adopt the pure, positive state ansatz in all of our reconstruction efforts below.

Neural network model.

While the quantum state (3) can in principle be inferred directly from a set of raw measurements (i.e. by inverting the measurement counts of each configuration to estimate $p(\bm{\sigma})$), such an approach is limited to small systems and very large datasets. In contrast, generative models used in unsupervised machine learning tasks can capture the structure of the distribution $p(\bm{\sigma})$, generalizing beyond a limited set of training samples. This results in a higher-fidelity reconstruction and a model size scaling polynomially in the system size (Fig. 2). Moreover, using a generative model rather than direct inference from the data enables automatic correction of this distribution for known measurement errors using a “noise layer” (see Fig. 1a and description below).

We parametrize $p(\bm{\sigma})$ with a generative model known as an RBM Ackley et al. (1985); Smolensky (1986), a stochastic neural network with two layers of binary units. The “visible” layer $\bm{\sigma}$ describes the atomic states of the Rydberg chain in the occupation number basis, while a hidden layer $\bm{h}$ captures correlations between visible units. The RBM defines the following probability distribution for the visible layer:

[TABLE]

where $Z_{\bm{\lambda}}$ is a normalization constant, and the real-valued network parameters are $\bm{\lambda}=\{\bm{W},\bm{b},\bm{c}\}$, with $\bm{W}$ being the weights connecting the two layers and $\bm{b}$ ($\bm{c}$) the visible (hidden) bias vectors. We use the visible layer of the RBM to define the projective measurement distribution $p(\bm{\sigma})$ of the pure state (3), resulting in an RBM wavefunction with positive amplitudes Torlai and Melko (2016): $\psi_{\bm{\lambda}}(\bm{\sigma})=\langle\bm{\sigma}|\psi_{\bm{\lambda}}\rangle=\sqrt{p_{\boldsymbol{\lambda}}(\bm{\sigma})}$. We have numerically verified that this RBM wavefunction can accurately describe states relevant to our experiment, with a number of parameters scaling only quadratically in system size (Fig. 2 and SI , Sec. IV), in accordance with recent scaling studies for quantum Ising ground states Sehayek et al. (2019). We point out that, although pure states with nontrivial phases (Torlai et al., 2018; Carleo and Troyer, 2017), as well as mixed state models (Torlai and Melko, 2018; Carrasquilla et al., 2019), could be applied using similar neural network models, measurements in other bases would be required.

Measurement process and noise layer.

Measurement data consists of a collection of $N$-bit strings $\bm{\tau}=(\tau_{1},\dots,\tau_{N}$), with $\tau_{j}=0,1$ indicating that atom $j$ was recorded as being in the ground $|g\rangle$ or Rydberg state $|r\rangle$ respectively Bernien et al. (2017). Such measurements are never perfect, and there are small measurement error probabilities $p(1|0)\sim 1\%$, $p(0|1)\sim 4\%$ Levine et al. (2018) for an atom in the ground state to be recorded as excited and vice-versa. These result in experimental data $\bm{\tau}$ that do not correspond to projective measurements. Instead, the measurement process can be described as a positive-operator valued measure (POVM) Nielsen and Chuang (2010) with measurement operators $\hat{\Pi}_{\bm{\tau}}=\sum_{\bm{\sigma}}p(\bm{\tau}|\bm{\sigma})|\bm{\sigma}\rangle\langle\bm{\sigma}|$, where $p(\bm{\tau}|\bm{\sigma})=\prod_{j=1}^{N}p(\tau_{j}|\sigma_{j})$ is the probability of the experimentalist recording $\bm{\tau}$ if the atoms are prepared in the state $|\bm{\sigma}\rangle$. The probability distribution sampled in the experiment is then $P_{\textnormal{exp}}(\bm{\tau})=\textnormal{Tr}\left[\hat{\rho}\hat{\Pi}_{\bm{\tau}}\right]$.

The experimental measurement process is incorporated into our model via a third binary layer, the so-called noise layer (Fig. 1a), which represents the observed POVM outcomes $\bm{\tau}$. The measurement error rates $p(\bm{\tau}|\bm{\sigma})$ are included as connections between the visible and noise layers Yichuan Tang et al. (2012), by assigning a probability $\tilde{p}_{\boldsymbol{\lambda}}(\bm{\bm{\tau}})=\sum_{\bm{\sigma}}p(\bm{\tau}|\bm{\sigma})p_{\boldsymbol{\lambda}}(\bm{\sigma})$ to the measurement result $\bm{\tau}$.

The full three-layer network is trained to learn parameters $\bm{\lambda}$ which maximize the log-likelihood of the recorded POVM outcomes under $\tilde{p}_{\boldsymbol{\lambda}}(\bm{\tau})$. During training, the noise layer prevents the parameters ${\bm{\lambda}}$ from fitting to spurious features in the data produced by measurement errors. This noise layer regularization significantly improves the fidelity between $|\psi_{\bm{\lambda}}\rangle$ and the state $\hat{\rho}$ underlying the data; numerical tests (Fig. 2) based on Lindbladian simulation of our experiment result in fidelities greater than 90% for the full many-body state at the end of the sweep, even when decoherence processes are included. All reconstructions presented below are obtained in this fashion.

Experimental reconstruction.

In the experiment, at fifteen subsequent time-steps $t$ (Fig. 1c), the sweep is halted and measurements $\bm{\tau}$ are sampled from the state $\hat{\rho}(t)$. At each time-step, a dataset of around 3,000 samples is collected and used to train a three-layer model with $2N=16$ hidden units. After training the networks, standard sampling methods can be applied to compute expectation values of observables, with a computational cost scaling polynomially in the network size SI . We consider in particular the connected correlation functions $\langle\hat{\sigma}_{i}^{\alpha}\hat{\sigma}_{j}^{\alpha}\rangle_{c}=\langle\hat{\sigma}_{i}^{\alpha}\hat{\sigma}_{j}^{\alpha}\rangle-\langle\hat{\sigma}_{i}^{\alpha}\rangle\langle\hat{\sigma}_{j}^{\alpha}\rangle$ for $\alpha=x,y,z$, and their spatial averages, $\bar{g}^{\alpha\alpha}(s)=\frac{1}{N-s}\sum_{i=1}^{N-s}\langle\hat{\sigma}_{i}^{\alpha}\hat{\sigma}_{i+s}^{\alpha}\rangle_{c}$.

In Fig. 3a-b, we verify that our reconstructions learn to represent their training sets, by examining their ability to accurately reproduce observables which are diagonal in the occupation number basis. The networks learn the strong two-body correlations $\langle\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}\rangle_{c}$ present in the experimental data. We compare the results of the reconstruction process to the exact solutions of a Lindblad master equation for the full many-body evolution. Our Lindbladian simulation predicts Rydberg excitation profiles in excellent agreement with experiment, but its significantly weaker correlations suggest our model for the sweep dynamics is partially incomplete.

Turning to experimentally inaccessible quantities (Fig. 3c-d), the reconstructions and simulation agree qualitatively in the temporal and spatial profiles of the transverse field $\langle\hat{\sigma}_{i}^{x}\rangle$ and its two-point correlation function, although the RBMs predict somewhat larger values in the ordered phase. Note that the distinct spatial variation of the transverse field correlations, a signature of quantum fluctuations captured in the approximate state (2), is reconstructed directly from our experimental data. Training on synthetic data SI indicates that a large portion of the disagreement between reconstruction and simulation is due to the discrepancy between our Lindbladian model and experiment evident in Fig. 3a-b, not the RBM model itself.

Beyond few-body observables, an important question is whether entanglement properties are reproduced accurately in reconstruction. From our RBMs, the Rényi entropy – which requires specialized or hardware-specific protocols to access directly in experiment Islam et al. (2015); Brydges et al. (2018) – may be extracted in a scalable fashion by applying a state-replication and swap procedure virtually Hastings et al. (2010); Torlai et al. (2018). In fact, for pure experimental states, positive-pure ansatzes such as the RBM wavefunction provide a lower bound on the mutual information defined by the Rényi entropy (Zhang et al. (2011), Grover and Fisher (2015), see also SI , Sec. VIII), regardless of the sign structure of the true state. We demonstrate a reconstruction of the mutual information defined by the Rényi entropy in Fig. 4, finding that the RBM values are in remarkable agreement with the results of numerical simulation. Reconstructions on experimental states of $N=9$ capture a similar buildup in the mutual information during the sweep predicted by Lindbladian simulation SI .

Conclusions.

In this Letter, we have demonstrated neural-network reconstruction of experimental quantum states from data produced by a programmable Rydberg-atom quantum simulator. By leveraging the real-positive nature of the ground state wavefunction expected from the effective Hamiltonian, we trained restricted Boltzmann machines on measurements in the occupation basis only. An additional noise layer was added to the standard RBM architecture to mitigate measurement errors. Once trained, the RBM was queried to produce a variety of observables not accessible in the original experimental setup, including the Rényi entropy - a basis independent measure of the quantum entanglement of the wavefunction.

Our approach can be integrated without alteration into existing platforms where a positive wavefunction ansatz is a valid approximation, such as Bose-Hubbard experiments and some non-frustrated quantum spin simulators Bakr et al. (2009); Weitenberg et al. (2011); Kaufman et al. (2016); Labuhn et al. (2016). Access to multiple measurement bases would allow enhanced certification of the reconstruction, by providing direct experimental access to observables which are informationally complete for local subsystems. Also, with access to multiple bases the RBM protocol can be easily adapted to reconstruct non-positive and complex wavefunctions Torlai et al. (2018). Identifying the minimal set of measurement bases and the optimal protocol to collect the statistics represents a crucial step towards reconstruction of quantum states prepared by fermionic quantum simulators and non-equilibrium dynamics Cheuk et al. (2015); Greif et al. (2016).

In conclusion, machine learning techniques offer a means of increasing the amount of useful information that can be extracted from experiments, especially when hardware constrains the quantity or quality of accessible measurements. They can be used to offload the burden of technically expensive – or fundamentally impossible – measurements from experimental platforms in a noise-resilient fashion. We expect experimentalists will profit from deeper integration of machine learning architectures with quantum devices.

I Approximate Eight-atom ground state

The full Rydberg Hamiltonian is

[TABLE]

At the end of the experimental sweep, the Hamiltonian has a positive detuning and a small transverse field: $\Delta>0,V_{nn}\gg\Delta\gg|\Omega|$; furthermore, interactions between sites separated by more than two lattice spacings may be neglected, as they are weak compared to the frequencies which characterize the sweep profile. In this regime the four-excitation states

[TABLE]

are degenerate under the classical part of the Hamiltonian $-\Delta\sum_{i=1}^{N}\hat{n}_{i}+\sum_{i<j}\frac{V_{nn}}{|i-j|^{6}}\hat{n}_{i}\hat{n}_{j}$. The ground state lies in the subspace spanned by these three states: adding or removing an excitation requires an energy penalty proportional to $V_{nn}$ or $\Delta$ respectively. This degeneracy is lifted by a nonzero transverse field, which couples the blockaded states at second order in $\Omega$ through the three-excitation subspace. Using perturbation theory, an effective Hamiltonian Sachdev (2011) $H_{\textnormal{eff}}$ may be constructed for the blockaded subspace, whose nonzero matrix elements are given by

[TABLE]

The corresponding ground state is $|\Psi\rangle=\frac{1}{\sqrt{2}}|e_{1}\rangle+\frac{1}{2}\left(|e_{2}\rangle+|e_{3}\rangle\right)$.

II Reconstruction Methods

II.1 Note on terminology

Below we discuss strategies for training on experimental data which has been corrupted by a fixed, known noise process. $\bm{\sigma}$ will denote the variables prior to corruption by measurement errors, while $\bm{\tau}$ will denote those which have been subjected to the noise channel – that is, for a fixed true value $\bm{\sigma}$, the noisy outputs are distributed according to $p(\bm{\tau}|\bm{\sigma})$. In our experiment, $\bm{\tau}$ are the only accessible variables, which yield the bitstrings recorded in each dataset. A model with parameters $\bm{\lambda}$ specifies a distribution $p_{\boldsymbol{\lambda}}(\bm{\sigma})$ over the uncorrupted variables $\bm{\sigma}$, and a corresponding corrupted distribution $\tilde{p}_{\boldsymbol{\lambda}}(\bm{\tau})=\sum_{\bm{\sigma}}p(\bm{\tau}|\bm{\sigma})p_{\boldsymbol{\lambda}}(\bm{\sigma})$.

II.2 Standard RBM training method

The standard training method involves fitting the RBM distribution $p_{\boldsymbol{\lambda}}(\bm{\sigma})=\frac{1}{Z_{\bm{\lambda}}}\sum_{\bm{h}}e^{\>\bm{h}^{\top}\bm{W}\bm{\sigma}+\bm{b}\cdot\bm{\sigma}+\bm{c}\cdot\bm{h}}$ directly to the experimental datasets; in other words, it assumes a noise-free source of data:

[TABLE]

The optimal parameters $\bm{\lambda}=\{\mathbf{W},\mathbf{b},\mathbf{c}\}$ for which the RBM best reproduces the measurement data are found by minimizing the negative log-likelihood

[TABLE]

of the RBM distribution $p_{\boldsymbol{\lambda}}$ averaged over the dataset $\mathcal{D}$ ($|\mathcal{D}|$ denotes the size of the dataset). The gradient of the log-likelihood cost function with respect to the trainable parameters $\bm{\lambda}$ may be written

[TABLE]

where $\langle\cdot\rangle_{p_{\boldsymbol{\lambda}}(\bm{\sigma})}$ denotes the expectation value with respect to the distribution $p_{\boldsymbol{\lambda}}(\bm{\sigma})$, and the effective energies

[TABLE]

are defined by $p_{\boldsymbol{\lambda}}(\bm{\sigma})=\frac{1}{Z_{\bm{\lambda}}}e^{\mathcal{E}_{\textnormal{eff}}(\bm{\sigma})}$.

The second term in the cost function gradient (12) is estimated using a batch of samples $\bm{\tau}_{i}$ of size $M$ drawn from the training set $\mathcal{D}$:

[TABLE]

Exact computation of the expectation value with respect to $p_{\boldsymbol{\lambda}}(\bm{\sigma})$ requires summing over a number of configurations which is exponential in the system size, and is therefore not tractable. It can also be approximated by drawing $M$ samples $\bm{\sigma}_{i}$ distributed according to $p_{\boldsymbol{\lambda}}(\bm{\sigma})$ and using the estimator

[TABLE]

In principle, samples which obey the model distribution $p_{\boldsymbol{\lambda}}(\bm{\sigma})$ can be generated by block Gibbs sampling Hinton (2012), which involves repeatedly sampling from the conditional distributions $p_{\boldsymbol{\lambda}}(\bm{\sigma}|\bm{h})$ and $p_{\boldsymbol{\lambda}}(\bm{h}|\bm{\sigma})$. Because of the restricted nature of the RBM graph – there are no intra-layer connections – the conditional distributions factorize and each unit in a given layer can be exactly sampled simultaneously. In pseudocode, starting from a ‘seed’ visible state $\bm{\sigma}_{1}$, the Gibbs sampling algorithm is:

for  $i$ in $[1,2,...,k]$ do

     Sample $\bm{h}_{i}$ from $p_{\boldsymbol{\lambda}}(\bm{h}|\bm{\sigma}_{i})$

     Sample $\bm{\sigma}_{i+1}$ from $p_{\boldsymbol{\lambda}}(\bm{\sigma}|\bm{h}_{i})$

end for

– the output of the algorithm is the visible state $\bm{\sigma}_{k+1}$, which will obey the model distribution $p_{\boldsymbol{\lambda}}(\bm{\sigma})$ for a sufficiently large number of sampling steps. In practice, the contrastive divergence algorithm Hinton (2002) is applied, where the visible state is seeded with samples from the training set, and only a small number of sampling steps $k$ is used. In practice, moderate values $k\sim 10$ are sufficient for training with stochastic gradient descent. Additional information about the RBM and its training can be found in Ref. Torlai and Melko (2019). An open-source software library for RBM reconstruction of generic wavefunctions is also available Beach et al. (2019).

II.3 Noise-regularized training method

In the case where the training set is known to be corrupted by a noise process $p(\bm{\tau}|\bm{\sigma})$, our goal is to learn a model $p_{\boldsymbol{\lambda}}(\bm{\sigma})$ whose corresponding noise-corrupted distribution $\tilde{p}_{\boldsymbol{\lambda}}(\bm{\tau})$ fits the observed data. We therefore define the corresponding log-likelihood cost function

[TABLE]

and train the network to minimize it on each dataset. The cost gradient takes a form nearly identical to that of the standard training method (12),

[TABLE]

The second term in the gradient update step is now computed not directly from the training set samples $\bm{\tau}\in\mathcal{D}$, but rather from the Bayesian posterior distribution

[TABLE]

which the RBM assigns to visible states $\bm{\sigma}$, given an observation $\bm{\tau}$ in the noisy training set.

This alteration to the cost gradient may be viewed as a regularization of the training based on prior knowledge of the sampling process. Regularization in machine learning generally refers to techniques for improving the generalization performance of a model trained on a particular data set to new datasets drawn from the ‘ground truth’ source. A typical regularization scheme like weight decay does not specify a priori how the in-sample and true distributions differ, and therefore typically requires some sort of validation process – testing the model on held-out data – to select good hyperparameters. In contrast, our regularization method is applied in a context where all accessible datasets are corrupted by the same noise process. This makes validation as a means of selecting regularization hyperparameters impossible – but if the noise process is known, this is no obstacle as there are no free hyperparameters to select.

In applying equation (17) to the unsupervised training of an RBM, both contributions to the gradient now require computation of expectation values over marginalized distributions $p_{\boldsymbol{\lambda}}(\bm{\sigma})$, $\tilde{p}_{\boldsymbol{\lambda}}(\bm{\sigma}|\bm{\tau})$ of the RBM, and are therefore intractable to compute exactly. As in the noise-free training case, this problem may be circumvented using the contrastive divergence method: the first term is estimated by repeated sampling from the conditional distributions $p_{\boldsymbol{\lambda}}(\bm{\sigma}|\bm{h}),p_{\boldsymbol{\lambda}}(\bm{h}|\bm{\sigma})$, while the second uses the same alternating sampling from the ‘data-clamped’ distributions $\tilde{p}_{\boldsymbol{\lambda}}(\bm{\sigma}|\bm{h},\bm{\tau}),\tilde{p}_{\boldsymbol{\lambda}}(\bm{h}|\bm{\sigma},\bm{\tau})=p_{\boldsymbol{\lambda}}(\bm{h}|\bm{\sigma})$. As noted above, $p_{\boldsymbol{\lambda}}(\bm{\sigma}|\bm{h}),p_{\boldsymbol{\lambda}}(\bm{h}|\bm{\sigma})$ are both efficiently computable due to the restricted structure of the RBM layers $\bm{\sigma},\bm{h}$. Similarly, $\tilde{p}_{\boldsymbol{\lambda}}(\bm{\sigma}|\bm{h},\bm{\tau})$ is efficiently computable if the error probabilities satisfy a weaker condition, namely factorizing over the uncorrupted variables:

[TABLE]

In this case, the clamped distribution may be computed explicitly as

[TABLE]

amenable to efficient block-Gibbs sampling.

Fig. 5 provides an intuitive way to understand the noise regularization – the corrupted variables $\bm{\tau}$ may be included as a third noise layer appended to the standard, two-layer RBM graph, with conditional probabilities depending on the $\bm{\sigma}$ layer only. These can be interpreted as effective biases for the noise layer, which depend on the uncorrupted variables – for example, the independent bit-flip errors used to model our Rydberg experiment may be written as

[TABLE]

For brevity, we will sometimes refer to RBMs trained with this regularization as ‘three-layer’ machines, as opposed to their ‘two-layer’ counterparts trained in the standard fashion. Similar graphical models known as Deep Belief Nets Hinton et al. (2006) have previously been used for unsupervised learning tasks, but with a different, layer-wise training algorithm that does not incorporate prior information; a gated RBM architecture similar to the three-layer machine has also been applied to Gaussian noise models in occluded images Yichuan Tang et al. (2012).

II.4 Sampling from trained RBMs

After an RBM has been trained, new configurations of the uncorrupted variables $\{\bm{\sigma}\}$ can be drawn from the distribution $p_{\boldsymbol{\lambda}}(\bm{\sigma})$ using the block-Gibbs sampling techniques discussed above. The expectation value of a generic observable $\hat{\mathcal{O}}$ in the state $\psi_{\bm{\lambda}}(\bm{\sigma})=\sqrt{p_{\boldsymbol{\lambda}}(\bm{\sigma})}$ can then be approximated with a Monte Carlo average over $n_{\mathrm{mc}}$ samples:

[TABLE]

where the “local estimate” of the observable is defined to be $\mathcal{O}_{L}(\bm{\sigma})=\sum_{\bm{\sigma}^{\prime}}\langle\bm{\sigma}|\hat{\mathcal{O}}|\bm{\sigma}^{\prime}\rangle\frac{\psi_{\bm{\lambda}}(\bm{\sigma}^{\prime})}{\psi_{\bm{\lambda}}(\bm{\sigma})}$. In the case of nontrivial noise processes, to sample from the corrupted distribution $\tilde{p}_{\boldsymbol{\lambda}}(\bm{\tau})$ one may first generate an uncorrupted batch $\{\bm{\sigma}\}$ of data and then sample once from the conditional distribution $p(\bm{\tau}|\bm{\sigma})$ for each uncorrupted configuration.

For an RBM with $N$ visible and $N_{h}$ hidden units, the times for training and Monte Carlo observable estimation scale as $\mathcal{O}(NN_{h})$, or in terms of the model complexity $\alpha=N_{h}/N$, as $\mathcal{O}(\alpha N^{2})$; note that the number of visible units is fixed by the system size. A universal approximation theorem Le Roux and Bengio (2008) guarantees that RBMs can represent any distribution over binary variables, although an exponentially large number of hidden units may be required in general. However, many quantum states relevant to experiment, such as ground states of paradigmatic Hamiltonians and some matrix product states, have been found to admit efficient descriptions Torlai et al. (2018); Torlai and Melko (2016); Deng et al. (2017); Glasser et al. (2018); Chen et al. (2018). In the present work with eight atoms, the Hilbert space is small enough that all amplitudes and expectation values can be computed exactly, providing a valuable check on our procedure. Such a benchmark quickly becomes impossible with current classical hardware when the number of atoms approaches $\sim 20$ for pure states, and at even smaller chain lengths for the exact evaluation of non-pure states.

III Training details

III.1 Methods

The reconstructions presented in this work were trained using the three-layer scheme detailed above on experimental datasets of $N\approx 3000$ samples each. Training was performed using stochastic gradient descent with a decayed learning rate, the gradients being estimated via contrastive divergence with $k=30$ sampling steps. Since the visible layers of our machines are relatively small, exact computation of the negative-log-likelihood was possible on each set. Hyperparameters for training were therefore selected by cross-validation on a randomly chosen experimental set; the same hyperparameters were used in training on all datasets. The reconstructions presented in the text were trained on the full datasets; RBMs were also trained on 90/10 splits of each dataset in order to verify that the out-of-sample negative log likelihood did not grow during training. Error bars on reconstructed observables were computed from their variation across these training subsets in the final epochs of training. We found it beneficial to train each machine with the error rates set to zero for the first epoch.

To check that the networks learned a consistent representation of the experimental data, we performed a scaling analysis of the number of hidden units $N_{h}$ of the RBM when training on experimental data. Increasing the number of hidden units, we found convergence of the observables and log-likelihood for $N_{h}\sim N$ (see Fig. 6 for examples). The reconstructions presented in this work used RBMs with $N_{h}=2N=16$.

III.2 Training on larger systems

As a test of the robustness of our reconstruction procedure, we also trained RBMs on a second set of Rydberg atom data sampled from a larger chain of $N=9$ atoms. The dynamics of this system is governed by a master equation identical in structure to that used for modeling the eight-atom data presented in the main text, but with slightly different detuning and Rabi frequency profiles, and different effective decoherence rates.

Fig. 7 compares the results of this reconstruction to predictions of the relevant Lindbladian model, as well as experimental values where appropriate. Without alteration of the training procedure, the RBMs reconstruct quantum dynamics, as manifested in the transverse field and mutual information, in good agreement with Lindbladian predictions. This is a key benefit conferred to the experimentalist by the RBM reconstruction method. Indeed, given previous knowledge regarding the properties of the quantum state prepared in the experiment, RBM reconstruction of experimentally inaccessible observables allows for rapid and inexpensive detection of errors in state preparation and manipulation.

IV Generalization capabilities

A generative model is of little use if it merely mimics the statistics of the training set. Successful machine learning applications are built upon the ability to generalize from a given dataset, extracting representations of the data that capture relevant features of the ground truth distribution from which it was sampled. This requires some structure in the data for the machine to learn, and the extent to which it succeeds in doing so depends upon the architecture of the machine as well as the size of the dataset.

For relatively small datasets such as those used in this work, it is natural to wonder whether the apparatus of machine learning is necessary at all. In particular, given access to the frequency distribution (FD)

[TABLE]

defined by a particular dataset $\mathcal{D}$ consisting of $N_{s}$ samples, one may define a naive frequency distribution reconstruction of a pure state corresponding to the data, which simply memorizes the training set:

[TABLE]

The FD state model can be computed and stored in a time linear in the size of the dataset, by building a lookup table that associates each observed bitstring $\bm{\tau}$ with its empirical probability in the dataset, and assigning probability zero to all other bitstrings. Such a model may then be used to produce Monte-Carlo estimates of desired observables, in the same fashion as for RBM states.

In general, the FD reconstruction approach cannot scale to high-entropy distributions – if $H_{2}$ is the second-order Renyi entropy of the ground-truth distribution $P_{\textnormal{GT}}(\bm{\tau})$, the fidelity $F\left(P_{\textnormal{FD}},P_{\textnormal{GT}}\right)=\sum_{\bm{\tau}}\sqrt{P_{\textnormal{FD}}(\bm{\tau})P_{\textnormal{GT}}(\bm{\tau})}$ between the frequency distribution and the ground truth obeys the inequality

[TABLE]

– for a proof, see Section (VII). In particular, if the measurement-basis entropy is proportional to the system size – as is the case in even some very simple states, such as a product state of spins not aligned with the measurement basis – the frequency-distribution fidelity will decay exponentially in system size. The ability to extract a modest number of physically relevant features is therefore essential for accurate state reconstruction from generic datasets of realistic size. However, our eight-atom system is small enough compared to the size of the datasets ($\sqrt{N_{s}}\sim 2^{L}$) that the FD approach is not a priori infeasible.

To quantify the performance of RBM and FD reconstructions in the small-system regime, we sampled synthetic datasets (in the occupation number basis) of size $N_{s}$ up to $10^{5}$ from ground states of the Rydberg Hamiltonian in equation (5), for a selection of system sizes up to $N=16$ atoms. Ground state wavefunctions were computed using the QuSpin exact diagonalization package Weinberg and Bukov (2017); the Hamiltonian parameters were constant throughout and chosen to place the system near the phase transition into $\mathbb{Z}_{2}$ state: $V_{nn}=30$MHz, $\Omega=2$MHz, $\Delta\approx 1$MHz. For each dataset, we computed the fidelity $F_{\textnormal{FD}}$ of the frequency distribution state onto the ground-truth Rydberg wavefunction; an RBM with $N_{h}=N$ hidden units was then trained on the same dataset, and its fidelity $F_{\textnormal{RBM}}$ onto the true state was also recorded. The RBMs were all trained with the hyperparameters described in section III, but with $k=10$ contrastive divergence steps. Fig. 8 plots the resulting fidelities achieved by both reconstructions as a function of system size – RBMs of fixed complexity achieve significantly higher fidelities for large systems, with small improvements even at $N=8$.

Another issue of practical relevance is model size: given a dataset $\mathcal{D}$ of a particular size $N_{s}$, how many parameters are required to store each trained model? For the RBM, the number of (real-valued) parameters required to specify the model completely is determined by the size of the bias vectors and weight matrix, $N\cdot N_{h}+N+N_{h}$, and therefore quadratic in the system size for a fixed model complexity $N_{h}/N$. For the FD model, the number of parameters is determined by the size of the lookup table, i.e. the number of unique samples present in the dataset, and therefore bounded above by the dataset size $N_{s}$.

In Fig. 9, the model sizes of the RBM and FD reconstructions from Fig. 8 are compared as a function of system size; for $N\gtrsim 8$ atoms the RBMs are a significantly more efficient (not to mention more accurate) description of the quantum state.

Finally, we note that even for small systems, generative models provide an additional advantage in state reconstruction from noisy data: in the presence of measurement errors, the FD model is not representative of the ground truth for any dataset size, and simply inverting the conditional probabilities will generally result in unphysical prior distributions. Denoising methods for cleaning noisy binary datasets prior to reconstruction Greig et al. (1989); Geman and Geman (1984); Batson and Royer (2019) may be applied, but a model training step is still required.

V Effects of Decoherence

For pure state reconstruction to be useful in near-term quantum simulators, realistic decoherence processes must be accounted for. Here, we provide a brief description of the Lindbladian master equation used in our modeling of the experiment, and discuss means of assessing the quality of pure state reconstructions in the presence of decoherence.

V.1 Master equation for the Rydberg machine

To account for decoherence processes quantitatively, we have used a Lindblad model, described in detail in Ref. Levine et al. (2018), which includes two jump operators $\tilde{\sigma}_{i}^{rg}=|g\rangle\langle r|,\tilde{\sigma}_{i}^{gg}=|g\rangle\langle g|$ to represent decay and dephasing processes acting on atom $i$. The time evolution of the full state is given by the master equation

[TABLE]

where $\hat{H}_{dis}=-\sum_{i=1}^{N}\delta_{i}\hat{n}_{i}$ is the static disorder Hamiltonian containing the Doppler shifts $\delta_{i}$, and $\gamma_{t}$, $t=rg,gg$ are decoherence rates estimated from single-atom measurements Levine et al. (2018) as $1/\gamma_{rg}=$80$\mu$s, $1/\gamma_{gg}=$40$\mu$s respectively. The Doppler shifts $\delta_{i}$ were assumed to be Gaussian-distributed with an rms width of $2\pi\cdot 43.5$kHz. Direct spontaneous decay processes from the Rydberg states, which occur over longer timescales, were neglected. Numerical solutions of the master equation (29) were performed using QuTiP Johansson et al. (2013), and observables were averaged over 100 disorder realizations $\{\delta_{i}\}$. Uncertainties in observables were computed from the standard error of the mean of these realizations. We note that the experiment has additional loss mechanisms, as well as imperfections in the laser sweep profile, which are not well characterized and not included in this Lindbladian model. We believe this accounts for the discrepancy with experimental correlation functions noted in the main text.

This master equation predicts a substantial loss in purity $\textnormal{Tr}\left[\hat{\rho}^{2}\right]$ for states produced at the end of the sweep (Fig. 10), whose detrimental effects on our pure-state reconstruction process we quantify below.

V.2 Reconstruction fidelities

To assess the quality of quantum state reconstruction, we consider the fidelity between two states $\hat{\rho},\hat{\sigma}$,

[TABLE]

which reduces to the norm of the overlap in the case where $\hat{\rho},\hat{\sigma}$ are pure states. An ideal state reconstruction $\hat{\sigma}$ of a mixed state $\hat{\rho}$ would yield $F(\hat{\rho},\hat{\sigma})=1$. For pure state reconstructions $\hat{\sigma}=|\psi_{\bm{\lambda}}\rangle\langle\psi_{\bm{\lambda}}|$, this is not possible if the true state $\hat{\rho}$ is non-pure. However, one may still seek an approximate reconstruction which reproduces the local reduced density operators of $\hat{\rho}$. In particular, specializing to the case of one-dimensional systems, we can consider contiguous subsystems formed from $s$ adjacent sites, $A^{(s)}_{i}=\{i,i+1,...,i+s-1\}$. Given two density operators $\hat{\rho},\hat{\sigma}$ for the global system $\mathcal{S}$ of size $N$, the reduced density operators which describe the subsystem in each state are obtained by tracing out the rest of the chain,

[TABLE]

Then we define a subsystem averaged fidelity as the spatial average of the fidelity between these local operators, over all subsystems of a particular size $s$:

[TABLE]

$F_{s}\left(\hat{\rho},\hat{\sigma}\right)$ is a measure of how well, on average, $\hat{\sigma}$ is able to reproduce the $s$-local physics of $\hat{\rho}$.

To examine the quality of the RBM states $\hat{\rho}_{\bm{\lambda}}=|\psi_{\bm{\lambda}}\rangle\langle\psi_{\bm{\lambda}}|$ in reproducing local density operators, we solved the master equation (29) for set of decay rates $\gamma_{rg}=\alpha\gamma_{rg}^{\textnormal{exp}},\gamma_{gg}=\alpha\gamma_{gg}^{\textnormal{exp}}$, with $\gamma_{rg}^{\textnormal{exp}},\gamma_{gg}^{\textnormal{exp}}$ denoting our estimates of the experimental values, and $\alpha$ a dimensionless parameterization of the overall decoherence strength. For each set of decoherence rates, the master equation was solved and synthetic data sampled from the resulting mixed states. Pure state RBMs were trained on each of these datasets, and the resulting averaged fidelities $F_{s}\left(\hat{\rho},\hat{\rho}_{\bm{\lambda}}\right)$ were computed.

As a representative example, Fig. 11 shows how the fidelities computed in the final state of the sweep vary as a function of the average Renyi entropy $\bar{s_{2}}$ of subsystems of a given size – one observes a roughly linear decay in the average fidelity with the averaged entropy. These numerical results suggest that pure state reconstruction techniques should focus on few-body operators, where the entropy build-up due to global decoherence process is limited in proportion to the system size.

VI Reconstruction improvement from noise layer regularization

Numerical experiments have demonstrated that noise layer regularization results in higher-fidelity pure state reconstruction when training on uniformly noisy data.

Fig. 12 compares fidelities achieved by regularized and unregularized RBMs, when trained on synthetic datasets subjected to the bitflip error channel described in the main text. The improvement is significant, especially in the ordered phase, where global state purity is lowest. It is important to note that in these experiments the noise process is known ahead of time and built into the three layer networks as in Fig. 5. We have also trained three-layer machines using incorrect values of the error rates on the same noisy synthetic data. Although the fidelity performance varied somewhat, depending most sensitively on $p(0|1)$, the quality of the regularized reconstructions is generally robust, and deep in the ordered phase all three-layer machines exhibited higher fidelities than their two layer counterparts on the corresponding datasets, for error rates with bounds set by single-atom measurements (Levine et al., 2018). Generically, of course, a sufficiently large mismatch between the true and assumed error rates will lead to decreased reconstruction fidelity. Future work will investigate more generally the task of selecting a regularization method for noisy quantum data.

Fig. 13 compares the predictions of these synthetically trained two- and three-layer machines (using the known noise values) for some of the observables discussed in the main text. We find that noise-layer training allows the RBMs to provide much tighter agreement in, for example, values of the transverse field and mutual information. Surprisingly, the three-layer machines actually produce poorer estimates of the transverse field correlator in the ordered phase, despite yielding two-body density operators with higher fidelities for all sampled states. A more detailed analysis of the ordered phase states reveals that regularized training does indeed produce better estimates of the one- and two-body expectation values $\langle\sigma_{i}^{x}\rangle$, $\langle\hat{\sigma}_{i}^{x}\hat{\sigma}_{i+1}^{x}\rangle$ for all sites $i$. However, at bonds (3, 4) and (5, 6), where quantum fluctuations are strongest, the three-layer improvement in the two-body expectation value is relatively small, while the reduction in the one-body expectation value is substantial. Upon computing the connected correlator $\langle\hat{\sigma}_{i}^{x}\hat{\sigma}_{i+1}^{x}\rangle-\langle\sigma_{i}^{x}\rangle\langle\sigma_{i+1}^{x}\rangle$, the overall effect is an overestimate of the true correlation.

VII Appendix: Proof of classical fidelity bound

Inequality (28) is obtained by bounding the probability of the most likely outcome using the Renyi entropy $H_{2}$ in the measurement basis. By definition, $H_{2}=-\log\sum_{\bm{\tau}}P_{\textnormal{GT}}(\bm{\tau})^{2}$, and $\sum_{\bm{\tau}}P_{\textnormal{GT}}(\bm{\tau})^{2}\geq\max_{\bm{\tau}}P_{\textnormal{GT}}(\bm{\tau})^{2}$, so $-\log\sum_{\bm{\tau}}P_{\textnormal{GT}}(\bm{\tau})^{2}\leq-2\log\max_{\bm{\tau}}P_{\textnormal{GT}}(\bm{\tau})$. Rearranging, $\max_{\bm{\tau}}P_{\textnormal{GT}}(\bm{\tau})\leq e^{-H_{2}/2}$. In particular, this bounds the probability of any event in the training set, so

[TABLE]

VIII Appendix: Renyi entropy bound from positive pure states

The $n$th order Renyi entropy of a quantum state $\hat{\rho}$ is defined as $S_{n}\left[\hat{\rho}\right]=\frac{1}{1-n}\log\textnormal{Tr}\hat{\rho}^{n}$.

Consider a system $S$ partitioned into subsets $A$ and $B$, and a density operator $\hat{\rho}$ defined on $S$; its reduced density operator in the $A$ subsystem is $\hat{\rho}_{A}=\textrm{Tr}_{B}\hat{\rho}$. Let $|i\rangle$, $|j\rangle$ denote orthonormal bases for $A,B$ respectively, so that the set of product states $|i,j\rangle$ forms an orthonormal basis for the full system $S$. Let $p_{i,j}$ be the probability assigned by $\hat{\rho}$ to the measurement outcome $i,j$: $p_{i,j}=\textnormal{Tr}\left(\hat{\rho}|i,j\rangle\langle i,j|\right)$. The positive-pure partner to the mixed state is defined as

[TABLE]

and the corresponding reduced density operator on $A$ is $\hat{\rho}^{P}_{A}=\textrm{Tr}_{B}|\Psi^{P}[\rho]\rangle\langle\Psi^{P}[\rho]|$.

**Theorem: ** For $n>1$, the Renyi entropies $S_{n}$ of the two density operators satisfy the inequality

[TABLE]

As a consequence, in the case of pure states $\hat{\rho}$, where the global Renyi entropy vanishes, the positive-pure partner provides a lower bound on the mutual information:

[TABLE]

We note that for the case of the $n=2$ Renyi entropy and pure states $\hat{\rho}$, this result has been obtained in previous work Zhang et al. (2011); Grover and Fisher (2015).

Proof: Choose an auxiliary system $R$ to purify $\hat{\rho}$: $\hat{\rho}=\textnormal{Tr}_{R}|\Psi\rangle\langle\Psi|$ for some pure state $|\Psi\rangle$ living in $S\otimes R$. If $|\alpha\rangle$ is an orthonormal basis for $R$, we can expand the larger pure state in the joint basis $|i,j,\alpha\rangle$ : $|\Psi\rangle=\sum_{i,j,\alpha}\Psi_{i,j}^{\alpha}|i,j,\alpha\rangle$ for some complex coefficients $\Psi_{i,j}^{\alpha}$.

In terms of these amplitudes, the reduced density operator of the mixed state on $A$ is

[TABLE]

and so

[TABLE]

with summation over all indices implied. The reduced density operator for positive-pure partner may be obtained from the definition above:

[TABLE]

whence

[TABLE]

(summation implied). Furthermore,

[TABLE]

and so by the Cauchy-Schwartz inequality,

[TABLE]

Therefore, writing $\mathbf{i}=(i_{1},...,i_{n}),\mathbf{j}=(j_{1},...,j_{n})$,

[TABLE]

Hence $-\log\textnormal{Tr}\left(\hat{\rho}^{P}_{A}\right)^{n}\leq-\log\textnormal{Tr}\left(\hat{\rho}_{A}\right)^{n}$, which means that for $n>1$, $S_{n}(\hat{\rho}^{P}_{A})\leq S_{n}(\hat{\rho}_{A})$.