{'fraction_non_alphanumeric': 0.03508644169009161, 'fraction_numerical': 0.010457478262136064, 'mean_word_length': 4.428955483513877, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Compressive gate set tomography

Flexible characterization techniques that identify and quantify experimental imperfections under realistic assumptions are crucial for the development of quantum computers. Gate set tomography is a characterization approach that simultaneously and self-consistently extracts a tomographic description of the implementation of an entire set of quantum gates, as well as the initial state and measurement, from experimental data. Obtaining such a detailed picture of the experimental implementation is associated with high requirements on the number of sequences and their design, making gate set tomography a challenging task even for only two qubits. In this work, we show that low-rank approximations of gate sets can be obtained from significantly fewer gate sequences and that it is sufficient to draw them randomly. Such tomographic information is needed for the crucial task of dealing with coherent noise. To this end, we formulate the data processing problem of gate set tomography as a rank-constrained tensor completion problem. We provide an algorithm to solve this problem while respecting the usual positivity and normalization constraints of quantum mechanics by using second-order geometrical optimization methods on the complex Stiefel manifold. Besides the reduction in sequences, we demonstrate numerically that the algorithm does not rely on structured gate sets or an elaborate circuit design to robustly perform gate set tomography. Therefore, it is more flexible than traditional approaches.


I. INTRODUCTION
The precise characterization of digital quantum devices is crucial for several reasons: (i) to obtain 'actionable advice' on how imperfections on their implementation can be reduced, e.g. by experimental control, (ii) to tailor applications to unavoidable device errors so that their effect can be mitigated, and (iii) to benchmark the devices for the comparison of different physical platforms and implementations. There is already a wide variety of protocols to characterize components of a digital quantum computing device with a trade off between the information gained about the system and the associated resource requirements and assumptions of the scheme [1,2].
One particular important requirement for practical characterization protocols for quantum gates is their robustness against errors in the state preparation and measurement (SPAM). There are two general approaches that SPAM-robustly characterize the implementation of entire gate sets of a quantum computer. On the low complexity side there is randomized benchmarking (RB) [3][4][5] and variants thereof [6], that typically aim at determining a single measure of quality for an experiment, though with the exception of RB tomography protocols [7][8][9][10]. Yet for the targeted improvement of individual quantum operations, protocols which provide more detailed information beyond mere benchmarking are crucial. This is the motivation of self-consistent gate set tomography (GST) [11][12][13][14][15][16][17][18]. GST estimates virtually all parameters describing a noisy implementation of a quantum computing device simultaneously from the measurements of many different gate sequences. This comprises tomographic estimates for all channels implementing the gates of the gate set, the initial state(s) and the measurement(s). Concomitant with the massive amount of inferred information and minimal assumptions, these protocols come with enormous resource requirements in terms of the necessary number of measurement rounds and the time and storage consumption of the classical post-processing. Standard GST, as described, e.g. by Nielsen et al. [18], uses many carefully designed gate sequences in the experiment and a sophisticated and challenging data processing pipeline in post. To arrive at physically interpretable estimates, i.e. completely positive and trace preserving (CPT) maps, additional postprocessing is required. The massive amount of specific data consumed by standard GST limits its practical applicability already for two-qubit gate sets.
It is worth noting that the focus of Nielsen et al. [18] and their implementation pyGSTi lies on so called long sequence GST, a method to improve an initial GST estimate by using gate sequences in which a building block is repeated many times. The resulting error amplification of the builing block is then used to significantly improve the accuracy of the GST estimate. In our work we focus on short sequences and and the problem of finding an initial estimate without assuming any prior knowledge on the gate set. We expect that progressively longer sequences can be added at the end of our optimization method, much in the same fashion as in pyGSTi, but we leave this to be pursued in future work.
The most important diagnostic information for a quantum computing device is often already contained in a lowrank approximation of the processes, states and measurements. Coherent errors are typically the ones that can be corrected by experimental control and are of interest for refining calibration models. The strength of incoherent noise on the other hand is arguably well-captured by average error measures as provided by RB outputs. Moreover, current fault tolerance thresholds often rely on worst-case error measures for which no good direct estimation technique exists [19][20][21][22] and coherent errors in particular hinder their indirect inference from average error measures [23][24][25]. For standard state and process tomography, it was realized that low-rank assumptions can crucially reduce the sample complexity, the required number of measurements and the post-processing complexity [9,[26][27][28][29][30][31][32][33][34][35] as well as improve the stability against imperfections in the measurements [36] using compressed sensing techniques [37][38][39].
In this work, we take a fresh look at the data processing problem of GST from a compressed sensing perspective and regard it as a highly-structured tensor completion problem. We develop a reconstruction method, called mGST, that exploits the geometric structure of CPT maps with low Kraus ranks. In numerical simulations we demonstrate that our structure-exploiting mGST approach (i) allows for maximal flexibility in the design of gate sequences, so that standard GST gate sequences and random sequences work equally well, and (ii) obtains low-rank approximations of the implemented gate set from a significantly reduced number of sequences and samples. This allows us to successfully perform GST with gate sets and sequences that are not amenable to the standard GST implementation pyGSTi [17,40]. As one example, while the sequence design of pyGSTi uses at least 907 specific sequences to reconstruct a two-qubit gate set, we numerically demonstrate low-rank reconstruction from 200 random sequences of maximal length = 7 with runtimes of less than an hour on a standard desktop computer. Thus, compressive GST significantly lowers the experimental resource requirements for maybe the most prominent use-cases of GST making it a tool that can be more easily and routinely applied. At the same time, for the default gate sets and sequences from the standard GST implementation, the novel algorithm matches state of the art results. The runtime and storage requirements of mGST still scales exponentially in the number of qubits as does the amount parameters of the gate set it identifies. This limits the feasibility of the classical post-processing of compressive GST to gate sets acting on only a few qubits without further assumptions. Nonetheless, we demonstrate that coherent errors and depolarizing noise of a 3-qubit gate set can be completely characterized, from as little as 128 sequences of length = 7 on desktop hardware in a few hours.
Our mGST reconstruction method relies on manifold optimization over complex Stiefel manifolds [41][42][43][44][45][46][47] in order to include the low-rank CPT constraints. Such constraints have been used in machine learning [46] and are also starting to become more popular in the quantum information literature, see e.g. [48][49][50]. In order to deal with the non-convex optimization landscape we adopt a second order saddle-free Newton method [51] to this setting. This involves the derivation of an analytic expression for geodesics, as well as an expression for the Riemannian Hessian in the respective product manifolds. Another important motivation for phrasing GST as a randomly subsampled tensor completion problem is to bring it closer to potential analytical recovery guarantees common for related tensor completion problems [52][53][54][55][56][57][58][59], opening up a new research direction.
Finally, being able to perform GST from random sequences enables one to use the same type of data for different increasingly refined characterization tasks from filtered RB [6], cross-entropy benchmarking (XEB) [60] and RB tomography [7][8][9] to GST. Unifying these approaches, random gate sequences can be regarded as the 'classical shadow' of a gate set from which many properties can be estimated efficiently [61]. Compressive gateset tomography provides more detailed diagnostic information and only requires to further increase the amount of data without changing the experimental instructions.
With randomized linear GST [62] and fast Bayesian tomography [63] related alternatives to tackle the GST data processing problem have been proposed. Here, the gates are assumed to be well-approximated by an a priori known unitary followed by a noise channel that is either linearized around the identity [62] or around a prior noise estimate [63]. This allows for a treatment of the outcome probabilities as approximately linear functions. The re-sulting scheme already works for random sequence data but comes at the expense of much stronger assumptions compared to the compressed sensing approach taken with mGST.
The rest of the paper is structured into three parts. In the following section, we formalize the data processing problem of GST as a constrained reconstruction problem. In Section III, we formulate the data processing problem as a geometric optimization task and and derive the mGST algorithm. In Section IV we demonstrate the performance of the novel algorithm in numerical simulations and compare our results with the standard GST processing pipeline of pyGSTi.

II. THE DATA PROCESSING PROBLEM OF GATE SET TOMOGRAPHY
In GST a quantum computing device is modeled as follows. The device is initialized with a state ρ ∈ S := {σ ∈ H : σ 0, Tr[σ] = 1} on a finite dimensional Hilbert space H = C d . Subsequently, a sequence of noisy operations from a fixed gate set (G i ) i∈ [n] can be applied, where we use the notation [n] := {1, 2, 3, . . . , n}. The noisy operations G i : L(H) → L(H) are CPT maps on L (H), the set of linear operators on H. We define [n] * := k=0 [n] k such that i ∈ [n] * defines a gate sequences with length of at most and associated CPT map G i := G i • · · · • G i1 , the concatenation of the gates in the sequence i. In the end, a measurement is performed described by a positive operator valued measure (POVM) with elements (E j ) j∈[n E ] , satisfying j E j = 1 and 0 E j 1 for all j ∈ [n E ]. The full description of the noisy quantum computing device is, thus, given by the triple of a quantum state, a physical gate set and a POVM. Let I ⊆ [n] * be the set of accessible gate sequences with n seq := |I| denoting the number of sequences. The probability of measuring outcome j upon applying a gate sequence i ∈ I is By (p i (X )) j := p j|i (X ) we denote the corresponding vector and, moreover, often omit the argument X . While being a fairly general description, this gate-set model relies on a couple of assumptions: (i) the physical system needs to be well-characterized by a Hilbert space of fixed dimension, (ii) the system parameters need to be time independent over different experiments, and (iii) a gate's action is independent of the gates applied before and after (Markovianity). There exist multiple descriptions of quantum computing devices within the gate-set model that yield the same measurement probabilities on all sequences. Below, we provide a more detailed description of this freedom in terms of gauge transformations. These are linear transformations that, when simultaneously applied to all gates, input state and POVM elements, leave the measurement statistics (2) invariant.
The task of GST is to infer the device's full description (1) from measured data. To this end, one estimates the output probabilities for a set of different sequences I ⊂ [n] * l by repeatedly performing the measurements of the corresponding sequences. Thus, we can state the dataprocessing problem of GST as follows.
Problem (GST data-processing). Let X be a gate set and I ⊂ [n] * l a set of sequences. Given empirical esti- Note that GST aims at solving an identification problem. That is, for sufficiently much data, find the unique device description of the device compatible with the data. In particular, the input data {p j|i } i∈I,j is required to uniquely single out the device description. This is related but distinct from the corresponding learning task to find a description that generalizes on unseen data.

A. Compressive gate set description
At the heart of our approach is to capture this dataprocessing problem as a highly structured tensor completion problem. The structure allows us to reduce the required size and structural assumption of the set I, in order to determine X . It is instructive to visualize the problem with tensor network diagrams. The gate set can be viewed as a tensor of five indices and the action of gate i on the initial state ρ can be visualized as where each leg represents an open index and the joining of legs represents summation of the corresponding indices; see e.g. [64,Chapter 5.1] for more information on the tensor network notation. Neglecting the finite statistics in estimating the probabilities, the GST dataprocessing problem can, thus, be rephrased as the problem of completing the translation-invariant matrix product state (MPS) [65][66][67] or tensor-train [68] from access to a couple of its entries. By the following assumptions one can introduce more structure. First, we assume the elements of the device description X to satisfy the physicality assumptions regarding normalization and positivity. Second, the assumption that they have lowrank approximations yields additional compressibility of X .
In more detail, physically implementable gate sets are completely positive and trace preserving if and only if they admit a Kraus decomposition [69], i.e. the i-th gate implementation can be written as for each i ∈ [n], where r K is the (maximum) Kraus rank of the CPT maps {G i }. We use the notation that K denotes the tensor containing all Kraus operators of all gates and K i contains the Kraus operators for gate i. In terms of tensor network diagrams the decomposition is represented as Moreover, the trace preservation constraints require K i viewed as a matrix in C r K d×d to be an isometry, i.e.
Constraints on a low Kraus rank r K can be naturally enforced in this parametrization by reducing the row dimension of K i . The initial state and POVM elements are constrained to be positive matrices, which we hence parameterize as with A j ∈ C r E ×d and B ∈ C d×rρ , where r E and r ρ are the matrix ranks. For the matrices A j to form a valid POVM, they have to satisfy a similar condition to the Kraus operators, while the initial state is of unit trace if With the physicality constraints incorporated, measurement outcome probabilities are given in terms of tensor network diagrams as Thus, we arrive at a compressive device description X c = (A, K, B) that considerably reduces the amount of parameters compared to the triple X = ((E j ) j∈[n E ] , (G i ) i∈[n] , ρ) when choosing small dimensions r ρ , r K , and r E . Correspondingly, we can adapt the the GST data processing task to demand only a compressive device description.
Problem (Compressive GST data processing). Let X be a gate set and I ⊂ [n] * l a set of sequences. Given empirical estimates {y j|i } i∈I,j∈[n E ] of {p j|i (X )} i∈I,j∈n E and ranks r ρ , r K , and r E , find the compressive device description X c = (A, K, B) of dimension r ρ , r K and r E respectively, so that the normalization constraints (7), (10), and (11) are satisfied.
As before, the set of sequences needs to be large enough so that this identification problem is well-defined. Again a desired compressive device description X c can only be determined up to gauge freedom. Note that for the identification problem to be well-defined, it is not required that the true gate set X that generated the data is of low-ranks itself. As one usually aims to implement unit rank states, unitary gates, and basis measurements, i.e., for r ρ = r K = r E = 1, it can be expected however that a compressive device description X c is often also a good approximation to the true gate set. Next, we discuss the gauge freedom more explicitly and introduce relevant error measures for gates sets X .

B. Gauge freedom and gate set metrics
So far we have not made explicit what 'finding a device description' actually means. What is well studied in the GST and RB literature [6,14,16,[70][71][72], is that without additional prior assumptions, there is a freedom in representing a device in the gate set model. In particular, this freedom needs to be considered when defining a metric for gate sets [70] w.r.t. which we want to recover the device description.
Gauge freedom refers to the following observation. The observable measurement probabilities p j|i of the form (4) are invariant under the transformation for any invertible superoperator T : L(H) → L(H) , where T † denotes the adjoint of T w.r.t. the Hilbert-Schmidt inner product. This invariance is also the wellknown gauge freedom of MPS [73].
If the gate set is universal and the initial state is pure then the gauge transformation T has to be either a unitary or anti-unitary channel. This statement can be seen as follows.
In our case the G i are constrained to be CPT. Hence, G → T −1 (G)T has to map CPT maps to CPT maps. Similarly, T −1 ρ has to be a density operator and {T † (E j )} j a valid POVM.
A more explicit condition on T can be obtained by considering gauge action on entire sequences. For all sequences i, we have where G i (ρ) is a positive operator if the gates G i are CPT. Now T −1 G i (ρ) has to be positive as well for all sequences i. Thus if the gate set is universal, the map T −1 has to be positive and trace preserving for all states. An analogous statement can be made for T † , by considering that T † G † i (E j ) has to be positive-definite for all POVM elements E j . This implies that T has to be a positive map as well.
It has been shown that any positive invertible map T with a positive inverse can be written either as T (ρ) = P U ρ U † P † or T (ρ) = P U ρ T U † P † for U ∈ GL(d, C) and P 0 [74, Theorem 2]. The condition that T needs to be trace preserving then yields P = 1, as can be seen from the Kraus decomposition. Hence, T is indeed either a unitary or anti-unitary channel.
We note that the map ρ → U ρ T U † is positive but not completely positive. However, it has the property that T −1 G i T is CPT whenever G i is CPT. This can be seen by observing that the Choi matrix of T −1 G i T is given by However, actual gate set implementations are noisy and hence not universal in the sense that they cannot prepare any pure state. Therefore, in practice, the gauge freedom can be larger [18]. For instance, if all gate implementations are given by unital channels then an additional freedom exists: depolarizing noise can be commuted through the circuit. Therefore, it can be distributed arbitrarily among initial state, gates, and measurement. Meaningful distance measures for gate sets should have the same gauge freedom as the GST data. The problem of finding gauge invariant distance has been studied by Lin et al. [70]. For individual gate sequences, any measure that compares only the ideal and observed outcome probabilities is naturally gauge invariant. The authors thus propose to use the total variation error, a natural error measure to compare probability distributions, for individual gate sequences. Let denote the probabilities of measuring the jth output of the POVM with elements E j after applying the sequence i of gates in G i to the state ρ. The total variation error for sequence i between two gate setsX = (Ê j ),Ĝ,ρ and X = ((E j ), G, ρ) is defined as The mean variation error (MVE) is defined as [70] MVE I (X , w.r.t. a set of sequences I, where i ∼ I means that i is drawn uniformly from I. Often, we omit the subscript I in the following. The MVE corresponds to taking the natural worst case error measure over the measurement outcomes (the total variation distance) and averaging it over the available gate sequences. Often I is chosen as the set of all gate sequences up to some length . Then the expectation value (18) contains a sum over exponentially many terms. However, since they are all non-negative, they can be estimated sampling efficiently via Monte Carlo sampling [70]. A closely related error measure is the mean squared error (MSE) which averages the squared deviation over all sequences and POVM elements.

III. GST DATA PROCESSING VIA RIEMANNIAN OPTIMIZATION
In the previous section, we defined the compressive GST data processing problem and introduced metrics for for the quality of reconstruction. We now turn to devising a concrete algorithm for the data processing problem. To this end, we formulate the reconstruction problem as a constraint optimization problem of a loss-function for the data fitting. A natural candidate for the loss-function is the MVE restricted to the set of measured sequences. As a proxy we instead minimize the MSE which depends smoothly on the gate set and is, therefore, more suitable for local optimization. In terms of the compressive device description, the MSE (19) can be written as Correspondingly, the compressive GST data processing problem can be cast as the constraint optimization problem: The constraints restrict the objective variables to embedded matrix manifolds. Therefore, algorithms for the optimization problem can be derived by generalizing standard optimization algorithms for functions on the Euclidean space to the geometric structure of these manifolds.

A. The complex Stiefel manifold
In order to formulate our main reconstruction algorithm we need to understand the matrix manifold that encompasses the physicality constraints mentioned in Section II A. We start by summarizing the elementary properties of these manifolds, to then derive a parametrization of geodesics and the Riemannian Hessian, thereby extending what was previously done for their real counterparts in Ref. [42]. For a comprehensive introduction to optimization on matrix manifolds we refer to the book by Absil, Mahony and Sepulchre [75].
Let (K l ) l∈[r] be the Kraus operators of a fixed gate. By stacking them along their row dimension to a new matrix K ∈ C dr×d , we can write the CPT constraint as K † K = 1. In the following we set D = dr. The set is called the D × d complex Stiefel manifold. This manifold is the set of isometries of the Euclidean space and contains the special cases of the sphere St(D, 1) and the unitary matrices U (D) = St(D, D). We regard it here as a submanifold of C dr×d . The tangent space of St(D, d) at K is given by The canonical inner product of ∆ 1 , ∆ 2 ∈ T K St(D, d) can be defined as with Γ = 1 − 1 2 KK † . Another choice is the standard Hilbert-Schmidt inner product of the embedding matrix space. However the advantage of the canonical inner product is that it weights all degrees of freedom on the tangent space equally. The Stiefel St(D, d) together with the metric given by (24) is a Riemannian manifold. The normal space is defined by The projector onto the normal space at position K is given by for X ∈ C dr×d and we can write the projector onto the tangent space at K as We wish to optimize the MSE over St(D, d). In analogy to the optimization over U(n) in [41], we will move along geodesics, which are the locally length minimizing curves. In Appendix A we show that within St(D, d), a geodesics starting at K t=0 ≡ K and going in the direction ∆ ∈ T K St(D, d) can be written as with Q,R given by the QR decomposition of (1−KK † )∆. Note thatK t | t=0 = ∆. Often simpler curves that just satisfy K 0 = K andK t | t=0 = ∆ are used instead of the geodesic in order to save computation time [75]. However, computing the exponential of the 2d-dimensional matrix in Eq. (27) provides no bottleneck in our scenario as the inversion of the 2nd 2 r K -dimensional Hessian is more costly (see Section IV F). In order to identify the Riemannian gradient and Hessian, we generalize results from the real case [42] to the complex case. Then we use the second order Taylor approximation of the objective function, which will be given below in terms of the MSE (19) along geodesics (see Appendix B). The same treatment can be applied to the POVM given by the matrices A j from the decomposition (9), where we define A as the matrix obtained from stacking the A j along their row dimension. The physicality constraint on A is then equivalent to A ∈ St(dn E , r E ) with n E being the number of POVM elements and r E their maximal rank.

B. The mGST estimation algorithm
With a better understanding of the underlying manifold structure we can now formulate a concrete optimization approach to tackle the estimation problem (21). The least squares cost function (19) is a polynomial of order at most the sequence length squared in the parameters of G, with a highly degenerate global minimum due to the gauge freedom. In analogy to the alternating minimization techniques which are successful for matrix product state completion [52,76,77] we alternate between updates on A, K and B. Each update would naively be done via a local optimization approach such as gradient descent. However, we observe that following the gradient direction on the respective manifolds is problematic around saddle points, which are frequently encountered in our optimization problem. In principle the gradient direction points away from saddle points, yet the norm of the gradient can be arbitrarily small. There are different approaches in the literature to deal with this problem. For instance information about the curvature can be included [51] or, if a saddle point is encountered, random update directions can be chosen to escape the area of vanishing gradient [44,78]. We find that the so-called saddle free Newton (SFN) method [51] yields considerably better results than first order methods. There the update direction is given by − |H| −1 g with H being the Hessian and g the gradient and the absolute value |H| define by spectral calculus. An instructive way to see why this leads to a speedup is to write the Hessian H as Since the vectors |v i form a basis, this can be interpreted as a rescaling of |g by |λ i | −1 in the directions |v i . As with the standard Newton method, this leads to a large rescaling if the curvature in a particular direction is small, resulting in large steps even close to the saddle point. Taking the absolute value of the eigenvalues then ensures that saddle points are repulsive. For numerical stability it is beneficial to introduce a damping term that offsets the eigenvalues of H that are very close to zero before the inversion. Algorithm 1 describes a single step of the damped saddle-free Newton method with damping parameter λ and is a generalization of the original SFN method [51] to manifolds.
Algorithm 1 is formulated in a way that is compatible with an update in Euclidean space as well as an update on the Stiefel manifold. In Euclidean space we update along the curve Y t (·, ·) : with the curve given by the geodesic (27) and T being the tangent bundle on St(D, d). The step size is determined by locally optimizing over the parameter t using standard gradient free optimizers. We derive an expression for the Hessian on the Stiefel manifold in Appendix B. In Appendix C, we also provide a detailed discussion and expressions for the optimization in complex Euclidean space.
Algorithm 2 describes the main mGST routine. It can be run with different choices of smooth objective functions and we use the MSE (19) by default. In our numerics we often find that optimizing the log-likelyhood function after the MSE can improve estimates (see Appendix D for a discussion).
The algorithm alternates updates on A, K and B. The updates are performed using the SFN update of Algorithm 1. For the updates on A and K, the SFN update is computed on the tangent spaces of the respective Stiefel manifolds, while for the unconstrained matrix B, the update is computed in complex Euclidean space.
The unit trace on ρ = BB † is the only physicality constraint that is not naturally enforced through the manifold structure. The trace value can be seen as a normalization constant that would be naturally fixed to 1 at an optimal point, if no noise due to the finite statistics is present in the data. One can perform an additional normalization step after each update B to exactly enforce the constraint. We observe no difference in the performance of the algorithm with or without explicitly enforcing this constraint.
In order to achieve good convergence, we perform the optimization with mGST in two consecutive steps: we start from a random initialization and perform a coarse grained optimization with a small batch size κ, i.e. only using κ many random gate sequences from I for each update step. The batching of data results in lower computation time for the derivatives and adds a factor of ran-domness to the optimization, which avoids getting stuck at suboptimal points to a certain degree. We terminate the first optimization loop when the objective function L I (A i , K i , B i |y) is smaller than an early stopping value δ, which is obtained from the data as follows.
For a number of m samples per sequence the outcome probabilities of each sequence for the true gate set are given by where k j|i is number of times outcome j is measured upon applying the gate sequence i. Due to Born's rule, k j|i is distributed according to the multinomial distribution M(m, (p 1|i , . . . , p n E |i )) with probabilities {p j|i } j and m trials. We estimate the expectation value of the objective function from the values y j|i . This provides us with a rough estimate for how low the objective function value can become, given the sample counts k j|i . Then we set the early stopping value to be twice that estimate, Hence, we require the objective function on the full data set to be close to its expectation value for the measured probabilities y j|i obtained from m samples. While computationally inexpensive the mini-batch stochastic optimization does not converge to an optimal point on the full data set I. In a second optimization loop, we initialize the mGST algorithm with the result from the first run and use all the data for the updates. Formally, we choose the batch size κ = |I| and, thereby, make the random batch selection obsolete. We perform these more costly update steps until the change in objective function reaches a desired relative precision , or a maximal number of iterations is exceeded. The first optimization run is initialized with a random gate set parameterized by A 0 , K 0 and B 0 (see Section II A). For the random initialization we make use of the Gaussian unitary ensemble (GUE). A matrix H belongs to the GUE if H = (M + M † )/2, where M is a complex Gaussian matrix, i.e, real and imaginary part of each M ij are independently drawn from N (0, 1), the normal distribution with zero mean and unit variance. In this case we we write H ∼ GUE. For A 0 and each gate in K 0 we take the first d columns of e iH with H ∼ GUE to obtain a random isometry K 0 . For B 0 we take a complex Gaussian matrix and normalize it such that Tr[B 0 † B 0 ] = 1.
Importantly, due to the nature of non-convex optimization, several initializations can be needed to converge to a satisfactory minimum.

IV. NUMERICAL ANALYSIS
In this section, we evaluate the performance of mGST in different scenarios in numerical simulations. In particular, we compare its performance to the state-of-the-art implementation for gate set tomography, pyGSTi [18], in the regimes where both methods can be applied.
For pyGSTi to be applicable one has to use structured gate sequences inspired by standard quantum process tomography. In Section IV A we evaluate the performance of mGST and pyGSTi on minimal measurement sequences and different models to find that mGST benefits from flexibility in the sequence design and a fully general model parametrization. Section IV B numerically validates the expected inverse square-root scaling of the reconstruction error with the number of measurement samples per sequence for different noise regimes. Section IV C numerically determines the required number of random sequences to accurately reconstruct simple and random gate set models with mGST for different Kraus ranks. In Section IV D we follow up with a numerical demonstration of unitary noise characterization for a three-qubit gate set using a priori knowledge in the initialization. Finally, in Sections IV E and IV F we discuss the choice of initialization and hyper-parameters, as well as the runtime of mGST.
For a model of n gates reconstructed from m measurements of sequence length , we validate the performance of mGST by computing the MVE (18) over all possible n sequences, or 10 4 random sequences of length if n > 10 4 . Usually m min(n , 10 4 ) and the MVE can be thought of as a generalization error on the predicted output probabilities of the gate set estimate. The gate sets studied in this section all use the same target initial state |0 0 | and computational basis measurement, although with different levels of noise applied to them. For instance, we often use global depolarizing noise, which acts on a quantum state ρ as ρ → (1 − p)ρ + p 1/d. For the numerics presented here, we use a maximum of 100 reinitializations (if not stated otherwise). A discussion of the required number of initializations is given in in Section IV E. A Python implementation of mGST and a short tutorial can be found on Github [79].

A. Gate set and measurement structure
We compare mGST and pyGSTi for the minimal number of sequences doable with each method and for gate sets of different conditionings, without using the compression capabilities of mGST yet. We find that mGST is more flexible in the sequence design and model parametrization, while generating estimators with lower mean variation errors in several regimes.
The traditional strategy for GST, akin to standard quantum process tomography, is to generate a frame for L(H), measure each gate in that frame and generate an estimate for each gate by applying the pseudo-inverse of the measurement operator.
This is particularly important for the first reconstruction step in pyGSTi where the sequences that generate the frame are called fiducials. The strategy of pyGSTi is to obtain an initial estimate via the pseudo-inverse, followed up by local optimization of a particular cost function [18]. In contrast, we perform mGST using random initializations and, thereby, not rely on designated fiducial sequences. Figure 1 compares mGST to pyGSTi focussing on the regime of very few gate sequences, showing what is needed in terms measurement effort to obtain low mean variation errors for different gate sets. In order to test pyGSTi in the regime of low sequence counts, we replace the 5 standard fiducial sequences with 2 or 3 fiducial sequences drawn uniformly at random, thereby reducing the total sequence number from 92 to 18 or 39 sequences. Since mGST is compatible with any sequences design, we use between 10 and 58 random sequences for mGST to explore the low sequence count region.
The first gate set we study is the so-called XYI model, the standard single qubit example in the pyGSTi package [40]. The XYI model consists of the identity gate, a π/2 X-rotation and a π/2 Y-rotation on the Bloch sphere, with initial state |0 0 | and measurement in the com-putational basis. Results for the ideal XYI-model can be seen on the left in Figure 1, with mGST and pyGSTi performing identically for 92 sequences. Comparing the results for 18 sequences we find that mGST does not converge on more than 50% of trials, which reflects in the median MVE being above 10 −1 , while pygsti achieves lower median errors. Comparing the 38 and 39 sequence medians however, we find that mGST yields lower error models than pyGSTi.
The subsequent models analyzed successively deviate from the simple unitary XYI-model and highlight the versatily of our manifold approach. Since the full CPT parametrization used in our optimization (21) is agnostic to any special gate set properties we expect it to perform well for all possible CPT maps as gate implementations. For instance, for random and specific non-Markovian channels. pyGSTi on the other hand uses a parametrization of Lindblad type and is therefore based on a more limited model space.
To illustrate this comparison, we perturb the XYImodel by adding amplitude damping noise to each gate. The amplitude damping channel can be written in terms of the Kraus operators K 1 = 1 0 0 Γ and , which arise e.g. from the Jaynes-Cummings model of a qubit system interacting with a quantized bosonic field [80]. How well mGST and pyGSTi perform on a model with Γ = 0.94 can be seen in the center left block of Figure 1 (XYI + Weak damping). We find generally similar performance, with mGST being a bit more accurate on 38 sequences, and a bit less accurate on 92 sequences with 10 6 samples. Increasing the interaction time between qubit and environment leads to memory effects and strong non-Markovianity of the amplitude damping channel at Γ = −0.6. This scenario is shown in the center right plot of Figure 1 and we see that while the accuracy of mGST is the same as before, the model parametrization of pyGSTi cannot fit the model with MVEs below 10 −2 , independent of the sequence or sample count.
For the last comparision (rightmost block in Figure 1) we look at the performance for random full Kraus rank channels. Each channels is constructed by drawing a Haar random d 3 × d 3 unitary and then taking its first d columns. The resulting d 3 × d matrix is an isometry and therefore constitutes a valid set of Kraus operators. Note that this construction is different from the previous construction of random channels via the Gaussian unitary ensemble. The results show that mGST can reconstruct these models from low sequence counts, while pyGSTi does not yield good estimators. Using the standard sequence design of 92 sequences, mGST and pyGSTi have identical accuracy again, suggesting that random channels are typically well within the model space of pyGSTi after all. These demonstrations show that mGST is indeed flexible in the sequence design with state of the art performance for arbitrary gate set implementations.

B. Number of samples per sequence
The probability associated to every sequence is estimated from a finite number of samples. Here, we study the resulting effect on the reconstruction accuracy as measured by the MVE more closely.
For a high number m of samples per sequence, each probability y j i in the objective function is estimated with an error of order 1/ √ m. Therefore, we expect the MVE to also decrease as 1/ √ m if the algorithm converges to the global minimum. This scaling was observed to hold true for pyGSTi [18]. In order to be able to compare the scaling of mGST directly to the one of pyGSTi, we use a standard pyGSTi setting: The gate set is the XYI-model (with π/2-rotations) and the gate sequences are the standard pyGSTi sequences for this model with a maximum sequence length of = 7.
We add noise to the gate set by varying the amount of depolarizing noise with strength p on each gate and also overrotating each gate by a random unitary. The random unitaries are given by e iγH with H ∼ GUE. In particular, this means that H can be bounded on average State preparation and measurement are assumed to be noise-free in this setup, however for a fixed sequence length the depolarizing noise per gate is equivalent to a global depolarizing channel applied before measurement, since it commutes with the unitary gates. Figure 2 depicts the resulting MVE-scaling of the reconstruction where data was generated using different numbers of samples per sequence m. We observe that mGST follows the expected scaling in m, matching the scaling of pyGSTifor different levels of unitary and depolarizing noise.
It is worth noting that in the plot on the left hand side of Figure 2 the blue lines that show the results for very low depolarizing noise indicate a lower error reached for pyGSTi. We observe this behaviour consistently for very low noise. This discrepancy originates from the choice of objective function. In the last step of its optimization procedure, pyGSTi uses a maximum likelihood estimator whereas mGST effectively minimizes a least squares objective function given by the MSE. In Appendix D we discuss this behaviour in more detail and show that while in ideally conditioned settings less accurate in MVE, mGST achieves lower values for the least squares objective function than pyGSTi. We further show that if the least squares estimate of mGST is used as a starting point for maximum likelihood estimation, the final MVE is reduced and precisely matches the results of pyGSTi.

C. Number of sequences
The arguably most challenging experimental requirement of GST is the number of measurement settings (sequences) that are required for a successful gate reconstruction. One of the main motivations of compressive GST is to employ structure constraints, i.e. to reduce the number of degrees of freedom of the reconstruction problem, in order to reduce the required number of measurements. Instead of reconstructing arbitrary quantum channels we aim at reconstructing low-rank approximations of the gate set elements. In addition, we expect that by using the mGST algorithm, compressive recovery is possible from already a 'few' randomly selected sequences. We here numerically demonstrate that this is indeed the case.
The top row of Figure 3 shows the median performance in MVE against the number of randomly chosen sequences for different Kraus ranks. On the left are the results for the single qubit XYI model as defined in Section IV A. On the right are the results for the XYICNOT gate set that is based on the identity, CNOT and Pauli-X and -Y rotations on each qubit individually, with rotation angle π/2.
We observe a phase transition in the MVE that indicates a minimal number of sequences that are required for the successful reconstructions of the gate sets. As expected, constraining the reconstruction to a lower Kraus rank indeed reduces the amount of required sequences in the reconstruction in most cases.
An intriguing exception is the r K = 1 reconstruction of the two qubit gate set that exhibits the worst reconstruction performance compared to higher rank constraints. We suspect that this is due to the optimization problem being more dependent on the initialization for r K = 1. In more general settings, it has been observed that the optimization over matrix-product states with fixed Kraus rank can be unstable and using rank-adaptive optimization techniques yield much better performance [82,83]. This motivates to use a slightly higher rank in the optimization than the expected rank of an effective approximation of the gate set. In accordance with this intuition, we find that it is beneficial to constrain the optimization to r K = 2 in order to achieve an accurate unit-rank approximation. The same effect is also observed in the single qubit example when taking a detailed look at the number of required initializations (see section IV E) , yet less pronounced. In the bottom row of Figure 3 we show the recovery rates for random unitary models, with the reconstruction now using a fixed Kraus rank of r K = 2.
Note that there are three sources of randomness present in the data, first the Haar-random unitary gates, then the random drawing of gate sequences and finally the random initialization of the algorithm. Each shade of green corresponds to one random gate set and the recovery rate tells us how many of the 10 random sequence sets lead to a successful reconstruction, given a budget of 33 initializations. We find that the random single qubit gate sets all have similar recovery rates, with a successful reconstruction possible from n seq = 20 to n seq = 30 sequences, and a high rate of recovery at n seq = 100 sequences. In the two qubit case a different picture emerges, where two of the random gate sets show a high recovery rate at n seq = 200 sequences (akin to the XYICNOT-model), while the least favorable random gate set was only recoverable at n seq = 500 sequences. This shows that random gate sets for two qubits can have very different conditioning.

Sequence number comparison to pyGSTi
Comparing the number of random sequences needed for mGST and the number of sequences for pyGSTi is not straightforward. The standard pyGSTi data-processing pipelines crucially relies on specific, fixed sequence construction. For this reason pyGSTi cannot be applied to the type of data that we use here. We can however compare the number of random sequences with the number of deterministic sequences that the standard implementation of pyGSTi uses. For the single qubit XYI-model, the minimal number of sequences given in the pyGSTi implementation is n seq = 92. This is significantly larger than the number of random sequences at which the phase transition of mGST in Figure 3 appears. However, the n seq = 92 sequences are overcomplete by design and we find that pyGSTi can also reconstruct the XYI model with n seq = 48 sequences. Yet we find the same sequence design not to be successful for three Haar-random single qubit gates, indicating that the choice of sequences is well-tailored to the XYI-model. The reduction in sequences becomes more pronounced for the two-qubit gate set studied in the top right of Figure 3. For this gate set, the minimal number of gate sequences that pyGSTi uses is n seq = 907, which is significantly larger than what mGST needs.

D. Characterizing unitary errors using prior knowledge
In the previous section we demonstrated compressive gate set tomography for one-and two-qubit gate sets using agnostic random initializations. A major obstacle in going beyond reconstructing entire two-qubits gate sets even compressively on desktop hardware is that besides run-time and storage also the number of required random initializations till proper convergence grows with the number of qubits. This situation can be remedied by using prior knowledge for the initialization. For example, if we have partial trust in the experimental realization, we can use the unitary target gate set for the initialization. To showcase the characterization of unitary errors using this prior knowledge, we take a three-qubit gate set that is the direct generalization of the previous two-qubit XYICNOT model, by adding the local Xand Y-rotations as individual gates to the 3rd qubit and adding a CNOT between qubits 2 and 3. We then apply a global random rotation to each gate individually, as well as depolarizing noise on each gate. From random sequences of fixed sequence length we can then, in theory, fit the noisy model perfectly via an r K = 1 approximation, as the depolarizing channels commute with the unitary gates and can be pulled into initial state or measurement. In Figure 4 we see that mGST is indeed able to precisely reconstruct the rotated gates, as shown by the average diamond norm error. We chose a comparatively high number of 10 5 samples per sequence to showcase that high accuracy can indeed be realized using this method: for instance, only 256 sequences are enough to achieve an average diamond norm distance of around 0.007 between the reconstructed unitary gates and the true unitary gates, which include overrotations. The fact that these overrotations were modeled as being global on all 3 qubits suggests that we can efficiently characterize unitary crosstalk as well.

E. Implementation details and calibration
We now provide more details on the simulations, the criteria for successful recovery and the required number of initializations. To simulate measurements on a gate sequence i, we first compute the outcome probabilities p j|i from Eq. (16) of the POVM elements according to the model gate set in question. Afterwards we draw m samples from the multinomial distribution M (m, (p 1|i , . . . , p n E |i )), where j p j|i = 1. Let k j|i be  (20). Therefore the computation time of the objective function and its derivatives scales linearly in κ, making a small batch size favorable. However it cannot be set too small, otherwise the update directions become highly erratic and no convergence is reached. A general rule of thumb is to set the batch size close to the number of free parameters in the model. Another hyper-parameter is the damping value λ for the saddle-free Newton method described in Algorithm 1. We find that a fixed value of λ = 10 −3 leads to the best results across the models tested.
Judging whether mGST recovers a gate set by looking at the attained objective function value can only be done if the set of measured sequences is informationally complete. Then there is a unique (up to gauge) global minimum in the least squares minimization problem and the minimum corresponds to the true gate set in the limit of infinitely many samples per sequence.
In Figure 5 we take a look at the correlation between the final least square objective function value L(X , y) and the mean variation error MVE(X , X ). We see that for a low number of sequences (10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20), a low objective function value does not imply a low MVE, yet for higher numbers of sequences, an objective function value below 10 −3 implies an MVE around 10 −2 . For sufficiently many sequences, the grey line indicating our success criterion clearly separates two clusters of points, meaning that no intermediate quality fits are found in our model space. In this sequence regime either the algorithm converges to a fit as good as the sample count allows, or it does not converge at all. Therefore, restarting the algorithm when an initialization turns out to be bad yields practically optimal results. A thorough analysis of the probability of obtaining an informationally complete set of random sequences is left for future work.
To give an intuition on how many initializations are required for mGST to converge, we can take a look at data from a modified XYI-model with gates {1, e i α 2 σy , e i α 2 σx }, simulating more difficult gate set conditioning. Figure 6 shows histograms for the number of reinitializations needed for convergence. The data combines the results for α between π/18 and π/2, with depolarizing noise of strength p = 0.001 on each gate and p = 0.01 on the initial state, as well as a maximum of 100 reinitializations.
We observe that in 48.4% of all cases convergence is reached on the first attempt, and in 90.4% of cases 4 or fewer reinitializations are required. The histogram indicates that the chance of needing multiple reinitializations rapidly decays and that rank-1 optimization is more sensitive to bad initializations compared to rank-4 optimization.

F. Runtime and scaling
In order to assess the runtime scaling of Algorithm 2 in the problem parameters we identify the two most time consuming steps as the computation of the second derivative terms in H and the diagonalization of H for the SFN update on the gates.
Recall that r K denotes the Kraus rank of the gate estimators, n the number of gates in the gate set, the number of gates per sequence and κ the number of sequences per batch. With ω being the exponent of matrix multiplication, the computation of the second derivative terms scales as O(κ 3 d 2ω ), while the eigendecomposition of H scales as O (2nr K d 2 ) 3 . A computationally less expensive variant is to not optimize over all variables of the full gate tensor at once, but over the individual gates one after another. The complexity of the eigendecompositions then reduces to O n(2r K d 2 ) 3 , which is beneficial for large gate sets. However, this approach also leads to slower convergence and we choose to optimize over the full tensor by default. Table I contains runtimes for the system sizes studied in our numerical experiments. We find that low-rank reconstruction of single qubit gate sets can be achieved within seconds and low rank 3-qubit gate sets can be done within a few hours on standard desktop hardware.

1Q
8s 1m 2Q 3m26s 31m 9h30m 3Q 5h30m Table I. Average runtimes for 1, 2 and 3 qubits with selected Kraus ranks for a single initialization on a modern 6-core desktop CPU. The one qubit scenario was taken to be the same as the top left of Figure 3 for α = π/2 and the two qubit scenario the same as in Figure 3. For the three-qubit runtime we tested a system of 4 random 3-qubit unitaries on 1000 sequences with maximum sequence length ≤ 8.

V. CONCLUSION AND OUTLOOK
We have revisited the data processing task of GST from a compressed sensing perspective regarding it as a highly structured and constrained tensor completion problem. In this formulation, we can naturally require the reconstructed gate set to be physical and, moreover, of low rank. Compressive gate set tomography, thus, aims at extracting considerably fewer parameters of the gate set.
At the same time we ave argued that the low-rank approximation to the implementation of a gate set contains the most valuable information about experimental imperfections for the practitioner.
The set of Kraus-operators of a low-rank gate can be regarded as isometries that make up the complex Stiefel manifold. This observation has motivated the solution of the compressive GST data processing problem via geometrical optimization on the respective product manifolds. We have devised the optimization algorithm mGST that performs an adapted saddle-free Newton method on the manifold. To this end, we have derived the Riemannian Newton equation, Hessian equation and geodesic curves.
In numerical experiments we have studied the performance of the mGST algorithm. We have compared it to pyGSTi, the state-of-the-art approach to the GST data processing problem, in settings where both algorithms can be applied. We have found that in these settings mGST matches the performance of pyGSTi, while offering a larger model space, more flexibility in the sequence design and allowing for low rank assumptions. Moreover, we have demonstrated numerically that making use of the low-rank constraints significantly reduces the required number of measured sequences and the runtime of the reconstruction algorithm for a standard single and two qubit model. Importantly, we have found that we can successfully reconstruct generic unitary channels and depolarizing noise of one-and two-qubit gate sets from random gate sequences. This reduces the demands of GST both for experiments and classical postprocessing: the data that compressive GST requires is virtually identical with the experimental data produced by randomized benchmarking experiments. The classical post-processing for a low-rank reconstruction of two qubit gate-sets takes only minutes even on desktop hardware. We expect that this can render compressive GST a routinely applied characterization protocol in the experimental design and calibration of gate sets.
Furthermore, compressive GST is moderately scalable to a few more qubits -especially if making use of prior information for the initialization. We demonstrated the tomography of unitary errors for three-qubit gate sets, using only a small number of random gate sequences, and with the post-processing still running on desktop hardware. We expect that even going slightly beyond three qubits is feasible by simply using more computing power. We leave it to future work to further tweak the numerical implementation in order to improve the scalability of the classical post-processing.
Finally, besides making GST more applicable and flexible in practice, our reformulation is motivated by bringing it closer to theoretical recovery guarantees quantifying a required and sufficient number of random sequences for accurate reconstruction. Regarding the data processing of GST as a translation-invariant matrix-productstate/tensor-train completion problem makes it more amenable to prove techniques from compressed sensing.
For example, establishing local convergence guarantees for mGST would allow one to quantify the assumptions on the experimental implementation that justify certain initialization of the algorithm. We hope that our work can serve as a foundation and inspiration in the quest of establishing mathematically rigorous guarantees for GST.

ACKNOWLEDGMENTS
We thank Lennart Bittel for helpful discussions and hints regarding optimization methods and their implementation and Markus Heinrich for a discussion on in-verses of positive maps on operator spaces. We are also particularly thankful to Kenneth Rudinger and Robin Blume-Kohout for valuable feedback on a previous version of the manuscript and help with finding an error in our code. The work of RB and MK has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within the Emmy Noether program (grant number 441423094) and by the German Federal Ministry of Education and Research (BMBF) within the funding program "quantum technologies -from basic research to market" via the joint project MIQRO (grant number 13N15522). IR acknowledges funding from the BMBF (DAQC) and the Einstein Foundation (Einstein Research Unit).

VI. APPENDIX
In this appendix, we provide the mathematical details required for the saddle-free Newton method within the Riemannian optimization framework, see Appendices A, B, and C. Moreover, we compare the dependence of the mGST MVE on the choice of of objective function (mean squared error vs. maximum likelihood) in Appendix D.
A. Geodesics on the Stiefel manifold Edelman, Arias, and Smith [42] derived the geodesic on the real Stiefel manifold by solving the respective geodesic equation. We now show that the simple generalization given in Eq. (27) is indeed the correct geodesic in the complex case. For a curve K t ≡ K(t) the general geodesic equation is [75,Chapter 5.4

, Proposition 5.3.2]
where the Christoffel symbol C Kt depends on the chosen metric. Here, we use the canonical metric with Γ = 1 − 1 2 K t K † t . Using the Einstein summation convention, the Christoffel symbol at K can be computed as where C k Kt is the k-th component of the Christoffel symbol at K t with respect to a basis {E k , E * k } k∈[Dd] on the ambient space C D×d .
We determine the derivatives of g needed for the Christoffel symbol by explicitly writing out g as from where we can find the derivatives by K and K * as With these derivatives we can calculate first and second derivatives of K t : We immediately see that K † tK1 = A 2 − R † R, which is Hermitian, as A is skew Hermitian. We will now show that K † tK2 = 0 starting with To simplify the last expression, we set and obtain K † tK2 = (U † 00 K † QU 11 + U † 10 Q † KU 01 )RA. From the series representation of the matrix exponential we gather that U 11 = 1 + R · X for some matrix X. Moreover U 10 = RX and U † 10 =X † R † for someX, leading to This shows that theK t lies in the normal space, leaving us with the terms in the geodesic equation (35) that depend only onK t : The last line follows fromK 2 A =K 2 and our previous observation that K † tK2 = 0, which implies that K † tK2 = 0 as well. The remaining termK † 2K 2 can be computed similarly to K † tK2 and we obtain where we used again that K † QR = R † Q † K = 0 in the last line.
We can now put all the terms obtained by multiplying Eq. (49) with K † t from the left together and find We see that these remaining terms are Hermitian and therefore the left hand side of Eq. (49) is in the normal space.

B. Complex Newton equation
In this section, we derive the Riemannian Hessian operator and solve the Hessian equation to obtain an update direction on the tangent space, which we can follow along the geodesic defined in Eq. (27). This can be done for each gate individually, or simultaneously over all gates, in which case we operate on the Cartesian product St(D, d) ×n of single Stiefel manifolds. We consider the latter case, whereby we obtain the single Stiefel Newton equation (Eq. (77)) as a byproduct. The method is based on the real case [42]. See also [43] for a recent treatment of second order optimization on the complex Stiefel manifold, where instead of following geodesics, each optimization step is done in Euclidean space followed by a projection onto the manifold.
First let us make a general observation that will be useful at several points.
is given by Proof. It is straightforward to see that X ∈ T K St(D, d) by applying the projector onto the tangent space: To show that X solves Eq. (3), we will use that K † ∆ is skew Hermitian, as well as the fact that Re Tr[HS] = 0 for any skew Hermitian matrix S and Hermitian matrix H. Plugging X = F * − KF T K into Eq. (3) we obtain Our goal is to simultaneously update all gates along the geodesic n i=1 K i (t) ∈ St(D, d) ×n , with the single Stiefel geodesics K i (t) being given by Eq. (27). We define the initial directions as ∆ i =K i (0). The first step to identify the Riemannian gradient and Hessian is to compute the second order Taylor series expansion of L in t at t = 0. Using (K i ) lm ≡ K ilm and Einstein notation we find L(K 1 ⊕ · · · ⊕ K n ; K * 1 ⊕ · · · ⊕ K * n ) = L| t=0 + 2 Re The solution for G in the canonical metric (24) is given by as per Lemma 3.
where K(t∆) satisfiesK(t∆)| t=0 = ∆ (see [42] for a discussion of the real case). The individual terms in Eq. (74) can be determined from our general Taylor approximation in Eq. (70), i.e. with i = j = 1, if we take K = K 1 . The term ∂L ∂K ∂ 2 K(t∆) ∂t 2 t=0 contains second derivatives of the geodesic given in Lemma 2, which we write out next.K(t) is given byK It follows using QR = (1 − KK † )∆ =: Π∆ and A = K † ∆ from the definition of the geodesic, thaẗ Putting the terms together, we arrive at The terms involving L KK and L K * K satisfy L KK [∆, Ω] = L KK [Ω, ∆] and L K * K [∆ * , Ω] = L K * K [Ω * , ∆], by the symmetry of second derivatives. Using this symmetry property we obtain the full Hessian (74), which turns out to be where it is helpful to note that Eq. (76) is related to Eq. (75) via a symmetrization of the L K term.
Theorem 5. Let vec(∆) be the row major vectorization of ∆ ∈ T K St(D, d). Furthermore, let T andL KK be defined by T vec(X) = vec(X T ) and L KK (∆, · ) =L T KK vec(∆). Then the solution ∆ of the linear equation in vec(∆) and vec(∆ * ) given by is the update direction along the geodesic given in Lemma 2 for the complex Newton method of a real function L at position K ∈ St(D, d).
Proof. The update direction ∆ for the standard Newton method [42] is determined through the equation which can be solved by rewriting the left hand side as Hess(∆, Ω) = f (∆), Ω K (for some yet to be determined f ) and setting Ω = P T (X) with arbitrary matrix X. This leads us to P T (f (∆)), in the second line the scalar product is extended from the canonical scalar product initially defined on T K St(D, d) to C D×d and we will use the same notation for both. The second line follows from the fact that any matrix X can be decomposed as X = P T (X) + P N (X) and from The same argument can be made for the L K * K term, leading to To identify f K (∆) we rewrite the second line in (76) as follows: by using Π = Π † , ΠΓ −1 = Γ −1 Π = Π, K † Γ −1 = 2K † . Eq. (81) implies P T (f (∆)) = −G and it remains to compute P T (f K (∆)), as P T (f KK (∆)) = f KK (∆) and P T (f K * K (∆ * )) = f K * K (∆ * ) (see Lemma 3). After a straightforward computation using Γ −1 = 1 + KK † , ΠΓ −1 = Π, P T (ΠZ) = P T (Z), as well as K † Π = ΠK = 0, we get with skew(A) = (A − A † )/2. Finally, by plugging Eqs. (82), (83) and (84) into Eq. (81), we obtain the Newton equation This is a linear equation in ∆ and ∆ * that can be solved via rewriting it as an equation in vec(∆). Using row-major vectorization with vec(AXB) = (A ⊗ B T ) vec(X), the matrices T andL KK defined by T vec(X) = vec(X T ) and L KK (∆, · ) =L T KK vec(∆), we arrive at the final equation for the single gate case Now for the simultaneous optimization over all gates on St(D, d) ×n , the Hessian as defined in Eq. (74) is determined by including all terms in Eq. (70). The Newton equation reads Hess(∆ 1 ⊕ · · · ⊕ ∆ n ; Ω 1 ⊕ · · · ⊕ Ω n ) = − We can use Lemma 3 again and obtain We are now faced with an equation of the type Ax + Bx * = c, a solution to which can be obtained by solving Eq. (88) is an equation on the tangent space and can be solved by finding a basis therein. However, in order to avoid a basis change at every step, we choose to solve it on the ambient space by setting ∆ i = P Ti (∆ i ) and ∆ * i = P * Ti (∆ * i ). The matrix equation for the update directions ∆ i is now given by where the sub matrices H G←∆ and H G←∆ * can be identified from Eq. (88). The update directions for the saddle-free Newton method 1 are calculated by applying to the right hand side of Eq. (90). In practice we use 1 2 (H + H † ) as the Hessian, since there exist more efficient methods for diagonalizing a Hermitian matrix compared to an arbitrary matrix.
where the higher order contribution r (dZ, dZ * ) satisfies meaning that derivatives of the objective function reduce to the derivatives of the tensor p. Taking the derivative of a tensor network w.r.t. one of its constituent tensors can be easily done in the pictorial representation by removing the respective tensor. For instance, ∂p i ∂K * can be calculated as shown in Figure 7, using the product rule. Care has to be taken for the order of open indices when removing a tensor. In practice we do not calculate the full tensor ∂p ∂K * of size n and only compute ∂p i ∂K * for i ∈ I, since usually |I| n . In Section IV B of the main text we saw that the MVE of the pyGSTi estimate is lower than the mean variation error of the mGST estimate. Here, we show that this performance difference is due to the choice of objective function.
In general, there are several potential reasons the observed discrepancy in performance: (i) mGST converges to a local minimum, (ii) mGST is near the global minimum but not fully converged, (iii) the final objective function used in the pyGSTi procedure is better suited to produce estimates with a low generalization error. Option (ii) can be easily ruled out by increasing the maximum number of steps and we find that due to the 2nd order optimization method used in mGST, convergence in objective function up to relative precisions of 10 −4 − 10 −6 is attained well before the step limit of 200 used to obtain the data. Option (i) seems unlikely, since we observe the correct scaling in the number of samples per sequence and the results of pyGSTi and mGST are fairly close. In Figure 8 we see that if we compare the least squares objective function values reached by pyGSTi and mGST, then mGST fits the data slightly better in this metric. This is to be expected, since pyGSTi performs maximum likelyhood estimation (MLE) which typically leads to small total variation errors. More specifically, in the final step of pyGSTi, the log-likelihood function logL I (A, K, B; y) := m i∈I j y j|i log p j|i (A, K, B) is maximized. In Figure 8 we show that if we take the estimated gate set given by mGST and use it as an initialization for maximum likelihood estimation, then the final MVE is of the same accuracy as results of pyGSTi. To perform this maximum likelihood estimation we use Algorithm 2 on the negative log-likelihood function defined in Eq. (95). This shows that option (iii) is the reason for the discrepancy between the accuracy of pyGSTi and mGST in the low noise and low sample size regime.