Estimation and uncertainty quantification for the output from quantum simulators

TL;DR

This paper develops methods for estimating and quantifying uncertainty in distributions generated by quantum simulators, using maximum entropy, Monte Carlo, and Bayesian techniques, with applications to quantum systems of multiple qubits.

Contribution

It introduces a novel combination of maximum entropy, Robbins Monro, and Langevin MCMC methods for quantum distribution estimation and uncertainty quantification.

Findings

Effective estimation of quantum distributions from moments

Unbiased gradient estimators enable convergent MCMC sampling

Methods demonstrated on simulated quantum system outputs

Abstract

The problem of estimating certain distributions over $\{0,1\}^d$ is considered here. The distribution represents a quantum system of $d$ qubits, where there are non-trivial dependencies between the qubits. A maximum entropy approach is adopted to reconstruct the distribution from exact moments or observed empirical moments. The Robbins Monro algorithm is used to solve the intractable maximum entropy problem, by constructing an unbiased estimator of the un-normalized target with a sequential Monte Carlo sampler at each iteration. In the case of empirical moments, this coincides with a maximum likelihood estimator. A Bayesian formulation is also considered in order to quantify posterior uncertainty. Several approaches are proposed in order to tackle this challenging problem, based on recently developed methodologies. In particular, unbiased estimators of the gradient of the log posterior are constructed and used within a provably convergent Langevin-based Markov chain Monte Carlo method. The methods are illustrated on classically simulated output from quantum simulators.

Figures

Equations

$$\rho =\sum _{j=1}^{2^d}{p}_j|{\psi }_j⟩⟨{\psi }_j|$$

$$O=\sum _{j=1}^{2^d}{o}_j|{\psi }_j⟩⟨{\psi }_j|$$

$$p(x;A)=\frac{1}{Z(A)}\exp (x^{T}Ax),Z(A)=\sum _{x\in \{0,1{\}}^d}\exp (x^{T}Ax),$$

$$m=\Pi (\varphi |\lambda ):=\sum _{x\in E^d}\varphi (x)\pi (x|\lambda ),$$

$$H(p)=-\sum _{x\in E^d}p(x)\log (p(x)),$$

$$q(x|\lambda )=\exp ({\lambda }^T\varphi (x)),$$

$$\Pi (m-\varphi |\lambda )=\frac{Q(m-\varphi |\lambda )}{Q(1|\lambda )}=0.$$

$$\mathbb{P}(Y|\lambda )=\prod _{i=1}^M\pi (y^{(i)}|\lambda ).$$

$${\nabla }_{\lambda }\log \mathbb{P}(Y|\lambda )=M\Pi (\hat{m}-\varphi |\lambda ),$$

$$\hat{m}=\frac{1}{M}\sum _{i=1}^M\varphi (y^{(i)}).$$

$$\Pi (\hat{m}-\varphi |\lambda )=0,$$

$$\mathbb{P}(\lambda |Y)\propto \mathbb{P}(Y|\lambda )\mathbb{P}(\lambda ),$$

$${\eta }_j(dx)={\Gamma }_j(dx)/{\Gamma }_j(1),$$

$${\eta }_J(\phi ):={\int }_U\phi (x){\eta }_J(dx).$$

$${\eta }_J(\phi )\approx \frac{1}{N}\sum _{i=1}^N\phi (x_J^i),$$

$${\gamma }_j(x)=\exp \left[\sum _{i=1}^j{\lambda }_i{\varphi }_i(x)\right].$$

$$G_j(x)={\gamma }_{j+1}(x)/{\gamma }_j(x).$$

$$({\eta }_j{M}_j)(dx):={\int }_U{\eta }_j(d{x}^{\prime })M_j(x^{\prime },dx)={\eta }_j(dx).$$

$${\eta }_j^N(\phi ):=\frac{1}{N}\sum _{i=1}^N\phi (x_j^i).$$

$$0<G_j(x)\le C.$$

$$\mathbb{E}|{\eta }_j^N(\phi )-{\eta }_j(\phi ){|}^p\le \frac{C(j)}{N^{p/2}}.$$

$${\eta }_j(G_j)=\frac{1}{{\Gamma }_j(1)}\sum _{x\in U}\frac{{\gamma }_{j+1}(x)}{{\gamma }_j(x)}{\gamma }_j(x)={\Gamma }_{j+1}(1)/{\Gamma }_j(1).$$

$${\Gamma }_J(1)=\prod _{j=0}^{J-1}{\eta }_j(G_j).$$

$${\eta }_j^N(G_j)=\frac{1}{N}\sum _{i=1}^N{G}_j(x_j^i).$$

$${\Gamma }_J^N(\phi ):=\prod _{j=0}^{J-1}{\eta }_j^N(G_j){\eta }_J^N(\phi ).$$

$$\frac{Q(m-\varphi |\lambda )}{Q(1|\lambda )}=0.$$

$${\lambda }^{k+1}={\lambda }^k+{\delta }_k\hat{F}({\lambda }^k).$$

$$c<q(x|\lambda )<c^{-1}.$$

$$\sum _{l=0}^{\infty }\mathbb{E}|{\Delta }_l|<\infty ,$$

$$\mathbb{E}Z=\sum _{l=0}^{\infty }\mathbb{E}{\Delta }_l=\sum _{l=0}^{\infty }{p}_l\frac{\mathbb{E}{\Delta }_l}{p_l}.$$

Full text

∗Oak Ridge National Laboratory, Oak Ridge, TN, 37831

+Department of Statistics and Applied Probability, National University of Singapore, Singapore

†School of Mathematics, University of Manchester, Manchester, UK, M13 4PL

$Oak Ridge National Laboratory, Oak Ridge, TN, 37831

1. Introduction and motivation

Quantum computing holds a great promise to enable future simulations of quantum systems that would otherwise suffer from the curse of dimensionality when simulated classically. However, we are currently entering the Noisy Intermediate Scale Quantum (NISQ) device era [40]. The hallmarks of NISQ era are relatively small (10-50 qubits) quantum computing devices and noisy qubits. Therefore, to succeed with quantum simulation on NISQ quantum computers, it is imperative to understand and fully characterize the effects of noise on the state of a quantum device.

