Quantum tomography via compressed sensing: error bounds, sample complexity and efficient estimators

In recent years, there has been amazing progress in studying complex quantum mechanical systems under controlled laboratory conditions [1]. Quantum tomography of states and processes is an invaluable tool used in many such experiments, because it enables a complete characterization of the state of a quantum system or process (see, e.g., [2–16]). Unfortunately, tomography is very resource intensive, and scales exponentially with the size of the system. For example, a system of n spin-1/2 particles (qubits) has a Hilbert space with dimension d = 2ⁿ, and the state of the system is described by a density matrix with d² = 4ⁿ entries.

Recently, a new approach to tomography was proposed: compressed quantum tomography, based on techniques from compressed sensing [17, 18]. The basic idea is to concentrate on states that are well approximated by density matrices of rank r ≪ d. This approach can be applied to many realistic experimental situations, where the ideal state of the system is pure, and physical constraints (e.g. low temperature or the locality of interactions) ensure that the actual (noisy) state still has low entropy.

This approach is convenient because it does not require detailed knowledge about the system. However, note that when such knowledge is available, one can use alternative formulations of compressed tomography, with different notions of sparsity, to further reduce the dimensionality of the problem [19]. We will compare these methods in section 6.2.

The main challenge in compressed tomography is how to exploit this low-rank structure, when one does not know the subspace on which the state is supported. Consider the example of a pure quantum state. Since pure states are specified by only O(d) numbers, it seems plausible that one could be reconstructed after measuring only O(d) observables, compared with O(d²) for a general mixed state. While this intuition is indeed correct [20–23], it is a challenge to devise a practical tomography scheme that takes advantage of this. In particular, one is restricted to those measurements that can be easily performed in the laboratory; furthermore, one then has to find a pure state consistent with measured data [24], preferably by some procedure that is computationally efficient (note that, in general, finding minimum-rank solutions is NP-hard, hence computationally intactable [25]).

Compressed tomography provides a solution that meets all these practical requirements [17, 18]. It requires measurements of two-outcome Pauli observables, which are feasible in many experimental systems. In total, it uses a random subset of m = O(rd log d) Pauli observables, which is just slightly more than the O(rd) degrees of freedom that specify an arbitrary rank r state. Then the density matrix ρ is reconstructed by solving a convex program. This can be done efficiently using general-purpose semidefinite program (SDP) solvers [26] or specialized algorithms for larger instances [27–29]. The scheme is robust to noise and continues to perform well when the measurements are imprecise or when the state is only close to a low-rank state.

Here, we follow up on [17, 18] by giving a stronger (and completely different) theoretical analysis, which is based on the restricted isometry property (RIP) [30–32]. This answers a number of questions that could not be addressed satisfactorily using earlier techniques based on dual certificates. Firstly, how large is the error in our estimated density matrix when the true state is full-rank with decaying eigenvalues? We show that the error is not much larger than the 'tail' of the eigenvalue spectrum of the true state. Secondly, how large is the sample complexity of compressed tomography, i.e. how many copies of the unknown state are needed to estimate its density matrix? We show that compressed tomography achieves nearly optimal sample complexity among all procedures using Pauli measurements and, surprisingly, the sample complexity of compressed tomography is nearly independent of the number of measurement settings m as long as m ⩾ Ω(rd poly log d).

In addition, we use numerical simulations to investigate the question: given a fixed time T during which an experiment can be run, is it better to do compressed tomography or full tomography, i.e. is it better to use a few measurement settings and repeat them many times or do all possible measurements with fewer repetitions? For the situations we simulate, we find that compressed tomography provides better accuracy at a reduced computational cost compared to standard maximum-likelihood estimation (MLE).

Finally, we provide two useful tools: a procedure for certifying the accuracy of low-rank state estimates and a very simple compressed sensing technique for quantum process tomography.

We now describe these results in more detail.

Theoretical analysis using RIP. We use a fundamental geometric fact: the manifold of rank-r matrices in $\mathbb {C}^{d\times d}$ can be embedded into O(rd poly log d) dimensions, with low distortion in the two-norm. An embedding that does this is said to satisfy the RIP [32]. In [30], it was shown that such an embedding can be realized using the expectation values of a random subset of O(rd poly log d) Pauli matrices. This implies the existence of the so-called 'universal' methods for low-rank matrix recovery: there exists a fixed set of O(rd poly log d) Pauli measurements, which has the ability to reconstruct every rank-r d × d matrix. Moreover, with high probability, a random choice of Pauli measurements will achieve this. (The earlier results of [17] placed the quantifiers in the opposite order: for every rank-r d × d matrix ρ, most sets of O(rd poly log d) Pauli measurements can reconstruct that particular matrix ρ.)

Intuitively, the RIP says that a set of random Pauli measurements is sensitive to all low-rank errors simultaneously. This is important, because it implies stronger error bounds for low-rank matrix recovery [31]. These bounds show that, when the unknown matrix ρ is full-rank, our method returns a (certifiable) rank-r approximation of ρ that is almost as good as the best such approximation (corresponding to the truncated eigenvalue decomposition of ρ).

In [31], these error bounds were used to show the accuracy of certain compressed sensing estimators, for measurements with additive Gaussian noise. Here, we use them to upper-bound the sample complexity of our compressed tomography scheme. (That is, we bound the errors due to estimating each Pauli expectation value from a finite number of experiments.) Roughly speaking, we show that our scheme uses O(r²d² log d) copies to characterize a rank-r state (up to constant error in trace norm). When r = d, this agrees with the sample complexity of full tomography. Our proof assumes a binomial noise model, but minor modifications could extend this result to other relevant noise models, such as multinomial, Gaussian or Poissonian noise.

Furthermore, we show an information-theoretic lower bound for tomography of rank-r states using adaptive sequences of single-copy Pauli measurements: at least Ω(r²d²/log d) copies are needed to obtain an estimate with constant accuracy in the trace distance. This generalizes a result from [33] for pure states. Therefore, our upper bound on the sample complexity of compressed tomography is nearly tight, and compressed tomography nearly achieves the optimal sample complexity among all possible methods using Pauli measurements.

Our observation that incomplete sets of observables are often sufficient to unambiguously specify a state gives rise to a new degree of freedom when designing experiments: when aiming to reduce statistical noise in the reconstruction, one can either estimate a small set of observables relatively accurately or else a large (e.g. complete) set of observables relatively coarsely. Our bounds (as well as our numerics) show that, remarkably, over a very large range of m the only quantity relevant for the reconstruction error is t, the total number of experiments performed. It does not matter over how many observables the repetitions are distributed. Thus, fixing t and varying m, the reduction in the number of observables and the increase in the number of measurements per observable have no net effect with regard to the fidelity of the estimate as long as m ⩾ Ω(rd poly log d). This finding holds independently of whether fidelity or trace distance is used to measure the reconstruction error, and we believe that it is plausible that a similar phenomenon occurs for other cost functions as well.

Certification. We generalize the technique of direct fidelity estimation (DFE) [33, 34] to work with low-rank states. Thus, one can use compressed tomography to get an estimated density matrix $\hat {\rho }$, and use DFE to check whether $\hat {\rho }$ agrees with the true state ρ. This check is guaranteed to be sound even if the true state ρ is not approximately low rank. Our extension of DFE may be of more general interest, since it can be used to efficiently certify any estimate $\hat \rho$ regardless of whether it was obtained using compressed sensing or not, as long as the rank r of the estimate is small (and regardless of the 'true' rank).

Numerical simulations. We compare the performance of several different estimators (methods for reconstructing the unknown density matrix). They include: constrained trace-minimization (also known as the matrix Dantzig selector), least squares with trace-norm regularization (or the matrix Lasso), as well as a standard MLE [35–37] for comparison.

We observe that our estimators outperform MLE in essentially all aspects, with the matrix Lasso giving the best results. The fidelity of the estimate is consistently higher using the compressed tomography estimators. Also, the accuracy of the compressed sensing estimates is (as mentioned above) fairly insensitive to the number of measurement settings m (assuming that the total time available to run the experiment is fixed). So by choosing m ≪ d², one still obtains accurate estimates, but with much faster classical post-processing, since the size of the data set scales as O(m) rather than O(d²).

It may be surprising to the reader that we outperform MLE, since it is often remarked (somewhat vaguely) that 'MLE is optimal'. However, MLE is a general-purpose method that does not specifically exploit the fact that the state is low-rank. Also, the optimality results for MLE only hold asymptotically and for informationally complete measurements [38, 39]; for finite data [40] or for incomplete measurements, MLE can be far from optimal.

