SYSTEMS AND METHODS FOR COMPRESSED SENSING MEASUREMENT OF LONG-RANGE CORRELATED NOISE

Abstract

A method for detecting a two-qubit correlated dephasing error includes accessing a signal of a quantum system, where the quantum system includes a plurality of qubits. Every qubit has a nonzero rate of dephasing and some qubits have a nonzero rate of correlated dephasing. The signal further includes information about a matrix that includes diagonal elements and off-diagonal elements. The off-diagonal elements of the matrix are 2 s-sparse. The method further includes performing randomized measurements of the off-diagonal elements of the matrix and recovering the matrix based on a direct measurement of the diagonal elements of the matrix.

Claims

1. A computer-implemented method for detecting a two-qubit correlated dephasing error, the method comprising: accessing a signal of a quantum system, the quantum system includes a plurality of qubits, wherein every qubit has a nonzero rate of dephasing, wherein some qubits have a nonzero rate of correlated dephasing, wherein the signal includes information about a matrix, the matrix including diagonal elements and off-diagonal elements, and wherein the off-diagonal elements of the matrix are 2s-sparse; performing randomized measurements of the off-diagonal elements of the matrix; and recovering the matrix based on a direct measurement of the diagonal elements of the matrix.

2. The computer-implemented method of claim 1, further comprising estimating a two-qubit correlated dephasing error based on the recovered matrix.

3. The computer-implemented method of claim 1, wherein performing the randomized measurements includes preparing entangled Greenberger-Horne-Zeilinger states (GHZ states) on random subsets of the plurality of qubits.

4. The computer-implemented method of claim 1, wherein the randomized measurements are based on noise spectroscopy and quantum sensing.

5. The computer-implemented method of claim 1, wherein the recovered matrix includes a restricted isometry property.

6. The computer-implemented method of claim 1, further comprising estimating a vector of decay rates that are dependent on the matrix.

7. The computer-implemented method of claim 6, wherein the recovered matrix is further based on the estimated decay rate.

8. The computer-implemented method of claim 6, further comprising estimating a relaxation time of the quantum system based on the estimated vector of decay rates.

9. The computer-implemented method of claim 6, further comprising estimating a decoherence time of the quantum system based on the estimated vector of decay rates.

10. The computer-implemented method of claim 1, further comprising detecting a long-range correlated dephasing error based on the recovered matrix.

11. A computer-implemented method for detecting a two-qubit correlated dephasing error, the method comprising: accessing a signal of a quantum system, the quantum system includes a plurality of qubits, wherein every qubit of the plurality of qubits has a nonzero rate of dephasing and wherein some qubits have a nonzero rate of correlated dephasing; dephasing entangled states of the plurality of qubits based on performing Ramsey spectroscopy using entangled states of random subsets of qubits of the plurality of qubits; measuring a linear function of a correlation matrix, where the correlation matrix corresponds to correlated Markovian dephasing between pairs of qubits of the plurality of qubits; generating a first vector and a second vector, wherein a plurality of elements of the first vector and of the second vector are randomly chosen; and estimating a decay rate based on the first vector and the second vector.

12. The computer-implemented method of claim 11, further comprising generating a restricted isometry property-based recovery matrix, wherein the recovery matrix is further based on the estimated decay rate.

13. The computer-implemented method of claim 12, further comprising detecting a long-range correlated dephasing error based on the recovery matrix.

14. The computer-implemented method of claim 11, further comprising estimating a relaxation time of the quantum system based on the estimated decay rate.

15. The computer-implemented method of claim 11, further comprising estimating a decoherence time of the quantum system based on the estimated decay rate.

16. A system for detecting a two-qubit correlated dephasing error, the system comprising: a processor; and a memory, including instructions stored thereon, which when executed by the processor, cause the system to: access a signal of a quantum system, the signal includes a plurality of qubits, wherein every qubit has a nonzero dephasing rate; and generate a matrix C in R.sup.n×n, wherein its off-diagonal part is 2s sparse and wherein the matrix C is based on the accessed signal.

17. The system of claim 16, wherein the instructions, when executed by the processor, further cause the system to: perform randomized measurements based on: estimating (b−a).sup.TC(b−a) for any a≠b in {0,1}n; choosing a sequence of random vectors r(1), . . . , r(m)in {1,0,−1}n; and estimating ΦC=(r.sup.(1){circumflex over ( )}TCr.sup.(1)), . . . , (r.sup.(m){circumflex over ( )}TCr.sup.(m), where Φ:R.sup.n×n.fwdarw.R.sup.m.

18. The system of claim 17, wherein the instructions, when executed by the processor, further cause the system to: estimate a decay rate Γ.sub.ab based on any a≠b in {0,1}.sup.n; and estimate a decoherence time of the quantum system based on the decay rate Γ.sub.ab.

19. The system of claim 18, wherein the instructions, when executed by the processor, further cause the system to: recover Matrix C based on: measuring a plurality of diagonal elements of matrix C directly; and determining an off-diagonal element of matrix C based on using .sub.1-regularized least-squares regression.

20. The system of claim 19, wherein the instructions, when executed by the processor, further cause the system to: detect a long-range correlated dephasing error based on the recovered matrix C.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0027] A better understanding of the features and advantages of the present disclosure will be obtained by reference to the following detailed description that sets forth illustrative embodiments, in which the principles of the present disclosure are utilized, and the accompanying drawings of which:

[0028] FIG. 1 is a diagram of an exemplary noisy intermediate-scale quantum information (NISQ) processor, in accordance with examples of the present disclosure;

[0029] FIG. 2 is a high level flow diagram of the computer-implemented method for detecting a two-qubit correlated dephasing error, in accordance with examples of the present disclosure;

[0030] FIG. 3 is a correlation graph of a noise model of the quantum system of FIG. 1, in accordance with examples of the present disclosure;

[0031] FIG. 4 illustrates a flow diagram of the computer-implemented method 600 for determining a two-qubit correlated dephasing error for the quantum system 100, such as the quantum system of FIG. 1, in accordance with examples of the present disclosure;

[0032] FIG. 5 is a graph that illustrates single qubit Ramsey spectroscopy, in accordance with examples of the present disclosure;

[0033] FIG. 6 is a correlation matrix C corresponding to the correlation graph of FIG. 3, in accordance with examples of the present disclosure;

[0034] FIG. 7A is a table that illustrates sample complexity of different methods for reconstructing the correlation matrix C, in accordance with examples of the present disclosure;

[0035] FIG. 7B is a table that illustrates sample complexity of different methods for reconstructing the correlation matrix C, in accordance with examples of the present disclosure;

[0036] FIGS. 8A-8C are graphs that illustrate scaling of the reconstruction error under various circumstances, in accordance with examples of the present disclosure;

[0037] FIG. 9A is a graph that illustrates trajectories of a random walk used in estimating an evolution time, in accordance with examples of the present disclosure;

[0038] FIG. 9B is a graph that illustrates the accuracy of a final estimate of the evolution time, in accordance with examples of the present disclosure;

[0039] FIGS. 10A-10C are graphs that illustrate an exponential decay of {tilde over (P)}.sub.ab(t) when there are state-preparation-and-measurement (SPAM) errors for a 3-qubit Greenberger-Horne-Zeilinger (GHZ) state, in accordance with examples of the present disclosure; and

[0040] FIG. 11 illustrates a flow diagram of the computer-implemented method for determining a two-qubit correlated dephasing error for the quantum system, such as the quantum system of FIG. 1, in accordance with examples of the present disclosure.

DETAILED DESCRIPTION

[0041] The present disclosure relates generally to the field of noisy intermediate-scale quantum information (NISQ) processors. More specifically, an aspect of the present disclosure provides detecting two-qubit correlated dephasing errors in NISQ devices.

[0042] Embodiments of the present disclosure are described in detail with reference to the drawings wherein like reference numerals identify similar or identical elements.

[0043] Although the present disclosure will be described in terms of specific examples, it will be readily apparent to those skilled in this art that various modifications, rearrangements, and substitutions may be made without departing from the spirit of the present disclosure. The scope of the present disclosure is defined by the claims appended hereto.

[0044] For purposes of promoting an understanding of the principles of the present disclosure, reference will now be made to exemplary embodiments illustrated in the drawings, and specific language will be used to describe the same. It will nevertheless be understood that no limitation of the scope of the present disclosure is thereby intended. Any alterations and further modifications of the novel features illustrated herein, and any additional applications of the principles of the present disclosure as illustrated herein, which would occur to one skilled in the relevant art and having possession of this disclosure, are to be considered within the scope of the present disclosure.

[0045] The notation x∝y means that x and y are proportional, i.e., x=cy where c is a “universal” constant. A universal constant is a quantity whose value is fixed once and for all, and does not depend on any other variable. Polylogarithmic factors are written in a compact way as follows: log.sup.cn=(log n).sup.c. Big-O notation, such as O(n.sup.3), denotes an asymptotic upper bound. Matrices are denoted by capital letters. For a vector v, ∥v∥.sub.p denotes the .sub.p norm (in cases where p=1,2, ∞). For a matrix M, ∥M∥.sub.F=(Σ.sub.jk|M.sub.jk|.sup.2).sup.1/2 denotes the Frobenius norm (i.e., the Schatten 2-norm, or the .sub.2 norm of the vector containing the entries of M), ∥M∥ denotes the operator norm (i.e., the Schatten ∞-norm), ∥M∥.sub.tr denotes the trace norm (i.e., the Schatten 1-norm, or the nuclear norm), and ∥M=Σ.sub.jk|M.sub.jk| denotes the norm of the vector containing the entries of M. Given an n×n matrix M, let uvec(M)=(M.sub.ij).sub.i<j denote the vector (of dimension n(n−1)/2) containing the entries in the upper-triangular part of M, excluding the diagonal. Statistical estimators are written with a hat superscript, e.g., {circumflex over (Γ)} is a random variable that represents an estimator for some unknown quantity Γ. For a random variable {circumflex over (Γ)}, ∥{circumflex over (Γ)}∥.sub.ψ.sub.2 denotes the sub-Gaussian norm, and ∥{circumflex over (Γ)}∥.sub.ψ.sub.1 denotes the subexponential norm.

[0046] The development of noisy intermediate-scale quantum information processors (NISQ devices) has the potential to advance many areas of computational science. An important problem is the characterization of noise processes in these devices, in order to improve their performance (via calibration and error correction), and to ensure correct interpretation of the results. One challenge is to characterize all of the noise processes that are likely to occur in practice, using some experimental procedure that is efficient and can scale up to large numbers of qubits.

[0047] Compressed sensing offers one approach to solving this problem, where one uses specially-designed measurements (and classical post-processing) to learn an unknown signal that has some prescribed structure. For example, the unknown signal can be a sparse vector or a low-rank matrix, the measurements can consist of random projections sampled from various distributions, and the classical post-processing can consist of solving a convex optimization problem (e.g., minimizing the .sub.1 or trace norm), using efficient algorithms. This approach has been used in several previous works on quantum state and process tomography, and estimation of Hamiltonians and Lindbladians.

[0048] From a theoretical perspective, one of the main challenges is to design measurements that have the mathematical properties needed for compressed sensing, and can be implemented efficiently on a quantum device. There has been a substantial amount of work in this area, which can be broadly classified into two approaches: “sparsity-based” and “low-rank” compressed sensing. For compressed sensing of “low-rank” objects (e.g., low-rank density matrices and quantum processes), there seems to be a few natural choices for measurements, including random Pauli measurements, and fidelities with random Clifford operations. For compressed sensing of “sparse” objects (e.g., sparse Hamiltonians, Lindbladians, or Pauli channels), however, the situation is more complicated, as a larger number of different measurement schemes and classical post-processing methods have been proposed, and the optimal type of measurement seems to depend on the situation at hand. This complicated state of affairs can be explained in part because “sparsity” occurs in a wider variety of situations than “low-rankness.”

[0049] Long-range correlated errors can severely impact the performance of NISQ (noisy intermediate-scale quantum) devices, and fault-tolerant quantum computation. Characterizing these errors is important for improving the performance of these devices, via calibration and error correction, and to ensure correct interpretation of the results. A compressed sensing method for detecting two-qubit correlated dephasing errors may be used for characterizing these errors, assuming only that the correlations are sparse (i.e., at most s pairs of qubits have correlated errors, where s<<n(n−1)/2, and n is the total number of qubits). In particular, the disclosed method can detect long-range correlations between any two qubits in the quantum system (i.e., the correlations are not restricted to be geometrically local).

[0050] The disclosed method is highly scalable, because the disclosed method requires as few as m∝s log n measurement settings (where ∝ denotes proportionality), and efficient classical post-processing based on convex optimization. In addition, when m∝s log.sup.4n, the disclosed method is highly robust to noise, and has sample complexity O(max(n, s).sup.2 log.sup.4(n)), which can be compared to conventional methods that have sample complexity O(n.sup.3). Thus, the disclosed method is advantageous when the correlations are sufficiently sparse, that is, when s≤O(n.sup.3/2/log.sup.2n). The disclosed method also performs well in numerical simulations on small quantum system sizes, and has some resistance to state-preparation-and-measurement (SPAM) errors. A feature of the disclosed method is a new type of compressed sensing measurement, which works by preparing entangled Greenberger-Horne-Zeilinger states (GHZ states) on random subsets of qubits, and measuring their decay rates with high precision.

[0051] Characterizing errors when building large-scale quantum devices is needed for calibrating, correcting the errors, and interpreting the result of the quantum devices correctly. For example, quantum simulations have only a limited simulation time. Metrology and quantum sensing may lead to reduced sensitivity, while error correction and fault tolerance require increased overhead.

[0052] Thus, there is a need to find methods for characterizing all the likely noise processes in quantum systems that are efficient and scalable.

[0053] Referring to FIG. 1, a diagram of an example quantum system 100, is shown. The quantum system 100 may include, for example, a NISQ (noisy intermediate-scale quantum) device (i.e., a NISQ processor). The term “noisy” refers to the fact that quantum processors are very sensitive to the environment and may lose their quantum state due to quantum decoherence (i.e., dephasing). The quantum system 100 includes qubits 102 which are linked by a coupler 104. The qubits are the quantum version of classic binary bits realized with a two-state quantum-mechanical system. The coupler 104, for example, may be implemented by connecting the two qubits to an intermediate electrical coupling circuit. The circuit might be a fixed element, such as a capacitor, or a controllable element. It is contemplated that the quantum system 100 may include any suitable number of qubits 102 and/or couplers 104.

[0054] The theory of “sparsity-based” quantum compressed sensing may be applied to a physical problem that is relevant to the development of quantum system 100 (e.g., a NISQ device), such as detecting long-range correlated dephasing errors. A simple model of correlated dephasing, is described by the following Markovian master equation:

[00001]$\frac{d\rho }{\mathrm{dt}}=ℒ\left(\rho \right)={.Math.}_{j,k=1}^n{c}_{\mathrm{jk}}\left(Z_k\rho Z_j-\frac{1}{2}\left\{Z_j{Z}_k,\rho \right\}\right).$

[0055] Here the quantum system 100 consists of n qubits, and Z.sub.j and Z.sub.k are Pauli σ.sub.z operators that act on the j'th and k'th qubits, respectively. The noise is then completely described by the correlation matrix C=(c.sub.jk)∈.sup.n×n. (See FIGS. 3 and 6.) Generalizations of this model with complex coefficients c.sub.jk, and an additional environment-induced Hamiltonian are considered.

[0056] Here, the diagonal elements c.sub.jj show the rates at which single qubits dephase, and the off-diagonal elements c.sub.jk show the rates at which pairs of qubits undergo correlated dephasing. Typically, the diagonal of C will be nonzero, while the off-diagonal part may be dense or sparse, depending on the degree of connectivity between the qubits and their environment.

[0057] This master equation describes a number of physically plausible scenarios, such as spin-½ particles coupling to a shared bosonic bath. It has also been studied as an example of how collective decoherence can affect physical implementations of quantum computation and quantum sensing.

[0058] This model of correlated dephasing is quite different from other models of crosstalk that are based on quantum circuits or Pauli channels. The disclosed model describes crosstalk that arises from the physical coupling of the qubits to their shared environment. The disclosed model has a different character from crosstalk that arises from performing imperfect two-qubit gates, or correlated errors that result when the physical noise processes are symmetrized by applying quantum operations such as Pauli twirls.

[0059] The disclosed model of correlated dephasing can be learned efficiently when the off-diagonal part of the correlation matrix C is sparse. It is assumed that C has at most s<<n(n−1)/2 nonzero entries above the diagonal, but these entries may be distributed arbitrarily; in particular, long-range correlated errors are allowed. (Note that physical constraints imply that C is positive semidefinite, hence C=C.sup.†)

