\section{Introduction}

Quantum information processing promises to dramatically advance physics and chemistry by providing efficient simulators for the Schr\"odinger or Dirac equations~\cite{lloyd_universal_1996,alan_qchem_2005,gerritsma_diracsim_2010}.
This is important because conventional methods are inefficient, scaling exponentially in the number of interacting subsystems.  Consequently, quantum simulations beyond a few tens of interacting particles are generally believed to be beyond the limitations of conventional supercomputers.   This inability to simulate large quantum systems means that important questions in condensed matter, such as the shape of the phase diagram for the Fermi--Hubbard model, remain open.
Analog quantum simulation raises the possibility that
special purpose \emph{analog devices} may be able to address such problems using current or near--future hardware~\cite{simon_simulation_2011,britton_simulation_2012,kim_simulation_2010}.  
A major objection to this avenue of inquiry is that analog simulators are not necessarily trustworthy~\cite{hauke_trust_2012,gogolin_boson_2013} and certification of them is not known to be efficient.
Without such certification, an analog simulator can at best only provide  hints about the answer to a given computational question.  A resolution to this problem is therefore essential if analog quantum simulators are to compete on an even footing with classical supercomputers.

An important first step towards a resolution is provided in~\cite{daSilva_practical_2011}, where it is shown that quantum systems with local time--independent Hamiltonians can be efficiently characterized given ensemble readout.  However, the method is not generally applicable, can be expensive and is not known to be either error robust or stable in cases where single shot measurements are used.
A number of machine learning and statistical inference methods  \cite{hentschel_machine_2010,hentschel_efficient_2011,sergeevich_characterization_2011,ferrie_how_2012,sergeevich_optimizing_2012,granade_robust_2012,lovett_differential_2013,svore_faster_2013} have been recently introduced to address similar problems in metrology or Hamiltonian learning.  In the context of Hamiltonian learning, such ideas have are known to be error--robust and lead to substantial reductions in the cost of high--precision Hamiltonian inference~\cite{granade_robust_2012}, albeit at the price of sacrificing the efficient scaling exhibited by~\cite{daSilva_practical_2011}.

We overcome these challenges by providing a robust method that can be used to characterize unknown Hamiltonians by unifying statistical inference with quantum simulation.  The key insight behind this is that Bayesian inference reduces the problem of Hamiltonian estimation to a problem in Hamiltonian simulation that can be efficiently solved using a trusted quantum simulator.  Our algorithm achieves this through the following steps. We begin by positing a Hamiltonian model for the system and a probability distribution over the parameters of the Hamiltonian model.  We then use a novel guess heuristic for the optimal experiment that adaptively chooses experiments based on the current uncertainty in the Hamiltonian.  The experiment is then performed and the trusted quantum simulator is used to efficiently compute the likelihood of the measurement outcome occurring if each hypothetical model were true.  These likelihoods are then used by the algorithm to update its knowledge of the Hamiltonian parameter via Bayes rule, resulting in an updated probability distribution, called the posterior distribution. This process is then repeated until the uncertainty in the unknown Hamiltonian parameters (as measured by the posterior variance) becomes sufficiently small.  This iterative process is depicted in \fig{flowchart}.








To make the problem concrete, we represent each hypothetical Hamiltonian $H_j$ by a vector of real numbers $\vec{x}_j \in \mathbb{R}^d$ such that $H_j = H(\vec{x}_j)$.  The Hamiltonian model is therefore specified by $H(\vec{x})$. 





\begin{figure}[t!]
\centering
\includegraphics[height=2.75in]{flowchart2.pdf}
\caption{\label{fig:flowchart} Flowchart for Hamiltonian learning algorithm.}
\end{figure}


We consider three classes of experiments that can be performed to infer the Hamiltonian, $H$, given an initial state $\ket{\psi}$ (typically a  pseudorandom state~\cite{emerson_pseudo-random_2003}): (a) Classical Likelihood Evaluation (CLE), (b) Quantum Likelihood Evaluation (QLE) and (c) Interactive Quantum Likelihood Evaluation (IQLE).  CLE is  the simplest of these experiments and is discussed in detail in~\cite{granade_robust_2012}.  It involves simply picking an experimental time $t$, and computing the likelihood $\Pr(D|\vec{x}_i)=|\bra{D}e^{-iH(\vec{x}_i)t}\ket{\psi}|^2$ using a classical computer, where $\vec{x}_i$ is a given set of Hamiltonian parameters and $D$ is the experimental outcome.  This function, known as the likelihood function, will not generally be efficiently computable on a classical computer because it involves quantum simulation.

