A Appendix A: Derivation of Cramer-Rao Bounds

# How to best sample a periodic probability distribution, or on the accuracy of Hamiltonian finding strategies

## Abstract

Projective measurements of a single two-level quantum mechanical system (a qubit) evolving under a time-independent Hamiltonian produce a probability distribution that is periodic in the evolution time. The period of this distribution is an important parameter in the Hamiltonian. Here, we explore how to design experiments so as to minimize error in the estimation of this parameter. While it has been shown that useful results may be obtained by minimizing the risk incurred by each experiment, such an approach is computationally intractable in general. Here, we motivate and derive heuristic strategies for experiment design that enjoy the same exponential scaling as fully optimized strategies. We then discuss generalizations to the case of finite relaxation times, .

Introduction. Measurement adaptive tomography has recently been suggested as an efficient means of performing partial quantum process tomography (11); (5). Little is known about optimal protocols when realistic experimental restrictions are imposed — as opposed to the case where one is allowed arbitrary quantum resources1. Indeed, even in the simplest examples, not even bounds have been given on the proposed protocols. Here, we give analytic bounds on both non-adaptive and adaptive estimation protocols for a Hamiltonian parameter estimation problem. Moreover, we derive estimation protocols which asymptotically achieve these bounds. Adaptive protocols are typically difficult to implement because a complex optimization problem must be solved after each measurement. We instead derive a heuristic that is easy to implement and achieves the exponentially improved asymptotic risk scaling of the optimal solution.

Within the nuclear magnetic resonance (NMR) community, similar concerns have motivated the examination of the use of maximum entropy (1) and maximum likelihood (3) methods for obtaining spectra. Recently, computational power has become available such as to make these methods feasible for use in analyzing non-uniform data obtained from high-dimensional NMR experiments (8). These studies have produced qualitatively similar strategies for how to best design experiments when each sample is expensive to collect.

The paper is organized as follows. First, we define the model Hamiltonian which we want to estimate the parameters of, along with our metric of success. Then we give both frequentist and Bayesian lower bounds on the risk derived from this metric. Finally, we derive strategies which achieve the asymptotic scaling of these bounds.

Problem statement. The model we consider is a qubit evolving under the Hamiltonian

 H=ω2σz.

Here is the unknown parameter whose value we want to ascertain. We make the problem dimensionless by assuming . An experiment consists of preparing a single known input state , evolving under the Hamiltonian for a controllable time and performing a measurement in the basis. We emphasize here that we are assuming strong projective measurements on individual copies of a quantum preparation, rather than weak measurements on physical ensembles such as those studied in NMR experiments.

The outcomes of the measurement we label , where and refer to and , respectively. An experiment design consists of a specification of the time that we evolve a qubit under before we measure. The likelihood function for a given experiment is then given by the Born rule and . Using our model Hamiltonian, we can express the likelihood more simply as:

 Pr(d|ω,t)=sin2(ω2t)dcos2(ω2t)1−d. (1)

Note that this model does not include noise. Below, we somewhat generalize this model by including limited visibility and a dephasing process.

If we desire an estimate of the true value , a commonly used figure of merit is the squared error loss:

 L(ω,^ω)=|ω−^ω|2.

The risk of an estimator, which is a function that takes data sets to estimates , is its expected performance with respect to the loss function:

 R(ω,^ω)=∑DPr(D|ω,T)L(ω,^ω(D,T)).

For squared error loss, the risk is also called the mean squared error (MSE).

Mean squared error lower bound. The difficulty here is that the random outcomes of the measurements are not identically distributed. In fact, since they depend on the measurement time, each one could be different. Although, asymptotic results exist for non-identically distributed random variables2, these results are derived for insufficient statistics, such as the sample mean. Moreover, we desire to provide computationally tractable heuristics that permit useful estimates with a finite number of samples.

Although it is quite difficult to obtain exact expressions for the risk for arbitrary measurement times, in some cases we have obtained an asymptotically tight lower bound. For unbiased estimators, we can appeal to the Cramer-Rao bound (4)

 R(ω,^ω)≥1I(ω), (2)

where

 I(ω)=−∑DPr(D|ω,T)∂2log(Pr(D|ω,T))∂ω2 (3)

is called the Fisher information. In our particular case, the Fisher information reduces to quite a simple form in

 I(ω)=N∑k=1t2k, (4)

which is conveniently independent of (a derivation is given in Appendix A). Thus, the mean squared error is lower bounded by

 R(ω,^ω)≥1∑Nk=1t2k. (5)

