Machine learning quantum states in the NISQ era

TL;DR

This paper reviews machine learning techniques, especially restricted Boltzmann machines, for reconstructing noisy quantum states on NISQ devices, highlighting recent experimental advances and future prospects.

Contribution

It provides a comprehensive overview of generative modeling for quantum state reconstruction, emphasizing practical applications on NISQ hardware and discussing future directions.

Findings

Successful reconstruction of classical and quantum states using RBMs

Application to experimental cold atom wavefunctions

Discussion of future NISQ-era quantum state tomography

Abstract

We review the development of generative modeling techniques in machine learning for the purpose of reconstructing real, noisy, many-qubit quantum states. Motivated by its interpretability and utility, we discuss in detail the theory of the restricted Boltzmann machine. We demonstrate its practical use for state reconstruction, starting from a classical thermal distribution of Ising spins, then moving systematically through increasingly complex pure and mixed quantum states. Intended for use on experimental noisy intermediate-scale quantum (NISQ) devices, we review recent efforts in reconstruction of a cold atom wavefunction. Finally, we discuss the outlook for future experimental state reconstruction using machine learning, in the NISQ era and beyond.

Figures

Equations

$$p(\mathbf{x})=P_{\text{data}}(\mathbf{x})=\frac{1}{\Vert \mathcal{D}\Vert }\sum _{{\mathbf{x}}_k\in \mathcal{D}}{\delta }_{\mathbf{x},{\mathbf{x}}_k}.$$

$${\mathbb{KL}}_{\mathbf{\lambda }}(P\Vert p)=\sum _{\mathbf{x}}P(\mathbf{x})\log \frac{P(\mathbf{x})}{p_{\mathbf{\lambda }}(\mathbf{x})},$$

$${\mathbb{KL}}_{\mathbf{\lambda }}(P\Vert p)\approx -\frac{1}{\Vert \mathcal{D}\Vert }\sum _{\mathbf{x}\in \mathcal{D}}\log p_{\mathbf{\lambda }}(\mathbf{x})-{\mathbb{H}}_{\mathcal{D}},$$

$$E(\mathbf{x})=-\sum _{ij}{W}_{ij}{x}_i{x}_j.$$

