Efficient Measurement of Quantum Dynamics via Compressive Sensing

Compressed Quantum Process Tomography A. Shabani(1) , R. L. Kosut(2) , and H. Rabitz(1)

arXiv:0910.5498v1 [quant-ph] 28 Oct 2009

(1)

Department of Chemistry, Princeton University, Princeton, New Jersey 08544 (2) SC Solutions, Sunnyvale, CA 94085

The characterization of a decoherence process is among the central challenges in quantum physics. A major difficulty with current quantum process tomography methods is the enormous number of experiments needed to accomplish a tomography task. Here we present a highly efficient method for tomography of a quantum process that has a small number of significant elements. Our method is based on the compressed sensing techniques being used in information theory. In this new method, for a system with Hilbert space dimension n and a process matrix of dimension n2 × n2 with sparsity s, the required number of experimental configurations is O(s log n4 ). This heralds a logarithmic advantage in contrast to other methods of quantum process tomography. More specifically, for q-qubits with n = 2q , the scaling of resources is O(sq) — linear in the product of sparsity and number of qubits. PACS numbers:

Introduction.— The characterization of quantum system dynamics is a demanding task in quantum control and quantum information sciences. In particular, identification of the dynamics of a quantum system embedded in an environment is an important task, for the main obstacle in manipulating systems in the quantum regime is ambient disturbances. Among many theories for the dynamics of an open quantum system [1], a completely positive trace preserving (CPTP) linear map is a formulation that captures the essence of the decoherence process. In the past decade, different methods have been developed to identify the CPTP map, the process matrix. These methods are known as quantum process tomography (QPT). Common to all QPT techniques is the requirement of preparing the system in some known initial state, measuring an observable at a later fixed time interval, repeating this for a different states and measurement protocols, but always for the same time interval. The latter step can be done indirectly by estimating the state of the decohered system, employed in Standard QPT [2, 3] and Ancilla-Assisted QPT [4], or directly by measuring the expectation values of certain observables, used in Direct QPT [5]. A major problem for QPT is the enormous number of unknown parameters to be identified. Recall that a CPTP map has the general representation [2] 2

S(ρ) =

n X

χαβ Γα ρΓ†β ,

α,β=1

X

χαβ Γ†β Γα = In

(1)

α,β

where n is the dimension of the system Hilbert space and {Γα } is a complete set of matrix orthonormal basis, Tr(Γ†β Γα ) = δαβ . The CP condition implies that the process matrix, χ, is positive semidefinite; the TP condition is enforced by the above linear equality. The process matrix 2 2 χ ∈ Cn ×n is thus constrained to the convex set, X χ ≥ 0, χαβ Γ†β Γα = In (2) α,β

With these two constraints on χ the total number of parameters is n4 − n2 . In the context of quantum computation, it is

24q − 22q for q qubits, yielding an exponential growth in the number of parameters to be estimated, and thus also in the distinct experimental configurations (initial state preparation and final measurement). Therefore the task of performing QPT, on even a small system, is very challenging. Here we develop a new method for tomography of process matrices with sparse or a few significant elements, the latter is referred to here as almost sparse; this will be made precise in the sequel. We demonstrate that the number of required experimental configurations scales linearly with the number of qubits, that is, O(s log(n4 )) = O(sq) for a process matrix χ with s nonzero or significant elements. The sparsity of a process matrix is a relevant assumption in different physical systems. A sparse χ is observed in several QPT experiments, NMR [6], ion trap [7] and linear optics [8]. As discussed in these works, the sparsity of the dynamics is related to there being a small number of sources of errors. Intuitively, a decoherence process originating from a few noise sub-dynamics yields a sparse process matrix. In some cases, the process can be well approximated by a sparse χ matrix, for instance the secular approximation in NMR [9]. In addition, as posited in [18], for an initially well designed system which is close to a desired unitary, the process matrix in a basis corresponding to the ideal unitary is almost sparse. Our highly efficient algorithm for QPT, termed compressed quantum process tomography (CQPT), is based on the recently developed idea of compressed sensing/sampling (CS) in information theory. CS provides techniques to recover a sparse vector from a relatively small number of measurements. Some introductory reviews of CS can be found in [10]. Here we briefly review the mathematical foundation of this method before presenting the CQPT formulation. Logarithmically efficient encoding.— With proper techniques QPT can be accomplished with drastically reduced number of experimental configurations. This requires an efficient way to encode the information of the sparse vector into a few measurement outputs. Existence of such an encoding was originally proved by Johnson and Lindenstrauss (JL) [11]. To

