Quantum process identification: a method for characterizing non-markovian quantum dynamics

The problem we address is to how to experimentally characterize the behavior of a quantum system of interest under some dynamical process, without assuming anything about the system's size or composition or the nature of the process. The process in question may be a repeatable procedure ${ \mathcal P }$ with definitive start and end, or it may be a continuous, ongoing dynamic. In the latter case we let ${ \mathcal P }$ denote evolution for some fixed time Δτ, where Δτ becomes the time resolution at which the dynamics are resolved.

Ultimately, all that can be known about a system comes from one's interactions with it. The kind of characterization we develop is relative to, and in terms of, the available ways of interacting with the system. In this work we consider only interactions in the form of 'experiments' in which the system is prepared, made to evolve under ${ \mathcal P }$, and subsequently measured. (In section 6 we will speculate how QPI may be extended to more complicated interactions.) Let ${ \mathcal I }$ be the set of ways the system can be prepared (initialized) and let ${ \mathcal M }$ be the set of ways it can be measured. While in principle there exists a continuum of possible initializations and possible measurements, in practice one can utilize only finitely many of these; thus we assume ${ \mathcal I }$ and ${ \mathcal M }$ are finite sets¹ . An experiment is a specified by a triple $(i,t,m)$ which specifies the initialization $i\in { \mathcal I }$, followed by t repetitions of ${ \mathcal P }$, followed by the measurement $m\in { \mathcal M }$. For notational convenience we suppose that each measurement has only two possible outcomes, YES or NO. (A measurement with k possible outcomes can be cast as k distinct YES/NO measurements.) Under these premises, the observable behavior of the system is the set of YES outcome probabilities for all possible experiments. The YES probability for experiment (i, t, m) will be denoted by ${F}_{i,m}^{(t)}$.

It is worth reiterating that the characterization of ${ \mathcal P }$ produced by QPI is completely relative to the set of initializations ${ \mathcal I }$ and set of measurements ${ \mathcal M }$ considered. ${ \mathcal I }$ and ${ \mathcal M }$ need not be tomographically complete for the system, but if they are not then the resulting characterization may not cover the full state space of the system. Furthermore, if an 'absolute' characterization is desired, ${ \mathcal I }$ and ${ \mathcal M }$ should include the procedures that define the reference frame. As a side result, QPI also yields a characterization of the initializations and measurements relative to one another.

For our purposes, characterization amounts to constructing a mathematical model of the process ${ \mathcal P }$, initializations ${ \mathcal I }$, and measurements ${ \mathcal M }$. We opt not to use the traditional Hilbert-space formulation of quantum mechanics, but instead the more general framework known as generalized probabilistic theory (GPT) [59–61]. GPT is a framework for a class of physical theories that includes quantum mechanics and classical mechanics as special cases. For us GPT has two appealing features: first, physical phenomena are described in operational terms. Rather than expressing phenomena in terms of complex operators that are not directly observable, GPT expresses phenomena in terms of the outcome probabilities of measurements that could be performed. This matches the task of experimental characterization extremely well. Secondly, the GPT formalism does not distinguish between quantum and classical behaviors; it treats all behaviors in a uniform way. This is useful because although the system of interest and the environment are nominally quantum, their behavior may have some purely classical aspects. GPT allows one to construct the most economical model of a process without choosing it to be quantum or classical. On the other hand, some symmetries required by quantum mechanics (such as complete positivity) are less conveniently expressed in GPT. We emphasize that QPI is not tied to the GPT formalism; if desired, QPI could be formulated entirely in terms of the traditional Hilbert space formalism.

One approach to representing non-Markovian dynamics is to construct dynamical maps that are not linear and/or not CPTP [36]. Another strategy is to go beyond time-local maps and use structures that explicitly encode temporal correlations, such as spectral densities [62, 63], autoregression parameters [19], or process tensors [20, 64]. QPI takes a different approach: the system's temporal evolution (whether Markovian not) is modeled by a fixed Markovian process involving both the system of interest and a sufficiently large abstract environment [65]. (This is analogous to Stinespring dilation of a CPTP channel [66, 67].) That is, an extended Markovian model is formulated to predict the system's observable (possibly non-Markovian) behavior. In technical terms, the central premise of QPI is that the system's observable behavior can be modeled to desired accuracy by a sufficiently large finite-dimensional, time-independent, time-local dynamical model² . Significantly, such a model can be inferred from measurements that nominally involve the system alone. This is important because the experimentalist often does not have access to, or even knowledge of, all the pertinent parts of the environment. By inferring a model on an extended state space, QPI implicitly characterizes the involvement of pertinent environmental degrees of freedom (though these degrees of freedom are abstract and might not be readily associated with particular physical degrees of freedom).

The only other assumption of QPI is that each experiment is statistically independent. One way of satisfying this second assumption is to use strong initializations that erase (as far as can be detected) any correlations between the system and the environment due to previous experiments. For example, the initialization procedure may include waiting a time much longer than the plausible memory lifetime of the environment. Note, however, that we allow the initialization procedure itself to correlate the system and environment. A weaker way of satisfying the independence assumption is to randomize the order of the experiments so that any residual correlations between the system and environment at the start of each experiment are sampled fairly. In this case, QPI infers a distribution of correlated initial states.

In QPI, a model of dimension d consists of a state vector ${s}_{i}\in {{\mathbb{R}}}^{d}$ for each initialization procedure $i\in { \mathcal I }$, a property vector³ ${p}_{m}\in {{\mathbb{R}}}^{d}$ for each measurement procedure $m\in { \mathcal M }$, and a transfer matrix $T\in {{\mathbb{R}}}^{d\times d}$ for the process ${ \mathcal P }$. (The relation between these model elements and conventional quantities is shown in table 1.) Then the YES probability of experiment (i, t, m) is