$$x_i=\{\begin{array}{ll}1,& \text{mbox}if{\sum }_j{W}_{ij}{x}_j>0\\ 0,& \text{mbox}otherwise.\end{array}$$

$$p_i=\frac{1}{1+e^{-\beta \Delta E_i}},$$

$$p(\mathbf{x})=\frac{1}{Z}{e}^{-\beta E(\mathbf{x})},Z=\sum _{\mathbf{x}}{e}^{-\beta E(\mathbf{x})}.$$

$$\Delta W_{ij}=\beta\Big{[}\langle x_{i}x_{j}\rangle_{\mathcal{D}}-\langle x_{i}x_{j}\rangle_{p(\bm{x})}\Big{]}\>.$$

$$E_{\mathbf{\lambda }}(\mathbf{v},\mathbf{h})=-\sum _{ij}{W}_{ij}{h}_i{v}_j-\sum _j{b}_j{v}_j-\sum _i{c}_i{h}_i,$$

$$p_{\mathbf{\lambda }}(\mathbf{v})=\sum _{\mathbf{h}}{p}_{\mathbf{\lambda }}(\mathbf{v},\mathbf{h})=\frac{1}{Z_{\mathbf{\lambda }}}\sum _{\mathbf{h}}{e}^{-E_{\mathbf{\lambda }}(\mathbf{v},\mathbf{h})}=\frac{1}{Z_{\mathbf{\lambda }}}{e}^{-{\mathcal{E}}_{\mathbf{\lambda }}(\mathbf{v})},$$

$${\mathcal{E}}_{\mathbf{\lambda }}(\mathbf{v})=-\sum _j{b}_j{v}_j-\sum _i\log \left(1+e^{{\sum }_j{W}_{ij}{v}_j+c_i}\right).$$

$$\begin{array}{rl}{\mathcal{C}}_{\mathbf{\lambda }}& =-\frac{1}{\Vert \mathcal{D}\Vert }\sum _{\mathbf{v}\in \mathcal{D}}\log p_{\mathbf{\lambda }}(\mathbf{v})=-\frac{1}{\Vert \mathcal{D}\Vert }\sum _{\mathbf{v}\in \mathcal{D}}{\mathcal{E}}_{\mathbf{\lambda }}(\mathbf{v})-\log Z_{\mathbf{\lambda }},\end{array}$$

$$\Delta\bm{\lambda}\propto-\nabla_{\bm{\lambda}}\mathcal{C}_{\bm{\lambda}}=\big{\langle}\nabla_{\bm{\lambda}}\mathcal{E}_{\bm{\lambda}}(\bm{v})\big{\rangle}_{\mathcal{D}}-\big{\langle}\nabla_{\bm{\lambda}}\mathcal{E}_{\bm{\lambda}}(\bm{v})\big{\rangle}_{p_{\bm{\lambda}}(\bm{v})}\>$$

$$\big{\langle}\nabla_{\bm{\lambda}}\mathcal{E}_{\bm{\lambda}}(\bm{v})\big{\rangle}_{\mathcal{D}}=\frac{1}{\|\mathcal{D}\|}\sum_{\bm{v}\in\mathcal{D}}\nabla_{\bm{\lambda}}\mathcal{E}_{\bm{\lambda}}(\bm{v})\>.$$

$$\begin{split}\big{\langle}\nabla_{\bm{\lambda}}\mathcal{E}_{\bm{\lambda}}(\bm{v})\big{\rangle}_{p_{\bm{\lambda}}(\bm{v})}&=\frac{1}{Z_{\bm{\lambda}}}\sum_{\bm{v}}p_{\bm{\lambda}}(\bm{v})\nabla_{\bm{\lambda}}\mathcal{E}_{\bm{\lambda}}(\bm{v})\\ &\approx\frac{1}{M}\sum_{k=1}^{M}\nabla_{\bm{\lambda}}\mathcal{E}_{\bm{\lambda}}(\bm{v}_{k})\>,\end{split}$$

$$p_{\mathbf{\lambda }}(\mathbf{v}|\mathbf{h})=\prod _j{p}_{\mathbf{\lambda }}(v_j|\mathbf{h}),p_{\mathbf{\lambda }}(\mathbf{h}|\mathbf{v})=\prod _i{p}_{\mathbf{\lambda }}(h_i|\mathbf{v}),$$

$${\mathbb{CD}}_k=\mathbb{KL}(P_{\text{data}}|p_{\mathbf{\lambda }})-\mathbb{KL}(p_{\mathbf{\lambda }}^{(k)}|p_{\mathbf{\lambda }}),$$

$$\Delta\bm{\lambda}\propto\big{\langle}\nabla_{\bm{\lambda}}\mathcal{E}_{\bm{\lambda}}(\bm{v})\big{\rangle}_{\mathcal{D}}-\big{\langle}\nabla_{\bm{\lambda}}\mathcal{E}_{\bm{\lambda}}(\bm{v})\big{\rangle}_{p^{(k)}_{\bm{\lambda}}(\bm{v})}\>.$$

$$\hat{\varrho }=\frac{e^{-\beta \hat{H}}}{\text{Tr}[e^{-\beta \hat{H}}]}=\sum _{\mathbf{\sigma }}{P}_{\beta }(\mathbf{\sigma })|\mathbf{\sigma }⟩⟨\mathbf{\sigma }|$$

$$H(\mathbf{\sigma })=-\sum _{⟨ij⟩}{\sigma }_i{\sigma }_j$$

$${\psi }_{\mathbf{\lambda }}(\mathbf{\sigma })=\sqrt{p_{\mathbf{\lambda }}(\mathbf{\sigma })}=\frac{1}{\sqrt{Z_{\mathbf{\lambda }}}}{e}^{-{\mathcal{E}}_{\mathbf{\lambda }}(\mathbf{\sigma })/2}.$$

$$\begin{array}{rl}{\mathcal{C}}_{\mathbf{\lambda }}& =\sum _{\mathbf{\sigma }}|\Psi (\mathbf{\sigma }){|}^2\log \frac{|\Psi (\mathbf{\sigma }){|}^2}{|{\psi }_{\mathbf{\lambda }}(\mathbf{\sigma }){|}^2}\\ & =-\frac{1}{\Vert \mathcal{D}\Vert }\sum _{{\mathbf{\sigma }}_k\in \mathcal{D}}\log p_{\mathbf{\lambda }}(\mathbf{\sigma })-{\mathbb{H}}_{\mathcal{D}},\end{array}$$

$$⟨\hat{\mathcal{O}}⟩=⟨{\psi }_{\mathbf{\lambda }}|\hat{\mathcal{O}}|{\psi }_{\mathbf{\lambda }}⟩=\sum _{\mathbf{\sigma }}{p}_{\mathbf{\lambda }}(\mathbf{\sigma }){\mathcal{O}}_{\mathbf{\sigma \sigma }},$$

$${\mathcal{O}}_L(\mathbf{\sigma })=\sum _{{\mathbf{\sigma }}^{\prime }}\frac{{\psi }_{\mathbf{\lambda }}({\mathbf{\sigma }}^{\prime })}{{\psi }_{\mathbf{\lambda }}(\mathbf{\sigma })}{\mathcal{O}}_{{\mathbf{\sigma }}^{\prime }\mathbf{\sigma }},$$

$$S_{\alpha }({\hat{\rho }}_A)=\frac{1}{1-\alpha }\log \text{Tr}({\hat{\rho }}_A^{\alpha }),$$

$$\begin{split}S_{2}(\hat{\rho}_{A})&=-\log\text{Tr}(\hat{\rho}_{A}^{2})\\ &=-\log\Big{[}\langle\Psi_{1}|\otimes\langle\Psi_{2}|\text{Swap}_{A}|\Psi_{1}\rangle\otimes|\Psi_{2}\rangle\Big{]}\>.\end{split}$$

$$\hat{H}=-\sum _{⟨ij⟩}{\hat{\sigma }}_i^z{\hat{\sigma }}_j^z-h\sum _i{\hat{\sigma }}_i^x.$$

$${\psi }_{\mathbf{\lambda \mu }}(\mathbf{\sigma })=\frac{1}{\sqrt{Z_{\mathbf{\lambda }}}}{e}^{-({\mathcal{E}}_{\mathbf{\lambda }}(\mathbf{\sigma })+i{\mathcal{E}}_{\mathbf{\mu }}(\mathbf{\sigma }))/2}.$$

$$\Psi ({\mathbf{\sigma }}^{\mathbf{b}})=\sum _{\mathbf{\sigma }}{\mathcal{U}}_{{\mathbf{\sigma }}^{\mathbf{b}}\mathbf{\sigma }}\Psi (\mathbf{\sigma }),$$

$$\begin{split}\mathcal{C}_{\bm{\lambda\mu}}&=-\frac{1}{\|\mathcal{D}\|}\sum_{\bm{\sigma^{b}}\in\mathcal{D}}\log|\psi_{\bm{\lambda\mu}}(\bm{\sigma}^{\bm{b}})|^{2}\\ &=-\frac{1}{\|\mathcal{D}\|}\sum_{\bm{\sigma^{b}}\in\mathcal{D}}\Bigg{[}\log\left(\sum_{\bm{\sigma}}\mathcal{U}_{\bm{\sigma}^{\bm{b}},\bm{\sigma}}{\psi}_{\bm{\lambda\mu}}(\bm{\sigma})\right)+c.c\Bigg{]}\>,\end{split}$$

Full text

I Introduction

We are entering the age where quantum computers with tens – or soon hundreds – of qubits are becoming available. These noisy intermediate-scale quantum (NISQ) devices Preskill (2018) are being constructed out of cold atoms Bernien et al. (2017), superconducting quantum circuits Kandala et al. (2017), trapped ions Zhang et al. (2017), and other quantum systems for which we have achieved an exquisite degree of control. NISQ devices will soon play an important role, since they are poised to surpass the ability of the world’s most powerful computers to perform exact simulations of them, ushering in the era of so-called quantum supremacy Terhal (2018).

Some of the first tasks for these NISQ devices will be as simulators (or emulators) of other highly-entangled quantum many-body systems. The goal is to supplant our current conventional computer simulation technology, such as exact diagonalization, quantum Monte Carlo, or tensor network methods. Efforts to produce quantum simulators of some of the most important physical systems, such as the fermionic Hubbard model Mazurenko et al. (2017), are progressing in earnest. However, with the advent of increasingly larger NISQ devices comes a paradox: how will we simulate the simulators? That is, how will we validate an intermediate-scale quantum device, confirming that it is producing the behavior it was designed for? Along with quantum supremacy comes the necessary breakdown of conventional tomography – the gold standard for quantum state reconstruction. We are left searching for imperfect alternatives.

The answer may lie in new data-driven approaches inspired by rapid advances in machine learning. A strategy for unsupervised learning, called generative modeling, has demonstrated the ability to integrate well with the data produced by NISQ devices. In industry applications, the goal of generative modeling is to reconstruct an unknown probability distribution $P(\bm{x})$ from a set of data $\bm{x}$ drawn from it. In the most powerful versions of generative modeling, the reconstructed probability distribution is represented approximately by a graphical model or neural network – the weights and biases serving as a parameterization of $P(\bm{x})$. After training, these generative models can be used to estimate the likelihood, or to produce samples, of new $\bm{x}$ in a way that generalizes and scales well.

This procedure can be extended to data produced by quantum devices, with the goal of reconstructing the quantum wavefunction (a complex generalization of a classical probability distribution). NISQ devices with single-site control are particularly suited to this data-driven approach, since they can produce projective measurements of the state of individual qubits. If a sufficient type and number of projective measurements can be obtained, industry-standard algorithms for unsupervised learning of the relevant probability distributions (produced according to the Born rule) can be used to reconstruct the underlying quantum state.

Such data-driven state reconstruction may play by different rules than Hamiltonian-driven discovery of quantum states. In fact, the latter consists of obtaining a quantum state underlying a microscopic model (i.e. a Hamiltonian), and it is a benchmark for quantum supremacy. Instead, the former assumes no knowledge of the Hamiltonian, but requires informationally-complete sets of measurement data on the quantum state. The question of how efficiently the data-driven approach scales for wavefunctions of various structures of interest to physicists, and how it compares to the more conventional Hamiltonian-driven approach, is largely unanswered.

The most obvious role for a quantum state reconstructed via generative modeling is to produce new physical observables. To be useful, this must be done in a tractable way that scales efficiently with increasing number of qubits, while generalizing well to unseen data. The observables in question may be inaccessible to the device, such as those encoded in a basis for which no projective measurement was taken, or those (such as Renyi entanglement entropies) that require elaborate technical setups Islam et al. (2015). Generative models are also capable of mitigating noise in the state preparation and measurement, a ubiquitous and defining condition in NISQ devices. Finally, the ability to off-load the production of various observables to a parameterized model frees experimentalists to focus solely on the production of high-quality projective measurements. It is this type of inelegant compromise that will allow machine learning techniques to contribute to the verification of quantum devices as they grow into the NISQ era and beyond.

In this paper, we review the development of generative modeling for quantum state reconstruction. Beginning with the classical treatment of probability distributions, we motivate the use of a restricted Boltzmann machine (RBM), and demonstrate its ability to parameterize the thermal distribution of data drawn from a classical Ising model. The same type of RBM is shown to faithfully reconstruct a real-positive wavefunction, and we demonstrate the production of non-trivial observables from the parameterized model. We then discuss extensions of standard RBMs to reconstruct complex wavefunctions and density matrices. We end this review with a discussion of recent efforts in reconstruction of a real-world Rydberg atom quantum simulator. Despite challenges in noisy state preparation and measurement, the demonstration of real experimental state reconstruction is a milestone for the use of machine learning in the NISQ era.

II Generative modeling

Let us begin by considering an unknown probability distribution $P(\bm{x})$ defined over the $2^{N}$-dimensional space of binary states $\bm{x}=(x_{1},\dots,x_{N})$, and a set of data $\mathcal{D}=\{\bm{x}_{k}\}$ distributed according to $P(\bm{x})$. Can we infer the features of such distribution, such as regularities and correlations, directly from the observation of the data? In other words, can we discover an approximate representation $p(\bm{x})\approx P(\bm{x})$ from the limited-size dataset $\mathcal{D}$? The simplest approach consists of approximating the unknown probability with the frequency distribution obtained by inverting the measurement counts in the dataset:

[TABLE]

The validity of this approximation depends on the size of the system $N$ the entropy of the distribution and the size $\|\mathcal{D}\|$ of the dataset. For most practical purposes however, it fails to generalize the features of $P(\bm{x})$ beyond the training set. In contrast, generative modeling aims to discover an approximation of the unknown distribution that captures the underlying structure and it is also capable of generalization.

The first ingredient is a compact representation of the probability distribution, i.e. a parametrization $p_{\bm{\lambda}}(\bm{x})$ in terms of a set of parameters $\bm{\lambda}$ whose number is much smaller than the size of the configuration space. Then, generative modeling consists of finding an optimal set of parameters $\bm{\lambda}^{*}$ such that the parametric distribution $p_{\bm{\lambda}^{*}}(\bm{x})$ mimics the unknown distribution $P(\bm{x})$ underlying the finite number of dataset samples. In practice, this search is carried out through an optimization procedure, where the distance between the two distributions is minimized with respect to the model parameters $\bm{\lambda}$. The distance between two probability distributions can be quantified by the Kullbach-Leibler (KL) divergence

[TABLE]

a non-symmetric statistical measure such that $\mathbb{KL}_{\bm{\lambda}}(P\,\|\,p_{\bm{\lambda}})>0$ and $\mathbb{KL}_{\bm{\lambda}}(P\,\|\,p_{\bm{\lambda}})=0$ if and only if $P=p_{\bm{\lambda}}$. By approximating the KL divergence with the measurement data, we obtain

[TABLE]

where $\mathbb{H}_{\mathcal{D}}$ is the dataset entropy. This quantity can be minimized iteratively by one of the many variants of the gradient descent algorithm. This procedure allows one to obtain a representation of the unknown distribution and generate new configurations that were not encountered in the learning stage. The most successful approach relies on the representation of $p_{\bm{\lambda}}(\bm{x})$ in terms of networks of artificial neurons.

II.1 Artificial neural networks

Artificial neural networks, the bedrock of modern machine learning and artificial intelligence LeCun et al. (2015), have a history spanning decades. Initially investigated to understand the process of human cognition, neural networks models are based on the idea that information (in the brain) has a distributed representation over a large collection of elementary units (neurons), and information processing occurs through the mutual interaction between neurons Rumelhart and McClelland (1986). The fundamental ingredients are: i) a set of neurons, each one applying a simple type of computation to the input signal it receives; ii) a set of interactions defined over a graph structure connecting the neurons; iii) an external environment providing a “teaching signal”; iv) a learning rule, i.e. a prescription for modifying the interactions according to the external environment.