In quantum mechanics the state of a $d$-qubit quantum device is generally described by a density matrix $\rho$ – a $2^{d}$-by-$2^{d}$ complex-valued positive semidefinite matrix with ${\rm Tr}\rho=1$. From the definition, it takes $2^{2d}-1$ real parameters to fully specify an arbitrary quantum state. If $O$ denotes a quantum observable, described by a $2^{d}$-by-$2^{d}$ Hermitian matrix ($O^{\dagger}=O$), then measuring $O$ yields one of its eigenvalues sampled from a distribution determined by $\rho$. In particular, the probability of the $i$th eigenvalue is $p_{i}=\langle\psi_{i}|\rho|\psi_{i}\rangle$, where $|\psi_{i}\rangle$ is the $i$th eigenvector in the bra-ket notation. Since there are at most $2^{d}$ eigenvalues, repeated measurements of any observable yields information about at most $2^{d}-1$ parameters of $\rho$. To estimate every parameter of the state one must measure at least $2^{d}+1$ independent observables. This approach is known as quantum state tomography [37, 24]. Its applicability is limited to small quantum devices due to exponential scaling in the number of measurements and related issues of solving exponentially large systems of equations. Bayesian methods also have been suggested [29, 5, 25, 20] to facilitate uncertainty quantification (UQ) in the context of quantum state tomography. However, the existing Bayesian quantum tomography solutions face similar scalability challenge as the number of parameters to infer still grows exponentially with the size of a quantum device.

Here we consider a fundamentally new approach, which provides a scalable solution to UQ for the output of quantum simulations. We will restrict our attention to a simplified scenario in which the density matrix and observable commute, hence are simultaneously diagonalizable, and the eigenvectors are known. In this way we are able to treat the quantum state, which is represented by the density matrix, as a classical probability distribution, given by the diagonal of the density matrix. Note that the standard quantum state tomography methods would still scale exponentially even for this simplified scenario. To avoid this problem we use the maximum entropy (MaxEnt) principle [27, 28] to model a complete probability distribution based on a selected subset of observed or exact moments. We then use well-established stochastic simulation algorithms [43, 30, 13, 14] to find the MaxEnt distribution with a cost which is polynomial in $d$. Adopting a MaxEnt ansatz for the likelihood provides a Bayesian posterior which (i) converges to the MaxEnt distribution as the number of observations tends to infinity, and (ii) can be simulated with a cost which is polynomial in $d$ using recent stochastic simulation methods [32, 42, 49, 45]. It is the topic of ongoing investigation to generalize the results to density matrices, which will be reported in a future publication. It is hoped that a polynomial cost in $d$ can be preserved in this case as well, for example combining the methodology here with that of [22, 8].

Assume that we know the basis in which the output state from a $d$-qubit quantum computer is diagonal. The density matrix in that basis reads,

[TABLE]

where $\{|\psi_{j}\rangle\}_{j=1}^{2^{d}}$ comprise a basis in the Hilbert space of the quantum computer $\mathbb{C}^{2^{d}}$. We also assume the observable $O$ that we measure is diagonal in the same basis i.e.

[TABLE]

where $\{o_{j}\}_{j=1}^{2^{d}}$ are the eigenvalues of $O$ and $p_{j}=\langle\psi_{j}|\rho|\psi_{j}\rangle$ is the probability of observing state a measurement outcome corresponding to the eigenvalue $o_{j}$.

Without loss of generality we can interpret the probability distribution $p$ as a classical probability density $p(x)$ over the binary random variable $x\in\{0,1\}^{d}$. The mathematical task at hand is then to find $p(x)$ consistent with the observed outcomes $\{o_{j}\}$. A canonical example of such density in our context is the Ising model

[TABLE]

which is parametrized by some matrix $A\in\mathbb{R}^{d\times d}$.

Given first and second moments of a distribution $\Pi$ on $\{0,1\}^{d}$, the MaxEnt distribution is the Ising model (3) which satisfies the moment constraints [11] . With finitely many observations, one can replace the moments with observed empirical moments, and the corresponding MaxEnt distribution is the same one which maximizes the likelihood under a MaxEnt ansatz for its parameterized form [34]. Through a Bayesian framework [21], this likelihood can be combined with a prior on the parameters, resulting in a posterior distribution, which quantifies the uncertainty in the distribution given the observations. We will develop and implement novel algorithms for the solution of these problems with a cost which is polynomial in $d$, based on recently developed technologies mentioned above. Now we proceed with a more rigorous mathematical setup, and a summary of the results of the paper.

2. Mathematical Setup and Summary of Results

Let $E=\{0,1\}$, and denote the state $x=[x_{1},\dots,x_{d}]\in E^{d}$. Consider a probability distribution $\Pi:\sigma(E^{d})\times\Xi\rightarrow[0,1]$, parameterized by $\lambda\in\Xi$, with density $\pi:E^{d}\times\Xi\rightarrow[0,1]$. The target $\Pi(\cdot|\lambda)$ represents the distribution over $d$ dependent qubits $X=[X_{1},\dots,X_{d}]^{T}$ given a particular value of $\lambda$, and the objective of this work is to identify $\lambda$. We will denote the un-normalized measure by $Q:\sigma(E^{d})\times\Xi\rightarrow[0,1]$, with density $q:E^{d}\times\Xi\rightarrow[0,1]$, so that $\Pi=\frac{1}{Z(\lambda)}Q$ and $\pi=\frac{1}{Z(\lambda)}q$, with $Z(\lambda)=Q(1|\lambda):=\sum_{x\in E^{d}}1q(x|\lambda)$.

Let $\phi_{j}:E^{d}\rightarrow\mathbb{R}$, for $j=1,\dots,J$, and define $\phi=[\phi_{1},\dots,\phi_{J}]^{T}:E^{d}\rightarrow\mathbb{R}^{J}$. It is assumed that, for a particular unknown value of $\lambda$, we are given either

• (P1)

$J$ exact moments $m=[m_{1},\dots,m_{J}]^{T}$:

[TABLE]

• (P2)