$$\begin{eqnarray}&&{F}_{i,m}^{(t)}={s}_{i}{T}^{t}{p}_{m}.\end{eqnarray} \tag{ 1 }$$

(We have adopted the convention in which states are row vectors and operators are applied on the right.) In terms of the matrices $S=[\begin{array}{ccc}{s}_{1}; & {s}_{2}; & \ldots \end{array}]$ and $P=[\begin{array}{ccc}{p}_{1} & {p}_{2} & \ldots \end{array}]$, we have

$$\begin{eqnarray}&&{F}^{(t)}={{ST}}^{t}P.\end{eqnarray} \tag{ 2 }$$

The goal of QPI is to find matrices S, P, and T that accurately predict F^(t) for all t = 0, 1, 2, ..., i.e. the observable behavior of the system as defined above. Here S and P are a characterization of the state preparations and measurements, respectively, while T characterizes the process ${ \mathcal P }$ itself. We note that these matrices can be determined only up to a similarity transformation, since for any invertible matrix G the mapping $S\to {SG}$, $P\to {G}^{-1}P$, $T\to {G}^{-1}{TG}$ leaves all probabilities unchanged. This indeterminacy of representation is known as 'gauge indeterminacy'. While gauge indeterminacy does lead to a minor theoretical conundrum regarding the assessment of process fidelity [68], it has (by definition) no experimental consequence and thus is not a limitation of QPI.

Table 1.  Elements of a QPI model and the corresponding objects in the Hilbert-space formulation of quantum mechanics.

QPI model element	Corresponding Hilbert-space object
State row vector s	Density operator ρ
Property vector p	Observable Hermitian operator A
Transfer matrix T	CPTP map ${ \mathcal E }$
Expectation value sTp	Expectation value $\mathrm{Tr}(A{ \mathcal E }(\rho ))$

In the case of a continuous-time process, QPI can be used to obtain an effective master equation. Let $\rho (\tau )\in {{\mathbb{R}}}^{d}$ denote the extended system-environment state at time τ with F^(τ) ≡ ρ(τ)P. From (2) we have $\rho (t{\rm{\Delta }}\tau )={{ST}}^{t}=\rho ((t-1){\rm{\Delta }}\tau )T$ where Δτ is the chosen time resolution. For sufficiently small Δτ, T ≈ 1 + LΔτ for some matrix L. Then ρ approximately obeys the master equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\rho (\tau )}{{\rm{d}}\tau }=\rho (\tau )L.\end{eqnarray} \tag{ 3 }$$

A key to solving the problem of characterizing non-Markovian quantum processes is a technique that was devised half a century ago in the context of discrete-time linear systems and is now a standard topic in introductory texts on systems engineering [21, 22]. Consider a classical input-output system whose output $x\in {{\mathbb{R}}}^{n}$ at each time $t\in {\mathbb{Z}}$ is a linear function of its input $u\in {{\mathbb{R}}}^{m}$ at the previous d times. Ho and Kalman [69] showed that the system can be described by a d-dimensional state model of the form

$$\begin{eqnarray}\begin{array}{rcl}s(t+1) & = & u(t)A+s(t)B\\ x(t) & = & s(t)C\end{array}\end{eqnarray} \tag{ 4 }$$

where $s\in {{\mathbb{R}}}^{d}$, $A\in {{\mathbb{R}}}^{m\times d}$, $B\in {{\mathbb{R}}}^{d\times d}$, and $C\in {{\mathbb{R}}}^{d\times n}$. Furthermore, they showed how to obtain the coefficients A, B, C from the observable quantities X(1), ..., X(2d) where X_i,j(t) is the value of the output x_j(t) in response to a unit impulse on the ith input (u_i) at time 0. In this case $X(t+1)={{AB}}^{t}C$. Let

$$\begin{eqnarray}H=\left[\begin{array}{ccc}X(1) & \cdots & X(d)\\ \vdots & & \vdots \\ X(d) & \cdots & X(2d-1)\end{array}\right]\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}{H}^{{\prime} }=\left[\begin{array}{ccc}X(2) & \cdots & X(d+1)\\ \vdots & & \vdots \\ X(d+1) & \cdots & X(2d)\end{array}\right].\end{eqnarray} \tag{ 6 }$$

The key fact here is that H and ${H}^{{\prime} }$ have related rank-d factorizations, H = LR and ${H}^{{\prime} }={LBR}$ where $L=[A;{AB};\,\cdots ;\,{{AB}}^{d-1}]$ and $R=[{CBC}\cdots {B}^{d-1}C]$. The Ho–Kalman result is that a model of the form (4) can be obtained by finding any rank-d factorization H = LR, taking the first m rows of L for A, the first n columns of R for C, and taking $B={L}^{+}{H}^{{\prime} }{R}^{+}$ where ${}^{+}$ denotes the Moore–Penrose pseudo-inverse.

The relevance of this result to QPI becomes clear upon noting that $X(t+1)={{AB}}^{t}C$ has the exactly same form as equation (2). From the perspective of process tomography, A describes the preparable states and C describes the measurable properties. When $d\gt \min (m,n)$ these do not span the state space of the system, yet a full reconstruction of the process is still possible. That is, observations of the system over a sufficiently long period of time can yield complete information about the system, even when the preparable states and measurable properties are not intrinsically complete [27]. The critical implication of this result for QPI is that the properties of a system of interest over time are sufficient to reconstruct an accurate model of all relevant dynamics, including the dynamics of any external degrees of freedom with which the system interacts.

