US11120359B2 - Phase estimation with randomized hamiltonians - Google Patents

Abstract

Existing methods for dynamical simulation of physical systems use either a deterministic or random selection of terms in the Hamiltonian. In this application, example approaches are disclosed where the Hamiltonian terms are randomized and the precision of the randomly drawn approximation is adapted as the required precision in phase estimation increases. This reduces both the number of quantum gates needed and in some cases reduces the number of quantum bits used in the simulation.

Description

CROSS REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional Application No. 62/819,301 entitled “PHASE ESTIMATION WITH RANDOMIZED HAMILTONIANS” and filed on Mar. 15, 2019, which is hereby incorporated herein by reference in its entirety.

FIELD

This application relates generally to quantum computing. In more detail, example approaches are disclosed where Hamiltonian terms are randomized and the precision of the randomly drawn approximation is adapted as the required precision in phase estimation increases.

SUMMARY

Existing methods for dynamical simulation of physical systems use either a deterministic or random selection of terms in the Hamiltonian. In this disclosure, example approaches are disclosed where the Hamiltonian terms are randomized and the precision of the randomly drawn approximation is adapted as the required precision in phase estimation increases. This reduces both the number of quantum gates needed and in some cases reduces the number of quantum bits used in the simulation.

Embodiments comprise randomizing phase estimation by replacing the Hamiltonian with a randomly generated one each time it is simulated. Further embodiments involve the use of randomization within an iterative phase estimation algorithm to select Hamiltonian terms for inclusion in the approximation as well as their ordering. Certain embodiments involve the use of importance functionals based on the significance of the term in the groundstate to determine whether it gets included in the randomly sampled Hamiltonian. Further embodiments involve the use of importance sampling based on the variational approximations to the groundstates, such as but not limited to, CISD states. Certain embodiments involve the use of adaptive Bayesian methods in concert with this process to quantify the precision of the Hamiltonian needed given the current uncertainty in the eigenvalue that the algorithm is estimating.

In this application, example methods for performing a quantum simulation using adaptive Hamiltonian randomization. The particular embodiments described should not be construed as limiting, as the disclosed method acts can be performed alone, in different orders, or at least partially simultaneously with one another. Further, any of the disclosed methods or method acts can be performed with any other methods or method acts disclosed herein. In particular embodiments, a Hamiltonian to be computed by the quantum computer device is inputted; a number of Hamiltonian terms in the Hamiltonian is reduced using randomization within a phase estimation algorithm; and a quantum circuit description for the Hamiltonian is output with the reduced number of Hamiltonian terms.

In certain embodiments, the reducing comprises selecting one or more random Hamiltonian terms based on an importance function; reweighting the selected random Hamiltonian terms based on an importance of each of the selected Hamiltonian random terms; and generating the quantum circuit description using the reweighted random terms. Some embodiments further comprise implementing, in the quantum computing device, a quantum circuit as described by the quantum circuit description; and measuring a quantum state of the quantum circuit. Still further embodiments comprise re-performing the method based on results from the measuring (e.g., using an iterative process). In some embodiments, the iterative process comprises computing a desired precision value for the Hamiltonian; computing a standard deviation for the Hamiltonian based on results from the implementing and measuring; and comparing the desired precision value to the standard deviation. Some embodiments further comprise changing an order of the Hamiltonian terms based on the reducing. Certain embodiments further comprise applying importance functions to terms of the Hamiltonian in a ground state; and selecting one or more random Hamiltonian terms based at least in part on the importance functions. Some embodiments comprise using importance sampling based on a variational approximation to a groundstate. Certain embodiments further comprise using adaptive Bayesian methods to quantify a precision needed for the Hamiltonian given an estimate of the current uncertainty in an eigenvalue.

Other embodiments comprise one or more computer-readable media storing computer-executable instructions, which when executed by a computer cause the computer to perform a method comprising inputting a Hamiltonian to be computed by the quantum computer device; reducing a number of Hamiltonian terms in the Hamiltonian using randomization within a phase estimation algorithm; and outputting a quantum circuit description for the Hamiltonian with the reduced number of Hamiltonian terms.

The foregoing and other objects, features, and advantages of the disclosed technology will become more apparent from the following detailed description, which proceeds with reference to the accompanying figures.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a quantum circuit for performing iterative phase estimation.

FIGS. 2-9 comprise graphs that show the average ground energy shift (compared to unsampled Hamiltonian), variance in ground energies over sampled Hamiltonians, average qubit requirement, and average number of terms in sampled Hamiltonians for Li₂, as a function of number of samples taken to generate the Hamiltonian and the value of the parameter p.

FIG. 10 is a flow chart showing an example method for implementing an importance sampling simulation method according to an embodiment of the disclosed technology.

FIG. 11 is a flow chart showing an example method for performing a quantum simulation using adaptive Hamiltonian randomization.

FIG. 12 illustrates a generalized example of a suitable classical computing environment in which aspects of the described embodiments can be implemented.

FIG. 13 shows an example of a possible network topology (e.g., a client-server network) for implementing a system according to the disclosed technology.

FIG. 14 shows another example of a possible network topology (e.g., a distributed computing environment) for implementing a system according to the disclosed technology.

FIG. 15 shows an exemplary system for implementing the disclosed technology.

FIG. 16 is a flow chart showing an example method for performing a quantum simulation using adaptive Hamiltonian randomization.

DETAILED DESCRIPTION I. Introduction

Not all Hamiltonian terms are created equally in quantum simulation. Hamiltonians that naturally arise from chemistry, materials and other applications often are composed of terms that are negligibly small. These terms are often culled from the Hamiltonian well before it hits the simulator. Other terms that are formally present in the Hamiltonian are removed, not because of their norm, but rather because they are not expected to impact the quantity of interest. For example, in quantum chemistry one usually selects an active space of orbitals and ignores any orbitals outside the active space. This causes many large terms to be omitted from the Hamiltonian.

This process, often called decimation, often involves systematically removing terms from the Hamiltonian and simulating the dynamics. The idea behind such a scheme is to remove terms in the Hamiltonian until the maximum shift allowed in the eigenvalues is comparable to the level of precision needed. For the case of chemistry, chemical accuracy sets a natural accuracy threshold for such simulations, but in general this precision requirement need not be viewed as a constant.

One of the example innovations of this disclosure is that, in iterative phase estimation, the number of terms taken in the Hamiltonian should ideally not be held constant. The reason why is that the high-order bits are mostly irrelevant when one is trying to learn, for example, a given bit of a binary expansion of the eigenphase. A much lower accuracy simulation can be tolerated than it can when learning a high-order bit. It then makes sense to adapt the number of terms in the Hamiltonian as iterative phase estimation proceeds through the bits of the phase estimation. Example embodiments of the disclosed technology provide a systematic method for removing terms and provides formal proofs that such processes need not dramatically affect the results of phase estimation nor its success probability.

One of the concepts behind the example advanced decimation procedure is that a form of importance sampling is used to estimate, a priori, which terms in the Hamiltonian are significant. These randomized Hamiltonians are then used within a simulation circuit to prepare approximate ground states. It is then shown that, using analysis that is reminiscent of that behind the Zeno effect or the quantum adiabatic theorem, that the errors in the eigenstate prepared at each round of phase estimation need not have a substantial impact on the posterior mean of the eigenphase estimated for the true Hamiltonian. This shows, under appropriate assumptions on the eigenvalue gaps, that this process can be used to reduce the time complexity of simulation and even under some circumstances reduce the space complexity by identifying qubits that are not needed for the level of precision asked of the simulation.

The disclosure proceeds by first reviewing iterative phase estimation and Bayesian inference, which are used to quantify the maximal error in the inference of the phase. The disclosure then proceeds to examine the effect of using a stochastic Hamiltonian on the eigenphases yielded by phase estimation in the simple case where a fixed, but random, Hamiltonian is used at each step of iterative phase estimation. The more complicated case is then examined where each repetition of e^−iHt in the iterative phase estimation circuit is implemented with a different random Hamiltonian. The theoretical analysis concludes by showing that the success probability is not degraded substantially if the eigenvalue gaps of the original Hamiltonian are sufficiently large. Further, numerical examples of this sampling procedure are shown and from that it can be concluded that the example sampling process for the Hamiltonian can have a substantial impact on the number of terms in the Hamiltonian and even in some cases the number of qubits used in the simulation.

II. Iterative Phase Estimation