$M$ independent and identically distributed (i.i.d.) noisy observations of this process $y^{(i)}\sim\Pi(\cdot|\lambda)$, for $i=1,\dots,M$.

As mentioned above, the aim in either case is to reconstruct $\Pi$ by identifying $\lambda$.

First, consider problem (P(P1)). Assume we seek the probability distribution $p$ which maximizes the entropy

[TABLE]

subject to the constraints (4), with $p$ replacing $\pi(\cdot|\lambda)$. This is in fact a convex optimization problem. Therefore, if there exists a probability distribution $p$ which satisfies the constraints (4), then the necessary first order optimality conditions are also sufficient for global optimality [7]. For this particular problem, the first order conditions can be simplified to $p(x)=\pi(x|\lambda)=\frac{1}{Z(\lambda)}q(x|\lambda)$, where

[TABLE]

and the vector of Lagrange multipliers $\lambda=[\lambda_{1},\dots,\lambda_{J}]^{T}$ is the solution to the system of equations

[TABLE]

This is referred to as the method of maximum entropy (MaxEnt). See Chapter 11.1 of [11] for the derivation. A density which is parametrized in the form (6) is more generally referred to as a MaxEnt ansatz [34].

In general, the system of equations (7) does not have a closed form solution, and it is a challenging problem to solve, particularly for large $d$, as one evaluation of the system would consist of $\mathcal{O}(2^{d})$ operations. The issue is known as the curse of dimensionality, and precludes the use of standard iterative algorithms for solving the problem. We propose to utilize a Robbins Monro type [43, 30] stochastic iterative algorithm for its solution, together with a sequential Monte Carlo (SMC) sampler [15] unbiased estimator of a system of equations equivalent to (7), $Q(m-\phi|\lambda)=0$, within each iteration of the Robbins Monro. It is well-known that SMC samplers are stable with $\mathcal{O}(d^{2})$ operations for such targets which can be evaluated pointwise (in $x$) [3], thus bypassing the curse of dimensionality. That is the first major contribution of the paper.

Now, consider problem (P(P2)), concerned with recovering $\pi(\cdot|\lambda)$ from the i.i.d. observations $Y:=\{y^{(i)}\}_{i=1}^{M}$. Following the MaxEnt ethos, we choose an appropriate set of functions $\phi$ as above (4), based on our assumptions of the underlying quantum system, and adopt the ansatz (6) for the likelihood of a single observation, given $\lambda$: $\mathbb{P}(y^{(i)}|\lambda)=\pi(y^{(i)}|\lambda)$. For example, one could use the set of first and second moments $\phi(x)={\rm vec}(\{x_{i}x_{j}\}_{1\leq i\leq j\leq d})$. The likelihood of $Y:=\{y^{(i)}\}_{i=1}^{M}$ given $\lambda$ is given by

[TABLE]

The maximizer of (8) is called the maximum likelihood estimator (MLE). It is equivalent to the maximizer of $\log\mathbb{P}(Y|\lambda)$. We have

[TABLE]

where

[TABLE]

The necessary first order optimality condition is vanishing gradient, i.e.

[TABLE]

which is clearly equivalent to to the MaxEnt solution (7), except with empirical moments replacing exact moments. This problem can therefore be solved using the same method described above.

This solution is still somewhat unsettling. Even if we are content with our choice of $\phi$, we have arrived at a point estimate for inference given a finite set of observations. The scrupulous Bayesian will not allow the story to end here. Indeed for the sake of uncertainty quantification (UQ), we require a posterior distribution over $\pi(\cdot|\lambda)$, or equivalently, a posterior distribution over $\lambda$. In order to accomplish this, we complete the picture with a prior on $\lambda$, giving

[TABLE]

where the constant of proportionality depends only upon the observations $Y$. The second major contribution of the paper concerns quantifying a posteriori uncertainty using this approach. Due to the fact that the likelihood $\mathbb{P}(Y|\lambda)$ itself has an intractable normalizing constant $Q(1|\lambda)^{M}$, this notoriously difficult problem is referred to as doubly-intractable [33, 35].

If we could at least evaluate a non-negative and unbiased estimator of the likelihood $\mathbb{P}(Y|\lambda)$, or equivalently (12), then there are algorithms which can be used, such as pseudo-marginal Markov chain Monte Carlo (MCMC) method [2, 1]. The recent work [31] proposes to construct unbiased but signed estimates of the target (12), and then they use such estimate within pseudo-marginal MCMC by considering the absolute value of the estimate in the acceptance ratio and then storing the sign of the accepted samples in order to correct the ultimate estimator. The work [48] proposes to use alternative unbiased estimators in this same framework. Note that earlier work also appeared which allows one to sample from such a target [33, 35], however these works require being able to sample from $\pi(\cdot|\lambda)$ for a given value of $\lambda$. This is possible in some cases, for example the Ising model (or Boltzmann machine) [41, 9], but not in general. The works [31, 48] are generally applicable.

We will take a new approach to doubly intractable simulation in the present work. A number of recently developed methods require only unbiased estimators of (9). These include the stochastic gradient Langevin dynamics (SGLD) [49], the zig-zag sampler [4], or the bouncy particle sampler [6, 39]. We propose to construct an unbiased estimator of (9) using the debiasing technique from the recent works [32, 42], and then use this unbiased estimator within such MCMC methods. This is the second major contribution of this work. In particular, we constrain our attention here to the SGLD.

The rest of the paper is organized as follows. In Section 3 the SMC sampler is introduced, which will be used as a crucial component within all subsequent algorithms. In particular, it will be shown how to construct consistent estimators of $\Pi(\varphi|\lambda)$ and unbiased estimators of $Q(\varphi|\lambda)$. In Section 4 we introduce the Robbins Monro algorithm, and present the first main result regarding its implementation using the SMC sampler. Section 5 concerns Bayesian posterior UQ in the context of finitely many observations. In Section 6 we simulate data from a toy model, and illustrate the algorithm from Section 4 for computing the MaxEnt solution using both exact moments and observed moments (where it coincides with the MLE), as well as the algorithm from Section 5 for UQ in the case of observations.

3. SMC samplers for estimation of integrals

