We establish methods for quantum state tomography based on compressed sensing. These methods are specialized for quantum states that are fairly pure, and they offer a significant performance improvement on large quantum systems. In particular, they are able to reconstruct an unknown density matrix of dimension d and rank r using O(rd log 2 d) measurement settings, compared to standard methods that require d 2 settings. Our methods have several features that make them amenable to experimental implementation: they require only simple Pauli measurements, use fast convex optimization, are stable against noise, and can be applied to states that are only approximately low-rank. The acquired data can be used to certify that the state is indeed close to pure, so no a priori assumptions are needed. We present both theoretical bounds and numerical simulations. The tasks of reconstructing the quantum states and processes produced by physical systems -known respectively as quantum state and process tomography [1] -are of increasing importance in physics and especially in quantum information science. Tomography has been used to characterize the quantum state of trapped ions [2] and an optical entangling gate [3] among many other implementations. But a fundamental difficulty in performing tomography on many-body systems is the exponential growth in the state space dimension. For example, to get a maximum-likelihood estimate of a quantum state of 8 ions, Ref. [2] required hundreds of thousands of measurements and weeks of post-processing. Still, one might hope to overcome this obstacle, because the vast majority of quantum states are not of physical interest. Rather, one is often interested in states with special properties: pure states, states with particular symmetries, ground states of local Hamiltonians, etc., and tomography might be more efficient in such special cases [4] . In particular, consider pure or nearly pure quantum states, i.e., states with low entropy. More precisely, consider a quantum state that is essentially supported on an r-dimensional space, meaning the density matrix is close (in a given norm) to a matrix of rank r, where r is small. Such states arise in very common physical settings, e.g. a pure state subject to a local noise process [20]. A standard implementation of tomography [5, 6] would use d 2 or more measurement settings, where d = 2 n for an nqubit system. But a simple parameter counting argument suggests that O(rd) settings could possibly suffice -a significant improvement. However, it is not clear how to achieve this performance in practice, i.e., how to choose these measurements, or how to efficiently reconstruct the density matrix. For instance, the problem of finding a minimum-rank matrix subject to linear constraints is NP-hard in general [7] . In addition to a reduction in experimental complexity, one might hope that a post-processing algorithm which takes as input only O(rd) ≪ d 2 numbers could be tuned to run considerably faster than standard methods. Since the output of the procedure is a low-rank approximation to the density operator and only requires O(rd) numbers be specified, it becomes conceivable that the run time scales better than O(d 2 ), clearly impossible for naive approaches using dense matrices. In this Letter, we introduce a method to achieve such drastic reductions in measurement complexity, together with efficient algorithms for post-processing. The approach further develops ideas that have recently been studied under the label of "compressed sensing". Compressed sensing [8] provides techniques for recovering a sparse vector from a small number of measurements [9] . Here, sparsity means that this vector contains only a few non-zero entries in a specified basis, and the measurements are linear functions of its entries. When the measurements are chosen at random (in a certain precise sense), then with high probability two surprising things happen: the vector is uniquely determined by a small number of measurements, and it can be recovered by an efficient convex optimization algorithm [8] . Matrix completion [10] [11] [12] is a generalization of compressed sensing from vectors to matrices. Here, one recovers certain "incoherent" low-rank matrices X from a small number of matrix elements X i,j . The problem of low-rank quantum state tomography bears a strong resemblance to matrix completion. However, there are important differences. We wish to use measurements that can be more easily implemented in an experiment than obtaining elements ρ i,j of density matrices. Previous results [10] [11] [12] cannot be applied to this more general situation. We would also like to avoid any unnatural incoherence assumptions crucial in prior work [10] . Our first result is a protocol for tomography that overcomes both of these difficulties: it uses Pauli measurements only, and it works for arbitrary density matrices. We prove that only O(rd log 2 d) measurement settings suffice. What is more, our proof introduces some new techniques, which both generalize and vastly simplify the previous work on matrix completion. We sketch the proof here; a more complete version appears in [25]. This provides the basic theoretical justification for our method of doing tomography. We then consider a number of practical issues. In a real experiment, the measurements are noisy, and the true state is only approximately low-rank. We show that our method is robust to these sources of error. We also describe ways to certify that a state is nearly pure without any a priori assumptions. Finally, we present fast algorithms for reconstructing the density matrix from the measurement statistics based on semidefinite programming -a feature not present in earlier methods for pure-state tomography [4] [5] [6] . These are adapted from algorithms for matrix completion [14] , and they are much faster than standard interior-point solvers. Reconstructing a low-rank density matrix for 8 qubits takes about one minute on an ordinary laptop computer. While our methods do not overcome the exponential growth in measurement complexity (which is provably impossible for any protocol capable of handling generic pure states), they do significantly push the boundary of what can be done in a realistic setting. Our techniques also apply to process tomography: to characterize an unknown quantum process E, prepare the Jamiołkowski state ρ E , and perform state tomography on ρ E . Our methods work when E can approximately be written as a sum of only a few Kraus operators, because this implies that ρ E has small rank. Matrix recovery using Pauli measurements. We consider the case of n spin-1/2 systems in an unknown state ρ [16]. An n-qubit Pauli matrix is of the form w = n i=1 w i , where The protocol proceeds as follows: choose 2 ] at random and measure the expectation values tr ρw(A i ). One then solves a convex optimization problem: minimize σ tr [17] subject to Theorem 1 (Low-rank tomography) Let ρ be an arbitrary state of rank r. If m = cdr log 2 d randomly chosen Pauli expectations are known, then ρ can be uniquely reconstructed by solving the convex optimization problem (1) with probability of failure exponentially small in c. The proof is inspired by, but technically very different from, earlier work on matrix completion [10] . Our methods are more general, can be tuned to give tighter bounds, and are much more compact, allowing us to present a fairly complete argument in this Letter. A more detailed presentation of this technique -covering the reconstruction of low-rank matrices from few expansion coefficients w.r.t. general operator bases (not just Pauli matrices or matrix elements) -will be published elsewhere [25]. Proof: Here we sketch the argument and explain the main ideas; detailed calculations are in the EPAPS supplement. Note that the linear constraints (1) depend only on the projection of ρ onto the span of the measured observables w(A 1 ), . . . , w(A m ). This is precisely the range of the "sampling operator" Indeed, the convex program can be written as min σ σ tr s.t. Rσ = Rρ. Evidently, the solution is unique if for all deviations ∆ := σ − ρ away from ρ either R∆ = 0 or ρ + ∆ tr > ρ tr . We will ascertain this by using a basic idea from convex optimization: constructing a strict subgradient Y for the norm. A matrix Y is a strict subgradient if ρ + ∆ tr > ρ tr + tr Y ∆ for all ∆ = 0. The main contribution below is a method for constructing such a Y which is also in the range of R. For then R∆ = 0 implies that ∆ is orthogonal to the range of R, thus tr Y ∆ = 0 and the subgradient condition reads ρ + ∆ tr > ρ tr . This implies uniqueness. (In fact, it is sufficient to approximate the subgradient condition in a certain sense). Let E be the projection onto the range of ρ, let T be the space spanned by those operators whose row or column space is contained in range ρ. Let P T be the projection onto T , P ⊥ T onto the orthogonal complement. Decompose ∆ = ∆ T +∆ ⊥ T , the parts of ∆ that lie in the subspaces T and T ⊥ . We distinguish two cases: . Case (i) is easier. In this case, ∆ is well-approximated by ∆ T and essentially we only have to show that the restriction A := P T RP T of R to T is invertible. Using a noncommutative large deviation bound (see EPAPS supplement), where κ = m/(dr) [17]. Hence the probability that A −  First, suppose such a Y exists. Then a simple calculation (see EPAPS) using the condition (ii) shows that R∆ = 0 indeed implies ρ + ∆ tr > ρ tr as hinted at above. This proves uniqueness in case (ii). The difficult part consists in showing that an almost-subgradient exists. To this end, we design a recursive process (the "golfing scheme" [25]) which converges to a subgradient exponentially fast. Assume we draw l batches of κ 0 rd Pauli observables independently at random (κ 0 will be chosen later). Define recursively X 0 = E, Y = Y l . Let R i be the sampling operator associated with the ith batch, and A i its restriction to T . Assume that in each run Denote the probability of this event not occurring by p 2 . Then the first part of (3), as soon as l ≥ log 2 (2d 2 √ r). We turn to the second part. Again using large-deviation techniques (EPAPS) we find P ⊥ T R i X i−1 ≤ 1/(4 √ r) X i−1 2 with some (high) probability (1 − p 3 ). Therefore: which is the second part of (3). Lastly, we have to bound the total probability of failure p f ≤ p 1 + p 2 + p 3 . Set κ 0 = 64µ(1 + ln(8dl)), which means that m = dr(ln d) 2 O(1) coefficients will be sampled in total. A simple calculation gives p f ≤ e −µ . This completes the proof of our main result. In the remaining space, we address the important aspects of resilience against noise, certified tomography, and numerical performance. Owing to space limitations, the presentation will focus on conceptual issues, with the details in [24]. Robustness to noise. Realistic situations will differ from the previous case in two regards. First, the true state ρ t may not be low-rank, but only well approximated by a state ρ of rank r: ρ t − ρ 2 ≤ ε 1 . Second, due to systematic and statistical noise, the available estimates for the Pauli expectations are not exactly tr ρ t w(a), but of the form tr ωw(a) for some matrix ω. Assume Rω −Rρ t 2 ≤ ε 2 (in practical situations, ε 2 may be estimated from the error bars associated with the individual Pauli expectation values [21]). In order to get an estimate for ρ t , choose some λ ≥ 1 and ε ≥ λ( d 2 /m)ε 1 + ε 2 , and solve the convex program min σ tr , subject to Rσ − Rω 2 ≤ ε. Observation 1 (Robustness to noise) Let ρ t be an approximately low-rank state as described above. Suppose m = cdr log 2 d randomly chosen Pauli expectations are known up to an error of ε as in (6), and let σ ⋆ be the solution of (6). Then the difference σ ⋆ − ρ t tr is smaller than O(ε √ rd). This holds with probability of failure at most 1/λ 2 plus the probability of failure in Theorem 1. The proof combines ideas from Ref. [12] with our argument above [19]. The main difference from the noise-free case is that, instead of using tr Y ∆ = 0, we must now work with | tr Y ∆| ≤ 2 Y 2 δ. With this estimate, Observation 1 follows from the noise-free proof, together with some elementary calculations (see EPAPS). We remark that the above bound is likely to be quite loose; based on related work involving the "restricted isometry property," we conjecture that the robustness to noise is actually substantially stronger than what is shown here [13] . Certified tomography of almost pure states. The preceding results require an a priori promise: that the true state ρ t is δ 1 -close to a rank-r state. However, when performing tomography of an unknown state, neither r nor δ 1 are known beforehand. There are a few solutions to this quandary. First, r and δ 1 may be estimated from other physical parameters of the system, such as the strength of the local noise [20]. Another approach is to estimate r and δ 1 from the same data that is used to reconstruct the state. When r = 1, this approach is particularly effective, in entirely assumption-free tomography: one can estimate δ 1 , using only O(d) Pauli expectation values. This is because δ 1 is related to the purity Tr ρ 2 , which has a simple closed-form expression in terms of Pauli expectation values. See EPAPS for details. We get: Observation 2 (Certified tomography) Assume that the unknown physical state is close to being pure. Then one can find a certificate for that assumption, and reconstruct the state with explicit guarantees on the reconstruction error, from O(cd log 2 d) Pauli expectation values. The probability of failure is exponentially small in c. Finally, when the state is approximately low-rank but not nearly pure (r > 1), one may perform tomography using different numbers of random Pauli expectation values m. When m is larger than necessary (corresponding to an over-estimate of r), we are guaranteed to find the correct density matrix. When m is too small, we find empirically that the algorithms for reconstructing the density matrix (i.e., solving the convex program (1)) simply fail to converge. A hybrid approach to matrix recovery. Here we describe a variant of our tomography method that makes the classical post-processing step (i.e., solving the convex program (1) to reconstruct the density matrix) faster. This method also uses random Pauli measurements, but they are chosen in a structured way. Any Pauli matrix is of the form w(u, v) = n . We choose a random subset S ⊂ {0, 1} n of size O(r polylog(d)), and then for all u ∈ S and v ∈ {0, 1} n , measure the Pauli matrix w(u, v). We call this the "hybrid method" because it is equivalent to a certain structured matrix completion problem. This fact implies that certain key computations in solving the convex program (1) can be implemented in time . However, the hybrid method is not covered by the strong theoretical guarantees shown earlier, though it does give accurate results in practice. For a more complete discussion, see the EPAPS supplement. Numerical results. We numerically simulated both the random Pauli and hybrid approaches discussed above. For both approaches, we used singular value thresholding (SVT) [14] . Instead of directly solving Eq. (6), SVT minimizes τ σ tr + σ 2 2 /2 subject to | tr(σ − ω)w(A i )| ≤ δ, which is a good proxy to Eq. (6) when τ dominates the second term; the programs are equivalent in the limit τ → ∞ (provided Eq. (6) has a unique solution) [14] . Estimating the second term for typical states suggests choosing 2τ r ≫ 1; we use τ = 5. To simulate tomography, we chose a random state from the Haar measure on a d × r dimensional system and traced out the r-dimensional ancilla, then applied depolarizing noise of strength γ. We sampled expectation values associated with randomly chosen operators as above, and added additional statistical noise (respecting Hermiticity) which was i.i.d. Gaussian with variance σ 2 and mean zero. We used SVT and quantified the quality of the reconstruction by the fidelity and the trace distance for various values of m, each averaged over 5 simulations. This dependence is shown in Fig. 1 . The reconstruction is remarkably high fidelity, despite severe undersampling and corruption by both depolarizing and statistical noise [26]. Using the hybrid method with 8 qubits on a rank 3 state plus γ = 5% depolarizing, and statistical noise strength σd = 0.1, we typically achieve 95% fidelity reconstructions in under 10 seconds on a modest laptop with 2 GB of RAM and a 2.2 GHz dual-core processor using MATLAB -even though 90% of the matrix elements remain unsampled. Increasing the number of samples only improves our accuracy and speed, so long as sparsity is maintained. Using truly randomly chosen Pauli observables (instead of the hybrid method) slightly increases the processing time due to the dense matrix multiplications involved: in our setup about one minute. However, this method achieves even better performance with respect to errors, as seen in Fig. 1 . The simulations above show that our method work for generic low rank states. Lastly, we demonstrate the functioning of the approach in the experimental context of the state ρ found in the 8 ion experiment of Ref. [2] . To exemplify the above results, we simulated physical measurements by sampling from the probability distribution computed using the Born rule applied to the reconstructed state ρ. This state is approximately low-rank, with 99% of the weight concentrated on the first 11 eigenvectors. The standard deviation per observable was 3/d. Fewer than 30% of all Pauli matrices were chosen randomly. From this information, a rank = 3 approximation σ with fidelity of 90.5% with respect to ρ was found in about 3 minutes on the aforementioned laptop. Discussion. We have presented new methods for low-rank quantum state tomography, which require only O(rd log 2 (d)) measurements, where r is the rank of the unknown density matrix and d is the Hilbert space dimension. Our methods are based on and further develop the new paradigm of compressed sensing, and in particular, matrix completion [10, 11] . We use measurements that are experimentally feasible, together with very fast classical post-processing. The methods perform well in practice, and are also supported by theoretical guarantees. It would be interesting to further flesh out the trade off between the need for measurements that can be performed easily in an experiment and the need for sparse matrices during the classical post-processing step. It is the hope that this work stimulates such further investigations. Acknowledgments. We thank E. Candès and Y. Plan for useful discussions.
He also presents another special case of the results announced in REF : the reconstruction of low-rank matrices from randomly sampled matrix elements.
Quantum state tomography via compressed sensing
Phys. Rev. Lett. 105, 150401 (2010)
Physics Medicine