Mathematically, the construction of a d-dimensional model from H is possible because H has rank d. What if A and C are already tomographically complete? In that case it ought to be sufficient to use just X(1) = AC for H and X(2) = ABC for ${H}^{{\prime} }$, as in QPT. In this case, the textbook Ho–Kalman method asks for more data than is actually necessary. The following Lemma, proved in appendix A, states that the number of time steps needed to ensure a complete characterization is determined by the number of latent (unmeasured) degrees of freedom only, not the total number of degrees of freedom:

Lemma 1. Let $X(t)$ be the impulse response function for an arbitrary linear, discrete-time, time-shift invariant system of dimension d. Let $r={\max }_{t}\mathrm{rank}\ X(t)$ denote the number of directly accessible dimensions and let $l\equiv d-r$ denote the number of latent dimensions. Let

$$\begin{eqnarray}{H}^{(t;k)}\equiv \left[\begin{array}{cccc}X(t) & X(t+1) & \cdots & X(t+k)\\ X(t+1) & & & \vdots \\ \vdots & & & X(t+2k-1)\\ X(t+k) & \cdots & X(t+2k-1) & X(t+2k)\end{array}\right].\end{eqnarray} \tag{ 7 }$$

For $k\geqslant l$, $\mathrm{rank}\ {H}^{(t;k)}=d$, i.e. ${H}^{(t;k)}$ captures the full dimensionality of the system. Furthermore a minimal ($d$-dimensional) model of the system can be constructed from ${X}^{(1)},\,\ldots ,\,{X}^{(2l+2)}$. In particular, if ${LR}$ is any rank factorization of ${H}^{(1;l)}$, then $X(t+1)={{AB}}^{t}C$ where $A$ is the first $m$ rows of $L,C$ is the first $n$ columns of $R$, and $B={L}^{+}{H}^{(2;l)}{R}^{+}$. Some $d$ dimensional systems can be reconstructed from fewer observations and/or observations at non-consecutive times, but some cannot. Thus $k\geqslant l$ is necessary to guarantee complete characterization of an unknown $d$ dimensional system.

The Ho–Kalman method is the basis of our approach to QPI. However, there are two important differences between a (classical) linear input-output system and a quantum process that require consideration:

• In the Ho–Kalman problem X_i,j(t) is (in principle) a directly observable quantity, whereas in QPI the measurable quantity ${F}_{i,m}^{(t)}$ is a probability. One consequence is that in QPI, many repetitions of an experiment are needed to obtain each value ${F}_{i,m}^{(t)}$. Even then, one only obtains an estimate ${\tilde{F}}_{i,m}^{(t)}$ whose precision is limited by the number of repetitions of the experiment. There is a large body of work on the extension of the Ho–Kalman approach to noisy data, but it considers additive or multiplicative Gaussian noise and is not applicable to QPI, where ${\tilde{F}}_{i,m}^{(t)}$ has a binomial distribution whose width depends on the (unknown) value of ${F}_{i,m}^{(t)}$.
• The observables in the Ho–Kalman method are all mutually compatible. In principle, the values ${X}_{i,j}(t)$ for all j and all t ≥ 1 can be obtained in a single experiment with initial state i. In contrast, measuring a quantum system generally disturbs it (an aspect of the uncertainty principle). Thus in QPI, not all measurable properties of the system can be measured at the same time, nor is the evolution of the system after a measurement the same as if the system had not been measured. Each ${F}_{i,m}^{(t)}$ must be estimated from a distinct set of experiments, in which the system evolves unobserved for t process repetitions and is then measured.

In the next section we present a method for QPI based on an extension of the Ho–Kalman approach to the quantum domain.

QPI is accomplished by performing a series of experiments, each repeated many times (figure 2). Each experiment is designated by a triple (i, t, m) which designates the ith way of initializing the system, followed by t repetitions of the process to be characterized, followed by the mth way of measuring the system. For each experiment (i, t, m) one records ${N}_{i,m}^{(t)}$, the number of times the experiment was performed, and ${Y}_{i,m}^{(t)}$, the number of times the YES outcome was observed. The frequency ${\tilde{F}}_{i,m}^{(t)}\equiv {Y}_{i,m}^{(t)}/{N}_{i,m}^{(t)}$ is the experimental estimate of ${F}_{i,m}^{(t)}$. (Throughout this work, an experimental estimate of a quantity Q is denoted by $\tilde{Q}$.)

Zoom In Zoom Out Reset image size

The set of experiments performed must be chosen with some care in order for QPI to be effective. If the experimental probabilities could be determined perfectly, then by lemma 1 ${F}^{(0)},\,\ldots ,\,{F}^{(2l+1)}$ would be sufficient to characterize any process of dimension $d\leqslant l+\mathrm{rank}\ {F}^{(0)}$. (Note that $\mathrm{rank}\ {F}^{(0)}$ is just the dimension of the space spanned by the initial states and measurements. If these are tomographically complete for the system of interest, $\mathrm{rank}\ {F}^{(0)}$ is just the system dimension.) However, some of the degrees of freedom may be slow to manifest in the sense that it takes many repetitions of the process for their impact on observable properties to become significant. In such cases the matrix ${H}^{(0;l)}$ is ill-conditioned and even a small amount of statistical error in the data yields an extremely poor characterization. For this reason it is advantageous to perform experiments covering a wide range of t values. As a general rule, an experiment with $2t$ process repetitions is twice as sensitive to the process parameters T as an experiment with t process repetitions [7].

The other main consideration is that the Ho–Kalman reconstruction technique requires the matrix H to have a generalized Hankel block structure:

$$\begin{eqnarray}H=\left[\begin{array}{cccc}{F}^{({i}_{1}+{j}_{1})} & {F}^{({i}_{1}+{j}_{2})} & \cdots & {F}^{({i}_{1}+{j}_{n})}\\ {F}^{({i}_{2}+{j}_{1})} & {F}^{({i}_{2}+{j}_{2})} & & {F}^{({i}_{2}+{j}_{n})}\\ \vdots & & & \vdots \\ {F}^{({i}_{m}+{j}_{1})} & {F}^{({i}_{m}+{j}_{2})} & \cdots & {F}^{({i}_{m}+{j}_{n})}\end{array}\right]\end{eqnarray} \tag{ 8 }$$

for some choice of integers i₁, ..., i_m and j₁, ..., j_n.

Our approach to satisfying these criteria is to use multiple 'flights' of experiments, where a flight is a set of experiments {(i, t, m)} with $2\tilde{l}+1$ consecutive t values and all values of i, m for each t (figure 3(a)). By lemma 1, each flight contains enough information to characterize a process of dimension at least $\tilde{l}+\mathrm{rank}\ {F}^{(0)}$. To efficiently cover a wide range of t values we use flights whose base t values increase exponentially, {0, 1, 2, 4, ..., 2ⁱ}. We then add flights as needed to support the Ho–Kalman structure, equation (8). The result is that the experiments performed are those with t values in the set

$$\begin{eqnarray}&&{ \mathcal T }={{ \mathcal T }}_{R}+{{ \mathcal T }}_{C},\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&{{ \mathcal T }}_{R}=\{0,1,2,4,\,\ldots ,\,{2}^{a}\}+\{0,1,2,\,\ldots ,\,\tilde{l}\}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{{ \mathcal T }}_{C}=\{0,1,2,4,\,\ldots ,\,{2}^{b}\}+\{0,1,2,\,\ldots ,\,\tilde{l}\}.\end{eqnarray} \tag{ 11 }$$

Here + is the sumset operation and $a,b,\tilde{l}$ are chosen integers. We found that such a pattern of experiments is generally necessary to obtain accurate models, especially for systems in which the non-Markovian aspects manifest only over long timescales. Notably, the characterization accuracy grows as $\max { \mathcal T }$, whereas the number of different experiments to be performed grows asymptotically as only ${(\mathrm{logmax}{ \mathcal T })}^{2}$.

Zoom In Zoom Out Reset image size

Once the experiments are performed, the data is analyzed to infer a model of the process ${ \mathcal P }$, the initial states, and the measurements. Our basic approach is to seek the model with the highest likelihood for the given observations. This is a challenging task for two reasons: (1) The size of the model (number of dimensions) is unknown. (2) The likelihood is extremely nonlinear in T: T appears with powers up to $\max { \mathcal T }$, which in practice can be on the order of 10³ or more. This creates an extremely unforgiving optimization landscape with many local (and very suboptimal) maxima.

Regarding the first challenge, a number of techniques have been devised to compare the likelihood of models with different numbers of parameters (e.g. the Akaike information criterion [70]) or to construct models whose complexity automatically scales with the complexity of the data (such as Dirichlet processes [71]). We use a simple but effective technique to first estimate the dimensionality of the process in question, which then determines the set of parameters to be estimated. The second challenge is essentially the same as that faced in gate set tomography. Building upon the progressive fitting approach described in [7], we devised a 4-stage model inference procedure, carefully honed over the course of our studies to reliably find concise, accurate process models. An overview of our procedure is given below; additional details are given in appendices B–D.

• 1.  
  Estimate the Model Dimension. First, the data is arranged into a matrix $\tilde{H}$ which is the experimental estimate of a Ho–Kalman matrix H. $\tilde{H}$ has the structure given in equation (8) with rows indexed by values from ${{ \mathcal T }}_{R}$ and columns indexed by values from ${{ \mathcal T }}_{C}$. This may be viewed as an arrangement of the data from all the flights as shown in figure 3. The effective dimension d is then estimated as the number of statistically significant non-zero singular values of $\tilde{H}$ [72], where significance is assessed by a chi-square test. In brief, the uncertainties in the experimental data are estimated and propagated to yield uncertainties in the singular values of $\tilde{H}$. This yields for each k = 1, 2, ... a threshold for the sum of squares of the k smallest singular values. If for some k the sum falls below the threshold, those singular values are discounted (deemed insignificant). The dimension d is taken to be the number of singular values that remain after the largest subset of singular values has been discounted. We note that this procedure does not determine the inherent dimension of a process so much as the dimension that is revealed by the data.
• 2.  
  Obtain an Initial Model. This stage is a variation of the Ho–Kalman method, modified to account for uncertainty in the experimental data. Let ${\tilde{H}}^{{\prime} }$ be the 'time shifted' version of $\tilde{H}$, i.e. the matrix obtained by increasing the t value of each element of $\tilde{H}$ by 1. First, the variance of each experimental quantity ${\tilde{F}}_{i,m}^{(t)}$ is estimated. Then a rank-d factorization $\tilde{H}\approx {LR}$ is found by minimizing ${\sum }_{i,m}{W}_{i,m}({LR}-\tilde{H}{)}_{i,m}^{2}$ where W_i,m is the inverse variance of ${\tilde{H}}_{i,m}$. The process matrix is estimated as the matrix T that minimizes ${\sum }_{i,m}{W}_{i,m}^{{\prime} }({LTR}-{\tilde{H}}^{{\prime} }{)}_{i,m}^{2}$ where ${W}_{i,m}^{{\prime} }$ is the inverse variance of ${\tilde{H}}_{i,m}^{{\prime} }$. The first few rows of L and columns of R may be taken as initial estimates of S and P, respectively.