The idea behind iterative phase estimation is based on the aim to build a quantum circuit that acts as an interferometer wherein the unitary one wishes to probe is applied in one of the two branches of the interferometer but not applied in the other. When the quantum state is allowed to interfere with itself at the end of the protocol the interference pattern reveals the eigenphase. This process allows the eigenvalues of U to be estimated within the standard quantum limit (i.e. the number of applications of U needed to estimate the phase within error ϵ is in Θ(1/ϵ²). If the quantum state is allowed to pass repeatedly through the interferometer circuit (or entangled inputs are used) then this scaling can be reduced to Θ(1/ϵ) which is known as the Heisenberg limit. Such a circuit is shown in schematic block diagram 100 of FIG. 1. In particular, FIG. 1 shows a quantum circuit for performing iterative phase estimation. M is the number of repetitions of the unitary U (not necessarily an integer), and θ is a phase offset between the ancilla |0

and |1

states.

The phase estimation circuit is easy to analyze in the case where U|ψ

=e^iϕ|ψ

. If U is repeated M times and θ is a phase offset then the likelihood of a given measurement outcome oϵ{0, 1} for the circuit in FIG. 1 with these parameters is

$\begin{array}{cc}\mathrm{Pr}\left(o❘\varphi ;M,\theta \right)=\frac{1+\left(-1\right)\mathrm{°cos}\left(M\left(\theta -\varphi \right)\right)}{2}.& \left(1\right)\end{array}$

There are many free parameters that can be used when designing iterative phase estimation experiments. In particular, the rules for generating M and θ for each experiment vary radically along with the methods used to process the data that comes back from these experiments. Approaches such as Kitaev's phase estimation algorithm, robust phase estimation, information theory phase estimation, or any number of approximate Bayesian methods, provide good heuristics for picking these parameters. In this disclosure, it is assumed that one does not wish to specify to any of these methods for choosing experiments, nor does one wish to focus on the specific data processing methods used. Nonetheless, Bayesian methods are relied on to discuss the impact that randomizing the Hamiltonian can have on an estimate of the eigenphase.

Bayes' theorem can be interpreted as giving the correct way to update beliefs about some fact given a set of experimental evidence and prior beliefs. The initial beliefs of the experimentalist are encoded by a prior distribution Pr(ϕ). In many cases, it is appropriate to set Pr(ϕ) to be a uniform distribution on [0, 2π) to represent a state of maximal ignorance about the eigenphase. However, in quantum simulation broader priors can be chosen if each step in phase estimation uses U_j=e^{−iHt j} and obeys U_j|ψ

=e^{−iE 0 t j} |ψ

for different t_j, since such experiments can learn E₀ as opposed to experiments with a fixed t which yield ϕ=E₀t mod 2π.

Bayes' theorem then gives the posterior distribution Pr(ϕ|o; ϕ, M) to be

$\begin{array}{cc}\mathrm{Pr}\left(\varphi ❘o;\varphi ,M\right)=\frac{\mathrm{Pr}\left(o❘\varphi ;M,\theta \right)\mathrm{Pr}\left(\varphi \right)}{\int \mathrm{Pr}\left(o❘\varphi ;M,\theta \right)\mathrm{Pr}\left(\varphi \right)d\varphi }.& \left(2\right)\end{array}$
Given a complete data set rather than a single datum, one has that

$\begin{array}{cc}\mathrm{Pr}\left(\overrightarrow{\varphi }❘\overrightarrow{o};\varphi ,\overrightarrow{M}\right)=\frac{\stackrel{}{{\prod }_j}\mathrm{Pr}\left(o_j❘\varphi ;M_j,{\theta }_j\right)\mathrm{Pr}\left(\varphi \right)}{\int {\prod }_j\mathrm{Pr}\left(o_j❘\varphi ;M_j,{\theta }_j\right)\mathrm{Pr}\left(\varphi \right)d\varphi }.& \left(3\right)\end{array}$
This probability distribution encodes the experimentalist's entire state of knowledge about ϕ given that the data is processed optimally.

It is not customary to return the posterior distribution (or an approximation thereof) as output from a phase estimation protocol. Instead, a point estimate for ϕ is given. The most frequently used estimate is the maximum a posteriori (MAP) estimate which is simply the d that has the maximum probability. While this quantity has a nice operational interpretation, it suffers from a number of deficiencies for purposes of this disclosure. The main drawback here is that the MAP estimate is not robust, in the sense that if two different values of ϕ have comparable likelihoods then small errors in the likelihood can lead to radical shifts in the MAP estimate. The posterior mean is a better estimate for this purpose, which formally is ∫ Pr(ϕ|{right arrow over (o)}; {right arrow over (M)}, {right arrow over (θ)})ϕdϕ. The posterior mean can be seen as the estimate that reduces the mean square error in any unbiased estimate of ϕ and thus it is well motivated. It also has the property that it is robust to small perturbations in the likelihood, which is a feature that is used below to estimate the impact on the results of a phase estimation experiment.

III. Errors in Likelihood Function A. Linear Combinations of Unitaries

Linear combination of unitary methods for quantum simulation have rapidly become a favored method for simulating Hamiltonian dynamics in quantum systems. Unlike Trotter decompositions, many of these methods do not necessarily yield a unitary approximation to the simulated quantum dynamics. This means that it is impossible to use Stone's theorem directly to argue that the linear combination of unitaries method implements e^{−i{tilde over (H)}t} in place of e^−iHt. In turn, since the standard analysis of the error propagation from Trotter decompositions to the estimated phase in iterative phase estimation fails because one cannot reason about the eigenvalues of {tilde over (H)} directly.

Here, to address this in part, a discussion is provided of the impact that such errors can have on the likelihood function for iterative phase estimation.

Lemma 1. Let V be a possibly non-unitary operation that can be non-deterministically performed by a quantum computer in a controlled fashion such that there exists a unitary U with ∥U−V∥≤δ<1. If one defines the likelihood function post-selected on V succeeding to be {tilde over (P)}(o|Et; M, θ) then

$P\left(o❘\mathrm{Et};M,\theta \right)-\tilde{P}\left(o❘\mathrm{Et};M,\theta \right)\le \frac{\delta }{1-\delta }.$
Proof. Let one assume that o=0. Then, one has that for input state |ψ

the error in the likelihood function output by iterative phase estimation is

$\begin{array}{cc}P\left(0❘\mathrm{Et};M,\theta \right)-\tilde{P}\left(0❘\mathrm{Et};M,\theta \right)=\mathrm{Tr}\left(\left[\frac{{\mathrm{Ue}}^{i\theta }\psi ⟩⟨\psi +e^{-i\theta }\psi ⟩⟨\psi U^{†}}{4}\right]-\left[\frac{{\mathrm{Ve}}^{i\theta }\psi ⟩⟨\psi +e^{-i\theta }\psi ⟩⟨\psi V^{†}}{4V\psi ⟩}\right]\right)\le \frac{1}{2}\mathrm{Tr}\left(U\psi ⟩⟨\psi -\frac{V\psi ⟩⟨\psi }{V\psi ⟩}\right)\le \frac{1}{2}\mathrm{Tr}\left(U\psi ⟩⟨\psi -V\psi ⟩⟨\psi \right)+\frac{1}{2}\mathrm{Tr}\left(V\psi ⟩⟨\psi -\frac{V\psi ⟩⟨\psi }{V\psi ⟩}\right)\le \frac{\delta }{2}+\frac{V}{2}\left(\frac{1}{1-\delta }-1\right)\le \frac{\delta }{2}+\frac{1+\delta }{2}\left(\frac{\delta }{1-\delta }\right)=\frac{\delta }{1-\delta }.& \left(4\right)\end{array}$

Since P(0|Et; M, θ)+P(1|Et; M, θ)=1 it follows that the same bound must apply for o=1 as well. Thus the result holds for any o as claimed.

This result, while straight forward, is significant because it allows the maximum errors in the mean of the posterior distribution to be propagated through the iterative phase estimation protocol. The ability to propagate these errors will ultimately allow us to show that iterative phase estimation can be used to estimate eigenvalues from linear combinations of unitary methods.

B. Subsampling Hamiltonians

The case where terms are sampled uniformly from the Hamiltonian is now considered. Let the Hamiltonian be a sum of L simulable Hamiltonians H_l, H

, Throughout, an eigenstate |ψ

of H and its corresponding eigenenergy E is considered. From the original, one can construct a new Hamiltonian

$\begin{array}{cc}{H}_{\mathrm{est}}=\frac{L}{m}\sum _{i=1}^m{H}_{{\ell }_i}& \left(5\right)\end{array}$
by uniformly sampling terms H_{l i} from the original Hamiltonian.

When one randomly sub-samples the Hamiltonian, errors are naturally introduced. The main question is less about how large such errors are, but rather how they impact the iterative phase estimation protocol. The following lemma states that the impact on the likelihood functions can be made arbitrarily small.

Lemma 2. Let

_i be an indexed family of sequences mapping {1, . . . . , m}→{1, . . . , L} formed by uniformly sampling elements from {1, . . . , L} independently with replacement and let {H_l:

=1, . . . , L} be a corresponding family of Hamiltonians with H=

. For |ψ

an eigenstate of H such that H|ψ

=E|ψ

and

$H_{\mathrm{samp}}=\frac{L}{m}\sum _{k=1}^m{H}_{{\ell }_i\left(k\right)}$
with corresponding eigenstate H_samp|ψ_i

=E_i|ψ_i

one then has that the error in the likelihood function for phase estimation vanishes with high probability over H_samp in the limit of large m:

$P\left(o❘\mathrm{Et};M,\theta \right)-P\left(o❘E_{i}t;M,\theta \right)\in O\left(\frac{\mathrm{MtL}}{\sqrt{m}}\sqrt{{𝕍}_{\ell }\left(⟨\psi H_{\ell }\psi ⟩\right)}\right)$
Proof. Because the terms H_{l i (k)} are uniformly sampled, each set of terms {

_i} is equally likely, and by linearity of expectation

[H_samp]=H, from which one knows that

_{i} [

ψ|H−H_est|ψ

]=0.

The second moment is easy to compute from the independence property of the distribution.

$\begin{array}{cc}{𝕍}_{\left\{i\right\}}\left(⟨\psi H_{\mathrm{samp}}\psi ⟩\right)=\frac{L^2}{m^2}{𝕍}_{\left\{i\right\}}\left(⟨\psi \sum _{k=1}^m{H}_{{\ell }_i\left(k\right)}\psi ⟩\right)=\frac{L^2}{m^2}\sum _{k=1}^m{𝕍}_{\left\{i\right\}}\left(⟨\psi H_{{\ell }_i\left(k\right)}\psi ⟩\right).& \left(6\right)\end{array}$
Since the different

_i are chosen uniformly at random the result then follows from the observation that

_{i}(

ψ|

|ψ

)=

(

ψ|

|ψ

).

From first order perturbation theory, one has that the leading order shift in any eigenvalue is O(

ψ|(H−H_samp)|ψ

) to within error O(L/√{square root over (m)}). Thus the variance in this shift is from Eq. (6)

$\begin{array}{cc}{𝕍}_{\left\{i\right\}}\left(⟨\psi \left(H-H_{\mathrm{samp}}\right)\psi ⟩\right)=\frac{L^2}{m}{𝕍}_{\ell }\left(⟨\psi H_{\ell }\psi ⟩\right).& \left(7\right)\end{array}$
This further implies that the perturbed eigenstate |ψ_i

has eigenvalue

$\begin{array}{cc}{H}_{\mathrm{samp}}{\psi }_i⟩=E{\psi }_i⟩+O\left(\frac{L}{\sqrt{m}}\sqrt{{𝕍}_{\ell }\left(⟨\psi H_{\ell }\psi ⟩\right)}\right)& \left(8\right)\end{array}$
with high probability over i from Markov's inequality. It then follows from Taylor's theorem and Eq. (1) that

$\begin{array}{cc}P\left(o❘\mathrm{Et};M,\theta \right)-P\left(o❘E_{i}t;M,\theta \right)\in O\left(\frac{\mathrm{MtL}}{\sqrt{m}}\sqrt{{𝕍}_{\ell }\left(⟨\psi H_{\ell }\psi ⟩\right)}\right),& \left(9\right)\end{array}$
with high probability over i.

This result shows that if one samples the coefficients of the Hamiltonian that are to be included in the sub-sampled Hamiltonian uniformly then one can make the error in the estimate of the Hamiltonian arbitrarily small. In this context, taking m→∞ does not cause the cost of simulation to diverge (as it would for many sampling problems). This is because once every possible term is included in the Hamiltonian, there is no point in sub-sampling and one may as well take H_samp to be II to eliminate the variance in the likelihood function that would arise from sub-sampling the Hamiltonian. In general, one needs to take mϵΩ(

(

ψ|

|ψ

)/(MtL)²) in order to guarantee that the error in the likelihood function is at most a constant. Thus, this shows that as any iterative phase estimation algorithm proceeds, that (barring the problem of accidentally exciting a state due to perturbation) one will be able to find a good estimate of the eigenphase by taking m to scale inverse quadratically with M.

IV. Bayesian Phase Estimation Using Random Hamiltonians

Theorem 3. Let E be an event and let P(E|θ) and P′(E|θ) for θϵ[−π, π) be two likelihood functions such that max_θ(|P(E|θ)−P′(E|θ)|)≤Δ and further assume that for prior P(θ) one has that min(P(E), P′(E))≥2Δ. One then has that

$\int \theta \left(P\left(\theta ❘E\right)-P^{\prime }\left(\theta ❘E\right)\right)d\theta \le \frac{5\pi \Delta }{P\left(E\right)}.\mathrm{and}\mathrm{if}$ $P\left(E❘\theta \right)=\prod _{j=1}^{N}P\left(E_j❘\theta \right)\mathrm{with}1-P^{\prime }\left(E_j❘\theta \right)/P\left(E_j❘\theta \right)\le \gamma \mathrm{then}$ $\int \theta \left(P\left(\theta ❘E\right)-P^{\prime }\left(\theta ❘E\right)\right)d\theta \le 5\pi \left({\left(1+\gamma \right)}^N-1\right).$
Proof. From the triangle inequality, one has that
|P(E)−P′(E)|=|∫P(θ)(P(E|θ)−P′(E|θ))dθ|≤Δ.  (10)
Thus it follows from the assumption that P′(E)≥2Δ that

$\left(11\right)$ $\begin{array}{c}P\left(\theta ❘E\right)-P^{\prime }\left(\theta ❘E\right)=P\left(\theta \right)\frac{P\left(E❘\theta \right)}{P\left(E\right)}-\frac{P^{\prime }\left(E❘\theta \right)}{P^{\prime }\left(E\right)}\\ \le P\left(\theta \right)\left(\frac{P\left(E❘\theta \right)}{P\left(E\right)}-\frac{P^{\prime }\left(E❘\theta \right)}{P\left(E\right)}+\frac{P^{\prime }\left(E❘\theta \right)}{P\left(E\right)}-\frac{P^{\prime }\left(E❘\theta \right)}{P^{\prime }\left(E\right)}\right)\\ \le P\left(\theta \right)\left(\frac{\Delta }{P\left(E\right)}+\frac{P^{\prime }\left(E❘\theta \right)}{P^{\prime }\left(E\right)-\Delta }-\frac{P^{\prime }\left(E❘\theta \right)}{P^{\prime }\left(E\right)}\right)\\ \le \Delta \left(\frac{P\left(\theta \right)}{P\left(E\right)}+\frac{2P^{\prime }\left(E❘\theta \right)}{P^{\prime 2}\left(E\right)}\right)\end{array}$
Thus one has that

$\begin{array}{cc}\begin{array}{c}\int \theta \left(P\left(\theta ❘E\right)-P^{\prime }\left(\theta ❘E\right)\right)d\theta \le \Delta \left(\frac{{⟨\theta ⟩}_{\mathrm{prior}}}{P\left(E\right)}+\frac{2{⟨\theta ⟩}_{\mathrm{posterior}}^{\prime }}{P^{\prime }\left(E\right)}\right),\\ \le \pi \Delta \left(\frac{1}{P\left(E\right)}+\frac{2}{P^{\prime }\left(E\right)}\right)\\ \le \pi \Delta \left(\frac{1}{P\left(E\right)}+\frac{2}{P\left(E\right)-\Delta }\right)\\ \le \frac{5\pi \Delta }{P\left(E\right)}.\end{array}& \left(12\right)\end{array}$
Now if one assumes that one has a likelihood function that factorizes into N experiments, one has that one can take

$\begin{array}{cc}\frac{\Delta }{P\left(E\right)}=\frac{\prod _{j=1}^{N}P\left(E_j❘\theta \right)-\prod _{j=1}^N{P}^{\prime }\left(E_j❘\theta \right)}{\prod _{j=1}^{N}P\left(E_j❘\theta \right)}=\prod _{j=1}^{N}1-\prod _{j=1}^N\frac{P^{\prime }\left(E_j❘\theta \right)}{P\left(E_j❘\theta \right)}.& \left(13\right)\end{array}$
From the triangle inequality

$\begin{array}{cc}\prod _{j=1}^{N}1-\prod _{j=1}^N\frac{P^{\prime }\left(E_j❘\theta \right)}{P\left(E_j❘\theta \right)}\le \prod _{j=1}^{N-1}1-\prod _{j=1}^{N-1}\frac{P^{\prime }\left(E_j❘\theta \right)}{P\left(E_j❘\theta \right)}+1-\frac{P^{\prime }\left(E_N❘\theta \right)}{P\left(E_N❘\theta \right)}{\left(1+\gamma \right)}^N\le \prod _{j=1}^{N-1}1-\prod _{j=1}^{N-1}\frac{P^{\prime }\left(E_j❘\theta \right)}{P\left(E_j❘\theta \right)}+{\gamma \left(1+\gamma \right)}^N.& \left(14\right)\end{array}$
Solving this recurrence relation gives

$\begin{array}{cc}\prod _{j=1}^{N}1-\prod _{j=1}^N{P}^{\prime }\left(E_j❘\theta \right)P\left(E_j❘\theta \right)\le {\left(1+\gamma \right)}^N-1.& \left(15\right)\end{array}$
Thus the result follows.

V. Shift in the Posterior Mean from Using Random Hamiltonians

In this section, the shift in the posterior mean of the estimated phase is analyzed assuming a random shift δ(ϕ) in the joint likelihood of all the experiments,
P′({right arrow over (o)}|ϕ;{right arrow over (M)},{right arrow over (θ)})=P({right arrow over (o)}|ϕ;{right arrow over (M)},{right arrow over (θ)})+δ(ϕ).  (16)
Here, P({right arrow over (o)}|ϕ; {right arrow over (M)}, {right arrow over (θ)}) is the joint likelihood of a series of N outcomes {right arrow over (o)} given a true phase ϕ and the experimental parameters {right arrow over (M)} and {right arrow over (θ)} for the original Hamiltonian. P′({right arrow over (o)}|ϕ; M, {right arrow over (θ)}) is the joint likelihood with a new random Hamiltonian in each experiment. By a vector like {right arrow over (M)}, the repetitions are meant for each experiment performed in the series; M_i is the number of repetitions in the i^th experiment.

First, one can work backwards from the assumption that the joint likelihood is shifted by some amount δ(ϕ), to determine an upper bound on the acceptable difference in ground state energies between the true and the random Hamiltonians. One can do this by working backwards from the shift in the joint likelihood of all experiments, to the shifts in the likelihoods of individual experiments, and finally to the corresponding tolerable differences between the ground state energies. Second, one can use this result to determine the shift in the posterior mean in terms of the differences in energies, as well as its standard deviation over the ensemble of randomly generated Hamiltonians.

A. Shifts in the Joint Likelihood

The random Hamiltonians for each experiment lead to a random shift in the joint likelihood of a series of outcomes
P′({right arrow over (o)}|ϕ;{right arrow over (M)},{right arrow over (θ)})=P({right arrow over (o)}|ϕ;{right arrow over (M)},{right arrow over (θ)})+δ(ϕ).  (17)
Assume that one would like to determine the maximum possible change in the posterior mean under this shifted likelihood. One can work under the assumption that the mean shift in the likelihood over the prior is at most |δ|≤P({right arrow over (o)})/2. The posterior is

$\begin{array}{cc}\begin{array}{c}{P}^{\prime }\left(\varphi ❘\overrightarrow{o};\overrightarrow{M},\overrightarrow{\theta }\right)=\frac{P^{\prime }\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)}{\int P^{\prime }\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)d\varphi }\\ =\frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)+\delta \left(\varphi \right)P\left(\varphi \right)}{\int \left(P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)+\delta \left(\varphi \right)P\left(\varphi \right)\right)d\varphi }\\ =\frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\theta \right)+\delta \left(\varphi \right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)+\stackrel{\_}{\delta }}.\end{array}& \left(18\right)\end{array}$
One can make progress toward bounding the shift in the posterior by first bounding the shift in the joint likelihood in terms of the shifts in the likelihoods of the individual experiments, as follows.
Lemma 4. Let P(o_j|ϕ; M_j, θ_j) be the likelihood of outcome o_j on the j^th experiment for the Hamiltonian H, and P′(o_j|ϕ; M_j, θ_j)=P(o_j|ϕ; M_j, θ_j)+ϵ_j(ϕ) be the likelihood with the randomly generated Hamiltonian H_j. Assume that N max_j(|ϵ_j(ϕ)|/P(o_j|ϕ, M_j, θ_j))<1 and |ϵ_j(ϕ)|≤P(o_j|ϕ, M_j, θ_j)/2 for all experiments j. Then the mean shift in the joint likelihood of all N experiments,
|δ|=|∫P(ϕ)(P′({right arrow over (o)}|ϕ;{right arrow over (M)},{right arrow over (θ)})−P({right arrow over (o)}|ϕ;{right arrow over (M)},{right arrow over (θ)}))dϕ|,
is at most