[0060] The fact that the disclosed model enables long-range correlations is quite powerful, and is relevant to a number of interesting settings in quantum information science. These settings include experimental NISQ devices, which are often engineered to have long-range interactions, in order to perform quantum computations more efficiently; the execution of quantum circuits on distributed quantum computers, where long-range correlations are generated when qubits are moved or teleported from one location to another; and quantum sensor arrays, where a detection event at an arbitrary location (j,k) in the array may be registered as a pairwise correlation between a qubit that is coupled to row j and a qubit that is coupled to column k in the array.

[0061] Referring to FIG. 2 a high level flow diagram of a computer-implemented method 200 for detecting a two-qubit correlated dephasing error is shown. A quantum system 100 such as a NISQ device is observed under noisy evolution 606 (FIG. 4). Random measurements 210 (e.g., Ramsey spectroscopy, FIG. 5) are made and a sparse correlation matrix C 400 (FIG. 6) is characterized based on the random measurements 210. The correlation matrix C 400 corresponds to correlated Markovian dephasing between pairs of qubits.

[0062] Referring to FIG. 3 an illustration of a correlation graph of a noise model 300 of the quantum system 100 of FIG. 1 is shown. The qubits 302 experience correlated Markovian dephasing. Non-zero c.sub.ij are shown, indicating correlated noise affecting the pairs of qubits connected by the lines 304.

[0063] FIG. 4 illustrates a computer-implemented method 600 for determining a two-qubit correlated dephasing error for the quantum system 100, such as the quantum system of FIG. 1. Initially, a general form of Ramsey spectroscopy using entangled states on multiple qubits is described. A simple method for directly measuring the correlation matrix C 400, by performing spectroscopy on every pair of qubits is described. This simple method serves as a baseline for measuring the performance of the disclosed compressed sensing method. It is contemplated that the disclosed methods may be performed by a processor and memory. The memory may include instructions thereon, which are executed by the processor.

[0064] Initially, a signal is accessed from a quantum system 100. The quantum system 100 contains a plurality of qubits, every qubit has a nonzero rate of dephasing, and some pairs of qubits have a nonzero rate of correlated dephasing. At step 602, vectors a and b whose elements are randomly chosen from {0,1} are generated. At step 604, the operation U.sub.ab prepares the state

[00002]$\frac{1}{\sqrt{2}}(.Math.a.Math.+.Math.b.Math.).$

At step bub the quantum system 100 evolves under dephasing noise for time t, represented by Et. At step 608, U.sub.ab.sup.† is applied and a computational basis measurement is performed. The probability of obtaining the outcome 0, P.sub.ab, decays exponentially (towards ½) as t increases. At step 610, matrix C 400 (FIG. 6) may be recovered by measuring the decay rate Γ.sub.ab for various a's and b's.

[0065] Here it is assumed that the entries in the matrix C 400 are real (i.e., with imaginary part equal to zero). This holds true in a number of cases, for instance, when the qubits are coupled to a bath at a high temperature. When C is complex, it can be handled using a generalization of the disclosed method. Note that physical constraints imply that C is positive semidefinite; hence C=C.sup.T in the real case, and C=C.sup.† in the complex case.

[0066] In aspects, the entangled states may undergo dephasing. The following procedure enables the measurement of certain linear functions of the correlation matrix C 400. This procedure is very general, and includes single- and two-qubit Ramsey spectroscopy as special cases. Consider an n-qubit state of the form

[00003]$\begin{array}{cc}.Math.{\psi }_{\mathrm{ab}}.Math.=\frac{1}{\sqrt{2}}(.Math.a.Math.+.Math.b.Math.)\in {\left({ℂ}^2\right)}^{.Math.n},& \left(\mathrm{Eq}.2.1\right)\end{array}$

where a, b∈{0,1}.sup.n, a≠b, |a=|a.sub.1, a.sub.2, . . . , a.sub.n and |b=|b.sub.1, b.sub.2, . . . , b.sub.n. By choosing a and b appropriately, one can make |ψ.sub.ab be a single-qubit |+ state, a two-qubit Bell state, or a many-qubit Greenberger-Horne-Zeilinger state (GHZ state) (while the other qubits are in a tensor product of |0 and |1 states).

[0067] For example, the state |ψ.sub.ab may be prepared, the state may be allowed to evolve for time t under the Lindbladian (Eq. 1.1). ρ(t) is the resulting density matrix. As can be seen in Eq. (1.2), the coherences (that is, the off-diagonal elements |a b|) of ρ(t) decay as exp(−Γ.sub.abt), where the decay rate Γ.sub.ab ∈ is defined so that (|ab|)=−Γ.sub.ab|ab|. This decay rate can be estimated by allowing the quantum system 100 to evolve for a suitable amount of time t, and then measuring in the basis

[00004]$\frac{1}{\sqrt{2}}(.Math.a.Math.+.Math.b.Math.).$

[0068] The decay rate Γ.sub.ab shows a certain linear function of the correlation matrix C, which can be written explicitly as follows. Let α.sub.i=(−1).sup.a.sup.i and β.sub.i=(−1).sup.b.sup.i denote the expectation values of Z.sub.i corresponding to the states |a and |b, respectively. In addition, define the vectors α=(α.sub.1, . . . , α.sub.n) and β=(β.sub.1, . . . , β.sub.n).

[00005]$\begin{array}{cc}{\Gamma }_{\mathrm{ab}}=-{.Math.}_{\mathrm{ij}}{c}_{\mathrm{ij}}\left[{\alpha }_i{\beta }_j-\frac{1}{2}{\alpha }_i{\alpha }_j-\frac{1}{2}{\beta }_i{\beta }_j\right]& \left(\mathrm{Eq}.2.2\right)\end{array}$$\begin{array}{cc}=2r^T\mathrm{Cr},& \left(\mathrm{Eq}.2.3\right)\end{array}$

[0069] Recall the definition of r, r=b−a, (Eq. 2.4) note that

[00006]$r=\frac{\alpha -\beta }{2},$

and the fact that C=C.sup.T is used to symmetrize the equation.