• 3.  
  Obtain an Improved Model. Although the previous stage uses all the data, it uses only the fact that $\tilde{H}$ and ${\tilde{H}}^{{\prime} }$ are related by a (single) factor of T. It does not use the fact that data from different flights should be related by higher powers of T. To incorporate this relationship into the model estimate, the data is organized into groups of flights or 'blocks' with exponentially increasing t values (figure 3(c)). Starting with the matrices L, T, R obtained in stage 2, a search is performed to find matrices A, T, B that best fit the data in all the blocks. The first few rows of A and columns of B may be taken as improved estimates of S and P, respectively. To ensure that the search finds a good solution and does not get stuck in a local extremum, the search is performed progressively, starting with the lowest-t blocks and gradually incorporating data from blocks with higher t values. The low-t blocks yield coarse estimates of T which are gradually refined using the data in higher-t blocks.
• 4.  
  Obtain a physically valid model estimate. The model produced by stage 3 is usually fairly accurate, but tends to make slightly unphysical predictions (i.e. probabilities outside the interval [0,1]). In this last stage, the model likelihood is maximized with the addition of penalty terms to encourage the validity of predicted probabilities. As in stage 3, the log-likelihood is approximated by a weighted sum of squared errors; however, this time the weights are not based on the (fixed) experimental frequencies but are based on the predicted probabilities. Also, in this step the data is not arranged into Ho–Kalman matrices; the model ${F}^{(t)}={{ST}}^{t}P$ (equation (2)) is fit directly to each experiment.

In our simulation studies, the model inference procedure described above proved to be both fast and reliable. In over 99% of the thousands of simulations we performed, it found a high-likelihood model without any intervention. In the few cases in which it did not, making minor adjustments to optimization parameters in stage 3 or 4 (e.g. penalty weights, convergence thresholds) invariably led to success. We found that all four stages were usually necessary to obtain accurate models. On a standard desktop computer, the time to complete all four stages ranged from just a few seconds for 7 dimensional processes to a few minutes for a 19-dimensional process.

We note that the exact arrangement of the data in steps 1–3 does not appear be critical. Beyond satisfying the Ho–Kalman structure, it appears to be important only that each row and column of $\tilde{H}$ cover both a wide range of t values and multiple runs of roughly l consecutive t values.

To illustrate the QPI procedure let us consider how Alice would use it to characterize a 'black box' version of the slot machine described in section 2. In this version the slot machine has a display screen, a lever, and buttons labeled I0, I1, and M. Pushing I0 or I1 turns the machine on. When the lever is pulled the machine makes a momentary noise but otherwise has no observable effect. Pressing M causes the screen to briefly display the number 0 or 1, after which the machine turns off. In terms of the QPI formalism the initialization procedures are ${ \mathcal I }=\{{\mathsf{I}}{\mathsf{0}},{\mathsf{I}}{\mathsf{1}}\}$. Pulling the lever performs an unknown process ${ \mathcal P }$. Measurement is accomplished by pressing M. This is treated as two procedures, M0 and M1, which test whether the outcome is 0 or 1, respectively. Thus ${ \mathcal M }=\{{\mathsf{M}}{\mathsf{0}},{\mathsf{M}}{\mathsf{1}}\}$.

The behavior of the machine, which Alice will attempt to discover, is as follows: If the machine is turned on using I0, the output after an even number of lever pulls is always 0; after an odd number of pulls it is 0 with probability 0.5 and 1 with probability 0.5. Likewise, if the machine is turned on using I1, the output after an even number of lever pulls is always 1; after an odd number of pulls it is 0 with probability 0.5 and 1 with probability 0.5. This behavior is expressed by the series of matrices F⁽⁰⁾, F⁽¹⁾, F⁽²⁾, ...with

$$\begin{eqnarray}{F}^{(t)}\equiv \left[\begin{array}{cc}{F}_{{\mathsf{I}}{\mathsf{0}},{\mathsf{M}}{\mathsf{0}}}^{(t)} & {F}_{{\mathsf{I}}{\mathsf{0}},{\mathsf{M}}{\mathsf{1}}}^{(t)}\\ {F}_{{\mathsf{I}}{\mathsf{1}},{\mathsf{M}}{\mathsf{0}}}^{(t)} & {F}_{{\mathsf{I}}{\mathsf{1}},{\mathsf{M}}{\mathsf{1}}}^{(t)}\end{array}\right]=\left\{\begin{array}{ll}\left[\begin{array}{cc}1 & 0\\ 0 & 1\end{array}\right] & t\ \mathrm{even}\\ \left[\begin{array}{cc}0.5 & 0.5\\ 0.5 & 0.5\end{array}\right] & t\ \mathrm{odd}\end{array}\right..\end{eqnarray} \tag{ 12 }$$

Alice seeks a model $(S,T,P)$ such that ${{ST}}^{t}P={F}^{(t)}$ for all t. The rows of $S\in {{\mathbb{R}}}^{2\times d}$ will express the states produced by I0 and I1 respectively in some d-dimensional state space. $T\in {{\mathbb{R}}}^{d\times d}$ will express the action of the lever in this state space; and the columns of $P\in {{\mathbb{R}}}^{d\times 2}$ will express the properties measured by M0 and M1 respectively.

To perform QPI Alice first estimates the probabilities ${F}_{i,m}^{(t)}$ for all $i\in \{{\mathsf{I}}{\mathsf{0}},{\mathsf{I}}{\mathsf{1}}\}$, all $m\in \{{\mathsf{M}}{\mathsf{0}},{\mathsf{M}}{\mathsf{1}}\}$, and for selected values of t. Each experiment involves initializing the machine using either I0 or I1 and pulling the lever some fixed number of times t before performing readout M. She repeats each experiment 100 times and records the fraction of times outcome m = {0, 1} occurs as ${\tilde{F}}_{i,m}^{(t)}\approx {F}_{i,m}^{(t)}$. Alice first obtains

$$\begin{eqnarray}{\tilde{F}}^{(0)}=\left[\begin{array}{cc}1 & 0\\ 0 & 1\end{array}\right].\end{eqnarray} \tag{ 13 }$$

Since ${\tilde{F}}^{(0)}\approx {F}^{(0)}$ she surmises that $\mathrm{rank}\ {F}^{(0)}=2$, i.e. the initializations and readout span a 2 dimensional space. Alice decides that she is willing to construct a model that can capture at least one latent dimension, $\tilde{l}=1$. (In practice she should probably consider a larger value, but here we use a small value for simplicity of presentation.) Each flight of experiments will involve $2\tilde{l}+1=3$ consecutive t values. She chooses the flights to have base values from 0 up to 2¹ + 2¹ = 4. The complete set of t values used in Alice's characterization is ${ \mathcal T }={{ \mathcal T }}_{R}+{{ \mathcal T }}_{C}=\{0,1,2,3,4,5,6\}$ where ${{ \mathcal T }}_{R}={{ \mathcal T }}_{C}=\{0,{2}^{0},{2}^{1}\}+\{0,1\}$. For the sake of discussion, we suppose Alice obtains

$$\begin{eqnarray}{\tilde{F}}^{(0)}={\tilde{F}}^{(2)}={\tilde{F}}^{(4)}={\tilde{F}}^{(6)}=\left[\begin{array}{cc}1 & 0\\ 0 & 1\end{array}\right],\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}{\tilde{F}}^{(1)}=\left[\begin{array}{cc}0.48 & 0.52\\ 0.53 & 0.47\end{array}\right],{\tilde{F}}^{(3)}=\left[\begin{array}{cc}0.49 & 0.51\\ 0.46 & 0.54\end{array}\right],{\tilde{F}}^{(5)}=\left[\begin{array}{cc}0.52 & 0.48\\ 0.57 & 0.43\end{array}\right].\end{eqnarray} \tag{ 15 }$$