2 embed a given set of p points in RN into a space of lower dimension Rm , while approximately preserving the relative distances between any two points. JL have shown the possibility of such an embedding if m > O(log(p)/ǫ2 ), where ǫ ∈ (0, 1) is the deviation of the relative distance between points. Since the early existence proof by JL, several constructive proofs have been given, that to our knowledge, all consider linear maps for embedding. In the context of CS, we are interested in extracting information from a vector (signal) x ∈ Rn by m linear measurements given by y = Φx where Φ is a m × n matrix. We say Φ satisfies the restricted isometry property (RIP) if it embeds the set of s-sparse vectors xs into a smaller space y ∈ Rm in the JL sense. A s-sparse vector is one with at most s non-zero entries. More rigorously, a matrix Φ obeys RIP of order s if there exists a constant δs ∈ (0, 1) such that [12] 2

2

2

(1 − δs ) kxT kℓ2 ≤ kΦxT kℓ2 ≤ (1 + δs ) kxT kℓ2

(3)

for all T -sparse vectors with T ≤ s. Convex decoding.— Is there a decoder mapping ∆ : Rm → n R that recovers the vector x from the extracted noisy information ky − Φxkℓ2 ≤ ǫ without knowing the sparsity pattern? The answer to this question is affirmative and related to the bounds on Gelfand widths [13, 14] which is beyond the scope of this paper to address. In the following we only mention the relevant results from CS. Define the decoder to be [15] ∆(y) = arg min kxkℓ1 , subject to ky − Φxkℓ2 ≤ ǫ.

(4)

x

Suppose Φ obeys √ RIP (3) of the order s. Cand´es [16] shows that if δ2s < 2 − 1, then ∃ constants C, D, kx − ∆(y)kℓ2

C ≤ √ kx − xs kℓ1 + Dǫ s

(5)

where xs is the best s-sparse approximation of x. The former tells us that the decoder is an ℓ1 -norm minimization of the unknown vector under measurement linear constraints with ℓ2 noise, thereby producing a convex optimization which can be implemented very efficiently. The latter determines the performance of this optimization recipe. For a physical signal, ideal sparsity, having a large number of elements equal to zero, is not a realistic circumstance, thus we employ the notion of almost sparsity, meaning that most of the parameters are below a threshold that can be approximated to zero. Inequality (5) guarantees that the signal x is recoverable with the same accuracy level as such an approximation, and with an ℓ2 -measurment noise level ǫ. Compressed QPT.— CS theory promotes the possibility of estimating the process matrix χ with a number of experiment configurations much less than n4 − n2 . In each experiment we initialize the system in some known state ρ and then an observable M is measured after a certain decoherence time. Suppose the system can be prepared in any state of the set {ρ1 , ..., ρk }, and further assume that we are able to measure a limited set of observables {M1 , ..., Mh }.

For a pair of (ρ, M ), the measurement output hM i will be Pn2 † yM,ρ = Tr(S(ρ)M ) = α,β=1 χαβ Tr(Γα ρΓβ M ) If the experiment is repeated for a total of m different configurations of (ρi , Mi ), i = 1, ..., m, the relation between measurement data y = (1/m)(yM1 ,ρ1 , . . . , yMm ,ρm )T and the true process matrix χ0 can be represented by a linear relation 4 y = Φvec(χ0 ). where vec(χ0 ) ∈ Cn is the vectorized form of the process matrix, and Φ is a m × n4 matrix of coefficients Tr(Γα ρΓ†β M )/m. (A factor 1/m is included for simplicity of the proofs.) As in (4), we assume ℓ2 -noisy measurements: ky − Φvec(χ0 )kℓ2 ≤ ǫ