Here it will be described how to estimate expectations $\Pi(\varphi|\lambda):=\sum_{x\in E^{d}}\varphi(x)\pi(x|\lambda)$, for bounded functions $\varphi:E^{d}\rightarrow\mathbb{R}$, and construct non-negative unbiased estimators of $Q(\varphi|\lambda)=\sum_{x\in E^{d}}\varphi(x)q(x|\lambda)$, for a given value of $\lambda$, using an auxillary variable technique.

Let $\eta_{0}=\Gamma_{0}$ be a known probability distribution over a general state space $U$, and let $\eta_{j}$ be some probability over $U$ for $j=0,\dots,J$, whose density can be evaluated up to a normalizing constant, but cannot be sampled from directly. Define the unnormalized measure associated to $\eta_{j}$ by $\Gamma_{j}$, with density $\gamma_{j}$, so that

[TABLE]

where $\Gamma_{j}(1)=\int_{U}\gamma_{j}(x)dx$.

The objective is to approximate expectations for bounded functions $\varphi:U\rightarrow\mathbb{R}$

[TABLE]

We cannot sample from this distribution, but we can obtain a convergent estimator

[TABLE]

using SMC samplers, as described below [15] (see also [26, 36, 10]).

Remark 3.1**.**

In the context of this paper the target $\eta_{J}=\Pi(\cdot|\lambda)$ will be considered in each inner iteration of the stochastic iterative algorithm. So $\gamma_{J}(x)=q(x|\lambda)$. To make the presentation more concrete, let $\gamma_{0}\propto 1$ and

[TABLE]

However, note that the choice of intermediate targets $\gamma_{j}$ to interpolate between an initial density $\gamma_{0}$, which can be sampled from exactly, and $\gamma_{J}$ are essentially arbitrary. For example, one could incorporate several constraints per iteration or anneal in between individual constraints.

3.1. Estimating the target

Define

[TABLE]

Let $M_{j}$ denote a Markov kernel such that

[TABLE]

Iterate Algorithm 1 and define

[TABLE]

Assumption 3.1**.**

There exists some $C>0$ such that for all $x\in U$

[TABLE]

Proposition 3.1**.**

Assume 3.1 holds. Then for $N>1$, $p\geq 1$,

[TABLE]

The constant $C(j)$ can be made to be independent of $j$ if an additional strong mixing condition is assumed on $M_{j}$, uniformly in $j$ [13].

3.2. Estimating the normalizing constant

Recall $\Gamma_{j}(1)=\sum_{x\in U}\gamma_{j}(x)$, and observe that

[TABLE]

Noting that $\Gamma_{0}(1)=\eta_{0}(1)=1$, therefore

[TABLE]

This is the normalizing constant of $\eta_{J}$.

Observe that using Algorithm 1 we can construct an estimator by

[TABLE]

Recall that for any $\varphi:U\rightarrow\mathbb{R}$ we have $\Gamma_{J}(\varphi)=\Gamma_{J}(1)\eta_{J}(\varphi)$ by definition (13). Define

[TABLE]

The following proposition is proven in [13]. The proof is reproduced in Appendix A for completeness, due to its importance to the results of this paper.

Proposition 3.2**.**

Assume 3.1 holds. Then $\mathbb{E}\Gamma_{J}^{N}(\varphi)=\Gamma_{J}(\varphi)$.

4. Stochastic approximation

Here it will be assumed that we are given values $m=\Pi(\phi|\lambda)\in\mathbb{R}^{J}$, i.e. Problem (P(P1)). As mentioned in Section 2, the MaxEnt approach entails solving the convex optimization problem (5) (which has domain $\{x;p(x)>0\}$) subject to the constraints $\sum_{x\in E^{d}}p(x)=1$, and (4). It is necessary and sufficient to solve the equation (7), reproduced below

[TABLE]

As mentioned above, a naive computation of the left-hand side of the above equation for a single value of $\lambda$ would incur a cost of $\mathcal{O}(2^{d})$, which makes the problem intractable for very small values of $d>15$ or so. It is however possible to obtain a noisy estimate using a stochastic simulation algorithm. Algorithms such as MCMC [23, 18] are designed to sample the state space in $\mathcal{O}(d)$ steps [44], and provide a consistent estimate of (17) for a given value of $\lambda$.

Returning to (17), observe that solving this equation is equivalent to solving $Q(m-\phi|\lambda)=0$. Assume we can simulate unbiased estimates $\widehat{F}(\lambda)$, i.e. estimates such that $\mathbb{E}\widehat{F}(\lambda)=Q(m-\phi|\lambda)$. The Robbins Monro algorithm [43] for stochastic fixed point iteration for finding a root of the equation $Q(m-\phi|\lambda)=0$ takes the form

[TABLE]

Under technical conditions, and for appropriate choice of $\delta_{k}$ such that $\sum_{k}\delta_{k}=\infty$ and $\sum_{k}\delta_{k}^{2}<\infty$, it can be shown that this algorithm converges to a fixed point of (17) [43, 30] without ever solving the equation exactly, even once.

We make the following assumption here, which is sufficient to ensure Assumption 3.1

Assumption 4.1**.**

There exists a $c>0$ such that for all $\lambda\in\Xi$ and $x\in E^{d}$

[TABLE]

Proposition 4.1**.**

Assume 4.1. Let $\widehat{F}(\lambda)$ be such that $\mathbb{E}\widehat{F}(\lambda)=Q(m-\phi|\lambda)$. If $\{\delta_{k}\}_{k=1}^{\infty}$ is chosen so that $\sum_{k}\delta_{k}=\infty$ and $\sum_{k}\delta_{k}^{2}<\infty$, then $\lambda^{k}\rightarrow\lambda^{*}$, where $\lambda^{k}$ is given by (18) and $\lambda^{*}$ is the solution to $Q(m-\phi|\lambda^{*})=0$, hence (17).

Consider now running the SMC sampler Algorithm 1 with the choice $\gamma_{0}\propto 1$ and $\gamma_{j}=\exp\{\sum_{i=1}^{j}\lambda_{i}\phi_{i}\}$ to construct the unbiased estimator (16) $\Gamma_{J}^{N}(m-\phi|\lambda)$ such that $\mathbb{E}\Gamma_{J}^{N}(m-\phi|\lambda)=Q(m-\phi|\lambda)$. With this unbiased estimator of $Q(m-\phi|\lambda)$, we can use Robbins Monro to solve the fixed point problem for $\lambda$ (17), hence the MaxEnt target $\pi(x|\lambda)$. The following theorem makes this statement precise.