[0070] Single-qubit Ramsey spectroscopy (FIG. 4) is a special case of this procedure, where a=(0, 0, . . . , 0) and b=(0, . . . , 0,1,0, . . . , 0) (where the 1 appears in the j'th position). Then |ψ.sub.ab is a |+ state on the j'th qubit, and Γ.sub.ab=2c.sub.jj tells the rate of dephasing on the j'th qubit.

[0071] Two-qubit generalized spectroscopy is another special case, where a=(0, 0, . . . , 0) and b=(0, . . . , 0,1,0, . . . , 0,1,0, . . . , 0) (where the l's appear in the j'th and k'th positions). Then |ψ.sub.ab is a maximally-entangled state on qubits j and k, and Γ.sub.ab=2(c.sub.jj+c.sub.jk+c.sub.kj+c.sub.kk) provides information about the rate of correlated dephasing on qubits j and k.

[0072] An advantage of the disclosed methods is that the methods enable the characterization of errors. The characterized errors may be used to improve quantum systems (and devices) to reduce errors. For example, errors may be mitigated by applying a gate opposite to the characterized error. In another example, there are many quantum systems 100 that include a master clock, that is configured to track the phase of all of the qubits, so if dephasing occurs, the system loses the ability to track the phase of the qubits. The disclosed methods enable mitigation of the dephasing error based on the characterized error.

[0073] FIG. 5 is a graph that illustrates single qubit Ramsey spectroscopy. Ramsey spectroscopy is a form of particle interferometry that uses the phenomenon of magnetic resonance to measure transition frequencies of particles. The plot 502 shows the overlap P.sub.+, which decays exponentially (towards ½) with decay rate Γ, as a function of time t. The inset block diagram 504 illustrates the Ramsey protocol, where a superposition of qubit states is prepared with the first Hadamard gate H, the quantum system 100 undergoes dephasing (represented by the noise channel ε.sub.t) for time t, and the second H followed by a measurement in the computational basis measures the overlap P.sub.+.

[0074] In aspects, the decay rate Γ.sub.ab may be estimated by:

[0075] 1. Choosing some evolution time t≥0 such that

[00007]$\frac{1}{2}\le {\Gamma }_{\mathrm{ab}}t\le 2.$

This can be done in various ways, for instance, by starting with an initial guess t=τ.sub.0 and performing binary search, using ˜|log(Γ.sub.abτ.sub.0)| trials of the experiment.

[0076] 2. Repeating the following experiment N.sub.trials times: (N.sub.trials may be set using equation (2.11) below)

[0077] (a) Preparing the state

[00008]$.Math.{\psi }_{\mathrm{ab}}.Math.=\frac{1}{\sqrt{2}}(.Math.a.Math.+.Math.b.Math.),$

allow the state to evolve for time t, then measure in the basis

[00009]$\frac{1}{\sqrt{2}}(.Math.a.Math.+.Math.b.Math.).$

[0078] Letting N.sub.+ and N.sub.− be the number of

[00010]$\frac{1}{\sqrt{2}}(.Math.a.Math.+.Math.b.Math.)\mathrm{and}\frac{1}{\sqrt{2}}(.Math.a.Math.+.Math.b.Math.)$

outcomes, respectively. Note that the probabilities of these outcomes are given by

[00011]$P_{+}=\frac{1}{2}\left(1+e^{-{\Gamma }_{\mathrm{ab}}t}\right)\mathrm{and}{P}_{-}=\frac{1}{2}\left(1-e^{-{\Gamma }_{\mathrm{ab}}t}\right).$

[0079] 3. Defining:

[00012]$\begin{array}{cc}\Delta =\frac{N_{+}-N_{-}}{N_{t\mathrm{rials}}}.& (\mathrm{Eq}.2.5\end{array}$

[0080] Note that Δ is an unbiased estimator for P.sub.+−P.sub.−, that is, (Δ)=P.sub.+−P.sub.−=e.sup.−Γ.sup.ab.sup.t. This motivates the definition of an estimator for Γ.sub.ab:

[00013]$\begin{array}{cc}{\hat{\Gamma }}_{ab}=-\frac{1}{t}\ln \left(\Delta \right).& \left(\mathrm{Eq}.2.6\right)\end{array}$

[0081] Next, some bounds on the accuracy of the estimator {circumflex over (Γ)}.sub.ab may be set. To do this, the notion of a sub-Gaussian random variable is introduced (e.g., a random variable whose moments and tail probabilities behave like those of a Gaussian distribution). A real-valued random variable X is sub-Gaussian if there exists a real number K.sub.2 such that, for all p≥1, ((|X|.sup.p)).sup.1/p≤K.sub.2√{square root over (p)}. (Eq. 2.7)

[0082] The sub-Gaussian norm of X, denoted ∥X∥.sub.ψ.sub.2, is defined to be the smallest choice of K.sub.2 in (2.7), i.e.,

[00014]$\begin{array}{cc}{.Math.X.Math.}_{{\psi }_2}=\underset{p\ge 1}{\sup }{p^{-1/2}\left(𝔼\left({.Math.X.Math.}^p\right)\right)}^{1/p}.& \left(\mathrm{Eq}.2.8\right)\end{array}$

[0083] In addition, it is known that X is sub-Gaussian if and only if there exists a real number K.sub.1 such that, for all t≥0, Pr[|X|>t]≤exp(1−t.sup.2/K.sub.1.sup.2). (Eq. 2.9)

[0084] The smallest choice of K.sub.1 in (Eq. 2.9) is equivalent to the sub-Gaussian norm∥X∥.sub.ψ.sub.2, in the following sense: there is a universal constant c such that, for all sub-Gaussian random variables X, the smallest choice of K.sub.1 in (Eq. 2.9) differs from ∥X∥.sub.ψ.sub.2 by at most a factor of c.

[0085] {circumflex over (Γ)}.sub.ab−Γ.sub.ab is a sub-Gaussian random variable, whose sub-Gaussian norm is bounded by

[00015]$\begin{array}{cc}{.Math.{\hat{\Gamma }}_{ab}-{\Gamma }_{ab}.Math.}_{{\psi }_2}\le C_0\frac{{\Gamma }_{ab}}{\sqrt{N_{\mathrm{trials}}}},& \left(\mathrm{Eq}.2.1\right)\end{array}$

[0086] where C.sub.0 is a universal constant (e.g., 1/√{square root over (N.sub.trials)}) In particular, the accuracy of {circumflex over (Γ)}.sub.ab can be controlled by setting N trials appropriately: for any δ>0 and ∈>0, where

[00016]$\begin{array}{cc}{N}_{\mathrm{trials}}\ge \frac{2}{{\delta }^2}\ln \left(\frac{2}{ϵ}\right),& \left(\mathrm{Eq}.2.11\right)\end{array}$

[0087] then {circumflex over (Γ)}.sub.ab satisfies the following error bound: with probability at least 1−∈, |{circumflex over (Γ)}.sub.ab−Γ.sub.ab|≤2 δe.sup.2Γ.sub.ab. (Eq. 2.12)

[0088] This shows that the error in {circumflex over (Γ)}.sub.ab is at most a small fraction of the true value of Γ.sub.ab, independent of the magnitude of Γ.sub.ab.

[0089] In aspects, this may derive an error bound that involves {circumflex over (Γ)}.sub.ab rather than Γ.sub.ab, and may be computed from the observed value of {circumflex over (Γ)}.sub.ab. To show such an error bound, the triangle inequality may be used to write |{circumflex over (Γ)}.sub.ab−Γ.sub.ab|≤2δe.sup.2(|{circumflex over (Γ)}.sub.ab|+{circumflex over (Γ)}.sub.ab−Γ.sub.ab|), and divide by (1−2δe.sup.2) to obtain:

[00017]$\begin{array}{cc}.Math.{\hat{\Gamma }}_{ab}-{\Gamma }_{ab}.Math.\le \frac{2\delta e^2}{1-2\delta e^2}.Math.{\hat{\Gamma }}_{ab}.Math..& \left(\mathrm{Eq}.2.13\right)\end{array}$

[0090] It remains to prove (Eq. 2.10) and (Eq. 2.12). First use Hoeffding's inequality to show that Δ is close to its expectation value: Pr[|Δ−(P.sub.+−P.sub.−)|≥δ]≤2 exp(−N.sub.trialsδ.sup.2/2). (Eq. 2.14)

[0091] When |Δ−(P.sub.+−P.sub.−)|≤δ, the error in {circumflex over (Γ)}.sub.ab may be bounded as follows: (using the fact that P.sub.+−P.sub.−=e.sup.−Γ.sup.ab.sup.t≥e.sup.−2 and

[00018]$t\ge \frac{1}{2{\Gamma }_{ab}})$

[00019]$\begin{array}{cc}\begin{array}{c}-{\hat{\Gamma }}_{ab}\le \frac{1}{t}\ln \left(P_{+}-P_{-}+\delta \right)\\ \le \frac{1}{t}\left[\ln \left(P_{+}-P_{-}\right)+\delta e^2\right]\\ \le -{\Gamma }_{ab}+2{\Gamma }_{ab}\delta e^2,\end{array}& \left(\mathrm{Eq}.2.15\right)\end{array}$$\begin{array}{cc}\begin{array}{c}-{\hat{\Gamma }}_{ab}\ge \frac{1}{t}\ln \left(P_{+}-P_{-}-\delta \right)\\ \ge \frac{1}{t}\left[\ln \left(P_{+}-P_{-}\right)-\delta e^2\right]\\ \ge -{\Gamma }_{ab}-2{\Gamma }_{ab}\delta e^2.\end{array}& \left(\mathrm{Eq}.2.16\right)\end{array}$

[0092] Thus:

Pr[|{circumflex over (Γ)}.sub.ab−Γ.sub.ab|≥2δe.sup.2Γ.sub.ab]≤2 exp(−N.sub.trialsδ.sup.2/2).  (Eq. 2.17)

[0093] (Eq. 2.17) implies (Eq. 2.10) and (Eq. 2.12).

[0094] Referring to FIG. 6, matrix C 400 corresponding to the correlation graph in FIG. 3 is shown. The diagonal elements correspond to single qubit dephasing 402. The off-diagonal elements indicate correlated dephasing noise 404. In aspects, a direct estimation of the correlation matrix C 400 may be performed. The correlation matrix C 400 may be estimated directly, by performing single-qubit spectroscopy to measure the diagonal elements c.sub.jj, and by performing two-qubit spectroscopy to measure the off-diagonal elements c.sub.jk. This method may be used as a baseline to measure the performance of the disclosed compressed sensing method.

[0095] For simplicity, the case where C is real is considered. Since C is positive definite, this implies that c.sub.jj≥0 and c.sub.jk=c.sub.kj.

[0096] First, the diagonal elements c.sub.jj, for j=1, n, are estimated as follows:

[0097] 1. Let a=(0, 0, . . . , 0) and b=(0, . . . , 0,1,0, . . . , 0) (where the 1 appears in the j'th position).

[0098] 2. Construct an estimate {circumflex over (Γ)}.sub.ab of the decay rate Γ.sub.ab=2c.sub.jj. Define g.sub.j={circumflex over (Γ)}.sub.ab/2.

[0099] To write this in a compact form, define diag(C)=(c.sub.11, . . . , c.sub.nn), (Eq. 2.18)

[0100] that is, the diagonal of C. Next, define g=(g.sub.1, . . . , g.sub.n), (Eq. 2.19)

[0101] This is viewed as an estimator for diag(C).

[0102] Next, the off-diagonal elements c.sub.jk, for 1≤j<k≤n, are estimated as follows:

[0103] 1. Let a=(0, 0, . . . , 0) and b=(0, . . . , 0, 1, 0, . . . , 0, 1, 0, . . . , 0) (where the 1's appear in the j'th and k'th positions).

[0104] 2. Construct an estimate {circumflex over (Γ)}.sub.ab of the decay rate Γ.sub.ab=2(c.sub.jj+2c.sub.jk+c.sub.kk) and define

[00020]$h_{jk}=\frac{1}{4}{\hat{\Gamma }}_{ab}-\frac{1}{2}{g}_j-\frac{1}{2}{g}_k.$

[0105] To write this in a compact form, define C′ to be the matrix C 400 with the diagonal entries replaced by zeros,

[00021]$\begin{array}{cc}{C}_{jk}^{\prime }=\{\begin{array}{cc}{c}_{jk}& \mathrm{ifj}\ne k,\\ 0& \mathrm{ifj}=k.\end{array}& \left(\mathrm{Eq}.2.2\right)\end{array}$

[0106] The “off-diagonal portion” of matrix C is 2 s sparse. An estimator for Ĉ′ is constructed, as follows:

[00022]$\begin{array}{cc}{\hat{C}}_{jk}^{\prime }=\{\begin{array}{cc}{h}_{\mathrm{jk}}& \mathrm{ifj}<k,\\ h_{\mathrm{kj}}& \mathrm{ifj}>k,\\ 0& \mathrm{ifj}=k.\end{array}& \left(\mathrm{Eq}.2.21\right)\end{array}$

[0107] The accuracy of these estimators may be analyzed as follows. Initially, two parameters δ.sub.1 and δ.sub.2 are chosen. Suppose that, during the measurement of the diagonal elements c.sub.jj, the decay rates Γ.sub.ab are estimated with accuracy ∥{circumflex over (Γ)}.sub.ab−Γ.sub.ab∥.sub.ψ.sub.2≤δ.sub.1Γ.sub.ab, (Eq. 2.22) and during the measurement of the off-diagonal elements c.sub.jk, the decay rates Γ.sub.ab are estimated with accuracy ∥−Γ.sub.ab∥.sub.ψ.sub.2≤δ.sub.2Γ.sub.ab. (Γ.sub.ab (Eq. 2.23).

[0108] Bounds of the form (Eq. 2.22) and (Eq. 2.23) can be obtained from equation (Eq. 2.10), by setting N.sub.trials˜1/δ.sub.1.sup.2 and N.sub.trials˜1/δ.sub.2.sup.2, respectively. Counting those trials of the experiment that are used to choose the evolution time t may be neglected, because the number of those trials grows only logarithmically with Γ.sub.ab. δ.sub.1 and δ.sub.2 may be different, because in many experimental scenarios, measurements of c.sub.jj take less time than measurements of c.sub.jk. Using (Eq. 2.22) and (Eq. 2.23), it can be shown that bounds on the accuracy of g.sub.j and h.sub.jk are:

[00023]$\begin{array}{cc}{.Math.g_j-c_{\mathrm{jj}}.Math.}_{{\psi }_2}\le {\delta }_1{c}_{\mathrm{jj}},& \left(\mathrm{Eq}.2.24\right)\end{array}$$\begin{array}{cc}{.Math.h_{jk}-c_{jk}.Math.}_{{\psi }_2}\le {\delta }_2{c}_{jk}+\frac{1}{2}\left({\delta }_1+{\delta }_2\right)\left(c_{\mathrm{jj}}+c_{kk}\right).& \left(\mathrm{Eq}.2.25\right)\end{array}$

[0109] These bounds on the sub-Gaussian norm imply bounds on the moments, such as |X|≤∥X∥.sub.ψ.sub.2 and (X.sup.2)≤2∥X∥.sub.ψ.sub.2.sup.2.

[0110] The results to bound the accuracy of the estimators g and Ĉ′ may be followed. For g, using the .sub.1 and .sub.2 norms, the bounds are:

(∥g−diag(C)∥.sub.1)≤δ.sub.1∥diag(C)∥.sub.1,  (Eq. 2.26)

(∥g−diag(C)∥.sub.2.sup.2)≤2δ.sub.1.sup.2∥diag(C)∥.sub.2.sup.2.  (Eq. 2.27)

[0111] For Ĉ′, using the .sub.1 vector norm and the Frobenius matrix norm, the bounds are:

[00024]$\begin{array}{cc}𝔼\left({.Math.{\hat{C}}^{\prime }-C^{\prime }.Math.}_{{\ell }_1}\right)\le 2{.Math.}_{j<k}\left({\delta }_2{c}_{\mathrm{jk}}+\frac{1}{2}\left({\delta }_1+{\delta }_2\right)\left(c_{\mathrm{jj}}+c_{\mathrm{kk}}\right)\right)={\delta }_2{.Math.}_{j\ne k}{c}_{\mathrm{jk}}+\frac{1}{2}\left({\delta }_1+{\delta }_2\right){.Math.}_{j\ne k}\left(c_{\mathrm{jj}}+c_{\mathrm{kk}}\right)\le \left(n-1\right)\left({\delta }_1+{\delta }_2\right).Math.\mathrm{diag}\left(C\right){.Math.}_{{\ell }_1}+{\delta }_2{.Math.C^{\prime }.Math.}_{{\ell }_1},& \left(\mathrm{Eq}.2.28\right)\end{array}$$\begin{array}{cc}𝔼\left({.Math.{\hat{C}}^{\prime }-C^{\prime }.Math.}_F^2\right)\le 2{.Math.}_{j<k}{2\left[{\delta }_2{c}_{\mathrm{jk}}+\frac{1}{2}\left({\delta }_1+{\delta }_2\right)\left(c_{\mathrm{jj}}+c_{\mathrm{kk}}\right)\right]}^2\le 2{.Math.}_{j<k}6\left[{\delta }_2^2{c}_{\mathrm{jk}}^2+\frac{1}{4}{\left({\delta }_1+{\delta }_2\right)}^2{c}_{\mathrm{jj}}^2+\frac{1}{4}{\left({\delta }_1+{\delta }_2\right)}^2{c}_{\mathrm{kk}}^2\right]\le 6{\delta }_2^2{.Math.C^{\prime }.Math.}_F^2+\frac{3}{2}{\left({\delta }_1+{\delta }_2\right)}^2{.Math.}_{j\ne k}\left(c_{\mathrm{jj}}^2+c_{\mathrm{kk}}^2\right)\le 3\left(n-1\right){\left({\delta }_1+{\delta }_2\right)}^2.Math.\mathrm{diag}\left(C\right){.Math.}_2^2+6{\delta }_2^2{.Math.C^{\prime }.Math.}_F^2.& \left(\mathrm{Eq}.2.29\right)\end{array}$

[0112] Note that, in the second step of (Eq. 2.29), the fact that for any real numbers a.sub.1, a.sub.2, a.sub.3, (a.sub.1+a.sub.2+a.sub.3).sup.2≤3(a.sub.1.sup.2+a.sub.2.sup.2+a.sub.3.sup.2) was used.

[0113] These are bounds on the expected error of g and {tilde over (C)}′. One can then use Markov's inequality to prove bounds that hold with high probability. For instance, using (Eq. 2.27), with probability at least 1−η, results in

[00025]$\begin{array}{cc}{.Math.g-\mathrm{diag}\left(C\right).Math.}_2\le \frac{1}{\sqrt{\eta }}\sqrt{2}{\delta }_1{.Math.\mathrm{diag}\left(C\right).Math.}_2,& \left(\mathrm{Eq}.2.3\right)\end{array}$

[0114] Using (Eq. 2.29), with probability at least 1−η, results in:

[00026]$\begin{array}{cc}{.Math.{\hat{C}}^{\prime }-C^{\prime }.Math.}_F\le {\frac{1}{\sqrt{\eta }}\left[3\left(n-1\right){\left({\delta }_1+{\delta }_2\right)}^2{.Math.\mathrm{diag}\left(C\right).Math.}_2^2+6{\delta }_2^2.Math.C^{\prime }.Math.\right]}^{1/2}\le \frac{1}{\sqrt{\eta }}\left[\sqrt{3}\sqrt{n-1}\left({\delta }_1+{\delta }_2\right){.Math.\mathrm{diag}\left(C\right).Math.}_2+\sqrt{6}{\delta }_2{.Math.C^{\prime }.Math.}_F\right].& \left(\mathrm{Eq}.2.31\right)\end{array}$

[0115] In fact, one can prove tighter bounds on the failure probability, by replacing Markov's inequality with a sharp concentration bound for sub-Gaussian random variables.

[0116] The bounds (Eq. 2.30) and (Eq. 2.31) give a rough sense of how well this estimator performs. In particular, these bounds show that the error in estimating C′ depends on the magnitude of diag(C), as well as the magnitude of C′. This is due to the fact that the procedure for estimating the off-diagonal matrix element c.sub.jk also involves the diagonal matrix elements c.sub.jj and c.sub.jk. These bounds may be used as a baseline to understand the performance of the disclosed compressed sensing estimators.

[0117] FIGS. 7A and 7B are tables illustrating sample complexity of different methods for reconstructing a correlation matrix C 400 (FIG. 6), of size n×n, with 2 s nonzero elements off the diagonal. In aspects, the naive method is to measure each element of C separately. “RIP-based” and “RIPless” refer to different analytical bounds on the accuracy of the reconstruction of matrix C 400 (FIG. 6). The asterisk (*) indicates that the results using the RIPless bound hold under a technical assumption that the diagonal of C does not have any unusually large elements. The different methods are parameterized in such a way that they reconstruct the diagonal of matrix C 400 (FIG. 6) up to an additive error of size δ∥diag(C)∥.sub.2, and they reconstruct the off-diagonal part of C up to an additive error of size δ∥C∥.sub.F. Each method makes use of single-qubit spectroscopy (with n experimental configurations or “measurement settings”), as well as multi-qubit spectroscopy (with m measurement settings, where m varies between m∝n.sup.2 and m∝s log n). The total sample complexity, shown in the tables of FIGS. 7A and 7B, include both single-qubit and multi-qubit spectroscopy. For the compressed sensing method, the number of measurement settings can be as low as m∝s log n, using the RIPless bound, but the best sample complexity is achieved when m∝s log.sup.4n, using the RIP-based bound. This sample complexity, O(max(n, s).sup.2 log.sup.4(n)), compares favorably with the sample complexity of the naive method, which is O(n.sup.3). Thus, Compressed sensing has an advantage over the naive method whenever s≤O(n.sup.3/2/log.sup.2n).

[0118] The performance of the disclosed compressed sensing method may be compared to the naive method where one measures each element of C independently. A rigorous comparison may be made based on the error bounds for each of these methods. These can be summarized as follows:

Naive:∥Ĉ−C∥.sub.F≤O(√{square root over (n)}(δ.sub.1+δ.sub.2)∥diag(C)∥.sub.2+δ.sub.2∥C′∥.sub.F),  (Eq. 1.4)

RIP-based:∥W.sup.(opt)−C∥.sub.F≤O(√{square root over (n)}δ.sub.1∥diag(C)∥.sub.2+δ.sub.2(√{square root over (n)}∥diag(C)∥.sub.2+√{square root over (2s)}∥C′∥.sub.F)),  (Eq. 1.5)

RIPless:∥W.sup.(opt)−C∥.sub.F≤O(s log.sup.5/2(n)[δ.sub.1∥diag(C)∥.sub.∞√{square root over (sn log)}(n)+δ.sub.2√{square root over (n)}∥diag(C)∥.sub.2+δ.sub.2√{square root over (2s)}∥C′∥.sub.F]).  (Eq. 1.6)

[0119] Here Ĉ is the estimate of C using the naive method, and W.sup.(opt) is the estimate of C using compressed sensing. diag(C) denotes the diagonal of C, and C′ denotes the off-diagonal part of C. m∝s log.sup.4n is set for the RIP-based bound, and m∝s log n for the RIPless bound. (2.30) may be used to bound the error in estimating diag(C), and (Eq. 2.31), (Eq. 5.48) and (Eq. 5.39) to bound the error in estimating C′. Finally, δ.sub.1 and δ.sub.2 are parameters that control the accuracy of the single-qubit and multi-qubit spectroscopy procedures. Using the disclosed method rigorous bounds on the sample complexity may be proven that are required for each method to reconstruct C with comparable accuracy.

[0120] For the compressed sensing method, the number of measurement settings can be as low as m∝s log n, using the RIPless bound, but the best sample complexity is achieved when m∝s log.sup.4n, using the RIP-based bound. This sample complexity, O(max(n, s).sup.2 log.sup.4(n)), compares favorably with the sample complexity of the naive method, which is O(n.sup.3). Thus, compressed sensing has an advantage over the naive method whenever s≤O(n.sup.3/2/log.sup.2n).

[0121] The estimation of the decay rate Γ.sub.ab=2r.sup.TCr may be performed with high precision (i.e., with small multiplicative error), by allowing the quantum system 100 to evolve for a time t˜1/Γ.sub.ab. Here the time t is chosen by an adaptive procedure that starts with an initial guess τ.sub.0 and converges (provably) after˜|log(Γ.sub.ab τ.sub.0)| steps (see FIG. 3). To achieve high precision, the time evolution operator is treated exactly.

[0122] The above-noted feature is helpful in experimental setups where measurements are time-consuming and the values of Γ.sub.ab span several orders of magnitude, thus making it difficult to determine the appropriate evolution time and accurately measure the decay rates using conventional methods. This feature can also be used as a standalone technique to estimate the relaxation (T.sub.1) and decoherence (T.sub.2) times of any quantum system 100.

[0123] Since the disclosed method relies on estimating exponential decay rates, it can be made at least partially robust to state-preparation-and-measurement (SPAM) errors (see FIGS. 9A and 9B). This follows some of the same approaches used in randomized benchmarking and gate set tomography. The estimated decay rate may be used to estimate a relaxation time of the quantum system 100. In aspects, a decoherence time of the quantum system 100 may be estimated based on the estimated decay rate.

[0124] The disclosed method enables a unitary evolution according to some Hamiltonian H.sub.S, and enables the matrix C 400 (FIG. 6) to have complex entries c.sub.jk. A similar approach can be used to estimate the Hamiltonian H.sub.S, as well as the imaginary part of c.sub.jk. The disclosed technology measures correlated dephasing and not Pauli errors. In aspects, the measurement may include Ramsey spectroscopy. In aspects the reconstruction algorithm may include, for example, convex optimization over a continuous domain.

[0125] The disclosed method assumes that the noise is sparse, meaning that the noise is supported on a subset S of the domain, where S is small, but unknown.

[0126] It is contemplated that the disclosed method for SPAM-robust estimation of sparse noise models, may be used to characterize other kinds of correlated noise processes in many-body quantum systems.

[0127] The disclosed technology includes an efficient method for learning the off-diagonal part of the correlation matrix C, under the assumption that it is sparse, i.e., the part of C that lies above the diagonal has at most s nonzero elements, where s<<n(n−1)/2. (Since C is Hermitian, it is sufficient to learn the part that lies above the diagonal; this then determines the part that lies below the diagonal.)

[0128] For simplicity, first consider the special case where the entries in the matrix C are real (i.e., with zero imaginary part), which occurs in a number of physical situations (for example, when the system is coupled to a bath at a high temperature.

[0129] The disclosed method consists of two steps: first, single-qubit Ramsey spectroscopy is performed in order to learn the diagonal elements of C; second, techniques from compressed sensing (e.g., random measurements, and .sub.1-minimization) are applied in order to recover the off-diagonal elements of C. In aspects, a long-range correlated dephasing error may be detected based on the recovered matrix.

[0130] In aspects, random measurements of the correlation matrix may be performed. Initially, each of the diagonal elements c.sub.jj, for j=1, . . . , n, are estimated using single-qubit Ramsey spectroscopy, as illustrated in FIG. 5. Let g=(g.sub.1, . . . , g.sub.n)∈.sup.n be the output of this procedure. g may be viewed as an estimate of a “sensing operator” that returns the diagonal elements of the matrix C 400,

g≈diag(C)=(c.sub.11,c.sub.22, . . . ,c.sub.nn).  (Eq. 3.1)

[0131] (Note that c.sub.jj≥0, since C is positive semidefinite.)

[0132] In order to estimate the off-diagonal part of C, a compressed sensing technique is used, which involves a certain type of generalized Ramsey measurement with random GHZ-type states, see FIG. 4. First, a parameter m is chosen, which can be roughly m˜s log n or m˜s log.sup.4n, which controls the number of different measurements. The particular choice of m may be motivated by a theoretical recovery guarantee. Now, for j=1, . . . , m, perform the following procedure:

[0133] 1. Choose vectors a, b∈{0,1}.sup.n uniformly at random. As in equation (2.4), define

r=b−a.  (Eq. 3.2)

[0134] 2. Prepare the state

[00027]$.Math.{\psi }_{\mathrm{ab}}.Math.=\frac{1}{\sqrt{2}}(.Math.a.Math.+.Math.b.Math.).$

This is a GHZ state on a subset of the qubits, with some bit flips. It can be created by preparing those qubits i where a.sub.i=b.sub.i in the state |a.sub.i, preparing those qubits i where a.sub.i≠b.sub.i in a GHZ state, and applying a Pauli X operator on those qubits i where a.sub.i>b.sub.i. (This requires a quantum circuit of depth ┌log.sub.2(n)┐+2.)

[0135] 3. Construct an estimate {circumflex over (Γ)}.sub.ab of the decay rate Γ.sub.ab=2r.sup.TCr. Define h.sub.j={circumflex over (Γ)}.sub.ab.

[0136] Let h=(h.sub.1, . . . , h.sub.m)∈.sup.m be the output of the above procedure. Again, h may be viewed as an estimate of a “sensing operator” Φ:.sup.n×n.fwdarw..sup.m that is applied to the matrix C 400 (FIG. 6),

h≈Φ(C)=[Φ.sub.j(C)].sub.j=1, . . . ,m  (Eq. 3.3)

Φ.sub.j(C)=2(r.sup.(j)).sup.TCr.sup.(j)  (Eq. 3.4)

[0137] where r.sup.(1), r.sup.(2), . . . , r.sup.(m)∈{1,0,−1}.sup.n are independent random vectors chosen from the same distribution as r (described above). Note that Φ.sub.j(C)≥0, since C is positive semidefinite. The factor of 2 is chosen to ensure that Φ has a certain isotropy property, as described herein.

[0138] In aspects, the correlation matrix C∈.sup.n×n may be reconstructed. C is positive semidefinite, due to physical constraints, and its off-diagonal part is sparse, with at most s<<n(n−1)/2 nonzero elements above the diagonal. In general, this sparsity constraint leads to an optimization problem that is computationally intractable. However, in this particular case, this problem can be solved using a strategy from compressed sensing: given g≈diag(C)∈.sup.n and h≈Φ(C) ∈.sup.m, matrix C 400 is recovered by solving a convex optimization problem, where the .sub.1 (vector) norm of the off-diagonal part of the matrix is minimized. This strategy succeeds when m≥c.sub.0s log n, and is highly robust to noise when m≥c.sub.0s log.sup.4n, where c.sub.0 is some universal constant.

[0139] Initially, the case where the measurements are noiseless, i.e., g=diag(C) and h=Φ(C) is considered. The following convex optimization problem is solved:

[00028]$\begin{array}{cc}\mathrm{Find}W\in {ℝ}^{n\times n}\mathrm{that}\mathrm{minimizes}{.Math.}_{i\ne j}.Math.W_{\mathrm{ij}}.Math.,& \left(\mathrm{Eq}.3.5\right)\end{array}$$\mathrm{such}\mathrm{that}:$$\begin{array}{cc}\mathrm{diag}\left(W\right)=g,& \left(\mathrm{Eq}.3.6\right)\end{array}$$\begin{array}{cc}\Phi \left(W\right)=h,& \left(\mathrm{Eq}.3.7\right)\end{array}$$\begin{array}{cc}W⪰0.& \left(\mathrm{Eq}.3.8\right)\end{array}$

[0140] Here, W0 means that W is positive semidefinite, which implies that W=W.sup.T. As a sanity check, note that W=C is a feasible solution to this problem. (Recall that C is positive semidefinite.)

[0141] In phase retrieval, one wishes to estimate an unknown vector x from measurements of the form |r.sup.Tx|.sup.2. In aspects, a PhaseLift algorithm works by “lifting” the unknown vector x to a matrix X=xx.sup.T, so that the problem becomes one of learning a rank-1 matrix X from measurements of the form r.sup.TXr; then one solves a convex relaxation of this problem. It is contemplated that in cases where the unknown vector x is sparse (“compressive phase retrieval”), variants of the PhaseLift algorithm (as well as other approaches) can also be used.

[0142] In the disclosed method, the unknown matrix C is almost always full-rank (because every qubit has a nonzero rate of dephasing). Physical constraints imply that C is positive semidefinite, so it can be factored as C=BB.sup.T, which is superficially similar to X=xx.sup.T; however, an important difference is that B is a square matrix, whereas x is a vector.

[0143] In aspects, the matrix C 400 may be reconstructed from noisy measurements. In the case where the measurements of g and h are noisy, the above convex optimization problem is modified, by relaxing the constraints (Eq. 3.6) and (Eq. 3.7). This leads to some technical complications, due to the fact that two variables are being reconstructed that have different characteristics: the diagonal part of C (which is not sparse), and the off-diagonal part of C (which is sparse).

[0144] To deal with these issues, two different ways of performing this reconstruction, when the measurements are noisy are: (1) simultaneous reconstruction of both parts of C, and (2) sequential reconstruction of the diagonal part of C, followed by the off-diagonal part of C. The former approach is arguably more natural, but the latter approach enables more rigorous analysis of the accuracy of the reconstruction.

[0145] Assuming the following bounds on the .sub.2 norms of the noise terms (u and v), that is,

g=diag(C)+u,∥u∥.sub.2≤∈.sub.1,  (Eq. 3.9)

h=Φ(C)+v,∥v∥.sub.2≤∈.sub.2.  (Eq. 3.10)

[0146] Simultaneous reconstruction of the diagonal and off-diagonal parts of C: the constraints (Eq. 3.6) and (Eq. 3.7) are relaxed in the simplest possible way, by replacing them with:

∥diag(W)−g∥.sub.2≤∈.sub.1,  (Eq. 3.11)

∥Φ(W)−h∥.sub.2≤∈.sub.2.  (Eq. 3.12)

[0147] This leads to a convex optimization problem that attempts to reconstruct both the diagonal part of C, which is not necessarily sparse, and the off-diagonal part of C, which is assumed to be sparse. (Note that W=C is a feasible solution to this problem.)

[0148] The reconstruction algorithm involves two different estimators (an .sub.1-regularized estimator for the off-diagonal part of C, and a least-squares estimator for the diagonal part of C). These two estimators are coupled together (through the constraint on Φ(W), and the positivity constraint W0).

[0149] Therefore, this method can behave quite differently, depending on whether the dominant source of error is g or h. When g is the dominant source of error, this method will behave like a least-squares estimator, whose accuracy scales according to the dimension n; when h is the dominant source of error, this method will behave like an .sub.1-regularized estimator, whose accuracy scales according to the sparsity s (neglecting log factors). From a theoretical point of view, this makes it more difficult to prove recovery guarantees for this method.

[0150] Sequential reconstruction of the diagonal part of C, followed by the off-diagonal part of C: In practice, one is often interested in the regime where g is known with high precision, and h is the dominant source of error. This is because measurements of g are relatively easy to perform, because they only require single-qubit state preparations and measurements; whereas measurements of h are more costly, because they require the preparation and measurement of entangled states on many qubits. So measurements of g can often be performed more quickly, and measurements on different qubits can be performed simultaneously in parallel; hence one can repeat the measurements more times, to obtain more accurate estimates of g.

[0151] In this regime, it is natural to try to recover the diagonal part of C directly from g, and then use .sub.1-minimization to recover only the off-diagonal part of C. This leads to a convex optimization problem which is arguably less natural, but it makes it easier to prove rigorous guarantees on the accuracy of the reconstruction of C. Starting from the convex optimization problem ((Eq. 3.5)-(Eq. 3.8)) for the noiseless case, and the last two constraints are relaxed to get:

[00029]$\begin{array}{cc}\mathrm{Find}W\in {ℝ}^{n\times n}\mathrm{that}\mathrm{minimizes}{.Math.}_{i\ne j}.Math.W_{\mathrm{ij}}.Math.,& \left(\mathrm{Eq}.3.13\right)\end{array}$$\mathrm{such}\mathrm{that}:$$\begin{array}{cc}\mathrm{diag}\left(W\right)=g,& \left(\mathrm{Eq}.3.14\right)\end{array}$$\begin{array}{cc}{.Math.\Phi \left(W\right)-h.Math.}_2\le {ϵ}_2+{ϵ}_1\sqrt{\mathrm{mn},}& \left(\mathrm{Eq}.3.15\right)\end{array}$$\begin{array}{cc}d\left(W,K_{+}\right)\le {ϵ}_1,\mathrm{and}W=W^T.& \left(\mathrm{Eq}.3.16\right)\end{array}$

[0152] Here K.sub.+={W′∈.sup.n×n|W′0} denotes the (real) positive semidefinite cone, and

[00030]$\begin{array}{cc}d\left(W,K_{+}\right)=\underset{W^{\prime }\in K_{+}}{\min }{.Math.W-W^{\prime }.Math.}_F& \left(\mathrm{Eq}.3.17\right)\end{array}$

[0153] to be the minimum distance from W to a point W′ in K.sub.+, measured in Frobenius norm; note that this is a convex function. While this convex optimization problem looks complicated, it follows from a simple underlying idea: since the diagonal elements of W are fixed by the constraint (Eq. 3.14), this is simply an .sub.1-regularized estimator for the sparse, off-diagonal part of W.

[0154] The attentive reader will notice two potential concerns with this approach. First, in general, C will not be a feasible solution to this convex optimization problem, since the diagonal elements of C will not satisfy (Eq. 3.14). However, C lies close to a feasible solution. To see this, let {tilde over (C)} be the matrix whose off-diagonal elements agree with C, and whose diagonal elements agree with g. Then C is within distance ∈.sub.1 of {tilde over (C)} (in Frobenius norm), and {tilde over (C)} is a feasible solution. To see this, check that {tilde over (C)} satisfies the constraints (Eq. 3.14), (Eq. 3.15) and (Eq. 3.16), since:

[00031]$\begin{array}{cc}\begin{array}{c}{.Math.\Phi \left(\tilde{C}\right)-\Phi \left(C\right).Math.}_2={.Math.\Phi \left(\mathrm{diag}\left(g-\mathrm{diag}\left(C\right)\right)\right).Math.}_2\\ \le {\left(m{.Math.g-\mathrm{diag}\left(C\right).Math.}_1^2\right)}^{1/2}\\ \le {ϵ}_1\sqrt{\mathrm{mn}.}\end{array}& \left(3.18\right)\end{array}$

[0155] (Here, {tilde over (C)} is written in a compact form, {tilde over (C)}=C+diag (g−diag(C)), where the diag(.Math.) notation has the following meaning: for a matrix M, diag (M) is the vector containing the entries M.sub.jj that lie along the diagonal of M; and for a vector v=(v.sub.1, . . . , v.sub.n), diag(v) is the diagonal matrix with v.sub.1, . . . , v.sub.n along the diagonal.)

[0156] Second, the optimal solution W may violate the positivity constraint (Eq. 3.8), making it un-physical/non-physical. (Similar issues can arise when performing quantum state and process tomography.) However, W can be easily corrected to get a physically-admissible solution. This follows because equation (3.16) shows that W lies within distance ∈.sub.1 of a physically-admissible solution W′0, and this solution W′ can be obtained by truncating the negative eigenvalues of W.

[0157] Finally, there are different ways of relaxing the positivity constraint (Eq. 3.8), and (Eq. 3.16) is not the strongest possible choice. For instance, a stronger constraint could be used than (Eq. 3.16), such as: d′(W,K.sub.+)≤∈.sub.1, where d′(W,K.sub.+)=min.sub.W′∈K.sub.+∥diag (W−W′)∥.sub.2. However, the constraint (Eq. 3.16) may be simpler to implement using numerical linear algebra software.

[0158] Since these are convex optimization problems, they can be solved efficiently (both in theory and in practice), for instance by using interior point algorithms. Nonetheless, some care is needed to ensure that these algorithms can scale up to solve very large instances of these problems. In particular, enforcing the positivity constraint (Eq. 3.8), and its relaxed version (Eq. 3.16), can be computationally expensive.

[0159] In aspects, the positivity constraint may be omitted. The matrix C 400 can be reconstructed by solving the convex optimization problem (Eq. 3.13)-(Eq. 3.16). This analysis holds even without the positivity constraint (Eq. 3.16). In aspects, it is easy to check that the positivity constraint is not used, and indeed, most of the theory of compressed sensing applies to all sparse signals, not just positive ones, although positivity can be helpful in certain situations.

[0160] This observation has a practical consequence: by omitting the positivity constraint (Eq. 3.16), one can make the convex optimization problem simpler, and thus easier to solve in practice (e.g., by using second-order cone programming, rather than semidefinite programming). One can then take the resulting solution, and project it onto the positive semidefinite cone, as is sometimes done in quantum state tomography, without increasing the error (in Frobenius norm). This technique is useful for scaling up the disclosed method to extremely large numbers of qubits.

[0161] Next, setting the Error Parameters ∈.sub.1 and ∈.sub.2 is described. To set the parameters ∈.sub.1 and ∈.sub.2 in equations (3.9) and (3.10): First, consider ∈.sub.1, which bounds u, the error in g. Note that, when g.sub.j is estimated using the above noted procedure, large-deviation bounds on u.sub.j are obtained. In particular, for some δ.sub.1>0, ∥u.sub.j∥.sub.ψ.sub.2≤δ.sub.1c.sub.jj, (Eq. 3.19)

[0162] where ∥.Math.∥.sub.ψ.sub.2 is the sub-Gaussian norm. This bound can be obtained from equation (2.10), by setting N.sub.trials˜1/δ.sub.1.sup.2. The trials of the experiment that are used to choose the evolution time t may be neglected in the count, because the number of those trials grows only logarithmically with Γ.sub.ab.

[0163] This implies that ∥u∥.sub.2.sup.2 is a subexponential random variable, whose subexponential norm is at most 2δ.sub.1.sup.2∥diag(C)∥.sub.2.sup.2. This implies that ∥u∥.sub.2 is bounded with high probability: for any τ≥1, Pr[∥u∥.sub.2≥τδ.sub.1∥diag(C)∥.sub.2]≤e.Math.exp(−τ.sup.2/2c), (Eq. 3.20)

[0164] where c>0 is some universal constant. τ may be chosen to be a sufficiently large constant, so that the failure probability is small. (Note that in some cases, one can prove stronger bounds, by taking advantage of the fact that the coordinates of u are independent, and using a Bernstein-type inequality. This bound is stronger when the diagonal elements of C satisfy∥diag(C)∥.sub.∞<<∥diag(C)∥.sub.2, i.e., when u has many independent coordinates with similar magnitudes.)

[0165] The above bound does not immediately tell how to set ∈.sub.1, because the bound depends on diag(C), which is not known exactly. Instead, a bound that depends on g is derived, which is known explicitly, and can be used to set ∈.sub.1. To do this, assume that δ.sub.1 is sufficiently small so that τδ.sub.1<1/4. With high probability:

[00032]$\begin{array}{cc}\begin{array}{c}{.Math.u.Math.}_2\le {\mathrm{\tau \delta }}_1{.Math.\mathrm{diag}\left(C\right).Math.}_2\\ \le {\mathrm{\tau \delta }}_1\left({.Math.g.Math.}_2+{.Math.u.Math.}_2\right)\\ \le \frac{{\mathrm{\tau \delta }}_1}{1-{\mathrm{\tau \delta }}_1}{.Math.g.Math.}_2\text{=:}{ϵ}_1,\end{array}& \left(3.21\right)\end{array}$

[0166] Using the triangle inequality, and some algebra shows how to set ∈.sub.1 so that equation (3.9) holds.

[0167] ∈.sub.1 may also be bounded in terms of diag(C), as follows:

[00033]$\begin{array}{cc}\begin{array}{c}{ϵ}_1\le \frac{{\mathrm{\tau \delta }}_1}{1-{\mathrm{\tau \delta }}_1}\left(.Math.\mathrm{diag}\left(C\right){.Math.}_2+{.Math.u.Math.}_2\right)\\ \le \frac{{\mathrm{\tau \delta }}_1}{1-{\mathrm{\tau \delta }}_1}\left(1+{\mathrm{\tau \delta }}_1\right){.Math.\mathrm{diag}\left(C\right).Math.}_2\\ \le 2{\mathrm{\tau \delta }}_1{.Math.\mathrm{diag}\left(C\right).Math.}_2.\end{array}& \left(\mathrm{Eq}.3.22\right)\end{array}$ [0168] In aspects, this bound may be used to analyze the accuracy of the estimate of matrix C 400. ∈.sub.2 bounds v, which is the error in h. The same approach as above may be used. When h.sub.j is estimated a bound on the sub-Gaussian norm of v.sub.j is obtained: for some δ.sub.2>0,

∥v.sub.j∥.sub.ψ.sub.2≤δ.sub.2Φ.sub.j(C).  (Eq. 3.23)

[0169] (This bound can be obtained from equation (2.10), by setting N.sub.trials˜1/δ.sub.2.sup.2. δ.sub.2 to may be different from δ.sub.1, because the measurements used to estimate h.sub.j are more costly than the measurements used to estimate g.sub.j, hence one may prefer to use different values for N.sub.trials in each case.)

[0170] This implies that ∥v∥.sub.2.sup.2 is a subexponential random variable, hence ∥v∥.sub.2 is bounded with high probability: for any τ≥1,

Pr[∥v∥.sub.2≥τδ.sub.2∥Φ(C)∥.sub.2]≤e.Math.exp(−τ.sub.2/2c),  (Eq. 3.24)

[0171] where c>0 is some universal constant, choose τ to be a sufficiently large constant, so that the failure probability is small. (Note that, for typical choices of Φ(.Math.), it is expected that ∥Φ(C)∥.sub.∞<<∥Φ(C)∥.sub.2. This implies that v has many independent coordinates with similar magnitudes. When this occurs, one can prove a stronger bound using a Bernstein-type inequality.

[0172] The above bound does not immediately show how to set ∈.sub.2, because the bound depends on Φ(C), which is not known exactly. Instead, derive a bound that depends on h, which is known explicitly, and can be used to set ∈.sub.2. To do this, it is assumed that δ.sub.2 is sufficiently small so that τδ.sub.2<1/4. With high probability:

[00034]$\begin{array}{cc}\begin{array}{c}{.Math.v.Math.}_2\le {\mathrm{\tau \delta }}_2{.Math.\Phi \left(C\right).Math.}_2\\ \le {\mathrm{\tau \delta }}_2\left({.Math.h.Math.}_2+{.Math.v.Math.}_2\right)\\ \le \frac{{\mathrm{\tau \delta }}_2}{1-{\mathrm{\tau \delta }}_2}{.Math.h.Math.}_2\text{=:}{ϵ}_2.\end{array}& \left(3.25\right)\end{array}$

[0173] This indicates how to set ∈.sub.2 so that equation (3.10) holds. Finally, ∈.sub.2 can also be bounded in terms of ∥C.sub.1, as follows:

[00035]$\begin{array}{cc}\begin{array}{c}{ϵ}_2\le \frac{{\mathrm{\tau \delta }}_2}{1-{\mathrm{\tau \delta }}_2}\left({.Math.\Phi \left(C\right).Math.}_2+{.Math.v.Math.}_2\right)\\ \le \frac{{\mathrm{\tau \delta }}_2}{1-{\mathrm{\tau \delta }}_2}\left(1+{\mathrm{\tau \delta }}_2\right){.Math.\Phi \left(C\right).Math.}_2\\ \le 2{\mathrm{\tau \delta }}_2{.Math.\Phi \left(C\right).Math.}_2\\ \le 2{\mathrm{\tau \delta }}_2\sqrt{m}{.Math.\Phi \left(C\right).Math.}_{\infty }\\ \le 4{\mathrm{\tau \delta }}_2\sqrt{m}{.Math.C.Math.}_{{\ell }_1.}\end{array}& \left(3.26\right)\end{array}$

[0174] In aspects, this bound may be used when analyzing the accuracy of the estimate of C.

[0175] Referring to FIGS. 8A-8C, graphs that illustrating scaling of the reconstruction error ∥W.sup.(opt)−C∥.sub.∞ under various circumstances are shown. Here W.sup.(opt) denotes the estimate of C obtained via compressed sensing. The solid lines are the average errors over 100 random instances of the problem, and the shaded region is their 95% confidence interval obtained by bootstrapping. FIG. 8A illustrates the reconstruction error as a function of the number of measurement settings m (assuming noiseless measurements) for various values of sparsity s (and n=64 qubits). The errors go through a phase transition whose location m.sub.c scales linearly with s. FIG. 8B illustrates the reconstruction error as a function of the number of measurement settings m (assuming noiseless measurements) for various numbers of qubits n (and sparsity s=12). The phase transition point m.sub.c scales logarithmically with n. FIG. 8C illustrates the reconstruction error as a function of the number of measurement settings m, for different values of added noise strength σ.sub.∈, with fixed parameters (n, s)=(64,12). The inset shows that the recovery errors (when m>m.sub.c, i.e., after the transition point) scale linearly with σ.sub.∈, as expected.

[0176] Numerical simulations are used to test the performance of the disclosed method on realistic system sizes, with different levels of sparsity, and when the data contain statistical fluctuations due to finite sample sizes. As shown in FIGS. 8A-8C, the disclosed method performs well overall.

[0177] The protocol for randomly chosen C matrices is numerically simulated. In these examples it is assumed that the diagonal elements of C are known, that is, ∈.sub.1=0 in Eq. (3.11). The convex optimization problem given by equations (3.5)-(3.7) is solved using, for example, CVXPY, a convex optimization package for Python®.

[0178] The case of noiseless measurements, corresponding to ∈.sub.2=0 in Eq. (3.12) is investigated. In FIG. 8A the recovery error is shown as a function of the number of measurements, m, for a fixed number of qubits, n, and various choices of the off-diagonal sparsity, s. The sharp transition in the recovery error as a function of m is evident. Moreover, as shown in the inset of FIG. 8A, the transition point m.sub.c, which is defined as the point where ∥C−W.sup.(opt)∥.sub.∞ drops below 0.25, scales linearly with s, consistent with the analytical results. In FIG. 8B s is fixed, n is varied, and the recovery error as a function of m is studied. A phase transition may be observed as m increases. In this case, m.sub.c scales polynomially with log(n) as suggested in the inset of FIG. 8B.

[0179] The effect of noisy measurements on the recovery error is investigated. Random C matrices may be generated, with a fixed number of qubits n and sparsity s. Noise is simulated by adding a random vector e, whose entries are independent Gaussian random variables with mean 0 and standard deviation σ.sub.∈, to measurement vector h. Eq. (3.7) may be replaced in the previous convex program with Eq. (3.12) and choose ∈.sub.2=√{square root over (m)}σ.sub.∈. The scaling of the reconstruction error ∥W.sup.(opt)−C∥.sub.∞ as a function of σ.sub.∈ is shown in FIG. 8C. The recovery error after the phase transition point scales linearly with σ.sub.∈, consistent with the proposed analytical bounds.

[0180] The convex optimization problem (Eq. 3.13)-(Eq. 3.16) is used to prove rigorous recovery guarantees that show that the optimal solution W.sup.(opt) is close to the true correlation matrix C 400 (FIG. 6), provided that m≥c.sub.0s log n (and with better robustness to noise, when m≥c.sub.0s log.sup.4n). Here, m is the dimension of the measurement vector h, s is the sparsity (the number of nonzero elements) in the off-diagonal part of the matrix C, and c.sub.0 is some universal constant.

[0181] In aspects, two different results are proven: a non-universal recovery guarantee, using the “RIPless” framework, as well as a universal recovery guarantee, using RIP-based techniques. Here, RIP refers to the “restricted isometry property,” a fundamental proof technique in compressed sensing. There are different advantages to the RIPless and RIP-based bounds: the RIPless bounds require slightly fewer measurements, while the RIP-based bounds are more robust when the measurements are noisy.

[0182] In aspects, two variants of the problem (Eq. 3.13)-(Eq. 3.16) are introduced: constrained .sub.1-minimization and the LASSO (“least absolute shrinkage and selection operator”). Generally speaking, recovery guarantees that hold for one of these problems can be adapted to the other one, with minor modifications.

[0183] To simplify the problem, start with the convex optimization problem (Eq. 3.13)-(Eq. 3.16). |Remove the positivity constraint d(W,K.sub.+)≤∈.sub.1. Change the objective function to sum over all i<j rather than all i≠j; since W is symmetric, this merely changes the objective function by a factor of 2. Finally, shift the variable W by subtracting away diag(g), so that its diagonal elements are all zero. Similarly, shift the measurement vector h to get h′=h−Φ(diag(g)). (Eq. 5.1)

[0184] This provides an equivalent problem:

[00036]$\begin{array}{cc}\mathrm{Find}W\in {ℝ}^{n\times n}\mathrm{that}\mathrm{minimizes}{.Math.}_{i<j}.Math.W_{\mathrm{ij}}.Math.,& \left(\mathrm{Eq}.5.2\right)\end{array}$$\mathrm{such}\mathrm{that}:$$\begin{array}{cc}\mathrm{diag}\left(W\right)=0,& \left(\mathrm{Eq}.5.3\right)\end{array}$$\begin{array}{cc}.Math.\Phi \left(W\right)-h^{\prime }{.Math.}_2\le {ϵ}_2+{ϵ}_1\sqrt{\mathrm{mn},}& \left(\mathrm{Eq}.5.4\right)\end{array}$$\begin{array}{cc}W=W^T.& \left(\mathrm{Eq}.5.5\right)\end{array}$

[0185] An operation diag(.Math.) that has two meanings: given an n×n matrix M, diag(M) returns an n-dimensional vector containing the diagonal elements of M; and given an n-dimensional vector v, diag(v) returns an n×n matrix that contains v along the diagonal, and zeros off the diagonal.

[0186] Define C′ to be the off-diagonal part of the correlation matrix C, that is, C′ is the matrix whose off-diagonal elements match those of C, and whose diagonal elements are zero.

C′=C−diag(diag(C)).  (Eq. 5.6)

[0187] h′ may be viewed as a measurement of C′, with additive error z, h′=Φ(C′)+z. (Eq. 5.7)

[0188] The solution W.sup.(opt) is an accurate estimate of C′. The error term z may be written in the form:

[00037]$\begin{array}{cc}\begin{array}{c}z=h^{\prime }-\Phi \left(C^{\prime }\right)\\ =h-\Phi \left(C\right)-\Phi \left(\mathrm{diag}\left(g-\mathrm{diag}\left(C\right)\right)\right)\\ =v-\Phi \left(\mathrm{diag}\left(u\right)\right),\end{array}& \left(5.8\right)\end{array}$

[0189] where u and v are the noise terms in (Eq. 3.9) and (Eq. 3.10). next, bound z using (Eq. 3.10) and (Eq. 3.18), ∥z∥.sub.2≤∈.sub.2+∈.sub.1√{square root over (mn)}. (Eq. 5.9)

[0190] It will be convenient to write

D={W∈.sup.n×n|W.sup.T=W,diag(W)=0}  (Eq. 5.10)

[0191] to denote the subspace of symmetric matrices whose diagonal elements are all 0. Let uvec: D.fwdarw..sup.n(n−1)/2 denote the linear operator that returns the upper-triangular part of W, uvec: W(W.sub.ij).sub.i<j. (Eq. 5.11)

[0192] Write Φ.sub.D: D.fwdarw..sup.m to denote the measurement operator Φ restricted to act on the subspace D. The problem (Eq. 5.2)-(Eq. 5.5) may be rewritten in a more concise form:

[00038]$\begin{array}{cc}\mathrm{Find}W\in D\mathrm{that}\mathrm{minimizes}{.Math.\mathrm{uvec}\left(W\right).Math.}_1,& \left(\mathrm{Eq}.5.12\right)\end{array}$$\mathrm{such}\mathrm{that}:$$\begin{array}{cc}{.Math.{\Phi }_D\left(W\right)-h^{\prime }.Math.}_2\le {ϵ}_2+{ϵ}_1\sqrt{\mathrm{mn}.}& \left(\mathrm{Eq}.5.13\right)\end{array}$

[0193] In aspects, a variant of the problem includes the LASSO:

[0194] Find W∈D that minimizes

[00039]$\begin{array}{cc}\frac{1}{2}.Math.{\Phi }_D\left(W\right)-h^{\prime }{.Math.}_2^2+\lambda .Math.\mathrm{uvec}\left(W\right){.Math.}_1.& \left(\mathrm{Eq}.5.14\right)\end{array}$

[0195] This can be viewed as a Lagrangian relaxation of the previous problem ((Eq. 5.12)-(Eq. 5.13)), or as an .sub.1-regularized least-squares problem. In addition, the convex optimization problems that were described earlier can also be relaxed into a LASSO-like form, in a similar way. The choice of the regularization parameter λ requires some care.

[0196] In aspects, the Regularization Parameter λ may be set by the following procedure. In general, the regularization parameter λ controls the relative strength of the two parts of the objective function in. When the noise in the measurement of h′ is strong, then λ must be set large enough to ensure that the .sub.1 regularization term still has the desired effect. However, if A is too large, it strongly biases the solution W.sup.(opt), making it less accurate.

[0197] An approach to setting λ, may be to have a goal that is to ensure that the solution W.sup.(opt) converges to (a sparse approximation of) the true correlation matrix C. To do this, λ must be set large enough to satisfy two constraints, which involve the noise in the measurement of h′ (see equation (4.1) and the equation below (4.2)). When these constraints are satisfied, the error in the solution W.sup.(opt) is bounded by equation (4.2). (Note that this error bound grows with 2L, hence one should choose the smallest value of λ that satisfies the above constraints.)

[0198] In order to set 2L, first, precise statements of the two constraints on λ are provided:

∥uvec(Φ.sub.D.sup.†z)∥.sub.∞≤λ,  (Eq. 5.15)

∥uvec(Φ.sub.D,T.sub.c.sup.†(I−P)z∥.sub.∞≤λ.  (Eq. 5.16)

[0199] Here, z is the noise term in the measurement of h′ in equation (5.7); Φ.sub.D.sup.†:.sup.m.fwdarw.D is the adjoint of the measurement operator Φ.sub.D; T⊂{(j,j′)|1≤j<j′≤n} is the support of (a sparse approximation of) the true correlation matrix C; Φ.sub.D,T is the sub-matrix of Φ.sub.D that contains those columns of Φ.sub.D whose indices belong to the set T; P is the projection onto the range of Φ.sub.D,T; T.sup.c is the complement of the set T; and Φ.sub.D,T.sub.c is the sub-matrix of Φ.sub.D that contains those columns of Φ.sub.D whose indices belong to the set T.sup.c.

[0200] In order to set λ, compute the quantities in equations (5.15) and (5.16), and to do this, some bounds on the noise z are needed. Two ways of obtaining such bounds include:

[0201] One straightforward way is as follows. Equation (5.8) can be used to write

z=v−Φ(diag(u)),  (Eq. 5.17)

[0202] where u and v are the noise terms in the measurements of g and h, respectively. Bounds on u and v (see equations (3.21) and (3.25)) imply bounds on z, via an elementary calculation. However, one can get better bounds on z by using a more sophisticated approach, starting with bounds on the sub-Gaussian norms of u.sub.j and v.sub.j, such as equations (3.19) and (3.23). Assume that g and h are measured using the procedures described previously. Then equations (3.19) and (3.23) provide bounds on the sub-Gaussian norms of u.sub.j and v.sub.j:

∥u.sub.j∥.sub.104 .sub.2≤δ.sub.1c.sub.jj,  (Eq. 5.18)

∥v.sub.j∥.sub.104 .sub.2≤δ.sub.2Φ.sub.j(C).  (Eq. 5.19)

[0203] Using these bounds, set λ so that it satisfies equations (5.15) and (5.16) with high probability. More precisely, set

[00040]$\begin{array}{cc}\lambda :=\left({ϵ}_1^{\mathrm{\prime ″}}.Math.4\sqrt{\mathrm{mn}}+{ϵ}_2^{\mathrm{\prime ″}}\right)4\sqrt{m}\left(1+\sqrt{2}\right)\sqrt{\ln \left(n\right)/c^{\prime }},& \left(\mathrm{Eq}.5.2\right)\end{array}$$\begin{array}{cc}\mathrm{where}{ϵ}_1^{\mathrm{\prime ″}}:=\frac{{\delta }_1}{1-{ϵ}_1^{″}}{.Math.g.Math.}_{\infty },{ϵ}_1^{″}:=2\sqrt{\ln \left(n\right)/c_0}{\delta }_1,& \left(\mathrm{Eq}.5.21\right)\end{array}$$\begin{array}{cc}{ϵ}_2^{\mathrm{\prime ″}}:=\frac{{\delta }_2}{1-{ϵ}_2^{″}}{.Math.h.Math.}_{\infty },{ϵ}_2^{″}:=2\sqrt{\ln \left(m\right)/c_0}{\delta }_2.& \left(\mathrm{Eq}.5.22\right)\end{array}$

[0204] Here assume that δ.sub.1 and δ.sub.2 are sufficiently small (for instance, δ.sub.1≲1/√{square root over (ln(n))} δ.sub.2≲1/√{square root over (ln(n))} to ensure that ∈″.sub.1<½ and ∈″.sub.2<½. Also, here c′ and c.sub.0 are universal constants (which are defined in the proof below).

[0205] Now let the correlation matrix C and measurement operator Φ.sub.D be fixed, and note that the noise term z is stochastic. With high probability (over the random realization of z), equations (5.15) and (5.16) are be satisfied; here the failure probability is at most (e/n.sup.3)+(e/m.sup.3)+(2e/n.sup.1+2√{square root over (2)}).

[0206] Finally, the following simple upper bounds on ∈′″.sub.1, ∈′″.sub.2 and λ (these follow from the definitions of ∈′″.sub.1 and ∈′″.sub.2, the definitions of g and h, and the definition of Φ(.Math.)):

∈′″.sub.1≤3δ.sub.1∥diag(C)∥.sub.28,  (Eq. 5.23)

∈′″.sub.2≤6δ.sub.2∥C.sub.1,  (Eq. 5.24)

λ≤O(√{square root over (mln(n))}(δ.sub.1∥diag(C)∥.sub.∞√{square root over (mn)}+δ.sub.2∥C.sub.1)).  (Eq. 5.25)

[0207] These bounds may be used to analyze the accuracy of the estimate of matrix C 400 (FIG. 6). Next, isotropy and incoherence of the measurement operator is described. The rows of the measurement operator Φ.sub.D have two properties, isotropy and incoherence, which play a fundamental role in compressed sensing. Let Q be the matrix (of size m by n(n−1)/2) that represents the action of Φ.sub.D (using the fact that the subspace D is isomorphic to .sup.n(n−1)/2); that is, Q and Φ.sub.D are related by the equation: Φ.sub.D(C)=Q.Math.uvec(C),∀C∈D. (Eq. 5.26)

[0208] The rows of Q are chosen independently at random, and each row has the form q=4uvec(rr.sup.T)∈.sup.n(n−1)/2, (Eq. 5.27)

[0209] where r is sampled from the distribution described in (Eq. 3.2). q is centered if it has mean (q)=0, and q is isotropic if its covariance matrix is the identity:

(qq.sup.T)=I.  (Eq. 5.28)

[0210] It is straightforward to check that q is centered and isotropic (up to a normalization

[00041]$\begin{array}{cc}\mathrm{factor}\mathrm{of}2),\mathrm{since}:\{\begin{array}{cc}𝔼\left[r_i{r}_j\right]=0& i<j\\ 𝔼\left[r_i{r}_j{r}_k{r}_l\right]=0& i<j\mathrm{and}k<l\mathrm{and}\left\{i,j\right\}\cap \left\{k,l\right\}=\varnothing \\ 𝔼\left[r_i{r}_j{r}_k{r}_l\right]=0& i<j\mathrm{and}k<l\mathrm{and}.Math.\left\{i,j\right\}\cap \left\{k,l\right\}.Math.=1\\ 𝔼\left[r_i{r}_j{r}_k{r}_l\right]=\frac{1}{4}& i<j\mathrm{and}k<l\mathrm{and}i=k\mathrm{and}j=l\end{array}.& \left(\mathrm{Eq}.5.29\right)\end{array}$

[0211] (Note that in the last line of Eq. 5.29, a case with i=l and j=k, as the requirements of i<j and k<l leads to a contradiction.)

[0212] In addition, q is incoherent with parameter μ>0 if, with probability 1, all of its coordinates are small: ∥q∥.sub.∞.sup.2≤μ. (Eq. 5.30)

[0213] In order for this to be useful for compressed sensing, one needs μ to be small, at most polylogarithmic in the dimension of q. For example, q is incoherent with parameter μ=16.

[0214] A non-universal recovery guarantee, using the “RIPless” framework is described herein, which in turn relies on the isotropy and incoherence properties described above.

[0215] Let C∈.sup.n×n be a correlation matrix, and let C″∈D be its off-diagonal part (see equation (5.6)). C′ is approximately s-sparse, i.e., there exists a matrix C.sup.(s)∈D that has at most s nonzero entries above the diagonal, and that approximates C′ in the (vector) .sub.1 norm, up to an error of size η.sub.s. This can be written compactly as:

∥uvec(C′−C.sup.(s))∥.sub.1≤η.sub.s,  (Eq. 5.31)

[0216] where uvec(.Math.) was defined in equation (5.11). (Recall that both C′ and C.sup.(s) are symmetric, with all zeros on the diagonal. Hence it suffices to consider those matrix entries that lie above the diagonal.)

[0217] Choose the measurement operator Φ at random (see equation (3.4)). Assume that m (the dimension of h) satisfies the bound:

[00042]$\begin{array}{cc}m\ge {\tilde{c}}_0\left(1+\beta \right).Math.4s\log \frac{n\left(n-1\right)}{2}.& \left(\mathrm{Eq}.5.32\right)\end{array}$

[0218] Here, {tilde over (c)}.sub.0 is a universal constant, and β>0 is a parameter that can be chosen freely by the experimenter. Note that m scales linearly with the sparsity s, but only logarithmically with the dimension n of the matrix C. This scaling is close to optimal, in an information-theoretic sense. This is one way of quantifying the advantage of compressed sensing, when compared with measurements that do not exploit the sparsity of C.

[0219] Measure g≈diag(C) and h≈Φ(C), and calculate h′ (see equation (5.1)). Let δ.sub.1 and δ.sub.2 quantify the noise in the measurements of g and h, as described in equations (5.18) and (5.19). Next, solve the LASSO problem, setting the regularization parameter λ according to (Eq. 5.20). Let W.sup.(opt) be the solution of this optimization problem.

[0220] Now the following recovery guarantee is noted, which shows that W.sup.(opt) provides a good approximation to C′, in both the Frobenius (.sub.2) norm, and the vector .sub.1 norm. Theorem 1: For any correlation matrix C satisfying the sparsity condition (5.31), with probability at least

[00043]$1-\frac{12}{n\left(n-1\right)}-6e^{-\beta }$

(over the random choice of the measurement operator Φ, with m set according to (Eq. 5.32)), the solution W.sup.(opt) is close to C′, with an error that is bounded by:

[00044]$\begin{array}{cc}{.Math.W^{\left(\mathrm{opt}\right)}-C^{\prime }.Math.}_F\le \sqrt{2}c\left(1+\alpha \right)\left[\frac{{\eta }_s}{\sqrt{s}}+\lambda \sqrt{s}\right],& \left(\mathrm{Eq}.5.33\right)\end{array}$$\begin{array}{cc}{.Math.W^{\left(\mathrm{opt}\right)}-C^{\prime }.Math.}_{{\ell }_1}\le 2c\left(1+\alpha \right)\left[{\eta }_s+\mathrm{\lambda s}\right],& \left(\mathrm{Eq}.5.34\right)\end{array}$

[0221] where c is a universal constant, and

[00045]$\alpha \le {\log }^{3/2}\left(\frac{n\left(n-1\right)}{2}\right).$

[0222] In these bounds, the first term upper-bounds the error that results from approximating C′ by a sparse matrix, and the second term upper-bounds the error due to noise in the measurements of g and h.

[0223] To make the second term more transparent, the second term can be combined with the bound on λ from equation (5.25):

λ≤O(√{square root over (m ln(n))}(δ.sub.1∥diag(C)∥.sub.∞√{square root over (mn)}+δ.sub.2∥C.sub.1)),  (Eq. 5.35)

[0224] where δ.sub.1 and δ.sub.2 quantify the noise in the measurements of g and h, as described earlier. Also, it is useful to consider the special case where C′ is exactly s-sparse, so η.sub.s=0, and where as few measurement settings as possible are used, by setting m∝s log n. In this case:

∥W.sup.(opt)−C′∥.sub.F≤O(log.sup.3/2(n)λ√{square root over (s)})  (Eq. 5.36)

≤O(√{square root over (s)} log.sup.3/2(n)√{square root over (m log(n))}.Math.[δ.sub.1∥diag(C)∥.sub.∞√{square root over (mn)}+δ.sub.2∥C.sub.1])  ((Eq. 5.37)

≤O(s log.sup.5/2(n)[δ.sub.1∥diag(C)∥.sub.∞√{square root over (sn log(n))}+δ.sub.2∥C.sub.1])  (Eq. 5.38)

≤O(s log.sup.5/2(n)[δ.sub.1∥diag(C)∥.sub.∞√{square root over (sn log(n))}+δ.sub.2√{square root over (n)}∥diag(C)∥.sub.2+δ.sub.2√{square root over (2s)}∥C′∥.sub.F]).  (Eq. 5.39)

[0225] This can be compared with the error bound (Eq. 2.31) for the naive method, and the RIP-based error bound (Eq. 5.48) for compressed sensing. Generally speaking, compressed sensing has an advantage over the naive method when s is small, and the RIPless bound is useful in the regime between m∝s log n and m∝s log.sup.4n, where the RIP-based bound does not apply. (When m is in the regime m∝s log.sup.4n or larger, the RIP-based bound applies, and provides better results than the RIPless bound.)

[0226] Next, a universal recovery guarantee, using an older approach based on the restricted isometry property (RIP) is described. This also relies on the isotropy and incoherence properties shown above.

[0227] First, the number of measurement settings is set as: m≥c.sub.0s log.sup.4n, (5.40) where s is the sparsity parameter, and co is some universal constant. (Note that m is slightly larger, by a poly(log n) factor, compared to the RIPless case.) Also, recall that D is the subspace of symmetric matrices whose diagonal elements are all 0 (see equation (5.10)), and Φ.sub.D is the measurement operator restricted to act on this subspace.

[0228] In aspects, with high probability (over the random choice of Φ.sub.D), the normalized measurement operator Φ.sub.D/√{square root over (m)} satisfies the RIP (for sparsity level 2s). To see this, recall the isotropy and incoherence properties shown above. These properties imply that the measurement operator Φ.sub.D is sampling at random from a “bounded orthonormal system.” Such operators are known to satisfy the RIP, via a highly nontrivial proof.

[0229] From this point onwards, let the measurement operator Φ.sub.D be fixed. Φ.sub.D is capable of reconstructing the off-diagonal parts of all sparse matrices C, i.e., Φ.sub.D can perform “universal” recovery.

[0230] As described above, let C∈.sup.n×n be a correlation matrix, and let C′∈D be its off-diagonal part (see equation (5.6)). Assume that C′ is approximately s-sparse, i.e., there exists a matrix C.sup.(s)∈D that has at most s nonzero entries above the diagonal, and that approximates C′ in the (vector) .sub.1 norm, up to an error of size η.sub.s. This can be written compactly as: ∥uvec(C′−C.sup.(s))∥.sub.1≤η.sub.s, (Eq. 5.41)

[0231] where uvec(.Math.) was defined in equation (5.11). (Recall that both C′ and C.sup.(s) are symmetric, with all zeros on the diagonal. Hence it suffices to consider those matrix entries that lie above the diagonal.)

[0232] Measure g≈diag(C) and h≈Φ(C), and calculate h′ (see equation (5.1)). Assume that the noise in the measurements of g and h is bounded by δ.sub.1 and δ.sub.2, as described in equations (3.19) and (3.23). Then solve the .sub.1-minimization problem in (5.12) and (5.13), setting the parameters ∈.sub.1 and ∈.sub.2 according to equations (3.21) and (3.25). Let W.sup.(opt) be the solution of this problem.

[0233] The following recovery guarantee, which shows that W.sup.(opt) provides a good approximation to C′, in the Frobenius (.sub.2) norm, and in the .sub.1 vector norm. There is one subtle point: the convex optimization problem in equations (5.12) and (5.13) uses the unnormalized measurement operator Φ.sub.D, while the error bounds apply to the normalized measurement operator Φ.sub.D/√{square root over (m)}. Hence, the noise is smaller by a factor of Vi in these error bounds.)

[0234] Theorem 2: With high probability (over the choice of the measurement operator Φ, with m set according to (Eq. 5.40)), for all correlation matrices C (satisfying the sparsity condition (Eq. 5.41)), the solution W.sup.(opt) satisfies the following error bounds:

[00046]$\begin{array}{cc}{.Math.W^{\left(\mathrm{opt}\right)}-C^{\prime }.Math.}_F\le c_1\frac{{\eta }_s}{\sqrt{s}}+c_2\left(\frac{{ϵ}_2}{\sqrt{m}}+{ϵ}_1\sqrt{n}\right),& \left(\mathrm{Eq}.5.42\right)\end{array}$$\begin{array}{cc}{.Math.W^{\left(\mathrm{opt}\right)}-C^{\prime }.Math.}_{{\ell }_1}\le c_1{\eta }_s+c_2\sqrt{s}\left(\frac{{ϵ}_2}{\sqrt{m}}+{ϵ}_1\sqrt{n}\right),& \left(\mathrm{Eq}.5.43\right)\end{array}$

[0235] where c.sub.1 and c.sub.2 are universal constants.

[0236] In these bounds, the first term upper-bounds the error that results from approximating C′ by a sparse matrix, and the second term upper-bounds the error due to noise in the measurements of g and h.

[0237] In order to apply these bounds, one needs to know the values of ∈.sub.1 and ∈.sub.2. These can be obtained from equations (3.22) and (3.26):

∈.sub.1≤O(δ.sub.1∥diag(C)∥.sub.2),  (Eq. 5.44)

∈.sub.2≤O(δ.sub.2√{square root over (m)}∥C.sub.1).  (Eq. 5.45)

[0238] Also, it is useful to consider the special case where C′ is exactly s-sparse, so η.sub.s=0, and where as few measurement settings as possible are used, by setting m∝s log.sup.4n. In this case:

∥W.sup.(opt)−C′∥.sub.F  (Eq. 5.46)

≤O(√{square root over (n)}δ.sub.1∥diag(C)∥.sub.2+δ.sub.2∥C∥.sub.1)  (Eq. 5.47)

≤O(√{square root over (n)}δ.sub.1∥diag(C)∥.sub.2+δ.sub.2(√{square root over (n)}∥diag(C)∥.sub.2+√{square root over (2s)}∥C′∥.sub.F)).  (Eq. 5.48)

[0239] This can be compared with the error bound (Eq. 2.31) for the naive method, and the RIPless error bound (Eq. 5.39) for compressed sensing. Generally speaking, compressed sensing has an advantage over the naive method when s is small, and the RIP-based bound has better scaling (as a function of n, s, diag(C) and C′) than the RIPless bound, although it requires m to be slightly larger.

[0240] The performance of the disclosed compressed sensing method, for a typical measurement scenario is described below. Both the accuracy of the method, and the experimental resources required to implement it are considered. The asymptotic scaling of the disclosed method is described, and compared to the naive method, direct estimation of the correlation matrix.

[0241] Overall, the disclosed compressed sensing method has asymptotically better sample complexity, whenever the off-diagonal part of the correlation matrix C is sufficiently sparse. In particular, for a system of n qubits, the disclosed method is advantageous whenever the number of correlated pairs of qubits, s, is at most O(n.sup.3/2) (ignoring log factors). These results are summarized in FIGS. 8A-8C.

[0242] Let C′ be the off-diagonal part of the correlation matrix C 400, that is, C′ is the matrix whose off-diagonal elements match those of C, and whose diagonal elements are zero. C′ has at most s nonzero elements above the diagonal (and, by symmetry, at most s nonzero elements below the diagonal). The goal is to estimate both C′ and diag(C).

[0243] Compressed sensing allows the possibility of adjusting the number of measurement settings, m, over a range from m∝s log n to m∝n.sup.2. (Note that m∝s log n is just slightly above the information-theoretic lower bound, while m∝n.sup.2 is the number of measurement settings used by the naive method.) Compressed sensing works across this whole range, but the error bounds vary depending on m. There are two cases: (1) For m∝s log.sup.4n or larger, both the RIP-based and RIPless error bounds are available, and the RIP-based error bounds are asymptotically stronger. (2) For m between m∝s log n and m∝s log.sup.4n, only the RIPless error bound is available.

[0244] To make a fair comparison between compressed sensing and the naive method, the accuracy of these methods needs to be quantified in a consistent way. This is a nontrivial task, because the error bounds for the different methods have different dependencies on the parameters n, s, diag(C) and C′; this can be seen by comparing equation (2.31), and Theorems 1 and 2 as described herein.

[0245] A simple way of quantifying the accuracy of all of these methods may be used: given some δ>0, each method is required to return an estimate Ĉ′ that satisfies

∥Ĉ′−C′∥.sub.F≤δ∥C∥.sub.F.  (Eq. 6.1)

[0246] The Frobenius matrix norm may be used, which is equivalent to the vector .sub.2 norm. C (rather than C′) is written on the right hand side of the inequality, in order to allow the recovery error to depend on both the diagonal and the off-diagonal elements of C.

[0247] In addition, each method returns an estimate g of diag(C) that satisfies

∥g−diag(C)∥.sub.2≤δ∥diag(C)∥.sub.2.  (Eq. 6.2)

[0248] For both compressed sensing as well as the naive method, g is obtained in the same way, by performing single-qubit spectroscopy as in (2.19), and the error in g satisfies the same bound (Eq. 2.30).

[0249] The cost of implementing each method using real experiments should be accounted for. This cost depends on a number of factors. One factor is the total number of experiments that have to be performed, often called the sample complexity. This is the number of measurement settings, times the number of repetitions of the experiment with each measurement setting. Another factor is the difficulty of performing a single run of the experiment. This involves both the difficulty of preparing entangled states (random n-qubit GHZ states for the compressed sensing method, and 2-qubit Bell states for the naive method), and the length of time that one has to wait in order to observe dephasing.

[0250] An advanced quantum information processor is considered, where n-qubit GHZ states are fairly easy to prepare (using O(logn)-depth quantum circuits), and dephasing occurs at low rates, so that the main cost of running each experiment is the amount of time needed to observe dephasing. In this scenario, it is reasonable to use the sample complexity as a rough measure of the total cost of implementing the compressed sensing method, as well as the naive method.

[0251] In aspects, the sample complexity for three methods of interest may be calculated: (1) the naive method with m∝n.sup.2, (2) compressed sensing with m∝s log.sup.4n, and (3) compressed sensing with m∝s log n. This method (2) outperforms the naive method whenever s≤O(n.sup.3/2/log.sup.2n), and method (3) outperforms the naive method whenever s≤O(n.sup.2/3/log.sup.2n).

[0252] In addition, for each of these methods, the number of samples where single-qubit spectroscopy is performed, and the number of samples where multi-qubit spectroscopy is performed. (Recall that all of these methods use single-qubit spectroscopy to estimate the diagonal of C, and then use multi-qubit spectroscopy to estimate the off-diagonal part of C.) Both of these numbers can be important: multi-qubit spectroscopy is more expensive to implement on essentially all experimental platforms, and requires more samples when s is large; but it is possible for single-qubit spectroscopy to dominate the overall sample complexity, when s is small.

[0253] A naive method with m∝n.sup.2 is presented. Two parameters, δ.sub.1 and δ.sub.2, to quantify the accuracy of the measurements are used, as in equations (2.24) and (2.25). An estimate of C′ is provided, whose error is bounded by equation (2.31): with probability at least 1−η,

[00047]$\begin{array}{cc}{.Math.{\hat{C}}^{\prime }-C^{\prime }.Math.}_F^2\le \frac{1}{\eta }\left[3\left(n-1\right){\left({\delta }_1+{\delta }_2\right)}^2{.Math.\mathrm{diag}\left(C\right).Math.}_2^2+6{\delta }_2^2{.Math.C^{\prime }.Math.}_F^2\right].& \left(\mathrm{Eq}.6.3\right)\end{array}$

[0254] For simplicity, η is set to be some universal constant, for example, η=0.001. Now, given any δ>0:

[00048]$\begin{array}{cc}\begin{array}{c}{.Math.{\hat{C}}^{\prime }-C^{\prime }.Math.}_F^2\le {\delta }^2{.Math.\mathrm{diag}\left(C\right).Math.}_2^2+\left({\delta }^2/n\right){.Math.C^{\prime }.Math.}_F^2\\ \le {\delta }^2{.Math.C.Math.}_F^2\end{array},& \left(\mathrm{Eq}.6.4\right)\end{array}$

[0255] by setting δ.sub.1=δ.sub.2=O(δ/√{square root over (n)}). This satisfies the requirement (Eq. 6.1).

[0256] In addition, one can easily check that the estimate g for diag(C) satisfies the requirement (Eq. 6.2).

[0257] Then the sample complexity is as follows (see the discussion preceding (Eq. 2.24) and (Eq. 2.25)): the method performs single-qubit spectroscopy on O(n/δ.sub.1.sup.2)=O(n.sup.2/δ.sup.2) samples, and multi-qubit spectroscopy on O(n.sup.2/δ.sub.2.sup.2)=O(n.sup.3/δ.sup.2) samples. Hence the total sample complexity is O(n.sup.3/δ.sup.2).

[0258] In aspects, compressed sensing may be performed with m∝s log.sup.4n. The solution W.sup.(opt) for the .sub.1-minimization problem in (Eq. 5.12) and (Eq. 5.13) satisfies the RIP-based bound in Theorem 2. Two parameters, δ.sub.1 and δ.sub.2, are used to quantify the accuracy of the measurements, as in equations (3.19) and (3.23).

[0259] The estimator W.sup.(opt) satisfies the following error bound (see Theorem 2, and equation (5.48)):

∥W.sup.(opt)−C′∥.sub.F≤O(√{square root over (n)}δ.sub.1∥diag(C)∥.sub.2+δ.sub.2(√{square root over (n)}∥diag(C)∥.sub.2+√{square root over (2s)}∥C′∥.sub.F)).  (Eq. 6.5)

[0260] Now, given any δ>0:

[00049]$\begin{array}{cc}\begin{array}{c}{.Math.W^{\left(\mathrm{opt}\right)}-C^{\prime }.Math.}_F\le \frac{1}{\sqrt{2}}\delta {.Math.\mathrm{diag}\left(C\right).Math.}_2+\frac{1}{\sqrt{2}}\delta {.Math.C^{\prime }.Math.}_F\\ \le \delta {.Math.C.Math.}_F\end{array},& \left(\mathrm{Eq}.6.6\right)\end{array}$

[0261] by setting δ.sub.1=O(δ/√{square root over (n)}) and δ.sub.2=O(δ/√{square root over (max(n,s))}). (Eq. 6.7)

[0262] This satisfies the requirement (Eq. 6.1). In addition, one can easily check that the estimate g for diag(C) satisfies the requirement (Eq. 6.2). Then the sample complexity is as follows (see the discussion following (Eq. 3.19) and (Eq. 3.23)): the method performs single-qubit spectroscopy on O(n/q)=O(n.sup.2/δ.sup.2) samples, and multi-qubit spectroscopy on

[0263] O(m/δ.sub.2.sup.2)≤O(smax(n,s)log.sup.4(n)/δ.sup.2) (Eq. 6.8) samples. Hence the total sample complexity is at most O(max(n,s).sup.2 log.sup.4(n)/δ.sup.2). (Eq. 6.9)

[0264] This is less than the sample complexity of the naive method, provided the off-diagonal part of the correlation matrix is sufficiently sparse, i.e., when s≤O(n.sup.3/2/log.sup.2n).

[0265] In aspects, compressed sensing may be performed with m∝s log n. The LASSO optimization problem whose solution W.sup.(opt) satisfies the RIPless bound in Theorem 1. Two parameters, δ.sub.1 and δ.sub.2, are used to quantify the accuracy of the measurements, as in equations (3.19) and (3.23). In the following, the diagonal elements of C satisfy a bound of the form

[00050]$\begin{array}{cc}{.Math.\mathrm{diag}\left(C\right).Math.}_{\infty }\le 0\left(\frac{1}{\sqrt{n}}{.Math.\mathrm{diag}\left(C\right).Math.}_2\right).& \left(\mathrm{Eq}.6.1\right)\end{array}$

[0266] This assumption may be used to get a stronger error bound for W.sup.(opt). The assumption (Eq. 6.10) indicates that none of the diagonal elements c.sub.jj is too much larger than the others. This is plausible for a quantum system 100 that consists of many qubits that are constructed in a similar way.

[0267] In order to make this intuition more precise, (Eq. 6.10) may be written in an equivalent form:

[00051]$\begin{array}{cc}\underset{1\le j\le n}{\max }\left(c_{\mathrm{jj}}^2\right)\le 0\left(\frac{1}{n}{.Math.}_{j=1}^n{c}_{\mathrm{jj}}^2\right),& \left(\mathrm{Eq}.6.11\right)\end{array}$

which indicates that the largest c.sub.jj.sup.2 is at most a constant factor larger than the average of all of the c.sub.jj.sup.2. Also, it is informative to consider how (Eq. 6.10) and (Eq. 6.11) compare to the (arguably more natural) assumption that

[00052]$\begin{array}{cc}\underset{1\le j\le n}{\max }.Math.c_{\mathrm{jj}}.Math.\le 0\left(\frac{1}{n}{.Math.}_{j=1}^n.Math.c_{\mathrm{jj}}.Math.\right).& \left(\mathrm{Eq}.6.12\right)\end{array}$

[0268] In fact, (Eq. 6.12) is actually a stronger assumption, in the sense that it implies (Eq. 6.10) and (Eq. 6.11), via the Cauchy-Schwarz inequality. The estimator W.sup.(opt) satisfies an error bound given by Theorem 1, and equation (5.39). Combining this with the assumption (Eq. 6.10):

∥W.sup.(opt)−C′∥.sub.F≤O(s log.sup.5/2(n)[δ.sub.1√{square root over (s log(n))}∥diag(C)∥.sub.2+δ.sub.2√{square root over (n)}∥diag(C)∥.sub.2+δ.sub.2√{square root over (2s)}∥C′∥.sub.F]).  (Eq. 6.13)

[0269] Now, given any S>0, ensure that

[00053]$\begin{array}{cc}\begin{array}{c}{.Math.W^{\left(\mathrm{opt}\right)}-C^{\prime }.Math.}_F\le \frac{1}{\sqrt{2}}\delta {.Math.\mathrm{diag}\left(C\right).Math.}_2+\frac{1}{\sqrt{2}}\delta {.Math.C^{\prime }.Math.}_F\\ \le {\delta }^2{.Math.C.Math.}_F^2\end{array},& \left(\mathrm{Eq}.6.14\right)\end{array}$$\begin{array}{cc}\mathrm{by}\mathrm{setting}{\delta }_1=O\left(\frac{\delta }{s^{3/2}{\log }^3\left(n\right)}\right),& \left(\mathrm{Eq}.6.15\right)\end{array}$$\begin{array}{cc}\mathrm{and}{\delta }_2=O\left(\frac{\delta }{s{\log }^{5/2}\left(n\right)\sqrt{\max \left(n,s\right)}}\right).& \left(\mathrm{Eq}.6.16\right)\end{array}$

[0270] This satisfies the requirement (Eq. 6.1).

[0271] In addition, one can easily check that the estimate g for diag(C) satisfies the requirement (Eq. 6.2). Then the sample complexity is as follows (see the discussion following (3.19) and (3.23)). The disclosed method performs single-qubit spectroscopy on O(n/δ.sub.1.sup.2)=O(ns.sup.3 log.sup.6(n)/δ.sup.2) (Eq. 6.17)

[0272] In aspects, the disclosed method further includes performing single-qubit spectroscopy samples, and multi-qubit spectroscopy on O(m/δ.sub.2.sup.2)≤O(s.sup.3 max(n,s)log.sup.6(n)/δ.sup.2) (Eq. 6.18) samples. Hence the total sample complexity is at most O(s.sup.3 max(n,s)log.sup.6(n)/δ.sup.2). (Eq. 6.19)

[0273] The total sample complexity is less than the sample complexity of the naive method, provided the off-diagonal part of the correlation matrix is sufficiently sparse, i.e., when s≤O(n.sup.2/3/log.sup.2n)

[0274] Referring to FIGS. 9A and 9B, an evolution time t is chosen. During the physical implementation of the measurements of the correlation matrix C, decay rates Γ.sub.ab may be estimated. To do this, quantum states |ψ.sub.ab are prepared, allowing the quantum states to evolve for some time t, and then measure the quantum states in an appropriate basis. This works well when t is chosen appropriately, so that Γ.sub.abt˜1.

[0275] For example, one method of estimating the decay rate may use an evolution time t such that

[00054]$\frac{1}{2}\le {\Gamma }_{ab}t\le 2.$

An initial guess for t (τ.sub.0) is made. Next, a “binary search” is performed. The binary search includes running a sequence of experiments, where one observes the dephasing of the state |ψ.sub.ab for some time t, and after each experiment, one adjusts the time t adaptively, multiplying and dividing by factors of 2, in order to get the “right amount” of dephasing. Generally, the sequence of experiments is about ˜|log(Γ.sub.abτ.sub.0)| experiments. More precisely, the following procedure may be performed:

[0276] 1. Fix some τ.sub.0>0; as the initial guess for the evolution time t.

[0277] 2. For r=1, 2, . . . , N.sub.trials, perform the following: (N.sub.trials are set according to equation (7.9) below)

[0278] (a) Set s.sub.0=0 and t.sub.0=2.sup.s.sup.0τ.sub.0 (as the initial guess for t).

[0279] (b) For j=0,1,2, . . . , N.sub.steps−1, do the following: (N.sub.steps are set according to equation (7.5) below)

[0280] i. Prepare the state

[00055]$.Math.{\psi }_{ab}.Math.=\frac{1}{\sqrt{2}}(.Math.a.Math.+.Math.b.Math.),$

allow the state to dephase for time then measure in the basis

[00056]$\frac{1}{\sqrt{2}}(.Math.a.Math.\pm .Math.b.Math.)$

[0281] ii. If the measurement returns

[00057]$\frac{1}{\sqrt{2}}(.Math.a.Math.+.Math.b.Math.),$

then set

[00058]$\begin{array}{cc}{s}_{j+1}=\{\begin{array}{c}{s}_{j+1}\mathrm{with}\mathrm{probability}\frac{e-1}{e+1},\\ s_j\mathrm{otherwise}.\end{array}& \left(\mathrm{Eq}.7.1\right)\end{array}$

[0282] If the measurement returns

[00059]$\frac{1}{\sqrt{2}}(.Math.a.Math.-.Math.b.Math.),$

then set s.sub.j+1=s.sub.j−1.

[0283] iii. Set t.sub.j+1=2.sup.s.sup.j+τ.sub.0. (This is the next guess for t.)

[0284] (c) Define τ.sub.r to be the value of s.sub.j from the last iteration of the loop, that is, ξ.sub.r=s.sub.N.sub.steps.

[0285] 3. Compute the average

[00060]$\xi =\frac{1}{N_{\mathrm{trials}}}{.Math.}_{r=1}^{N_{\mathrm{trials}}}{\xi }_r.$

Return {circumflex over (t)}=2.sup.86 τ.sub.0. (This is the estimate for t.)

[0286] This procedure can be described intuitively as follows. The inner loop of this procedure (the loop indexed by j) can be viewed as a kind of stochastic gradient descent, which behaves like a random walk on real numbers of the form t=2.sup.sτ.sub.0 (s∈) (see the dashed curves in FIG. 9A).

[0287] In aspects, this random walk has a single basin of attraction at a point t*=2.sup.s*τ.sub.0 that satisfies Γ.sub.abt*≈1, that is, s*≈−log.sub.2(Γ.sub.abτ.sub.0) The random walk converges to this point: with high probability, the sequence s.sub.0, s.sub.1, s.sub.2, . . . will reach the point s* after O(|s*|)=O(|log(Γ.sub.abτ.sub.0)|) steps; after that point, the sequence will remain concentrated around s*, with exponentially-decaying tail probabilities (see FIG. 9B). This claim is made precise in equations (7.5) and (7.6).

[0288] Finally, the outer loop of this procedure (the loop indexed by r) computes an estimate ξ of s*, by averaging over several independent trials (see the solid curves in FIG. 9A). This then yields an estimate {circumflex over (t)} of t*. The required number of trials, and the accuracy of the resulting estimate {circumflex over (t)}, are analyzed in equations (7.9) and (7.10).

[0289] Referring to FIGS. 9A and 9B, plots illustrating the evolution time t are shown. FIG. 9A shows the trajectories of the random walk. Start with an initial guess τ.sub.0, that can either be shorter or longer than the optimal time. The evolution time is then stochastically halved or doubled over N.sub.steps iterations, according to the algorithm described herein. This procedure is repeated (dashed curves) and the outcome is averaged (solid curves) to obtain an estimate {circumflex over (t)}. The region in which

[00061]$\frac{1}{2}<{\Gamma }_{\mathrm{ab}}\hat{t}<2$

is shaded in green. FIG. 9B illustrates the accuracy of the final estimate t, as a function of the number of steps N.sub.steps and the starting point τ.sub.0. The green shading shows the region where {circumflex over (t)} satisfies the bound

[00062]$\frac{1}{2}<{\Gamma }_{\mathrm{ab}}\hat{t}<2.$

At the boundary of the green region, N.sub.steps scales logarithmically with |Γ.sub.abτ.sub.0|, as predicted by Eq. (7.5).

[0290] To choose t the random walk may be used. The variables s.sub.j, which are related to the t.sub.j via the identity t.sub.j=2.sup.s.sup.jτ.sub.0. It is easy to see that s.sub.0=0, s.sub.j+1∈{s.sub.j,s.sub.j−1,s.sub.j+1}, and

[00063]$\begin{array}{cc}\begin{array}{c}𝔼\left[s_{j+1}.Math.s_j\right]=s_j+\frac{1}{2}\left(1+e^{-{\Gamma }_{\mathrm{ab}}{t}_j}\right)\frac{e-1}{e+1}-\frac{1}{2}\left(1-e^{-{\Gamma }_{\mathrm{ab}}{t}_j}\right)\\ =s_j+\frac{1}{e+1}\left(e^{1-{\Gamma }_{\mathrm{ab}}{t}_j}-1\right).\end{array}& \left(\mathrm{Eq}.7.2\right)\end{array}$

[0291] Hence the sequence s.sub.0, s.sub.1, s.sub.2, . . . can be viewed as the trajectory of a random walk on a 1-dimensional chain, beginning at s.sub.0, with transition probabilities that vary along the chain. The expected behavior of the random walk can be bounded as follows:

[00064]$\begin{array}{cc}𝔼\left[s_{j+1}.Math.s_j\right]\{\begin{array}{c}\ge s_j+\mu \mathrm{when}{\Gamma }_{\mathrm{ab}}{t}_j\le \frac{1}{\sqrt{2}},\\ \le s_j-\mu \mathrm{when}{\Gamma }_{\mathrm{ab}}{t}_j\ge \sqrt{2},\\ \approx s_j+\frac{e-1}{e+1}\mathrm{when}{\Gamma }_{\mathrm{ab}}{t}_j\ll 1,\\ \approx s_j-\frac{1}{e+1}\mathrm{when}{\Gamma }_{\mathrm{ab}}{t}_j\gg 1,\end{array}& \left(\mathrm{Eq}.7.3\right)\end{array}$

[0292] where μ is a numerical constant,

μ=0.09.  (Eq. 7.4)

[0293] Hence the random walk will tend to converge towards some integer s* such that t*=2.sup.s*τ.sub.0 satisfies

[00065]$\frac{1}{\sqrt{2}}\le {\Gamma }_{\mathrm{ab}}{t}^{*}\le \sqrt{2}.$

[0294] Note that the expected position of the random walk moves towards s* at a rate that is lower-bounded by μ. So, in order to go from s.sub.0=0 to s*≈−log.sub.2(Γ.sub.abτ.sub.0), the random walk will take roughly

[00066]$\frac{1}{\mu }.Math.s^{*}.Math.=\frac{1}{\mu }.Math.{\log }_2\left({\Gamma }_{\mathrm{ab}}{\tau }_0\right).Math.$

steps. It is easy to see that the stationary distribution of the random walk is centered around s*, with exponentially decaying tails; hence, once the walk reaches s*, it will remain concentrated around that point.

[0295] Set the parameter N.sub.steps so that, with high probability, after N.sub.steps steps, the random walk will converge. Given an upper bound h on the magnitude of s*, i.e., |s*|≤h, or equivalently 2.sup.−h≤Γ.sub.abτ.sub.0≤2.sup.h. The random walk may be performed for a number of steps

[00067]$\begin{array}{cc}{N}_{\mathrm{steps}}=\frac{h}{\mu }+\eta ,& \left(\mathrm{Eq}.7.5\right)\end{array}$

[0296] where η≥0. Here,

[00068]$\frac{h}{\mu }$

is (an upper bound on) the expected number of steps needed to reach s*. Take an additional η steps to ensure that the walk does indeed reach s* with high probability (the probability of failure decreases exponentially with η).

[0297] After N.sub.steps steps, the final position of the walk is close to s*, with exponentially decaying tail probabilities: for any ≥1,

[00069]$\begin{array}{cc}\mathrm{Pr}[.Math.s_{N_{\mathrm{steps}}}-s^{*}.Math.\ge \ell .Math.s_0=0]\le \frac{16}{{\mu }^2}\exp \left(-\frac{\mu \left(\mu +1\right)}{8}\ell +\frac{{\mu }^2}{4}\right)+2\exp \left(-\frac{\mu }{16}\min \left\{\frac{\mathrm{\mu \eta }}{h},1\right\}\left(\ell +\mathrm{\mu \eta }\right)\right).& \left(\mathrm{Eq}.7.6\right)\end{array}$

[0298] In particular, when

[00070]$\eta \ge \frac{h}{\mu },$

this bound can be slightly simplified:

[00071]$\begin{array}{cc}\mathrm{Pr}\left[.Math.s_{N_{\mathrm{steps}}}-s^{*}.Math.\ge \ell |s_0=0\right]\le \frac{16}{{\mu }^2}\exp \left(-\frac{\mu \left(\mu +1\right)}{8}\ell +\frac{{\mu }^2}{4}\right)+2\exp \left(-\frac{\mu }{16}\left(\ell +\mathrm{\mu \eta }\right)\right).& \left(\mathrm{Eq}.7.7\right)\end{array}$

[0299] The bound (Eq. 7.6) is proved using martingale techniques.

[0300] In aspects, the parameter N.sub.trials may be set, and an error bound on the estimator {circumflex over (t)} may be derived. The bound (Eq. 7.6) implies that the random variables ξ.sub.r=s.sub.N.sub.steps are sub-exponential, and their sub-exponential norms are bounded by some constant ∥ξ.sub.r∥.sub.ψ.sub.1≤K. Hence their average ξ satisfies a Bernstein-type concentration inequality: for every δ≥0,

[00072]$\begin{array}{cc}\mathrm{Pr}\left[.Math.\xi -s^{*}.Math.\ge \delta \right]\le 2\exp \left(-c\min \left\{\frac{{\delta }^2}{K^2},\frac{\delta }{K}\right\}{N}_{trials}\right),& \left(\mathrm{Eq}.7.8\right)\end{array}$

[0301] where c>0 is a universal constant.

[0302] Now, for any ∈>0, set

[00073]$\begin{array}{cc}{N}_{\mathrm{trials}}=\frac{1}{c}\max \left\{\frac{K}{\delta },\frac{K^2}{{\delta }^2}\right\}\log \left(\frac{2}{ϵ}\right).& \left(\mathrm{Eq}.7.9\right)\end{array}$

[0303] The following error bound on the estimator {circumflex over (t)}=2.sup.ξτ.sub.0: with probability ≥1−∈, provides |ξ−s*|<δ, which implies that

2.sup.−0.5−δ≤Γ.sub.ab√{square root over (t)}<2.sup.0.5αδ,  (Eq. 7.10)

[0304] Assuming δ<½, this implies that

[00074]$\frac{1}{2}<{\Gamma }_{ab}\hat{t}<2,$

as desired.

[0305] When the disclosed measurement protocols are implemented in an experiment, errors may occur during state preparation and measurement (SPAM errors). The effect of these errors on estimating the decay rates Γ.sub.ab may be investigated. Let ρ.sub.0 and E.sub.0 denote the noiseless initial state and observable of interest, respectively.

[00075]$\begin{array}{cc}{\rho }_0=E_0=\frac{1}{2}(.Math.a.Math..Math.a.Math.+.Math.b.Math..Math.b.Math.+.Math.b.Math..Math.a.Math.+.Math.a.Math..Math.b.Math.).& \left(8.1\right)\end{array}$

[0306] Error channels ε.sub.s and ε.sub.m that act on state preparation and measurement operations may be considered as:

{tilde over (ρ)}=ε.sub.s(ρ.sub.0)=ρ.sub.0+δρ  (Eq. 8.2)

{tilde over (E)}=ε.sub.m(E.sub.0)=E.sub.0+δE,  (Eq. 8.3)

[0307] where ∥δρ∥.sub.tr≤∈.sub.s and ∥δE∥≤∈.sub.m, and ∈.sub.s and ∈.sub.m are small parameters. The outcome of the protocol is now given by

{tilde over (P)}.sub.ab(t)=Tr[{tilde over (E)}ε.sub.t({tilde over (ρ)})],  (Eq. 8.4)

[0308] where ∈.sub.t=exp(t) is the evolution under the correlated dephasing noise (1.1).

[0309] The protocol is robust against these kinds of errors, and for short times t the decay of {tilde over (P)}ab(t) is still dominated by Γ.sub.ab. Using Eqs. (8.2) and (8.3):

Tr[{tilde over (E)}ε.sub.t({tilde over (ρ)})]=Tr[E.sub.0ε.sub.t(ρ.sub.0)]  (Eq. 8.5)

+Tr[E.sub.0ε.sub.t(δρ)]+Tr[δEε.sub.t(ρ.sub.0)]

+Tr[δEε.sub.t(δρ)]

[0310] The first term is the outcome without errors:

[00076]$\begin{array}{cc}\mathrm{Tr}\left[E_0{ℰ}_t\left({\rho }_0\right)\right]=\frac{1+e^{-t{\Gamma }_{ab}}}{2}.& \left(\mathrm{Eq}.8.6\right)\end{array}$

[0311] Find the effect of errors on the second and third terms, by considering the effect of ε.sub.t on ρ.sub.1 and E.sub.1. Specifically:

Tr[E.sub.0ε.sub.t(δρ)]=η.sub.s+ζ.sub.se.sup.−Γ.sup.ab.sup.t,  (Eq. 8.7)

Tr[δEε.sub.t(ρ.sub.0)]=η.sub.m+ζ.sub.me.sup.−Γ.sup.ab,  (Eq. 8.8)

[0312] where η.sub.s,m and ζ.sub.s,m are constants that are determined by δρ and δE, for s and m, respectively. Therefore, these terms decay at the same rate as the first case and do not affect the exponential decay. However, the last term can, in principle, contain different decay rates and can cause deviation from a single exponential decay. Bound the rate at which R(t)=Tr[δEε.sub.t(δρ)] grows:

[00077]$\begin{array}{cc}.Math.\dot{R}\left(t\right).Math.=.Math.\frac{\partial }{\partial t}\mathrm{Tr}\left[\delta E{ℰ}_c\left(\delta \rho \right)\right].Math.& \left(\mathrm{Eq}.8.9\right)\end{array}$$\begin{array}{cc}=.Math.\mathrm{Tr}\left[\delta Eℒ\left(\mathrm{\delta \rho }\right)\right].Math.\le 2{ϵ}_m{ϵ}_s\left(n+s\right),& \left(\mathrm{Eq}.8.1\right)\end{array}$

[0313] Therefore:

[00078]$\begin{array}{cc}{\tilde{P}}_{ab}\left(t\right)=\frac{1}{2}\left(1+{\eta }_s+{\eta }_m+\left(1+{\zeta }_s+{\zeta }_m\right)e^{-{\Gamma }_{\mathrm{ab}}t}\right)+R\left(t\right),& \left(\mathrm{Eq}.8.12\right)\end{array}$

[0314] The deviations from a single exponential decay are attributed to R(t). Using Eqs. (8.11) and (8.12) the decay rate of {tilde over (P)}.sub.ab(t) is dominated by Γ.sub.ab for evolution times t≲1/(2∈.sub.s∈.sub.m(n+s)).sup.−1.

[0315] Referring to FIGS. 10A-10C, graphs illustrating the exponential decay of {tilde over (P)}.sub.ab(t) when there are SPAM errors is shown. The decay of a 3-qubit GHZ state is simulated. A correlation matrix C is selected with uniform single-qubit dephasing rate (c.sub.ii=γ.sub.0) and nearest-neighbor correlations c.sub.i,i±1=γ.sub.0/4. FIGS. 10A-10C show the decay of {tilde over (P)}.sub.ab(t) under different noise strengths: (FIG. 10A) Δ=0.01, (FIG. 10B) Δ=0.05, (FIG. 10C) Δ=0.1. The dashed lines show the decay with no SPAM errors, and solid lines show the decay with SPAM errors from randomly-sampled error channels. The solid lines resemble the dashed lines for short evolution times t.

[0316] Numerical simulations may be used to investigate the effect of SPAM errors on estimates of the decay rate Γ.sub.ab. Simulate SPAM errors by applying random error channels ε.sub.s and ε.sub.m, whose strengths are controlled by a parameter Δ. Compare the decay of {tilde over (P)}.sub.ab(t) with and without SPAM errors, for different values of Δ. Observe that, for short times t, the decay with SPAM errors matches the decay without SPAM errors, i.e., the decay rate is dominated by Γ.sub.ab, see FIGS. 10A-10C. This is consistent with the theoretical analysis.

[0317] The disclosed method for learning sparse correlated dephasing noise (FIGS. 3 and 4) can be extended to the most general case of the master equation (1.1), where the matrix C is complex, and there is an additional Hamiltonian term H.sub.s.

[0318] In aspects, the decay rates may be complex. The complete dynamics imposed by the environment on the quantum system 100 can have a coherent evolution in addition to the decay. The evolution of the quantum system

[00079]$100\frac{d\rho }{\mathrm{dt}}=ℒ\left(\rho \right)$

is in general given by the Lindblad generator

[00080]$\begin{array}{cc}ℒ\left(\rho \right)=-i\left[\rho ,H_s\right]+{.Math.}_{l,m}{c}_{lm}\left(Z_l\rho Z_m-\frac{1}{2}\left\{Z_l{Z}_m,\rho \right\}\right),& \left(\mathrm{Eq}.9.1\right)\end{array}$$\begin{array}{cc}\mathrm{where}{H}_s={.Math.}_{l,m}{r}_{lm}{Z}_l{Z}_m& \left(\mathrm{Eq}.9.2\right)\end{array}$

[0319] is sometimes called the (generalized) Lamb shift Hamiltonian, and C=(c.sub.lm) is now a complex matrix. Matrix C may be decomposed as C=V+iT, (Eq. 9.3)

[0320] where the real and imaginary parts are separated into V=(v.sub.lm), a real symmetric matrix, and T=(t.sub.lm), a real skew-symmetric matrix, respectively. Moreover, the Lamb shift Hamiltonian can be encoded in the symmetric matrix R=(r.sub.lm).

[0321] In aspects, the operator acts on the off-diagonal matrix elements |ab|. This involves “decay rates” that depend on R, V and T in a simple way, although these “decay rates” are now complex rather than real. Specifically,

(|ab|)=(−Γ.sub.ab+iΩ.sub.ab)|ab|,  (Eq. 9.4)

[0322] where Γ.sub.ab and Ω.sub.ab are real numbers that capture the decay and oscillations of the matrix element, respectively. The convention is defined for states |a and |b: Z.sub.i|a=α.sub.i|a, and similarly for |b and β.sub.i. Therefore, separating the real and imaginary parts of (|ab|) it is found that:

[00081]$\begin{array}{cc}{\Gamma }_{ab}=-{.Math.}_{l,m}{v}_{lm}\left({\alpha }_l{\beta }_m-\frac{1}{2}{\alpha }_l{\alpha }_m-\frac{1}{2}{\beta }_l{\beta }_m\right)=\frac{1}{2}{\left(\alpha -\beta \right)}^{T}V\left(\alpha -\beta \right),& \left(9.5\right)\end{array}$

[0323] The fact that V is symmetric can be used. Similarly:

[00082]$\begin{array}{cc}{\Omega }_{ab}=-{.Math.}_{l,m}{r}_{lm}\left({\alpha }_l{\alpha }_m-{\beta }_l{\beta }_m\right)+{.Math.}_{l,m}{v}_{lm}\left({\alpha }_l{\beta }_m-\frac{1}{2}{\alpha }_l{\alpha }_m-\frac{1}{2}{\beta }_l{\beta }_m\right)=-\left({\alpha }^{T}R\alpha -{\beta }^{T}R\beta \right)+\frac{1}{2}\left({\alpha }^{T}T\beta -{\beta }^{T}T\alpha \right).& \left(\mathrm{Eq}.9.6\right)\end{array}$

[0324] R and T, are symmetric and skew-symmetric matrices, respectively. Therefore, the coherences, that is, the |ab| elements of the density matrix ρ(t), exhibit both oscillations (at a frequency nab) and exponential decay (at a rate Γ.sub.ab) as a function of t. Note that it is possible to measure both Γ.sub.ab and Ω.sub.ab, by estimating the matrix elements of ρ(t) at different times t, using the same types of Ramsey experiments described herein. In particular, one can extract these parameters using standard techniques for analyzing spectral lines in atomic physics. Here, the squared magnitude of the Fourier transform of the measurement time series is a Lorentzian function, the center of the Lorentzian peak provides the oscillation frequency, and the width of the peak provides the decay rate.