(6)

From CS theory we prove the main result. Theorem 1 Let χ0 (s) denote an approximation to χ0 obtained by selecting the s largest elements of χ0 and setting the remaining elements to zero. If the measurement matrix 4 Φ ∈ Cm×n is drawn randomly from a probability distribution that satisfies the concentration inequality in (13) below, then ∃ constants C1 , C2 , C3 , C4 > 0 such that the solution χ⋆ to the convex optimization problem, minimize kvec(χ)kℓ1 subject to ky − Φvec(χ)kℓ2 ≤ ǫ, χ is CPTP (2)

(7)

satisfies, C1 kvec(χ⋆ − χ0 )kℓ2 ≤ √ kvec(χ0 (s) − χ0 )kℓ1 + C2 ǫ (8) s with probability ≥ 1 − 2e−mC3 provided that, m ≥ C4 s log(n4 /s)

(9)

If for a q-qubit QPT with n = (2q )2 the conditions of the theorem are satisfied, then the number of measurements m scales linearly with sq, specifically, m ≥ C4 s(8q log 2 − log s) = O(sq). If the measurements are noise free (ǫ = 0) and the process matrix is actually s-sparse (χ0 = χ0 (s)), then the right hand side of Eq.(8) is zero. Consequently, under the conditions of the theorem, perfect recovery is attained, i.e., χ⋆ = χ0 . Otherwise the solution tends to the best s-sparse approximation of the almost sparse process matrix. Unlike the standard CS form (4), CQPT (7) is constrained by the CPTP conditions (2). Since the true process matrix χ0 and the solution χ⋆ both obey the physical constraints of CPTP, we find that no changes are needed in the proofs presented in [16]. Proof. — The results in [13, 16] show that the RIP condition (3) holds with high likelihood if Φ satisfies (13). Suppose that for each experimental configuration we initialize the system randomly in a state ρ ∈ {ρi:1...k1 } and then measure an observable M randomly chosen from {Mj:1...k2 }. The corresponding matrix Φ then has independent rows with correlated elements of each row since they are functions of the same M and ρ. This differs from the randomization schemes in the CS

3 literature where all elements of Φ can be independently and randomly selected. Although Φ is a random matrix, because it is constructed from quantum states and observables of a finite dimensional system, it is bounded, i.e., 2

2

2

ℓ kxkℓ2 ≤ E kΦxkℓ2 ≤ u kxkℓ2 , ∀x ∈ RN

(10)

for some u, ℓ ∈ R. Here E denotes the expectation value. Next we apply, Hoeffding’s concentration inequality Let v1 , ..., vm be independent bounded random variables such that vi falls P in the interval [ai , bi ] with probability one. Then for S = i vi and any t > 0 we have 2

/

P

i (bi −ai )

2

2

/

P

i (bi −ai )

2

Pr {S − E(S) ≥ t} ≤ e−2t

Pr {S − E(S) ≤ −t} ≤ e−2t

(11)

In our problem vi is kφi xk2ℓ2 and S = kΦxk2ℓ2 . From the above inequalities and the relation (10) we find, ∀t+ , t− > 0 and ∀x n o 2 Pr S − u kxkℓ2 ≥ t+ ≤ Pr {S − E(S) ≥ t+ } 2

≤ e−2t+ /

P

i (bi −ai )

2

o n Pr S − l kxk2ℓ2 ≤ −t− ≤ Pr {S − E(S) ≤ −t− } 2

≤ e−2t− /

P

i (bi −ai )

2

(12)