Theorem 4.1** (Stochastic approximation with SMC).**

Assume 4.1 holds, and let $\lambda^{k}$ be generated by the iteration (18), where $\widehat{F}(\lambda):=\Gamma_{J}^{N}(m-\phi|\lambda)$ is independently generated from an SMC sampler algorithm at each iteration, as defined in (16). If $\{\delta_{k}\}_{k=1}^{\infty}$ is chosen so that $\sum_{k}\delta_{k}=\infty$ and $\sum_{k}\delta_{k}^{2}<\infty$ then $\lambda^{k}\rightarrow\lambda^{*}$, where $\lambda^{*}$ is a solution of $Q(m-\phi|\lambda^{*})=0$, hence (17).

Proof.

By Proposition 3.2 $\mathbb{E}\widehat{F}(\lambda)=Q(m-\phi|\lambda)$. The conclusion is an immediate consequence of Proposition 4.1. ∎

Assume the evaluation of $\phi_{j}$ is $\mathcal{O}(1)$. Then the evaluation of $\gamma_{j}$ is $\mathcal{O}(j)$, and since there will be $J$ steps, the total cost of the algorithm will be $\mathcal{O}(J^{2})$. Since $\lambda\in\mathbb{R}^{J}$, one may require up to $\mathcal{O}(J)$ steps of Robbins Monro. Therefore, a conservative estimate on the cost is $\mathcal{O}(J^{3})$. Recall that $J$ is the number of constraints. For $p^{\rm th}$ moments, there are $J={{d+p-1}\choose{p}}\lesssim d^{p}$ constraints, yielding a conservative estimate on the final cost of $\mathcal{O}(d^{3p})$. The constant should be rather small however, and this is of course a vast improvement on $\mathcal{O}(2^{d})$ for even the best linear algorithm for solution with exact evaluation of (17).

5. Uncertainty Quantification

In the previous section we arrive at a single optimal value $\lambda^{*}$, and a single MaxEnt target likelihood, i.e. uncertainty is not quantified. Here we venture to quantify uncertainty by putting a prior on $\lambda$, and computing (sampling from) the posterior $\mathbb{P}(\lambda|Y)$ (12). This leads to a so-called doubly-intractable target, as described above. In this section a new method is introduced to sample from the target. In particular, in subsection 5.1 we introduce the debiasing method mentioned above, and describe how to construct unbiased (signed) estimators of $\Pi(\varphi|\lambda)$ using estimator(s) from the SMC sampler. In subsection 5.2 we introduce the SGLD algorithm and present the second main result of the paper, regarding its implementation using the debiased estimators. The section is concluded with a description in subsection 5.3 of how to separate an explicit error distribution, which is more expensive and is left to future work.

5.1. Debiasing

Some MCMC methods, such as SGLD [49], as well as some recently developed piecewise-deterministic Markov processes (PMDP) [12, 4, 6, 39], require only an unbiased estimator of the gradient of the log-likelihood for implementation. Others, such as pseudo-marginal MCMC, require an unbiased and non-negative estimator of the likelihood itself. It is generally much easier to construct unbiased estimators which are signed [32, 42]. The work [31] shows how one can utilize the pseudo-marginal MCMC even with a signed estimator. However, one may easily argue that it is more natural to use a signed and unbiased estimator within SGLD or one of the PDMP algorithms. We consider the latter strategy here.

Suppose that we have a collection of random variables $\{\Delta_{l}\}_{l=0}^{\infty}$, where

[TABLE]

and we are interested in estimating $\mathbb{E}Z$, where

[TABLE]

Note that (19) guarantees that $\mathbb{E}Z<\infty$. If we draw $L\sim{\bf p}$, where ${\bf p}=[p_{0},\dots]$, then it is clear that

[TABLE]

is an unbiased estimator of $\mathbb{E}Z$. Furthermore, if we instead draw $L\sim{\bf\tilde{p}}$, where ${\bf\tilde{p}}$ is defined such that $p_{\ell}=\sum_{l\geq\ell}\tilde{p}_{l}$, then one can see that

[TABLE]

is again an unbiased estimator:

[TABLE]

The exchange of the infinite sum and expectation is allowed by (19) and Fubini theorem.

Assumption 5.1**.**

Let (19) hold and assume the probability ${\bf p}$ is chosen such that $\sum_{l=0}^{\infty}\frac{\mathbb{E}\Delta_{l}^{2}}{p_{l}}<\infty$.

Proposition 5.1**.**

For either estimator $r\in\{s,t\}$, let $\widehat{Z}^{i}_{r}$ be i.i.d. draws of the estimator (21) or (22), for $i=1,\dots,N$. Under assumption 5.1, $\mathbb{E}\widehat{Z}_{r}^{i}=\mathbb{E}Z$ and we have,

[TABLE]

with $\Sigma=\sum_{l=0}^{\infty}\frac{\mathbb{E}\Delta_{l}^{2}}{p_{l}}<\infty$. The expected cost is given by

[TABLE]

where $C(\Delta_{l})$ is the cost to obtain the realization $\Delta_{l}$.

See [32, 42, 31, 46] and references therein for proof and further discussion of such approaches. In particular, this approach will allow us to transform consistent estimators to unbiased estimators for use in the algorithms mentioned above. The recent literature on PDMP is quite extensive, and so we defer consideration of these algorithms to future work and focus on SGLD, which will be described below in generality. First we illustrate how to construct such de-biased estimator of (9).

5.1.1. Debiased SMC sampler estimator

We now describe one way to construct an unbiased estimator of (9) directly, using this de-biasing trick, without multiplying through by $Q(1|\lambda)$. In particular, the single term estimator (21) will be considered.