Later we show that this bound becomes exponentially suppressed when we include noise in our model. In general, this quantity is dependent on the true parameter .

The Bayesian solution considers the average of the risk, called the Bayes risk, with respect to some prior :

 r(π,^ω)=∫R(ω,^ω)π(ω)dω.

As in references (5); (11), we choose a uniform prior for . Then, the final figure of merit is the average mean squared error:

 r(^ω)=∫R(ω,^ω)dω.

The goal is to find a strategy which minimizes this quantity. Although there exist Bayesian generalizations of the Cramer-Rao bound (6), ours is independent of and thus remains unchanged by integrating equation (5) over the parameter space:

 r(^ω)≥1∑Nk=1t2k. (6)

Note also that, in general, Bayesian Cramer-Rao bounds require fewer assumptions to derive than the standard (frequentist) bound. Although they are the same for this model, they differ for a more general model considered later. In broad strokes, the difference in practice between Bayesian and frequentist methods is averaging versus optimization. Below we demonstrate a heuristic strategy which draws from both methods to achieve the goal of determining the measurement times which give the lowest possible achievable bound on the Bayes risk (6).

Looseness of the Cramer-Rao bound. As useful as the Bayesian Cramer-Rao lower bound (6) is, it is simple to see that it is not always achievable. We can obtain a lower bound by considering the best protocol we could possibly hope for in any two-outcome experiment. In such a protocol, one bit of experimental data provides exactly one bit of certainty about the parameter . If we learn the bits of in sequence, at each step , our risk is upper bounded by the worst-case where all the remaining bits of are either all 0 or all 1. In either case, the error incurred by estimating a point between the two extremes is given by , leading to the best possible MSE after measurements being , even though we can make a smaller Cramer-Rao bound by choosing times that grow faster than this exponential function. Note that this risk is achievable via the standard phase estimation protocol (2), but that this protocol requires quantum resources which are not part of our model.

Examples. Let us consider a couple of examples for which the lower bound can be further simplified. First, consider the case when all the measurement times are the same. This is by far the simplest case, since the outcomes become identically distributed. Recall . Then, the measurement time should be less then the first Nyquist time, , or the data will be consistent with more than one . That is, for (but less than , say), the likelihood function will have two equally likely maxima. We minimize the risk, then, by choosing . Then, the maximum likelihood estimator (MLE), for example, will be asymptotically efficient (9) achieving the Cramer-Rao lower bound

 r(^ωMLE)=2π2N+O(1N2).

Now consider a uniform grid of times. Since , we should choose the Nyquist sampling rate: . Then, for any estimator using data collected at these measurement times, the Cramer-Rao bound gives

 Missing or unrecognized delimiter for \left

Again, the maximum likelihood estimator will be asymptotically efficient. However, since the likelihood function will have many local maxima, the maximum likelihood estimator is non-trivial to find as gradient methods are not guaranteed to work. Bayesian estimators were derived in (11), where simulations yielded risk scaling which is asymptotically efficient.

Note that since we are considering a uniform spacing of times, we can apply a Fourier estimation technique without worrying about spectral aliasing introduced by non-uniformity (10). That is, we apply the discrete Fourier transform and estimate the peak of the power spectrum. Since the resolution in the frequency domain is , we expect the Bayes risk to be

 r(^ωFourier)=1π2N2.

The sampling theorem requires that we sample from a deterministic function, not a probability distribution. In practice, this condition is often approximately satisfied by sampling some stable statistic such as the mean value of the distribution at each time. This can be achieved by measuring at the same time until a sufficiently accurate estimate of the mean at that time is obtained, then repeating this for many other times. But as we have shown, this method can be quadratically improved by performing every single measurement at a different time.

Exponentially achievable lower bound. It has been shown that Bayesian adaptive solutions lead to risk decreasing exponentially with the number of measurements (11). However, these results are given by fits to numerical data. Here, we give an analytic lower bound on the risk of these protocols.

The local (in time) Bayesian adaptive protocol can be described as follows: (1) begin with a uniform prior and determine the first measurement time which minimizes the average (over the two possible outcomes) variance of the posterior distribution; (2) perform a measurement at , record the outcome , and update the distribution via Bayes’ rule; (3) repeat step (1) replacing the current prior with the current posterior. Note that the expected variance in the posterior is the Bayes risk. Thus, the protocol attempts to minimize the risk assuming the next measurement is the last. Strategies that are local in this sense are called a greedy strategies, as opposed to strategies which attempt to minimize the risk over all future experiments.