A worst-case estimation of m[ai , bi ]/ kxk2ℓ2 is found from, [min(φαβ (ρ, M )), max(φαβ (ρ, M ))]. The choice of t+ = 2 2 (δ + 1 − u) kxkℓ2 and t− = (ℓ − 1 + δ) kxkℓ2 in the above inequalities yields o n 2 2 2 2 2 Pr | kΦxkℓ2 − kxklm | ≥ δ kxkℓ2 ≤ 2e−2m(δ+ǫ) /(wu −wl ) 2

(13)

with ǫ = min{1 − u, ℓ − 1}. To insure that t+ , t− > 0, we need 1 − δ < ℓ ≤ u < 1 + δ. Since the observable M can be scaled by any real number, a sufficient condition is u/ℓ < (1 + δ)/(1 − δ)

(14)

Based on [13, Thm.5.2] and [16, Thm.1.2], since the random matrix Φ satisfies the concentration inequality (13) for a given 0 < δ < 1, the results in (7) follow. This completes the proof. Simulation results.— QPT is performed by solving (7) with noise-free measurments (ǫ = 0) where the system is designed to be a quantum Fourier transform (QFT) with unitary representation Uqft ∈ Cn×n . This is corrupted by interaction with an unknown nE -dimensional environment via the total con˜ with the evolution stant Hamiltonian, H = InE ⊗ Hqft + γ H time interval normalized to one. The interaction Hamiltonian ˜ is randomly selected, then normalized to kHk ˜ = 1; γ is H true then the interaction magnitude. The system S (1) and associated process matrix χtrue are extracted via the partial trace over the environment. The values of γ selected for the simulations are in the range [0, 1.25]. These induce various fidelity values as will be indicated.

We present results for both a 2-qubit (240 parameters) and 3-qubit (4032 parameters) system. The estimates from (7) are obtained in the singular value decomposition (SVD)-basis [Γα in (1)] of the ideal QFT where u/l ≈ 1.3 ensuring δ ≈ 0.13. The process matrix of the ideal unitary in this basis is maximally sparse with the single non-zero 11-element equal to n [18]. The environmental interactions make the process matrix almost sparse as defined in (7). For example, Fig. 1 shows 3D bar plots of the real and imaginary elements of the process matrices in the SVD-basis of the ideal QFT. Column 1 shows the true system with γ = 1.25 which gives a low channel fidelity with respect to the ideal unitary of fchn (Uqft , S true ) = 0.732 [21]. Column 2 shows an estimate using m = 36 configurations as will be described below. Both the true process matrix and the estimate exhibit the expected large 11-element. The remaining elements are much smaller by comparison. To evaluate the estimation quality we use the worst-case pure state fidelity between the true system S true and the estimated system S est [22]. In this case, fwc (S true , S est ) = 0.90 (see 4th row of Fig. 3, p = 6). Considering the paucity of resources and relatively low value of channel fidelity, this is a surprisingly good estimate. This result can be understood by examining Fig. 2, which shows the absolute value, sorted by magnitude, of the true process matrix elements relative to the 11-element. The contributions of the approximately 200 out of 256 elements with magnitudes below the 0.01 level have very little effect. This is the typical pattern observed of an almost sparse process matrix expressed in the SVD basis of the desired unitary. As we have observed from the simulation study to follow, the resulting estimates behave like the best s-sparse reconstruction as predicted by the theory (7). This is why we get a relatively good (0.90) worst-case fidelity. 4 To form the measurement matrix Φ ∈ Cm×n , we randomly generated p inputs of the form |ψ1 i|ψ2 i and randomly selected p single qubit observables from the Pauli set {XI, Y I, ZI, IX, IY, IZ}, thus, m = p2 . For 2-qubits (n = 4, 240 process parameters) we solved (7) for all 12 combinations of m ∈ {4, 16, 36} and 4 instances of random environment interactions giving rise to the following channel fidelities: fchn (Uqft , S true ) ∈ {1.0, 0.988, 0.895, 0.736}. Fig. 3 shows the results. The black bars show fchn (Uqft , S true ), the true channel fidelities between the ideal unitary Uqft and the true system which is exactly (χsvd )11 /n. The gray bars in the left plot are the estimated channel fidelity fchn (Uqft , S est ) = (χ⋆svd )11 /n where χ⋆svd is an estimate of χsvd obtained from the solution of (7) in the SVD basis of the ideal unitary. Each white bar shows the worst-case fidelity fwc (S true , S est ). Remarkably, even for the low channel fidelity of 0.732 (4th row), the worst-case fidelities never fall below 0.90. In other words, the process matrix estimates are reasonably good for very few experiment resources. This is borne out even more for the 3-qubit example, particularly in comparison with the number of variables to be estimated in the process matrix: 4,032. Fig. 4 shows 3 rows corresponding to all the 12 combinations of channel fidelities