Let $N_{l}=N_{0}2^{2l}$. Replace step (ii) from the SMC samplers algorithm with $\hat{x}_{j}^{i}=x_{j-1}^{i}$, i.e. omit the resampling step. This algorithm was named annealed importance sampling (AIS) in [36]. The estimator $\frac{1}{N}\sum_{i=1}^{N}\prod_{j=1}^{J}G_{j-1}(x_{j-1}^{i})\varphi(x_{J}^{i})$ is an unbiased estimator of $Q(\varphi|\lambda)$, and the individual terms are independent: $x_{j}^{i}\perp x_{j}^{i^{\prime}}$ for $i\neq i^{\prime}$ and all $j=1,\dots,J$. This means in particular that if we construct

[TABLE]

from two independent replications, one with $N$ total particles and one with $N/4$ total particles, then two estimators have the same distribution; in particular, the same expectation. So we have the following proposition.

Proposition 5.2**.**

Define the coupled estimator $\Delta_{l}:=\tilde{\eta}^{N_{l}}_{J}(\varphi)-\tilde{\eta}^{N_{l}/4}_{J}(\varphi)$, with particles coming from the same realization of the algorithm described above with $N_{l}$ particles. Then $\mathbb{E}\widehat{Z}_{s}=\eta_{J}(\varphi)$, with $\widehat{Z}_{s}$ given by (21). The factor $4$ here relates to the specific choice of $N_{l}$ and can be replaced with $N_{l+1}/N_{l}$ in general.

Proof.

Following from the explanation above, the equation (20) takes the form $\eta_{J}(\varphi)=\sum_{l=1}^{\infty}\mathbb{E}\Delta_{l}$, so the result follows from Proposition 5.1. ∎

Observing that in this case $\mathbb{E}\Delta_{l}^{2}\propto N_{l}^{-1}$ and $C(\Delta_{l})\propto N_{l}$, we can observe that this is the so-called sub-canonical case [42]. That is, if the rate of convergence of $\mathbb{E}\Delta_{l}^{2}$ were faster (or rate of growth of expected cost slower), then we could choose $p_{l}$ in such a way that both asymptotic variance and cost were finite. As this is not the case, we instead choose $p_{l}\propto N_{l}^{-1}l\log(l)$, as proposed in [42] (see also [16]), which results in finite variance and expected cost which is finite with high probability.

5.2. SGLD

Given a symmetric positive definite matrix $P(\lambda)\in\mathbb{R}^{J\times J}$ the following SDE keeps $p(\lambda)$ invariant [19, 21]

[TABLE]

where for $j=1,\dots,J$

[TABLE]

For constant $P(\lambda)$, it was proposed in [49] to replace $\nabla\log p(\lambda)$ by an unbiased estimator $\widehat{\nabla\log p(\lambda)}$ (i.e. $\mathbb{E}\widehat{\nabla\log p(\lambda)}=\nabla\log p(\lambda)$) and approximate the SDE using Euler Maruyama method with stepsize $\delta_{n}\propto 1/n$, similar to Robbins Monro, as follows

[TABLE]

For the moment, assume $P(\lambda)=1$, and $\widehat{\nabla\log p(\lambda)}=\nabla\log p(\lambda)$, and define the measure corresponding to $p(\lambda)$ by $\mu$. Furthermore, define the empirical measure corresponding to $K$ steps of the constant $\delta$-discretized version of (27) by $\mu^{K}_{\delta}(\varphi)=\frac{1}{K}\sum_{n=1}^{K}\varphi(\lambda^{n})$ and define the limiting diffusion (25) sampled at the same $K$ points $i\delta$, $i=1,\dots,K$, by $\mu^{K}$. Then we have

[TABLE]

The discretization error $\delta^{2}$ comes from the standard Euler-Maruyama strong convergence rate for SDE with constant diffusion. Define the auto-covariance

[TABLE]

where $\lambda_{\infty}^{n}$ is a sample from (25) at time $\delta n$ and $\bar{\lambda}_{\infty}=\int\lambda\mu(d\lambda)$. One expects $\rho_{n}\approx(1-\delta)^{n}$. Therefore the integrated auto-covariance IACT$=2\sum_{n=0}^{\infty}\rho_{n}\approx\delta^{-1}$, and the sample error is asymptotically (after burn-in) proportional to IACT$/K\approx(\delta K)^{-1}$. Balancing the error terms, one finds that $\delta\propto K^{-1/3}$ is optimal, leading to MSE of $K^{-2/3}$. It is slightly worse than the MSE of $K^{-1}$ for $K$ steps of standard MCMC methods. This simple argument suggests that the optimal choice for $\delta_{n}$ is $\delta_{n}\propto n^{-1/3}$, which is not the same as for SGD. In fact, this result is made rigorous in the work [45, 47] (see also references therein), leading to the following proposition.

Proposition 5.3**.**

Let $P(\lambda)=1$, define the measure corresponding to $p(\lambda)$ by $\mu$, choose $\delta_{n}\propto n^{-1/3}$, and let $\lambda^{n}$ be generated by (27). Then for suitably regular $\varphi$, one has

[TABLE]

Furthermore, this choice of $\delta_{n}$ is asymptotically optimal in terms of cost, and for $K$ steps, one has MSE $\propto K^{-2/3}$.

Let $p(\lambda)=\mathbb{P}(\lambda|Y)\propto\mathbb{P}(Y|\lambda)\mathbb{P}(\lambda)$, and recall (9), so that we have

[TABLE]

Assume that $\nabla\log\mathbb{P}(\lambda)$ can be computed exactly. Now define an unbiased estimator of (30) using the method of Proposition 5.2, i.e. draw $L\sim{\bf p}$, and let

[TABLE]

Theorem 5.1**.**

Let $P(\lambda)=1$, define the measure corresponding to $p(\lambda)=\mathbb{P}(\lambda|Y)$ by $\mu$, choose $\delta_{n}\propto n^{-1/3}$, and let $\lambda^{n}$ be generated by (27), with $\widehat{\nabla\log p(\lambda^{n})}$ defined as in (31). Then for suitably regular $\varphi$, one has

[TABLE]

Furthermore, this choice of $\delta_{n}$ is asymptotically optimal in terms of cost, and for $K$ steps, one has MSE $\propto K^{-2/3}$.

Proof.