Having obtained experimental estimates of F^(t) for a number of t values, Alice proceeds to perform the four steps of data analysis:

• 1.  
  Estimate the model dimension. Alice arranges nearly all the data into a Ho–Kalman matrix

  $$\begin{eqnarray}\tilde{H}=\left[\begin{array}{ccc}{\tilde{F}}^{(0)} & {\tilde{F}}^{(1)} & {\tilde{F}}^{(2)}\\ {\tilde{F}}^{(1)} & {\tilde{F}}^{(2)} & {\tilde{F}}^{(3)}\\ {\tilde{F}}^{(2)} & {\tilde{F}}^{(3)} & {\tilde{F}}^{(4)}\\ {\tilde{F}}^{(3)} & {\tilde{F}}^{(4)} & {\tilde{F}}^{(5)}\end{array}\right].\end{eqnarray} \tag{ 16 }$$

  (The column that would contain ${\tilde{F}}^{(6)}$ is excluded because step 2 will require all the t values to be increased by one, and no data was taken for t = 7.) Using the procedure detailed in appendix B, Alice computes the singular values of $\tilde{H}$, the cumulative sum of squares of the singular values from smallest to largest, and the corresponding chi-square thresholds:

  

  Since the sum of squares falls below the threshold for d ≥ 4 but above the threshold for d ≤ 3, Alice concludes d = 3.
• 2.  
  To obtain an initial estimate, Alice shifts the entries in $\tilde{H}$ by one lever pull to obtain

  $$\begin{eqnarray}{\tilde{H}}^{{\prime} }=\left[\begin{array}{ccc}{\tilde{F}}^{(1)} & {\tilde{F}}^{(2)} & {\tilde{F}}^{(3)}\\ {\tilde{F}}^{(2)} & {\tilde{F}}^{(3)} & {\tilde{F}}^{(4)}\\ {\tilde{F}}^{(3)} & {\tilde{F}}^{(4)} & {\tilde{F}}^{(5)}\\ {\tilde{F}}^{(4)} & {\tilde{F}}^{(5)} & {\tilde{F}}^{(6)}\end{array}\right].\end{eqnarray} \tag{ 18 }$$

  She numerically obtains a rank-3 factorization ${LR}\approx \tilde{H}$ that minimizes the total weighted error. The first 2 rows of L yield an initial estimate of S while the first two columns of R yield an initial estimate of P. She obtains an initial estimate of T by minimizing total weighted error of the approximation ${LTR}\approx {\tilde{H}}^{{\prime} }$, keeping L and R fixed. Her initial model estimate is

  $$\begin{eqnarray}S=\left[\begin{array}{ccc}-0.659 & 0.705 & -0.010\\ -0.657 & -0.709 & -0.016\end{array}\right],\ T=\left[\begin{array}{ccc}1.00 & 0.015 & 0.013\\ 0.005 & -0.064 & -0.839\\ 0.006 & -1.19 & 0.057\end{array}\right],\ \tilde{P}=\left[\begin{array}{cc}-0.761 & -0.758\\ 0.706 & -0.708\\ -0.016 & 0.009\end{array}\right].\end{eqnarray} \tag{ 19 }$$