The first artificial neuron capable of computation, the perceptron, was proposed by Frank Rosenblatt as early as 1957 Rosenblatt (1958). Based on the previous work of McCullogh and Pitts McCulloch and Pitts (1943), the perceptron was capable of discriminating different classes of input patterns, a process called supervised learning. It was later shown that a single layer perceptron is only capable of learning linearly separable functions Minsky and Papert (1969), and since no learning algorithms were known for multi-layer perceptrons, the model was abandoned, leading to a decrease in both popularity and fundings of neural networks (called the first AI winter). The first resurgence of the field took place more than a decade later, with the invention of the backpropagation algorithm Rumelhart et al. (1986) and the Boltzmann machine (BM) Ackley et al. (1985). The latter, was directly built on the connection between cognitive science and statistical mechanics, established by the works of condensed matter physicists William Little Little (1974); Little and Shaw (1978) and John Hopfield Hopfield (1982).

II.1.1 The Hopfield model

The Hopfield network, introduced in 1982 as a model for associative memories Hopfield (1982), was inspired by the concept of emergence in condensed matter physics, where complex behaviors effectively emerge from the mutual interactions of a large number of degrees of freedom. In this context, Hopfield formulated a physics-inspired model of cognition for the task of recovering a corrupted memory. By regarding a memory as a state $\bm{x}$ containing $N$ bits of information, the corresponding network consists of $N$ binary neurons fully connected with symmetric weights (or interactions), described by an energy function

[TABLE]

Each neuron in the network carries out the computation of $\sum_{j}W_{ij}x_{j}$, and update itself according to the following rule:

[TABLE]

Since the energy difference between the two possible states of the $i$-th neuron is $\Delta E_{i}=\sum_{j}W_{ij}x_{j}$, the dynamics resulting from the asynchronous update of the neurons monotonically minimizes the total energy. Therefore, given an initial state, the network evolves in time by following the above equation of motion until a stable configuration (i.e. a local energy minimum) is found. In the context of associative memories, given a set of desired memory states $\{\bar{\bm{x}}_{1},\bar{\bm{x}}_{2},\dots\}$, there exists a learning rule to modify the interactions $\bm{W}$ in such a way that these states become local minima in the energy landscape Hopfield (1982). Thus, if the network is initialized to a corrupted memory $\bar{\bm{x}}_{k}+\bm{\delta}$ which is sufficiently close to the true state (i.e. small $\bm{\delta}$), the network is able to recover the correct memory simply by evolving with its equations of motion.

II.1.2 The Boltzmann machine

