Stochastic Estimation of Dynamical Variables

TL;DR

The paper introduces STEADY, a robust stochastic estimation method for accurately determining the Hamiltonian of quantum systems, leveraging all measurement data and scalable to information-theoretic limits.

Contribution

It presents a novel stochastic gradient descent-based estimator that efficiently and robustly estimates complete Hamiltonians of quantum systems, including non-unitary parameters.

Findings

Performance scales at information-theoretic limits

Robust to state preparation and measurement errors

Applicable to piecewise-differentiable Hamiltonians

Abstract

Estimating the parameters governing the dynamics of a system is a prerequisite for its optimal control. We present a simple but powerful method that we call STEADY, for STochastic Estimation algorithm for DYnamical variables, to estimate the Hamiltonian (or Lindbladian) governing a quantum system of a few qubits. STEADY makes efficient use of all measurements and its performance scales as the information-theoretic limits for such an estimator. Importantly, it is inherently robust to state preparation and measurement errors. It is not limited to evaluating only a fixed set of possible gates, rather it estimates the complete Hamiltonian of the system. The estimator is applicable to any Hamiltonian that can be written as a piecewise-differentiable function and it can easily include estimators for the non-unitary parameters as well. At the heart of our approach is a stochastic gradient descent over the difference between experimental measurement and model prediction.

Equations

$$|\dot{\psi }(t)⟩=-iH\left(\mathbf{d}(t)\right)|\psi (t)⟩.$$

$${\tilde{H}}_{ij}(\mathbf{\sigma },\mathbf{h};\mathbf{d})=h_{ij}+\stackrel{[}{k}=1]D\sum {\sigma }_{ijk}{d}_k,$$

$$\begin{array}[]{c}\tilde{H}(\bm{\alpha},\bm{\beta};\bm{d})=\stackrel{{\scriptstyle[}}{{k}}=1]{M}{\sum}a_{k}A_{k},\\ \text{where }a_{k}=\stackrel{{\scriptstyle[}}{{l}}=1]{D}{\sum}\alpha_{kl}d_{l}+\beta_{k}.\end{array}$$

$$|{\psi }_i⟩={\mathrm{e}}^{-iH({\mathbf{d}}_i)T}|0⟩.$$

$$C(\mathbf{\omega })=\frac{1}{P}\sum _{i=1}^P\text{dist}\left({\hat{\mathbf{p}}}_i,{\tilde{\mathbf{p}}}_i(\mathbf{\omega })\right).$$

$$\hat{\mathbf{\omega }}=\arg \underset{\mathbf{\omega }}{\min }C(\mathbf{\omega }).$$

$$V(\mathbf{\omega })=\frac{1}{P_v}\sum _{i=1}^{P_v}\text{dist}\left({\mathbf{p}}_i,{\tilde{\mathbf{p}}}_i\left(\mathbf{\omega }\right)\right),$$

$$\text{var}({\omega }_l)\ge \frac{1}{P\times S\times ⟨{\mathcal{I}}_l^{\text{samp}}⟩}.$$

$$\mathbb{E}\left(({\omega }_l-{\omega }_{l,0}{)}^2\right)\ge \frac{{\left(1-\frac{\partial b_l}{\partial {\omega }_l}\right)}^2}{P\times S\times ⟨{\mathcal{I}}_l^{\text{samp}}⟩}+b_l^2,$$

$$H(\mathbf{d})=$$

$$+$$

$$\text{where }{\delta }_i={\kappa }_i{d}_i$$

$$\mathbf{G}(\mathbf{\omega })\approx \mathbf{G}({\mathbf{\omega }}_0)+\mathbf{H}({\mathbf{\omega }}_0).(\mathbf{\omega }-{\mathbf{\omega }}_0),$$

$$\hat{\mathbf{\omega }}\approx {\mathbf{\omega }}_0-\mathbf{H}({\mathbf{\omega }}_0{)}^{-1}.\mathbf{G}({\mathbf{\omega }}_0).$$

$$\mathbb{E}\left(\left({\hat{\omega }}_i-{\omega }_{0i}\right)\left({\hat{\omega }}_j-{\omega }_{0j}\right)\right)\approx$$

$$\text{where }{\mathcal{I}}_{\mathbf{C}}({\mathbf{\omega }}_0)=$$

$$A_{ij}=$$

$$B_{ij}=$$

$$l(\omega )=\log \left(\prod _{i=1}^P\left(S!\prod _{k=1}^{2^Q}\frac{{\tilde{p}}_{ik}(\mathbf{\omega }{)}^{{\hat{p}}_{ik}S}}{({\hat{p}}_{ik}S)!}\right)\right),$$

$${{\mathcal{I}}_{\mathbf{F}}}_{ij}=$$

$$=$$

$$\approx$$

$$C(\mathbf{\omega })+{\lambda }_k\Vert \mathbf{\omega }{\Vert }_1$$

$${\rho }_{\text{init}}=(1-Qs)|0⟩⟨0|+s\sum _{q=1}^Q|0\dots 01_{q}0\dots 0⟩⟨0\dots 01_{q}0\dots 0|,$$

$$\mathcal{S}=(1-Qs)I+{\mathcal{S}}^{\prime }$$

$$C_s(\mathbf{\omega })$$

$$=\frac{1}{P}\sum _{i=1}^P\text{dist}\left({\hat{\mathbf{p}}}_i+s{\mathcal{S}}^{\prime }{\hat{\mathbf{p}}}_i+\mathcal{O}(s^2),{\tilde{\mathbf{p}}}_i(\mathbf{\omega })\right)$$

$$=C(\mathbf{\omega })+\frac{1}{P}\sum _{i=1}^{P}s{\mathcal{S}}^{\prime }{\hat{\mathbf{p}}}_i.\left({\hat{\mathbf{p}}}_i-{\tilde{\mathbf{p}}}_i(\mathbf{\omega })\right)$$

$$+\mathcal{O}(s^2).$$

$$C_s({\mathbf{\omega }}_0+\mathbf{b})=$$

$$+\frac{1}{P}\sum _{i=1}^{P}s{\mathcal{S}}^{\prime }{\hat{\mathbf{p}}}_i.\left({\hat{\mathbf{p}}}_i-{\tilde{\mathbf{p}}}_i(\mathbf{\omega })\right)$$

$$+\delta s^2+\text{higher orders in s or ∣b∣},$$

$$C_s({\mathbf{\omega }}_0+\mathbf{b})$$

$$=$$

Full text

A common task in physics and engineering is the control of a system, where the control pulses sent to the system pass through a complex transfer function before they effect a useful change to the state of the system. There are two overarching prerequisites for good control: learning the dynamical law that governs the system (the goal of disciplines like experimental design and parameter estimation) and, consecutively, the derivation of control pulses for the given system (broadly covered by optimal control theory). Advances in these areas are crucial for applications in quantum information science, where the precise control of well-characterized quantum systems will form the basis for quantum computers.

Here we present STEADY, a conceptually simple but performant method for approaching the parameter estimation problem for dynamical variables. We can model a piece of quantum hardware with a Hamiltonian (or Lindbladian) $\tilde{H}(\bm{\omega};\bm{d})$ which depends on the parameters to be estimated $\bm{\omega}$ and on the control pulses $\bm{d}(t)$. Our goal becomes finding the value for $\bm{\omega}$ that leads to an $\tilde{H}$ that (for any value of the control pulses $\bm{d}$) most closely mimics the dynamical law $H$ governing the real hardware. As is commonly done in parameter estimation, we do this by searching for a value of $\bm{\omega}$ that minimizes some measure of distance between $\tilde{H}$ and $H$.

Our contribution follows in the rich traditions of stochastic methods and compressed sensing: instead of performing full process tomography on the hardware which would be extremely time consuming, we run a relatively small number of random control pulses on it and study its response. For each control pulse we sample the final state of the system (for instance by projectively measuring the qubits in the computational basis). We then estimate the difference between this experimental measurement and the prediction based on the $\tilde{H}({\bm{\omega}})$ model. This measure of “difference” is stochastic, as it uses only a small finite sample of possible control drives.