In QLE experiments, a trusted quantum simulator is used to ameliorate these problems.  It does so by estimating $\Pr(D|\vec{x}_i)$ to be the fraction of times outcome $D$ occurs in a sufficiently large set of simulated experiments, which is efficient if $\Pr(D|\vec{x}_i)$ is only polynomially small.  
  This approach allows a complex quantum simulator, such as a fault tolerant quantum computer, to act as a certifier for an analog quantum simulator.  A trusted quantum simulator could also be constructed using a bootstrapping protocol wherein a smaller trusted analog simulator is the certifier.  This is possible if a compressed simulation scheme~\cite{kraus_compressed_2011} for the dynamics exists.

The Loschmidt echo famously shows that, for complex quantum systems, two nearly identical Hamiltonians will typically generate evolutions that diverge exponentially after a short time, before saturating at an exponentially small overlap~\cite{Haa06}.  This means that QLE will often be restricted to short evolution times  to guarantee efficiency (which is undesirable~\cite{granade_robust_2012}).  We resolve this by using IQLE experiments, which are described in~\fig{models}.  These experiments are reminiscent of the Hahn echo experiments commonly used in magnetic resonance and experimental quantum information processing~\cite{hahn_spin_1950}.
An IQLE experiment swaps the state of the unknown quantum system with that of a trusted quantum simulator then inverts the evolution based on a guessed Hamiltonian $H_{-}$.  The measurement in IQLE is always assumed to be in an orthonormal basis that has $\ket{\psi}$ as an element.  This produces $\Pr(D|\vec{x}_i)=|\bra{D}e^{iH_{-} t}e^{-iH(\vec{x}_i)t}\ket{\psi}|^2$