For some choices of measurement times, including those given by the protocol above, the posterior will be approximately normally distributed3. This is guaranteed in the asymptotic limit, but the posterior distribution near its peak is also remarkably well approximated by a Gaussian after as few as 15 reasonably chosen measurements (we found a uniform grid to be sufficient for “warming up” to the Gaussian approximation). Thus, we approximate the current distribution (at given some sufficiently long measurement record ) as

 Pr(ω|D)=1√2πσ2e−(ω−μ)22σ2,

with some arbitrary mean and variance implied by . The expected posterior variance (which is equal to the Bayes risk) of the probability distribution of the next measurement is

 r(t)=σ2(1+t2σ2sin(μt)2−et2σ2+cos(μt)2), (7)

(derived in Appendix B) which oscillates with frequency within an envelope . Asymptotically, the minimum risk will approach the minimum of the envelope for all , but will be a lower bound on the risk otherwise. This minimum occurs at with a risk of , which is also the variance of the updated probability distribution since both outcomes are equally probable at . Thus, at each measurement step we reduce the risk by . Thus, the risk scales exponentially as and is achieved at measurement times which scale as

 tk∼1σ(1−e−1)k/2≈1.26kσ.

These times are guaranteed to be optimal only in the asymptotic limit. For finite numbers of samples, we suggest two simple heuristics. First, we suggest the use of exponentially increasing times, where the base of the exponent is optimized offline, followed by the use of the maximum likelihood estimator for these times. Second, we suggest a simpler adaptive scheme based on the assumption that the distribution remains Gaussian after each measurement. Making use of this normality assumption, we only need update equations for the mean and variance of the distribution over . In deriving the update equations, we also take into account the oscillations of the expected Bayes risk by finding the nearest achievable minima to the one given by the lower bound. We provide the update equations in Appendix C.

Generalization to finite . In practice, we will have to consider not only experimental restrictions but also noise and relaxation processes. Processes which do not affect the quantum state can be effectively modeled by random bit-flip errors occurring with probability . Processes which do affect the quantum state (decoherence) are modeled by an exponential decay of phase coherence4 with characteristic time . Since the state being measured lies in the -plane of the Bloch sphere, this loss of phase coherence manifests as an exponential decaying envelope being applied to the original likelihood (1). The model is thus fully specified by the likelihood function

 Pr(0|ω,t,η,T2)=η⎛⎜⎝e−tT2cos2(ω2t)+1−e−tT22⎞⎟⎠+1−η2. (8)

The Cramer-Rao bound is now given by

 R(ω,^ω)≥⎛⎜ ⎜⎝N∑k=1t2kη2sin2(ωtk)e2tkT2−η2cos2(ωtk)⎞⎟ ⎟⎠−1. (9)

Note that unlike the Cramer-Rao bound (5) for the noiseless case, the above bound is not independent of and thus we must appeal to the Bayesian Cramer-Rao bound so that the measurement times can be chosen independently of the true parameter. However, the Bayesian bound turns out to be very loose. A sharper bound is given by first upper bounding each term in the denominator to give

 r(^ω)≥1η2∑Nk=1t2ke−2tkT2.

The noise term (or visibility) simply gives a constant reduction in the achievable accuracy. The relaxation process provides a more interesting dynamic as we see that the gains from longer times are exponentially suppressed. In other words, strategies are restricted to explore . We can thus do no better than

 r(^ω)≥e2Nη2T22. (10)

The adaptive strategy discussed above can be generalized to include noise and relaxation but the expressions are more lengthy (see Appendix B). To illustrate the performance of our adaptive strategy, we simulate the adaptive strategy along with offline strategies using identical times (), linearly spaced times () and exponentially sparse times (). For each strategy, we perform simulations for experiments consisting of different numbers of samples , up to , and repeat each such simulation to obtain an estimate of the Bayes risk for that strategy and experiment size. In Fig. 1, we present the results of these simulations for the noiseless case, and for the cases and .

Note that in all cases, the adaptive strategy achieves exponential scaling until the times selected reach . At that point, the risk will then scale linearly if the remaining measurement times are . However, if the protocol continues to select larger measurement times, the information gained from those measurements will tend to zero and the risk will remain constant.

Summary and conclusions. By using the Cramer-Rao bound along with analytic expressions for the variance of each posterior distribution, we have motivated a heuristic method for choosing experiment designs that asymptotically admits exponentially small error scaling in the number of measurements. For finite measurements, we have relied on numerical simulation to demonstrate that this scaling is well-achieved even for . Numerical simulations for finite , moreover, have suggested that we can enjoy exponential scaling of the risk until the measurement times saturate the bound, at which point the risk scaling switches to the asymptotic scaling of . In both cases, the heuristics used to design experiments are quite computationally tractable, thus motivating the utility of our heuristics to actual experimental practice.