This leads to a number of properties that make STEADY perform particularly well. First, the distance measure that we use is differentiable with respect to $\omega$ which lets us use efficient (stochastic) gradient descent to rapidly find optimal values for the parameters being estimated. Moreover, the stochastic nature of our estimator leads to much lower resource requirements. We avoid doing full tomography, which greatly reduces the number of necessary measurements, while only modestly increasing the number of steps required by the (now stochastic) gradient descent. Simultaneously, the stochasticity lets us surpass the error floor otherwise imposed by the finite number of measurements performed when sampling the final states of the system. By using pulse sequences of varying lengths, our method becomes inherently insensitive to state preparation and measurement (SPAM) errors. The fact that dynamical variables—Hamiltonians and Lindbladians—are local quantities means that our estimates are generally sparse descriptions of the noise in a system, in contrast to an estimate that reconstructs a finite-time evolution of the noise. This allows regularized estimators to be used that avoid overfitting and lead to good estimates with surprisingly few data. In fact, our estimator can be restated as a maximum likelihood estimator which naturally approaches the information-theoretic limit of the Cramér-Rao bound.

Background

In the field of quantum computing there is a long history of achievements in both parameter estimation and control theory: the design of precise control schemes Peirce et al. (1988); Dawson and Nielsen (2005); Krotov and Feldman (1983); Khaneja et al. (2005); Caneva et al. (2011) for the preparation of states Tannor and Rice (1985); Law and Eberly (1996), unitary operations Unanyan et al. (1998); Mischuck and Mølmer (2013); Heeres et al. (2015, 2017a); Krastanov et al. (2015), and even complete quantum channels Bacon et al. (2001); Shen et al. (2017) has permitted advances in metrology, chemistry, communication, and—more recently—quantum computation. These control schemes are, in return, informed by the precise description of the system, obtained through various tomographic measurements Leibfried et al. (1996); Chuang and Nielsen (1997); Poyatos et al. (1997).

One can simplify the control problem by considering only a discrete set of gates instead of studying the continuous control that a complete knowledge of the Hamiltonian would provide. As long as the available set of basic gates generates the unitary group, there is a known efficient compilation procedure Dawson and Nielsen (2005); Parent et al. (2017); Amy et al. (2017); Svore et al. (2017). The performance bottleneck is in estimating the exact behavior of the given gates. Techniques like process tomography Leibfried et al. (1996); Chuang and Nielsen (1997); Poyatos et al. (1997) have existed for a while, but they are susceptible to state preparation and measurement (SPAM) errors Merkel et al. (2013). Gate set tomography Merkel et al. (2013); Blume-Kohout et al. (2013); Greenbaum (2015) mostly circumvents issues of SPAM, by requiring the preparation of only one type of state (e.g. the ground state) and only one type of measurement (e.g. in the computational basis). However gate set tomography is still susceptible to what Greenbaum (2015) calls “intrinsic SPAM errors”: state preparation errors for the ground state (e.g. a finite temperature of the system) and any imperfections in the projective measurement in the computational basis. Finally, there are tools like randomized benchmarking Emerson et al. (2005); Knill et al. (2008) that quantify average error rates of quantum processes (instead of the entire channel) without SPAM, and so can be used only as a benchmark, not as a tool for direct calibration (though extensions of the idea do allow this Kimmel et al. (2014)).

On the other end of the spectrum are control schemes that compute new control pulses for every single unitary operation (instead of compiling them out of the predetermined “universal” set of gates). This type of “continuous” control provides for quantum circuits with an order of magnitude smaller depth Heeres et al. (2017b); Boixo et al. (2018); Neill et al. (2018), however, it is also computationally more difficult. A popular approach to it is the use of gradient-based methods like GRAPE Khaneja et al. (2005). The fidelity of a given operation (the difference between the desired operation and the operation actually implemented by the control pulse) is computed as a function of the control pulse parameters. The fidelity is a differentiable function and its gradient with respect to the control pulse parameters is also computed, thus permitting efficient gradient descent, leading to a locally optimal control pulse. A large body of work is available discussing how to avoid getting stuck in local optima Brif et al. (2010); Zahedinejad et al. (2015); Palittapongarnpim et al. (2016). Recent development in the use of reinforcement learning has even provided for gradient-free techniques robust to noise Chen et al. (2014); Bukov et al. (2017); Niu et al. (2018).

However, most of the gradient-based control techniques require precise knowledge of the Hamiltonian in order to provide high-fidelity control drives. Some adaptive techniques get around this problem by switching to an in-situ method when they are near the optimal pulse. When the imprecisions in the model Hamiltonian start dominating, they forsake the model and start measuring the pulse fidelity experimentally through process tomography Egger and Wilhelm (2014); Wu et al. (2018) or randomized benchmarking Kelly et al. (2014). However, this makes the optimization much less rapid as the gradient is not directly available anymore and techniques like downhill simplex become necessary. Such optimization techniques are limited by the statistical error in the fidelity estimation. Recently, an elegant workaround based on “simultaneous perturbation stochastic approximation” was suggested in the ACRONYM method Ferrie and Moussa (2015), breaking through this statistical error floor. Impressive improvements have been seen in experiments following variations of these methods Rol et al. (2017). In either case, however, many additional experimental samples are required, but are then discarded after the current iteration of the optimizer. Moreover, all of these additional measurements are done for the sake of designing one specific gate with exquisite precision, but they do not contribute to estimating the Hamiltonian of the system and are forsaken when later on one tries to design another gate.

Here we focus on the precise estimation of the Hamiltonian itself, which can later be used in any control scheme. We suggest STEADY, a simple, but powerful approach that exploits the entirety of the measurement data in the Hamiltonian estimation, reaching the fidelity limits imposed by information theory. Similarly to ACRONYM, we use stochastic techniques to surpass the statistical error floor. Similarly to gate set tomography we are inherently insensitive to SPAM errors. Moreover, borrowing ideas from randomized benchmarking, our method can circumvent even the intrinsic SPAM errors (like finite temperature in the preparation of the ground state). The random pulses, informationally incomplete measurements, sparse models, and regularized estimators that we use take advantage of compressed sensing methods for state and process tomography Gross et al. (2010); Flammia et al. (2012) to improve accuracy. Lastly, STEADY estimates the complete Hamiltonian or Lindbladian (or other dynamical models, e.g. a stochastic master equation) of a system, not just one gate or a set of gates.

A recent independent preprint (Flurin et al., 2018) by Flurin et al. presents techniques similar to STEADY’s with significant differences in the design choices. While the prescription for gathering experimental data and the comparison of the model to that data are very similar in both approaches, Flurin et al. use a recurrent neural net (RNN) to model the dynamics, while we use physical models with very general parametrization. Flurin’s black box approach is promising for models that are particularly difficult to differentiate, like stochastic master equations, however modern autodifferentiation frameworks (Innes, 2018) enable the use of STEADY as well. Furthermore, the RNN size is expected to grow exponentially with the number of qubits in order to simulate the quantum dynamics. That exponential cost is explicitly present in STEADY. It would be interesting to examine whether the RRN, after hyperparameter optimization, would find a sparse representation of the dynamics, similarly to STEADY’s use of Hamiltonian and Lindbladian generators which are explicitly a sparser representation of the otherwise dense superoperator.

In the following we specify the formalism we use to describe the protocol and demonstrate that it reaches the information theoretical limit in estimation fidelity. We discuss the effect of the intrinsic SPAM errors and how to circumvent them. Experimental design techniques that further improve the fidelity of our estimator are described. We briefly discuss the effects of parameter drift, non unitary errors, and nonlinearities in the Hamiltonian (as a function of the control pulses). Together with this manuscript we also provide a software package based on a popular differentiable programming framework Abadi et al. (2015) that implements our techniques for various models including unitary or non-unitary evolution.

Problem statement

A system of $Q$ qubits (a $2^{Q}$-dimensional Hilbert space) is controlled by a Hamiltonian $H(\bm{d})$ (that is itself a function of time-dependent control pulse $\bm{d}(t)$ set by the experimentalist). $\bm{d}$ is a $D$-dimensional real vector, where $D$ is the number of control parameters available to the experimentalist. The evolution of an initial state $|\psi\rangle$ will then be expressed as