• 3.  
  To obtain a better estimate, Alice arranges her data into three blocks:

  $$\begin{eqnarray}{D}_{0}=\left[\begin{array}{cc}{\tilde{F}}^{(0)} & {\tilde{F}}^{(1)}\\ {\tilde{F}}^{(1)} & {\tilde{F}}^{(2)}\\ {\tilde{F}}^{(2)} & {\tilde{F}}^{(3)}\\ {\tilde{F}}^{(3)} & {\tilde{F}}^{(4)}\end{array}\right],\qquad {D}_{1}=\left[\begin{array}{cc}{\tilde{F}}^{(1)} & {\tilde{F}}^{(2)}\\ {\tilde{F}}^{(2)} & {\tilde{F}}^{(3)}\\ {\tilde{F}}^{(3)} & {\tilde{F}}^{(4)}\\ {\tilde{F}}^{(4)} & {\tilde{F}}^{(5)}\end{array}\right],\qquad {D}_{2}=\left[\begin{array}{cc}{\tilde{F}}^{(2)} & {\tilde{F}}^{(3)}\\ {\tilde{F}}^{(3)} & {\tilde{F}}^{(4)}\\ {\tilde{F}}^{(4)} & {\tilde{F}}^{(5)}\\ {\tilde{F}}^{(5)} & {\tilde{F}}^{(6)}\end{array}\right]\end{eqnarray} \tag{ 20 }$$

  Alice then uses the numerical search procedure described in appendix C to find $A\in {{\mathbb{R}}}^{8\times 3}$, $B\in {{\mathbb{R}}}^{3\times 4}$, and $T\in {{\mathbb{R}}}^{3\times 3}$ such that AB ≈ D₀, ATB ≈ D₁, and ${{AT}}^{2}B\approx {D}_{2}$. For the starting values of A, T, B she uses L,T, and the first four columns of P from the previous step. This yields

  $$\begin{eqnarray}S=\left[\begin{array}{ccc}-0.659 & 0.705 & -0.014\\ -0.657 & -0.709 & 0.003\end{array}\right],\ T=\left[\begin{array}{ccc}1.00 & 0.015 & 0.013\\ 0.005 & -0.064 & -0.834\\ 0.006 & -1.19 & 0.057\end{array}\right].\ P=\left[\begin{array}{cc}-0.761 & -0.758\\ 0.706 & -0.708\\ -0.017 & 0.011\end{array}\right].\end{eqnarray} \tag{ 21 }$$
• 4.  
  Alice improves her estimate even further using the procedure described in appendix D to directly minimize the total error of the model ${{ST}}^{t}P={F}^{(t)}$ over all $t\in { \mathcal T }$, using the values of S, T, P obtained from step 3. This yields her final estimated model

  $$\begin{eqnarray}S=\left[\begin{array}{ccc}-0.659 & 0.701 & -0.025\\ -0.657 & -0.712 & 0.000\end{array}\right]\ T=\left[\begin{array}{ccc}1.00 & 0.010 & 0.009\\ 0.005 & -0.061 & -0.840\\ 0.006 & -1.19 & 0.055\end{array}\right]\ P=\left[\begin{array}{cc}-0.764 & -0.755\\ 0.706 & -0.708\\ -0.020 & 0.013\end{array}\right].\end{eqnarray} \tag{ 22 }$$

For comparison, an exact minimal model (S, T, P) for this system is

$$\begin{eqnarray}{S}^{{\prime} }=\left[\begin{array}{ccc}-0.5 & -0.5 & -0.707\\ 0.5 & 0.5 & -0.707\end{array}\right],\qquad {T}^{{\prime} }=\left[\begin{array}{ccc}-1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{array}\right],\qquad {P}^{{\prime} }=\left[\begin{array}{cc}-0.5 & 0.5\\ -0.5 & 0.5\\ -0.707 & -0.707\end{array}\right].\end{eqnarray} \tag{ 23 }$$

Alice's final model does not appear to be close to this. However, recall that models are meaningful only up to similarity transformations. A suitable similarity transformation on Alice's final model yields

$$\begin{eqnarray}\begin{array}{c}S=\left[\begin{array}{ccc}-0.486 & -0.498 & -0.707\\ 0.511 & 0.493 & -0.707\end{array}\right],\,T=\left[\begin{array}{ccc}-1.003 & 0.004 & 0.000\\ -0.004 & 0.997 & 0.000\\ 0.006 & 0.000 & 1.000\end{array}\right]\,P=\left[\begin{array}{cc}-0.503 & 0.503\\ -0.497 & 0.497\\ -0.710 & -0.704\end{array}\right].\end{array}\end{eqnarray} \tag{ 24 }$$

which clearly approximates equation (23). As it turns out, this model is accurate only for t ≲ 10, a consequence of the fact that Alice tested only small values of t and repeated each experiment relatively few times.

In quantum computing technologies, operations on qubits are typically accomplished by applying control pulses that temporarily modify the Hamiltonian to induce a unitary rotation of one or more qubits. In many implementations the amount of rotation is proportional to the pulse energy. A miscalibration of the pulse intensity or duration results in over- or under-rotation of the qubits, which generally contributes to computational error [73].

In this example the process of interest is an intended π rotation of a qubit about the y axis, i.e. a bit flip. We suppose however that the intensity of the control pulse drifts sinusoidally over time, such that the actual angle of rotation produced by the pulse at time t (where each pulse takes one unit of time) is

$$\begin{eqnarray}&&{\theta }_{t}=\pi +\epsilon \sin {\rm{\Omega }}t.\end{eqnarray} \tag{ 25 }$$

We choose  = 0.01 and Ω = 0.02, which yields a small, slow oscillation of the qubit about the intended state in the x–z plane. This oscillation cannot be described by any fixed Markovian process involving only the qubit. In fact, the error can be described as phase modulation, a process which has an infinite number of sidebands. In other words, this process technically has an infinite dimension. However, relatively few dimensions are needed to accurately reproduce typical experimental data.

The evolution of the qubit under this process was simulated for two different initial states, $| +{\rangle }_{z}$ and $| +{\rangle }_{x}$. For each initial state, tomographic measurements of the qubit were simulated at 304 times ranging from t = 0 to t = 1035, forming 57 flights of length 12. Each measurement was repeated 10⁴ times. In nearly all realizations, model inference yielded an 11-dimensional model. Figure 4 shows the average characterization error as a function of t.