$\stackrel{\_}{\delta }\le 2\sum _{j=1}^N{\max }_{\varphi }\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ;M_j,{\theta }_j\right)}P\left(\overrightarrow{o}\right).$
Proof. One can write the joint likelihood in terms of the shift ϵ_j(ϕ) to the likelihoods of each of the N experiments in the sequence, P′(o_j|ϕ; M_j, θ_j)=P(o_j|ϕ; M_j, θ_j)+ϵ_j(ϕ). The joint likelihood is P′({right arrow over (o)}|ϕ; {right arrow over (M)}, {right arrow over (θ)})=Π_j=1 ^N (P(o_j|ϕ; M_j, θ_j)+ϵ_j(ϕ)), so

$\begin{array}{cc}\begin{array}{c}\log P^{\prime }\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)=\log \left(\prod _{j=1}^N\left(P\left(o_j❘\varphi ;M_j,{\theta }_j\right)+{ϵ}_j\left(\varphi \right)\right)\right)\\ =\sum _{j=1}^N\log P\left(o_j❘\varphi ,M_j,{\theta }_j\right)+\\ \log \left(1+\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ,M_j,{\theta }_j\right)}\right)\\ =\log P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)+\\ \sum _{j=1}^N\log \left(1+\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ,M_j,{\theta }_j\right)}\right)\end{array}& \left(19\right)\end{array}$
This gives one the ratio of the shifted to the unshifted joint likelihood,

$\begin{array}{cc}\frac{P^{\prime }\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)}{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)}=\exp \left[\sum _{j=1}^N\log \left(1+\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ,M_j,{\theta }_j\right)}\right)\right].& \left(20\right)\end{array}$
One can then, for example, linearize and simplify this using inequalities for the logarithm and exponential. By finding inequalities which either upper or lower bound both these functions, one can upper bound on |δ(ϕ)| in terms of the unshifted likelihoods P({right arrow over (o)}|ϕ; {right arrow over (M)}, {right arrow over (θ)}) and P(o_j|ϕ; M_j, θ_j), and the shift in the single-experiment likelihood ϵ_j(ϕ).

The inequalities that one can use to sandwich the ratio are
1−|x|≤exp(x), for x≤0   (21)
exp(x)≤1+2|x|, for x<1
−2|x|≤log(1+x), for |x|≤1/2
log(1+x)≤|x|, for x∈

In order for all four inequalities to hold, one must have that N max_j(|ϵ_j(ϕ)|/P(o_j|ϕ, M_j, θ_j))<1 (for the exponential inequalities) and |ϵ_j(ϕ)|≤P(o_j|ϕ, M_j, θ_j)/2 for all j (for the logarithm inequalities). Using them to upper bound the ratio of the shifted to the unshifted likelihood,

$\begin{array}{cc}\frac{P^{\prime }\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)}{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)}=\frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)+\delta \left(\varphi \right)}{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)}\le \exp \left[\sum _{j=1}^N\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ,M_j,{\theta }_j\right)}\right]\le 1+2\sum _{j=1}^N\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ,M_j,{\theta }_j\right)};& \left(22\right)\\ \delta \left(\varphi \right)\le 2P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)\sum _{j=1}^N\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ,M_j,{\theta }_j\right)}.& \end{array}$
On the other hand, using them to lower bound the ratio,

$\begin{array}{cc}\frac{P^{\prime }\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)}{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)}=\frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)+\delta \left(\varphi \right)}{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)}\ge \exp \left[-2\sum _{j=1}^N\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ;M_j,{\theta }_j\right)}\right]\ge 1-2\sum _{j=1}^N\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ;M_j,{\theta }_j\right)};& \left(23\right)\\ \delta \left(\varphi \right)\ge -2P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)\sum _{j=1}^N\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ;M_j,{\theta }_j\right)}.& \end{array}$
The upper and lower bounds are identical up to sign. This allows one to combine them directly, so one has

$\begin{array}{cc}\delta \left(\varphi \right)\le 2P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)\sum _{j=1}^N\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ;M_j,{\theta }_j\right)}.& \left(24\right)\end{array}$
From this, one can find an upper bound on the mean shift over the posterior, |δ|, since by the triangle inequality

$\begin{array}{cc}\stackrel{\_}{\delta }=\int \delta \left(\varphi \right)P\left(\varphi \right)d\varphi \le 2\sum _{j=1}^N\int \left(\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ;M_j,{\theta }_j\right)}P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)\right)d\varphi \le 2\sum _{j=1}^N\underset{\varphi }{\max }\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ;M_j,{\theta }_j\right)}P\left(\overrightarrow{o}\right).& \left(25\right)\end{array}$

So one has a bound on the shift in the joint likelihood in terms of the shifts in the likelihoods of individual experiments. These results allow one to bound the shift in the posterior mean in terms of the shifts in the likelihoods of the individual experiments ϵ_j(ϕ).