The two major limitations of the Hopfield model are the tendency for the network to get trapped into local minima, and its restricted capacity. Nevertheless, it suggested the important connection between cognitive science and statistical physics, which was was further strengthened with the invention of the BM by Ackley, Hinton and Sejnowski in 1985 Ackley et al. (1985). Similarly to the Hopfield model, the BM consists of a set of $N$ binary neurons interacting with energy given in Eq. 4. However, in order to allow the system to escape local minima, the neural network is placed at thermal equilibrium at some inverse temperature $\beta=1/T$. The update rule becomes now stochastic, with the $i$-th neuron activating ($x_{i}=1$) with probability

[TABLE]

where $\Delta E_{i}$ is once again the energy difference between its two internal states. As the temperature goes to zero ($\beta\rightarrow\infty)$, one recovers the Hopfield model with a deterministic dynamics minimizing energy. Instead, for a stochastic dynamics at finite temperature, the network minimizes the free energy instead, equilibrating to the canonical Boltzmann distribution:

[TABLE]

The BM is one of the simplest examples of generative models. In fact, the set of interactions can be considered as tunable parameters, resulting into a parametric distribution $p_{\bm{\lambda}}(\bm{x})$ (with $\bm{\lambda}=\bm{W}$). Then, the interactions can be modified following an unsupervised learning procedure in order for the network distribution to mimic an unknown probability distribution underlying a given set of data points $\mathcal{D}=\{\bm{x}\}$. By minimizing the statistical divergence between the data and model distribution (Eq. 3), one obtains the following learning rule for the parameters Ackley et al. (1985),

[TABLE]

In the positive phase, the weight $W_{ij}$ is increased according to the average value of $x_{i}x_{j}$ over the data points in $\mathcal{D}$, corresponding to a traditional Hebbian learning (i.e. “neurons that fire together wire together”). This term effectively lowers the energy of all configurations that are compatible with the dataset, thus increasing their probability. In contrast, in the negative phase, the same process occurs with the reverse sign, decreasing the probability of configurations generated by the BM when running freely at thermal equilibrium. Clearly, as the two averages coincide, the BM distribution reproduces the dataset and there is not net change in the parameter. Otherwise, the network is trying to unlearn configurations generated at equilibrium that lead to an imbalance with respect to the data. It is interesting to note how this learning and unlearning process had been already proposed in an ad-hoc way by Hopfield to eliminate spurious minima in his model of associate memories Hopfield et al. (1983).

The major limitation of this network is the structure of the energy function, allowing the BM to capture only pairwise correlation in the data (e.g. it cannot learn the XOR function Hinton and Sejnowski (1986)). The simplest way to increase the reach of its representational capabilities is to introduce an auxiliary set of neurons which do not appear in the input space of the data. The full network is then divided as $\bm{x}=(\bm{v},\bm{h})$, where $\bm{v}$ are called visible units, corresponding to the degrees of freedom in the dataset, and $\bm{h}$ are called hidden units Ackley et al. (1985). In order to derive a learning rule for the network parameters, one needs to eliminate the hidden degrees of freedom so that an explicit distribution over the visible neurons $p_{\bm{\lambda}}(\bm{v})=\sum_{\bm{h}}p_{\bm{\lambda}}(\bm{v},\bm{h})$ can be obtained. Therefore, to attain a tractable marginal distribution, one can restrict the interactions between neurons on different layers, resulting in the famous restricted Boltzmann machine (RBM).

II.2 Restricted Boltzmann machines

The RBM, originally introduced by Smolensky under the name of Harmonium Smolensky (1986), is a probabilistic graphical model with energy

[TABLE]

where we have added bias terms (i.e. magnetic fields) $\bm{b}$ and $\bm{c}$ for the visible and hidden layer respectively. The set of tunable parameters is now $\bm{\lambda}=(\bm{W},\bm{b},\bm{c})$ (Fig. 1b). The marginal distribution, obtained by tracing out the hidden neurons, can be calculated analytically

[TABLE]

where we set the inverse temperature to $\beta=1$ and we introduced the new energy function

[TABLE]

The energy $\mathcal{E}_{\bm{\lambda}}(\bm{v})$ defines an effective system consisting of the visible neurons only. We can see that the energy contains two terms: a mean-field contribution, proportional to the visible bias $\bm{b}$ and a non-linearity containing correlations between visible neurons at all orders. The particular structure of such an effective energy allows the RBM distribution $p_{\bm{\lambda}}(\bm{v})$ to be a universal function approximator of discrete distributions Le Roux and Bengio (2008). This means that, given a large enough number of hidden neurons, any function of discrete binary variables can be approximated to arbitrary precision. However, in the worst case scenario, the number of hidden neurons may grow exponentially with the visible layer.

II.2.1 Unsupervised learning

The goal of unsupervised learning is to discover a set of parameters so that the RBM distribution mimics the unknown distribution underlying a dataset $\mathcal{D}=\{\bm{v}_{1},\bm{v}_{2},\dots\}$ of visible samples. The cost function, given by the KL divergence from Eq. 3, is

[TABLE]

where we have omitted the constant entropy term $\mathbb{H}_{\mathcal{D}}$. The learning rule for the RBM parameters is obtained by taking the gradient of $\mathcal{C}_{\bm{\lambda}}$

[TABLE]

where the gradients $\nabla_{\bm{\lambda}}\mathcal{E}_{\bm{\lambda}}(\bm{v})$ are straightforwards to calculate. Similar to the regular BM, the gradient contains two competing terms, driven respectively by the data and the RBM distribution. The first term (the positive phase) is trivial to compute, being an average over the data:

[TABLE]

Conversely, the calculation of the negative phase is in general intractable. It needs to be approximated using a Monte Carlo simulation,

[TABLE]

where the configurations $\bm{v}_{k}$ are drawn from a Markov chain running on the distribution $p_{\bm{\lambda}}(\bm{v})$.

The sampling stage to estimate the negative phase, which is the computational bottleneck of the training, is aided by the restricted nature of the RBM graph. In fact, because of that, neurons in a given layer are conditionally independent of one another. That is, due to the lack of intra-layer connections, the conditional probabilities for the neurons in one layer conditioned on the current state of the other factorize over each individual neuron,

[TABLE]

and can be easily calculated analytically Fischer and Igel (2012). When running the Markov chain to collect the statistic in Eq. 15, one can sample the state of each neuron in one layer simultaneously using the above conditional probabilities, alternating between visible and hidden layer. This sampling strategy is called block Gibbs sampling.

II.2.2 Training by contrastive divergence

The calculation of the negative phase, even if carried out using block Gibbs sampling, is still computationally intensive. In fact, at each training iteration, the RBM needs to reach its equilibrium distribution $p_{\bm{\lambda}}(\bm{v})$ before collecting the statistics for the negative phase calculation. Furthermore, the gradient in Eq. 13 can display a large variance, being the difference of two averages computed from two different distributions. A solution to both these issues is to consider a different cost function. Namely, the contrastive divergence (CD) between the data and the RBM after a sequence of $k$ block Gibbs sampling steps is Hinton (2002),

[TABLE]

where $p^{(k)}_{\bm{\lambda}}$ is the probability distribution of the visible layer after $k$ steps. The new update from the gradient of $\mathbb{CD}_{k}$ becomes Bengio and Delalleau (2009); Fischer and Igel (2011); Carreira-Perpiñán and Hinton (2005)

[TABLE]