Acknowledgements. We thank Miriam Diamond for assistance in testing and developing the simulation software. CF thanks Josh Combes for helpful discussions. This work was financially supported by NSERC and CERC.

## Appendix A Appendix A: Derivation of Cramer-Rao Bounds

In this Appendix, we show that for the simple model represented by the likelihood function presented in equation (1), the Fisher information given by (3) reduces to the form claimed in (4). To show this, we first note that the likelihood for a vector of observations at times is given by a product of the likelihoods for each individual measurement,

 Pr(D|ω,T)=∏kPr(dk|ω,tk).

Thus, the log-likelihood function is simply a sum over the individual log-likelihoods. Since the derivative operator commutes with summation, we obtain that

 ∂2∂ω2logPr(D|ω,T)=∑k∂2∂ω2logPr(dk|ω,tk).

This in turn implies that the Fisher information for a vector of measurements is given by the sum for each measurement of that measurement’s Fisher information.

To calculate the single-measurement Fisher information, we find the second derivative of the log-likelihood for a single measurement is given by

 ∂2∂ω2logPr(dk|ω,tk)=t2k(2dk−1)(1−2dk+cos(ωtk))((2dk−1)cos(ωtk)−1)2.

Thus, we find that the single-measurement Fisher information is given by

 I(ω|tk) =−∑dk∈{0,1}Pr(dk|ω,tk)∂2∂ω2logPr(dk|ω,tk) Missing or unrecognized delimiter for \left =t2k.

We conclude that , as claimed.

For the model with finite and limited visibility, given by the likelihood function (8), we can follow the same logic. We find the second derivative of (8) with respect to gives us

 ∂2∂ω2logPr(dk|ω,tk)=ηt2k⋅(2dk−1)(η(1−2dk)+etkT2cos(ωtk))(η(1−2dk)cos(ωtk)+etkT2)2.

The expected value of this derivative then gives us the Fisher information for a single measurement in the finite- model,

 I(ω|tk)=η2t2ksin2(ωtk)e2tkT2−η2cos2(ωtk).

Taking the sum of this information then produces the Cramer-Rao bound given in (9).

## Appendix B Appendix B: Asymptotic Scaling of the Bayes Risk

In this Appendix, we derive expressions for posterior distributions under the assumption of a normally-distributed prior, and then apply these expressions to show the asymptotic scaling of the Bayes risk. We also derive update rules that allow for expedient implementation of the greedy algorithm described in the main text.

Under the assumption of a normally-distributed prior, all prior information about the parameter can be characterized by the mean and variance of the prior distribution. Thus, we shall write our priors as to reflect the assumption of normality. Then, the probability of obtaining a datum at time given such prior information is then given by

 Pr(d|t;μ,σ2)=∫∞−∞Pr(d|t,ω)Pr(ω|μ,σ2)dω=14(2−(2d−1)(1+e2iμt)e−12t(σ2t+2iμ)).

Applying Bayes’ rule then produces the posterior distribution

 Pr(ω|d,t;μ,σ2)=Pr(ω|μ,σ2)Pr(d|t,ω)Pr(d|t;μ,σ2)=√2πe−(μ−ω)22σ2((1−2d)cos(tω)+1)σ(2−(2d−1)(1+e2iμt)e−12t(σ2t+2iμ)).

The mean and variance of this distribution are given by:

 E[ω|d,t;μ,σ2] =2((2d−1)e−12σ2t2(σ2tsin(μt)−μcos(μt))+μ)2−(2d−1)(1+e2iμt)e−12t(σ2t+2iμ) V[ω|d,t;μ,σ2] =μ2+σ2−2((2d−1)e−12σ2t2(σ2tsin(μt)−μcos(μt))+μ)2−(2d−1)(1+e2iμt)e−12t(σ2t+2iμ) −2(2d−1)σ2teiμt(σ2tcos(μt)+2μsin(μt))(2d−1)(1+e2iμt)−2e12t(σ2t+2iμ)

To chose optimal times, we wish to pick so as to minimize the expected value over of the variance, where this expectation is taken over possible data. Based on the previous expressions, we find that

 Missing or unrecognized delimiter for \left

in agreement with Equation (7).

