When quantum tomography goes wrong: drift of quantum sources and other errors

The goal of quantum-state tomography [1] is to give a statistically reliable estimate of a quantum state ρ. Two further questions may come to mind: (i) what is the purpose of that estimate ρ? and (ii) why or when are we correct in giving an estimate of just one quantum state?

There are at least two answers to the first question: our experiment may be aimed at producing a particular state, say, a cluster state, and we may just want to verify how close ρ is to the desired state. But that answer provides really just an intermediate goal. The ultimate goal is always to use the desired state for some particular quantum information processing task. So we could say that the goal of producing an estimate ρ is to be able to predict the future performance in a particular protocol of one or more unmeasured quantum system(s) produced by the same source.

Now there is a nice statistical method for ranking different models according to their ability to predict future measurement results (not on how well they fit the past data!), based on the Akaike information criterion (AIC) [2]. That criterion was developed entirely within a classical context, but it ought to apply to quantum-state estimation, too. We show this is true, even though we will point out some interesting differences between classical and quantum statistics.

The motivation behind the second question is as follows. Since we do not have full control over all physical quantities relevant to the quantum-state generation process (for example, even the best laser suffers from phase diffusion; and there are always spatially and temporally fluctuating magnetic and electric fields), the quantum states produced by a quantum source are not all identical. A possible description of the individual states of M systems k = 1,...,M would be a sequence {ρ_k,k = 1,...,M} where each ρ_k+1 is a little different from the previous one (even with entanglement or correlation between the different systems, we can define ρ_k by tracing out all the other systems). So, why would we use just a single estimate ρ in this case? One aspect of the answer is, of course, that we have no way of estimating each individual ρ_k. A more positive answer is that multiple measurements of a given observable $\hat {O}$ only yield estimates of average quantities such as $\left \langle O\right \rangle =\overline {\mathrm {Tr}\,\rho _k\hat {O}}$ or $p_n=\overline {\mathrm {Tr}\,\rho _k| O_n\rangle \!\langle O_n |}$, where the average is over those k on which $\hat {O}$ was measured, and where $\left | O_n \right \rangle$ denotes an eigenstate of $\hat {O}$. These averages being linear in ρ_k are determined by a single density matrix, namely the average density matrix $\rho =\overline {\rho _k}$. This simple picture has been made much more rigorous by Renner in [3]. He showed that the crucial ingredient (missing in the simple picture) is permutation invariance. That is, if we randomly permute the sequence of quantum systems, and then trace out some subset, the joint state of the remaining systems is to a good approximation independently and identically distributed (i.i.d.). In our context this means that as long as the quorum of observables is measured in a random order, then to a good approximation any one of the remaining unmeasured systems can be described by a single density matrix ρ. We now discuss what may go wrong if we measure the observables constituting a quorum in a nonrandom order.

It is much easier to measure a given observable from the quorum many times in a row, before switching to measurement of the next observable. Such a procedure is standard practice, but it voids Renner's proof, and so it may be that there is not a single density matrix that can be validly assigned to the remaining unmeasured quantum systems.

Let us introduce this problem with a simple example. Given an ensemble of 3N ≫ 1 qubits that—we assume!—are identically and independently prepared, we want to estimate their density matrix. So we divide them into three equal and sequential groups, and measure σ_x on samples 1,...,N, σ_y on samples N + 1,...,2N, and σ_z on the last N. Now, if the samples are indeed identically prepared in some state ρ, then we can safely perform the measurements in this order—the state ρ^⊗3N is invariant under permutations, so all orderings are equivalent. But if the source is drifting over time, the first N copies are best described by a mean density matrix $\bar {\rho }_{1}$, while the second and third sets of N qubits are best described by (possibly different) average states $\bar {\rho }_{2}$ and $\bar {\rho }_{3}$, respectively.