B. Shift in the Posterior Mean

One can use the assumption that |δ|≤P({right arrow over (o)})/2 to bound the shift in the posterior mean.

Lemma 5. Assuming in addition to the assumptions of Lemma 4 that |δ|≤P({right arrow over (o)})/2, the difference between the posterior mean that one would see with the ideal likelihood function and the perturbed likelihood function is at most

$\stackrel{\_}{\varphi }-{\stackrel{\_}{\varphi }}^{\prime }\le 8\underset{\varphi }{\max }\left(\sum _{j=1}^N\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ,M_j,{\theta }_j\right)}\right){\stackrel{\_}{\varphi }}^{\mathrm{post}}.$
Proof. One can approach the problem of bounding the difference between the posterior means by bounding the point-wise difference between the shifted posterior and the posterior with the original Hamiltonian,

$\begin{array}{cc}P\left(\varphi ❘\overrightarrow{o};\overrightarrow{M},\overrightarrow{\theta }\right)-P^{\prime }\left(\varphi ❘\overrightarrow{o};\overrightarrow{M},\overrightarrow{\theta }\right)=\frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)}-\frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)+\delta \left(\varphi \right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)+\stackrel{\_}{\delta }}.& \left(26\right)\end{array}$
As a first step, one can place an upper bound on the denominator of the shifted posterior,

$\begin{array}{cc}\left(P\left(\overrightarrow{o}\right)+\stackrel{\_}{\delta }\right)\text{:}\begin{array}{c}\frac{1}{P\left(\overrightarrow{o}\right)+\stackrel{\_}{\delta }}=\frac{1}{P\left(\overrightarrow{o}\right)}\sum _{k=0}^{\infty }{\left(\frac{-\stackrel{\_}{\delta }}{P\left(\overrightarrow{o}\right)}\right)}^k\\ =\frac{1}{P\left(\overrightarrow{o}\right)}-\frac{\stackrel{\_}{\delta }}{{P\left(\overrightarrow{o}\right)}^2}+\frac{{\stackrel{\_}{\delta }}^2}{{P\left(\overrightarrow{o}\right)}^3}\sum _{k=0}^{\infty }{\left(\frac{-\stackrel{\_}{\delta }}{P\left(\overrightarrow{o}\right)}\right)}^k\\ \le \frac{1}{P\left(\overrightarrow{o}\right)}+\frac{2\stackrel{\_}{\delta }}{{P\left(\overrightarrow{o}\right)}^2}=\frac{1+2\stackrel{\_}{\delta }/P\left(\overrightarrow{o}\right)}{P\left(\overrightarrow{o}\right)},\end{array}& \left(27\right)\end{array}$
where in the two inequalities the assumption that |δ|≤P({right arrow over (o)})/2 was used. Using this, the point-wise difference between the posteriors is at most

$\begin{array}{cc}\frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)}-\frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)+\delta \left(\varphi \right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)+\stackrel{\_}{\delta }}\le \frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)}-\frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)+\stackrel{\_}{\delta }}+\frac{\delta \left(\varphi \right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)+\stackrel{\_}{\delta }}\le \frac{2\stackrel{\_}{\delta }P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)}{{P\left(\overrightarrow{o}\right)}^2}+\frac{\delta \left(\varphi \right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)}\left(1+\frac{2\stackrel{\_}{\delta }}{P\left(\overrightarrow{o}\right)}\right)\le \frac{2\stackrel{\_}{\delta }P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)}{{P\left(\overrightarrow{o}\right)}^2}+\frac{2\delta \left(\varphi \right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)},& \left(28\right)\end{array}$
again using |δ|≤P({right arrow over (o)})/2. With this, one can bound the change in the posterior mean,

$\begin{array}{cc}\begin{array}{c}\stackrel{\_}{\varphi }-{\stackrel{\_}{\varphi }}^{\prime }\le \int \varphi ❘P\left(\varphi ❘\overrightarrow{o};\overrightarrow{M},\overrightarrow{\theta }\right)-P^{\prime }(\varphi \overrightarrow{o};\overrightarrow{M},\overrightarrow{\theta })d\varphi \\ \le \frac{2}{P\left(\overrightarrow{o}\right)}\int \varphi \left(\frac{\stackrel{\_}{\delta }P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)}+\delta \left(\varphi \right)P\left(\varphi \right)\right)d\varphi \\ \le \frac{2}{P\left(\overrightarrow{o}\right)}\int \varphi \delta \left(\varphi \right)P\left(\varphi \right)d\varphi +\frac{2\stackrel{\_}{\delta }}{P\left(\overrightarrow{o}\right)}\int \varphi \left(\frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)}\right)d\varphi \\ \le \frac{2}{P\left(\overrightarrow{o}\right)}\left(\int \varphi \delta \left(\varphi \right)P\left(\varphi \right)d\varphi +{\stackrel{\_}{\varphi }}^{\mathrm{post}}\stackrel{\_}{\delta }\right)\end{array}& \left(29\right)\end{array}$

Now, the bounds from Lemma 4 allow one to bound the shift on the posterior mean in terms of the shifts in the likelihoods of individual experiments,

(ϕ),

$\begin{array}{cc}\stackrel{\_}{\varphi }-{\stackrel{\_}{\varphi }}^{\prime }\le \frac{2}{P\left(\overrightarrow{o}\right)}\left(\int \varphi \delta \left(\varphi \right)P\left(\varphi \right)d\varphi +{\stackrel{\_}{\varphi }}^{\mathrm{post}}\stackrel{\_}{\delta }\right)\le \frac{2}{P\left(\overrightarrow{o}\right)}\left(2\underset{\varphi }{\max }\sum _{j=1}^N\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ,M_j,{\theta }_j\right)}P\left(\overrightarrow{o}\right)\int \varphi \left(\frac{P\left(\overrightarrow{o}❘\varphi ;\overrightarrow{M},\overrightarrow{\theta }\right)P\left(\varphi \right)}{P\left(\overrightarrow{o}\right)}\right)d\varphi +{\stackrel{\_}{\varphi }}^{\mathrm{post}}\stackrel{\_}{\delta }\right),& \left(30\right)\end{array}$
where in the last step, one can multiple and divide by P({right arrow over (o)}). This is

$\begin{array}{cc}\begin{array}{c}\stackrel{\_}{\varphi }-{\stackrel{\_}{\varphi }}^{\prime }\le \frac{2}{P\left(\overrightarrow{o}\right)}\left(2\underset{\varphi }{\max }\left(\sum _{j=1}^N\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ,M_j,{\theta }_j\right)}\right)P\left(\overrightarrow{o}\right){\stackrel{\_}{\varphi }}^{\mathrm{post}}+{\stackrel{\_}{\varphi }}^{\mathrm{post}}\stackrel{\_}{\delta }\right)\\ \le 8\underset{\varphi }{\max }\left(\sum _{j=1}^N\frac{{ϵ}_j\left(\varphi \right)}{P\left(o_j❘\varphi ,M_j,{\theta }_j\right)}\right){\stackrel{\_}{\varphi }}^{\mathrm{post}}.\end{array}& \left(31\right)\end{array}$

C. Acceptable Shifts in the Phase

A further question is what the bound on the shift in the posterior mean is in terms of shifts in the phase.

Theorem 6. If the assumptions of Lemma 5, for all j and x

$\in \left[-\pi ,\pi \right)P\left(o_j|\theta ;x,{\theta }_j\right)=\frac{1+{\left(-1\right)}^{o_j}\cos \left(M_j\left({\theta }_j-x\right)\right)}{2},$
for each of the N experiments, one has that the eigenphases used in PE {ϕ′_j=1, . . . N} and the eigenphase the true Hamiltonian ϕ obey |ϕ−ϕ′_j≤Δϕ and additionally P(o_j|ϕ, M_j, θ_j)ϵΘ(1) then one has that the shift in the posterior mean of the eigenphase that arises from inaccuracies in the eigenvalues in the intervening Hamiltonians obeys

$\stackrel{\_}{\varphi }-{\stackrel{\_}{\varphi }}^{\prime }\le 8\underset{\varphi }{\max }\left(\sum _{j=1}^N\frac{M_j}{P\left(o_j❘\varphi ;M_j,{\theta }_j\right)}\right)\Delta \varphi .$
Furthermore if Σ_j M_jϵO(1/ϵ_ϕ) and P(o_j|ϕ; M_j, θ_j)ϵΘ(1) for all j then

$\stackrel{\_}{\varphi }-{\stackrel{\_}{\varphi }}^{\prime }\in O\left(\frac{\Delta \varphi }{{ϵ}_{\varphi }}\right)$
Proof. One can express the shift in the posterior mean in terms of the shift in the phase applied to the ground state, Δϕ, by bounding ϵ_j(ϕ) in terms of it. Recall that the likelihood with the random Hamiltonian is
P′(o _j |ϕ;M _j,θ_j)=P(o _j |ϕ;M _j,θ_j)+ϵ_j(ϕ),  (32)
where the unshifted likelihood for the j^th experiment is P(o_j|ϕ; M_j, θ_j)=½(1+(−1)^{o j} cos(M_j(ϕ−θ_j)). Thus,
|ϵ_j(ϕ)|=½|cos(M _j(ϕ+Δϕ−θ_j)−cos(M _j(ϕ−θ_j)|≤M _j|Δϕ|,  (33)
using the upper bound on the derivative sin(x)≤|x|. In sum, one has that the error in the posterior mean in the posterior mean is at most

$\begin{array}{cc}\stackrel{\_}{\varphi }-{\stackrel{\_}{\varphi }}^{\prime }\le 8\underset{\varphi }{\max }\left(\sum _{j=1}^N\frac{M_j}{P\left(o_j|\varphi ,M_j,{\theta }_j\right)}\right){\stackrel{\_}{\varphi }}^{\mathrm{post}}\mathrm{\Delta \varphi }.& \left(34\right)\end{array}$
The result then follows from the fact that the absolute value of the posterior mean is at most π if the branch [−π, π) is chosen.

VI. Shift in the Eigenphase with a New Random Hamiltonian in Each Repetition

One can reduce the variance in the applied phase by generating a different Hamiltonian in each repetition. However, this is not without its costs: it can be viewed either as leading to a failure probability in the evolution or more generally to an additional phase shift.

The reason this reduces the variance is somewhat complex to formalize mathematically—it comes in when you compute the variance of |Δϕ|. Instead of just having the indices for a single Hamiltonian, the variance is over the indices of M_j Hamiltonians. Because of this, it only scales as M

οϕ_est] instead of M²

[ϕ_est] as it usually would (there is an underlying variance in ϕ_est). The cost is that, by reducing the variance in the phase in this way, one causes an additional shift in the phase. If one did not do this across multiple steps, one would be adiabatic all the way with the same wrong Hamiltonian. Instead, this means that one is only approximately adiabatic at the cost of the variance being lower by a factor M_j. Since the additional shift is also linear in [M_j], this can lead to an improvement. It generally requires that the gap be small, and that λ_j∝∥H_j−H_j-1∥ be small.

One then has a competition between the standard deviation being M_j

[ϕ_est] or √{square root over (M_j)}

[ϕ_est], and this new shift which is linear in M_j. So depending on the gap, it might be hard to get a rigorous bound showing that this is better, and that one should not just stick with the higher variance from a single Hamiltonian.

A. Failure Probability of the Algorithm

For phase estimation, one can reduce the variance of the estimate in the phase by randomizing within the repetitions for each experiment. For the j^th experiment with M_j repetitions (recall that M_j is not necessarily an integer), one divides into [M_j] repetitions.

Within each repetition, one can randomly generate a new Hamiltonian H_k. Each Hamiltonian H_k has a slightly different ground state and energy than all the others.

The reason this reduces the variance in the estimated phase is that the phases between repetitions are uncorrelated whereas for the single-Hamiltonian case, the variance in the phase exp(−iMϕ_est) is

[Mϕ_est]=M²

[ϕ_est], when one simulates a different random Hamiltonian in each repetition (and estimate the sum of the phases, as exp(−iΣ_k=1 ^M ϕ_k,est)), the variance is

[Σ_k=1 ^M ϕ_k,est]=Σ_k=1 ^M

[ϕ_k,est].

By evolving under a different random instantiation of the Hamiltonian in each repetition, the variance in the phase is quadratically reduced; the only cost is that the algorithm now has either a failure probability (of leaving the ground state from repetition to repetition, e.g. in the transition from the ground state of H_k−1 to the ground state of H_k) or an additional phase shift compared to the true sum of the ground state energies. The first case is simpler to analyze: it is shown in Lemma 7, provided that the gap is sufficiently small, that the failure probability can be made arbitrarily small. One can do this by viewing the success probability of the algorithm as the probability of remaining in the ground state throughout the sequence of [M_j] random Hamiltonians. In the second case, one can prove in Lemma 8 a bound on the difference between eigenvalues if the state only leaves the ground space for short intervals during the evolution.

Lemma 7. Consider a sequence of Hamiltonians {H_k}_k=1 ^M, M>1. Let γ be the minimum gap between the ground and first excited energies of any of the Hamiltonians, γ=min_k(E₁ ^k−E₀ ^k). Similarly, let λ=max_k∥H_k−H_k−1∥ be the maximum difference between any two in the sequence. The probability of leaving the ground state when transferring from H_i to H₂ through to H_M in order is at most 0<ϵ<1 provided that

$\frac{\lambda }{\gamma }<\sqrt{1-\exp \left(\frac{\log \left(1-c\right)}{M-1}\right)}.$
Proof. Let |ψ_i ^k

be the i^th eigenstate of the Hamiltonian H_k and let E_i ^k be the corresponding energy. Given that the algorithm begins in the ground state of H₁ (|ψ₀ ¹

), the probability of remaining in the ground state through all M steps is
|

ψ₀ ^M|ψ₀ ^M-1

. . .

ψ₀ ²|ψ₀ ¹

|².  (35)
This is the probability of the algorithm staying in the ground state in every segment. One can simplify this expression by finding a bound for ∥ψ₀ ^k

−|ψ₀ ^k−1

|². Let λ_kV_k=H_k−H_k−1, where one can choose λ_k such that ∥V_k∥=1 to simplify the proof. Treating λ_kV_k as a perturbation on H_k−1, the components of the shift in the ground state of H_k−1 are bounded by the derivative

$\begin{array}{cc}\frac{\partial }{\partial \lambda }⟨{\psi }_0^{k-1}|{\psi }_{\ell }^k⟩=\frac{⟨{\psi }_{\ell }^{k-1}V_k{\psi }_0^{k-1}⟩}{E_{\ell }^{k-1}-E_0^{k-1}}& \left(36\right)\end{array}$
multiplied by λ=max |λ_k|, where the maximization is over both k as well as perturbations for a given k. Using this,

$\begin{array}{cc}{⟨{\psi }_{\ell }^k|{\psi }_0^{k-1}⟩}^2\le {\lambda }^2\frac{{⟨{\psi }_{\ell }^{k-1}V_k{\psi }_0^{k-1}⟩}^2}{{\left(E_{\ell }^{k-1}-E_0^{k-1}\right)}^2}\le {\lambda }^2\frac{{⟨{\psi }_{\ell }^{k-1}V_k{\psi }_0^{k-1}⟩}^2}{{\gamma }^2}.& \left(37\right)\end{array}$
This allows one to write |ψ₀ ^k+1

=(1+δ₀)|ψ₀ ^k

+

|

, where |

≤λmax_k

$\frac{⟨{\psi }_{\ell }^k\left|V_k\right|{\psi }_0^k⟩}{E_0^k-E_{\ell }^k}.$
Letting V_k|ψ₀ ^k

=κ_k|ϕ_k), where one can again choose κ_k such that |ϕ_k

is normalized,

$\begin{array}{cc}\begin{array}{c}{{\psi }_0^k⟩-{\psi }_0^{k-1}⟩}^2={\mathrm{<figure-callout\; id="0"\; label="\delta "\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">\delta </figure-callout>}}_0^2+\sum _{\ell >0}{\delta }_{\ell }^2\le {\mathrm{<figure-callout\; id="0"\; label="\delta "\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">\delta </figure-callout>}}_0^2+\frac{{\lambda }^2}{{\gamma }^2}\sum _{\ell }{⟨{\psi }_{\ell }^{k-1}V_k{\psi }_0^{k-1}⟩}^2\\ ={\mathrm{<figure-callout\; id="0"\; label="\delta "\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">\delta </figure-callout>}}_0^2+\frac{{\lambda }^2}{{\gamma }^2}{\kappa }_k^2\sum _{\ell }{⟨{\psi }_{\ell }^k|{\varphi }_k⟩}^2\\ ={\mathrm{<figure-callout\; id="0"\; label="\delta "\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">\delta </figure-callout>}}_0^2+\frac{{\lambda }^2}{{\gamma }^2}{\kappa }_k^2.\end{array}& \left(38\right)\end{array}$
One can solve for δ₀ ² in terms of

since (1+δ₀)²+

=1. Since √{square root over (1−x)}≥1−x for xϵ[0, 1],

$\begin{array}{cc}{\mathrm{<figure-callout\; id="0"\; label="\delta "\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">\delta </figure-callout>}}_0^2={\left(1-\sqrt{1-\sum _{\ell >0}{\delta }_{\ell }^2}\right)}^2\le {\left(\sum _{\ell >0}{\delta }_{\ell }^2\right)}^2\le \sum _{\ell >0}{\delta }_{\ell }^2& \left(39\right)\end{array}$
since

. Finally, returning to ∥ψ₀ ^k

−|ψ₀ ^k−1

|², since κ_k≤1 (this is true because ∥V_k∥=1), the difference between the ground states of the two Hamiltonians is at most

$\begin{array}{cc}{{\psi }_0^k⟩-{\psi }_0^{k-1}⟩}^2\le \frac{2{\lambda }^2}{{\gamma }^2}& \left(40\right)\end{array}$

This means that the overlap probability between the ground states of any two adjacent Hamiltonians is

${⟨{\psi }_0^{k+1}|{\psi }_0^k⟩}^2\ge 1-\frac{{\lambda }^2}{{\gamma }^2}.$
Across M segments (M−1 transitions), the success probability is at least

${\left(1-\frac{{\lambda }^2}{{\gamma }^2}\right)}^{M-1}.$
If one wishes for the failure probability to be at most some fixed 0<ϵ<1, one must have

$\begin{array}{cc}{\left(1-\frac{{\lambda }^2}{{\gamma }^2}\right)}^{M-1}>1-ϵ\frac{\lambda }{\gamma }<\sqrt{1-\exp \left(\frac{\log \left(1-ϵ\right)}{M-1}\right)}.& \left(41\right)\end{array}$

If one can only prepare the ground state |ψ₀

of the original Hamiltonian, the success probability has an additional factor |

ψ₀ ¹|ψ₀

|². In this case, one can apply Lemma 7 with ∥H−H₁∥ included in the maximization for λ. Further, since γ=min_k(E₁ ^k−E₀ ^k)≤E₁−E₀−2λ|

ψ₀|H_k−H|ψ₀

≤E₁−E₀−2λ, where E₁−E₀ is the gap between the ground and first excited states of H, one needs

$\begin{array}{cc}\frac{\lambda }{E_1-E_0-2\lambda }<\sqrt{1-\exp \left(\frac{\log \left(1-ϵ\right)}{M}\right)}.& \left(42\right)\end{array}$
Provided that this occurs, one stays in the ground state of each Hamiltonian throughout the simulation with probability 1−ϵ. In this case, the total accumulated phase is

$\begin{array}{cc}\frac{\begin{array}{c}⟨{\mathrm{<figure-callout\; id="0"\; label="\psi "\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">\psi </figure-callout>}}_0^{\lceil M_j\rceil }e^{-{\mathrm{iH}}_{\lceil M_j\rceil }\Delta t}{\psi }_0^{\lceil M_j\rceil -1}⟩\dots ⟨{\psi }_0^2e^{-{\mathrm{iH}}_2\Delta t}{\psi }_0^1⟩\\ ⟨{\psi }_0^1e^{-{\mathrm{iH}}_1\Delta t}{\psi }_0⟩\end{array}}{⟨{\psi }_0^{\lceil M_j\rceil }|{\mathrm{<figure-callout\; id="0"\; label="\psi "\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">\psi </figure-callout>}}_0^{\lceil M_j\rceil -1}⟩\dots ⟨{\psi }_0^2|{\mathrm{<figure-callout\; id="0"\; label="\psi "\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">\psi </figure-callout>}}_0^1⟩⟨{\psi }_0^1|{\mathrm{<figure-callout\; id="0"\; label="\psi "\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">\psi </figure-callout>}}_0⟩}=\exp \left(-i\sum _{k=1}^{\lceil \mathrm{Mj}\rceil }{E}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">k</figure-callout>}}^0\Delta t\right),\mathrm{where}\Delta t=M_{j}t/\lceil M_j\rceil .& \left(43\right)\end{array}$

B. Phase Shifts Due to Hamiltonian Errors

One can generalize the analysis of the difference in the phase by determining the difference between the desired (adiabatic) unitary and the true one. Evolving under M random Hamiltonians in sequence, the unitary applied for each new Hamiltonian H_k is

$\begin{array}{cc}{U}_k=\exp \left(-{\mathrm{iH}}_k\Delta t\right)=\sum _{\ell }{\psi }_{\ell }^k⟩⟨{\psi }_{\ell }^ke^{-{\mathrm{iE}}_{\ell }^k\Delta t},& \left(44\right)\end{array}$
while the adiabatic unitary one would ideally apply is

$\begin{array}{cc}{U}_{k,\mathrm{ad}}=\sum _{\ell }{\psi }_{\ell }^{k+1}⟩⟨{\psi }_{\ell }^ke^{-{\mathrm{iE}}_{\ell }^k\Delta t}.& \left(45\right)\end{array}$
The difference between the two is that true time evolution U_k under H_k applies phases to the eigenstates of H_k, while the adiabatic unitary U_k,ad applies the eigenphase, and then maps each eigenstate of H_k to the corresponding eigenstate of H_k+1. This means that if the system begins in the ground state of H₁, the phase which will be applied to it by the sequence (U_{[M j ],ad}U_{[M j ]−1,ad} . . . U_2,adU_1,ad is proportional to the sum of the ground state energies of each Hamiltonian in that sequence. By comparison, U_{[M j ]}U_{[M j ]−1} . . . U₂U₁ will include contributions from many different eigenstates of the different Hamiltonians H_k.

One can bound the difference between the unitaries U_k and U_k,ad as follows.

Lemma 8. Let P₀ ^k be the projector onto the ground state of H_k, |ψ₀ ^k

, and let the assumptions of Lemma 7 hold. The difference between the eigenvalues of

$U_k{P}_0^k=\exp \left(-{\mathrm{iH}}_k\Delta t\right)P_0^k={\sum }_{\ell }{\psi }_{\ell }^k⟩⟨{\psi }_{\ell }^ke^{-{\mathrm{iE}}_{\ell }^k\Delta t}{P}_0^k\mathrm{and}$ $U_{k,\mathrm{ad}}{P}_0^k={\psi }_0^{k+1}⟩⟨{\psi }_0^ke^{-{\mathrm{iE}}_{\ell }^k\Delta t}{P}_0^k,$
where Δt is the simulation time, is at most

$\left(U_k-U_{k,\mathrm{ad}}\right)P_0^k\le \frac{2{\lambda }^2}{{\left(\gamma -2\lambda \right)}^2}.$
Proof. First, one can expand the true unitary using the resolution of the identity Σ_p|ψ_p ^k+1

ψ_p ^k+1|, the eigenstates of the next Hamiltonian, H_k+1:

$\begin{array}{cc}{U}_k=\sum _{p,\ell }{\psi }_p^{k+1}⟩⟨{\psi }_p^{k+1}{\psi }_{\ell }^k⟩⟨{\psi }_{\ell }^ke^{-{\mathrm{iE}}_{\ell }^k\Delta t}.& \left(46\right)\end{array}$
Let

=

ψ_p ^k+1|ψ_l ^k

for p≠

and 1+Δ_pp=

_p ^k+1|ψ_p ^k

when p=

. In a sense, one is writing the new eigenstate |ψ_p ^k+1

as a slight shift from the state |

): this e reason that one chooses

ψ_p ^k+1|ψ_p ^k

=1+Δ_pp. Using this definition, one can continue to simplify U_k, as

$\begin{array}{cc}{U}_k=\sum _p\left(1+{\Delta }_{\mathrm{pp}}\right){\psi }_p^{k+1}⟩⟨{\psi }_p^ke^{-{\mathrm{iE}}_p^k\Delta t}+\sum _{p\ne \ell }{\Delta }_{p\ell }{\psi }_p^{k+1}⟩⟨{\psi }_{\ell }^ke^{-{\mathrm{iE}}_{\ell }^k\Delta t}.& \left(47\right)\end{array}$
One is now well-positioned to bound ∥(U_k−U_k,ad)P₀ ^k∥. Noting that U_k,ad exactly equals the 1 in the first sum in U_k.

$\begin{array}{cc}{\left(U_k-U_{k,\mathrm{ad}}\right)P_0^k=\sum _p{\Delta }_{p\ell }{\psi }_p^{k+1}⟩⟨{\psi }_{\ell }^ke^{-{\mathrm{iE}}_{\ell }^k\Delta t}{\psi }_0^k⟩⟨{\psi }_0^k=\underset{\psi ⟩}{\max }\sum _{p\ell }{\Delta }_{p\ell }{\psi }_p^{k+1}⟩⟨{\psi }_{\ell }^ke^{-{\mathrm{iE}}_{\ell }^k\Delta t}{\psi }_0^k⟩⟨{\psi }_0^k\psi ⟩}^2={\sum _p{\Delta }_{p0}{e}^{-{\mathrm{iE}}_0^k\Delta t}}^2\le \sum _p{{\Delta }_{p0}}^2.& \left(48\right)\end{array}$

The final step in bounding ∥U_k−U_k,ad∥ is to bound |

|=

ψ_p ^k+1|

similarly to how one bounded

. For p≠

,

is given by

$\begin{array}{cc}{{{\Delta }_{p\ell }}^2=⟨{\psi }_p^{k+1}{\psi }_{\ell }^k⟩}^2\le {\lambda }^2\frac{{⟨{\psi }_p^kV_k{\psi }_{\ell }^k⟩}^2}{{\left(E_p^k-E_{\ell }^k\right)}^2}.& \left(49\right)\end{array}$
So, as with the bounds on |δ₀|² and Σ_l>0|δ_l|² in Lemma 7, ∥(U_k−U_k,ad)P₀ ^k∥ is upper bounded by

$\begin{array}{cc}\sum _p{{\Delta }_{p0}}^2=\sum _p{⟨{\psi }_p^{k+1}{\psi }_0^k⟩}^2\le \sum _p{\lambda }^2\frac{{⟨{\psi }_p^kV_k{\psi }_0^k⟩}^2}{{\left(E_p^k-E_0^k\right)}^2}\le \frac{2{\lambda }^2}{{\left(\gamma -2\lambda \right)}^2},& \left(50\right)\end{array}$
which completes the proof.
Theorem 9. Consider a sequence of Hamiltonians {H_k}_k=1 ^M, M>1. Let γ be the minimum gap between the ground and first excited energies of any of the Hamiltonians, γ=min_k(E₁ ^k−E₀ ^k). Similarly, let λ=max_k∥H_k−H_k−1∥ be the maximum difference between any two Hamiltonians in the sequence. The maximum error in the estimated eigenphases of the unitary found by the products of these M Hamiltonians is at most

${\varphi }_{\mathrm{est}}-{\varphi }_{\mathrm{true}}\le \frac{2M{\lambda }^2}{{\left(\gamma -2\lambda \right)}^2},$
with a probability of failure of at most e provided that

$\frac{\lambda }{\gamma }<\sqrt{1-\exp \left(\frac{\log \left(1-ϵ\right)}{M-1}\right)}.$
Proof. Lemma 8 gives the difference between eigenvalues of U_kP₀ ^k and U_k,adP₀ ^k. Across the entire sequence, one has

$\begin{array}{cc}U_M{P}_0^M\dots U_k{P}_0^k\dots U_1{\mathrm{<figure-callout\; id="0"\; label="P"\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">P</figure-callout>}}_0^1-U_{M,\mathrm{ad}}{P}_0^M\dots U_{k,\mathrm{ad}}{P}_0^k\dots U_{1,\mathrm{ad}}{\mathrm{<figure-callout\; id="0"\; label="P"\; filenames="US11120359-20210914-D00009.png"\; state="\{\{state\}\}">P</figure-callout>}}_0^1\le \frac{2M{\lambda }^2}{{\left(\gamma -2\lambda \right)}^2}.& \left(51\right)\end{array}$
This is the maximum possible difference between the accumulated phases for the ideal and actual sequences, assuming the system leaves the ground state for at most one repetition at a time.

The probability of leaving the groundstate as part of a Landau-Zener process instigated by the measurement at adjacent values of the Hamiltonians is, under the assumptions of Lemma 7, that the failure probability occurring at each projection is ϵ if

$\begin{array}{cc}\frac{\lambda }{\gamma }<\sqrt{1-\exp \left(\frac{\log \left(1-ϵ\right)}{M-1}\right)},& \left(52\right)\end{array}$
thus the result follows trivially from these two results.

VII. Importance Sampling

The fundamental idea behind the example approach to decimation of a Hamiltonian is importance sampling. This approach has already seen great use in coalescing, but one can use it slightly differently here. The idea behind importance sampling is to reduce the variance of the mean of a quantity by reweighting the sum. Specifically, one can write the mean of N numbers F(j) as

$\begin{array}{cc}\frac{1}{N}\sum _{j}F\left(j\right)=\sum _{j}f\left(j\right)\frac{F\left(j\right)}{\mathrm{Nf}\left(j\right)},& \left(53\right)\end{array}$
where ƒ(j) is the importance of a given term. This shows that one can view the initial unweighted average as the average of a reweighted quantity x_j/(ƒ(j)N). While this does not have an impact on the mean of x_j it can dramatically reduce the sample variance of the mean and thus is widely used in statistics to provide more accurate estimates of means. The optimal importance function to take in these cases is ƒ(j)∝|x_j|. In such cases, it is straightforward to see that the variance of the resulting distribution is in fact zero if the sign of the x_j is constant. In such cases a straight forward calculation shows that this optimal variance is

_{ƒ opt} =(

(|F|))²−(

(F))²  (54)

The optimal variance in (54) is in fact zero if the sign of the numbers is constant. While this may seem surprising, it becomes less mysterious when one notes that in order to compute the optimal importance function you need the ensemble mean that you would like to estimate. This would defeat the purpose of importance sampling in most cases. Thus if one wants to glean an advantage from importance sampling for Hamiltonian simulation then it is important to show that one can use it even with an inexact importance function that can be, for example, computed efficiently using a classical computer.

It is now shown how this robustness holds below.

Lemma 10. Let F:

_N

with

(F)=N⁻¹Σ_j=0 ^N-1 F(j) be an unknown probability distribution that can be sampled from and let {tilde over (F)}:

_N

be a known function such that for all j, |{tilde over (F)}(j)|−|F(j)|=δ_j with |δ_j|≤|F(j)|/2. If importance sampling is used with an importance function ƒ(j)=|{tilde over (F)}(j)|/Σ_k|{tilde over (F)}(k)| then the variance obeys

${𝕍}_f\left(F\right)=\frac{1}{N^2}\sum _{j}f\left(j\right){\left(\frac{F\left(j\right)}{\mathrm{Nf}\left(j\right)}\right)}^2-{\left(𝔼\left(F\right)\right)}^2\le \frac{4}{N^2}\left(\sum _k{\delta }_k\right)\left(\sum _jF\left(j\right)\right)+{𝕍}_{f_{\mathrm{opt}}}\left(F\right)$
Proof. The proof is a straight forward exercise in the triangle inequality once one uses the fact that |δ_j|≤|F(j)|/2 and the fact that 1/(1−|x|)≤1+2|x| for all xϵ[−½, ½]:

$\begin{array}{cc}\begin{array}{c}{𝕍}_f\left(F\right)=\frac{1}{N^2}\left(\sum _kF\left(k\right)+{\delta }_k\right)\left(\sum _j\frac{F^2\left(j\right)}{F\left(j\right)+{\delta }_j}\right)-{\left(𝔼\left(F\right)\right)}^2\\ \le \frac{1}{N^2}\left(\sum _kF\left(k\right)+{\delta }_k\right)\left(\sum _j\frac{F^2\left(j\right)}{F\left(j\right)-{\delta }_j}\right)-{\left(𝔼\left(F\right)\right)}^2\\ \le \frac{1}{N^2}\left(\sum _kF\left(k\right)+{\delta }_k\right)\left(\sum _jF\left(j\right)+2{\delta }_j\right)-{\left(𝔼\left(F\right)\right)}^2\\ =\frac{1}{N^2}\left(\sum _k{\delta }_k\right)\left(\sum _jF\left(j\right)+2{\delta }_j\right)+\\ \frac{1}{N^2}\left(\sum _kF\left(k\right)\right)\left(2\sum _j{\delta }_j\right)+{\left(𝔼\left(F\right)\right)}^2-{\left(𝔼\left(F\right)\right)}^2\\ \le \frac{4}{N^2}\left(\sum _k{\delta }_k\right)\left(\sum _jF\left(j\right)\right)+{𝕍}_{f_{\mathrm{opt}}}\left(F\right).\end{array}& \left(55\right)\end{array}$

This bound is tight in the sense that as max_k|δ_k|→0 the upper bound on the variance converges to (

(|F|))²−(

(F))² which is the optimal attainable variance.

In applications such as quantum chemistry simulation, what one wants to do is minimize the variance in Lemma 2. This minimum variance could be attained by choosing ƒ(j)∝|

ψ|H_j|ψ

|. However, the task of computing such a functional is at least as hard as solving the eigenvalue estimation problem that one wants to tackle. The natural approach to take is to take inspiration from Lemma 10 and instead take ƒ(j)∝|

{tilde over (ψ)}|H_j|{tilde over (ψ)}

| where |{tilde over (ψ)}

is an efficiently computable ansatz state such as a CISD state. In practice, however, the importance of a given term may not be entirely predicted by the ansatz prediction. In which case a hedging strategy can be used wherein for some ρϵ[0, 1], ƒ(j)∝(1−ρ)

H_j

+ρ∥H_j∥. This strategy allows you to smoothly interpolate between importance dictated by the magnitude of the Hamiltonian terms as well as the expectation value in the surrogate for the groundstate.

VIII. Numerical Results

This disclosure has shown that it is possible to use iterative phase estimation using a randomized Hamiltonian. To show how effective example embodiments can be, two examples of diatomic molecules, dilithium and hydrogen chloride, are considered. In both cases, the molecules are prepared in a minimal STO6G basis and use CISD states found by variationally minimizing the groundstate energy over all states within 2 excitations away from the Hartree Fock state. One can then randomly sample Hamiltonians terms and then examine several quantities of interest including the average groundstate energy, the variance in the groundstate energies and the average number of terms in the Hamiltonian. Interestingly, one can also look at the number of qubits present in the model. This can differ because some randomly sampled Hamiltonians will actually choose terms in the Hamiltonian that do not couple with the remainder of the system. In these cases, the number of qubits required to represent the state can in fact be lower than the total number that would be ordinarily expected.

One can see in FIGS. 2-5 and FIGS. 6-9 that the estimates of the ground state energy varies radically with the degree of hedging used. It is found that if ρ=1 for both cases then in all cases, one has a very large variance in the groundstate energy, as expected since importance sampling has very little impact in that case. Conversely, it is found that is one takes ρ=0, one can maximally privilege the importance of Hamiltonian terms from the CISD state which leads to very concise models but with shifts in groundstate energies that are on the order of 1 Ha for even 10⁷ randomly selected terms (some of which may be duplicates). If one instead use a modest amount of hedging (ρ=2×10⁻⁵) then one notices that the shift in the ground state energy is minimized assuming that a shift in energy of 10%; of chemical accuracy or 0.1 mHa is acceptable for Hamiltonian truncation error. For dilithium, this represents a 30% reduction in the number of terms in the Hamiltonian; whereas for HCl this reduces the number of terms in the Hamiltonian by a factor of 3. Since the cost of a Trotter-Suzuki simulation of chemistry scales super-linearly with the number of terms in the Hamiltonian this constitutes a substantial reduction in the complexity.

One can also note that for the case of dilithium, the number of qubits needed to perform the simulation varied over the different runs. In contrast Chlorine showed no such behavior. This difference arises from the fact that dilithium requires six electrons that reside in 20 spin orbitals. Hydrogen Chloride consists of eighteen electrons and the example Fock space consists of 20 spin orbitals also. As a result, nearly every spin orbital will be relevant in that which explains why the number of spin orbitals needed to express dilithium to a fixed degree of precision changes whereas it does not for HCl. This illustrates that embodiments of the disclosed randomization procedure can be used to help select an active space for a simulation on the fly as the precision needed in the Hamiltonian increases through a phase estimation procedure.

FIGS. 2-9 comprise graphs 200, 300, 400, 500, 600, 700, 800, and 900 that show the average ground energy shift (compared to unsampled Hamiltonian), variance in ground energies over sampled Hamiltonians, average qubit requirement, and average number of terms in sampled Hamiltonians for Li₂, as a function of number of samples taken to generate the Hamiltonian and the value of the parameter ρ. A term in the Hamiltonian H_α is sampled with probability ρ_α∝(1−ρ)

H_α

+ρ∥H_α∥, where the expectation value is taken with the CISD state.

IX. Example Embodiments

In this section, example methods for performing the disclosed technology are disclosed. The particular embodiments described should not be construed as limiting, as the disclosed method acts can be performed alone, in different orders, or at least partially simultaneously with one another. Further, any of the disclosed methods or method acts can be performed with any other methods or method acts disclosed herein.

FIG. 10 is a flow chart 1000 showing an example method for implementing an importance sampling simulation method according to an embodiment of the disclosed technology.

FIG. 11 is a flow chart 1100 showing an example method for performing a quantum simulation using adaptive Hamiltonian randomization.

FIG. 16 is a flow chart 1600 showing an example method for performing a quantum simulation using adaptive Hamiltonian randomization.

At 1610, a Hamiltonian to be computed by the quantum computer device is inputted.

At 1612, a number of Hamiltonian terms in the Hamiltonian is reduced using randomization within a phase estimation algorithm.

At 1614, a quantum circuit description for the Hamiltonian is output with the reduced number of Hamiltonian terms.

In certain embodiments, the reducing comprises selecting one or more random Hamiltonian terms based on an importance function; reweighting the selected random Hamiltonian terms based on an importance of each of the selected Hamiltonian random terms; and generating the quantum circuit description using the reweighted random terms. Some embodiments further comprise implementing, in the quantum computing device, a quantum circuit as described by the quantum circuit description; and measuring a quantum state of the quantum circuit. Still further embodiments comprise re-performing the method based on results from the measuring (e.g., using an iterative process). In some embodiments, the iterative process comprises computing a desired precision value for the Hamiltonian; computing a standard deviation for the Hamiltonian based on results from the implementing and measuring; and comparing the desired precision value to the standard deviation. Some embodiments further comprise changing an order of the Hamiltonian terms based on the reducing. Certain embodiments further comprise applying importance functions to terms of the Hamiltonian in a ground state; and selecting one or more random Hamiltonian terms based at least in part on the importance functions. Some embodiments comprise using importance sampling based on a variational approximation to a groundstate. Certain embodiments further comprise using adaptive Bayesian methods to quantify a precision needed for the Hamiltonian given an estimate of the current uncertainty in an eigenvalue.

Other embodiments comprise one or more computer-readable media storing computer-executable instructions, which when executed by a computer cause the computer to perform a method comprising inputting a Hamiltonian to be computed by the quantum computer device; reducing a number of Hamiltonian terms in the Hamiltonian using randomization within a phase estimation algorithm; and outputting a quantum circuit description for the Hamiltonian with the reduced number of Hamiltonian terms.

The method can comprise selecting one or more random Hamiltonian terms based on an importance function: reweighting the selected random Hamiltonian terms based on an importance of each of the selected Hamiltonian random terms; and generating the quantum circuit description using the reweighted random terms. The method can further comprise causing a quantum circuit as described by the quantum circuit description to implemented by the quantum computing device; and measuring a quantum state of the quantum circuit. The method can further comprise computing a desired precision value for the Hamiltonian; computing a standard deviation for the Hamiltonian based on results from the implementing and measuring: comparing the desired precision value to the standard deviation; and re-performing the reducing based on a result of the comparing.

Another embodiment is a system, comprising a quantum computing system; and a classical computing system configured to communicate with and control the quantum computing system. In such embodiments, the classical computing system is further configured to: input a Hamiltonian to be computed by the quantum computer device; reduce a number of Hamiltonian terms in the Hamiltonian using randomization within an iterative phase estimation algorithm; and output a quantum circuit description for the Hamiltonian with the reduced number of Hamiltonian terms. The classical computing system can be further configured to: select one or more random Hamiltonian terms based on an importance function; reweight the selected random Hamiltonian terms based on an importance of each of the selected Hamiltonian random terms; and generate the quantum circuit description using the reweighted random terms. The classical computing system can be further configured to cause a quantum circuit as described by the quantum circuit description to be implemented by the quantum computing device; and measure a quantum state of the quantum circuit. Still further, the classical computing system can be further configured to compute a desired precision value for the Hamiltonian; compute a standard deviation for the Hamiltonian based on results from the implementing and measuring; compare the desired precision value to the standard deviation; and re-perform the reducing based on a result of the comparing. In still further embodiments, the classical computing system can be further configured such that, as part of the randomization, one or more unnecessary qubits are omitted.

X. Example Computing Environments

FIG. 12 illustrates a generalized example of a suitable classical computing environment 1200 in which aspects of the described embodiments can be implemented. The computing environment 1200 is not intended to suggest any limitation as to the scope of use or functionality of the disclosed technology, as the techniques and tools described herein can be implemented in diverse general-purpose or special-purpose environments that have computing hardware.

With reference to FIG. 12, the computing environment 1200 includes at least one processing device 1210 and memory 1220. In FIG. 12, this most basic configuration 1230 is included within a dashed line. The processing device 1210 (e.g., a CPU or microprocessor) executes computer-executable instructions. In a multi-processing system, multiple processing devices execute computer-executable instructions to increase processing power. The memory 1220 may be volatile memory (e.g., registers, cache, RAM, DRAM, SRAM), non-volatile memory (e.g., ROM, EEPROM, flash memory), or some combination of the two. The memory 1220 stores software 1280 implementing tools for performing any of the disclosed techniques for operating a quantum computer to perform Hamiltonian randomization as described herein. The memory 1220 can also store software 1280 for synthesizing, generating, or compiling quantum circuits for performing any of the disclosed techniques.

The computing environment can have additional features. For example, the computing environment 1200 includes storage 1240, one or more input devices 1250, one or more output devices 1260, and one or more communication connections 1270. An interconnection mechanism (not shown), such as a bus, controller, or network, interconnects the components of the computing environment 120). Typically, operating system software (not shown) provides an operating environment for other software executing in the computing environment 1200, and coordinates activities of the components of the computing environment 1200.

The storage 1240 can be removable or non-removable, and includes one or more magnetic disks (e.g., hard drives), solid state drives (e.g., flash drives), magnetic tapes or cassettes, CD-ROMs, DVDs, or any other tangible non-volatile storage medium which can be used to store information and which can be accessed within the computing environment 1200. The storage 1240 can also store instructions for the software 1280 implementing any of the disclosed techniques. The storage 1240 can also store instructions for the software 1280 for generating and/or synthesizing any of the described techniques, systems, or quantum circuits.

The input device(s) 1250 can be a touch input device such as a keyboard, touchscreen, mouse, pen, trackball, a voice input device, a scanning device, or another device that provides input to the computing environment 1200. The output device(s) 1260 can be a display device (e.g., a computer monitor, laptop display, smartphone display, tablet display, netbook display, or touchscreen), printer, speaker, or another device that provides output from the computing environment 1200.

The communication connection(s) 1270 enable communication over a communication medium to another computing entity. The communication medium conveys information such as computer-executable instructions or other data in a modulated data signal. A modulated data signal is a signal that has one or more of its characteristics set or changed in such a manner as to encode information in the signal. By way of example, and not limitation, communication media include wired or wireless techniques implemented with an electrical, optical, RF, infrared, acoustic, or other carrier.

As noted, the various methods and techniques for performing Hamiltonian randomization, for controlling a quantum computing device, to perform circuit design or compilation/synthesis as disclosed herein can be described in the general context of computer-readable instructions stored on one or more computer-readable media. Computer-readable media are any available media (e.g., memory or storage device) that can be accessed within or by a computing environment. Computer-readable media include tangible computer-readable memory or storage devices, such as memory 1220 and/or storage 1240, and do not include propagating carrier waves or signals per se (tangible computer-readable memory or storage devices do not include propagating carrier waves or signals per se).

Various embodiments of the methods disclosed herein can also be described in the general context of computer-executable instructions (such as those included in program modules) being executed in a computing environment by a processor. Generally, program modules include routines, programs, libraries, objects, classes, components, data structures, and so on, that perform particular tasks or implement particular abstract data types. The functionality of the program modules may be combined or split between program modules as desired in various embodiments. Computer-executable instructions for program modules may be executed within a local or distributed computing environment.

An example of a possible network topology 1300 (e.g., a client-server network) for implementing a system according to the disclosed technology is depicted in FIG. 13. Networked computing device 1320 can be, for example, a computer running a browser or other software connected to a network 1312. The computing device 1320 can have a computer architecture as shown in FIG. 12 and discussed above. The computing device 1320 is not limited to a traditional personal computer but can comprise other computing hardware configured to connect to and communicate with a network 1312 (e.g., smart phones, laptop computers, tablet computers, or other mobile computing devices, servers, network devices, dedicated devices, and the like). Further, the computing device 1320 can comprise an FPGA or other programmable logic device. In the illustrated embodiment, the computing device 1320 is configured to communicate with a computing device 1330 (e.g., a remote server, such as a server in a cloud computing environment) via a network 1312. In the illustrated embodiment, the computing device 1320 is configured to transmit input data to the computing device 1330, and the computing device 1330 is configured to implement a technique for controlling a quantum computing device to perform any of the disclosed embodiments and/or a circuit generation/compilation/synthesis technique for generating quantum circuits for performing any of the techniques disclosed herein. The computing device 1330 can output results to the computing device 1320. Any of the data received from the computing device 1330 can be stored or displayed on the computing device 1320 (e.g., displayed as data on a graphical user interface or web page at the computing devices 1320). In the illustrated embodiment, the illustrated network 1312 can be implemented as a Local Area Network (LAN) using wired networking (e.g., the Ethernet IEEE standard 802.3 or other appropriate standard) or wireless networking (e.g. one of the IEEE standards 802.11a, 802.11b, 802.11g, or 802.11n or other appropriate standard). Alternatively, at least part of the network 1312 can be the Internet or a similar public network and operate using an appropriate protocol (e.g., the HTTP protocol).

Another example of a possible network topology 1400 (e.g., a distributed computing environment) for implementing a system according to the disclosed technology is depicted in FIG. 14. Networked computing device 1420 can be, for example, a computer running a browser or other software connected to a network 1112. The computing device 1420 can have a computer architecture as shown in FIG. 12 and discussed above. In the illustrated embodiment, the computing device 1420 is configured to communicate with multiple computing devices 1430, 1431, 1432 (e.g., remote servers or other distributed computing devices, such as one or more servers in a cloud computing environment) via the network 1412. In the illustrated embodiment, each of the computing devices 1430, 1431, 1432 in the computing environment 1400 is used to perform at least a portion of the Hamiltonian randomization technique and/or at least a portion of the technique for controlling a quantum computing device to perform any of the disclosed embodiments and/or a circuit generation/compilation/synthesis technique for generating quantum circuits for performing any of the techniques disclosed herein. In other words, the computing devices 1430, 1431, 1432 form a distributed computing environment in which aspects of the techniques for performing any of the techniques as disclosed herein and/or quantum circuit generation/compilation/synthesis processes are shared across multiple computing devices. The computing device 1420 is configured to transmit input data to the computing devices 1430, 1431, 1432, which are configured to distributively implement such as process, including performance of any of the disclosed methods or creation of any of the disclosed circuits, and to provide results to the computing device 1420. Any of the data received from the computing devices 1430, 1431, 1432 can be stored or displayed on the computing device 1420 (e.g., displayed as data on a graphical user interface or web page at the computing devices 1420). The illustrated network 1412 can be any of the networks discussed above with respect to FIG. 13.

With reference to FIG. 15, an exemplary system for implementing the disclosed technology includes computing environment 1500. In computing environment 1500, a compiled quantum computer circuit description (including quantum circuits for performing any of the disclosed techniques as disclosed herein) can be used to program (or configure) one or more quantum processing units such that the quantum processing unit(s) implement the circuit described by the quantum computer circuit description.

The environment 1500 includes one or more quantum processing units 1502 and one or more readout device(s) 1508. The quantum processing unit(s) execute quantum circuits that are precompiled and described by the quantum computer circuit description. The quantum processing unit(s) can be one or more of, but are not limited to: (a) a superconducting quantum computer; (b) an ion trap quantum computer; (c) a fault-tolerant architecture for quantum computing; and/or (d) a topological quantum architecture (e.g., a topological quantum computing device using Majorana zero modes). The precompiled quantum circuits, including any of the disclosed circuits, can be sent into (or otherwise applied to) the quantum processing unit(s) via control lines 1506 at the control of quantum processor controller 1520. The quantum processor controller (QP controller) 1520 can operate in conjunction with a classical processor 1510 (e.g., having an architecture as described above with respect to FIG. 12) to implement the desired quantum computing process. In the illustrated example, the QP controller 1520 further implements the desired quantum computing process via one or more QP subcontrollers 1504 that are specially adapted to control a corresponding one of the quantum processor(s) 1502. For instance, in one example, the quantum controller 1520 facilitates implementation of the compiled quantum circuit by sending instructions to one or more memories (e.g., lower-temperature memories), which then pass the instructions to low-temperature control unit(s) (e.g., QP subcontroller(s) 1504) that transmit, for instance, pulse sequences representing the gates to the quantum processing unit(s) 1502 for implementation. In other examples, the QP controller(s) 1520 and QP subcontroller(s) 1504 operate to provide appropriate magnetic fields, encoded operations, or other such control signals to the quantum processor(s) to implement the operations of the compiled quantum computer circuit description. The quantum controller(s) can further interact with readout devices 1508 to help control and implement the desired quantum computing process (e.g., by reading or measuring out data results from the quantum processing units once available, etc.)

With reference to FIG. 15, compilation is the process of translating a high-level description of a quantum algorithm into a quantum computer circuit description comprising a sequence of quantum operations or gates, which can include the circuits as disclosed herein (e.g., the circuits configured to perform one or more of the procedures as disclosed herein). The compilation can be performed by a compiler 1522 using a classical processor 1510 (e.g., as shown in FIG. 12) of the environment 1500 which loads the high-level description from memory or storage devices 1512 and stores the resulting quantum computer circuit description in the memory or storage devices 1512.

In other embodiments, compilation and/or verification can be performed remotely by a remote computer 1560 (e.g., a computer having a computing environment as described above with respect to FIG. 12) which stores the resulting quantum computer circuit description in one or more memory or storage devices 1562 and transmits the quantum computer circuit description to the computing environment 1500 for implementation in the quantum processing unit(s) 1502. Still further, the remote computer 1500 can store the high-level description in the memory or storage devices 1562 and transmit the high-level description to the computing environment 1500 for compilation and use with the quantum processor(s). In any of these scenarios, results from the computation performed by the quantum processor(s) can be communicated to the remote computer after and/or during the computation process. Still further, the remote computer can communicate with the QP controller(s) 1520 such that the quantum computing process (including any compilation, verification, and QP control procedures) can be remotely controlled by the remote computer 1560. In general, the remote computer 1560 communicates with the QP controller(s) 1520, compiler/synthesizer 1522, and/or verification tool 1523 via communication connections 1550.

In particular embodiments, the environment 1500 can be a cloud computing environment, which provides the quantum processing resources of the environment 1500 to one or more remote computers (such as remote computer 1560) over a suitable network (which can include the internet).

XI. Concluding Remarks

This application has shown that iterative phase estimation is more flexible than it previously may have been thought and that the number of terms in the Hamiltonian be randomized at each step of iterative phase estimation without substantially contributing to the underlying variance of an unbiased estimator of the eigenphase. It was further shown numerically that by using such strategies for sub-sampling the Hamiltonian terms that one can perform a simulation using fewer Hamiltonian terms than would be necessary for traditional approaches require. These reductions in the number of terms directly impact the complexity of Trotter-Suzuki based simulation and indirectly impact qubitization and truncated Taylor series simulation methods because they also reduce the 1-norm of the vector of Hamiltonian terms.

Having described and illustrated the principles of the disclosed technology with reference to the illustrated embodiments, it will be recognized that the illustrated embodiments can be modified in arrangement and detail without departing from such principles. For instance, elements of the illustrated embodiments shown in software may be implemented in hardware and vice-versa. Also, the technologies from any example can be combined with the technologies described in any one or more of the other examples. It will be appreciated that procedures and functions such as those described with reference to the illustrated examples can be implemented in a single hardware or software module, or separate modules can be provided. The particular arrangements above are provided for convenient illustration, and other arrangements can be used.

Claims (20)

The invention claimed is:

1. A method of operating a quantum computing

device, comprising:

inputting a Hamiltonian to be computed by the quantum

computing device;

reducing a number of Hamiltonian terms in the Hamiltonian using randomization within a phase estimation algorithm; and

outputting a quantum circuit description for the Hamiltonian with the reduced number of Hamiltonian terms.

2. The method of claim 1, wherein the method is performed by a classical computer.

3. The method of claim 1, wherein the reducing comprises:

selecting one or more random Hamiltonian terms based on an importance function;

reweighting the selected random Hamiltonian terms based on an importance of each of the selected random Hamiltonian terms;

generate the quantum circuit description using the reweighted random terms.

4. The method of claim 3, further comprising:

implementing, in the quantum computing device, a quantum circuit as described by the quantum circuit description;

measuring a quantum state of the quantum circuit.

5. The method of claim 4, further comprising re-performing the method of claim 4 based on results from the measuring.

6. The method of claim 5, herein the re-performing is performed based on an iterative process.

7. The method of claim 6, wherein the iterative process comprises:

computing a desired precision value for the Hamiltonian;

computing a standard deviation for the Hamiltonian based on results from the implementing and measuring; and

comparing the desired precision value to the standard deviation.

8. The method of claim 1, further comprising changing an order of the Hamiltonian terms based on the reducing.

9. The method of claim 1, further comprising:

applying importance functions to terms of the Hamiltonian in a ground state; and

selecting one or more random Hamiltonian terms based at least in part on the importance functions.

10. The method of claim 1, further comprising:

using importance sampling based on a variational approximation to a groundstate.

11. The method of claim 1, further comprising:

using adaptive Bayesian methods to quantify a precision needed for the Hamiltonian given an estimate of the current uncertainty in an eigenvalue.

12. One or more computer-readable media storing computer-executable instructions, which when executed by a computer cause the computer to perform a method, the method comprising:

inputting a Hamiltonian to be computed by a quantum computer device;

reducing a number of Hamiltonian terms in the Hamiltonian using randomization within a phase estimation algorithm; and

outputting a quantum circuit description for the Hamiltonian with the reduced number of Hamiltonian terms.

13. The one or more computer-readable media of claim 12, wherein the method further comprises:

selecting one or more random Hamiltonian terms based on an importance function;

reweighting the selected random Hamiltonian terms based on an importance of each of the selected random Hamiltonian terms; and

generating the quantum circuit description using the reweighted random terms.

14. The one or more computer-readable media of claim 13, wherein the method further comprises:

causing a quantum circuit as described by the quantum circuit description to be implemented by the quantum computing device; and

measuring a quantum state of the quantum circuit.

15. The one or more computer-readable media of claim 14, wherein the method further comprises:

computing a desired precision value for the Hamiltonian;

computing a standard deviation for the Hamiltonian based on results from the implementing and measuring;

comparing the desired precision value to the standard deviation; and

re-performing the reducing based on a result of the comparing.

16. A system, comprising:

a quantum computing system; and

a classical computing system configured to communicate with and control the quantum computing system, the classical computing system being further configured to:

input a Hamiltonian to be computed by the quantum computing system;

reduce a number of Hamiltonian terms in the Hamiltonian using randomization within an iterative phase estimation algorithm; and

output a quantum circuit description for the Hamiltonian with the reduced number of Hamiltonian terms.

17. The system of claim 16, wherein the classical computing system is further configured to:

select one or more random Hamiltonian terms based on an importance function;

reweight the selected random Hamiltonian terms based on an importance of each of the selected random Hamiltonian terms; and

generate the quantum circuit description using the reweighted random terms.

18. The system of claim 17, wherein the classical computing system is further configured to:

cause a quantum circuit as described by the quantum circuit description to be implemented by the quantum computing device; and

measure a quantum state of the quantum circuit.

19. The system of claim 18, wherein the classical computing system is further configured to:

compute a desired precision value for the Hamiltonian;

compute a standard deviation for the Hamiltonian based on results from the implementing and measuring;

compare the desired precision value to the standard deviation; and

re-perform the reducing based on a result of the comparing.

20. The system of claim 16, wherein the classical computing system is further configured to omit one or more unnecessary qubits based on the reduced number of Hamiltonian terms.

Images