Zoom In Zoom Out Reset image size

Leakage refers to the loss of a qubit, or more precisely, the transition of a qubit to a non-computational state. This process is usually modeled as irreversible and incoherent, in which case the qubit has no effect on any subsequent operation. More realistically, qubit operations create and manipulate superpositions of computational and non-computational states [26]. In such cases a qubit can effectively leave and later return to the computational subspace, interfering with the computational evolution. The particular details of the qubit's behavior will depend strongly on the physics of the physical implementation of the qubit and the control mechanisms employed. For this example we use the 3-level anharmonic oscillator model which is applicable to popular qubit technologies including superconducting flux qubits and trapped ions [74]. In this model, the control pulse drives not only the intended $| 0\rangle \leftrightarrow | 1\rangle$ transition, but also off-resonant transitions between $| 1\rangle$ and the next excited state $| 2\rangle$. The interaction Hamiltonian is

$$\begin{eqnarray}&&H={\rm{\Omega }}(t)\left(| 0\rangle \langle 1| +\sqrt{2}| 1\rangle \langle 2| +{\rm{h}}.{\rm{c}}.\right)+{\rm{\Delta }}| 2\rangle \langle 2| ,\end{eqnarray} \tag{ 26 }$$

where Ω(t) is the control amplitude and Δ is the detuning of the non-computational state. For simulations we used Gaussian control pulses of width 0.25 time steps (truncated to have 1 time step duration) and chose an excited state detuning Δ = 20 yielding a system in which the leakage in and out of the computational space is small but not negligible. For these parameter values the total probability of state $| 2\rangle$ oscillates with a period of nearly 30 time steps, with maximal value of a few percent. Meanwhile, the computational state slowly precesses about the intended state.

The evolution of this 3-level system was simulated for four different initial qubit states ($| +{\rangle }_{z}=| 0\rangle$, $| -{\rangle }_{z}=| 1\rangle$, $| +{\rangle }_{x}$, and $| +{\rangle }_{y}$). For each initial state, tomographic measurements of the qubit were simulated at 196 times ranging from t = 0 to t = 1029, forming 57 flights of length 6. Each measurement was repeated 10⁴ times. In nearly all cases, the model inference yielded a 7-dimensional model. The average characterization error for this process is shown in figure 5.

Zoom In Zoom Out Reset image size

Another potential source of error in qubit devices is unintended coupling between qubits, or between a qubit and a nearby impurity. A simple model for such coupling is the isotropic spin exchange Hamiltonian

$$\begin{eqnarray}&&H=\gamma \left({\sigma }_{x}\otimes {\sigma }_{x}+{\sigma }_{y}\otimes {\sigma }_{y}+{\sigma }_{z}\otimes {\sigma }_{z}\right).\end{eqnarray} \tag{ 27 }$$

Taking γ = 0.01 to model weak coupling, the evolution of a qubit and impurity spin was simulated for three initial qubit states ($| +{\rangle }_{x}$, $| +{\rangle }_{y}$, and $| +{\rangle }_{z}$). The impurity is always initially in the state $| +{\rangle }_{z}$. For each initial state, tomographic measurements were simulated at 64 times ranging from t = 0 to t = 1030, forming 12 flights of length 7. Each measurement was repeated 10⁴ times. In nearly all cases, model inference yielded a 7-dimensional model. This may be compared with the nominal dimension 4² = 16 of a two qubit system. The average characterization error for this process is shown in figure 6.

Zoom In Zoom Out Reset image size

For our final case study, we considered a ${}^{31}{\rm{P}}$ donor in ²⁸Si with ²⁹Si impurities. The nuclear spin of the ${}^{31}{\rm{P}}$ defines a qubit. The qubit is manipulated by a time-varying magnetic field which also affects the impurities. Meanwhile, the ³¹P and ²⁹Si impurities interact via hyperfine coupling to the donor electron. A detailed model for this system was formulated in [75]. We simulated the behavior of the qubit and impurities in response to a sequence of pulses that each nominally produce a π rotation of the qubit (bit flip). For these simulations, the donor was imagined to be at the center of a ${(5\mathrm{nm})}^{3}$ cube with a ²⁹Si concentration of 800 ppm, which translates to 6 impurities in the qubit volume. The impurities were distributed randomly throughout the volume; in all, 200 random configurations were simulated. The background magnetic field B_z was 1.5 T, the amplitude B_x of the driving magnetic field was 0.15 T, and the ambient temperature was 250 mK. The total length of each pulse sequence was $2.9\,\mu {\rm{s}}$. (For further details see [75].)

For each configuration of the impurities, the evolution of the donor and impurities was simulated for 4 different initial qubit states ($| +{\rangle }_{x}$, $| +{\rangle }_{y}$, $| +{\rangle }_{z}$, and $| -{\rangle }_{z}$) and with the impurities initially in thermal equilibrium. For each initial state, tomographic measurements of the qubit were simulated at 272 different time steps ranging from t = 0 to t = 1033, forming 57 flights of length 10. Each measurement was repeated 10³ times. In most realizations the inferred model dimension was 19 or 20, which is much less than the nominal dimension of the ³¹P + ²⁹Si system (7 nuclear spins + 1 electron spin = 8 spins with a nominal dimension of ${\left({2}^{8}\right)}^{2}=65536$). Thus QPI reveals that the effective interaction between the qubit and the impurities involves only a small subspace of the available Hilbert space, an insight that is not readily apparent from the underlying model. The average characterization error for this system is shown in figure 7.

Zoom In Zoom Out Reset image size