For an amusing (albeit extreme) example, consider a situation where the first N copies are best described by $\bar {\rho }_{1} = | +\rangle \!\langle + |$, the second group by $\bar {\rho }_{2}=| +i\rangle \!\langle +i |$, and the third by $\bar {\rho }_{3}=| 0\rangle \!\langle 0 |$. The measurement outcomes in this case are not random at all: every single measurement (of σ_x, σ_y, σ_z) will yield eigenvalue +1. Linear inversion tomography will yield a radically nonpositive state

$$\begin{equation} \hat{\rho}_{\mathrm{tomo}} = \left(\begin{array}{@{}cc@{}} 1 & \frac{1+\mathrm{i}}{2} \\ \frac{1-\mathrm{i}}{2} & 0 \end{array}\right), \end{equation} \noindent \tag{ 1 }$$

and maximum likelihood estimation (MLE) yields the projector onto $\hat {\rho }_{\mathrm {tomo}}$'s positive eigenspace. Although both estimates are plausible answers to 'what single matrix best fits the observed data?', neither one of them is of any predictive use at all! The source is drifting so rapidly and drastically that this set of 3N samples really tells us almost nothing about future observations. This is the simplest and best conclusion at which our data analysis should arrive.

This is a rather extreme and contrived example of experimental drift (below we will discuss a more common type of nonrandom experiment where the above cycle of measurements is repeated once: so we measure σ_x on the first N/2 copies, then σ_y, then σ_z, and then σ_x,σ_y,σ_z again, each on N/2 sequential copies). More realistic examples show similar behavior, though. The statistics given above are actually more consistent with a different (and still plausible) mechanism: when the measurement apparatus is 'rotated' to perform a different measurement, the experimenter inadvertently 'rotates' the samples as well. A particularly naïve version of this could occur with photon polarization, where one way to physically rotate a polarizer is for the experimentalist to simply rotate his own frame of reference (e.g., by lying down). Such a passive rotation obviously fails to change the relative orientation of samples and apparatus. More realistic examples occur when similar quantum gate devices are used to (i) prepare states (e.g. EPR states) and (ii) implement measurements. In quantum process tomography, this sort of pitfall is well known; it violates the conditions for complete positivity of processes, and causes negative eigenvalues just as in our example above [4].

All of these failures are examples of a single phenomenon: sample-apparatus correlation. In process tomography, this is usually explained by correlation between the system and its environment. In state tomography, there is no environment per se, but if the state of the kth sample is (in any way) correlated with the behavior of the measurement apparatus (e.g. with what measurement it is oriented to perform), then tomography goes wrong. Experimental drift is a simple and easy to understand example: the sample state is correlated with time, and if the apparatus setting is also allowed to vary with time, then there will be sample-apparatus correlation. As noted above, this can be eliminated by explicitly randomizing the order of measurements, so that while the samples are still time-dependent, the apparatus is not. Other kinds of sample-apparatus correlation are not so easy to remedy.

In the example given above, the extremity of the data—and the fact that the linear inversion estimate is radically negative—are a dead giveaway. On the other hand, linear inversion can produce negative estimates even with ideal data [5, 6] because of statistical fluctuations. The raison d'etre of MLE is to fix this negativity, but by constraining the estimate to positive states, MLE also hides the tell-tale signature of failed tomography. Moreover, negative estimates are not (in general) a reliable symptom even of drastic experimental drift. If the drifting states in the example above were a bit more mixed—e.g. $\bar {\rho }_{k}' = \frac{1}{2}\bar {\rho }_{k} + \frac {1}{4}{\mathbbm{1}}$—then linear inversion and MLE would yield identical and positive density matrices. But, just as in the original example, those estimates would be useless and not predictive.

Fortunately, there is a general solution to this problem. It elegantly generalizes the observation (made above) that a radically negative $\hat {\rho }_{\mathrm {tomo}}$ should trigger skepticism. It can also diagnose drift in the absence of negativity if the data are sufficiently rich. It is called model selection.

The core principle is that, when tomography fails:

• 1.  
  The standard model for tomography—i.i.d. samples described by a single density matrix—is bad.
• 2.  
  Some other model will be better.
• 3.  
  We can quantify 'bad' and 'better', and use the results to decide whether our tomography went wrong.

Clearly, putting this into practice requires that we come up with alternative models to describe the data. Model design is more of an art than a science. Here, we demonstrate alternative models for some simple and relevant problems, and leave the rich problems of general and optimal alternative-model design to future work. Instead, we focus on model selection, which means determining whether (i) the standard tomographic model is pretty good, or (ii) some other model (e.g. a drifting source model) is better.

To accomplish this, we propose, as we mentioned above, to use the AIC [2]. Widely used outside of physics [7, 8], the AIC is relatively unknown within the physics community. However, it has been applied in astrophysics [9], entanglement verification [10] and quantum state estimation [11–13]. Its function is to quantify (by assigning a real number) how well a given model describes the data from a given experiment. The AIC's absolute value is not meaningful, but the relative AIC values for multiple different models have a deep and useful meaning (see following section for a more detailed discussion of the AIC, its meaning, and its derivation). Their simplest use is to rank all the different models, and thus to identify (a) which is the best, and (b) how significantly 'worse' the others are.

The AIC assigns a number Ω_k to each model k, given by⁵

$$\begin{equation} \Omega_k:=\ln {\cal L}_k -K_k, \end{equation} \tag{ 2 }$$

where ${\cal L}_k$ is the likelihood of model k—or, if model k has adjustable parameters (as is usually the case), the maximum of the likelihood over all those parameters—and K_k is the number of independent model parameters used in model k to fit the data⁶. The larger the AIC (Ω_k) is, the higher the model is ranked. While Ω_k's absolute value is meaningless, the difference Δ = Ω_k − Ω_k' represents (roughly speaking) the weight of evidence in favor of k over k', measured in bits. So, for example, if we want to report a weighted average of the two models, the ratio of the weights assigned to models k and k' should be $w_k/w_{k'} = \exp (\Omega _k-\Omega _{k'})$.

The AIC's simple form admits a simple interpretation: fitting the data better (higher likelihood) is good, but extra parameters are bad. Additional parameters must justify their existence by improving the likelihood (a measure of goodness-of-fit) by at least a factor of e. This helps to prevent overfitting. Adding adjustable parameters will always improve a model's fit—but a good fit to past data is not a guarantee that the model will accurately predicting future measurements.

Example. If we measure each of 3N qubits, measuring $\hat {O}_j$ on qubits j for j = 1,...,3N, then the best possible fit to the data is to assume that each qubit j just happened to be in the appropriate eigenstate of $\hat {O}_j$ so that the probability of the observed data is $\mathcal {L} = 1$!

Intuitively, this 'explanation' is absurd. The AIC quantifies that intuition; that model requires a huge (O(N)) number of parameters, and the resulting penalty will overwhelm its higher likelihood, ensuring that its AIC is far worse than that of simpler models.

To apply the AIC to our example of section 1.2, we need an alternative model (the 'standard model' just uses a single density matrix for all 3N qubits). A simple alternative that describes experimental drift (as well as some other forms of sample-apparatus correlation) is to use one density matrix for each of the three groups of samples. This alternative model will always fit the data at least as well, but it may use more parameters⁷. The AIC ranks both models, and quantifies how much better one is than the other. We perform and analyze this calculation for our single-qubit example (where just two models are sufficient) in section 2.1, and address more complicated variations on this theme—with multiple alternative models—in section 2.3.

To conclude this (long) introduction, we note that the appearance of maximum likelihoods in the AIC does not imply any privileged role for MLE estimation of states or any other physical quantities. The likelihood is a central concept in statistics, and appears in almost every method. In the AIC, it is used specifically to quantify goodness-of-fit, and (obviously) the AIC balances this quantity against another (model complexity). Moreover, the AIC is used only to rank different models. There is no implicit requirement that the highest-ranked model must be chosen exclusively (in fact, a common strategy is to average over high-ranked models), and even if the 'best' model is chosen, we remain free to analyze that model without MLE (e.g. via Bayesian averaging).

