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

Intuitively, if a density operator has small rank, then it should be easier to estimate from experimental data, since in this case only a few eigenvectors need to be learned. We prove two complementary results that confirm this intuition. Firstly, we show that a low-rank density matrix can be estimated using fewer copies of the state, i.e. the sample complexity of tomography decreases with the rank. Secondly, we show that unknown low-rank states can be reconstructed from an incomplete set of measurements, using techniques from compressed sensing and matrix completion. These techniques use simple Pauli measurements, and their output can be certified without making any assumptions about the unknown state. In this paper, we present a new theoretical analysis of compressed tomography, based on the restricted isometry property for low-rank matrices. Using these tools, we obtain near-optimal error bounds for the realistic situation where the data contain noise due to finite statistics, and the density matrix is full-rank with decaying eigenvalues. We also obtain upper bounds on the sample complexity of compressed tomography, and almost-matching lower bounds on the sample complexity of any procedure using adaptive sequences of Pauli measurements. Using numerical simulations, we compare the performance of two compressed sensing estimators—the matrix Dantzig selector and the matrix Lasso—with standard maximum-likelihood estimation (MLE). We find that, given comparable experimental resources, the compressed sensing estimators consistently produce higher fidelity state reconstructions than MLE. In addition, the use of an incomplete set of measurements leads to faster classical processing with no loss of accuracy. Finally, we show how to certify the accuracy of a low-rank estimate using direct fidelity estimation, and describe a method for compressed quantum process tomography that works for processes with small Kraus rank and requires only Pauli eigenstate preparations and Pauli measurements.