[TABLE]

To accurately predict this dynamics we must learn the Hamiltonian $H$ as a map $H:\bm{d}\rightarrow H(\bm{d})$ in order to be able to control the quantum hardware using a control pulse $\bm{d}$. We introduce a parameterized model for the Hamiltonian, $\tilde{H}(\bm{\omega};\bm{d})$, in which case the problem becomes finding the values of all parameters in the array $\bm{\omega_{0}}$ for which $H(\bm{d})=\tilde{H}(\bm{\omega_{0}};\bm{d})$ for all $\bm{d}$. We will also discuss the case where the model $\tilde{H}$ cannot exactly represent the reality of $H$, as well as cases where non-unitary evolution is non-negligible.

With some a priori knowledge of the physical system, an experimentalist might be able to deduce an approximation of $\bm{\omega_{0}}$, however the experimentalist could also run experiments on the hardware to learn successively better approximations of $\bm{\omega_{0}}$. An experimentalist can run a control pulse $\bm{d}(t)$ and then project and measure the final state of the system in the computational basis. The details of how the pulses are chosen and how the measurement data is used distinguishes the various approaches to Hamiltonian estimation and calibration like process tomography Leibfried et al. (1996); Poyatos et al. (1997); Chuang and Nielsen (1997), gate set tomography Blume-Kohout et al. (2013); Greenbaum (2015), and randomized benchmarking Emerson et al. (2005); Knill et al. (2008); Kimmel et al. (2014).

Unlike most estimation techniques, we work at the level of control pulses, without hiding them behind a set of precompiled gates. Moreover, our protocol is inherently untroubled by SPAM errors, as no special states or measurements are necessary, besides preparing the ground state and performing measurements in the computational basis, just like in gate set tomography. We improve even further by partially circumventing the intrinsic SPAM errors found in gate set tomography Greenbaum (2015).

Methods

To describe our protocol, we will first consider only Hamiltonians that are linear in the control pulse. This is an appropriate description when some prior knowledge of the relationship between the control fields and the Hamiltonian is known, but more general mappings could be addressed as well, as we discuss below. The model for such a Hamiltonian in its most general form would be (in index notation for $i,j\in[1..2^{Q}]$)

[TABLE]

where $h_{ij}$ and $\sigma_{ijk}$ belong to arrays of complex numbers representing the parameters $\bm{\omega}$ that need to be learned. The numbers $h_{ij}$ form the elements of a $2^{Q}\times 2^{Q}$ matrix and can be interpreted as the drift Hamiltonian of the system, and the numbers $\sigma_{ijk}$ (forming a $2^{Q}\times 2^{Q}\times D$ array) can be interpreted as the list of control Hamiltonians (one for each control parameter in $\bm{d}$, where $\bm{d}$ could be time dependent). Hermiticity can be ensured if the parametrization is done in terms of pairs symmetric and antisymmetric real matrices. There are $2^{2Q}\times(D+1)$ real parameters to be learned in this case.

Such a high level of parametrization might be unnecessary for a well studied system. In such a case one can list the fixed known Hermitian operators that are part of the Hamiltonian in a large list $\{A_{1},A_{2},\dots,A_{M}\}$ (where we have $M$ such possible operators) and only parameterize how these operators are summed together in

[TABLE]

Here $\alpha_{kl}$ and $\beta_{k}$ are arrays of real numbers representing the parameters $\bm{\omega}$ that need to be learned. The $M$-dimensional real vector $\bm{\beta}$ represents all the drift components of the Hamiltonian, while the $M\times D$ real matrix $\bm{\alpha}$ represents the linear mixing/crosstalk between drives. There are $M\times(D+1)$ real parameters to be learned in this case. For the next few paragraphs we will consider only this low-degree-of-freedom parametrization, and later on we will discuss when one might want to use the more general approach. We stress that our approach is not limited to Hamiltonians linear in the control pulse parameters—any functional dependence can be used in place of the two discussed above. In particular, we discuss the estimation of a Lindbladian later in this manuscript, and even the measurement backaction parameters of a stochastic master equation can be studied using STEADY.

Experiments in this setup proceed by performing some form of regression on data that has been gathered from the hardware in order to find an approximation of $\bm{\omega_{0}}$. To gather the data we suggest the following approach. Begin by generating a large number $P$ of random control pulses on the classical computer controlling the hardware. For simplicity we will initially consider only constant pulses of fixed duration $T$ where each random pulse is taken from a normal distribution of unit variance centered on zero, but our approach is equally easy to apply to random time-dependent pulses sampled from other distributions with fixed variance. This provides us with a list $\{\bm{d}_{1},\dots\bm{d}_{P}\}$ of control pulses. Each control pulse is run on the hardware, initialized to the ground state, resulting in a final state

[TABLE]

Experimentalists cannot exactly read the components of $|\psi_{i}\rangle$ in a given basis (e.g. the computational basis), rather they can only estimate them through projective measurements. Specifically, our protocol requires the experimentalist to run each pulse $S$ times in order to repeatedly sample through projective measurements. For each $\bm{d}_{i}$ (and corresponding $|\psi_{i}\rangle$) this provides a vector of estimated Born rule probabilities $\hat{\bm{p}}_{i}$. With an infinite number of samples $S$, and the assumption that the unknown Hamiltonian parameters do not drift, the estimate would converge to be exactly the Born rule probability vector $\bm{p}_{i}:=\lim_{S\to\infty}\hat{\bm{p}}_{i}$ (the $k$-th component of $\bm{p}_{i}$ is $\left|\langle k|\psi_{i}\rangle\right|^{2}$, where $|k\rangle$ enumerates the computational basis). Similarly, the Born rule probability for a given drive $\bm{d}_{i}$ predicted by the model Hamiltonian $\tilde{H}(\bm{\omega};\bm{d}_{i})$ will be denoted $\tilde{\bm{p}}_{i}(\bm{\omega})$ (its $k$-th component is $\left|\langle k|\mathrm{e}^{-i\tilde{H}(\bm{\omega};\bm{d}_{i})T}|0\rangle\right|^{2}$). We note that when the model contains the true Hamiltonian then $\bm{p}_{i}=\tilde{\bm{p}}_{i}(\bm{\omega_{0}})$.

We can define the “distance” between the measured estimate for the population and the predicted population, averaged over the $P$ random pulses:

[TABLE]

Our estimator $\hat{\bm{\omega}}$ for $\bm{\omega_{0}}$ is the minimum of this distance measure (called the “cost function” from here on):

[TABLE]

If the distance function $\text{dist}()$ is the cross entropy, i.e. $\text{dist}(\bm{a},\bm{b})=-\bm{a}.\log(\bm{b})$, then our estimator is a maximum-likelihood estimator. For most of the numerical examples, the distance function we use is the mean squared error $\text{dist}(\bm{a},\bm{b})=\left(\bm{a}-\bm{b}\right)^{2}$, which is simpler, but in practice leads to the same estimate. The cost function is differentiable, which permits us to run automated stochastic gradient descent optimizers in the search for $\bm{\omega_{0}}$. In practice, we also augment this cost function with a regularization cost for $\bm{\omega}$ to avoid overfitting and enhance convergence, as discussed in the supplement. Ideally, we would have the sum run over all possible control pulses, but this is unfeasible in finite time. Stochastic gradient descent, where only a small number $P$ of random pulses is used, is what enables our method thanks to its strong guarantees of convergence to the same minimum.

Results

In the next few paragraphs we study the performance of this method. In order to check our susceptibility to over-fitting or convergence failures, we also introduce a validation cost function (which would be unavailable to the experimentalist, but is available on our simulated “mock” hardware)

[TABLE]