\begin{figure}
    \[
      \xymatrix @*=<0em> @R 1em @C 1.5em {
                        &                                      &                &                          & \ustick{t, \vec{x}_{-}}                 &        &                &  \\
                        & \qw                                  & *=<0em>{\times} \qw \qwx[1] & \qw & \gate{e^{-i H_{-} t}} \cwx[-1] & *=<1.8em,1.4em>{\xy ="j","j"-<.778em,.322em>;{"j"+<.778em,-.322em> \ellipse ur,_{}},"j"-<0em,.4em>;p+<.5em,.9em> **\dir{-},"j"+<2.2em,2.2em>*{},"j"-<2.2em,2.2em>*{} \endxy} \POS ="i","i"+UR;"i"+UL **\dir{-};"i"+DL **\dir{-};"i"+DR **\dir{-};"i"+UR **\dir{-},"i" \qw & \rstick{D} \cw &  \\
    \lstick{\ket{\psi}} & \gate{e^{-i H(\vec{x}_0) t}} \cwx[1] & *=<0em>{\times} \qw         & \qw                      & \qw                                   & \qw    & \qw            &  \\
                        & \dstick{t}                           &                &                          &                                       &        &                &
      }
    \]
    \caption{\label{fig:models}IQLE. The upper register is the trusted simulator, while the lower register is the system under study.  QLE is conceptually similar but with $H_{-}=0$.
\end{figure}




Although the Loschmidt echo may also seem to be problematic for IQLE experiments, we exploit it in our guess heuristic for $H_{-}$.  We call this heuristic the ``particle guess heuristic'' (PGH), which chooses $H_{-}:= H(\vec{x}_{-})$ by sampling $\vec{x}_{-}$ from the prior probability $\Pr({\vec{x}})$, which describes our current knowledge of the Hamiltonian parameters.  The set of parameters, $\vec{x}_{-}$, is called a \emph{particle} because it is described by a Dirac-delta distribution over parameter space.  The time $t$ is chosen by drawing a second particle $\vec{x}_{-}' \ne \vec{x}_{-}$ and setting $t=1/\|\vec{x}_{-}'-\vec{x}_{-}\|_2$.  As the uncertainty in the estimated parameter shrinks, the PGH adaptively picks longer times to ensure that informative experiments continue to be chosen as certainty about the unknown parameters increases.  The PGH also causes $e^{-iH(\vec{x}_i)t}$ to result in substantially different likelihoods for $\vec{x}_i$ that are within one standard deviation of the prior's mean, which we show in the appendix is optimal for certain learning problems.  

IQLE experiments with two outcomes also ensure that $\Pr(D|\vec{x})$ will not be exponentially small (with high probability) for $H$ an affine transformation acting on $\vec{x}$, since
\begin{align}
|\bra{\psi}e^{iH_{-} t}e^{-iHt}\ket{\psi}|&\ge 1-2\|H-H_{-}\|_2t \nonumber\\
&\ge 1- O(\|\vec{x}-\vec{x}_{-}\|_2t).
\end{align}
If the prior distribution has converged to a unimodal distribution centered near the correct Hamiltonian (this is typical for Bayesian inference of non-degenerate learning problems~\cite{granade_robust_2012}) then $\|\vec{x}-\vec{x}_{-}\|_2\in \Theta(1/t)$.  This means that if we use a POVM with two elements: $\ket{\psi}\!\!\bra{\psi}$ and its orthogonal compliment $\openone-\ket{\psi}\!\!\bra{\psi}$ then we expect (a) neither probability will be exponentially small if $\vec{x}_{-}$ and $\vec{x_j}$ are near the mean and (b) $\Pr(\psi|\vec{x_j})$ will typically be exponentially small for $H(\vec{x_j})$ that differ substantially from the correct Hamiltonian.  The PGH therefore leads to IQLE experiments that rapidly eliminate incorrect hypotheses about the correct Hamiltonian.


The measurement outcomes yielded by the experiments are \emph{immediately processed} using Bayesian inference, as described in~\fig{flowchart}.  This immediate processing allows our algorithm to adaptively choose experiments based on its current knowledge of the correct Hamiltonian.
The state of knowledge is represented by a distribution that is called, previous to the next update step, the \emph{prior}.  In the cases we consider, the initial prior distribution before any data is observed is taken to be uniform.  This encodes a state of maximum ignorance about the correct $\vec{x}$.
The prior distribution is updated as measurement outcomes are recorded using Bayes' rule, which gives the proper way of computing the probability of each $\vec{x}_j$ being correct given the observed data and the prior.  It states that if datum $D$ is recorded then
\begin{equation}
\Pr(\vec{x}_j|D) \propto \Pr(D|\vec{x}_j)\Pr(\vec{x}_j),\label{eq:bayes}
\end{equation}
up to a normalization factor and
 $\Pr(\vec{x}_j|D)$ is called the posterior distribution.  

Eq.~\eq{bayes} can be efficiently computed (for a polynomial number of $H_j$) only if the likelihood function $\Pr(D|\vec{x}_j)$
is tractable.  QLE and IQLE experiments allows $\Pr(D|\vec{x}_j)$ to be efficiently estimated, which removes the main obstacle to using Bayesian methods to learn the correct $\vec{x}$.


\begin{figure}[t!]
\includegraphics{complete_scaling.pdf}
\caption{The quadratic loss  plotted as a function of the number of inversion experiments for Ising models on the complete graph. The shaded areas show a $50\%$ confidence interval for the quadratic loss.}\label{fig:complete_scale}
\end{figure}

A secondary problem is that \emph{exact} computations of the update rule are intractable in practice because an infinite number of Hamiltonians  could potentially describe the system; hence, a probability
distribution over Hamiltonians cannot be exactly represented on either a classical or quantum computer.  This problem can be addressed by using the sequential Monte Carlo (SMC) approximation \cite{huszar_adaptive_2011,granade_robust_2012,doucet_tutorial_2011}, which approximates the probability distribution using a weighted sum of particles (Dirac delta functions). Each particle corresponds to a particlar $\vec{x}_j$ and is a hypothesis about the correct Hamiltonian parameters $\vec{x}$.  SMC assigns a weight $w_j$ to each particle that represents the probability of that hypothesis.  The weights are normalized such that $\sum_{j}w_j =1$.  The update rule for the probability distribution under the SMC approximation then becomes
$w_j \mapsto \Pr(D|\vec{x}_j)w_j$, followed by normalization.
If necessary, a resampling step is used after updating to ensure that the inference procedure remains stable, as discussed in~\cite{granade_robust_2012} and in \app{bayes}.

The algorithm then iteratively updates the weights $w_i$ and positions $\vec{x}_i$ of the sequential Monte Carlo particles representing the distribution over Hamiltonians $\Pr(H | D)$, conditioned on the data recorded at each step. In this way, the full state of knowledge at each step is iteratively carried forward, and is used to heuristically design future experiments according to the PGH.  Subsequent updates will then refine this estimate of the unknown Hamiltonian parameter until the uncertainty of the estimated Hamiltonian is sufficiently small, as measured by the trace of the posterior covariance matrix.

\begin{figure}[t!]
\centering
\includegraphics{exp_scale.pdf}
\caption{The median decay exponent for the quadratic loss as a function of the number of parameters in the Ising model, $d$.}\label{fig:exp_scale}
\end{figure}

Now that we have discussed how our algorithm works, we will proceed to assess its cost.  We will show that the cost of Hamiltonian inference on a fixed number of IQLE experiments is exponentially smaller than the cost of using CLE.  This is significant because CLE  gives the best known methods for some problems~\cite{granade_robust_2012}.  

A natural measure of the cost is the number of quantum simulations needed to estimate the Hamiltonian parameters.  The total cost is therefore,
\begin{equation}
{\rm{Cost}} = N_{\rm steps}(\delta) \times {\rm{Cost}(\rm update;\epsilon)}.\label{eq:costbasic}
\end{equation}
Here $N_{\rm steps}$ is the number of updates needed to make the uncertainty less than $\delta$ and ${\rm{Cost}(\rm update;\epsilon)}$ is the number of samples from the trusted simulator that are needed to update the particle weights using Eq.~\eq{bayes} within error $\epsilon$ in the $1$--norm.
We show in \app{stability} that, with high probability,
${\rm{Cost}(\rm update;\epsilon)}$ scales as
\begin{equation*}
\frac{|\{\vec{x}_i\}|}{\epsilon^2}\left(\mathbb{E}_{D|H}\left[\frac{\max_k \Pr(D|\vec{x}_k)(1-\Pr(D|\vec{x}_k))}{\left( {\sum_k \Pr(D|\vec{x}_k)\Pr(\vec{x}_k)}\right)^2}\right]\right).
\end{equation*}
This implies that the update process will be efficient if the number of particles required is small and the resultant probability distribution is not too flat.  That is, $|\{\vec{x}_i\}|\in O({\rm poly}(n))$ and  $\sum_k \Pr(D|H(\vec{x}_k))\Pr(H(\vec{x}_k))\in O(1/{\rm poly}(n))$, where $n$ is the number of interacting systems.  It has been shown that SMC algorithms require a number of particles that scales sub--exponentially in $d$~\cite{beskos_stability_2011}, which itself may not be a function of $n$.  This means that in practice, a small number of particles will typically be required.  The robustness of the algorithm to sampling errors is discussed in~\cite{ferrie_likelihood-free_2013} as well as in \app{stability}, so relatively large $\epsilon$ can be tolerated.








If the posterior distribution has converged to a unimodal distribution such that $\vec{x}$ is within a fixed distance from the mean, then the PGH and~\eq{bayes} ensure that $\mathbb{E}_{H_{-}} [|\!\bra{\psi}e^{iH_{-} t}e^{-iHt}\ket{\psi}|^2] \in \Theta(1)$ since $t\in \Theta(|\vec{x}-\vec{x}_{-}|^{-1})$.  If a two outcome measurement is used then  Markov's inequality implies that  $\sum_k \Pr(D|H(\vec{x}_k))\Pr(H(\vec{x}_k))\in \Theta(1)$ with high probability. 
By a trivial generalization of this argument, it is clear that a super--polynomial reduction in the cost of performing~\eq{bayes} relative to CLE is obtained with high probability for IQLE experiments if $d\in O({\rm poly}(n))$ and the  \emph{effective} number of outcomes, $\sum_j \Pr(j|\vec{x}_k)^{-2}$, is at most $O({\rm poly}(n))$ for each $\vec{x}_k$.

In contrast, QLE experiments may not lead to a super--polynomial separation in the cost estimates for generic Hamiltonians and large $t$ because $\sum_k \Pr(D|H(\vec{x}_k))\Pr(H(\vec{x}_k))\in 2^{-\Theta(n)}$ with high probability for complex quantum systems~\cite{ududec_equilibration_2013,Haa06}.  This  can be rectified by choosing small $t$ as per~\cite{daSilva_practical_2011}, but such QLE experiments will be much less informative~\cite{granade_robust_2012}.







If a fixed number of updates are required, then the previous discussion and~\eq{costbasic} suggest that IQLE will provide an exponential advantage over CLE.  If inference within a fixed error tolerance, $\delta$, is required then the cost estimate is much more challenging.
 Each two-outcome measurement yields at most one bit of information about $H$ per measurement hence $N_{\rm steps}(\delta) \in \Omega(d\log_2(1/\delta))$.  For most models of interest, $d$ is polynomial (or even constant) in $n$ and hence a small number of updates should typically suffice.  It is, however, unclear whether this lower bound is tight; hence, we turn to numerical evidence to show that our algorithm can efficiently learn Hamiltonians in certain cases.






Consider the problem of learning $H({\vec{x}})$ using IQLE experiments for an Ising model with no transverse field:
\begin{equation}\label{eq:Ising}
 H(\vec{x}) = \sum_{(i,j)\in G} x_{i,j}\ \sigma_z^{(i)}\sigma_z^{(j)},
\end{equation}
where $G$ is the edge set of an interaction graph on $n$ qubits. Unless otherwise specified, we take $x_{i,j}\in [-1/2,1/2]$ uniformly at random.  We take the initial state for the evolution to be $\ket{\psi}=\ket{+}^{\otimes n}$.  We choose this Hamiltonian not only because it is physically relevant~\cite{richerme_trapped-ion_2013}, but also for numerical expediency, since the learning process require the algorithm to perform thousands of simulated evolutions of the initial state.  All measurements are performed in the eigenbasis of $X^{\otimes n}$.  Restricting the measurements to two outcomes is unnecessary for these experiments because IQLE and the PGH concentrates $\Pr(D|\vec{x}_i)$ over a small number of outcomes for this Hamiltonian.

\fig{complete_scale} shows that the quadratic loss (a generalization of the mean--squared error for multiple parameters) shrinks exponentially with the number of experiments performed; however, the rate at which the error decreases slows as the number of qubits $n$ increases.  This is expected because $d=n(n-1)/2$ for the case of a complete interaction graph, which implies that the learning problem becomes more difficult as $n$ increases.  The data for interactions on the line is similar and is presented in \app{line} and QLE data is given in \app{QLE}.



\begin{figure}[t!]
\centering
\includegraphics[width=0.95\linewidth]{corner.pdf}
\caption{An approximate $1$ parameter Ising model on the complete graph.  The thin lines give the best fits to the exponential decays, which scale as $e^{-0.07N }$ and $e^{-0.23N}$ for $n=4$ and as $e^{-0.029N}$ and $e^{-0.23N}$ for $n=6$ qubits.\label{fig:corner}}
\end{figure}


The rate at which the learning process slows as $n$ increases is investigated in
\fig{exp_scale}.  We examine the slowing of the learning problem by fitting the quadratic loss, $\delta$, in each experiment to $Ae^{-\gamma N_{\rm steps}}$.  The median decay exponent, which is the median of the values of $\gamma$ attained for a set of experiments with constant $n$, measures how rapidly the algorithm learns the unknown parameters. \fig{exp_scale}  shows that these decay constants scale as $O(1/d)$ for the complete graph, and provides weaker evidence for the line.  This implies that $N_{\rm steps}(\delta)=O(d\log(1/\delta))$ for this Hamiltonian, which implies that the inference is efficient.  Similarly, the PGH implies that the total simulation time needed (for fixed $|\{\vec{x}_i\}|$) scales as $N_{\rm steps}\delta^{-1}\approx \delta^{-3}$, which is relevant in cases where the cost of a simulation is dominated by the evolution time.

Although $d=n(n-1)/2$ or $d=n-1$ in the examples considered above, $d$ can be approximately independent of $n$ in some cases.
An example of this behavior is given in \fig{corner}, where we consider the case where each of the $x_{i,j}$ is approximately the same value chosen uniformly on $[0,100]$, but with small normally distributed fluctuations with mean $0$ and variance $10^{-4}$.  This causes the learning problem to be effectively one--dimensional initially, and then transition to $d=n$ when the small fluctuations  need to be identified to learn the Hamiltonian parameters within a fixed accuracy.  The transition from a single-- to a multi--parameter learning problem happens at $\delta \approx d\times 10^{-4} \approx 10^{-3}$, which coincides with the point when the slope in \fig{corner} changes.  This emphasizes that the cost of Hamiltonian estimation using our method only implicitly depends on $n$ through $d$.  In fact, the difference in the observed scaling of $\gamma$ is approximately a factor of $2.5$, which is what would be expected if  $\gamma\propto 1/d$.



In conclusion, we have shown that Bayesian inference combined with the SMC approximation provides an ideal way to leverage a (potentially non--universal) quantum
simulator to characterize an unknown or unreliable quantum system.  We provide theoretical evidence that shows that the update rule, which is at the heart of the learning algorithm, can be performed efficiently using quantum resources.  We then illustrate the practicality of the algorithm and show that it is capable of learning unknown Ising couplings with surprisingly few experiments even in the presence of sampling errors.  We will show elsewhere that the algorithm is highly resilient to depolarizing noise and other forms of noise that can be introduced via a noisy swap gate



















\begin{acknowledgements}
  We acknowledge Troy Borneman for suggesting bootstrapping and also Krysta Svore and Allan Geller for useful feedback and discussion.  The numerical experiments performed here used SciPy, F2Py and QInfer~\cite{SciPy2001,peterson_f2py:_2009,ferrie_qinfer_2012}.  This work was supported by funding from USARO-DTO, NSERC, CERC and CIFAR.  CF was supported in part by NSF Grant Nos.~PHY-1212445 and PHY-1005540.
\end{acknowledgements}