The resulting CD training consists of initializing the RBM to a random sample from the dataset $\mathcal{D}$ and using the visible state after $k$ steps of block Gibbs sampling to evaluate the negative phase.

Once the gradient of the cost function is calculated, the parameters $\bm{\lambda}$ are updated with a gradient descent algorithm. The simplest one, called stochastic gradient descent, uses a random set of data to evaluate the positive phase and performs the update $\Delta\bm{\lambda}=-\eta\nabla_{\bm{\lambda}}\mathcal{C}_{\bm{\lambda}}$ where $\eta$ is the step-size of the update, also called the learning rate. The total number of data samples used for the update is called the batch size. Other algorithms can be used to speed up the convergence Sutskever et al. (2013) and tune the learning rate in an adaptive way Zeiler (2012); Kingma and Ba (2014). Furthermore, an additional term should be added to the cost function to help generalization, i.e. avoid the overfitting of data points. A common choice is weight decay Krogh and Hertz (1992) regularization, which penalizes large value of the weights. We refer the reader to Ref Hinton (2012) for more details on the practical training of RBMs and a description of the various training hyperparameters (and how to choose them).

III Quantum state reconstruction

Image that an experimental NISQ apparatus in the laboratory containing $N$ qubits is prepared in some quantum state of interest, described by a density operator $\hat{\varrho}$. Because of the practical limitations imposed by the hardware, measurements of properties of interest might be costly, or not technically possible. It is then highly desirable to be able to reconstruct the quantum state $\hat{\varrho}$ from simple, experimentally feasible measurements.

The traditional approach for reconstructing a quantum state from measurement data is called quantum state tomography (QST) Vogel and Risken (1989); Ježek et al. (2003); Banaszek et al. (2013). A typical procedure consists of maximum-likelihood reconstruction of a density operator parametrized as $\hat{\rho}\propto\hat{T}^{\dagger}\hat{T}$ James et al. (2001), where $\hat{T}$ is a tri-diagonal hermitian matrix, enforcing the positive semi-definite requirement on $\hat{\rho}$. Such procedures assume no a priori phase structure to the quantum state, or even whether it is necessarily pure. Such “full” QST therefore typically scales exponentially. Given this, full QST can only be effectively carried out for systems with a relatively small number of particles or qubits Häffner et al. (2005). In general however, physical quantum states – such as ground states of local Hamiltonians – possess large degrees of structure. This often makes it possible to obtain a compact representation with resources scaling polynomially with the size of the Hilbert space. The most notable example are matrix product states (MPS), which have been used to successfully reconstruct quantum states outside the reach of full QST Cramer et al. (2010); Lanyon et al. (2017). However, so-called MPS tomography inherits the intrinsic limitations of the MPS representation, namely the restriction to one dimensional systems and low-entangled states, which limits the reconstruction of short-time dynamics, for example. The inherent structure of a quantum state can also be exploited in alternative ways, such as in permutationally invariant QST Tóth et al. (2010); Moroder et al. (2012) or compressed sensing Gross et al. (2010).

In this section, we overview a machine learning-based approach to QST, and show that unsupervised learning of generative models provides a very natural framework for reconstructing quantum many-body states. As described in the last section, RBMs offer a generative modeling framework that is conceptually interpretable in the context of statistical physics. In addition, they have been more widely explored in applications in classical and quantum state reconstruction than any other generative model. We start by considering the simplest case of reconstructing a thermal state in the classical limit, and proceed with increasing complexity to the case of pure quantum wavefunctions and finally density operators.

III.1 Classical limit

We start with the reconstruction of a physical system at thermal equilibrium, and consider the classical limit where the Hamiltonian under consideration is diagonal in the measurement basis $\{|\bm{\sigma}\rangle\}$ . The density operator we aim to reconstruct simply reduces to

[TABLE]

where $P_{\beta}(\bm{\sigma})=e^{-\beta H(\bm{\sigma}^{z})}/Z_{\beta}$ is the classical Boltzmann distribution in the canonical ensemble and $Z_{\beta}$ its partition function. State reconstruction is inherently a classical problem here, corresponding to the unsupervised learning of the distribution $P_{\beta}(\bm{\sigma})$. A simple yet non-trivial example is given by the Ising model, where $N$ spins interact with Hamiltonian

[TABLE]

with the sum running over nearest neighbours on a lattice. In two dimensions, the spin system displays ferromagnetic order at low temperature and a high-temperature disordered state, separated by a continuous phase transition.

As first demonstrated in Ref. Torlai and Melko (2016), different RBMs can be trained on datasets containing spin configurations at different temperatures across the phase diagram, generated by importance sampling the partition functions using Monte Carlo simulations Hastings (1970). The quality of the reconstruction can be assessed by comparing expectation values of thermodynamics observables generated by the RBM with the exact values calculated on the datasets. In Fig. 2 we report such a comparison for the magnetization and specific heat, with a varying number of hidden units in the RBM. While the magnetization converges very quickly – since it it explicitly encoded in the dataset – a larger number of hidden units is required to accurately reproduce the specific heat, particularly in the presence of large fluctuations at the critical point. Finally, we point out the curious observation that the quality of the reconstruction does not obviously improve for deep versions of the RBMs Morningstar and Melko (2017), such as deep belief networks Hinton et al. (2006) or deep Boltzmann machines Salakhutdinov and Hinton (2009).

III.2 Positive wavefunctions

We now turn to quantum states described by pure density operators $\hat{\varrho}=|\Psi\rangle\langle\Psi|$, where the wavefunction has representation $|\Psi\rangle=\sum_{\bm{\sigma}}\Psi(\bm{\sigma})|\bm{\sigma}\rangle$ with coefficients $\Psi(\bm{\sigma})=\langle\bm{\sigma}|\Psi\rangle$ in the measurement basis $\{|\bm{\sigma}\rangle\}$. In addition, we assume for now that the pure state $|\Psi\rangle$ has a real and positive representation in this basis, $\Psi(\bm{\sigma})\in\mathbb{R}$ and $\Psi(\bm{\sigma})>0\>\>\forall|\bm{\sigma}\rangle$. Under this assumption, valid for example for ground states of so-called “stoquastic" Hamiltonians Bravyi et al. (2008), the wavefunction $|\Psi\rangle$ is uniquely characterized by the probability distribution underlying a set of projective measurements, given by the Born rule $P(\bm{\sigma})=|\Psi(\bm{\sigma})|^{2}$. The inherently probabilistic nature of quantum mechanics provides a simple and natural way to define a representation of a pure and positive quantum state in terms of an RBM Torlai et al. (2018),

[TABLE]

Note that, since RBMs are universal approximators of any discrete probability distribution, provided the number of hidden units in the network is sufficiently large, the RBM wavefunction $\psi_{\bm{\lambda}}(\bm{\sigma})$ is capable of representing any positive quantum state to arbitrary accuracy.

Because of the positivity of the target state, quantum state reconstruction in this case is equivalent to conventional RBM unsupervised learning. Upon minimization of the KL divergence between the projective measurement distribution and the RBM distribution,

[TABLE]

the RBM wavefunction approximates the target state $\psi_{\bm{\lambda}}\sim\Psi$ as desired.

III.2.1 Measurement of physical observables