This expected variance, which describes our risk incurred by measuring at a given , is bounded below by an envelope . A pair of examples of the envelope and achievable risk is illustrated in Figure 2.

Note that the envelope is minimized by . Moreover, the expected variance saturates the lower bound at intervals in of , but the width of the envelope’s minimum grows as , so that as more measurements are performed, the bound becomes a good approximation for the minimum achievable risk. Thus, in the asymptotic limit of large numbers of experiments, we have that the risk at scales with each step as the minimum of the envelope,

 E(^t,σ2)σ2=1−e−1≈0.632.

We conclude that in the asymptotic limit, the risk decays as , where is the number of measurements performed.

## Appendix C Appendix C: Update Equations for μ, σ2

In this Appendix, we state without derivation the update rules for and after obtaining a measurement result from an experiment performed at time , under the assumption of an normal prior. For the simple model described by Equation (1),

 E[ω|d] =μ−π(2d−1)σ2(−1)k(2k−1)exp(−π2σ2(1−2k)28μ2)2μ (11) V[ω|d] =σ2−π2(1−2d)2σ4(1−2k)2exp(−π2σ2(1−2k)24μ2)4μ2, (12)

where is used to pick the intersection of and to the minimum of , as described in Appendix B.

For the finite- model, the updated mean and variance are given by

 E[ω|d] =μ+π(2d−1)(−1)k(2k−1)σ2exp(−(π−2πk)(−2πkσ2T2+4μ+πσ2T2)8μ2T2)2μ (13) V[ω|d] =σ2−π2(2d−1)2(2k−1)2σ4exp(−(π−2πk)(−2πkσ2T2+4μ+πσ2T2)4μ2T2)4μ2, (14)

where in this case,

 k=round⎡⎢ ⎢⎣μ−μ√4σ2T22+1+πσ2T22πσ2T2⎤⎥ ⎥⎦.

### Footnotes

1. As in the standard phase estimation protocol. See e.g. (2).
2. The frequentist reference is (7), while a useful Bayesian reference is (12).
3. This is true asymptotically and higher order corrections can be used if required (12).
4. We do not include amplitude damping in our model since our populations remain equal throughout evolution and thus only manifests as a contribution to .

### References

1. J. C. J. Barna, E. D. Laue, M. R. Mayger, J. Skilling, and S. J. P. Worrall. Exponential sampling, an alternative method for sampling in two-dimensional NMR experiments. Journal of Magnetic Resonance (1969), 73(1):69–77, June 1987.
2. Andrew M. Childs, John Preskill, and Joseph Renes. Quantum information and precision measurement. Journal of Modern Optics, 47(2):155–176, 2000.
3. Roger A. Chylla and John L. Markley. Theory and application of the maximum likelihood principle to NMR parameter estimation of multidimensional NMR data. Journal of Biomolecular NMR, 5, April 1995.
4. Thomas M. Cover and Joy A. Thomas. Elements of Information Theory. Wiley-Interscience, second edition, July 2006.
5. Christopher Ferrie, Christopher E Granade, and D. G Cory. Adaptive Hamiltonian estimation using Bayesian experimental design. AIP Conf. Proc. 1443:165–173, 2011.
6. Richard D. Gill and Boris Y. Levit. Applications of the van Trees Inequality: A Bayesian Cramér-Rao Bound. Bernoulli, 1(1/2), 1995.
7. Bruce Hoadley. Asymptotic Properties of Maximum Likelihood Estimators for the Independent Not Identically Distributed Case. The Annals of Mathematical Statistics, 42(6), 1971.
8. Sven G. Hyberts, Haribabu Arthanari, and Gerhard Wagner. Applications of Non-Uniform sampling and processing. Topics in Current Chemistry, July 2011. PMID: 21796515.
9. E. L. Lehmann and George Casella. Theory of Point Estimation. Springer, second edition, September 1998.
10. Mark W. Maciejewski, Harry Z. Qui, Iulian Rujan, Mehdi Mobli, and Jeffrey C. Hoch. Nonuniform sampling and spectral aliasing. Journal of Magnetic Resonance, 199(1):88–93, July 2009.
11. Alexandr Sergeevich, Anushya Chandran, Joshua Combes, Stephen D. Bartlett, and Howard M. Wiseman. Characterization of a qubit Hamiltonian using adaptive measurements in a fixed basis. Physical Review A, 84:052315, 2011.
12. Ruby C. Weng. A Bayesian Edgeworth expansion by Stein’s Identity. Bayesian Analysis, 5:741–764, 2010.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters