Quantum process tomography via completely positive and trace-preserving projection
George C. Knee Eliot Bolduc Jonathan Leach Erik M. Gauger



[email: ]gk@physics.org
Department of Physics, University of Warwick, Coventry, CV4 7AL, UK

SUPA, Institute of Photonics and Quantum Sciences, Heriot-Watt University, David Brewster Building, Edinburgh, EH14 4AS, UK

SUPA, Institute of Photonics and Quantum Sciences, Heriot-Watt University, David Brewster Building, Edinburgh, EH14 4AS, UK

SUPA, Institute of Photonics and Quantum Sciences, Heriot-Watt University, David Brewster Building, Edinburgh, EH14 4AS, UK

                                           


We present an algorithm for projecting superoperators onto the set of completely positive, trace-preserving maps. When combined with gradient descent of a cost function, the procedure results in an algorithm for quantum process tomography: finding the quantum process that best fits a set of sufficient observations. We compare the performance of our algorithm to the diluted iterative algorithm as well as second-order solvers interfaced with the popular cvx package for matlab, and find it to be significantly faster and more accurate while guaranteeing a physical estimate. 

Quantum process tomography via completely positive and trace-preserving projection
    Erik M. Gauger
    December 30, 2023
==================================================================================



§ INTRODUCTION

As experimental quantum information science progresses, researchers are increasingly turning to the characterisation of quantum processes (see Fig <ref>a) rather than the quantum states which they act upon. For example, the fault-tolerance theorem <cit.> (which underpins the feasibility of large scale quantum computers) places a requirement on the accuracy of operations performed on quantum bits. In fact, the theorem relies on a worst-case analysis with respect to state preparation (before the operation) and with respect to measurement (after the operation) in judging their accuracy. Fortunately, to evaluate such metrics one does not need to actually prepare and measure the worst case: A complete description of the process can be inferred from a suitable set of preparations and measurements and subsequently any property of interest may then be calculated. Such inference is known as quantum process tomography (QPT).  QPT has been performed in many different physical settings <cit.>: the challenge is producing an estimate of the process in a reasonable amount of time which matches the data as closely as possible and is consistent with a quantum mechanical model of the experiment. We present algorithms for reconstructing the process subject to appropriate constraints, via a subroutine that implements a composite projection onto the set of quantum channels– see Fig <ref>b.

In Sec. <ref>, we describe the maximum likelihood approach to process tomography, along with the Choi representation of processes and individual projections onto important constraint sets. In Sec. <ref> we describe our main results, including the composite projection onto quantum channels, and two algorithms that exploit it to estimate a quantum process from data. Our algorithms are benchmarked against existing approaches in Sec. <ref>, and conclusions are drawn in Sec <ref>. Further details are given in Appendices. 





§ PROCESS TOMOGRAPHY


A quantum process is a linear `superoperator' which acts on an input density operator ρ to produce an output density operator ρ':

    ℰ(ρ)=ρ'.

Here ρ and ρ' are order 2 tensors (i.e. operators), which we assume (for simplicity) to act on Hilbert spaces of equal dimension d. ℰ is an order 4 tensor specified by d^4-d^2 parameters, where d is itself exponential in the number of subsystems. Assumption-free QPT is therefore necessarily very expensive: cheaper alternatives are made possible by relying on prior knowledge (e.g. compressive sensing <cit.>) or by estimating a single summary statistic (as in randomised benchmarking <cit.>), but such alternatives may be insufficient to decide whether ℰ  meets fault tolerance requirements <cit.>. It is therefore important to reduce the costs associated with assumption-free QPT while maintaining accuracy.



 §.§ Estimation of processes
The probability of observing an outcome corresponding to E_j (a positive measurement operator), when the quantum process has transformed some input state ρ_i, is

    p_ij = Tr( ℰ(ρ_i)E_j).

For each input state, N_i samples are drawn and the number of outcomes n_ij that correspond to each outcome j are recorded so that ∑_j n_ij=N_i.  	In order to produce an estimate for a quantum process from such a set of observations, a cost function (or measure of `un-fitness') f(ℰ) is chosen and then minimised over a set of candidate processes. The minimisation is handled by, for example, gradient descent steps of the cost function. Maximum-likelihood (ML) estimation is a principled and prominent choice, and one that we make here in order to showcase our general method. The likelihood of ℰ given the `data' n_ij is
ℒ = ∏_ij p_ij^n_ij
up to an irrelevant constant. According to the principle of ML, we seek ℰ^*, the quantum process that maximises ℒ. This approach is justified heuristically by considering ℒ as a measure of agreement between the model p_ij and the data n_ij, assuming a multinomial statistical model. ML enjoys the key property of asymptotic efficiency: as the number of trials N_i becomes large, under mild conditions <cit.> the variance of the ML estimator is lower than any unbiased estimator, saturating the Cramér-Rao bound <cit.>. Alternative estimators may be preferable for finite N_i <cit.>, but will generally perform worse asymptotically in terms of their precision and accuracy. In the limit where n_ij are all large, the likelihood may be approximated by a Gaussian function and the problem reduces to minimising the cost function f(ℰ)→∑_ij(p_ij-n_ij)^2 <cit.>, which may be solved by linear inversion <cit.>. This is an unjustified simplification in general, but may be transformed into an effective heuristic, as we will show below.  Regardless of the choice of cost function, our approach is flexible in that it is not tied (as some other methods are <cit.>) to a particular choice of (preparation and/or) measurement operators. Proceeding with the (convex) cost function

    f(ℰ)   =-lnℒ=-∑_ij n_ijln p_ij,

we will introduce below an algorithm minimising f subject to appropriate constraints. The minimiser of f is also the maximiser of ℒ due to the monotonicity of the logarithm.  Unconstrained tomographic procedures usually produce unphysical quantum states or processes <cit.>, often due to an implicit and idealised assumption of zero noise. A physical estimate is paramount if one is to use it in any further theoretical analysis, such as calculating purity, fidelity, entanglement, expectation values and so on. We will now introduce the constraints on ℰ and our approach to enforcing them, which uses a certain representation of ℰ.


 §.§ Choi representation

It is well known that quantum processes, considered as superoperators with dimension d× d × d × d may be represented as d^2× d^2-dimensional operators on a Hilbert space ℋ_in⊗ℋ_out. These are known as Choi operators, defined as:

    C_ℰ = ∑_ij |i⟩⟨ j  | ⊗ℰ(|i⟩⟨ j  |)

for an orthonormal set of kets {|i⟩} which form a basis of both the input and output Hilbert spaces. This relation implies that ℰ may be inferred by applying it to one part of a d^2-dimensional system prepared in a single maximally entangled state, and performing quantum state tomography on the output: a technique known as ancilla-assisted QPT <cit.>.  Here we continue with the standard approach which uses multiple input states in d dimensions, and considers the complete dataset as a whole (rather than simply performing state tomography on each output state independently <cit.>). One can easily verify that ℰ(ρ)=Tr_ℋ_in[(ρ^T⊗𝕀)C_ℰ] which implies that the forward model can be rewritten as

    p_ij = Tr([ ρ_i^T  ⊗ E_j]C_ℰ).

Now the cost function has a new argument f(ℰ)→ f(C_ℰ). The gradient of f with respect to this new argument follows from operator calculus 

    ∇ f (C_ℰ) =∂  f (C_ℰ)/∂ C_ℰ= -∑_ijn_ij/p_ijρ_i ⊗ E_j^T.

These expressions can be vectorised to enable faster evaluation on a computer – see Appendix <ref>.



 §.§ Constraints and projections

Just as some d× d operators fail to qualify as proper density matrices (quantum states), ℰ must satisfy certain constraints in order for it to represent a proper quantum process. Firstly, it should be completely positive (CP). This means that it must preserve the positive semidefiniteness of an arbitrary input state when acting on only a part of a larger space of arbitrary dimension: ρ≽ 0⇒[ℰ⊗𝕀_d'](ρ)≽ 0   ∀ d'.  A theorem due to Choi <cit.> states that the complete positivity of ℰ is equivalent to the positive semidefiniteness of C_ℰ [A positive semidefinite matrix is denoted by the expression C≽ 0, meaning that C belongs to the subset of Hermitian matrices with nonnegative eigenvalues.]:

    ℰ∈CP⇔ C_ℰ≽0.

Since CP is a convex set, we can find the closest CP map to any given Hermitian superoperator by using a projection of its Choi matrix onto the set of positive semidefinite operators. This is achieved via

    𝒞𝒫(𝒞_ℰ)=Vmax(D,0)  V^†

with 𝒞_ℰ=VDV^†: i.e. through removing negative terms in the eigendecomposition of 𝒞_ℰ. 𝒞𝒫 is then a projection orthogonal under the Frobenius norm <cit.>.
Secondly, it is common to assume the process satisfies the property of trace preservation (TP): Tr(ℰ(ρ))=Tr(ρ)  ∀ρ. In the Choi operator picture, this is equivalent to the condition 

    Tr_ℋ_out(C_ℰ)=𝕀,

and the TP set is therefore defined by an affine constraint <cit.>. In Appendices  <ref> and <ref>, we (respectively) explain why the TP constraint makes QPT a non-trivial extension of state tomography, and discuss relaxing the constraint to trace non-increasing (probabilistic) processes. We recast Eq. (<ref>) as M[C_ℰ] = [𝕀], defining

    M = ∑_k 𝕀⊗⟨ k | ⊗𝕀⊗⟨ k|

The vec[·] operation reshapes its argument into a vector, and we denote the inverse operation by vec^-1[·]. 
The orthogonal projection onto TP is the solution to

    min_C1/2 ||[C]-[C_ℰ]||^2   s.t.   M[C] = [𝕀]

given by


    𝒯𝒫(C_ℰ)    =^-1[[C_ℰ]- M^†(MM^†)^-1(M[C_ℰ]-[𝕀])]
       =^-1[[C_ℰ]- 1/d M^† M [C_ℰ] +1/dM^†[𝕀]]
    .


In the second line, we used MM^†=d 𝕀⊗𝕀; it is useful to note that M^† M = ∑_ij𝕀⊗ |j⟩⟨ i|⊗𝕀⊗|j⟩⟨ i|. 


§ COMPOSITE PROJECTION


It is not obvious that the ability to separately project into the CP and TP sets of superoperators enables projection into the intersection set CPTP.  While 𝒞𝒫 and 𝒯𝒫 are separately projections, their sequential action will not generally result in a matrix belonging to CPTP (unlike the equivalent sequential actions for projection onto quantum states <cit.>). This may be intuited with the help of the geometrical picture in Fig. <ref>b.
Repeated averaged projections, however, defined (at iteration k of the outer loop of our algorithm) by

    H^0 = C^k_ℰ;      H^l+1 = 1/2(𝒯𝒫(H^l)+𝒞𝒫(H^l))

will converge (in linear time <cit.>) to a point in CPTP – see Fig. <ref>b. Dykstra's alternating projection algorithm is a superior alternative, which we choose because it achieves projection 𝒞𝒫𝒯𝒫(C_ℰ) onto the closest point in the intersection of the sets <cit.>. A similar method was used in Ref <cit.> to solve a feasibility problem where input and output states are given and an exact but non-unique solution is sought. By contrast the ML approach (which we adopt) finds the best possible fit to some given, noisy data and treats the more realistic situation where an exact solution cannot be expected. Pseudo-code for this subroutine (showing our use of a robust stopping criterion due to Birgin and Raydan <cit.>) is given below. 




 §.§ Proposed algorithm: pgdB

Returning to the main problem of QPT, we are now ready to apply the principle of projected gradient descent. Such an approach has recently been shown to offer speed benefits when applied to quantum state tomography <cit.>. Applying it to QPT implies taking a single gradient descent step before running an alternating projection <cit.> under 𝒞𝒫 and 𝒯𝒫 until this `inner loop' converges, before taking another gradient descent step and repeating. Formally the ML-QPT problem is written

C_ℰ∈ℂ^d^2×d^2f=-∑_ij n_ij ln( Tr([ ρ_i^T  ⊗E_j]C_ℰ) )

C_ℰ≽0  (CP constraint)
Tr_ℋ_out(C_ℰ)=𝕀_d×d  (TP constraint),

which our proposed algorithm solves via the iterative update rule:

    C^0_ℰ   =𝕀_d^2× d^2/d;
         C_ℰ^k+1   =𝒞𝒫𝒯𝒫(C_ℰ^k-1/μ∇ f(C_ℰ^k)).

Here μ is a step-size metaparameter. At each iteration a step is taken in the locally downhill direction of f, with the aim of increasing the likelihood, before a projection ensures that the constraints are satisfied (up to any desired tolerance) – see Fig <ref>c. Pseudocode of the full algorithm pgdb is given below, and shows our implementation of backtracking (Armijo line search) to improve convergence times <cit.>.




It is important to point out that, although by our choice of f the considered problem is convex (meaning there are no local optima which are not global optima) there may nevertheless be regions (where p_ij=0 for some i,j) where f has an undefined gradient and where the likelihood is stationary <cit.>. These subspaces correspond to those processes which take ρ_i to a state orthogonal to E_j. The probability of a randomly chosen process lying within these sub-dimensional regions is low. However, an iteration of our algorithm might encounter such a region, causing the algorithm to stall: this is especially the case since projections tend to result in rank deficient matrices. Our solution to this issue is presented in Appendix <ref>. 

Note that our approach can be adapted to solve other constrained optimisation problems over quantum processes: whenever the cost function f is convex and continuously differentiable the projected gradient descent procedure is guaranteed to converge to an optimal solution <cit.>.



 §.§ Linear inversion with a final projection (LIFP)
 To show that our CPTP projection has applications in QPT other than for ML, we present a further algorithm based on linear inversion of a vectorised form of Eq. (<ref>), (p_ij) =A(C_ℰ), where the rows of A are composed of [ρ_i⊗ E_j^T]^T. To allow for a comparison with the manifestly physical approach of the ML methods, we upgrade naïve linear inversion so as to force it to also produce physical estimates — simply by applying our novel, composite CPTP projection just once. We call the resultant estimator Linear Inversion with a Final Projection:

    C^LIFP_ℰ:= 𝒞𝒫𝒯𝒫(  vec^-1 [A^+n]),

where A^+=A^†(AA^†)^-1 is the pseudo-inverse of A <cit.> and n=[n_ij]. This approach is in the same spirit as the `quick and dirty' state tomography approach used by Kaznady and James <cit.>, suggested as an alternative to full-blown optimization. For full details, see Appendix <ref>.



§ NUMERICAL BENCHMARKING


Unlike other approaches that require an orthonormal operator basis <cit.> ours is flexible enough to deal with any informationally complete set of preparations and measurements. To test the performance of our proposed algorithm, we considered the reconstruction of qudit processes, with the d^2 preparations ρ_i=|ψ_i⟩⟨ψ_i| projecting onto the following pure states:

    |ψ_i⟩ = 
    
    |j⟩    j=1… d
    |j⟩+|k⟩/√(2)    j=1… d,k>j
    |j⟩+√(-1)|k⟩/√(2)    j=1… d,k>j

and a Positive Operator Valued Measure (POVM) formed with the 2d^2 elements

    E_j = 
    ρ_i/d^2     i=1… d^2 
    𝕀/d^2-ρ_i/d^2     i=1… d^2

which resolve to the identity ∑_jE_j=𝕀. This choice of preparations and measurements is a very practical one, when compared with the use of symmetric informationally complete POVMs (SIC-POVMs) <cit.> or ancillas <cit.> which require much greater experimental effort. Appendix <ref> shows how the choice of measurement and preparation operators affects the conditioning of the optimisation problem. We ignore systematic errors in state preparation and measurement (so-called SPAM errors): gate set tomography <cit.> is a recently proposed solution for a self-consistent characterization in the presence of such imperfections.

We generated random, purity >90% CPTP maps C_ℰ_true in d=2,…,8, simulating multinomial random data with sample size N_i=N ranging from 10 to 10^9 (see Appendix <ref>). Throughout our simulations we normalise the data such that ∑_j n_ij=1. Although this of little consequence for Maximum Likelihood (ML) and other optimisation approaches (the cost function is merely scaled), it is an important step for Linear Inversion and LIFP, since the latter approach is based on treating the data as probabilities – hence the need to normalise them into frequencies rather than raw counts.

We also simulated the infinite data case by setting the frequencies equal to the theoretical probabilities N→∞⇒ n_ij=p_ij.  Although not a realistic situation this is a useful device to compare various iterative QPT algorithms as the limit of many observations is approached, where the reconstruction error is dominated by numerical errors rather than noise in the data. We then ran our algorithm (as well as benchmark algorithms) to find C_est. We recorded running times and accuracies for each. The figure of merit for accuracy was taken to be the J distance, defined as:

    J(ℰ_est,ℰ_true)= || 𝒞_ℰ_est-𝒞_ℰ_true||_tr / 2d.

Here ||·||_tr is the trace norm, or sum of singular values, and J(ℰ_est,ℰ_true)∈[0,1] is related to the average probability of distinguishing the two processes <cit.>. Simulations were run using matlab R2017a and cvx v2.1 on a Intel Core i7-4790 3.6 GHz with 8 MB L3 cache. 



 §.§ Benchmark 1 : Diluted iterations

The state-of-the-art algorithm for ML-QPT is the technique of diluted iterations (dia), devised by Fiurasek and Hradil in 2001 <cit.>, who adapted an algorithm for quantum state tomography known as diluted Rρ R <cit.>. The main idea is to exploit an extremal equation obeyed by the ML-CPTP map:

    C_ℰ=Λ^-1(C_ℰ)∇ f (C_ℰ)C_ℰ∇ f (C_ℰ)Λ^-1(C_ℰ),

where Λ(C_ℰ)=(Tr_ℋ_out(C_ℰ))^1/2⊗𝕀 is derived by incorporating the TP constraint through a Lagrange multiplier. This equation is treated as an iterative update rule. The algorithm has been been applied to real data, mostly with Fock-space truncated optical (i.e. continuous variable) systems, according to the formulation by Anis and Lvovsky <cit.>. Cooper et al. reported that the algorithm ran for 6.5 hours (with a machine precision stopping rule) to reconstruct a d=6 conditional state engineering process <cit.>. Fedorov et al. performed QPT of a beamsplitter, revealing the Hong-Ou-Mandel effect <cit.>. The algorithm typically requires diluting with a parameter 0 < ϵ≤ 1 to prevent f from oscillating:

    R^k    = ϵ∇ f(C_ℰ^k)+(1-ϵ)𝕀,
    Λ^k    = (Tr_ℋ_out(R^kC_ℰ^kR^k))^1/2⊗𝕀,
    
    C_ℰ^k+1   =(Λ^k)^-1R^kC_ℰ^kR^k(Λ^k)^-1.

In our application of dia, we used a crude optimisation of ϵ at each iteration to prevent overshoots. Pseudocode is shown below:





 §.§ Benchmark 2: mosek

The problem (<ref>) is straightforward to enter into the cvx modelling environment <cit.>, which can then call one of a number of general-purpose solvers.  cvx solves problems featuring logarithmic cost functions with a successive approximation method. We found SeDuMi failed to find a solution and that sdpt3 succeeded but was slow. We therefore settled on the commercial solver mosek: it usually produced good results but sometimes failed in higher dimensions. We estimate the failure probability to be as high as 15%, see Appendix <ref>.



§ RESULTS AND CONCLUSIONS
 

Our results are summarised in Fig. <ref>: The CPTP-projection-based algorithms are not only faster than existing approaches, but significantly more accurate. This is especially the case as d increases but is true even for d=4 (corresponding, for example, to a 2-qubit gate such as the controlled-NOT). 

We confirmed that due to the iterative nature of pgdB, higher accuracies can be achieved by adjusting the stopping criterion and other metaparameters, at the expense of a longer algorithm runtime.  LIFP is a quick and effective method suited to situations with low noise and low d, and sacrificing a small amount of accuracy, typically less than an order of magnitude, with respect to the iterative methods. Appendix <ref> includes results corresponding to reconstruction of full rank processes, which show qualitatively similar trends. 

Our results indicate that CPTP projection is a valuable and versatile tool for QPT, holding great promise for applications where time and accuracy are both important. matlab code is available at <https://github.com/geoknee/CPTPprojection>.


 G.C.K was supported by the Royal Commission for the Exhibition of 1851 and thanks Yanbao Zhang for discussions. E.M.G. thanks the Royal Society of Edinburgh and the Scottish Government for support. We thank Yen-Huan Li for useful feedback on this document. 
 



§ VECTORIZING THE PROBLEM


Our forward model may be vectorised as
p := [p_ij]=A[C_ℰ]
upon defining

    A = 
    (
    [ [ρ_1^*⊗ E_1^†]^†;                ⋮; [ρ_1^*⊗ E_N^†]^†;                ⋮; [ρ_N^*⊗ E_N^†]^† ])=
    (
    [ [ρ_1⊗ E_1^T]^T;              ⋮; [ρ_1⊗ E_N^T]^T;              ⋮; [ρ_N⊗ E_N^T]^T ]).

Now f(C_ℰ)=-n^T ln(p) and ∇ f(C_ℰ) = - A^†η with η_ij = n_ij/p_ij and η=[η_ij] and so on. These facts follow from elementary matrix calculus along with the identity Tr(A^† B) ≡[A]^†[B]. Shang et al. present a method to speedup calculation of the gradient when A has a tensor product structure <cit.>: unfortunately this does not help in our case because we have only a single tensor product ρ_i^T⊗ E_j and multiple tensor products are needed to show an advantage. In a slight abuse of notation, above we have denoted the vectorization of a matrix M by [M_ij]. where M_ij are the matrix elements of M.


§ TRACE PRESERVATION CONSTRAINT



Counter to a common misconception, the celebrated `Choi-Jamiolokowski isomorphism' <cit.> does not imply a 1:1 correspondence between CPTP maps and density operators C_ℰ/Tr(C_ℰ) in a higher dimensional space. As discussed in the main paper, 
some such density operators do not correspond to CPTP maps <cit.>. See for example Refs <cit.>. As an example of why a stronger condition is important, consider the Choi matrix C_=diag(0.1,0.1,0.1,1.7). This is proportional to a perfectly valid, classically correlated 2-qubit state. It is manifestly positive, and has trace equal to 2. Note however, that Tr_2(C_)=diag(0.2,1.8)≠𝕀. This means that it does not represent a trace preserving map. Note that the given example actually increases the trace of some states and decreases the trace of others so is neither in TP nor in TNI (discussed below).



§ PROJECTION ONTO THE SET OF TRACE NON-INCREASING PROCESSES


It will be possible to relax our algorithmic projections in such a way as to search in larger spaces, i.e. supersets of TP. This may be important for certain applications: in fact, whenever there is loss or some other non deterministic process such as a measurement, CPTP maps (or `channels') are replaced by CPTNI maps (or  `operations') which are trace non-increasing. The trace of the output density matrix corresponds to the probability of success of the map. CP maps admit a Kraus representation <cit.>:

    ℰ(ρ)=∑_i K_iρ K_i^†.

The difference in trace between an input and an output state is given by

    Tr(ℰ(ρ))-Tr(ρ)   =Tr(∑_i K_iρ K_i^†)-Tr(ρ)
       =Tr([𝒴 -𝕀]ρ),

where we used the cyclic property of the trace and defined 𝒴=∑_i K_i^† K_i. Since ρ is an arbitrary positive operator, it is clear that trace preservation is equivalent to

    𝒴=𝕀    (TP).

 In fact, it is always possible to diagonalise 𝒴 with a similarity transformation 𝒱. Let λ_i be the eigenvalues of 𝒴 <cit.>. Then
 
    Tr(ℰ(ρ))-Tr(ρ) = Tr(diag[λ_1-1,λ_2-1,…]𝒱^†ρ𝒱).

So a general CP map will change the trace of an input state depending on the eigenvalues of 𝒴, and on the projection of the state onto the eigenvectors of 𝒴.
We can consider the case where the process has a uniform and known success probability p:

    𝒴=p𝕀    (US_p).

The projection onto this set is easy: simply alter the TP projection given in Eq. (<ref>) by taking [𝕀]→ p[𝕀]. 
 If we demand that the trace is non-increasing for an arbitrary input state, we have

    𝒴≼𝕀    (TNI).


A commonly used approach to reconstruct TNI operations involves introducing of a fictitious `shelving' state (thereby extending the space d→ d+1 <cit.>). One then performs a TP reconstruction in the larger space before projecting out the fictitious state. 

We present an alternative here that avoids (the possibly very expensive downside of) having to increase the size of the Hilbert space. By the definition of the Choi matrix, 

    Tr_ℋ_out(C_ℰ)=𝒴.

It is straightforward to project onto TNI (which is a convex superset of TP). The correlations χ in the Choi matrix are implicitly defined through

    C_ℰ = Tr_ℋ_out(C_ℰ)⊗𝕀+𝕀⊗Tr_ℋ_in(C_ℰ)+χ.

We can achieve the TNI projection via 


    𝒯𝒩ℐ(C_ℰ)     =
    C_ℰ+1/d[𝒢(Tr_ℋ_out(C_ℰ))-Tr_ℋ_out(C_ℰ)]⊗𝕀.


where we introduced 𝒢, defined and computed via

    𝒢(X)    = min_Tr_ℋ_out(B)≼𝕀||Tr_ℋ_out(B)⊗𝕀-X⊗𝕀||
       = min_Tr_ℋ_out(B)≼𝕀||Tr_ℋ_out(B)-X||· ||𝕀||
       = min_Tr_ℋ_out(B)≼𝕀||Tr_ℋ_out(B)-X||
       = min_E≼𝕀||E-X||
       = Vmin(D,1)  V^†

with X=VDV^† <cit.>. To see that we have the orthogonal projection, consider that for any B∈ TNI we may write B =Tr_ℋ_out(B)⊗𝕀+𝕀⊗Tr_ℋ_in(B)+χ_B. Now


    ||B-ℰ||    = ||  Tr_ℋ_out(B-C_ℰ)⊗𝕀+𝕀⊗Tr_ℋ_in(B-C_ℰ)+χ_B-χ||
       ≤ ||  Tr_ℋ_out(B-C_ℰ)⊗𝕀||+||𝕀⊗Tr_ℋ_in(B-C_ℰ)||+||χ_B-χ||
       ≤ ||  Tr_ℋ_out(B-C_ℰ)⊗𝕀||≤min_Tr_ℋ_out(B)≼𝕀 ||  Tr_ℋ_out(B)⊗𝕀-Tr_ℋ_out(C_ℰ)⊗𝕀||
       = min_Tr_ℋ_out(B)≼𝕀 ||  Tr_ℋ_out(B)-Tr_ℋ_out(C_ℰ)||· ||𝕀||=min_Tr_ℋ_out(B)≼𝕀1/d||  Tr_ℋ_out(B)-Tr_ℋ_out(C_ℰ)||
       = ||𝒯𝒩ℐ(C_ℰ)-C_ℰ||,


and therefore 𝒯𝒩ℐ(C_ℰ) is the closet element in TNI to C_ℰ. We verified this numerically using cvx.
Using the techniques shown in the main text, this result enables tomography of quantum operations belonging to CPTNI without expanding the reconstruction space.


§ STALLING



As mentioned in the main text, it is possible for pgdB to stall if a situation with p_ij=0 is encountered, and in this section we discuss a simple solution allowing the algorithm to continue. Importantly, our proposed solution alerts the user to when the stalling fix has been applied, turning such a rare occurrence into a heralded event. It is important to note that other potential solutions also exist, but a comprehensive study of their merits and drawbacks (including bounding any systemic inaccuracies that arise) are beyond the scope of the present work.

One straightforward solution to keeping the algorithm running when a p_ij=0 situation arises is to modify the cost function and gradient at each iteration with the heralded conditioning step:

    p_ij=max(p_ij,ϵ)

for ϵ some small parameter. It is clearly possible that the algorithm never encounters p_ij<ϵ; in that case our conditional modification does not occur and we proceed toward the ML solution. In the case that the modification is triggered, the algorithm is capable of raising a flag, warning the user to interpret the resultant solution accordingly – it would also be possible for the flag to trigger a restart of the algorithm with a random initial guess. In theory, the ML solution cannot have any p_ij=0 (since that is a zero of the likelihood, which is positive semidefinite). The exception is when p_ij=n_ij=0, trivial cases excluded by taking the convention 0^0=1. Therefore 

    p_ij=
    
    p_ij,       p_ij>0
    ϵ,       p_ij=0

is sufficient to keep the algorithm going without altering the turning points of the objective function. In practice, machines with finite precision cannot calculate the logarithm of very small numbers, necessitating the stronger conditioning step given in Eq. (<ref>). The implication is that, while the ML solution cannot be at p_ij=0, it can have p_ij<ϵ. Therefore with low probability the conditioning step can affect the accuracy of the reconstructed solution (but this possibility is heralded). The probability of this occurring is related to the Kraus rank of the process – see Fig. <ref>, which shows that the conditioning step prevents the stalling problem with no impact on the reconstruction error in most cases. 


§ LINEAR INVERSION WITH A FINAL PROJECTION ONTO CPTP



The most straightforward approach to inverting the forward model is arguably the approach of linear inversion. Here, the estimate is 

    C^LI_ℰ:=  vec^-1 [A^+n],

where A^+=A^†(AA^†)^-1 is the pseudo-inverse of A <cit.>, and A is defined in Eq. (<ref>). This estimator minimises the distance in Frobenius norm between the data and the forward model – a familiar notion of `least squares' fitting. This notion of discrepancy is distinct from the negative likelihood used in the main paper. Note however, that in the special case where the likelihood is a Gaussian function, or indeed when there is no noise at all, the two notions of discrepancy coincide <cit.>. Linear Inversion is a non-iterative approach so is likely to be fast in low dimensions, although the number of floating point operations will be O(d^12) and will therefore not scale as well as pgdB (which is O(d^8)). Furthermore without the explicit expression of the CP and TP constraints, there is no guarantee of a physical estimate. 

We implemented this estimator using the matlab backslash operator \, which avoids explicit calculation of the pseudo-inverse and instead uses a method of Gaussian elimination. We projected the estimate onto the space of Hermitian matrices C^LI_ℰ→(C^LI_ℰ+(C^LI_ℰ)^†)/2 in order to investigate their subsequently real eigenspectrum.  We found that estimates were unphysical with overwhelming probability – see Fig. <ref>. The severity of the problem increases with shot noise. It is therefore not possible to use the standard analytical tools of infidelity or J distance to meaningfully compare the accuracy of linear inversion with ML. It is worth noting that the ML approach provides estimates with smallest eigenvalue no less than some negative quantity which can be set arbitrarily small. 


To allow for a comparison with the manifestly physical approach of ML, we upgraded linear inversion to force it to also produce physical estimates — simply by reapplying the novel, composite CPTP projection derived in the main paper just once. We call the resultant estimator Linear Inversion with a Final Projection:

    C^LIFP_ℰ:= 𝒞𝒫𝒯𝒫(  vec^-1 [A^+n]),

once more forgoing explicit calculation of A^+ in favour of the matlab backslash operator \. This approach is very much in the same spirit as the `quick and dirty' state tomography approach used by Kaznady and James <cit.>. We supplement the results in the main paper with an investigation in d=7, shown in Fig. <ref>. The results confirm the idea that linear inversion with a final projection is a quick and effective method when the number of trials is taken to infinity, but that the loss in accuracy for finite N is exacerbated in higher dimensions. In fact, the projection is not necessary when there is zero noise. When noise is negligible, linear inversion provides a very high precision estimate due to its non-iterative nature, but our results suggest that ML approaches this situation faster as multinomial noise is decreased. Otherwise, in the physically relevant scenario of finite N a significant loss in precision is clear to see. The running time of Linear Inversion with a Final Projection can be favourable in low dimensions but will become much slower than the other methods as d increases.



§ CONDITION NUMBER OF A



Fig <ref> shows the condition number of A rising with d, for the `minimal' choice of preparations and measurements made in the main text. A high condition number means a slower runtime and a less accurate result in general. The condition number might be brought down by considering over-complete preparations and measurements, or SIC-POVMs <cit.>.  The density of A is typically around 10 percent or less, meaning efficiency savings for sparse operations with matlab. 

The computational complexity of the gradient-based algorithms is dominated by the matrix-vector product required to calculate the gradient. The complexity is 𝒪(d^4n_combs) where n_combs is the number of combinations of preparations and measurement outcomes. Our choice above implies n_combs=2d^4. 
An over-complete set would impact running time: for example, choosing the eigenvectors of the generalised Gell-Mann matrices <cit.> implies n_combs=(d^3-d)^2, and raises the complexity from 𝒪(d^8) to 𝒪(d^10).  



§ FAILURE OF CVX



Depending on dimensionality and statistical noise, we found the experimental method of successive approximation used by cvx <cit.> to be unreliable. This is a known issue and to be expected due to the heuristic nature of the method. Fig. <ref> shows the failure rate to be as high as 15%. The data for mosek used in Figs. <ref>, <ref> and <ref> are post-selected on the algorithm succeeding.



§ GENERATION OF FULL-RANK AND QUASI-PURE RANDOM CPTP MAPS


Random CPTP maps may be uniformly randomly generated according to a prescription by Bruzda et al. <cit.>:

    Y    = (Tr_ℋ_out[XX^†])^1/2,
    
    B   =(𝕀⊗ Y^-1/2)XX^†(𝕀⊗ Y^-1/2),

where X is a d^2× M complex random matrix with entries distributed normally. The Kraus rank of the CPTP map is M. We generated d^2 rank-1 matrices B_i and formed their convex combination C_ℰ_true=∑_iP_iB_i with an exponentially decaying probability distribution P_i with  ∑_i P_i^2=0.9. The resulting ensemble is CPTP (by convexity) and has purity Tr[C_ℰ^2]/d^2 no less than 90%. We call this the `quasi-pure' ensemble.


Using the first part of the above procedure we generated an alternative data set corresponding to full rank processes. The results are shown in Fig. <ref>, where the lower purity of this ensemble (compared to the quasi-pure ensemble of Fig. <ref>) is apparent from the relative closeness of a typical Choi matrix to diag(1,1,… 1)/d.  Our proposed algorithm is still superior to existing methods, although less dramatically so. 