4 est

χsvd

2

0

0

0

0

−0.2

−0.2

Imag

0.2

0.2

fchn = 0.988

2

fchn = 0.952

4

fchn = 0.895

4

Fidelities vs. Configurations

fchn = 0.732

Real

χtrue svd

1 0.5 0 1 0.5 0 1 0.5 0 1 0.5 0

4

16 Configurations (m)

36

16×16

FIG. 3: Ideal: 2-qubit QFT (χ ∈ C ). Estimates from (7) in SVD basis. Black bars: true compared to ideal fchn (S true , Uqft ). Gray bars: estimate compared to ideal fchn (S est , Uqft ). White bars: estimate compared to true fwc (S est , S true ).

FIG. 1: Ideal: 2-qubit QFT (χ ∈ C16×16 ). Real and imaginary process matrix elements for channel fidelity fchn (Uqft , S true ) = 0.732 (γ = 1.25). Column 1 True process matrix, Column 2 Estimated process matrix from m = 36 input/observable configurations for worst-case fidelity fwc (S true , S est ) = 0.90.

Fidelities vs. Configurations 1 fchn = 0.988

0

−2

0.5 0

10

1 −4

10

0

10

1

2

10

10

3

10

fchn = 0.896

Magnitude

10

Sorted Index

0

0.5 0

fchn (Uqft , S est ) ∈ {0.988, 0.896, 0.734} and configurations m ∈ {4, 16, 64, 256, 1024}. In all the rows the accuracy considerably improves for m = 256. In the second row, for example, the estimated channel fidelity (gray bar) is 0.896 which to three digits is exactly the true value (black bar) with a high worst-case fidelity (white bar) of 0.964; ; the latter of course is unknown except in simulations. Considering the very few relative measurements, this is clearly producing a useful process matrix estimate. Discussion and Conclusion.— In this letter, we introduced the scheme of compressed quantum process tomography (CQPT) for efficient estimation of sparse or almost process matrices. The logarithmic advantage of the required physical resources emanate in the theory from the randomness in system initialization and the final measurement. In practice, as long as the available set in the laboratory is sufficiently large, then the results will hold. Roughly, if the available resources engender a Φ matrix whose m rows are approximately independent, then there is an s-sparse approximation

2

1 fchn = 0.734

FIG. 2: Absolute values of the 256 process matrix elements of 16×16 χtrue sorted by relative magnitude (with respect to the svd ∈ C 11-element) for fchn (Uqft , S true ) = 0.732.

0.5

4

16

64 Configurations (m)

256

1024

FIG. 4: Ideal: 3-qubit QFT (χ ∈ C64×64 ). Estimates from (7) in SVD basis. Black bars: true compared to ideal fchn (S true , Uqft ). Gray bars: estimate compared to ideal fchn (S est , Uqft ). White bars: estimate compared to true fwc (S est , S true ).

which will be obtained by solving (7), where s, m are related via (8)-(9). Current CS theory, including the theory presented here, can be viewed as providing the “scaffold” upon which rest the many reported dramatic results, e.g., the many articles, tutorials, testimonials, and software on the compressive sensing website [10]. Applications to quantum systems are just beginning, e.g., for QPT as developed here and in [18], for low rank quantum state estimation via matrix completion in [19], and for ghost-imaging in [20]. For quantum system identification in general, compressive sensing compels the necessity for further theory development

5 to tighten the bounds on scaling laws, to pursue highly efficient convex computational algorithms, and expand the horizon, if possible, to include Hamiltonian parameter estimation and quantum metrology, Acknowledgement. — RLK and HR were supported by DARPA Grant FA9550-09-1-0710. We thank Sina Jafarpour for useful discussions.

[1] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, Oxford, 2002). [2] M.A. Nielsen and I.L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, UK, 2000). [3] J. F. Poyatos et al., Phys. Rev. Lett. 78, 390 (1997);I. L. Chuang and M. A. Nielsen, J. Mod. Opt. 44, 2455 (1997). [4] G. M. DAriano and P. Lo Presti, Phys. Rev. Lett. 86, 4195 (2001); J. B. Altepeter et al., Phys. Rev. Lett. 90, 193601 (2003). [5] M. Mohseni and D. A. Lidar, Phys. Rev. Lett. 90, 193601 (2006) [6] A. M. Childs et al., Phys. Rev. A 64, 012314 (2001);Boulant:03 N. Boulant et al., Phys. Rev. A 67, 042322 (2003);Weinstein:04 Y. S. Weinstein et al., J. Chem. Phys. 121, 6117 (2004). [7] M. Riebe et al., Phys. Rev. Lett. 97, 220407 (2006). [8] M. W. Mitchell et al., Phys. Rev. Lett. 91, 120402 (2003);O’Brien:04 J. L. O’Brien et al., Phys. Rev. Lett. 93, 080502 (2004);Nambu:05 Y. Nambu et al., Phys. Rev. Lett. 94,

010404 (2005). [9] R. r. Ernst et al., Principles of Nuclear Magnetics Resonance in One and Two Dimensions (Oxford University Press, Oxford, 1994). [10] E. J. Cand´es and M. B. Wakin, IEEE Sig. Proc. Mag. Mar., 21 (2008);E. J. Cands and T. Tao, IEEE Inform. Theory. News. Dec., 14 (2008);D. L. Donoho, IEEE Inform. Theory. News. Dec., 18 (2008); Compressive Sensing Website http://dsp.rice.edu/cs. [11] W. B. Johnson and J. Lindenstrauss, Contemporary Mathematics 26 189 (1984). √ P [12] For a vector x, kxkℓ2 = x† x and kxkℓ1 = i |xi |. [13] R. Baraniuk et al., Constr. Approx. 28 253 (2008). [14] A. Cohen et al., J. Amer. Math. Soc. 22 211 (2009). [15] E. J. Cand´es and T. Tao, IEEE Trans. Inform. Theory, 51 4203 (2004). [16] E. J. Cand´es, Compte Rendus de l’Academie des Sciences, Paris, Serie I, 346 589 (2008). [17] E. J. Cand´es textitet al., Comm. Pure Appl. Math, 59 1207 (2006). [18] R. L. Kosut, arXiv:0812.4323v1[quant-ph], (2008). [19] D. Gross et al. arXiv:0909.3304 (2009). [20] W. Gong and S. Han arXiv:0910.4823v1[quant-ph], (2009). [21] In the SVD-basis of any unitary P U , (χsvd )11 = nfchn (U, S) with fchn (U, S) = (1/n2 ) µ |Tr(U † Sµ )|2 [18]. [22] Worst-case fidelity for pure state inputs between two systems S1 , S2 is fwc (S1 , S2 ) = ‚ ”2 “‚p p ‚ ‚ † † minψ† ψ=1 ‚ S1 (ψψ ) S2 (ψψ )‚ where k·k∗ is ∗ the nuclear norm, the sum of the singular values [2].