where $P_{v}$ is the size of $\{\bm{d}_{1}^{v},\dots\bm{d}^{v}_{P_{v}}\}$, a new set of random “validation” control pulses sampled from a unit-variance distribution. This validation function does not suffer from the statistical noise inherent to finite $S$: it is a sample estimate of the cost function over pulses, and it is the expected value over measurements. Moreover $V$ is non-negative and $V(\bm{\omega_{0}})=0$, hence $V(\bm{\omega})=\mathcal{O}\left(\left(\bm{\omega}-\bm{\omega_{0}}\right)^{2}\right)$, given our choice of distance function (see supplementary materials). We will keep the validation set $\{\bm{d}_{i}^{v}\}$ the same in all comparisons, even if we change the size, variance, or anything else related to the training set $\{\bm{d}_{1},\dots\bm{d}_{P}\}$. Similarly, for the validation function we always use a pulse of unit duration, even if we decide to use different duration pulses for the cost function.

The validation function is defined so that small values of $V(\bm{\hat{\omega}})$ imply good predictive power of the empirically reconstructed model. Our validation function was chosen for this reason and for its close connection to the cost function, but in the supplementary materials we demonstrate that other more common measures of fidelity have the same scaling.

For most of this manuscript the distance function in $C$ and $V$ is the mean squared error and other choices are discussed in the appendix. Unless specified otherwise, numerical results are given for $T=1$, and a simulated system of $Q=3$ qubits, driven by Pauli drives and $\sigma_{i}^{+}\sigma_{i+1}^{-}$ nearest neighbor exchange interactions, where $\bm{\omega}$ gives the relative strengths of each drive. Details are provided in the supplementary materials.

Statistical Sampling Errors

If we could obtain a perfect estimate of the populations $\{\bm{p_{i}}\}_{i=1..P}$ (i.e. if we could have $S=\infty$), then even a small data set (a small $P$) would be sufficient to perfectly estimate the Hamiltonian parameters. The only issue would be ensuring our system is not under-constrained by having $P\gtrsim M\times(D+1)$, and regularizing the gradient descent procedure to ensure we are not stuck in a valley of the cost function. Indeed, when $S=\infty$, we rapidly converge to $C(\bm{\omega})\approx V(\bm{\omega})\approx 10^{-16}$, a floor imposed by the floating point precision.

However, in a realistic case we could run the quantum hardware only a finite number of times $P\times S$, where a finite $S$ will incur a statistical error on $\{\hat{\bm{p}}_{i}\}_{i=1..P}$. Given that obtaining $\hat{\bm{p}}_{i}$ is a multinomial sampling procedure from the distribution $\bm{p}_{i}$, we can expect an error $\propto\frac{1}{\sqrt{S}}$, which would cause an error floor of $C(\bm{\omega})\propto\frac{1}{S}$. We would need to increase $S$ in order to get a better estimate of $\bm{\omega_{0}}$, as can be seen from the behavior of $C(\bm{\omega})$ in Fig. 2a. From that figure one could think that increasing $S$ is important while increasing $P$ is a waste of resources, however, the minimized value of $C(\bm{\omega})$ is a bad proxy of the quality of our parameter estimation given that it is inherently plagued by the $\frac{1}{S}$ statistical error. This is why we have introduced the validation function V(\bm{\omega})=\mathcal{O}\bigl{(}(\bm{\omega}-\bm{\omega_{0}})^{2}\bigr{)}, which in Fig. 2b. shows that the total amount of information $P\times S$ is the important resource expended in parameter estimation. The precision of our estimate scales as $\frac{1}{P\times S}$, i.e. inversely with the total amount of measurements we take from the hardware. It is inconsequential how we group the data (more pulses or better estimation of the result from fewer pulses) as long as $P\times S$ is kept constant and $P$ is sufficiently large to constrain the system. In fact, as long as we have good regularization that ensures convergence of the gradient descent procedure, even $S=1$ (where $\hat{\bm{p}}_{i}$ becomes a binary vector) performs just as well.

On first sight it can be counter intuitive that the validation function continues to improve even when the actual cost function reaches a floor, but this is similar to the difference between the standard deviation of a distribution (taking the role of $C$) and the possibly much lower standard error in the estimator of the mean of that distribution (taking the role of $V$). Such behavior is typical of stochastic optimizers and can be seen in ACRONYM Ferrie and Moussa (2015) as well. A linearized example of this can be seen in the supplementary material. A more rigorous understanding of this effect can be presented in terms of the Cramér–Rao bound Cramér (2016). The variance of each component $\omega_{l}$ of our estimator for $\bm{\omega_{0}}$ is bounded by $\frac{1}{\mathcal{I}_{l}(\omega_{l})}$, where $\mathcal{I}_{l}$ is the corresponding component of the diagonal of the Fisher information matrix for our measurements (without loss of generality, we assume the parameter vector $\bm{\omega}$ is chosen so that no two distinct components have nonzero correlation). Moreover, $\mathcal{I}_{l}=P\times S\times\langle\mathcal{I}_{l}^{\text{samp}}\rangle$, where $\mathcal{I}_{l}^{\text{samp}}$ is the Fisher information for the estimator of a single projective measurement for a given pulse and $\langle\dots\rangle$ denotes an average over all control pulses $\{\bm{d_{1}},\dots\bm{d_{P}}\}$. The projective measurement is equivalent to sampling once the multinomial distribution described by $\bm{p}_{i}$, hence $\mathcal{I}_{l}^{\text{samp}}=\stackrel{{\scriptstyle[}}{{k}}=1]{2^{Q}}{\sum}p_{k}\left(\frac{\partial p_{k}}{\partial\omega_{l}}\right)^{2}$, where $p_{k}=\left|\langle k|\psi_{i}\rangle\right|^{2}$. For the purposes of statistical errors we can already observe that for every component $\omega_{l}$ of $\bm{\omega}$ we have

[TABLE]

This confirms our observations from Fig. 2, and proves the efficiency and unbiasedness of our estimator.

Intrinsic State Preparation and Measurement Errors

The discussion from the preceding paragraphs did not consider the effects of imperfect state preparation and measurements (SPAM). Unlike with process tomography, we do not need to prepare initial states spanning the whole Hilbert space, nor do we need to make projections in anything but the computational basis. As such, STEADY is not susceptible to the usual SPAM errors, an advantage we share with gate set tomography. However, both gate set tomography and STEADY suffer if the ground state is not properly cooled or if the measurement in the computational basis is not perfect (effects called “intrinsic SPAM” by Greenbaum (2015)). Below we study the error floor caused by the intrinsic SPAM and describe how our protocol is able to deal with it. Our model for intrinsic SPAM in the numerical examples is parameterized by a single parameter $s$, which is defined as the probability that any given qubit might have flipped from $|0\rangle$ to $|1\rangle$ during state preparation. We also choose the probability that a projective measurement is incorrectly reported as its opposite to be this same value $s$. For numerical simplicity second order effects are neglected, i.e. two qubits cannot flip at the same time (during preparation or measurement). The detailed exact expression for the intrinsic SPAM model we use is also given in the supplementary materials.

As one can see from Fig. 3(b), increasing the intrinsic SPAM probability (the per-qubit error rate for initial state preparation or final projective measurement) leads to a breakdown of the $\frac{1}{P\times S}$ scaling and an intrinsic error floor that cannot be surpassed by simply increasing the available data. The Cramér-Rao bound is not immediately useful in explaining this effect, because we have initially phrased it in terms of an unbiased estimator for $\bm{\omega_{0}}$. However, the intrinsic SPAM errors contribute to a bias in our estimator, in which case the bound becomes

[TABLE]

where $\mathbb{E}(\dots)$ is an expectation value over all possible sets of measurements, and $\bm{b}$ is the bias of the estimator of $\bm{\omega_{0}}$.

In Fig. 3(d) we see the appearance of an error floor in the fidelity/validation function. We observe $V\propto s^{2}$, which unsurprisingly points to a bias that grows with the intrinsic SPAM error $\left|\bm{b}\right|\propto s$ (behavior that can be observed by leading order expansions of the validation function, as done in the supplementary materials).

An inspiration from randomized benchmarking leads to a way to surpass this error floor, namely that longer control pulses would suffer from fixed SPAM errors, but greater and greater coherent bias-induced errors. Anywhere that $\bm{b}$ appears, it is multiplied by $T$, leading to $\mathrm{e}^{-i\tilde{H}(\bm{\omega}+\bm{b};\bm{d})T}=\mathrm{e}^{-i\tilde{H}(\bm{\omega};\bm{d})T}+\mathcal{O}\left(\left|\bm{b}\right|T\right)$. Doing the same Taylor expansion as described above leads us to $\left|\bm{b}\right|\propto\frac{s}{T}$. We do observe in Fig. 3(c) that we rapidly return to the much lower statistical error by using longer pulses. Equivalently, instead of using longer pulses, we can use pulses of greater average power. Thanks to this, STEADY provides for addressing not only statistical errors (by taking more samples), but also intrinsic SPAM errors (by using longer/stronger pulses).

A figure with a more exhaustive numerical report over various pulse lengths and SPAM error rates is presented in the supplementary materials.

Another way to fight the effects of SPAM is to include them in the model of the dynamics, and estimate the parameters governing SPAM together with the unitary parameters. Later in the manuscript we discuss such approaches of extending the model to describe more general dynamics like non-linear drives and non-unitary effects.

Optimal Control and Experimental Design

The usual goal of estimating the parameters describing one’s hardware would be to permit high-fidelity open-loop control. There is a rich history of methods designed for that purpose Law and Eberly (1996); Mischuck and Mølmer (2013); Krastanov et al. (2015); Heeres et al. (2017a); Leung et al. (2017); Mavadia et al. (2017). Lately, optimal control techniques have been used with great success. The “differentiable programming” toolkits Abadi et al. (2015) that have enabled the rapid development of STEADY is readily applicable to the reverse problem of optimal control Leung et al. (2017): we “freeze” the model parameters and now optimize with respect to the control drives, while the minimization target is not the distance between a measured and predicted state, but rather between the desired and predicted state.

However, in the context of our work, there is a more exciting application of optimal control, that would permit parameter estimation at much lower resource/time cost. As we have described in previous sections, the fidelity of any estimator is limited by the Fisher information contained in our measurements. Here we suggest a relatively straightforward optimization procedure to increase the available Fisher information. This procedure will perform well in practice, and as we have shown, it continues to saturate the scaling behavior of the Cramér–Rao Bound, but it is not globally optimal. As our task is a multiparameter estimation problem and parameters that couple to noncommuting observables cannot be measured simultaneously without disturbance, the optimal estimation scheme is unlikely to be as simple as the use of random control pulses in STEADY (we refer the reader to Refs. (Pezzè et al., 2017; Matsumoto, 2002; Yang et al., 2018) for an overview of the literature on optimal multiparameter estimation). Instead of using random control pulses we can also run gradient descent to find the control pulses that maximize the determinant of the Fisher information of the measurements (also known as D-optimal design in the field of experimental design). This optimization procedure does not need to converge with great precision to the true maximum, as we are interested in the gross gain obtained from switching from random pulses to optimized pulses (the minor gain from precisely finding the maximum is negligible in this context). Fig 4 demonstrates how STEADY can halve the required number of measurements while obtaining the same fidelity, thanks to this careful design of control pulses. We have kept the average power of the control pulses fixed to a single unit, both in the case of random and optimized pulses. In our mock tests this technique has been able to provide similar gains in the fidelity of the estimator over a large range of values for $P$ or $S$.

Discussion

There are a number of technical details and caveats that need to be addressed.

First, the stochastic gradient descent could have trouble reaching the optimum due to difficult-to-traverse valleys in the cost landscape. In practice, annealing the learning rate in an adaptive Nesterov-momentum Adam optimizer Nesterov (1983); Kingma and Ba (2014), together with annealing of a 1-norm regularizer (following MacKay’s “empirical Bayes” trick MacKay (1992)) was sufficient for robust performance of our estimator. These choices were informed by a hyper-parameter optimization study. We provide more details about this in the supplement.

Many hardware systems could be plagued by slow drifts in the Hamiltonian. Approaches to solving this issue span from estimating the drift (i.e. just making it one of the parameters describing the dynamics) to repeating the estimation procedure at regular intervals. Both techniques are readily applicable to our method, additionally opening the possibility to use efficient transfer learning Pratt (1993) if necessary.

It is important to also note that our approach for circumventing the intrinsic SPAM errors assumes that our model is capable of exactly describing the actual dynamics of the system. There are two cases in which this might not be true. On one hand, we might have coherent errors, due to small corrections to $H$ that are not described by the map $\bm{d}\mapsto\tilde{H}(\bm{d})$. Such problems can be straightforwardly addressed by parameterizations of $\tilde{H}$ that include quadratic and higher terms.

A special type of error that can never be captured by the Hamiltonian formulation above would be an inherently non-unitary error such as decay or dephasing. The only way to fully capture the effects of such non-unitary dynamics is to include them in the model used by the estimator, and the most natural approach is to simulate the evolution of the system with a master equation. This leads naturally to the Lindbladian version of STEADY, which follows the same principles as the Hamiltonian version detailed above. Fig. 5 demonstrates that including the non-unitary dynamics in our model permits us to reach the statistical error floor that was unreachable with an incomplete Hamiltonian model. A more detailed discussion in terms of Fisher Information content for this Lindbladian case is provided in the appendix. Even more general dynamical laws can be implemented in STEADY as well: whether for use in classical mechanics, other dynamical systems, or in the rich field of continuous weak measurements of quantum systems (Flurin et al., 2018; Tsang, 2009; Cortez et al., 2017).

Conclusion

We have introduced STEADY, a method for Hamiltonian or Lindbladian estimation that reaches the information-theoretic performance limits. The method is inherently insensitive to general SPAM errors plaguing approaches like process tomography and can even circumvent the intrinsic SPAM errors (e.g. errors in the preparation of the ground state). Working at the Hamiltonian/Lindbladian level gives us greater control than what methods restricted to sets of pre-compiled gates provide, letting us use optimal control techniques when manipulating the system. This versatility permits us to use well known techniques like D-optimal experimental design to further improve the fidelity of our estimator.

There are many avenues that remain open for further exploration and extension. Because STEADY uses an inherently sparse description of the noise in the system, it is conceivable that the method could be made scalable. At the moment, this is precluded by the need to simulate the full quantum dynamics, but methods to speed this up could be employed such as using non-universal quantum circuits instead of random pulses, or using a well-calibrated quantum device together with STEADY to calibrate a new device. There is much room for additional theoretical work on the performance of STEADY, especially in regards to the family of models that we consider. Here we mainly considered a linear coupling to drives, but it should be possible to give a complete theoretical analysis of performance for more general classes of functions. It would also be interesting to include weak measurements in a stochastic master equation formalism and attempt to infer the measurement parameters as well. One additional open problem is to connect the performance of STEADY to other popular error metrics for quantum gates so as to facilitate an analysis of fault tolerant capabilities for quantum computation Kueng et al. (2016). Lastly, experimental implementation of these ideas will undoubtedly lead to further improvements in the method.

Acknowledgements.

This work would not have been possible without the contributions of the Python, Jupyter, Matplotlib, Numpy, Qutip, and Julia open source projects and the Yale HPC team. We acknowledge support from the ARL-CDQI (W911NF-15-2-0067, W911NF-18-2-0237), ARO (W911NF-18-1-0020, W911NF-18-1-0212), ARO MURI (W911NF-16-1-0349), AFOSR MURI (FA9550-14-1-0052, FA9550-15-1-0015), DOE (DE-SC0019406), NSF (EFMA-1640959), and the Packard Foundation (2013-39273). This work was also supported by the Australian Research Council through the Centre of Excellence in Engineered Quantum Systems CE170100009, and by the US Army Research Office grant numbers W911NF-14-1-0098 and W911NF-14-1-0103.

SK performed the majority of the software development and numerical tests. The design, evaluation, and repeated iterative improvements of the method were jointly done between all four authors.

Supplementary Materials for Stochastic Estimation of Dynamical Variables

The software we provide can be used both for running our method of stochastic Hamiltonian estimation and in reverse for performing gradient-based control (à la GRAPE). We admit both constant-in-time drive pulses and time-dependent drive pulses. For Hamiltonian dynamics we use an integrator that directly computes the evolution operator through diagonalization of the Hamiltonian. For Lindbladian dynamics we use an RK4 integrator. All integrators are implemented as fully differentiable operations, i.e. their gradients with respect to any parameters (control pulse parameters, parameters of the Hamiltonian or the Lindbladian) are computable analytically, to be used for parameter estimation (e.g. our stochastic Hamiltonian estimation method), experimental design (e.g. maximizing Fisher information), or control. Later in this appendix we describe the implementation in more details.

For most of the numerical tests of our method we use a simulated hardware with $Q=3$ qubits governed by a Hamiltonian containing Pauli drives ($\sigma_{x}$, $\sigma_{y}$, and $\sigma_{z}$) for each of the three qubits (marked in many of the figures as $X1$ through $Z3$) as well as nearest neighbor exchange interactions $\sigma_{i}^{+}\sigma_{j}^{-}+h.c.$ (marked in plots as pairs $12$, $23$, and $31$, depending on which pair of qubits they correspond to). For exact numerical values (which are on the order of 1, generated randomly), consult the referenced source code. We have $D=12$ drives for each of the $M=12$ components of the Hamiltonian. Moreover, there is no mixing between the drives (i.e. the true $\bm{\alpha}$ is diagonal). The $\sigma_{z}$ and the nearest neighbor contributions to the Hamiltonian are present even in the absence of drives, i.e. the corresponding $\bm{\beta}$ coefficients are non-zero; the rest of the $\bm{\beta}$ vector ($\sigma_{x}$ and $\sigma_{y}$) is zero.

The source code for this project is available at github.com/Krastanov/hamiltonian_estimation.

Appendix A The

Choice of Validation Function: How to Evaluate the Quality of Parameter Estimation

Our parameter estimation procedure is inherently stochastic, which is both a blessing and a curse. It is what permits us to use the entirety of the measurements performed on the hardware and reach the information-theoretic limits of precision. For numerical simplicity we choose a simple cost function—the mean squared error between measured and predicted populations—however we do observe consistent behavior independent of the particular choice of stochastic validation function (see Fig. 6a).

However, there might be a “gauge” degree of freedom in the way we have parameterized the Hamiltonian, which does not affect the actual dynamics of the system. As such, there might be parameters whose value will neither affect the measurement data we gather, nor will be of any consequence when designing control pulses down the road. Our choice of “indirect” validation function permits us to disregard this degree of freedom, as the validation function is sensitive only to the prediction of our model, not to the particular parameterization we have used in our model.

However, it would be instructive to observe how some of these “unimportant” parameters behave in our estimator. For instance, the test system used in much of this manuscript has a Hamiltonian of the form

[TABLE]

where $\kappa_{i}$ denote the coupling strengths for each drive pulse component, and $\varepsilon_{i}$ and $\eta_{i}$ are the strengths of Hamiltonian components that are present even in the absence of drives. Our estimator tries to model this Hamiltonian as $\tilde{H}(\bm{\alpha},\bm{\beta};\bm{d})=\stackrel{{\scriptstyle[}}{{k}}=1]{M}{\sum}a_{k}A_{k}$, where $a_{k}=\stackrel{{\scriptstyle[}}{{l}}=1]{M}{\sum}\alpha_{kl}d_{l}+\beta_{k}$ and the set of possible $A_{k}$ is $\{\sigma_{z1},\sigma_{z2},\sigma_{z3},\sigma_{x1},\sigma_{x2},\sigma_{x3},\sigma_{y1},\sigma_{y2},\sigma_{y3},\sigma_{1}^{+}\sigma_{2}^{-}+h.c.,\sigma_{2}^{+}\sigma_{3}^{-}+h.c.,\sigma_{3}^{+}\sigma_{1}^{-}+h.c.,\}$. We can see that the “true” $\bm{\beta}$ is zero for the subset of operators $\{\sigma_{x1},\sigma_{x2},\sigma_{x3},\sigma_{y1},\sigma_{y2},\sigma_{y3}\}$, and the “true” $\bm{\alpha}$ is a diagonal matrix. However, in the measurement data used by the estimator, there is nothing that defines the orientation of the x axis in the xy plane. If we rotate all the qubits by the same angle around the z axis, the physics of the system (and the measured data) will not change. In other words, for this particular Hamiltonian, the change of xy basis generator $\sigma_{z1}\otimes\sigma_{z2}\otimes\sigma_{z3}$ commutes with the zero drive Hamiltonian $H(\bm{0})$, hence non-diagonal $\bm{\alpha}$ such as the one in Fig. 6c would predict the same evolution as the diagonal one.

In the case of a microwave superconducting circuit implementation, this freedom in the parameter would be due the fact that the master oscillator used in mixing the control signals has an arbitrary initial phase (this phase defines our choice of x and y axis in the rotating frame for all qubits).

More generally, as long as the random pulses we use in our stochastic cost function are representative of the pulses that will be used to control the hardware, such “gauge” degrees of freedom are inconsequential: if they do not affect our cost function, they will not affect the result of a computation running on the hardware either.

Appendix B The Choice of Cost Function, Estimator Efficiency, and the Cramér–Rao

Bound

In the main text we invoked the Cramér–Rao bound in order to explain the $\frac{1}{P\times S}$ fidelity scaling that we obtain in the case of unbiased estimation (before considering the effects of SPAM). That bound is a general information inequality that expresses a lower bound on the variance. Besides verifying numerically that we follow the bound (as done in the main text), it would be instructive to derive the variance of our particular estimator and compare it to this bound. This would help inform and defend our choice of $\text{dist}(\hat{\bm{p}},\tilde{\bm{p}})$, which until now has been based on computational convenience.

Let $\bm{G}(\bm{\omega})$ be the gradient of $C(\bm{\omega})$, and $\bm{H}(\bm{\omega})$ be the Hessian matrix. For our estimator we have $\bm{G}(\bm{\hat{\omega}})=0$ where $\bm{\hat{\omega}}$ is our estimate for the true value $\bm{\omega_{0}}$. Expanding $\bm{G}(\bm{\omega})$ around $\bm{\omega_{0}}$ we get

[TABLE]

which leads us to

[TABLE]

This permits us to estimate variances and covariances for the parameters

[TABLE]

If we pick a maximal likelihood estimator (i.e. if the distance function is the cross entropy $\text{dist}(\hat{\bm{p}},\tilde{\bm{p}})=\tilde{\bm{p}}.\log(\hat{\bm{p}})$), then we have $A=B=\bm{\mathcal{I}_{F}}$ and $\bm{\mathcal{I}_{C}}=\bm{\mathcal{I}_{F}}$, where $\bm{\mathcal{I}_{F}}$ is the Fisher information matrix, therefore proving that we would saturate the Cramér–Rao bound and have a fully efficient estimator. In practice, one would need to consider the numerical stability of gradient descent as well: the estimator would be run with a number of different distance functions, to see empirically which one proves most reliable numerically. This led us to use the mean squared error in our examples, but other choices might be more performant in other settings (Fig. 7).

For completeness we also derive the explicit form of the Fisher information matrix used in arguments in the main text. First we denote the log-likelihood function

[TABLE]

where $\tilde{p}_{ik}(\bm{\omega})$ is the predicted population in the $k$-th state for the $i$-th pulse in the training set $\{\bm{d}_{1},\dots,\bm{d}_{P}\}$ and $\hat{p}_{ik}$ is the measured population for that state and pulse (by taking $S$ samples). This leads us to an expression for the Fisher information

[TABLE]

where $\langle\dots\rangle$ denotes average over all pulses in the training set. Also, $\sum_{k=1}^{2^{Q}}\frac{1}{\tilde{p}_{k}}\frac{\partial\tilde{p}_{k}}{\partial\omega_{i}}\frac{\partial\tilde{p}_{k}}{\partial\omega_{j}}$ happens to be the Fisher information for a single measurement of a single pulse, i.e. the Fisher information is additive.

Appendix C Convergence, Overfitting, and Model Errors

In the previous discussion we neglected issues of convergence, overfitting, and model errors. We use two complementary tools to fight these problems. Firstly, overfitting or being stuck in a valley of the cost function can both be avoided with an annealed regularization cost applied to the parameters.

The specific annealed and regularized cost function that we use as a function of iteration $k$ is of the form

[TABLE]

where $C(\bm{\omega})$ is given by Eq. 5 and the $\lambda_{k}\geq 0$ are annealed as $k$ increases, meaning that our estimator is essentially a LASSO-type estimator Tibshirani (1996).

From a purely practical point of view this avoids the initial steps of the gradient descent going in a wildly unphysical direction of pathological values for the parameters. More importantly, minimizing the 1-norm of the parameters provides for a sparser realistic parametrization. For particularly difficult programs we observed that empirical tricks akin to MacKay’s regularization annealing schedule (MacKay, 1992), where the variance in the regularization is forced to follow the variance in the cost enables reliable convergence.

The second tool involves extending the dynamics permitted in the model. Throughout the main text and in the following sections we discuss the biasing effect of various types of model errors (e.g. intrinsic SPAM not being included in the simplest models, or the effect of neglected parasitic couplings or other coherent errors in Fig. 8, or incoherent noise that can not be expressed in a Hamiltonian formalism). A quick fix solution we discussed is making the estimator more sensitive to parameter errors than to model errors by, e.g., using stronger pulses that accentuate the parameter errors. However, the more powerful solution, as long as it is computationally feasible, is to extend the model to include these otherwise neglected dynamics. Both of these approaches are described in the following sections.

Lastly, it is important to note the interaction between the regularization approach and the model extension approach. Permitting many more degrees of freedom in the model also leads to a higher risk of overfitting and convergence issues, hence leading to the need for annealed regularization. For particularly difficult problems one can even envision extending STEADY to a “hyper parameter” estimator where a discrete optimization algorithm evaluates a family of models by the Akaike information criterion (Akaike, 1974).

Appendix D Intrinsic SPAM

As described in the main text, we used the following model for the intrinsic state preparation error

[TABLE]

and the following model for the measurement error: $\bm{p}_{\text{meas}}=\mathcal{S}\boldsymbol{p}$, where

[TABLE]

and the only non-zero components of $\mathcal{S^{\prime}}$ are $\mathcal{S^{\prime}}_{ij}=s$ where the Hamming distance between the binary representations of $i$ and $j$ is 1.

To find the bias $\bm{b}$ in our estimator $\hat{\bm{\omega}}$ we will need to take a more careful look at the minimum of $C_{s}(\bm{\omega})$ ($C_{s}(\bm{\omega})$ is what $C(\bm{\omega})$ becomes in the presence of SPAM errors as described below). By definition of $\bm{b}$, that minimum will be at $\bm{\omega_{0}}+\bm{b}$. For simplicity let us consider only measurement errors, which would cause the vector $\hat{\bm{p}}_{i}$ to become $\hat{\bm{\pi}}_{i}=\mathcal{S}\hat{\bm{p}}_{i}$, where $\mathcal{S}$ is an almost-diagonal stochastic matrix, with off-diagonal components on the order of $\mathcal{O}\left(s\right)$, i.e. $\mathcal{S}=I+s\mathcal{S}^{\prime}+\mathcal{O}(s^{2})$, where $\mathcal{S}^{\prime}$ contains constant components on the order of unity. $\mathcal{S}$ describes the chance that a measurement of a given state is wrongly reported as another state. While the state preparation errors are more complicated to express, because they happen before the dynamical evolution of the state, linearity guarantees that a similar treatment would work for an appropriately chosen ensemble of states. This leaves us with

[TABLE]

Evaluated at $\bm{\omega_{0}}+\bm{b}$ it results in

[TABLE]

where $\delta$ is a positive real number. Moreover, we have $C_{s}(\bm{\omega_{0}}+\bm{b})-C(\bm{\omega_{0}})\sim\alpha s^{2}+\mathcal{O}(s^{3})$ (because $s\mapsto\underset{\bm{\omega}}{\min}\left(C_{s}(\bm{\omega})\right)$ has its minimum at $s=0$) and $C(\bm{\omega_{0}}+\bm{b})-C(\bm{\omega_{0}})\sim\beta|\bm{b}|^{2}+\mathcal{O}(|\bm{b}|^{3})$ (because $\bm{b}\mapsto\underset{\bm{\omega}}{\min}\left(C(\bm{\omega}+\bm{b})\right)$ has its minimum at $\bm{b}=0$), where $\alpha$ and $\beta$ are positive real numbers. Therefore

[TABLE]

Given that we are interested in the regime where the bias overwhelms the statistical error, we can take the limit $S\to\infty$ which results to leading order in $\left(\alpha-\delta\right)s^{2}-\beta|\bm{b}|^{2}-\gamma s|\bm{b}|=0$, where $\gamma$ is a positive real number. This leads to $|\bm{b}|\propto s$ and the observed error floor of $V\propto|\bm{b}|^{2}\propto s^{2}$ (see Fig. 3d).

This error floor is unsurmountable by simply increasing the number of acquired measurements. However, if we redo the expansion while keeping record of $T$, we are left with $|\bm{b}|\propto\frac{s}{T}$. This immediately suggests a way to decrease the bias of our estimator: simply use longer (or more powerful) control pulses. Fig. 3c demonstrates the improvements due to this approach.

Appendix E Effects of

Non-unitary Decay on Information Content

As we have seen in the rest of the manuscript, having an incomplete model, one that is incapable of expressing the entire dynamics, would cause bias in our estimator. This was observed both in the case of intrinsic SPAM, and in the case of Lindbladian dynamics. This issue can be addressed in some cases by making the estimator more sensitive to estimator errors (by using longer pulses). When this fails one can instead extend the model to include the missing dynamics. Both of these approaches were discussed in the main text.

However, the good performance of the Lindbladian model estimator from Fig. 5 can be counterintuitive. Taken to the extreme, strong decay would cause all of the information about the unitary evolution to leak out to the environment before the measurement. This extreme example seems to imply that the performance of the estimator should drop at extremely high decay parameters, which we do not observe. Similarly to the discussion in A, this stems from our choice of validation function; that is, we evaluate the quality of the predictions we can make about our system, and we do not evaluate directly how precise each parameter is estimated. Hence, when we are interested in the quality of predictions, we do not need to worry about parameters that do not affect the dynamics strongly. In the case of strong decays, all other parameters become unimportant, and this is why we do not see a drop in performance. The same phenomenon that makes parameters hard to estimate also makes them inconsequential to the dynamics of the system.

This can be observed in the Fisher information. Consider for simplicity a two level system. The excited state Born probability is $p(\omega)$ in the case of unitary evolution. In the presence of decay, a first order approximation for that same probability is $p(\omega)\mathrm{e}^{-\Gamma T}$. The Fisher information with respect to $\omega$ for a single measurement is then $\bm{{\mathcal{I}_{F}}}_{\omega}=\frac{\mathrm{e}^{-\Gamma T}}{p}p^{\prime 2}(\omega)+\frac{\mathrm{e}^{-2\Gamma T}}{1-p\mathrm{e}^{-\Gamma T}}p^{\prime 2}(\omega)=p^{\prime 2}(\omega)\frac{\mathrm{e}^{-\Gamma T}}{p\left(1-p\mathrm{e}^{-\Gamma T}\right)}$, which in the limiting cases is $\bm{{\mathcal{I}_{F}}}_{\omega}\underset{\Gamma\rightarrow\infty}{\sim}\frac{p^{\prime 2}(\omega)\mathrm{e}^{-\Gamma T}}{p(\omega)}$ and $\bm{{\mathcal{I}_{F}}}_{\omega}\underset{\Gamma\rightarrow 0}{\sim}\frac{p^{\prime 2}(\omega)}{p(\omega)}+\mathcal{O}(\Gamma T)$. Similarly, the Fisher information with respect to $\Gamma$ is $\bm{{\mathcal{I}_{F}}}_{\Gamma}=p^{2}T^{2}\mathrm{e}^{-2\Gamma T}\left(\frac{1}{p\mathrm{e}^{-\Gamma T}}+\frac{1}{1-p\mathrm{e}^{-\Gamma T}}\right)$, with limiting behavior $\bm{{\mathcal{I}_{F}}}_{\Gamma}\underset{\Gamma\rightarrow\infty}{\sim}p^{2}T^{2}\mathrm{e}^{-\Gamma T}$ and $\bm{{\mathcal{I}_{F}}}_{\Gamma}\underset{\Gamma\rightarrow 0}{\sim}\frac{p}{1-p}T^{2}+\mathcal{O}(\Gamma T)$.