By discovering a set of parameters that successfully minimizes the cost function, the RBM builds an internal representation of the unknown target wavefunction and can be sampled to compute expectation values of physical observables. If the observable $\hat{\mathcal{O}}$ is diagonal in the measurement basis, its expectation value reduces to a thermal average with respect to the RBM distribution

[TABLE]

which can be approximated by a Monte Carlo average using block Gibbs sampling. Calculations of diagonal observables provide a direct verification of the quality of the training, since the expectation values can be compared with those calculated directly on the training dataset.

More interestingly, the RBM allows one to estimate average values of observables which are off-diagonal in the measurement basis. For this case, the expectation value reduces to the average $\langle\hat{\mathcal{O}}\rangle=\langle\mathcal{O}_{L}(\bm{\sigma})\rangle_{p_{\bm{\lambda}}(\bm{\sigma})}$, where

[TABLE]

is the so-called local estimate of the observables Becca and Sorella (2017). Provided the matrix representation of $\hat{\mathcal{O}}$ is sufficiently sparse in the measurement basis (i.e. the number of off-diagonal elements that are non-zero scales polynomially with $N$), its expectation value can be efficiently estimated with Monte Carlo.

Another important quantity amenable to calculation with an RBM is the entanglement of a subsystem $A$, which for pure states is quantified by the Renyi entropy Rényi (1961)

[TABLE]

with $\hat{\rho}_{A}=\text{Tr}_{A^{\perp}}|\Psi\rangle\langle\Psi|$ the reduced density matrix of $A$. For the case of $\alpha=2$, the entanglement entropy can be calculated by considering two identical replicas of the original system, and computing the overlap between the states with and without the configurations of subregion $A$ swapped between the replicas

[TABLE]

Calculations of the entanglement entropy via this procedure has been successfully carried out in numerical simulations using different flavors of quantum Monte Carlo algorithms Hastings et al. (2010); Kallin et al. (2011); Zhang et al. (2011); Stéphan et al. (2013). In the experimental context, the measurement of the swap operator has been performed for ultracold bosons in optical lattices by means of quantum interference Islam et al. (2015). Interestingly, the above implies that same measurement can be efficiently implemented using the RBM wavefunction, where instead of replicating the experimental NISQ hardware in the laboratory, one can first train the RBM to learn the experimental wavefunction, and then perform the calculation of the swap operator by replicating the neural network Torlai et al. (2018).

III.2.2 Reconstructing quantum spins on a lattice

As an example, we review a numerical experiment for the quantum reconstruction of the ground state of the transverse-field Ising model, with Hamiltonian

[TABLE]

This spin system undergoes a quantum phase transition between a ferromagnetic state for a small value of the transverse field $h$, and a paramagnetic state for large $h$. Measurement data in the $\{|\bm{\sigma}^{z}\rangle\}$ basis can be generated with standard methods White (1992); Evertz (2003). Similar to the case of the classical Ising model above, different RBMs are trained at different values of the transverse field, and then sampled to generate expectation values of observables Torlai et al. (2018). Fig. 3 shows the reconstruction of the average diagonal and off-diagonal (transverse) magnetizations for the quantum Ising model on a square lattice, and the entanglement entropy for the one-dimensional chain, calculated using the swap operator between replicated copies of the neural network.

III.3 Complex wavefunctions

The assumption of a pure and positive quantum state enables RBM reconstruction with a favorable scaling with respect to the number of particles in the system. In general however, experimental quantum states might violate this assumption, containing a sign or a phase structure, where the coefficients of the wavefunction can be both positive and negative, or complex-valued $\Psi(\bm{\sigma})=|\Psi(\bm{\sigma})|e^{i\phi(\bm{\sigma})}$. A sign structure often appears in ground states of non-stoquastic Hamiltonians, such as quantum spins with competing interactions on frustrated lattices, or fermions. In this case, data from a single measurement basis is clearly not sufficient to fully capture the quantum state, since the corresponding probability distribution $P(\bm{\sigma})=|\Psi(\bm{\sigma})|^{2}$ does not contain any fingerprints of the sign structure. Thus, reconstruction of the quantum state requires measurement in additional bases.

The first step for generalizing the RBM reconstruction to complex-valued wavefunctions is to define an appropriate neural-network parametrization of the quantum state. The most straightforward way consists of adding a phase factor to the positive RBM wavefunction defined in the previous section, $\psi_{\bm{\lambda\mu}}(\bm{\sigma})=\sqrt{p_{\bm{\lambda}}(\bm{\sigma})}e^{i\theta_{\bm{\mu}}(\bm{\sigma})}$. There is a large amount of freedom in choosing the nature of the phase function $\theta_{\bm{\mu}}(\bm{\sigma})$ in term of additional network parameters $\bm{\mu}$, and it needs not to be restricted to generative models. In fact, any feedforward neural network, such as convolutional networks Krizhevsky et al. (2012), could be used to this end. Another powerful way to adapt the RBM to quantum states is by using complex-valued weights and biases Carleo and Troyer (2017). In this review we will use an additional RBM to capture the phases, leading to the following neural-network wavefunction Torlai et al. (2018)

[TABLE]

Note that the generation of configurations in the reference basis corresponds to sampling the distribution $|\psi_{\bm{\lambda}}(\bm{\sigma})|^{2}=p_{\bm{\lambda}}(\bm{\sigma})$ which does not depend on the phases and can be then carried out by using block Gibbs sampling on the RBM with parameters $\bm{\lambda}$.

III.3.1 Learning the phase structure

The reconstruction of a phase structure requires performing additional measurements in bases different than the reference one where the RBM wavefunction is expressed. This involves applying a unitary transformation $\hat{\mathcal{U}}$ to the quantum state,

[TABLE]

where $|\bm{\sigma}^{\bm{b}}\rangle=|\sigma_{1}^{b_{1}},\dots,\sigma_{N}^{b_{N}}\rangle$ and $b_{j}$ identifies a particular choice of local basis for the $j$-th degree of freedom. The corresponding probability distribution after the measurement, $P(\bm{\sigma}^{\bm{b}})=|\Psi(\bm{\sigma}^{\bm{b}})|^{2}$, contains partial information on the phases and can be used to reconstruct the complex state. In general, such a unitary transformation consists of a collection of independent rotations of the local Hilbert spaces. The number and the type of rotations required to extract sufficient information to learn a phase depends on the structure of the specific quantum state under reconstruction.

Given a dataset $\mathcal{D}=\{\bm{\sigma}^{\bm{b}}\}$ of measurements in different bases, the RBM reconstruction can be realized by minimizing the total KL divergence in all bases,

[TABLE]

where we have omitted the constant entropy term. By taking the gradients with respect to the parameters one obtains:

[TABLE]

[TABLE]

where the averages over the quasi-probability distribution $Q^{\bm{b}}(\bm{\sigma})=\mathcal{U}_{\bm{\sigma}^{\bm{b}},\bm{\sigma}}\psi_{\bm{\lambda\mu}}(\bm{\sigma})$ are calculated directly on the measurement data. Since the negative phase does not depend on the phase parameters $\mu$, standard CD training can be directly applied here. A detailed derivation of the gradients can be found in Ref. Giacomo Torlai (2018).

III.4 Density operators