From these results, one can extract some lessons about how to use compressed tomography. Compressed tomography involves two separate ideas: (1) measuring an incomplete set of observables (i.e. choosing m ≪ d²) and (2) using trace minimization or regularization to reconstruct low-rank solutions. Usually, one does both of these things. Now, suppose that the goal is to reconstruct a low-rank state using as few samples as possible. Our results show that one can achieve this goal by doing (2) without (1). At the same time, there is no penalty in the quality of the estimate when doing (1), and there are practical reasons for doing it, such as reducing the size of the data set to speed up the classical post-processing.

Quantum process tomography. Finally, we adapt our method to perform tomography of processes with a small Kraus rank. Our method is easy to implement, since it requires only the ability to prepare eigenstates of Pauli operators and measure Pauli observables. In particular, we require no entangling gates or ancillary systems for the procedure. In contrast with [19], our method is not restricted to processes that are element-wise sparse in some known basis, as discussed in section 6.2. This is an important advantage in practice, because while the ideal (or intended) evolution of a system may be sparse in a known basis, it is often the case that the noise processes perturbing the ideal component are not sparse, and knowledge of these noise processes is key to improving the fidelity of a quantum device with some theoretical ideal.

1.1. Related work

1.2. Notation and outline

We describe our compressed tomography scheme in detail. Firstly, we describe the measurement procedure and the method for reconstructing the density matrix. Then we prove error bounds and analyze the sample complexity.

2.1. Random Pauli measurements

Consider a system of n qubits, and let d = 2ⁿ. Let $\mathcal {P}$ be the set of all d² Pauli operators, i.e. matrices of the form P = σ₁⊗···⊗σ_n where σ_i∈{I,σ^x,σ^y,σ^z}.

Our tomography scheme works as follows. Firstly, choose m Pauli operators, P₁,...,P_m, by sampling independently and uniformly at random from $\mathcal {P}$. (Alternatively, one can choose these Pauli operators randomly without replacement [65], but independent sampling greatly simplifies the analysis.) We will use t copies of the unknown quantum state ρ. For each i∈[m], take t/m copies of the state ρ, measure the Pauli observable P_i on each one and average the measurement outcomes to obtain an estimate of the expectation value Tr(P_iρ). (We will discuss how to set m and t later. Intuitively, to learn a d × d density matrix with rank r, we will set m ∼ rd log⁶ d and t ∼ r²d² log d.)

To state this more concisely, we introduce some notation. Define the sampling operator to be a linear map $\mathcal {A}:\: \mathbb {H}^d \rightarrow \mathbb {R}^m$ defined for all i∈[m] by

$$\begin{equation} (\mathcal{A}(\rho))_i = \sqrt{\frac{d}{m}} \mathrm{Tr}(P_i\rho) . \end{equation} \noindent \tag{ 1 }$$

(The normalization is chosen so that $\mathbbm {E} \mathcal {A}^*\mathcal {A} = \mathcal {I}$, where $\mathcal {I}$ denotes the identity superoperator, and $\mathcal {A}^*$ is the adjoint of $\mathcal {A}$. That is, $\mathcal {A}^*:\: \mathbb {R}^m \rightarrow \mathbb {H}^d$ is the complex-conjugate transpose of the linear map $\mathcal {A}$.) We can then write the output of our measurement procedure as a vector

$$\begin{equation} y = \mathcal{A}(\rho) + z, \end{equation} \noindent \tag{ 2 }$$

where z represents statistical noise due to the finite number of samples. (More generally, if the measurements contain systematic or adversarial noise, these can also be represented by z. There are various error bounds for this situation, although they are necessarily more pessimistic than the ones we show here [66].)

2.2. Reconstructing the density matrix

We now show two methods for estimating the density matrix ρ. Both are based on the same intuition: find a matrix $X \in \mathbb {H}^d$ that fits the data y while minimizing the trace norm $\Vert {X}\Vert _{\mathrm {tr}}$, which serves as a surrogate for minimizing the rank of X. In both cases, this amounts to a convex program, which can be solved efficiently.

(We mention that at this point we do not require that our density operators are normalized to have unit trace. We will return to this point later in section 5.)

The first estimator is obtained by constrained trace-minimization (also known as the matrix Dantzig selector):

$$\begin{equation} \hat{\rho}_{\mathrm{DS}} = \arg\min_X \Vert{X}\Vert_{\mathrm{tr}} \ \ \mbox{s.t.} \ \Vert{\mathcal{A}^*(\mathcal{A}(X)-y)}\Vert \leqslant \lambda. \end{equation} \tag{ 3 }$$

Note that the constraint bounds the operator norm of $\mathcal {A}^*(\mathcal {A}(X)-y)$. The parameter λ should be set according to the amount of noise in the data y; we will discuss this later.

The second estimator is obtained by least-squares linear regression with trace-norm regularization (also known as the matrix Lasso):

$$\begin{equation} \hat{\rho}_{\mathrm{Lasso}} = \arg\min_X \textstyle\frac{1}{2} \Vert{\mathcal{A}(X)-y}\Vert_2^2 + \mu \Vert{X}\Vert_{\mathrm{tr}}. \end{equation} \tag{ 4 }$$

Again the regularization parameter μ should be set according to the noise level; we will discuss this later.

One additional point is that we do not require positivity of the output in our definition of the estimators (3) and (4). One can add this constraint (since it is convex) and the conclusions below remain unaltered. We will explicitly add positivity later on when we do numerical simulations, and discuss any trade-offs that result from this.

2.3. Error bounds

In previous works on compressed tomography [17, 18], error bounds were proved by constructing a 'dual certificate' [67] (using convex duality to characterize the solution to the trace-minimization convex program). Here we derive stronger bounds using a different tool, known as the RIP. The RIP for low-rank matrices was first introduced in [32], and was recently shown to hold for random Pauli measurements [30].

We say that the sampling operator $\mathcal {A}$ satisfies the rank-r RIP if there exists some constant 0 ⩽ δ_r < 1 such that, for all $X \in \mathbb {C}^{d\times d}$ with rank r,

$$\begin{equation} (1-\delta_r) \Vert{X}\Vert_{\rm F} \leqslant \Vert{\mathcal{A}(X)}\Vert_2 \leqslant (1+\delta_r) \Vert{X}\Vert_{\rm F}. \end{equation} \tag{ 5 }$$

For our purposes, we can assume that X is Hermitian. Note that this notion of RIP is analogous to that used in [19], except that it holds over the set of low-rank matrices, rather than the set of sparse matrices.

The random Pauli sampling operator (1) satisfies RIP with high probability, provided that m ⩾ Crd log⁶ d (for some absolute constant C); this was recently shown in [30]. We note, however, that this RIP result requires m to be larger by a poly log d factor compared to previous results based on dual certificates. Although m is slightly larger, the advantage is that when combined with the results of [31], this immediately implies strong error bounds for the matrix Dantzig selector and the matrix Lasso.

To state these improved error bounds precisely, we introduce some definitions. For the rest of section 2, let C, C₀, C₁, C'₀ and C'₁ be fixed absolute constants. For any quantum state ρ, we write ρ = ρ_r + ρ_c, where ρ_r is the best rank-r approximation to ρ (consisting of the largest r eigenvalues and eigenvectors), and ρ_c is the residual part. Now we have the following:

Theorem 1. Let $\mathcal {A}$ be the random Pauli sampling operator (1) with m ⩾ Crd log⁶ d. Then, with high probability, the following holds:

Let $\hat {\rho }_{\mathrm {DS}}$ be the matrix Dantzig selector (3), and choose λ so that $\Vert {\mathcal {A}^*(z)}\Vert \leqslant \lambda$. Then

$$\begin{eqnarray*} &&\Vert{\hat{\rho}_{\mathrm{DS}} - \rho}\Vert_{\mathrm{tr}} \leqslant C_0 r\lambda + C_1 \Vert{\rho_c}\Vert_{\mathrm{tr}}. \end{eqnarray*}$$

Alternatively, let $\hat {\rho }_{\mathrm {Lasso}}$ be the matrix Lasso (4), and choose μ so that $\Vert {\mathcal {A}^*(z)}\Vert \leqslant \mu /2$. Then

$$\begin{eqnarray*} &&\Vert{\hat{\rho}_{\mathrm{Lasso}} - \rho}\Vert_{\mathrm{tr}} \leqslant C'_0 r\mu + C'_1 \Vert{\rho_c}\Vert_{\mathrm{tr}}. \end{eqnarray*}$$

In these error bounds, the first term depends on the statistical noise z. This, in turn, depends on the number of copies of the state that are available in the experiment; we will discuss this in the next section. The second term is the rank-r approximation error. It is clearly optimal, up to a constant factor.

Proof. These error bounds follow from the RIP as shown by Theorem 2.1 in [30], and a straightforward modification of Lemma 3.2 in [31] to bound the error in trace norm rather than Frobenius norm (this is similar to the proof of Theorem 5 in [66]).

The modification of Lemma 3.2 in [31] is as follows⁶. (For the remainder of this section, equation numbers of the form (III.x) refer to [31].) In the case of the Dantzig selector, let $H = \hat {\rho }_{\mathrm {DS}} - \rho$ be the difference of the matrix Dantzig selector and the quantum state. Following equation (III.8) and decomposing H = H₀ + H_c exactly as in [31], we can obtain the following bound:

$$\begin{equation} \Vert{H}\Vert_{\mathrm{tr}} \leqslant \Vert{H_0}\Vert_{\mathrm{tr}} + \Vert{H_c}\Vert_{\mathrm{tr}} \leqslant 2\Vert{H_0}\Vert_{\mathrm{tr}} + 2\Vert{\rho_c}\Vert_{\mathrm{tr}}, \end{equation} \tag{ 6 }$$

where we used the triangle inequality and equation (III.8). Then, at the end of the proof, we write

$$\begin{eqnarray} \Vert{H_0}\Vert_{\mathrm{tr}} &\leqslant& \sqrt{2r} \Vert{H_0}\Vert_{\rm F} \leqslant \sqrt{2r} \Vert{H_0+H_1}\Vert_{\rm F} \nonumber\\ & \leqslant& C_1 4\sqrt{2} r \lambda + C_1 2\sqrt{2} \delta_{4r} \Vert{\rho_c}\Vert_{\mathrm{tr}}, \end{eqnarray} \tag{ 7 }$$

where we used Cauchy–Schwarz, the fact that H₀ and H₁ are orthogonal, the bound on $\Vert {H_0+H_1}\Vert _{\mathrm { F}}$ following equation (III.13) and equation (III.7). Combining (6) and (7) gives our desired error bound. The error bound for the Lasso is obtained in a similar way; see section III.G in [31].   □

2.4. Sample complexity

Here we bound the sample complexity of our tomography scheme; that is, we bound the number of copies of the unknown quantum state ρ that are needed to obtain our estimate up to some accuracy. What we show, roughly speaking, is that $t = O\bigl ((\frac {rd}{\varepsilon })^2 \log d\bigr )$ copies are sufficient to reconstruct an estimate of an unknown rank r state up to accuracy ε in the trace distance. For comparison, note that when r = d, and one does full tomography, O(d⁴/ε²) copies are sufficient to estimate a full-rank state with accuracy ε in trace distance⁷.

To make this claim precise, we need to specify how we construct our data vector y from the measurement outcomes on the t copies of the state ρ. For the matrix Dantzig selector, suppose that

$$\begin{equation} t \geqslant 2 C_4 (C_0 r/\varepsilon)^2 d(d+1) \log d \end{equation} \tag{ 8 }$$

for some constants C₄ > 1 and ε ⩽ C₀. (For the matrix Lasso, substitute C'₀ for C₀ in these equations.) We construct an estimate of $\mathcal {A}(\rho )$ as follows: for each i∈[m], we take t/m copies of ρ, measure the random Pauli observable P_i on each of the copies and use this to estimate Tr(P_iρ). Then let y be the resulting estimate of $\mathcal {A}(\rho )$, and let $z = y - \mathcal {A}(\rho )$. Everything else is defined exactly as in theorem 1.

Theorem 2. Given $t = O\bigl ((\frac {rd}{\varepsilon })^2 \log d\bigr )$ copies of ρ as in equation (8) and measured as discussed above, the following holds with high probability over the measurement outcomes.

Let $\hat {\rho }_{\mathrm {DS}}$ be the matrix Dantzig selector (3), and set λ = ε/(C₀r) for some ε > 0. Then

$$\begin{eqnarray*} &&\Vert{\hat{\rho}_{\mathrm{DS}} - \rho}\Vert_{\mathrm{tr}} \leqslant \varepsilon + C_1 \Vert{\rho_c}\Vert_{\mathrm{tr}}. \end{eqnarray*}$$

Alternatively, let $\hat {\rho }_{\mathrm {Lasso}}$ be the matrix Lasso (4), and set μ = ε/(C'₀r). Then

$$\begin{eqnarray*} &&\Vert{\hat{\rho}_{\mathrm{Lasso}} - \rho}\Vert_{\mathrm{tr}} \leqslant \varepsilon + C'_1 \Vert{\rho_c}\Vert_{\mathrm{tr}}. \end{eqnarray*}$$

Proof. Our claim reduces to the following question: if we fix some value of λ > 0, how many copies of ρ are needed to ensure that the measurement data y satisfy $\Vert {\mathcal {A}^*(y-\mathcal {A}(\rho ))\Vert } \leqslant \lambda$? Then one can apply theorem 1 to obtain an error bound for our estimate of ρ.

Let t be the number of copies of ρ. Say we fix the measurement operator $\mathcal {A}$, i.e. we fix the choice of the Pauli observables P₁,...,P_m. (The measurement outcomes are still random, however.) For i∈[m] and j∈[t/m], let B_ij∈{1,−1} be the outcome of the jth repetition of the experiment that measures the ith Pauli observable P_i. Note that $\mathbbm {E} B_{ij} = \mathrm {Tr}(P_i\rho )$. Then construct the vector $y \in \mathbb {R}^m$ containing the estimated expectation values (scaled by $\sqrt {d/m}$):

$$\begin{equation} y_i = \sqrt{\frac{d}{m}} \cdot \frac{m}{t} \sum_{j=1}^{t/m} B_{ij}, \quad i \in [m]. \end{equation} \tag{ 9 }$$

Note that $\mathbbm {E} y = \mathcal {A}(\rho )$.

We will bound the deviation $\Vert {\mathcal {A}^*(y-\mathcal {A}(\rho ))\Vert }$, using the matrix Bernstein inequality. First, we write

$$\begin{equation} \mathcal{A}^*(y) = \sqrt{\frac{d}{m}} \sum_{i=1}^m P_i y_i = \frac{d}{t} \sum_{i=1}^m \sum_{j=1}^{t/m} P_i B_{ij}, \end{equation} \tag{ 10 }$$

and also

$$\begin{equation} \mathcal{A}^*\mathcal{A}(\rho) = \frac{d}{m} \sum_{i=1}^m P_i \mathrm{Tr}(P_i\rho). \end{equation} \tag{ 11 }$$

We can now write $\mathcal {A}^*(y-\mathcal {A}(\rho ))$ as a sum of independent (but not identical) matrix-valued random variables:

$$\begin{equation} \mathcal{A}^*(y-\mathcal{A}(\rho)) = \sum_{i=1}^m \sum_{j=1}^{t/m} X_{ij}, \quad X_{ij} = \frac{d}{t} P_i \bigl[B_{ij} - \mathrm{Tr}(P_i\rho)\bigr]. \end{equation} \tag{ 12 }$$

Note that $\mathbbm {E} X_{ij} = 0$ and $\Vert {X_{ij}}\Vert \leqslant 2d/t =: R$. Also, for the second moment we have

$$\begin{eqnarray} &&\mathbbm{E}\bigl(X_{ij}^2\bigr) = \mathbbm{E}\left(\frac{d^2}{t^2} I \bigl[B_{ij} - \mathrm{Tr}(P_i\rho)\bigr]^2\right)= \frac{d^2}{t^2} I \bigl[1 - \mathrm{Tr}(P_i\rho)^2\bigr]. \end{eqnarray} \tag{ 13 }$$

Then we have

$$\begin{eqnarray} &&\sigma^2 := \Bigl\Vert \sum_{ij} \mathbbm{E}(X_{ij}^2) \Bigr\Vert = \sum_{ij} \frac{d^2}{t^2} \bigl[1 - \mathrm{Tr}(P_i\rho)^2\bigr] \leqslant t \cdot \frac{d^2}{t^2} = \frac{d^2}{t}. \end{eqnarray} \tag{ 14 }$$

Now the matrix Bernstein inequality (Theorem 1.4 in [68]) implies that

$$\begin{eqnarray} &&\fl \Pr\bigl[\Vert{\mathcal{A}^*(y-\mathcal{A}(\rho))}\Vert \geqslant \lambda\bigr] \leqslant d\cdot \exp\left( -\frac{\lambda^2/2}{\sigma^2 + (R\lambda/3)}\right) \leqslant d\cdot \exp\left( -\frac{t\lambda^2/2}{d(d+1)} \right) \end{eqnarray} \tag{ 15 }$$

(where we assumed that λ ⩽ 1).

For the matrix Dantzig selector, we set λ = ε/(C₀r), and we get that, for any t ⩾ 2C₄λ⁻²d(d + 1)log d = 2C₄(C₀r/ε)²d(d + 1)log d,

$$\begin{eqnarray} &&\Pr\left[\Vert{\mathcal{A}^*(y-\mathcal{A}(\rho))}\Vert \geqslant \frac{\varepsilon}{C_0 r}\right] \leqslant d\cdot \exp(-C_4 \log d)= d^{1-C_4}, \end{eqnarray} \tag{ 16 }$$

which is exponentially small in C₄. Plugging into theorem 2 completes the proof of our claim. A similar argument works for the matrix Lasso.   □

How good are the sample complexities of our algorithms? Here we go a long way toward answering this question by proving nearly tight lower bounds on the sample complexity for generic rank r quantum states using single-copy Pauli measurements. Previous work on single-copy lower bounds has only treated the case of pure states [33].

Roughly speaking, we show the following, which we will make precise at the end of the section.

Theorem 3 The number of copies t needs to grow at least as fast as $\Omega \bigl (r^2d^2/\log d\bigr )$ to guarantee a constant trace-norm confidence interval for all rank-r states.

The argument proceeds in three steps. First, we fix our notion of risk to be the minimax risk, meaning we seek to minimize our worst-case error according to some error metric such as the trace distance. We want to know how many copies of the unknown state we need to make this minimax risk an arbitrarily small constant. For a fixed set of two-outcome measurements, say Pauli measurements, we then show that some states require many copies to achieve this. In particular, these states have the property that they are globally distinguishable (their trace distance is bounded from below by a constant) but their (Pauli) measurement statistics look approximately the same (each measurement outcome is close to unbiased). The more such states there are, the more copies we need to distinguish between them all using solely Pauli measurements. Finally, we use a randomized argument to show that, in fact, there are exponentially many such states. This yields the desired lower bound on the sample complexity.

Let Σ be some set of density operators. We want to put lower bounds on the performance of any estimation protocol for states in Σ. (We do not initially restrict ourselves to states with low rank.)

We assume that the protocol has access to t copies of an unknown state ρ∈Σ on which it makes measurements one by one. At the ith step, it has to decide which observable to measure. Let us restrict ourselves to two-outcome measurements, which can be described by positive operator-valued measures (POVMs) $\{\Pi _i, \mathbbm {1}-\Pi _i\}$, where each Π_i satisfies $0\leqslant \Pi _i \leqslant \mathbbm {1}$ and may be chosen from a set $\mathcal {P}$. (We do not initially restrict ourselves to Pauli measurements, either.) We allow the choice of the ith observable to depend on the previous outcomes. We refer to the random variables which jointly describe the choice of the ith measurement and its random outcome as Y_i. At the end, these are mapped to an estimate $\hat \rho (Y_1, \dots , Y_t)\in \Sigma$.

In other words, an estimation protocol is specified by a set of functions

$$\begin{eqnarray*} &&\Pi_i : Y_1, \dots, Y_{i-1} \mapsto \mathcal{P}, \\ &&\hat\rho : Y_1, \dots, Y_t \mapsto \Sigma. \end{eqnarray*}$$

Consider a distance measure $\Delta :\Sigma \times \Sigma \to \mathbb {R}$ on the states in Σ. (For example, this could be the trace distance or the infidelity; we need not specify.) Suppose that we are given t copies of the unknown state ρ, and the maximum deviation we wish to tolerate between our estimate $\hat {\rho }$ and the true state ρ is given by  > 0. Now define the minimax risk

$$\begin{equation} M^*(t,\epsilon) = \inf_{\langle \hat \rho, \Pi_i \rangle} \sup_{\rho\in\Sigma} \Pr\bigl[ \Delta(\hat\rho(Y), \rho) > \epsilon\bigr], \end{equation} \tag{ 17 }$$

where the infimum is over all estimation procedures $\langle \hat \rho , \Pi _i\rangle$ on t copies with estimator $\hat \rho$ and choice of observables given by Π_i. That is, we are considering the 'best' protocol to be the one whose worst-case probability of deviation is the least.

The next lemma shows that if there are a large number of states in Σ which are far apart (by at least ) and whose statistics look nearly random for all measurements in $\mathcal {P}$, then the minimax risk is lower-bounded by a function that decreases linearly with the number of copies t. Hence, to achieve a small minimax risk, the number of copies t must be large.

Lemma 1. Assume that there are states $\rho _1, \dots , \rho _s\subset \Sigma$ such that

$$\begin{eqnarray*} &\forall i\neq j: \Delta(\rho_i, \rho_j) \geqslant \epsilon ,\\ &\forall i,\forall\, \Pi \in \mathcal{P}: \bigl|\mathrm{tr} \rho_i \Pi - \textstyle\frac12\bigr| \leqslant \alpha . \end{eqnarray*}$$

Then the minimax risk as defined in (17) fulfills

$$\begin{eqnarray*} &&M^*(t,\epsilon)>1-\frac{4\alpha^2 t + 1}{\log s}. \end{eqnarray*}$$

Proof. Let X be a random variable taking values uniformly in [s]. Let $Y_1, \dots , Y_t$ be random variables, Y_i describing the outcome of the ith measurement carried out on ρ_X. Define

$$\begin{eqnarray*} &&\hat X(Y) := \arg \min_i \Delta(\hat\rho(Y),\rho_i). \end{eqnarray*}$$

Then

$$\begin{equation} \Pr\bigl[ \Delta(\hat\rho(Y),\rho_i) > \epsilon \bigr] \geqslant \Pr\bigl[ \hat X(Y) \neq X \bigr] . \end{equation} \tag{ 18 }$$

Now combine Fano's inequality

$$\begin{eqnarray*} &&H(X|\hat X)\leqslant 1+ \Pr[ \hat X \neq X ] \log s, \end{eqnarray*}$$

the data processing inequality

$$\begin{eqnarray*} &&I(X;\hat X(Y))\leqslant I(X;Y), \end{eqnarray*}$$

in terms of the mutual information I(X;Y ): = H(X) − H(X|Y ), and the fact that H(X) = log s to obtain

$$\begin{eqnarray*} \Pr[ \hat X(Y) \neq X ] &\geqslant \frac{H(X|\hat X) - 1}{\log s} \\ &= \frac{H(X)-I(X;\hat X) - 1}{\log s} \\ &= 1-\frac{I(X;\hat X) + 1}{\log s} \\ &\geqslant 1-\frac{I(X;Y) + 1}{\log s} \\ &= 1-\frac{H(Y) - H(Y|X) + 1}{\log s} \\ &\geqslant 1-\frac{t - \frac{1}{s} \sum_{i=1}^s H(Y|X=i) + 1}{\log s}. \end{eqnarray*}$$

Let h(p) be the binary entropy and recall the standard estimate

$$\begin{eqnarray*} &&h(1/2\pm\alpha) \geqslant (1-4\alpha^2). \end{eqnarray*}$$

Combine that with the chain rule [69, Theorem 2.5.1]

$$\begin{eqnarray*} H(Y|X=i) &= \sum_{j=1}^t H(Y_j | Y_{j-1}, \dots, Y_1, X=i) \\ &\geqslant t (1-4\alpha^2). \end{eqnarray*}$$

The advertised bound follows.   □

We now show how to construct a set of states, each with rank r, that satisfy the conditions of lemma 1 in terms of the trace distance and the set of Pauli measurements. The following lemma shows how, given such a set of states, we can enlarge it by one. We will then apply this lemma repeatedly, to get a total of s < e^crd states.

Lemma 2. For any $0< \epsilon < 1-\frac {r}{d}$, let $\rho _1, \dots , \rho _s$ be normalized⁸ rank-r projections on $\mathbb {C}^d$, where s < e^c()rd and c() is specified below. Then there exists a normalized rank-r projection ρ such that

$$\begin{equation} \forall i \in [s]: \textstyle \frac{1}{2}\|\rho - \rho_i\|_{\mathrm{tr}} \geqslant \epsilon, \end{equation} \tag{ 19 }$$

$$\begin{equation} \forall P_k \neq \mathbbm{1}: \bigl|\mathrm{tr}\bigl[\textstyle \frac{1}{2}(\mathbbm{1}\pm P_k) \rho\bigr]-\textstyle \frac{1}{2}\bigr| \leqslant \alpha. \end{equation} \tag{ 20 }$$

Here, $\alpha ^2 = O\bigr (\frac {\log d}{rd}\bigr )$, the P_k are n-qubit Pauli operators, and

$$\begin{eqnarray*} &&c(\epsilon) = \frac{\ln(8/\pi)}{2 r d}+ \frac{1}{32}\Bigl[\left(1-\frac{r}{d}\right)-\epsilon\Bigr]^2. \end{eqnarray*}$$

Proof. Let ρ₀ be some normalized rank-r projection and choose ρ according to⁹

$$\begin{eqnarray*} &&\rho = O \rho_0 O^T \end{eqnarray*}$$

for a Haar-random O∈ SO(d). Here we use the special orthogonal group SO(d) because the analysis becomes marginally simpler than if we use a unitary group.

To check (19), choose i∈[s] and define R_i to be the projector onto the range of ρ_i. Also define the function

$$\begin{eqnarray*} &&f: O \mapsto \|\rho_i - O \rho_0 O^T \|_{\mathrm{tr}}. \end{eqnarray*}$$

We can bound the magnitude of f using the pinching inequality:

$$\begin{eqnarray*} &&\fl f(O) \geqslant \|\rho_i \!-\! R_i O \rho_0 O^T R_i\|_{\mathrm{tr}} + \|R_i^\perp O \rho_0 O^T R_i^\perp\|_{\mathrm{tr}} \geqslant \mathrm{tr}(\rho_i) \!-\!\mathrm{tr}(R_i O \rho_0 O^T R_i) + \mathrm{tr}(R_i^\perp O \rho_0 O^T R_i^\perp) \\ &&\fl \hphantom{f(O)} = 1 +\mathrm{tr}\bigl[(R_i^\perp - R_i)O \rho_0 O^T\bigr] . \end{eqnarray*}$$

From this we can bound the expectation value of f over the special orthogonal group:

$$\begin{eqnarray*} \mathbbm{E}\bigl[f(O)\bigr] &\geqslant& 1+ \mathrm{tr}\bigl[(R_i^\perp - R_i) \big( \mathbbm{E} O \rho O^T \big)\bigr]\\ &=& 1+ \frac{1}{d} \mathrm{tr}\bigl[(R_i^\perp - R_i) \mathbbm{1}\bigr] \\ &=& 1 + \frac{d-2r}{d}=2\Bigl(1-\frac{r}{d}\Bigr). \end{eqnarray*}$$

Next we get an upper bound of $\frac {4}{\sqrt {r}}$ on the Lipschitz constant of f with respect to the Frobenius norm:

$$\begin{eqnarray*} |f(O+\Delta) - f(O)| &\leqslant& \|(O+\Delta) \rho_0 (O+\Delta)^T - O \rho_0 O^T\|_{\mathrm{tr}} \\ &\leqslant& \|O \rho_0 \Delta^T\|_{\mathrm{tr}} + \|\Delta\rho_0 O^T\|_{\mathrm{tr}} + \|\Delta \rho_0 \Delta^T\|_{\mathrm{tr}}\\ &=&2\|\Delta \rho_0\|_{\mathrm{tr}} +\mathrm{tr}(\Delta \rho_0 \Delta^T) \\ &\leqslant& 2\sqrt r \|\Delta \rho_0\|_{\rm F} + \mathrm{tr}(\rho_0 \Delta \Delta^T)\\ & \leqslant& 2\sqrt r \|\Delta\|_{\rm F} \|\rho_0\| + \frac{1}{r} \|\Delta\| \|\Delta^T\|_{\mathrm{tr}}\\ &\leqslant &\frac{2}{\sqrt r} \|\Delta\|_{\rm F} + \frac{2}{r} \sqrt{r} \|\Delta\|_{\rm F} \leqslant \frac{4}{\sqrt r} \|\Delta\|_{\rm F} , \end{eqnarray*}$$

where we use ∥Δ∥ ⩽ 2 in the last line, which follows from the triangle inequality and the fact that any Δ can be written as a difference Δ = O' − O for O'∈ SO(d).

From these ingredients, we can invoke Lévy's lemma on the special orthogonal group [70, Theorem 6.5.1] to obtain that, for all t > 0,

$$\begin{eqnarray*} &&\Pr\bigl[ \| \rho_i - \rho \|_{\mathrm{tr}} < 2(1-r/d)- t \bigr] \leqslant \mathrm{e}^{c_1} \exp\Bigl(- \frac{c_2 t^2 rd}{16}\Bigr), \end{eqnarray*}$$

where the constants are given by $c_1 = \ln \bigl (\sqrt {\pi /8}\bigr )$ and c₂ = 1/8. Now choose t = 2(1 − r/d) − 2 and apply the union bound to obtain

$$\begin{eqnarray*} &&\fl \Pr\bigl[(19) \textrm{ does not hold}\bigr] <\, \mathrm{e}^{c(\epsilon) rd} \Pr\bigl[ \| \rho_i - \rho \|_{\mathrm{tr}} < 2(1-r/d)- t \bigr]\\ &&\fl \hphantom{\Pr\bigl[(19) \textrm{ does not hold}\bigr]}\leqslant \exp\mathclose{}\bigg(rd\Bigl[c(\epsilon)-\frac{c_2 t^2 }{16}\Big]+ c_1\bigg)=1. \end{eqnarray*}$$

The upper bound on follows from the requirement that t > 0. This shows that ρ indeed satisfies equation (19).

Now we move on to equation (20). For any non-identity Pauli matrix P_k, define a function

$$\begin{eqnarray*} &&f: O \mapsto \mathrm{tr}(P_k O \rho_0 O^T) . \end{eqnarray*}$$

Clearly, we have $\mathbbm {E}[f(O)] = 0$. We again wish to bound the rate of change so that we can use Lévy's lemma, so we compute

$$\begin{eqnarray*} &&(\mathrm{d}_O f)(\Delta) = \mathrm{tr}\bigl(\rho_0 O^T P_k \Delta\bigr) + \mathrm{tr}\bigl(P_k O \rho_0 \Delta^T\bigr) = \mathrm{tr}\bigl[(\rho_0 O^T P_k + \rho_0^T O^T P_k^T) \Delta \bigr], \end{eqnarray*}$$

which implies that

$$\begin{eqnarray*} &&\|\nabla f(O)\|_2 = \|\rho_0 O^T P_k + \rho_0^T O^T P_k^T \|_{\rm F} \leqslant \frac{2}{\sqrt{r}} . \end{eqnarray*}$$

Lévy's lemma then gives for all t > 0

$$\begin{eqnarray*} &&\Pr\bigl[ \vert{ \mathrm{tr} P_k \rho }\vert > t \bigr] \leqslant \mathrm{e}^{c_1} \exp\biggl(- \frac{c_2 t^2 rd}{4}\biggr). \end{eqnarray*}$$

Choosing t = 2α, and α² = 4 ln(d⁴π/8)/(rd), then the union bound gives us

$$\begin{eqnarray*} \fl \Pr\bigl[ (20) \textrm{ does not hold}\bigr] < &\, d^2 \Pr\bigl[ |\mathrm{tr} P_k \rho| > 2\alpha\bigr] \leqslant \exp\Bigl(2\ln d -\frac{c_2(2\alpha)^2rd}{4}+ c_1\Bigr) = 1, \end{eqnarray*}$$

from which the lemma follows.   □

We remark that a version of lemma 2 continues to hold even if we can adaptively choose from as many as 2^O(n) additional measurements which are globally unitarily equivalent to Pauli measurements.

We now combine the two previous lemmas, to get a precise version of theorem 3. (To summarize: we use lemma 2 repeatedly to construct a large set of states that are difficult to distinguish using Pauli measurements; then we use lemma 1 to lower-bound the minimax risk for estimating an unknown state using t copies; hence, if an estimation procedure has a small minimax risk, t must be large.)

Theorem 4 Fix $\epsilon \in (0, 1-\frac {r}{d})$ and δ∈[0,1). Then, in order for an estimation procedure to satisfy M*(t,) ⩽ δ, the number of copies t must grow like

$$\begin{eqnarray*} &&t = \Omega\biggl(\frac{r^2 d^2}{\log d}\biggr), \end{eqnarray*}$$

where the implicit constant depends on δ and .

Here we sketch how the technique of DFE, introduced in [33, 34] for pure states, can be used to estimate the fidelity between a low-rank estimate $\hat \rho$ and the true state ρ. The only assumption we make is that $\hat \rho$ is a positive semidefinite matrix with $\mathrm {Tr}(\hat \rho ) \leqslant 1$ and $\mathrm {rank}(\hat \rho ) = r$. No assumption at all is needed on ρ. In fact, we also do not assume that we obtained the estimate $\hat \rho$ from any of the estimators which were discussed previously. Our certification procedure works regardless of how one obtains $\hat \rho$, and so it even applies to the situation where $\hat \rho$ was chosen from a subset of variational ansatz states, as in [61].

Recall that the main idea of DFE is to take a known pure state |ψ〉〈ψ| and from it define a probability distribution $\Pr (i)$ such that by estimating the Pauli expectation values Tr(ρP_i) and suitably averaging it over $\Pr (i)$ we can learn an estimate of 〈ψ|ρ|ψ〉. In fact, one does not need to do a full average; simply sampling from the distribution a few times is sufficient to produce a good estimate.

For non-Hermitian rank-1 matrices |ϕ_j〉〈ϕ_k| instead of pure states, we need a very slight modification of DFE. Following the notation in [33], we simply redefine the probability distribution as $\Pr (i) = \vert {\langle {\phi _j}|P_i|{\phi _k}\rangle }\vert ^2/d$ and the random variable X = Tr(P_iρ)/〈ϕ_k|P_i|ϕ_j〉. It is easy to check that $\mathbb {E}(X) = \langle {\phi _j}|\rho |{\phi _k}\rangle$ and the variance of X is at most one. Then all of the conclusions in [33, 34] hold for estimating the quantity 〈ϕ_j|ρ|ϕ_k〉. In particular, we can obtain an estimate to within an error ± with probability 1 − δ by using only single-copy Pauli measurements and $O(\frac {1}{\epsilon ^2\delta } + \frac {d \log (1/\delta )}{\epsilon ^2})$ copies of ρ.

Our next result shows that by obtaining several such estimates using DFE, we can also infer an estimate of the mixed state fidelity between an unknown state ρ and a low-rank estimate $\hat \rho$, given by

$$\begin{equation} F(\rho,\hat\rho) = \bigl[\mathrm{Tr} \sqrt{G}\bigr]^2, \end{equation} \tag{ 21 }$$

where for brevity we define $G = \sqrt {\hat \rho }\rho \sqrt {\hat \rho }$. (Note that some authors use the square root of this quantity as the fidelity. Our definition matches [33].) The asymptotic cost in sample complexity is far less than the cost of initially obtaining the estimate $\hat \rho$ when r is sufficiently small compared to d.

Theorem 5. Given a state estimate $\hat \rho$ with $\mathrm {rank}(\hat \rho ) = r$, the number of copies t of the state ρ required for estimating $F(\rho ,\hat \rho )$ to within ±ε with probability 1 − δ using single-copy Pauli measurements satisfies

$$\begin{equation} t=O\left(\frac{r^5}{\varepsilon^4}\bigr[d \log(r^2/\delta) + r^2/\delta\bigl]\right) . \end{equation} \tag{ 22 }$$

Proof. The result uses the DFE protocol of [33, 34], modified as mentioned above, where the states |ϕ_j〉 are the eigenstates of $\hat \rho$.

Expand $\hat \rho$ in its eigenbasis as $\hat \rho = \sum _{j=1}^r \lambda _j |{\phi _j}\rangle \!\langle {\phi _j}|$. Then we have

$$\begin{equation} G = \sum_{j,k=1}^r \sqrt{\lambda_j \lambda_k} \langle{\phi_j}| \rho |{\phi_k}\rangle|{\phi_j}\rangle\!\langle{\phi_k}| . \end{equation} \tag{ 23 }$$

For all 1 ⩽ j ⩽ k ⩽ r, we use DFE to obtain an estimate $\hat {g}_{jk}$ of the matrix element 〈ϕ_j|ρ|ϕ_k〉, up to some additive error _jk that is bounded by a constant, $\vert {\epsilon _{jk}}\vert \leqslant \epsilon _0$.

If each estimate is accurate with probability 1 − 2δ/(r² + r), then by the union bound the probability that they are all accurate is at least 1 − δ. The total number of copies t required for this is

$$\begin{equation} t=O\left(\frac{r^2}{\epsilon_0^2}\bigr[d \log(r^2/\delta) + r^2/\delta\bigl]\right) . \end{equation} \tag{ 24 }$$

Let $\hat {g}_{jk} = \hat {g}_{kj}^*$, and let

$$\begin{equation} \hat{G} = \sum_{j,k=1}^r \sqrt{\lambda_j \lambda_k} \hat{g}_{jk} |{\phi_j}\rangle\!\langle{\phi_k}| \end{equation} \tag{ 25 }$$

be our estimate for G. Finally, let $\hat {G}^+ = [\hat {G}]_+$ be the positive part of the Hermitian matrix $\hat {G}$, and let $\hat {F} = \bigl [\mathrm {Tr}\sqrt {\hat {G}^+}\bigr ]^2$ be our estimate of the fidelity $F(\rho ,\hat \rho )$. Note that we may assume that $\hat F \leqslant 1$, since if it were larger, we can only improve our estimate by just truncating it back to 1.

We now bound the error of this fidelity estimate. We can write $\hat {G}$ as a perturbation of G, $\hat {G} = G+E$, where

$$\begin{equation} E = \sum_{j,k=1}^r \sqrt{\lambda_j \lambda_k} \epsilon_{jk} |{\phi_j}\rangle\!\langle{\phi_k}| . \end{equation} \tag{ 26 }$$

First note that the Frobenius norm of this perturbation is small,

$$\begin{equation} \Vert{E}\Vert_{\rm F} = \biggl(\sum_{j,k} \lambda_j \lambda_k \vert{\epsilon_{jk}}\vert^2 \biggr)^{1/2} \leqslant \epsilon_0 \biggl\vert\sum_j \lambda_j \biggr\vert = \epsilon_0 . \end{equation} \tag{ 27 }$$

(If $\hat \rho$ is subnormalized, then this last equality becomes an inequality.)

Next observe that

$$\begin{eqnarray*} &&\vert{F-\hat F}\vert = \Bigl\vert \mathrm{Tr}(\sqrt{G} + \sqrt{\hat{G}^+}) \mathrm{Tr}(\sqrt{G} - \sqrt{\hat{G}^+}) \Bigr\vert \leqslant 2 \Delta , \end{eqnarray*}$$

where we define

$$\begin{equation} \Delta = \bigl\vert\mathrm{Tr}(\sqrt{G} - \sqrt{\hat{G}^+})\bigr\vert . \end{equation} \noindent \tag{ 28 }$$

Using the reverse triangle inequality, we can bound Δ in terms of the trace norm

$$\begin{equation} \Delta \leqslant \left\Vert\sqrt{G} - \sqrt{\hat{G}^+}\right\Vert_{\mathrm{tr}} . \end{equation} \noindent \tag{ 29 }$$

Using [71, Theorem X.1.3] in the first step, we find that

$$\begin{eqnarray*} &&\left\Vert\sqrt{G} - \sqrt{\hat{G}^+}\right\Vert_{\mathrm{tr}} \leqslant \Bigl\Vert\sqrt{\vert{G -\hat{G}^+}\vert}\,\Bigr\Vert_{\mathrm{tr}} \leqslant \sqrt{r}\sqrt{\bigl\Vert G -\hat{G}^+\bigr\Vert_{\mathrm{tr}}}, \end{eqnarray*} \noindent$$

where the second bound follows from the Cauchy–Schwarz inequality on the vector of eigenvalues of $\vert {G -\hat {G}^+}\vert$. Using the Jordan decomposition of a Hermitian matrix into a difference of positive matrices X = [X]₊ − [X]₋, we can rewrite $G-\hat {G}^+$ as

$$\begin{equation} G-\hat{G}^+ = G-[G+E]_+ = -E - [G+E]_- . \end{equation} \noindent \tag{ 30 }$$

Then by the triangle inequality and positivity of G,

$$\begin{equation} \Vert{G-\hat{G}^+}\Vert_{\mathrm{tr}} \leqslant \Vert{E}\Vert_{\mathrm{tr}} + \Vert{[G+E]_-}\Vert_{\mathrm{tr}} \leqslant 2 \Vert{E}\Vert_{\mathrm{tr}}. \end{equation} \noindent \tag{ 31 }$$

Using the standard estimate $\Vert {E}\Vert _{\mathrm {tr}} \leqslant \sqrt {r}\Vert {E}\Vert _{\mathrm { F}}$, we find that

$$\begin{equation} \Delta \leqslant r^{3/4} \sqrt{2 \epsilon_0} . \end{equation} \noindent \tag{ 32 }$$

This gives our desired error bound, albeit in terms of ₀ instead of the final quantity ε. To get our total error to vanish, we take $\varepsilon = 2 r^{3/4} \sqrt {2 \epsilon _0}$, which gives us the final scaling in the sample complexity.   □

As a final remark, we note that by computing Δ for the special case $\hat \rho = \rho = \mathbbm {1}/r$, $E = \epsilon _0 \mathbbm{1}/r$, we find that

$$\begin{equation} \Delta = \sqrt{1+r \epsilon_0}-1 , \end{equation} \tag{ 33 }$$

so the error bound for this protocol is tight with respect to the scaling in ₀ (and hence ε). However, we cannot rule out that there are other protocols that achieve a better scaling. Also, it seems that the upper bound for the current protocol with respect to r could potentially be improved by a factor of r with a more careful analysis.

We numerically simulated the reconstruction of a quantum state given Pauli measurements and the estimators from equations (3) and (4). Here we compare these estimates to those obtained from a standard maximum likelihood estimate (MLE) as a benchmark.

As mentioned previously, all of our estimators have the advantage that they are convex functions with convex constraints, so the powerful methods of convex programming [72] guarantee that the exact solution can be found efficiently and certified. We take full advantage of this certifiability and do simulations which can be certified by interior-point algorithms. This way we can separate out the performance of the estimators themselves from the (sometimes heuristic) methods used to compute them on very large scale instances.

For our estimators, we will explicitly enforce the positivity condition. That is, we will simulate the following slight modifications of equations (3) and (4) given by

$$\begin{equation} \hat{\rho}_{\mathrm{DS}}' = \arg\min_{X\succeq0}\, \mathrm{Tr} X \ \ \mbox{s.t.} \ \Vert{\mathcal{A}^*(\mathcal{A}(X)-y)}\Vert \leqslant \lambda \end{equation} \tag{ 34 }$$

and

$$\begin{equation} \hat{\rho}_{\mathrm{Lasso}}' = \arg\min_{X\succeq0}\, \frac{1}{2} \Vert{\mathcal{A}(X)-y}\Vert_2^2 + \mu \mathrm{Tr} X . \end{equation} \tag{ 35 }$$

Moreover, whenever the trace of the resulting estimate is less than 1, we renormalize the state, .

5.1. Setting the estimator parameters λ and μ

5.2. Time needed to switch measurement settings

5.3. Other simulation parameters

We consider the following ensembles of quantum states. We initially choose n = 5 qubit states from the invariant Haar measure on $\mathbb {C}^{2^n}$. Then we add noise to the state by applying independent and identical depolarizing channels to each of the n qubits. Recall that the depolarizing channel with strength γ acts on a single qubit as

$$\begin{equation} \mathcal{D}_\gamma(\rho) = \gamma \frac{\mathbbm{1}}{2} + (1-\gamma) \rho; \end{equation}\noindent \tag{ 36 }$$

that is, with probability γ the qubit is replaced by a maximally mixed state and otherwise it is left alone. Our simulations assume very weak decoherence, and we use γ = 0.01.

We use two measures to quantify the quality of a state reconstruction $\hat \rho$ relative to the underlying true state ρ, namely the (squared) fidelity $\Vert \sqrt {\hat \rho }\sqrt {\rho }\Vert ^2_{\mathrm {tr}}$ and the trace distance $\frac {1}{2} \Vert {\rho - \hat \rho }\Vert _{\mathrm {tr}}$.

We used the interior-point solver SeDuMi [26] to do the optimizations in equations (34) and (35). Although these algorithms are not suitable for large numbers of qubits, they produce very accurate solutions, which allows for a more reliable comparison between the different estimators. The data and the code used to generate the data for these plots can be found with the source code for the arXiv version of this paper. For larger numbers of qubits, one may instead use specialized algorithms such as singular value thresholding (SVT) or fixed-point continuation (FPCA) to solve these convex programs [27–29]; however, the quality of the solutions depends somewhat on the algorithm, and it is an open problem to explore this in more detail.

For the MLE, we used the iterative algorithm of [73], which has previously been used in e.g. [74]. This method converged on every example we tried, so we did not have to use the more sophisticated 'diluted' method of [75] that guarantees convergence.

For the purposes of our simulations, we sampled our Pauli operators without replacement.

5.4. Results and analysis

The data are depicted in figure 1. The three subfigures (a)–(c) use three different values for the total sampling time T in increasing order. We have plotted both the fidelity and the trace distance (inset) between the recovered solution and the true state.

Zoom In Zoom Out Reset image size

Several features are immediately apparent. First, we see that both of the compressed sensing estimators consistently outperform MLE along a wide range of values of m, the number of Paulis sampled, even when we choose the optimal value of m just for the MLE. Once m is sufficiently large (but still m ≪ d²) all of the estimators converge to a reasonably high-fidelity estimate. Thus, it is not just the compressed sensing estimators which are capable of reconstructing nearly low-rank states with few measurement settings, but rather this seems to be a generic feature of many estimators. However, the compressed sensing estimators seem to be particularly well suited for the task.

When the total amount of time T is very small (figure 1(a)), then there is a large advantage of using compressed sensing. In this regime, if we insist on making an informationally complete measurement, there is barely enough time to make one measurement per setting, so the measurement results are dominated by statistical fluctuations, i.e. the signal-to-noise ratio is poor. Compressed sensing techniques—measuring an incomplete set of observables, and regularizing the solution to favor low-rank states—lead to much better results. However, note that in this situation, it is very difficult for any method, including compressed sensing, to achieve extremely high-fidelity state reconstructions (with fidelity greater than 95%, say).

As the total amount of time available increases, all of the estimators of course converge to higher fidelity estimates. Interestingly, for the compressed sensing estimators, there seems to be no trade-off whatsoever between the number of measurement settings and the fidelity, once T and m are sufficiently large. That is, the curve is completely flat as a function of m above a certain cutoff. This is most notable for the matrix Lasso, which consistently performs at least as well as the matrix Dantzig selector, and often better. These observations are consistent with the bounds proven earlier.

The flat curve as a function of m is especially interesting, because it suggests that there is no real drawback to using small values of m. Smaller values of m are attractive from the computational point of view because the runtime of each reconstruction algorithm scales with m. This makes a strong case for adopting these estimators in the future, and at the very least more numerical studies are needed to investigate how well these estimators perform more generally.

We draw the following conclusion from these simulations. The best performance for a fixed value of T is given by the matrix Lasso estimator of equations (4) and (35) in the regime where m is nearly as small as possible. Here the fidelity is larger than the other estimators (if only by a small or negligible amount when T is large), but using a small value for m means smaller memory and processing time requirements when doing the state reconstruction. MLE consistently underperforms the compressed sensing estimators, but still seems to 'see' the low-rank nature of the underlying state and converges to a reasonable estimate even when m is small.

Compressed sensing techniques can also be applied to quantum process tomography. Here, our method has an advantage when the unknown quantum process has a small Kraus rank, meaning that it can be expressed using only a few Kraus operators. This occurs, for example, when the unknown process consists of unitary evolution (acting globally on the entire system) combined with local noise (acting on each qubit individually, or more generally, acting on small subsets of the qubits).

Consider a system of n qubits, with Hilbert space dimension d = 2ⁿ. Let $\mathcal {E}$ be a completely positive trace-preserving (CP) map from $\mathbb {C}^{d\times d}$ to $\mathbb {C}^{d\times d}$, and suppose that $\mathcal {E}$ has Kraus rank r. Using compressed sensing, we will show how to characterize $\mathcal {E}$ using m = O(rd²log d) settings. (For comparison, standard process tomography requires d⁴ settings, since $\mathcal {E}$ contains d⁴ independent parameters in general.) Furthermore, our compressed sensing method requires only the ability to prepare product eigenstates of Pauli operators and make Pauli measurements, and it does not require any ancilla qubits.

We remark that, except for the ancilla-assisted method, just the notion of 'measurement settings' for process tomography does not capture all of the complexity because of the need to have an input to the channel. Here we define one 'input setting' to be a basis of states from which the channel input should be sampled uniformly. Then the total number of settings m is the sum of the number of measurement settings (Paulis, in our case) and input settings. This definition justifies the claim that the number of settings for our compressed process tomography scheme is m = O(rd² log d).

The analysis here focuses entirely on the number of settings m. We forgo a detailed analysis of t, the sample complexity, and instead leave this open for future work.

Note that there is a related set of techniques for estimating an unknown process that is element-wise sparse with respect to some known, experimentally accessible basis [19]. These techniques are not directly comparable to ours, since they assume a greater amount of prior knowledge about the process, and they use measurements that depend on this prior knowledge. We will discuss this in more detail at the conclusion of this section.

6.1. Our method

First, recall the Jamiołkowski isomorphism [76]: the process $\mathcal {E}$ is completely and uniquely characterized by the state

$$\begin{eqnarray*} &&\rho_{\mathcal{E}} = (\mathcal{E}\otimes\mathcal{I})(|{\psi_0}\rangle\langle{\psi_0}|), \end{eqnarray*}$$

where $|{\psi _0}\rangle = \frac {1}{\sqrt {d}} \sum _{j=1}^{d} |{j}\rangle \otimes |{j}\rangle$. Note that when $\mathcal {E}$ has Kraus rank r, the state $\rho _{\mathcal {E}}$ has rank r. This immediately gives us a way to do compressed quantum process tomography: first prepare the Jamiołkowski state $\rho _{\mathcal {E}}$ (by adjoining an ancilla, preparing the maximally entangled state |ψ₀〉, and applying $\mathcal {E}$); then do compressed quantum state tomography on $\rho _{\mathcal {E}}$; see figure 2.

Zoom In Zoom Out Reset image size

We now show a more direct implementation of compressed quantum process tomography that is equivalent to the above procedure, but does not require an ancilla. Observe that in the above procedure, we need to estimate expectation values of the form

$$\begin{eqnarray*} &&\mathrm{Tr}\bigl((P_{A}\otimes P_{B}) \rho_{\mathcal{E}}\bigr) = \mathrm{Tr}\bigl((P_{A}\otimes P_{B}) (\mathcal{E}\otimes\mathcal{I}) (|{\psi_0}\rangle\langle{\psi_0}|)\bigr), \end{eqnarray*}$$

where P_A and P_B are Pauli matrices. By using the Kraus decomposition, it is straightforward to derive the equivalent expression

$$\begin{equation} \mathrm{Tr}((P_{A} \otimes P_{B}) \rho_{\mathcal{E}}) = \frac{1}{d} \mathrm{Tr}(P_{A} \mathcal{E}(\overline{P_{B}})), \end{equation} \tag{ 37 }$$

where the bar denotes complex conjugation in the standard basis.

We now show how to estimate the expectation value (37). Let λ_j and |ϕ_j〉 denote the eigenvalues and eigenvectors of $\overline {P_B}$. Then we have

$$\begin{eqnarray*} &&\mathrm{Tr}((P_A\otimes P_B) \rho_{\mathcal{E}}) = \frac{1}{d} \sum_{j=1}^{d} \lambda_j \mathrm{Tr}(P_A \mathcal{E}(|{\phi_j}\rangle\langle{\phi_j}|)). \end{eqnarray*}$$

To estimate this quantity, we repeat the following experiment many times, and average the results: choose j∈[d] uniformly at random, prepare the state |ϕ_j〉, apply the process $\mathcal {E}$, measure the observable P_A, and multiply the measurement result by λ_j. (See figure 3.) In this way, we learn the expectation values of the Jamiołkowski state $\rho _{\mathcal {E}}$ without using an ancilla. We then use compressed quantum state tomography to learn $\rho _{\mathcal {E}}$, and from this we recover $\mathcal {E}$. Note that since the approach described here can also be applied to the case that $\cal E$ is not trace-preserving, it applies to detector tomography as well.

Zoom In Zoom Out Reset image size

6.2. Related work

In this work, we have introduced two new estimators for compressed tomography: the matrix Dantzig selector and the matrix Lasso (equations (3) and (4)). We have proved that the sample complexity for obtaining an estimate that is accurate to within in the trace distance scales as $O(\frac {r^2 d^2}{\epsilon ^2} \log d)$ for rank-r states and that for higher rank states, the additional error is proportional to the truncation error. This error scaling is optimal up to constant multiplicative factors, and requires measuring only O(rd poly log d) Pauli expectation values, a fact we proved using the RIP [30]. We also proved that our sample complexity upper bound is within poly log d of the sample complexity of the optimal minimax estimator, where the risk function is given by a trace norm confidence interval. We showed how a modification of DFE can be used to unconditionally certify the estimate using a number of measurements which is asymptotically negligible compared to obtaining the original estimate. We numerically simulated our estimators and found that they outperform MLE, giving higher fidelity estimates from the same amount of data. Finally, we generalized our method to quantum process tomography using only Pauli measurements and preparation of product eigenstates of Pauli operators.

There are many interesting directions for future work. One major open problem is to switch the focus from two-outcome Pauli measurements to alternative measurements which are still experimentally feasible. For example, measurements in a local basis have 2ⁿ outcomes and are not directly analyzable using our techniques. It would be very interesting to give an analysis of our estimators from the perspective of such local basis measurements. One difficulty, however, is that something like the RIP is not likely to hold in this case, so we will need additional techniques.

Another direction is to find better methods for deriving error bars or confidence regions for compressed tomography. In principle, one can use the error bounds in this paper to estimate confidence regions (based on mild assumptions about the experimental setup), and then use DFE to check (unconditionally, i.e. without making any assumptions) whether the true state lies within that confidence region. However, while the error bounds in this paper have roughly the right asymptotic scaling (with the rank r and dimension d), in practice one may want to know finer properties of the confidence region, such as its shape [63, 64].

A related direction is to find methods for rank-based model selection (or, more simply, estimating the rank of the unknown state) [80]. In principle, one can always 'guess and check' the rank, using DFE; but a more sophisticated approach may perform better in practice.

Some other open problems are to tighten the gap that remains between the sample complexity upper and lower bounds, and to try to prove optimality with respect to alternative criteria, in addition to minimax risk. For example, it would be interesting to find a useful notion of average case optimality.

On the numerical side, some of the open questions are the following. First, it would be very interesting to carry out a more detailed numerical study of the performance of our estimators. While they have clearly outperformed MLE in the simulations reported here, there is no question that this is a narrow range of parameters on which we have tested these estimators. It would be interesting to do additional comparative studies between these and other estimators to see how robust these performance enhancements are. It would also be very interesting to study fast first-order solvers such as [27–29] which work particularly well for approximately low-rank matrices and which could compute estimates on a large number of qubits (ten or more).

The success of the estimators we studied depends on being able to find good values for the parameters λ and μ. While we have used simple heuristics to pick particular values, a detailed study of the optimal values for these parameters could only improve the quality of our estimators. Moreover, MLE seems to enjoy the same 'plateau' phenomenon, where the quality of the estimate is insensitive to m above a certain cutoff. This leads us to speculate that this is a generic phenomenon among many estimators, and that perhaps there are even better choices for estimators than the ones we benchmark here.

We thank R Blume-Kohout, M Kleinmann and T Monz for discussions, B Brown for some preliminary numerical investigations and S Glancy for helpful comments on a draft of the paper. We additionally thank B Englert for hosting the workshop on quantum tomography in Singapore where some of this work was completed. STF was supported by the IARPA MQCO program. DG was supported by the Excellence Initiative of the German Federal and State Governments (grant ZUK 43). Contributions to this work by NIST, an agency of the US government, are not subject to copyright laws. Although this paper refers to specific computer software, NIST does not endorse commercial products. JE was supported by EURYI, the EU (Q-Essence) and BMBF (QuOReP).