This follows directly from Propositions 5.2 and 5.3. ∎

5.2.1. Reimann manifold method

It is assumed in [38] that Proposition 5.3 holds also for deterministic evaluable $P(\lambda)$, yielding a Reimann manifold SGLD (RMSGLD), but this is not proven. In any case, based on the above heuristic one would expect a different scaling of $\delta_{n}$ to be optimal, since the strong rate of convergence is reduced to $\delta$, so balancing terms gives $\delta\propto K^{-1/2}$ and an MSE of $K^{-1/2}$. Nevertheless, here we further postulate that (32) holds indeed for the case in which we have an unbiased estimator $\widehat{g}(\lambda)\approx g(\lambda)$ and a non-negative unbiased estimator $\widehat{P}(\lambda)\approx P(\lambda)$.

Let $P(\lambda)=Q(1|\lambda)=\Gamma_{J}(1|\lambda)$, and let the target be $p(\lambda)=\mathbb{P}(\lambda|Y)\propto\mathbb{P}(Y|\lambda)\mathbb{P}(\lambda)$, so that

[TABLE]

where ${\bf 1}\in\mathbb{R}^{J}$ denotes the vector of all $1$s, and

[TABLE]

where $\widehat{m}=\frac{1}{M}\sum_{k=1}^{M}\phi(y^{(k)})$. Then using the notations of Section 3, (26) takes the form

[TABLE]

This can be unbiasedly approximated by

[TABLE]

as can $P(\lambda)$ by simply $\widehat{P}(\lambda)=\Gamma_{J}^{N}(1|\lambda)$.

This approach would allow direct simulation without de-biasing, since only unbiased estimates of the unnormalized measure $Q$ are required, which can be constructed directly from SMC samplers, as described in Proposition 3.2. Nonetheless, this is not pursued further here due to the expected higher cost of an optimal implementation. We note that the cost of these algorithms can be improved using the multilevel Monte Carlo method [17], but this is deferred to future work.

5.3. Separating an explicit error distribution

Before concluding the section, we mention an extension of the model considered in the rest of the paper which provides deeper UQ. Suppose that we do not have i.i.d. observations $y^{(i)}\sim\pi(\cdot|\lambda)$, but in fact we have observations of the form

[TABLE]

where b^{\epsilon}=[b_{1}^{\epsilon},\dots,b_{d}^{\epsilon}]^{T},\leavevmode\nobreak\ b_{i}^{\epsilon}\sim B(\epsilon)\leavevmode\nobreak\ {\rm i.i.d.}\leavevmode\nobreak\ for $i=1,\dots,d$, and $B(\epsilon)$ denotes a Bernoulli with probability $\epsilon$. In other words, each $b_{i}^{\epsilon}$ is independently $1$ with probability $\epsilon$ and [math] otherwise. This model explicitly accounts for error in the observations, i.e. with probability $\epsilon$ any qubit may have been observed wrong.

We have $\mathbb{P}(d|y,\lambda)=\epsilon^{|d-y|}(1-\epsilon)^{1-|d-y|}$. The likelihood of a single observation $d$ is given by

[TABLE]

and so we can see that again we have an integral with respect to $\pi(\cdot|\lambda)$. Let $D=[d^{(1)},\dots,d^{(M)}]$. Due to the independence between observations and errors, we have

[TABLE]

where we defined $\varphi_{d^{(i)}}(y):=\epsilon^{|d^{(i)}-y|}(1-\epsilon)^{1-|d^{(i)}-y|}$.

Therefore

[TABLE]

One can define

[TABLE]

An unbiased estimator of this function can be constructed from M independent SMC samplers, one for each factor. Despite being massively parallel, this is a significant computational burden, particularly for large $M$. Investigation of efficient algorithms for solving this problem is the topic of ongoing research.

6. Simulations

6.1. Maximum likelihood: exact marginals

Recall the form of $\widehat{F}(\lambda)$, for a fixed $\lambda$ from (4.1). Using the machinery from subsection 3.2 we could construct such an unbiased estimator, for any sequence of intermediate targets $\gamma_{j}(x|\lambda)$, and any $J$, provided that $\gamma_{J}(x|\lambda)=q(x|\lambda)$. The natural choice is of course to let

[TABLE]

as above, and this choice will be used here. The MCMC kernels $M_{j}$ are chosen as Metropolis-Hastings kernels with the reversible random walk proposal $q(\lambda,\lambda^{\prime})$ defined by $\lambda^{\prime}_{i}={\rm mod}(\lambda_{i}+1,2)$ with probability $\beta$ and $\lambda^{\prime}_{i}=\lambda_{i}$ otherwise, for $i=1,\dots,J$. In other words, $q(\lambda,\lambda^{\prime})=|\lambda-\lambda^{\prime}|$ is a vector of independent Bernoullis with probability $\beta$, and is hence symmetric, with tunable parameter $\beta$ determining the step-size. We note that a Gibbs sampler can also be used here. Now we let $\delta_{n}=1/n$, and iterate (18) until convergence.

In order to test the accuracy of our methodology, we simulate the data as follows. Choose $a_{ij}\in\mathbb{R}$, for $i=1,\dots,d$ and $j=i,\dots,d$, and let $A_{ij}=a_{ij}$ for $j\geq i$ and $A_{ij}=a_{ji}$ otherwise. Now consider the true target

[TABLE]