However, the validation function is to first order $V\propto\left(p(\omega_{0})\mathrm{e}^{-\Gamma_{0}T}-p(\omega)\mathrm{e}^{-\Gamma T}\right)^{2}$, where $\omega_{0}$ and $\Gamma_{0}$ are the true values of the parameters. Hence $V\propto\mathrm{e}^{-2\Gamma_{0}T}\left(p^{\prime}(\omega_{0})(\omega_{0}-\omega)-Tp(\omega_{0})(\Gamma_{0}-\Gamma)\right)^{2}\propto\mathrm{e}^{-2\Gamma_{0}T}\left(p^{\prime 2}(\omega_{0})\sigma_{\omega}^{2}+T^{2}p^{2}(\omega_{0})\sigma_{\Gamma}^{2}\right)$, where $\sigma_{\omega}$ and $\sigma_{\Gamma}$ are the variances in the estimators for each of the parameters (we neglect correlations). Due to the Cramér–Rao bound both variances scale as the inverse of the corresponding Fisher information, which leaves us with $V\underset{\Gamma\rightarrow\infty}{\sim}p(\omega_{0})\mathrm{e}^{-\Gamma T}$ (see Fig. 9). Therefore, as the decay rate goes higher and leaves us with less and less available information per measurement, it also causes less diversity in the final measurements (high probability that the final state is the ground state), leading to small (“good”) values for the validation function. In practice, this effect becomes noticeable only at impractically high values for $\Gamma$: while we do precisely predict in such cases how the state decays, this is not of particular use for devising control protocols for the system. Generally a system with strong decays is not particularly useful as a quantum hardware, unless we use “engineered dissipation” control schemes or control schemes employing virtual states like STIRAP (Unanyan et al., 1998), which are beyond the scope of this manuscript.

Appendix F Variance of Parameter Estimator in Linear Least Squares

Here we give a concrete analytical example of cases where the cost function has a floor dictated by $P$, while the validation function has a floor dictated by $S\times P$. The description here is generic; it does not refer to a model of quantum dynamics.

The model we are learning is $f(\omega;x)$ parameterized by $\omega$. We denote by $(y_{i},x_{i})$ the pairs of (population estimate, pulse) that we are learning from (we have $P$ such pairs). For each pair we used $S$ repetitions in order to estimate $y_{i}$, the population, for the given $x_{i}$, the pulse.

The true population is $y_{i}^{\circ}=y_{i}-e_{i}$ (given that we are performing multinomial sampling, we approximate $e_{i}$ as drawn from a normal distribution with $\sigma_{0}=\frac{\sqrt{p}}{\sqrt{S}}$).

We perform parameter estimation by minimizing least squares $C(\omega)=\frac{1}{P}\stackrel{{\scriptstyle[}}{{1}}]{P}{\sum}\left(y_{i}-f(\omega,x_{i})\right)^{2}=\frac{1}{P}\stackrel{{\scriptstyle[}}{{1}}]{P}{\sum}\left(y_{i}^{\circ}+e_{i}-f(\omega,x_{i})\right)^{2}$ and the validation function is $V(\omega)=\frac{1}{P}\stackrel{{\scriptstyle[}}{{1}}]{P}{\sum}\left(y_{i}^{\circ}-f(\omega,x_{i})\right)^{2}$.

For the purpose of this example, consider linear least squares: the model we are fitting to has the parameters $\omega=a,b$, and $f(a,b;x)=a+bx$. $a$ and $b$ denote the parameter values which minimize $C$. However, the data are governed by the “true” model $y^{\circ}=f(a^{\circ},b^{\circ};x)=a^{\circ}+b^{\circ}x$.

The minimum of $C(\omega)$ is reached at $b=\frac{\hat{\sigma}_{xy}}{\hat{\sigma}_{xx}}$ and $a=\bar{y}-b\bar{x}$, where:

• •

a bar denotes the typical estimator of an average value;

• •

$a=\bar{y}-b\bar{x}=\bar{y}^{\circ}+\bar{e}-b\bar{x}$;

• •

$\hat{\sigma}_{xx}$ denotes the estimator of the variance of $x$;

• •

$\hat{\sigma}_{xy}$ is the estimator of the covariance of $x$ and $y$.

The estimated covariance can be calculated as,

[TABLE]

Substituting the expressions for $y_{i}^{\circ}$ and $a$ in $V$ and factoring out the $(x_{i}-\bar{x})$ term leaves us with:

[TABLE]

Hence the expectation value of $V$ is $\frac{p}{PS}$.

Appendix G Software Implementation

The estimators we implement are all fully differentiable programs, that can run “batches” of measurement data. They can run on both CPUs and GPUs, as well as any other optimized tensor processing units supported by Tensorflow, like Google’s TPU chips. The code is documented inline and extremely short (<50 lines of code per estimator). We implement the following estimators (all of them can be run either as stochastic parameter estimators where measurements are provided as training data, or as optimal control optimizers where the parameters are fixed but the control pulses are optimized for a given target).

• •

StateProbabilitiesPaulied solves Schroedinger’s equation for

[TABLE]

and where $\{A_{1},\dots,A_{M}\}$ are fixed in advance. Control pulses for this solver are constant in time (with implied duration of $T=1$) provided as arrays of shape $P\times D$. The time-evolution operator is calculated through diagonalization of the Hamiltonian.

• •

StateProbabilities solves Schroedinger’s equation for

[TABLE]

with explicitly ensured Hermiticity (see the source code for details). Control pulses for this solver are constant in time (with implied duration of $T=1$) provided as arrays of shape $P\times D$. The time-evolution operator is calculated through diagonalization of the Hamiltonian.

• •

Non-linear drives can be straightforwardly added to the above solvers by adding higher-order terms like $\sigma_{ijkl}d_{k}d_{l}$, etc.

• •

StateProbabilitiesTimeDep solves Schroedinger’s equation in either of the previous two forms, but for time dependent control pulses $d_{k}$. Control pulses are piece-wise constant of implied total duration $T=1$ and a predetermined number of timesteps $\Theta$. The control pulses are provided as arrays of shape $P\times\Theta\times D$. The time evolution is computed through evaluating the $\Theta$ consecutive time-evolution operators. Each operator is computed either, as above, through diagonalization, or through faster, but less precise, Taylor expansion of the matrix exponential (in Horner’s form).

• •

StateProbabilitiesTimeDepLindblad solves Lindblad’s Master equation for time dependent control pulses $d_{k}$. The Hamiltonian part is provided in either of the two forms discussed above. The control pulses are provided in the same format as above. The non-unitary evolution is modeled by fixed predetermined collapse operators $\{L_{1},\dots,L_{C}\}$, however the strengths of each collapse operator $\{c_{1},\dots,c_{C}\}$ is a parameter optimized by the estimator. The full equation being modeled is

[TABLE]

The integrator uses Euler’s method.

• •

StateProbabilitiesTimeDepLindbladRK4 works as above but it uses the RK4 method for numerical integration.

• •

Unlike gradients, Hessians do not permit fast backpropagation methods for their calculation (and hence are not well supported by differentiable programming frameworks yet). This led us to writing a simple numerical procedure for estimating the Hessians needed for the calculation of the Fisher information in AvgFIStateProbabilitiesPaulied. It is a proof-of-concept that works only on the first of the estimators described above, but it is the tool that permitted us to run experimental-design optimizations as described in the main text. It is implemented in a separate Jupyter notebook.