When the purity of the quantum state of interest cannot be assumed, one needs to reconstruct the full density operator, $\varrho(\bm{\sigma},\bm{\sigma}^{\prime})=\langle\bm{\sigma}|\hat{\varrho}|\bm{\sigma}^{\prime}\rangle$. Similar to the case of a pure state, before handling the reconstruction we require a representation of the density matrix in terms of a set of network parameters, $\rho_{\bm{\lambda\mu}}(\bm{\sigma},\bm{\sigma}^{\prime})$, i.e. a neural density operator (NDO). However, in contrast with an RBM wavefunction, the construction of a NDO has more stringent requirements, namely the Hermitian condition $\hat{\rho}_{\bm{\lambda\mu}}=\hat{\rho}^{\dagger}_{\bm{\lambda\mu}}$ and the positive semi-definite condition $\hat{\rho}_{\bm{\lambda\mu}}\geq 0$. One way to enforce the latter directly into the neural network representation consists of adding a set of auxiliary degrees of freedom that purifies the mixed state of the physical system Torlai and Melko (2018).

III.4.1 Latent space purification

For any mixed quantum state, it is always possible to introduce a set of variables $\bm{\alpha}$ in such a way that the quantum state of the composite system is pure Benenti et al. (2004). In the context of neural networks, we can introduce a RBM wavefunction for the enlarged system

[TABLE]

and obtain a NDO by tracing out the auxiliary variables:

[TABLE]

By embedding the auxiliary units in the latent space of the neural network, it is possible to perform this trace analytically Torlai and Melko (2018)

[TABLE]

Here we defined

[TABLE]

and

[TABLE]

capturing, respectively, the correlations within the system, and the correlations between the system and the environment. The new parameters $\bm{V}_{\bm{\lambda}/\bm{\mu}}$ encode the degree of mixing of the state of the physical system – they are identically zero for a pure state.

The cost function for the quantum reconstruction of a NDO is given by

[TABLE]

and its gradients can be easily calculated analytically Giacomo Torlai (2018); Torlai and Melko (2018). Similarly to the case of a pure state, all the gradients can be evaluated directly on the training data (provided the appropriate unitary rotations are applied to the state). The exception is the term involving the partition function (the negative phase), which is approximated by the CD algorithm using a finite step of block Gibbs sampling (equivalent to sampling the distribution $\rho_{\bm{\lambda\mu}}(\bm{\sigma},\bm{\sigma})$). Given that the purification through the latent space of a RBM architecture generates a physical density operator ($\hat{\rho}_{\bm{\lambda\mu}}\geq 0$), this type of ansatz is also suitable for the simulation of quantum dynamics of open systems, which was recently explored in various numerical experiments Hartmann and Carleo (2019); Vicentini et al. (2019); Nagy and Savona (2019).

When evaluating the gradients of the cost function in Eq. 38, the NDO needs to be transformed back into the reference basis by appropriate unitary transformations related to the measurement basis, $\hat{\rho}^{\bm{b}}_{\bm{\lambda\mu}}=\hat{\mathcal{U}}_{\bm{b}}\hat{\rho}_{\bm{\lambda\mu}}\hat{\mathcal{U}}^{\dagger}_{\bm{b}}$. This rotation has to be carried out explicitly and it is then only feasible as long as $\hat{\mathcal{U}}$ acts non-trivially on a sufficiently small number of degrees of freedom. This limitation can be circumvented by avoiding the parametrization of the quantum state directly, and using instead a generative model to represent the probability distribution underlying the measurement outcomes of a set of information complete set of positive-operator valued measures (POVM) Carrasquilla et al. (2019).

III.5 Reconstruction of experimental wavefunctions

We have shown that RBMs trained with unsupervised learning offer a versatile approach to quantum state reconstruction of many-body systems. In this Section, we turn to the case of RBM reconstruction of experimental data from NISQ hardware.

III.5.1 Noise mitigation

One of the major obstacles in reconstructing quantum states from real experiments is the presence of measurement errors. In practice, when performing measurements on a system prepared in the quantum state $\hat{\varrho}$, one obtains measurement outcomes $\bm{\tau}$ which do not correspond to projective measurements, but are instead described by a POVM $\hat{\Pi}^{(\bm{\tau})}=\sum_{\bm{\sigma}}p(\bm{\tau}\>|\>\bm{\sigma})|\bm{\sigma}\rangle\langle\bm{\sigma}|$, where the distribution $p(\bm{\tau}\>|\>\bm{\sigma})$ is the probability of recording the outcome $|\bm{\tau}\rangle$ given the actual measurement $|\bm{\sigma}\rangle$. The probability distribution underlying a set of measurement data is then given by $P(\bm{\tau})=\text{Tr}[\Pi^{(\bm{\tau})}\hat{\varrho}]$. Assuming the rate of measurement errors $p(\bm{\tau}\>|\>\bm{\sigma})$ are known from the experiment, it is possible to incorporate the noisy measurements in the RBM architecture in such a way that the neural network learns the de-noised distribution, corresponding to the projective measurements prior to the noisy measurement process Torlai et al. (2019).

We show this so-called noise regularization for the case of a pure and positive RBM wavefunction $\psi_{\bm{\lambda}}(\bm{\sigma})=\sqrt{p_{\bm{\lambda}}(\bm{\sigma})}$, where $p_{\bm{\lambda}}(\bm{\sigma})$ now identifies the probability distribution over the de-noised visible layer $|\bm{\sigma}\rangle$ (before imperfect measurements). Rather than feeding the data directly to this layer, the RBM can be enlarged by adding a third (noise) layer encoding the variable $\bm{\tau}$, with underlying probability distribution $\widetilde{p}_{\bm{\lambda}}(\bm{\tau})=\sum_{\bm{\sigma}}p(\bm{\tau}\>|\>\bm{\sigma})p_{\bm{\lambda}}(\bm{\sigma})$. The weights connecting the noisy and the de-noised layer can be crafted in such a way that the resulting conditional distribution $p(\bm{\tau}\>|\>\bm{\sigma})$ matches the experimental values. The quantum reconstruction is performed by minimizing the negative log-likelihood with respect to the noisy data

[TABLE]

with the following update rule,

[TABLE]

The training is realized similarly to the case of ideal measurements, with the only difference being that the learning in the positive phase, instead of being calculated directly from the data, is driven by the Bayesian posterior distribution

[TABLE]

Following this training procedure, the trained RBM wavefunction $\psi_{\bm{\lambda}}(\bm{\sigma})=\sqrt{p_{\bm{\lambda}}(\bm{\sigma})}$ approximates the quantum state prior to the application of the measurement errors. A similar strategy based on autoencoders has also been put forward, and applied for the neural-network quantum reconstruction of experimental photonic states Macarone Palmieri et al. (2019).

III.5.2 Application to a Rydberg-atom quantum simulator

Finally, we summarize a recent experiment where RBM quantum reconstruction was applied to real data from a NISQ simulator. Specifically, the experimental system consists of an array of cold Rydberg atoms Bernien et al. (2017); Endres et al. (2016), one of the highest-quality platforms for programmable simulation of Ising-like quantum spins Schauß et al. (2015); Labuhn et al. (2016); Guardado-Sanchez et al. (2018). In the experiment, 87Rb atoms are individually trapped by optical tweezers in a defect-free array. The atomic ground state $|g\rangle$ is coupled to an highly excited Rydberg state $|r\rangle$ by a uniform laser drive, and the atoms interact with a Van der Waals potential, resulting into the Hamiltonian