We return to tomography of single qubits, where we measure σ_x on the first N qubits, then σ_y on the next N, and σ_z on the last N qubits. Denote the three thusly observed averages by $X:=\left \langle \sigma _x\right \rangle _{{\mathrm { obs}}}$, $Y:=\left \langle \sigma _y\right \rangle _{{\mathrm { obs}}}$, and $Z:=\left \langle \sigma _z\right \rangle _{{\mathrm { obs}}}$. In order to calculate likelihoods, we need the frequencies of having observed spin up (+) and down (−), respectively. They are given in terms of these averages by

$$\begin{equation} f_{\pm}^{(x)}=\frac{1\pm X}{2}, \end{equation} \tag{ 3a }$$

$$\begin{equation} f_{\pm}^{(y)}=\frac{1\pm Y}{2}, \end{equation} \tag{ 3b }$$

$$\begin{equation} f_{\pm}^{(z)}=\frac{1\pm Z}{2}. \end{equation} \tag{ 3c }$$

$$\begin{equation} \ln{\cal L}^{(x)}=Nf_+^{(x)}\ln(f_+^{(x)})+Nf_-^{(x)}\ln(f_-^{(x)}) =-NH\left(\frac{1+X}{2}\right), \end{equation} \tag{ 4 }$$

with H(.) the Shannon entropy. The same story holds for the next two sets of measurements, and so there is always a perfect fit to the data when we use the 'alternative model' with three density matrices, and that model needs three independent parameters. We conclude that the AIC assigns the following ranking to the alternative model:

$$\begin{equation} \Omega_{\mathrm{a}}=-N\left\{ H\left(\frac{1+X}{2}\right)+H\left(\frac{1+Y}{2}\right)+H\left(\frac{1+Z}{2}\right)\right\}-3. \end{equation} \noindent \tag{ 5 }$$

The performance of the 'standard model' depends on the value of just one number. If

$$\begin{equation} R^2:=X^2+Y^2+Z^2\leqslant 1, \end{equation} \noindent \tag{ 6 }$$

there is a single maximum likelihood density matrix $\bar {\rho }$ (with purity $\mathrm {Tr}\,\bar {\rho }^2=(R^2+1)/2$) that describes the whole measurement perfectly, just as the alternative model does. The standard model also needs three parameters in this case, and the maximum likelihood is also the same as for the alternative model. So, in this case there is no real difference between the two models—we could pick $\bar {\rho }_{1}=\bar {\rho }_{2}=\bar {\rho }_{3}=\bar {\rho }$—and we have Ω_s = Ω_a. There is no reason to reject the standard model when R ⩽ 1.

Now let us suppose that R > 1. We have then the choice between two descriptions:

• 1.  
  Alternative model. We describe each of the three measurements by their own density matrix. The maximum likelihood estimates of those three states satisfy

  $$\begin{equation} \mathrm{Tr}\,\bar{\rho}_{1}\sigma_x=X, \end{equation} \tag{ 7a }$$

  $$\begin{equation} \mathrm{Tr}\,\bar{\rho}_{2}\sigma_y=Y, \end{equation} \tag{ 7b }$$

  $$\begin{equation} \mathrm{Tr}\,\bar{\rho}_{3}\sigma_z=Z. \end{equation} \noindent \tag{ 7c }$$

  Three independent parameters are needed for this model. ( $\bar {\rho }_{1},\bar {\rho }_{2},\bar {\rho }_{3}$ are underdetermined, of course, but for the purpose of finding the maximum likelihood ${\cal L}_{\mathrm {a}}$ the information suffices.)
• 2.  
  Standard model. We use one density matrix to describe all three measurements together. The maximum likelihood estimate of that state will be pure. There is no known method to compute it exactly, but a generally good approximation is given by

  $$\begin{equation} \mathrm{Tr}\,\bar{\rho}_{\mathrm{s}}\sigma_x=X/R, \end{equation} \tag{ 8a }$$

  $$\begin{equation} \mathrm{Tr}\,\bar{\rho}_{\mathrm{s}}\sigma_y=Y/R, \end{equation} \tag{ 8b }$$

  $$\begin{equation} \mathrm{Tr}\,\bar{\rho}_{\mathrm{s}}\sigma_z=Z/R, \end{equation} \noindent \tag{ 8c }$$

  and this state's likelihood is a strict (but generally pretty tight) lower bound on the maximum likelihood for the standard model. Two independent parameters are needed in this model⁸.

The reason we end up with a pure maximum likelihood state in the standard model is that the single matrix fitting the data perfectly lies outside the set of physical states (it has a negative eigenvalue), and the closest physical state lies on the boundary [5]. In the case of qubits, this means a pure state. More precisely, if the unphysical best-fit matrix $\tilde {\rho }$ is written in its diagonal form, $\tilde {\rho }=\sum _{k=+,-}\lambda _k| \psi _k\rangle \!\langle \psi _k |$, with λ₊ > 1 and λ₋ < 0, then the maximum likelihood estimate would be $\bar {\rho }_{\mathrm {s}}=| \psi _+\rangle \!\langle \psi _+ |$. The latter state has the properties (8), as can be easily verified by explicit calculation.

Thus, when R > 1 the alternative model fits the data better but uses one more parameter than does the standard model. We can calculate the maximum likelihoods analytically in each of the two models, and thus obtain the relative AIC score of the two models:

$$\begin{eqnarray} \Omega_{\mathrm{s}}-\Omega_{\mathrm{a}}&=&1+N\sum_{M=X,Y,Z} \frac{1}{2}\ln \frac{1-M^2/R^2}{1-M^2} +\frac{M}{2}\ln\frac{(R+M)(1-M)}{(R-M)(1+M)}. \end{eqnarray} \tag{ 9 }$$

We accept the standard model as consistent iff Ω_s ⩾ Ω_a. This will happen only if R is sufficiently close to 1. If we expand R around 1, we can Taylor expand the right-hand side of (9) as

$$\begin{equation} \Omega_{\mathrm{s}}-\Omega_{\mathrm{a}}\approx 1-N\sum_{M=X,Y,Z}\frac{(R-1)^2M^2}{2(1-M^2)}, \end{equation} \tag{ 10 }$$

provided (R − 1)² ≪ (1 − M²) for M = X,Y,Z. That is, with this proviso, the standard model is consistent only when

$$\begin{equation} (R-1)\leqslant \frac{C}{\sqrt{N}}, \end{equation} \tag{ 11 }$$

with the constant C given by

$$\begin{equation} C=\frac{1}{\sqrt{\sum_M M^2/2(1-M^2)}}. \end{equation} \tag{ 12 }$$

The dependence of the condition (11) on N agrees with the simple idea that it is sufficient for R to be less than about a standard deviation or two above 1 for the standard model to still apply, and that standard deviation, of course, decays like $1/\sqrt {N}$ for N → ∞.

The implementation of tomography in the previous example is probably too simple and too obviously wrong for it to have been applied in an actual experiment. The straightforward improvement to measure each of σ_x,σ_y,σ_z in two separate blocks will allow one to detect drift. Let us denote the six observed averages by $X_{1,2}:=\left \langle \sigma _x\right \rangle _{{\mathrm { obs 1,2}}}$, $Y_{1,2}:=\left \langle \sigma _y\right \rangle _{{\mathrm { obs 1,2}}}$, and $Z_{1,2}:=\left \langle \sigma _z\right \rangle _{{\mathrm { obs 1,2}}}$. Drift can be detected by comparing the pairs of estimates X_1,2 with each other, Y_1,2 with each other, and Z_1,2 with each other. The AIC works as follows: we need again at least two different models for describing the data. One will be the standard model, with one density matrix describing all six measurements. This density matrix will be determined by the three averages (X₁ + X₂)/2, etc. The alternative model may consist of two independent density matrices (with six parameters in total) or of two density matrices that are not independent with either four or five parameters in total. Let us test the AIC in a simulation of data generated by single-qubit states of the form

$$\begin{equation} \rho_{{\rm actual}}=p| \psi_\phi\rangle\!\langle \psi_\phi | +(1-p){\mathbbm{1}}/2, \end{equation} \tag{ 13 }$$

where the pure state $\left | \psi _\phi \right \rangle$ depends on an angle ϕ, which we assume to undergo a random walk, and with the following meaning:

$$\begin{eqnarray} \left\langle \psi_\phi \right|\sigma_x\left| \psi_\phi \right\rangle&=&\cos\phi,\nonumber\\ \left\langle \psi_\phi \right|\sigma_y\left| \psi_\phi \right\rangle&=&\sin\phi,\\ \left\langle \psi_\phi \right|\sigma_z\left| \psi_\phi \right\rangle&=&0.\nonumber \end{eqnarray} \noindent \tag{ 14 }$$

For p we take the value p = 0.9. We perform in total 3000 measurements, divided into six groups of 500, in which we measure σ_x,σ_y,σ_z,σ_x,σ_y,σ_z in that order.

In figures 1 and 2 we plot two qualitatively different cases. In the first case the diffusion of ϕ is so fast that it leads to noticeably different values of X_1,2 and Y_1,2. The AIC in this case gives a very clear preference for the alternative model of using two density matrices with five parameters in total (only the expectation value of σ_z does not change over the course of the experiment). In the second case the drift over the course of the experiment is small enough so that the standard model is still the best, even though there is some drift, and even though the more complicated model does, of course, fit the data slightly better.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Consider now a tomographically complete measurement on 9N copies of two qubits, where on the first N pairs of qubits we measure σ_x on both qubits independently, then on the next N pairs we measure σ_x on one qubit and σ_y on the other, then on the third set of N pairs we measure σ_x on the one and σ_z on the other ... until on the last (ninth) set of N pairs we measure σ_z on both qubits. The first measurement is described by three independent averages that are obtained from measuring σ_x on both qubits independently:

$$\begin{equation} XX:=\left\langle\sigma_x\otimes\sigma_x\right\rangle_{{\rm obs,1}}, \end{equation} \tag{ 15a }$$

$$\begin{equation} IX:=\left\langle{\mathbbm{1}}\otimes\sigma_x\right\rangle_{{\rm obs, 1}}, \end{equation} \tag{ 15b }$$

$$\begin{equation} XI:=\left\langle\sigma_x\otimes{\mathbbm{1}}\right\rangle_{{\rm obs, 1}}. \end{equation} \tag{ 15c }$$

$$\begin{equation} XY:=\left\langle\sigma_x\otimes\sigma_y\right\rangle_{{\rm obs, 2}}, \end{equation} \tag{ 16a }$$

$$\begin{equation} IY:=\left\langle{\mathbbm{1}}\otimes\sigma_y\right\rangle_{{\rm obs, 2}}, \end{equation} \tag{ 16b }$$

$$\begin{equation} XI':=\left\langle\sigma_x\otimes{\mathbbm{1}}\right\rangle_{{\rm obs, 2}}. \end{equation} \tag{ 16c }$$

Writing down all different averages obtained from this particular experiment, we find nine quantities that are measured once, and six other quantities that are measured thrice. It becomes now much harder to judge when all the differences between those different estimates of the same quantities are, in total, statistically significant or not. That is, the generalization of the ad hoc method that worked fine for a single qubit, becomes troublesome. This, of course, becomes exponentially worse for more than two qubits.

On the other hand, the AIC can be applied straightforwardly to various alternative models. It is sufficient to find just one alternative model superior to the standard model in order to have succeeded in diagnosing an inconsistency in our tomographic experiment. Of course there is a large multitude of alternative models, but one can be guided in searching for such models by looking for those estimates of the same quantities that are the least consistent.