I. INTRODUCTION
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 in 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][3][4][5][6][7][8][9][10][11][12][13][14][15][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 n , and the state of the system is described by a density matrix with d 2 = 4 n 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 VI B.
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 2 ) for a general mixed state.While this intuition is indeed correct [20][21][22][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 lab; furthermore, one then has to find a pure state consistent with measured data [24], preferably by some procedure that is computationally efficient (note that finding minimum-rank solutions is NP-hard in general [25]).
Compressed tomography provides a solution that meets all of these practical requirements [17,18].It requires measurements of 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 generalpurpose SDP solvers [26], or specialized algorithms for larger instances [27][28][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 lowrank state.
Here we follow up on Refs.[17,18] by giving a stronger (and completely different) theoretical analysis, which is based on the restricted isometry property (RIP) [30][31][32].This answers a number of questions that could not be addressed satisfactorily using the earlier techniques based on dual certificates.First, 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.Second, 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 the 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, so 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.
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 C d×d can be embedded into O(rd poly log d) dimensions, with small distortion in the 2-norm.An embedding that does this is said to satisfy the restricted isometry property (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 so-called "universal" methods for low-rank matrix recovery: there exists a fixed set of O(rd poly log d) Pauli measurements, that 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 fullrank, 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 Ref. [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 2 d 2 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 2 d 2 / log d) copies are needed to obtain an estimate with constant accuracy in the trace distance.This generalizes a result from Ref. [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, when 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, so long as m ≥ Ω(rd poly log d).
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 ρ, and use DFE to check whether ρ 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 ρ regardless of whether it was ob-tained 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 (a.k.a. the matrix Dantzig selector), least squares with trace-norm regularization (a.k.a. the matrix Lasso), as well as a standard maximum likelihood estimation (MLE) [35][36][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 are (as mentioned above) fairly insensitive to the number of measurement settings m (assuming the total time available to run the experiment is fixed).So by choosing m d 2 , one still obtains accurate estimates, but with much faster classical post-processing, since the size of the data set scales like O(m) rather than O(d 2 ).
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 2 ), and (2) using trace minimization or regularization to reconstruct low-rank solutions.Usually one does both of these things.Now, suppose the goal is to reconstruct a lowrank 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 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 to Ref. [19], our method is not restricted to processes that are element-wise sparse in some known basis, as discussed in Section VI B. 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.
However, recently there has been a flurry of work which seeks to transcend the standard MLE methods and improve on them in various ways.Our contributions can also be seen in this context.
One way in which alternatives to MLE are being pursued is through what we call full rank methods.Here the idea is somewhat antithetical to ours: the goal is to output a full rank density operator, rather than a rank deficient one.This is desirable in a context where one cannot make the approximation that rare events will never happen.Blume-Kohout's hedged MLE [53] and Bayesian mean estimation [52] are good examples of this type of estimator, as are the minimax estimator of Ref. [54] and the so-called Max-Ent estimators [55][56][57][58].The latter are specifically for the setting where the measurement data are not informationally complete, and one tries to minimize the bias of the estimate by maximizing the entropy along the directions where one has no knowledge.
By contrast, our low rank methods do not attempt to reconstruct the complete density matrix, but only a rank-r approximation, which is accurate when the true state is close to low-rank.From this perspective, our methods can be seen as a sort of Occam's Razor, using as few fit parameters as possible while still agreeing with the data [59].Furthermore, as we show here and elsewhere [17], informationally incomplete measurements can still provide faithful state reconstructions up to a small truncation error.
One additional feature of our methods is that we are deeply concerned with the feasibility of our estimators for a moderately large number of qubits (say, [10][11][12][13][14][15].In contrast to most of the existing literature, we adopt the perspective that it is not enough for an estimator to be asymptotically efficient in the number of copies for fixed d.We also want the scaling with respect to d to be optimal.We specifically take advantage of the fact that many states and processes are described by low rank objects to reduce this complexity.In this respect, our methods are similar to tomographic protocols that are tailored to special ansatz classes of states, such as those recently developed for use with permutation-invariant states [60], matrix product states [61] or multi-scale entangled states [62].
Our error bounds are somewhat unique as well.Most prior work on error bounds used either standard resampling techniques or Bayesian methods [42,43,45,48,50,52]. Very recently, Christandl & Renner and Blume-Kohout independently derived two closely related approaches for obtaining confidence regions that satisfy or nearly satisfy certain optimality criteria [63,64].Especially these latter approaches can give very tight error bounds on an estimate, but they can be computationally challenging to implement for large systems.The error bounds which most closely resemble ours are of the "large deviation type"; see for example the discussion in Ref. [38].This is true for the new improved error bounds, as well as the original bounds proven in Refs.[17,18].These types of bounds are much easier to calculate in practice, which agrees with our philosophy on computational complexity, but may be somewhat looser than the optimal error bounds obtainable through other more computationally intensive methods such as those of Refs.[63,64].

B. Notation and Outline
We denote Pauli operators by P or P i .We define [n] = {1, . . ., n}.The norms we use are the standard Euclidean vector norm x 2 , the Frobenius norm X F = Tr(X † X), the operator norm X = λ max (X † X) and the trace norm X tr = Tr|X|, where |X| = √ X † X.The unknown "true" state is denoted ρ and any estimators for ρ are given a hat: ρ.The expectation value of a random variable X is denoted EX.We denote by H d the set of d × d Hermitian matrices.
The paper is organized as follows.In Section II we detail the estimators and error bounds, then upper bound the sample complexity.In Section III we derive lower bounds on the sample complexity.In Section IV we find an efficient method of certifying the state estimate.In Section V we detail our numerical investigations.We show how our scheme can be applied to quantum channels in Section VI and conclude in Section VII.

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

A. Random Pauli Measurements
Consider a system of n qubits, and let d = 2 n .Let P be the set of all d 2 Pauli operators, i.e., matrices of the form Our tomography scheme works as follows.First, choose m Pauli operators, P 1 , . . ., P m , by sampling independently and uniformly at random from 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 6 d and t ∼ r 2 d 2 log d.) To state this more concisely, we introduce some notation.Define the sampling operator to be a linear map The normalization is chosen so that EA * A = I, where I denotes the identity superoperator and A * is the adjoint of A. We can then write the output of our measurement procedure as a vector where z represents statistical noise due to the finite number of samples, or even due to an adversary.

B. 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 ∈ H d that fits the data y while minimizing the trace norm X 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 V.) The first estimator is obtained by constrained traceminimization (a.k.a. the matrix Dantzig selector): 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 (a.k.a. the matrix Lasso): 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), (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 tradeoffs that result from this.

C. Error Bounds
In previous work on compressed tomography [17,18], error bounds were proved by constructing a "dual certificate" [66] (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 restricted isometry property (RIP).The RIP for lowrank matrices was first introduced in Ref. [32], and was recently shown to hold for random Pauli measurements [30].
We say that the sampling operator A satisfies the rankr restricted isometry property if there exists some constant 0 ≤ δ r < 1 such that, for all X ∈ C d×d with rank r, For our purposes, we can assume that X is Hermitian.
Note that this notion of RIP is analogous to the one used in Ref. [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 6 d (for some absolute constant C); this was recently shown in Ref. [30,Theorem 2.1].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 Ref. [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 II, let C, C 0 , C 1 , C 0 and C 1 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 A be the random Pauli sampling operator (1) with m ≥ Crd log 6 d.Then, with high probability, the following holds: Let ρDS be the matrix Dantzig selector (3), and choose Alternatively, let ρLasso be the matrix Lasso (4), and choose µ so that 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 Ref. [30], and a straightforward modification of Lemma 3.2 in Ref. [31] to bound the error in trace norm rather than Frobenius norm (this is similar to the proof of Theorem 5 in Ref. [67]).
The modification of Lemma 3.2 in Ref. [31] is as follows1 .(For the remainder of this section, equation numbers of the form (III.x) refer to Ref. [31].)In the case of the Dantzig selector, let H = ρDS −ρ.Following equation (III.8),we can get the following bound: where we used the triangle inequality and equation (III.8).Then, at the end of the proof, we write: where we used Cauchy-Schwarz, the fact that H 0 and H 1 are orthogonal, the bound on H 0 + H 1 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 Ref. [31].

D. 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 ( rd ε )2 log d 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 4 /ε 2 ) copies are sufficient to estimate a full-rank state with accuracy ε in trace distance 2 .
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 for some constants C 4 > 1 and ε ≤ C 0 .(For the matrix Lasso, substitute C 0 for C 0 in these equations.)We construct an estimate of A(ρ) 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 A(ρ), and let z = y − A(ρ).Everything else is defined exactly as in Theorem 1.
Theorem 2. Given t = O ( rd ε ) 2 log d copies of ρ as in Eq. ( 8) and measured as discussed above, then the following holds with high probability over the measurement outcomes: Let ρDS be the matrix Dantzig selector (3), and set λ = ε/(C 0 r) for some ε > 0. Then Alternatively, let ρLasso be the matrix Lasso (4), and set µ = ε/(C 0 r).Then 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 satisfies A * (y − A(ρ)) ≤ λ? Then one can apply Theorem 1 to get an error bound for our estimate of ρ.
Let t be the number of copies of ρ.Say we fix the measurement operator A, i.e., we fix the choice of the Pauli observables P 1 , . . ., 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 j'th repetition of the experiment that measures the i'th Pauli observable P i .Note that EB ij = Tr(P i ρ).Then construct the vector y ∈ R m containing the estimated expectation values (scaled by d/m): Note that Ey = A(ρ).
We will bound the deviation A * (y −A(ρ)) , using the matrix Bernstein inequality.First we write and also We can now write A * (y − A(ρ)) as a sum of independent (but not identical) matrix-valued random variables: Note that EX ij = 0 and X ij ≤ 2d/t =: R. Also, for the second moment we have Then we have Now the matrix Bernstein inequality (Theorem 1.4 in [68]) implies that (where we assumed λ ≤ 1).
For the matrix Dantzig selector, we set λ = ε/(C 0 r), and we get that, for any t which is exponentially small in C 4 .Plugging into Theorem 1 completes the proof of our claim.A similar argument works for the matrix Lasso.

III. LOWER BOUNDS
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 (Imprecise formulation).The number of copies t needs to grow at least as fast as Ω r 2 d 2 / log d 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 the protocol has access to t copies of an unknown state ρ ∈ Σ on which it performs measurements one by one.At the ith step, it has to decide which observable to measure.Let us restrict ourselves to binary POVM measurements {Π i , 1 − Π i }, where each Π i satisfies 0 ≤ Π i ≤ 1 and may be chosen from a set 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 ρ(Y 1 , . . ., Y t ) ∈ Σ.
In other words, an estimation protocol is specified by a set of functions Consider a distance measure ∆ : Σ × Σ → R on the states in Σ. (For example, this could be the trace distance or the infidelity; we need not specify.)Suppose that the maximum deviation we wish to tolerate between our unknown state and the estimate is given by > 0. Now define the minimax risk where the infimum is over all estimation procedures ρ, Π i on t copies with estimator ρ 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 P, then the number of copies t must be large too to avoid having a large minimax risk.Lemma 4. Assume there are states ρ 1 , . . ., ρ s ⊂ Σ such that Then the minimax risk as defined in (17) fulfills Now combine Fano's inequality the data processing inequality in terms of the mutual information I(X; Y ) := H(X) − H(X|Y ), and the fact that Let h(p) be the binary entropy and recall the standard estimate Combine that with the chain rule [69, Theorem 2.5.1]: The advertised bound follows.
By applying the following lemma s < e crd times, we can randomly create a set of s states each with rank r that satisfy the conditions of Lemma 4 in terms of the trace distance and the set of Pauli measurements.
Lemma 5.For any 0 < < 1 − r d , let ρ 1 , . . ., ρ s be normalized 3 rank-r projections on C d , where s < e c( )rd and c( ) is specified below.Then there exists a normalized rank-r projection ρ such that: Here, α 2 = O log d rd , the P k are n-qubit Pauli operators, and Proof: Let ρ 0 be some normalized rank-r projection and choose ρ according to 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 We can bound the magnitude of f using the pinching inequality: From this we can bound the expectation value of f over the special orthogonal group: Next we get an upper bound of 4 √ r on the Lipschitz constant of f with respect to the Frobenius norm: 3 Normalized to have trace 1. 4 For the duration of this proof, the letter O denotes an element of SO(d) instead of the asymptotic big-O notation.
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 get that for all t > 0, where the constants are given by c 1 = ln π/8 and c 2 = 1/8.Now choose t = 2(1 − r/d) − 2 and apply the union bound to obtain Pr (19) does not hold The upper bound on follows from the requirement that t > 0. This shows that ρ indeed satisfies Eq. ( 19).Now we move on to Eq. ( 20).For any non-identity Pauli matrix P k , define a function Clearly, we have 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 Lévy's Lemma then gives for all t > 0 Pr |trP k ρ| > t ≤ e c1 exp − c 2 t 2 rd 4 .
We remark that a version of Lemma 5 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.
Combining the two previous lemmas yields a precise formulation of Theorem 3.
Theorem 6 (Precise version of Theorem 3).Fix ∈ (0, 1 − r d ) and δ ∈ [0, 1).Then for our bound to allow for M * ( ) ≤ δ we must have that the number of copies t of ρ grows like where the implicit constant depends on δ and .

IV. CERTIFYING THE STATE ESTIMATE
Here we sketch how the technique of direct fidelity estimation (DFE), introduced in Refs.[33,34] for pure states, can be used to estimate the fidelity between a low rank estimate ρ and the true state ρ.The only assumption we make is that ρ is a positive semidefinite matrix with Tr(ρ) ≤ 1 and rank(ρ) = r.No assumption at all is needed on ρ.In fact, we also do not assume that we obtained the estimate ρ from any of the estimators which were discussed previously.Our certification procedure works regardless of how one obtains ρ, and so it even applies to the situation where ρ 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 Ref. [33], we simply redefine the probability distribution as Pr(i) = | φ j |P i |φ k | 2 /d and the random variable X = Tr(P i ρ)/ φ k |P i |φ j .It is easy to check that E(X) = φ j |ρ|φ k and the variance of X is at most one.Then all of the conclusions in Refs.[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( 1 ) 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 ρ, given by where for brevity we define G = √ ρρ √ ρ. (Note that some authors use the square root of this quantity as the fidelity.Our definition matches Ref. [33].)The asymptotic cost in sample complexity is far less than the cost of initially obtaining the estimate ρ when r is sufficiently small compared to d.
Theorem 7. Given a state estimate ρ with rank(ρ) = r, the number of copies t of the state ρ required for estimating F (ρ, ρ) to within ±ε with probability 1 − δ using single-copy Pauli measurements satisfies Proof: The result uses the DFE protocol of Refs.[33,34], modified as mentioned above, where the states |φ j are the eigenstates of ρ.Expand ρ in its eigenbasis as ρ = r j=1 λ j |φ j φ j |.Then we have For all 1 ≤ j ≤ k ≤ r, we use direct fidelity estimation to obtain an estimate ĝjk of the matrix element φ j |ρ|φ k , up to some additive error jk that is bounded by a constant, If each estimate is accurate with probability 1 − 2δ/(r 2 + 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 Let ĝjk = ĝ * kj , and let Ĝ = r j,k=1 be our estimate for G. Finally, let Ĝ+ = [ Ĝ] + be the positive part of the Hermitian matrix Ĝ, and let F = Tr Ĝ+ 2 be our estimate of the fidelity F (ρ, ρ).Note that we may assume that F ≤ 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 Ĝ as a perturbation of G, Ĝ = G + E, where First notice that the Frobenius norm of this perturbation is small, (If ρ is subnormalized, then this last equality becomes an inequality.)Next observe that where we define Using the reverse triangle inequality, we can bound ∆ in terms of the trace norm Using [71, Thm.X. 1.3] in the first step, we find that where the second bound follows from the Cauchy-Schwarz inequality on the vector of eigenvalues of |G − Ĝ+ |.Using the Jordan decomposition of a Hermitian matrix into a difference of positive matrices Then by the triangle inequality and positivity of G, Using the standard estimate E tr ≤ √ r E F , we find This gives our desired error bound, albeit in terms of 0 instead of the final quantity ε.To get our total error to vanish, we take ε = 2r 3/4 √ 2 0 , which gives us the final scaling in the sample complexity.
As a final remark, we note that by computing ∆ for the special case ρ = ρ = 1/r, = 0 1/r, we find that so the error bound for this protocol is tight with respect to the scaling in 0 (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.

V. NUMERICAL SIMULATIONS
We numerically simulated the reconstruction of a quantum state given Pauli measurements and the estimators from Eqs. ( 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 Eqs. ( 3) and ( 4) given by ρ DS = arg min X 0 TrX s.t.A * (A(X) − y) ≤ λ , (34) and ρ Lasso = arg min Moreover, whenever the trace of the resulting estimate is is less than one, we renormalize the state, ρ ← ρ/Tr(ρ).
A. Setting the Estimator Parameters λ and µ From Thm. 1 we know roughly how we should choose the free parameters λ and µ, modulo the unknown constants C i , for which we will have to make an empirical guess.We guess that the log factor is an artifact of the analysis, and that the state is nearly pure.Then we could choose λ ∼ µ ∼ d/ √ t.For the matrix Dantzig selector of Eq. ( 34), we specifically chose λ = 3d/ √ t for our simulations.However, this did not lead to very good performance for the matrix Lasso of Eq. ( 35) when m was large.We chose instead µ = 4m/ √ t, which agrees with λ when m ∼ d, as it can be for nearly pure states.
We stress that very little effort was made to optimize these choices for the parameters.We simply picked a few integer values for the constants in front (namely 2, . . ., 5), and chose the constant that worked best upon a casual inspection of a few data points.We leave open the problem of determining what the optimal choices are for λ and µ.

B. Time Needed to Switch Measurement Settings
Real measurements take time.In a laboratory setting, time is a precious commodity for a possibly non-obvious reason.An underlying assumption in the way we typically formulate tomography is that the unknown state ρ can be prepared identically many times.However, the parameters governing the experiment often drift over a certain timescale, and this means that beyond this characteristic time, it is no longer appropriate to describe the measured ensemble by the symmetric product state ρ ⊗t .
To account for this characteristic timescale, we introduce the following simple model.We assume that the entire experiment takes place in a fixed amount of time T .In a real experiment, this is the timescale beyond which we cannot confidently prepare the same state ρ, or it is the amount of time it takes for the student in the lab to graduate, or for the grant funding to run out, whichever comes first.Within this time limit, we can either change the measurement settings, or we can take more samples.We assume that there is a unit cost to taking a single measurement, but that switching to a different measurement configuration takes an amount of time c.
At the one extreme, the switching cost c to change measurement settings could be quite large, so that it is only barely possible to conduct enough independent measurements that tomography is possible within the allotted time T .In this case, we expect that compressed tomography will outperform standard methods, since these methods have not been optimized for the case of incomplete data.At the other extreme, the switching cost c could be set to 0 (or some other small value), in which case we are only accounting for the cost of sampling.Here it is not clear which method is better, for the following reason.Although the standard methods are able to take a sufficient amount of data in this case, each observable will attain comparatively little precision relative to the fewer observables measured with higher precision in the case of compressed tomography for the same fixed amount of time T .One of our goals is to investigate if there are any tradeoffs in this simple model.
For the relative cost c between switching measurement settings and taking one sample, we use c = 20, a number inspired by the current ion trap experiments done by the Blatt group in Innsbruck.However, we note that the conclusions don't seem to be very sensitive to this choice, as long as we don't do something outrageous like c > t, which would preclude measuring more than even one observable.

C. Other Simulation Parameters
We consider the following ensembles of quantum states.We initially choose n = 5 qubit states from the invariant Haar measure on 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 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 ρ relative to the underlying true state ρ, namely the (squared) fidelity √ ρ√ ρ 2 tr and the trace distance 1  2 ρ − ρ tr .We used the interior-point solver SeDuMi [26] to do the optimizations in Eqs.(34) and (35).Although these algorithms are not suitable for large numbers of qubits, they produce very accurate solutions, which allows 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 SVT, TFOCS, or FPCA to solve these convex programs [27][28][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 Ref. [73], which has previously been used in e.g.Ref. [74].This method converged on every example we tried, so we did not have to use the more sophisticated "diluted" method of Ref. [75] that guarantees convergence.
For the purposes of our simulations, we sampled our Pauli operators without replacement.

D. Results and Analysis
The data are depicted in Fig. 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.
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 2 ) 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 wellsuited for the task.
When the total amount of time T is very small (Fig. 1a), then there is a large advantage to using compressed sensing.In this regime, there is barely enough time to make one measurement per setting if we insist on making an informationally complete measurement, so the measurement statistics aren't even Gaussian.Thus, the fluctuations are so large that trying to fit a density operator to these data just leads to poor results because you wind up fitting only noise.While the advantage of compressed sensing is clear in this regime, it is not applicable when we are interested in extremely high-fidelity state reconstructions (namely greater than 95%).
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 tradeoff 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  and c) are for a total measurement time of T = 41k, 80k, and 270k respectively, in units where the cost to measure a single copy is one unit of time, and the cost to switch measurement settings is c = 20 units.The three estimators plotted are the matrix Dantzig selector (Eq.(34), blue), the matrix Lasso (Eq.( 35), red), and a standard MLE (green).Each bar shows the mean and standard deviation from 120 Haar-random 5-qubit states with 1% local depolarizing noise.Our estimators consistently outperform MLE, even after optimizing the MLE over m.See the main text for more details.
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 Eqs. ( 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.

VI. PROCESS TOMOGRAPHY
Compressed sensing techniques can also be applied to quantum process tomography.Here, our method has an advantage when the unknown quantum process has small Kraus rank, meaning 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 n .Let E be a completely positive (CP) map from C d×d to C d×d , and suppose that E has Kraus rank r.Using compressed sensing, we will show how to characterize E using m = O(rd 2 log d) settings.(For comparison, standard process tomography requires d 4 settings, since E contains d 4 independent parameters in general.)Furthermore, our compressed sensing method requires only the ability to prepare product eigenstates of Pauli operators and perform 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 2 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 elementwise 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.

A. Our Method
First, recall the Jamio lkowski isomorphism [76]: the process E is completely and uniquely characterized by the state Note that when E has Kraus rank r, the state ρ E has rank r.This immediately gives us a way to do compressed quantum process tomography: first prepare the Jamio lkowski state ρ E (by adjoining an ancilla, preparing the maximally entangled state |ψ 0 , and applying E); then do compressed quantum state tomography on ρ E ; see Figure 2.
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 Tr (P A ⊗ P B )ρ E = Tr (P A ⊗ P B )(E ⊗ I)(|ψ 0 ψ 0 , where P A and P B are Pauli matrices.By using the Kraus decomposition, it is straightforward to derive the equivalent expression

Tr((P
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 P B .Then we have

Tr((P
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 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 lkowski state ρ E without using an ancilla.We then use compressed quantum state tomography to learn ρ E , and from this we recover E.

B. Related Work
Our method is somewhat different from the method described in Ref. [19].Essentially the difference is that our method works for any quantum process with small Kraus rank, whereas the method of Shabani et al. works for a quantum process that is elementwise sparse in a known basis (provided this basis is experimentally accessible in a certain sense).The main advantage of the Shabani et al. method is that it can be much faster: for a quantum process E that is s-sparse (i.e., has s nonzero matrix elements), it requires only O(s log d) settings.The main disadvantage is that it requires more prior knowledge about E, and is more difficult to apply.While it has been demonstrated in a number of scenarios, there does not seem to be a general recipe for designing measurements that are both experimentally feasible and effective in the sparse basis of E.
To clarify these issues, we now briefly review the Shabani et al. method.We assume that we know a basis Γ = {Γ α | α ∈ [d 2 ]} in which the process E is s-sparse.For example, when E is close to some unitary evolution U , one can construct Γ using the SVD basis of U .This guarantees that, if E contains no noise, it will be perfectly sparse in the basis Γ.However, in practice, E will contain noise, which need not be sparse in the basis Γ; any non-sparse components will not in general be estimated accurately.The success of the Shabani et al. method therefore rests on the assumption that the noise is also sparse in the basis Γ.Although this assumption has been verified in a few specific scenarios, it seems less clear why it should hold in general.By comparison, our method simply assumes that the noise is described by a process matrix that is low rank; this can be rigorously justified for any noise process that involves only local interactions or few-body processes.
The other complication with the Shabani et al. method concerns the design of the state preparations and measurements.On one hand, these must satisfy the RIP condition for s-sparse matrices over the basis Γ; on the other hand, they must be easy to implement experimentally.This has been demonstrated in some cases, by using states and measurements that are "random" enough to show concentration of measure behaviors, but also have tensor product structure.However, these constructions are not guaranteed to work in general for an arbitrary basis Γ.
We leave open the problem of doing a comparative study between these and other methods [77,78], akin to Ref. [79].

VII. CONCLUSION
In this work, we have introduced two new estimators for compressed tomography: the matrix Dantzig selector and the matrix Lasso (Eqs.( 3) and (4).)We have proved that the sample complexity for obtaining an estimate that is accurate to within in the trace distance scales like O( r 2 d 2 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 direct fidelity estimation 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 open questions that remain.On the theoretical side, one open problem is of course to tighten the gap that remains between the sample complexity upper and lower bounds.Another open problem is to try to prove optimality with respect to alternative criteria other than minimax risk.For example, it would be interesting to find a useful notion of average case optimality.
One major open problem is to switch focus from twooutcome Pauli measurements to alternative measurements which are still experimentally feasible.For example, measurements in a local basis have 2 n 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.
On the numerical side, some of the open questions are the following.First, it would be very interesting to do 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 Refs.[27][28][29] which could compute estimates on a large number of qubits (10 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.

20 FIG. 1 :
FIG.1: Fidelity and trace (inset) as a function of m, the number of sampled Paulis.Plots a), b) and c) are for a total measurement time of T = 41k, 80k, and 270k respectively, in units where the cost to measure a single copy is one unit of time, and the cost to switch measurement settings is c = 20 units.The three estimators plotted are the matrix Dantzig selector (Eq.(34), blue), the matrix Lasso (Eq.(35), red), and a standard MLE (green).Each bar shows the mean and standard deviation from 120 Haar-random 5-qubit states with 1% local depolarizing noise.Our estimators consistently outperform MLE, even after optimizing the MLE over m.See the main text for more details.

FIG. 2 :
FIG.2: Compressed quantum process tomography using an ancilla.The quantum circuit represents a single measurement setting, where one measures the observable PA on the system and PB on the ancilla.
Let X be a random variable taking values uniformly in[s].Let Y 1 , . . ., Y t be random variables, Y i describing the outcome of the ith measurement performed on ρ X .Define Proof: Compressed quantum process tomography, implemented directly without an ancilla.Here one prepares a random eigenstate of PB, applies the process E, and measures the observable PA on the output.