[TABLE]

Here, $\Omega$ is the Rabi frequency, $V_{nn}$ is the nearest-neighbour interaction, $\hat{n}_{i}$ is the occupation number operator for the Rydberg state and $\Delta$ is the laser detuning. For small “transverse field” $\Omega$ and large negative detuning $\Delta$, the ground state of $\hat{H}$ is approximately a product state with all atoms in the ground state $|g\rangle$. Conversely, when the detuning is large and positive, depending on the spacing between the atoms, the simulator can be driven into different broken-symmetry states Bernien et al. (2017).

After initializing the simulator in the fiducial state where all atoms are prepared in $|g\rangle$, the ground state of Hamiltonian (42) in the broken-symmetry phase is obtained by an adiabatic ramp in the laser detuning $\Delta$. Assuming perfect adiabatic evolution, the instantaneous state of the system can always be gauged to be positive. Therefore, under this assumption, measurements in the occupation number basis $|\bm{\sigma}\rangle$ (i.e. eigenstates of $\hat{n}_{j}$) are sufficient to reconstruct the quantum state. The real experimental state is generally expected to be described by a mixed density operator, due to the unavoidable decoherence. Nevertheless, if the purity of the state remains sufficiently high, a pure-state RBM approximation is capable of correctly capturing properties defined over local subsystems Torlai et al. (2019). As described in the previous Section, measurement errors are mitigated using a noise layer regularization.

The experiment is performed with $N=8$ atoms for the transition into the $\mathbb{Z}_{2}$ ordered phase. The adiabatic sweep is halted at subsequent times (corresponding to an increase laser detuning $\Delta$), where a collection of about 3,000 measurements is taken. Each measurement consists of a bit-string $\bm{\tau}$ recording whether each single atom was measured in the ground or in the Rydberg state. Each dataset is then independently input to an RBM, and used to discover an optimal pure state approximation by incorporating the measurement errors with the noise layer. Once trained, directly sampling the RBM gives access to observable diagonal in the measurement basis. As an example, we show in Fig. 4a the average Rydberg population at different detunings during the sweep. We can compare the values generated by the RBM results with the experimental data (providing direct verification) and with full simulations of the Limbladian master equation. The RBM can then be used to sample expectation values of observables not accessible in the experimental apparatus. In particular, we show the average transverse field (Fig. 4b) and the Renyi mutual information $I_{2}(s)=S_{2}(\hat{\rho}_{s}^{A})+S_{2}(\hat{\rho}_{s}^{B})-S_{2}(\hat{\rho})$ (Fig. 4c), showing an overall good agreement with the results from simulations of open system dynamics.

IV Conclusions and outlook

Modern machine learning has provided us with generative modeling techniques that are perfectly suited for the emerging landscape of NISQ hardware. Stochastic neural networks, such as RBMs and their cousins Carrasquilla et al. (2019); Rocchetto et al. (2018); Macarone Palmieri et al. (2019), are heuristically known to provide good quality state reconstructions for intermediate-scale and noisy data. In that sense, their adoption to quantum state reconstruction on devices of tens, hundreds, or even thousands of qubits should come as no surprise.

As discussed in this review, the systematic development of RBM theory for use in quantum state reconstruction is becoming well understood from a formal standpoint. Parallel to theoretical and algorithmic advancements, a crucial role is also played by the development of related open source software Beach et al. (2018); Carleo et al. (2019), easily accessible to experimentalists. However, many fundamental questions still remain to be answered if such machine learning techniques are to become fully integrated with NISQ hardware.

First, as evident in this review, the most well-studied cases involve wavefunctions that are real and positive – mathematically equivalent to probability distributions. There, standard generative models such as RBMs can be employed with little alteration from their original industry-motivated design. Under the assumption of purity, recent work has demonstrated the efficiency of modern algorithms for unsupervised learning in approximate state reconstruction. In this case, RBMs in particular have shown their utility in producing accurate and scalable estimators for physical observables not directly available from the original data set; i.e. they generalize well. Of particular interest is the basis-independent Renyi entanglement entropy, which can be measured directly from a trained RBM using a scalable algorithm involving replication of the model wavefunction. This is perhaps the most striking example of a measurement that is resource-intensive experimentally Islam et al. (2015), but relatively simple to implement in the trained generative model.

Real and positive wavefunctions occupy a special place in the landscape of physically-interesting states; for example, they are the ground states of stoquastic Hamiltonians. However, a large proportion of quantum states under study theoretically and experimentally cannot be assumed to have this significant simplification. As we have discussed in detail, in the case of complex wavefunctions, state reconstruction is possible with RBMs (and other generative models). What is required first is a convention to parameterize the phase, e.g. in additional hidden layers, or as complex weights Carleo and Troyer (2017). Then, measurements in more than one basis are needed to train the parameters encoding the phase of the wavefunction. Given this strategy, experimental NISQ wavefunctions, such as cold-atom implementations of the fermionic Hubbard model Mazurenko et al. (2017) or other interesting many-body Hamiltonians, may conceivably be reconstructed in the near future.

Herein lies one frontier for state reconstruction with machine learning. In the quest to construct a NISQ-compatible generative modeling method, foremost is the question of scaling of the number of measurement bases required for informational completeness. Very little is known theoretically about this scaling for wavefunctions of interest to NISQ simulators; in the pure case, the number of basis required to learn a $N$-qubit state could range from 1 (see above), to a number that grows exponentially in $N$ (see e.g. Ref. Ma et al. (2016)). This wide range of possibilities leaves open many questions about the learnability of quantum states. For example, for what other typical physical wavefunctions is the number of bases tractable in the context of generative models (RBMs or otherwise)? Also, how does the target wavefunction structure affect the number of measurements required in each basis? Finally, what is the relationship between these numbers and the scaling of the RBM parameters required for a desired representational accuracy? An entire field related to the study of how efficient learning relates to the sign or entanglement structure of a quantum state still lies ahead.

Moving away from pure states, the ability to represent density matrices suggests the possibility that machine-learning reconstruction can be expressed as an approximate re-formulation of more traditional quantum state tomography. The same scaling questions apply as above for the context of generic complex wavefunctions (albeit with the possibility of significant further roadblocks to scaling). A reformulation of the problem in the language of informationally-complete POVMs, briefly mentioned here Carrasquilla et al. (2019), offers the tantalizing possibility of scaling improvements, at the cost of (experimentally) more complicated measurements. Finally, success with the mixed-state density matrix formulation suggests that a re-imagination of process tomography as an unsupervised learning problem could also be in store. With further development along these lines, generative modeling is poised to breach beyond the realm of NISQ simulators, to become a tool for gate-based architectures in the near future.

Looking forward, it is clear that today’s hardware is a necessary stepping stone to the more powerful quantum technologies of the future. As these devices continue to grow, they will develop in lock-step with powerful classical algorithms, to aid in all stages of state preparation, measurement, verification, error correction, and more. With the dawn of artificial intelligence as the most powerful classical computing paradigm of a generation, it stands to reason that machine learning of quantum many-body states will play a critical role in the NISQ era and beyond.

Acknowledgements

The Flatiron Institute is supported by the Simons Foundation. R.G.M. is supported by NSERC of Canada, a Canada Research Chair, and the Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported through Industry Canada and by the Province of Ontario through the Ministry of Research & Innovation.