Similarly, let $\Lambda$ be defined as $\Lambda_{ij}=\lambda_{ij}$ for $j\geq i$ and $\Lambda_{ij}=\lambda_{ji}$ otherwise, for coefficients $\{\lambda_{ij}\}_{\{i=1,\dots,d,\leavevmode\nobreak\ j=i,\dots,d\}}$. We will identify $a,\lambda$ with vectors in $\mathbb{R}^{J}$, where $J=d(d+1)/2$. Notice this is of MaxEnt form with $\phi(x)={\rm vec}[\{x_{i}x_{j}\}_{\{i=1,\dots,d,\leavevmode\nobreak\ j=i,\dots,d\}}]$, i.e. first and second moments are observed. We choose $d=10$, let $N=2d$, take $d$ MCMC steps with $\beta=0.6$, and define $\delta_{n}=\delta^{1}_{n}1_{n\leq 2d}+\delta^{2}_{n}1_{n>2d}$, where $\delta^{1}_{n}=\epsilon n_{0}/(n_{0}+n)$, $\delta^{2}_{n}=\Gamma^{N}_{J}(m-\phi|\lambda^{2d})^{-1}\epsilon n_{0}/(n_{0}+n)$, $n_{0}=5d$, and $\epsilon=1$. For $n\leq 2d$ the direction of descent is taken as $m-\eta_{J}^{N}(\phi)$, while for $n>2d$ the full estimator (16) is used $\widehat{F}(\lambda^{n})=\Gamma_{J}^{N}(m-\phi|\lambda^{n})$, so that this is an unbiased estimator of $Q(m-\phi|\lambda)$. This construction prevents large fluctuations in $\Gamma^{N}_{J}(1|\lambda^{n})$ during early outer iterations from precluding stabilization of the algorithm. The truth, reconstruction, and convergence plot are given in Fig. 1.

6.2. Maximum likelihood: observations

In this section we explore maximum likelihood estimation with finitely many observations and empirical moments. We look at the cases of $M=1000$ and $M=50$ observations. The results are shown in Figure 2. The top two panels correspond to $M=1000$ (left) and $M=50$ (right). The good news is that we obtain the MLE, up to the bias arising due to the observational error, quite rapidly: within several hundred iterations for $M=1000$ and several tens for $M=50$. This is observed in the middle two panels. Notice from the plot for $M=50$(right) that the misfit error eventually starts to increase once the algorithm begins to fit noise. This is typical and can be expected. When the algorithm is stopped the error is quite significant. In the top two panels, we observe that the reconstruction for $M=50$ does not look similar to the true $A$ (shown in the left panel of Figure 1), however for $M=1000$ it looks quite acceptable. Just for a sanity check, observe the bottom two panels which show the observed moments for $M=50$ on the left and predicted moments using the reconstruction (top right) on the right. The agreement is quite decent, and reasonable given the noise level of the observations. The point is that the error in the observations translates to a much larger error in the parameters $\lambda$.

6.3. Uncertainty Quantification using de-biasing

Here we consider sampling from the posterior, with density $\mathbb{P}(\lambda|Y)$. The non-informative improper prior $\mathbb{P}(\lambda)\propto 1$ is adopted.

As an initial experiment we ran the SGD algorithm using

[TABLE]

with $L\sim{\bf p}$ as described in subsection 5.1, the numerator constructed as in (24), and ${\bf p}$ chosen as described in subsection 5.1.1, with the same prototype problem except with $d=4$, and with exact moments. The results are shown in Figure 3 left. The middle and right panels show the results of running the algorithm with the original estimator, and the simple biased but consistent estimator $\eta^{N}_{J}(m-\phi|\lambda)$, respectively. Note the algorithm with the biased but consistent drift appears to be convergent also. As mentioned earlier, this is done in practice, but we are not aware of theoretical results verifying this. Therefore we will proceed with the SGLD using this unbiased estimator, i.e. we iterate (27) with $P=1$ and $\delta_{n}=n^{-1/3}$, and drift given by (37), with $L\sim\bf{p}$.

The full SGLD MCMC is run with $M=1000$ observations on the $d=4$ qubit system. It takes roughly two hours on a laptop to obtain $K=10^{6}$ samples. The pairwise marginals on the diagonal ${\rm diag}(\Lambda)$ are illustrated in Figure 4. The histograms are constructed by discarding the first half of the samples for burn-in, and resampling the remaining samples according to their weights $\{\delta_{n}\}$. The true value is also indicated in red in the plots. While the truth is indeed in the region of high probability, and the mean is reasonably close to the truth, the spread of the posterior is still quite significant, and more than one would hope for with $M=1000$ observations. This is consistent with the results of Figure 2, where it is observed that the error in the coefficients is amplified quite a bit in comparison to the error in the moments. It is however not feasible to compute the pushforward distribution from the coefficients to the moments here, since such estimation for a single moment would require a separate large sample simulation (e.g. by MCMC or SMC) for each of the $K=10^{6}$ samples. We look at the same simulation with $M=10^{6}$ observations, and the results are plotted in Figure 5. The spread is much tighter, and also there is a strong correlation between the diagonal elements in this case.

7. Summary and Future Directions

A first attempt has been made to estimate the MaxEnt distribution associated to outputs of a quantum computer, and quantify the posterior uncertainty. This task has lead to a novel use of the SMC sampler for construction of an unbiased estimator of the intractable drift for use in a Robbins Monro stochastic approximation algorithm, which can be used to efficiently compute the MaxEnt distribution for exact moments or the MLE distribution for observed moments. The Bayesian formulation yields to a doubly-intractable target. For its solution we use de-biased SMC sampler estimators of the log-likelihood within a SGLD algorithm to construct consistent posterior expectations. All the algorithms are provably convergent, and numerical simulations are provided using (classically) simulated data from a toy model. This exercise reveals that the posterior uncertainty in the distribution can be significantly amplified with respect to the uncertainty in the moments arising from having finitely many observations. Topics for further research include: (0) the extension of the current framework to density matrices, (i) the incorporation of an explicit error distribution for the observations within the model, (ii) exploration of more complex targets, i.e. more qubits, (iii) MLMC acceleration of SGLD, (iv) investigation of PDMP MCMC methods as alternatives to SGLD, (v) other debiasing strategies and opportunites for acceleration, and of course (vi) investigation of outer problems, such as model selection.

Acknowledgements: This work is supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research (ASCR) quantum algorithm teams program, under field work proposal number ERKJ333.

Appendix A

The proof of Proposition 3.2 is given below.

Proof.

Define

[TABLE]

and observe that the iterates of the algorithm of section 3.1 can be rewritten in the concise form

[TABLE]

where we recall the definition (15) for the empirical measure $\eta^{N}_{j-1}$.

Notice that

[TABLE]

By the law of iterated expectations we have

[TABLE]

Iterating in this way, it is clear that

[TABLE]

Noting that $\Gamma_{0}=\eta_{0}$, the right-hand side is exactly (39).

∎