Fitting quantum noise models to tomography data

The presence of noise is currently one of the main obstacles to achieving large-scale quantum computation. Strategies to characterise and understand noise processes in quantum hardware are a critical part of mitigating it, especially as the overhead of full error correction and fault-tolerance is beyond the reach of current hardware. Non-Markovian effects are a particularly unfavourable type of noise, being both harder to analyse using standard techniques and more difficult to control using error correction. In this work we develop a set of efficient algorithms, based on the rigorous mathematical theory of Markovian master equations, to analyse and evaluate unknown noise processes. In the case of dynamics consistent with Markovian evolution, our algorithm outputs the best-fit Lindbladian, i.e., the generator of a memoryless quantum channel which best approximates the tomographic data to within the given precision. In the case of non-Markovian dynamics, our algorithm returns a quantitative and operationally meaningful measure of non-Markovianity in terms of isotropic noise addition. We provide a Python implementation of all our algorithms, and benchmark these on a range of 1- and 2-qubit examples of synthesised noisy tomography data, generated using the Cirq platform. The numerical results show that our algorithms succeed both in extracting a full description of the best-fit Lindbladian to the measured dynamics, and in computing accurate values of non-Markovianity that match analytical calculations.


Introduction
A key challenge in developing medium to large-scale quantum architectures is error mitigation [1]. Unfortunately, noise adversely affects all stages of quantum computation, in particular the manipulation of states through quantum gates and continuous-time quantum processes. On large scale devices, error correcting codes will be used to suppress noise and achieve fully fault-tolerant computation. But in the current NISQ era, the overhead of full fault tolerance is prohibitive, placing it beyond reach of current or near-term hardware. It is thus important to understand and characterise the underlying noise dynamics in current quantum hardware, both in order to inform hardware design, and to prepare error correction and mitigation protocols optimised for the specific noise in the apparatus.
Various methods have been devised to evaluate and analyse noise in quantum dynamics, making different assumptions on, and providing different information about, the underlying noise processes [2]. One way to understand the underlying noise model in a device is to look for compatible Markovian 1 evolutions. Knowing the Lindbladian which best approximates the generator of the physical process may help to calibrate error mitigation techniques.
However, noise in quantum devices may substantially deviate from a memoryless dynamics, so that no compatible Markovian description exists. This can be problematic: strategies for mitigating or correcting non-Markovian noise are currently less established than schemes designed for memoryless errors. For example, the fault tolerance thresholds that have been calculated for non-Markovian noise models [5,6,7] have worse estimates than traditional threshold theorems. Therefore, methods to give precise quantification of Markovian and non-Markovian dynamics are of considerable interest. A number of approaches to benchmark non-Markovianity have been devised; we give a brief overview in subsection 1.2. However, some of these procedures are impractical to compute, requiring either complete knowledge of the master equation, or an unbounded number of full process tomography snapshots. Other methods are only capable of providing onesided witnesses of non-Markovianity, failing even in principle to identify all non-Markovian effects, and providing no information about the underlying dynamics in the Markovian case. Thus, finding a theoretically well-defined, but also feasible and low-resource procedure to assess Markovian and non-Markovian noise in quantum devices is desirable.
In this paper we present two methods, both built on convex optimisation programmes, to characterise and quantify noise processes in quantum systems. The first is an efficient algorithm to compute the Lindblad generator of the best-fit Markovian master equation to the measured dynamics. This can be used to certify Markovian evolution within any desired level of error tolerance, and the form of the resulting Lindbladian gives insight into the noise processes present. Alternatively, if the distance between the memoryless channel generated by the Lindbladian and the experimentally measured channel is significant (in comparison to the error rate of the tomographic reconstruction), this difference constitutes an insightful quantity to evaluate non-Markovian dynamics. The second algorithm calculates a quantitative and operationally meaningful measure of non-Markovianity, first proposed in ref. [8], in terms of the minimal amount of isotropic noise to be added to the generator of an hermiticity-and trace-preserving map, close to the tomographic snapshot, in order to "wash out" memory effects and render the evolution Markovian. 2 Conveniently, this task is shown to have an efficient solution in the required precision for fixed Hilbert space dimension in the form of a semi-definite integer programme [9]. A key strength of our approach is that the algorithms we have developed do not require any a priori assumptions on the underlying dynamics, nor any access to or characterisation of the environment, but only a small number of tomographic snapshots -a single one would already suffice.
We extend the approach of refs. [8,9] in a number of important ways. We change the semi-definite integer programme into a convex optimisation algorithm that searches for the best-fit Lindbladian generator in a neighbourhood of a given size around the input. This allows us to include an error tolerance parameter that can be tuned at any desired level, making the scheme robust with respect to inaccuracies of the tomographic measurements, and applicable to real-world data with statistical errors.
We also fully generalise the set of input operators, whereas ref. [8] assumed that the input channel had a non-degenerate spectrum in order to ensure the uniqueness of the generator. We lift this limitation in two ways. First, we show that we can perturb any operator with multi-dimensional eigenspaces into an arbitrarily close matrix having non-degenerate spectrum and at the same time retaining the hermiticity-preserving property. This guarantees the uniqueness of the generator while preserving the outcome of the non-Markovianity measure by adjusting the error tolerance parameter accordingly. Second, we consider the more physically relevant case, where we have a non-degenerate quantum channel which arises as a perturbation from a degenerate channel. This is crucial for characterising noise in quantum computing devices, where the channels of interest are typically noisy versions of quantum gates having degenerate eigenvalues. This task is delicate, due to the sensitivity of the Lindblad form to perturbation in the presence of multi-dimensional eigenspaces. To deal with this issue, we develop a mathematical approach to reconstruct these subspaces by leveraging techniques from matrix perturbation theory [10], and implemented this as a set of algorithms that serve as a pre-processing phase for the convex optimisation task. Finally, we extend the algorithm to the case where we have a sequence of tomographic snapshots. This offers a significantly more sensitive test of non-Markovianity, or a considerably more precise fit to an underlying Markovian master equation, with only a linear increase in the number of measurement settings required.
We accompany our theoretical analysis with a Python implementation of all the algorithms, which we benchmark numerically on simulated tomographic data in Cirq [11]. The numerical results show that our algorithms successfully identify Markovian-compatible dynamics for a range of 1-and 2-qubit examples, both for noisy quantum gates with degenerate spectrum and non-degenerate quantum channels. The numerics also confirm that our algorithms are able to compute accurate values for the non-Markovianity measure of noisy, non-Markovian quantum channels, which we show to be consistent with a calculation of this measure done by hand.
It should be noted that throughout this paper we restrict our attention for simplicity to time-independent Markovian noise. This is reasonable for characterising errors in individual gates or low-depth circuits on current hardware, given the short timescales involved [12] over which noise processes are unlikely to vary significantly. However, it is likely the approach presented in the current work can be extended to encompass the more general case of time-dependent Markovian noise, which would allow assessment of errors in longer quantum circuits, over timescales in which the noise processes may fluctuate. This is an important topic for follow-up work.
The paper is organised as follows. In the remainder of this section we give an overview of our main results, and of prior work in the field. In Section 2 we introduce our notation and cover preliminaries on convex optimisation and matrix perturbation theory. The precise notion of Markovianity is defined in Section 3, where we also present the non-Markovianity measure first defined in [8]. In Section 4 we discuss how we extend the approach from [8] to all quantum channels, in particular ones having degenerate spectrum. Our algorithms are presented in Section 5, and finally results from numerical simulations are discussed in Section 6.

Main results
Our main result is a set of algorithms, implemented in Python and benchmarked on simulated data in Cirq, which characterise and quantify the noise processes in a quantum system from one (or a small number of) tomographic snapshot(s). In the case of Markovian or near-Markovian dynamics, Algorithm 1 computes a full description of the Lindbladian master equation that best fits the measured data, without requiring any a priori assumptions on the form of the master equation. Conversely, Algorithm 2 calculates an operationally meaningful measure of non-Markovianity, corresponding to the minimum amount of white noise that has to be added in order to "wash out" the memory effects and render the dynamics Markovian.
More precisely, consider a d 2 × d 2 matrix M resulting from process tomography on the dynamics of a d-dimensional quantum system, together with an error parameter ε which accounts for statistical errors and other sources of uncertainty in the tomographic data, by setting the maximum distance from M to look for a compatible Markovian map.
For simplicity, assume for the moment that M has non-degenerate spectrum. (Note that a significant part of this work concerns lifting this assumption, so that the algorithms accept arbitrary matrices as input. See Section 4 and subsection 5.1.) By formulating a convex optimisation programme whose constrains are exactly the necessary and sufficient conditions of a Markovian generator (cfr. Section 3 for an in-depth explanation), we retrieve the closest Lindbladian L ( m) to the complex branch L m of the matrix logarithm L 0 = log M . By iterating the convex optimisation task over the set of branches L m for m ∈ {−m max , −m max +1, . . . , 0, . . . , m max −1, m max } ×d 2 , we keep the Lindbladian L ( m) whose generated Markovian map exp{L ( m)} is the closest to M .
If no Lindbladian is found within the ε-ball, with Algorithm 2 we compute the non-Markovianity measure of M in the sense of [8], denoted by µ min . This value quantifies the minimal amount of isotropic noise to be added to the generator H of an hermiticity-and trace-preserving map in the ε-neighborhood of M in order to obtain a Markovian evolution; in other words, µ min is the smallest scalar µ such that L = H − µ ω ⊥ is a valid Lindbladian. This is equivalent to (cfr. Section 3) with the constraint that there exists a branch L m such that where Γ is the involution from the transfer matrix representation in the elementary basis, to the Choi representation (more details in subsection 2.1), X = (H ) Γ , ω ⊥ is the projection onto the orthogonal complement of the maximally entangled state and δ is the maximal distance for H from L m representing -in some approximation -the radius of the ε-ball around M under the matrix logarithm (more details in subsection 5.2).
We present here the pseudo-code of a simplified version of the two core algorithms. They will be re-stated in complete form in subsection 5.2.

Algorithm 1: Retrieve best-fit Lindbladian
j=1 m j P j (branches of the matrix log L 0 ) Run convex optimisation program on variables µ and X: minimize if H is not null then Output: hermiticity-and trace-preserving channel generator H , non-Markovianity parameter µ min end We have tested our approach on simulated tomographic data from the Cirq platform for the following examples: the X Pauli gate, a 1-qubit depolarising channel, a 1-qubit unital quantum channel, the 2-qubit ISWAP gate, and a 2-qubit non-Markovian depolarizing CZ-channel. These examples cover 1-and 2-qubit Markovian and non-Markovian channels, as well as channels with degenerate eigenspaces. The algorithm was successful in analysing all these examples. For the nearly-Markovian examples, the algorithm correctly identified a Markovian evolution within the error tolerance regime and extracted the best-fit Lindblad generator. For Markovian channels with degenerate spectrum, we also investigated how the number of random samples from the pre-processing phase affects the output.
In Fig. 1a we show how the distance of the closest Markovian map to the snapshot varies with increasing number of iterations. The optimal output is reliably achieved within approximately 3000 runs. See subsection 5.2 for a full discussion of the randomised pre-processing stages to our convexoptimisation algorithms.
(a) The distance between the tomography result and the Markovian channel constructed by Algorithm 1 against random samples for the simulated ISWAP channel.
(b) Numerics for the non-Markovianity measure from Algorithm 2 for a simulated unital quantum channel with the error tolerance parameter ε in the range (0.025, 0.609). For the unital channel, Fig. 1b illustrates the output of the algorithm for different values of ε. A µ-value is output when an hermiticity-and trace-preserving genenerator is first found within the ε-neighbourhood (at ε = 0.025). The µ-value decreases continuously until a discontinuity at ε = 0.58 indicates that a valid Lindbladian generator is retrieved, without any need of white noise addition. In subsection 6.1.4 we show that the numerical measure precisely matches the value calculated by hand.

Related work on assessing non-Markovian noise
The nature of non-Markovian noise has been investigated from a variety of perspectives and with a number of different approaches (cfr. review pa-pers [13,14,15] and the introduction of [16]). One of the principal ways to quantifying non-Markovianity is based on divisibility [17,18], i.e., the property of a channel encoding evolution from time t 0 to t 1 to be implemented as a concatenation of channels from time t 0 to t and t to t 1 for any t ∈ [t 0 , t 1 ]. Indeed, a channel which is Markovian, in the Lindblad sense, is also divisible. However, the converse is not necessarily true. 3 The original measure of non-Markovianity [8], on which this work builds, determines whether the observed tomographic data is consistent with a timeindependent Markovian master equation. If not, it provides a quantitative measure of how far the observed dynamics is from the closest Markoviant trajectory. This task (and also that of determining finite divisibility) was shown to be NP-hard in general [9,19], but efficiently (classically) computable for any fixed Hilbert space dimension.
Other methods for detecting and measuring non-divisibility and non-Markovianity are based on checking monotonicity of quantities that are known to decrease under completely positive maps such as quantum correlation (see [20], known as the RHP measure), quantum coherence [21], quantum relative entropy [22] or the quantum mutual information [23].
Another approach affirms that a non-Markovian map is one that allows information [24], e.g. the Fisher information [25], to flow from the environment to the system. Non-Markovianity can also be quantified by considering the change in the distinguishability of pairs of input states [26]. The observation of non-monotonic behaviour of channel capacity [27], the geometrical variation of the volume of the set of physical states [28], ensembles of Lindblad's trajectories [29] and deep-neural-network learning [30] are among the many alternative strategies to quantify dynamics with memory effects.
A limitation of many of these measures is that they provide one-sided witnesses of non-Markovianity, but cannot show that the dynamics is Markovian, or find the master equation consistent with or closest to the observed dynamics.

Notation and preliminaries
We denote elementary basis vectors by |e j = (0, . . . , 1, 0, . . . , 0) T with 1 in the j-th position. The maximally entangled state is |ω = d j=1 |e j , e j / √ d and ω ⊥ = 1 − |ω ω| is the projection onto its orthogonal complement. We write F to denote the flip operator interchanging the tensor product of elementary basis vectors, i.e. F|e j , e k = |e k , e j . We will use the Frobenius norm on matrices, defined by M F = j,k |m jk | 2 , and the 1-norm, M 1 = j,k |m jk |. They are both submultiplicative, that is, they satisfy AB ≤ A B for all square matrices A and B.

Channel representations and ND 2 matrices
We will consider quantum channels of finite dimension only, i.e. completely positive and trace preserving (CPT) linear operators acting on the space of d × d matrices, where d = 2 N . To represent a channel T as a d 2 × d 2 matrix T , we will adopt the elementary basis representation: This is sometimes called the transfer matrix of the channel in the elementary basis. The corresponding representation |v ∈ C d ⊗ C d of a d × d matrix V on which the channel acts is then In this representation, the action of the channel on a matrix becomes matrixvector multiplication, and the composition of channels corresponds to the product of their respective matrix representations.
In order to formulate necessary and sufficient conditions for the generator of a Markovian evolution, we will also make use of another representation, the Choi-matrix (or Choi-representation), defined as Conveniently, the two representations are directly related through the Γinvolution [8], acting on the elementary basis as |e j , e k e , e m | Γ := |e j , e e k , e m | .
Explicitly, we have The Choi-representation is very useful to investigate the hermiticity-preserving property of quantum channels, i.e., quantum maps T such that T (X † ) = T (X) † for all X. Indeed we have Lemma 1. T is hermiticity-preserving ⇐⇒ τ is hermitian.
We will use the terms hermiticity-preserving and Choi-hermitian interchangeably. Using the terminology from ref. [31], we define the following two properties Definition 2 (defective matrix). A matrix is defective if the geometric multiplicity of some eigenvalue is strictly less than its algebraic multiplicity. That is, the matrix is not diagonalizable.

Definition 3 (derogatory matrix).
A matrix is derogatory if some eigenvalue has geometric multiplicity strictly larger than one.
For brevity, we call non-defective, non-derogatory matrices ND 2 matrices. These matrices are diagonalisable by a unique choice of normalised eigenvectors. Crucial for our approach is the consequent structure of the complex matrix logarithm of an ND 2 matrix T , given by an infinite number of branches indexed by a vector m = (m 1 , . . . , m d 2 ) ∈ Z d 2 . The 0-branch of the matrix logarithm of a diagonalisable matrix T = d 2 j=1 λ j |r j j |, with {λ 1 , . . . λ d 2 } being the eigenvalues of T and j , r j the respective left and right eigenvectors such that j |r k = δ jk , is given by the m-branch is then Modified hermitian adjoint The hermitian adjoint operation on a matrix in its vector representation |v in the elementary basis, is given by |v † := F|v * . We call vectors such that |v = |v † = F|v * self-adjoint, and we say that two vectors v and w are hermitian-related if |w = |v † = F|v * (and equivalenty |v = |w † = F|w * ). This terminology is unusual with respect to the conventional hermitian conjugation operation on vectors, but the definition adopted here for d 2 -dimensional vectors exactly corresponds to the usual hermitian adjoint of their corresponding d×d matrices on which T acts. We will conversely denote the standard hermitian conjugation for a matrix A by A H .

Convex optimisation programmes
At the heart of our algorithms are convex optimisation programs which either retrieve the closest Lindbladian to the matrix logarithm, or the smallest non-Markovianity parameters (both objects explained into detail in Section 3). These convex optimisations over a (scalar or vector) variable x have the general form [32] standard form where f 0 , f 1 , . . . , f n are convex functions. A fundamental property of convex optimisation problems is that any locally optimal point is also globally optimal.
A class of convex optimisation problems, called second-order cone programs, takes the form where f , g, b j 's and c j 's are vectors (not necessarily of the same dimension), F and A j 's are matrices and d j 's are scalars. The condition A j x + b j 2 ≤ c j , x + d j is called second-order cone constraint. Important for our work is the fact that minimization objectives for the Frobenius norm over matrix variables can be converted into a second order cone program via epigraph formulation.
To numerically implement convex optimisation programme, we make use of the Python library CVXPY [33,34].

Matrix perturbation theory
In order to account for a parameter ε reflecting the error tolerance with respect to the inaccuracy of the tomographic measurement, we make use of the following techniques and results in perturbation theory.
We will call eigenspace of A a subspace of a matrix A related to a single eigenvalue, and simple invariant subspace an invariant subspace whose corresponding eigenvalues are all distinct from the ones of its complement subspace.

Perturbation of eigenvalues
Given an n-dimensional eigenspace of an operator T with respect to the eigenvalue λ, an analytic perturbation of the form T (ε) = T + ε T (1) + ε 2 T (2) + . . . will create a λ-cluster of eigenvaluesλ 1 , . . . ,λ r (also referred to as λ-group in [35]) with multiplicity n 1 , . . . , n r so that r j=1 n j = n. The branches of perturbed eigenvalues are analytic functions with respect to the perturbation parameter ε. In our work, we deal with ND 2 matrices only, such that every eigenvalue has multiplicity 1; this also implies that the perturbation of a multi-dimensional eigenspace will remain a simple subspace, and more generally that any subspace of the perturbed operator T (ε) is simple. This will be important in order to leverage Theorem 4, given below.

Perturbation of invariant subspaces
To establish whether a given operator is a noisy implementation of a quantum Markovian map, we will make use of and slightly adapt the matrix perturbation framework developed in [10].
Let S be an n-dimensional invariant subspace of an d 2 × d 2 -dimensional operator M and let the columns of U 1 be a set of orthonormal vectors spanning it. Then, we extend it into a unitary matrix (U 1 , U 2 ), where the columns of U 2 are a basis of the orthogonal complement S ⊥ of S. Under this basis transformation, the matrix representation of M with respect to (U 1 , U 2 ) is given by since U H 2 M U 1 = 0 d 2 −n,n . It is usually more convenient to represent M in a block-diagonal form: for U 1 , U 2 , M 1 , M 2 as above, there always exist matrices Z 1 and Z 2 such that (cfr. Section V, Thm 1.5 of [10]) It follows immediately that the orthogonal complement of S is also a simple invariant subspace. The analysis on perturbed subspaces is based on the separation function for a simple invariant subspace S of M , which is an expression to quantify the distance between M 1 and M 2 . Now back to our task, we consider some hermiticity-preserving matrix Λ and a perturbation operator E such that M = Λ+E is the quantum channel which we want to either identify as or discriminate from a Markovian process.
Let Λ have an n-dimensional eigenspace with respect to the eigenvalue λ, and let U 1 be the matrix whose orthonormal columns span the n-dimensional invariant subspace of M with respect to the corresponding λ-cluster, producing the spectral resolution of M as per eq. (14). Under this basis, we write the representation of the perturbation matrix E in the block-form Then we have as a special case of Theorem 2.8 in Chapter V of Stewart and Sun [10], where here we consider the backward perturbation F = −E such that Λ = M + F and then use the property for matrix norms −A = A , for some submultiplicative matrix norm · , then there exists a uniqe matrix P with P < 2 · E 21 /γ such that the columns of form a basis of the λ-eigenspace of Λ.

Quantum Markovian channels and embedding problem
Processes which do not retain any memory of their previous evolution are called Markovian and satisfy the Markov property: given a sequences of points in time t 1 < t 2 <, . . . , < t n−1 < t n , a stochastic process X t taking values on a countable space has the Markov property if Extending this notion, a quantum Markov process is described by a oneparameter semi-group giving rise to a continuous sequence of completely positive and trace preserving (CPT) channels. The generators of this type of semi-group is called Lindbladian and must take the well-known Lindblad form [3,4], where H is hermitian, G = (g αβ ) is positive semi-definite and {F α } α are orthonormal operators. The first term on the RHS is the Hamiltonian part and describes the unitary evolution of the density operator, while the second term represents the dissipative part of the process. By diagonalising G as Γ := U † GU = diag(γ α ) and defining the so-called jump operators J α := √ γ α β u βα F β , we can re-write eq. (21) in diagonal form, The question whether a given quantum map M is compatible with a Markovian process, in the sense that there exists a memoryless evolution that at a certain time is equal to M, has been investigated from different perspectives, e.g. in the context of complexity [9], channel divisibility [17], regarding spectrum [36] and toward the goal of achieving a quantum advantage [37], and it is sometimes referred to as the embedding problem. A method to ascertain whether a given channel is compatible with a Markovian dynamics has been developed in ref. [8], which provides three properties for L that are necessary and sufficient to satisfy eq. (21) (where L is the elementary basis representation of L). These are (iii) ω| L = 0|, which corresponds to the trace-preserving property.
In this work, we restrict the analysis to time-independent Lindbladian generators, thus their corresponding quantum channel at time t is given by T (t) = e L t . We will call quantum embeddable any map whose matrix logarithm admits at least one complex branch satisfying these properties.
Interpreting the above conditions, we observe that they impose a rigid structure on the operator in matrix form, since they involve both eigenvalues and eigenvectors. In particular, from the hermiticity-preserving condition we note that, if λ is an eigenvalue of L and |v the corresponding eigenvector, then it follows that: Thus λ * and |v † = F|v * are an eigenvalue and the corresponding eigenvector of L too. This implies an important property of Linbladians in their elementary basis representation L : complex eigenvalues must necessarily come in complexconjugate pairs λ, λ * and have the same multiplicity. Moreover, the eigenspace of λ must admit a set of basis vectors whose hermitian conjugates span the eigenspace of λ * . If λ is real, then it must necessarily admit a set of vectors spanning its eigenspace which either come in hermitian-related pairs or are self-adjoint. For a Markovian channel T = e L the conditions for eigenvalues that are either complex or positive and for the vectors spanning their subspaces follow exactly the same rules. On the other hand, negative eigenvalues of T must have even multiplicity (implying that no non-degenerate negative eigenvalue can occur); the eigenspace of a negative eigenvalue must admit a basis of hermitian-related pairs of vectors. This particular structure for Markovian maps and their Lindbladian generators will be exploited in subsection 4.2 to reconstruct the eigenspace of an originally degenerate eigenvalue from a set of perturbed eigenvectors.
Conditions (i), (ii) and (iii) will run throughout our analysis, and we will implement them in the algorithms in Section 5 as constraints of a convex optimisation problem. We will also sometimes need to express these conditions in the Choi representation. By expanding condition (iii) as we can re-formulate it as This is equivalent to Tr 1 (L ) Γ = 0 in any matrix norm. A measure of non-Markovianity for hermiticity-and trace-preserving ND 2 quantum channels M was introduced in [8] in terms of white noise addition. More precisely, given a set of Choi-hermitian generators {L m } m of M with ω| L m = 0|, we define as the non-Markovianity parameter. Then µ min is the smallest value µ such that L = L m − µ ω ⊥ is a Lindbladian generator for some L m and is a measure of Markovianity for M .

Multi-dimensional eigenspaces
The analysis in [8] is based on the assumption that all eigenvalues are nondegenerate. Indeed, as pointed out there, the subset of ND 2 matrices is dense in the matrix set with respect to the Zariski topology, whose closed sets are the roots of the resultant of the characteristic polynomial and its derivative. Note that the Zariski topology is weaker than the metric topology, and hence any Zariski dense set is also dense in the metric topology [39]. Working with an ND 2 matrix is needed in order to deal with a set of eigenvectors where each of them is unique (up to a scalar factor). This also ensures that the matrix logarithm is unique, up to complex branches. Conversely, if M is not an ND 2 matrix, there is then a continuous freedom in the choice of eigenbasis, and thus uncountably infinitely many different matrix logarithms (not just the countable infinity of complex branches of the logarithm). This is due to the fact that, in case the matrix is diagonalizable but has an eigenvalue which is not simple, the corresponding degenerate eigenspace allows for infinite number of choices of basis vectors. In the more general case, when M is not diagonalizable, the Jordan canonical form again admits an uncountably infinite number of choices of generalized eigenvectors [40].

Perturbation of hermiticity-preserving matrices
In this section we show that it is always possible, given a defective or derogatory Choi-hermitian matrix, to produce an arbitrarily close ND 2 matrix that preserves the hermiticity-preserving property. Formally, Theorem 5 (ND 2 matrices are dense in the Choi-hermitian matrix set).
Let M be an hermiticity-preserving matrix, either defective or derogatory (or both). Then for any there exists an hermiticity-preserving This result allows us to resolve the problem of the freedom of basis choice and reduce to a unique principle branch of the matrix logarithm. In real tomographic data, eigenvalues will invariably be non-degenerate. But we must nonetheless deal with this case, as there may be channels with degenerate eigenvalues within the statistical error of the data. This is particularly important in the case of noisy implementations of unitary quantum gates for quantum computation, which often have degenerate eigenvalues.
For the proof, we will make use of the argument presented in ref. [41] showing that the set of non-defective, non-derogatory matrices is dense to prove that the same is true when we restrict to the subset of hermiticitypreserving matrices. In particular, we will apply the following definition and results.
Definition 6 (Resultant of two polynomials). Let p(x) = a n x n + · · · + a 1 x + a 0 and q(x) = b m x m + · · · + b 1 x + b 0 be two polynomial in x of degree n and m, respectively. The resultant of p and q, denoted by Res(p, q), is then given by the determinant of the (n + m) × (n + m) Sylvester matrix of p and q.
We will need the following two results: where {r j } j and {s k } k are the roots of p and q, respectively We note that the characteristic polynomial p of M has multiple roots if and only if p and its derivative p , have at least one common root; according to Corollary 8, this means that Res(p, p ) = 0. Note that Res(p, p ) is a polynomial in α; because of the fundamental theorem of algebra, we then know that it can only have a finite number of roots or that it is equal to 0 for all α. We can immediately rule out the second case, since we know that for α = 1 we have M = E, which has no degenerate eigenvalues. This implies that there exists some value α min so that for every 0 For our specific task, Theorem 5 implies that we can choose a perturbation matrix E and a perturbation parameter α so that we can obtain a ND 2 matrix from any input channel, such that any (possibly defective or derogatory) quantum embeddable channel will be recognized as compatible with Markovian dynamics by our algorithm even after being perturbed, given a tolerance parameter ε ≥ α ( M F + E F ). Analogously, our algorithm running for the perturbation of a map with a Markovianity parameter µ will retrieve a channel with a parameter smaller or equal to µ. Therefore, without loss of generality, we can substitute any matrix with a sufficiently close ND 2 matrix and take its logarithm instead.
A simple way to construct a Choi-hermitian ND 2 perturbation operator E is to first consider a matrix D with dim D = dim M , diagonal in the elementary basis, with the constraint where d (j,k) is the diagonal element in the (j · √ d + k)-th row and column. Then, E := D + FD * F is hermiticity-preserving and ND 2 .

Reconstructing perturbed eigenspaces
The converse situation is when we are given an ND 2 matrix M which may come from the perturbation of a quantum embeddable operator having degenerate subspaces. For instance perturbations of many of the standard unitary gates in quantum computation, such as Pauli gates. This is a delicate situation since in the general case the hermiticity-preserving structure will be broken under perturbation due to the instability of the basis of multi-dimensional eigenspaces already discussed. Thus, when looking for the closest Lindbladian, the convex optimisation approach will in general retrieve a Lindbladian whose matrix exponential is very distant from the original unperturbed operator M .
To illustrate this argument, consider the X-Pauli gate and restrict our attention to the hermiticity-preserving condition (noting that the closest hermiticity-preserving matrix will always be closer than the closest Lindbladian since the latter imposes more constraints). The operator X has a twofold degenerate eigenvalue 1 with eigenspace span{(1, 1, 1, 1); (1, −1, −1, 1)} and another two-dimensional eigenspace span {(1, −1, 1, −1); (1, 1, −1, −1)} with respect to eigenvalue −1. We denote these vectors by w 1 , w 2 , w 3 , w 4 , respectively. Observe that the first two are self-adjoint, while w 3 and w 4 are hermitian-related. Now consider a perturbation along w 1 , w 2 . Define the self-adjoint eigenvectors |w 5 = (|w 3 + |w 4 )/ w 3 + w 4 and |w 6 6 }, and its logarithm log(X + E) has eigenvalues ε, −ε, iπ − ε, iπ + ε (up to first order in ε) with respect to the same eigenbasis. At this point, if we look for the closest hermiticity-preserving operator, since all eigenvectors are self-adjoint we obtain a matrix having again the same eigenbasis and keeping the real part of the eigenvalues of log(X + E), i.e., ε, −ε, −ε, ε. Clearly, the exponential of this matrix is close to the identity map and not the expected Pauli-X gate. The same will apply for any complex branch of log(X + E) where we can add 2πi mod k to any eigenvalue of log(X + E).
If we instead consider a perturbation of the same magnitude but along hermitian-related vectors of the eigenspace of -1, say, E = ε |w 1 w 1 | − |w 2 w 2 | + |w 3 w 3 | − |w 4 w 4 | , then log(X + E) has again eigenvalues 1 + ε, 1−ε, −1+ε, −1−ε but this time with respect to eigenbasis In this case, the complex branch log(X + E) − 2πi |w 4 w 4 | has eigenvalues ε, −ε, iπ − ε, −iπ + ε and its closest hermiticity-preserving map has eigenvalues ε, −ε, iπ, −iπ. As expected, taking the exponential of this matrix will give a map very close to X. This example highlights both the importance of reconstructing a pair of hermitian-related eigenvectors for the eigenvalue -1 as well as searching over complex branches of the matrix logarithm.
Our strategy is to reconstruct a compatible hermiticity-preserving structure for those eigenvalues of M that are close to each others and that presumably come from a perturbation of a unique degenerate eigenvalue. Hence we will look for a new basis of eigenvectors that we will interchange with the actual eigenbasis, creating a new operator R on which to run the convex optimisation problem to retrieve the closest Lindbladian to log R.
We first discuss the single-qubit case and then generalise the approach for multi-qubit quantum channels.

One-qubit case
For single-qubit channels we have the following possiblilities: (a) one pair of close eigenvalues, (b) three close eigenvalues, (c) two different pairs, or (d) all four eigenvalues are close.
Consider case (a) with a pair of eigenvalues that is close to the real negative axis and where w 1 and w 2 are the corresponding eigenvectors. Assume that they come from a real negative 2-fold degenerate eigenvalue. We seek a new pair of hermitian-related eigenvectors {v, v † } such that span{v, v † } = span{w 1 , w 2 }. Thus we want to find coefficients α, β, µ, ν such that |v = α|w 1 The solution will parametrise a set of compatible hermitian-related eigenvectors that we will interchange with the vectors w 1 and w 2 .
If the two close eigenvalues are near the positive real axis, then if we assume they come from a real eigenvalue we have an additional option: a pair of two self-adjoint eigenvectors {v 1 , v 2 } spanning the eigenspace of w 1 and w 2 . In other words, we look for coefficients α, β, µ, ν such that |v 1 = α|w 1 + β|w 2 and |v 2 = µ|w 1 + ν|w 2 with |v 1 = ±|v † 1 and |v 2 = ±|v † 2 . Again, we should implement any possible solution {v 1 , v 2 } as a new basis of eigenvectors related to the pair of eigenvalues.
The third option for case (a) is a pair of complex eigenvalues not close to the real axis. If this was originally a unique, two-fold degenerate complex eigenvalue λ, the hermiticity-preserving condition implies a second two-dimensional eigenspace with respect to eigenvalue λ * ; this will then be case (c). Now consider case (b), where three eigenvalues of M are all close. In order to represent the perturbation of a quantum embeddable channel, they cannot originate from a 3-fold degenerate complex eigenvalue, since it is not possible to pair three eigenvalues each with a complex partner in a 4dimensional space. They cannot originate from a real negative eigenvalue either, since the same argument will apply for the logarithm of M , which must also be hermiticity-preserving. M may instead be compatible with a Markovian dynamics if the eigenvalues are close to the real positive axis. Denoting by w 1 , w 2 and w 3 the corresponding eigenvectors, we want to substitute them with a set {v, v † , z} with z = z † and span{v, v † , z} = span{w 1 , w 2 , w 3 } reconstructing an original unperturbed 3-dimensional eigenspace. Alternatively, we should find a new eigenbasis of three self-adjoint vectors.  In case (c) we say that we have two pairs of different eigenvalues. As we discussed in case (a), if one of this pair stems from a 2-fold degenerate and complex eigenvalue λ, then by the hermiticity-preserving condition the other pair should be close to λ * . We should thus find a basis v 1 , v 2 for the eigenspace of λ and v 3 , v 4 for the eigenspace of λ * such that v 3 = v † 1 and v 4 = v † 2 . Now consider pairs close to the real axis. If both can be associated with a real negative eigenvalue, then the basis of each 2-dimensional subspace corresponding to one pair of eigenvalues should be chosen to be hermitianrelated vectors as in case (a). If one pair is close to a real positive number and the other to a negative one, than the eigenspace of the positive eigenvalue can be spanned either by an hermitian-related pair of eigenvectors or two self-adjoint vectors. The last option in case (c) is that both pairs comes from two different 2-fold degenerate real positive eigenvalues. In this case, again each eigenspace can be spanned by a pair of hermitian-related vectors or two self-adjoint vectors.
If all four eigenvalues are close and we presume that they come from a single 4-dimensional eigenspace, then M must necessarily be the (perturbed) identity channel up to a real scalar factor. We check if this is close enough according to the error tolerance parameter.
In Table 1 we summarize the above described scenarios for the unperturbed operator.

Multi-qubit case
Assume that the multi-qubit channel has a cluster of n eigenvalues that are close to each others, presumably stemming from an n-dimensional eigenspace with respect an eigenvalue λ. We denote this subspace of M by S λ and we distinguish three possible cases.
If λ is complex, then as previously discussed the hermiticity-preserving condition implies the existence of a second n-dimensional subspace with respect to the eigenvalue λ * , S λ * . This remark already tells us that n ≤ d 2 /2; indeed, for the single-qubit case we have ruled out the possibility of a 3dimensional eigenspace for a complex eigenvalue. If we do not identify a second cluster of m eigenvalues that are close to λ * , then the channel M is not quantum embeddable. Otherwise, we look for a basis {v 1 , . . . , v n } of S λ such that {v † 1 , . . . , v † n } is a basis for S λ * . Given the eigenvectors {w 1 , . . . , w n } of the cluster of eigenvalues of M related to the perturbed eigenvalues originated from λ, and {z 1 , . . . , z n } related to λ * , we hence seek n vectors of the form |v := α 1 |w 1 + · · · + α n |v n such that v † ∈ span{z 1 , . . . , z n }. The set of m independent vectors satisfying the condition will be chosen as a new eigenbasis of S λ and their hermitian counterparts as the eigenbasis of S λ * . The algorithm to retrieve the closest Lindbladian should then run over all feasible solutions of this eigenbasis problem.  Table 2: Structure of an n-dim eigenspace of a multi-qubit hermiticity-preserving operator. Here we abbreviate "hermitian related" by h.r. and "self-adjoint" by s.h.
If λ is negative, we recall that in this case the relevant constraint is that log M is Choi-hermitian. The related perturbed subspace of M then needs n/2 vectors v 1 , . . . v n such that {v 1 , v † 1 , . . . , v n/2 , v † n/2 } is a basis of S λ . This implies that n must necessarily be even, and indeed in the discussion of the one-qubit case we ruled out the possibility of having 3 eigenvalues stemming from a degenerate negative eigenvalue. Thus, again denoting by {w 1 , . . . , w n } the eigenvectors of the perturbed eigenvalues of λ, we seek m/2 vectors |v := α 1 |w 1 + · · · + α n |v n such that v † ∈ S λ . The set of vectors given by the n/2 pairs {v, v † } will be chosen as the new eigenbasis. Again, the algorithm for the Lindbladian should run for all feasible sets of hermitian related eigenbasis for the subspace S λ .
The remaining option is a real, positive unperturbed eigenvalue λ generating a cluster of n perturbed eigenvalues of M . Here we have an additional freedom due to the possibility of admitting self-adjoint eigenvectors. More precisely, for each p = 0, 1, . . . , n/2, we seek an eigenbasis of p pairs of hermitian-related vectors and n − 2p self-adjoint vectors. Again, in principle all sets of vectors with this structure should be used for the algorithms, for all values of p.
We schematically illustrate this analysis in Table 2.

Approximate hermiticity-preserving eigenspace structure
When perturbations "mix" eigenspaces in a way such that it is not anymore possible to obtain an exact choice of vectors according to the prescription given above, we should search for vectors that are close to satisfying those conditions. The motivation comes from the following analysis based on the relation between the set of vectors spanning an unperturbed eigenspace and its perturbed version, where we will make full use of the tools presented in subsection 2.3. We show that the hermiticity-preserving structure is a stable property with respect to perturbations. Formally, If λ is complex, then there exist a set of basis vectors {w i } n i=1 spanning the right invariant subspace of M with respect to the λ-cluster and a set of basis vectors {w i } 2n i=n+1 spanning the right invariant subspace of M with respect to the λ * -cluster such that where X 2 is the submatrix of the basis transformation (T 1 , X 2 ), (X 1 , T 2 ) H for the spectral resolution of M for the subspace of the λ-cluster and its complement subspace, E 21 = T H 2 ET 1 , Y 2 is the submatrix of the basis transformation (S 1 , Y 2 ), (Y 1 , S 2 ) H for the spectral resolution of M for the subspace of the λ * -cluster and its complement subspace, In principle, every basis set satisfying the bounds of Theorem 9 should be used as a new eigenbasis of R for the Lindbladian algorithms. However, due to the fact that perturbation terms E 12 , E 11 and E 22 cannot be evaluated precisely, in practice the expression in eq. (31) can only be estimated according to a guess on the perturbation matrix E reflected by the tolerance parameter ε. Nevertheless, when an exact hermiticity-preserving structure cannot be retrieved, this theorem motivates the approach of searching for a basis set made of pairs of eigenvectors that are close to being hermitianrelated, instead of performing a brute-force search over all basis transformations in the subspace.
Proof of Theorem 9. In this proof we will make use of the vector 1-norm and matrix 1-norm only, and denote them simply by · to ease the notation. Assume that Λ is an hermiticity-preserving operator such that Λ = M − E. Let {w j } n j=1 be the set of eigenvectors of M related to the cluster of n eigenvalues {λ j } j stemming from an n-degenerate eigenvalue λ < 0 of Λ, and arranged in columns forming a d 2 × n matrix W 1 . Let {w j } d 2 j=n+1 be the remaining eigenvectors with respect to the remaining eigenvalues {µ j } j arranged in the columns forming the d 2 × (d 2 − n) matrix W 2 . We can write Consider a new basis of orthonormal vectors for the subspace of the λ-cluster {u j } j given by |u j = n k=1 ς jk |w k , j = 1 . . . , n, arranged in columns to form a matrix U 1 , together with a Under this basis transformation, we write and With this choice of basis and matrix representation, we can apply Theorem 4 in subsection 2.3: the eigenspace of λ is spanned by the columns {v j } j of the d 2 × n matrix V 1 = U 1 + Z 2 P for a unique operator P with Now, since Λ is hermiticity-preserving, then there are n/2 vectors of the form |ṽ = n j=1 α j |v j such thatṽ † is in the same subspace, i.e., |ṽ † = n j=1 β j |v j for some {β j } j . If the α's and β's are a solution for this relation, then for a k := j α j ς jk and b k : The first term in eq. (39) is 0 by construction. Moreover, the 1-norm is invariant under hermitian conjugation, i.e.
As discussed above, if λ > 0 then self-adjoint vectors should be taken into account. Assume the same structure as in eqs. (34) and (35) again with |u j = n k=1 ς jk |w k , j = 1 . . . , n orthonormal vectors and with the columns of V 1 = U 1 + Z 2 P spanning the eigenspace of λ. Assume that there exists a set of coefficients {α j } j such that |ṽ = j α j |v j andṽ † =ṽ. Then a bound follows in the same fashion as eq. (43), where again a k = j α j ς jk . An analogous bound can be derived for a complex eigenvalue λ, where the condition should also account for a partner eigenspace with respect to λ * . In this case, we have two different basis choices to consider. Let the columns of T 1 span the subspace of the λ-cluster and let (T 1 , T 2 ) be unitary and such that Similarly, let the columns of S 1 be a basis of the subspace for the λ * -cluster so that with (S 1 , S 2 ) unitary. Then the columns of L 1 = T 1 + X 2 P span the eingenspace of λ and those of R 1 = S 1 + Y 2 Q span that of λ * . Additionally, assume that the eigenspace of λ admits a basis of vectors whose hermitian counterparts span the eigenspace of λ * , that is, there exist coefficients α's and β's such that j (α j | j ) † = j β j |r j . Now, let {w k } n k=1 denote the eigenvectors of the λ-cluster, and let {w k } 2n k=n+1 denote those of the λ * -cluster and write |t j = j ς jk |w k and |s j = j ς j,n+k |w n+k . Then where this time a k := j α j ς jk and b k := j β j ς j,n+k .

Algorithms for retrieval of best-fit Lindbladian and computing non-Markovianity measure µ min
In this section we present convex optimisation based algorithms to verify whether, given a channel M , there exists a Markovian channel in the εneighborhood (ε can be understood as the error-tolerance parameter which captures the error in the tomographic data as well as our maximal tolerance in the amount of approximation error). Alternatively, when such a channel does not exist, we look for the generator of an hermiticity-and trace-preserving channel in the δ-Ball of log M (induced again by the ε-ball of M ) with the smallest non-Markovianity parameter µ [8].
Since the M we are interested in arise as tomographic data, and ND 2 matrices form a dense set in the set of matrices, we assume that the input for the algorithms is ND 2 . However, following the argument of subsection 4.1 it is straightforward to include the possibility of inputs which are not ND 2 . This would just require an additional step to check the input, and if it was not ND 2 perturb it to a matrix that is.
Our main algorithm, which provides a recipe for evaluating non-Markovianity, is given in Algorithm 3. It takes as input an estimate of the channel, M , a precision parameter, p, an accuracy parameter ε, and an integer, r, which determines how many random basis choices will be tested in the case of degenerate eigenvalues. If a Markovian channel, T , exists within the εneighbourhood of M the algorithm returns the Lindbladian L satisfying T = e L . If no such Markovian channel exists, the algorithm returns the non Markovianity parameter µ min .
The first step in the algorithm is to determine whether the input matrix has eigenvalues that may originate from a perturbation of degenerate eigenspaces. If the matrix has no such eigenvalues, we run the convex optimisation algorithm directly on the input. Conversely, if the matrix has eigenvalues that may originate from perturbed degenerate eigenspaces, before running the convex optimisation we need to produce a matrix with perturbed degenerate eigenspaces having a suitable hermiticity-preserving structure, as discussed in subsection 4.2. The task is performed by Algorithm 4, Algorithm 5 and Algorithm 6, outlined in subsection 5.1.
Once these matrices have been constructed we run the convex optimisation Algorithm 1 on them to determine whether or not there exists a Markovian channel in the ε-neighbourhood. If a Markovian channel is found, the algorithm terminates and outputs the Lindbladian, along with the distance between the input channel and the resulting Markovian channel. If no Markovian channel is retrieved, the algorithm calls Algorithm 2 which calculates the non-Markovianity parameter µ min (see subsection 5.2).
Finally, in subsection 5.3 we present a version of the algorithm extended for the case of a sequence of snapshots from a given quantum process. We look for a Lindbladian whose generated evolution passes ε-close to each snapshot at the appropriate time.

Pre-processing algorithms
The first step in the pre-processing algorithm is to construct sets of real positive, real negative, and complex degenerate eigenvalues. If there are no degenerate eigenvalues then the algorithm calls the convex optimisation algorithms to find L or µ min (see subsection 5.2). If all the eigenvalues are real, positive and equal (to within precision p) then the pre-processing algorithm outputs that the channel is consistent with the identity map.
In all other cases the next step of pre-processing is to construct bases for the degenerate eigenspaces which have a compatible hermiticity-preserving structure (see the analysis in subsection 4.2).
The algorithm to construct the basis for a degenerate eigenvalue which is either complex or negative real is given in Algorithm 4. It takes as input the transfer matrix of the channel, M , and two sets of integers. The first set of integers, degenerate_set, contains the indices of the eigenvectors corresponding to the degenerate eigenvalue. The second set of integers, conj_degenerate_set, contains the indices of the eigenvectors corresponding to the conjugate degenerate eigenvalue. 4 The idea behind Algorithm 4 is that for each set of eigenvectors w 1 , w 2 ,...,w n associated with a degenerate eigenvalue λ and eigenvectors u 1 , u 2 ,...,u n associated with λ * we solve the equation: for α, β,...γ, which allows us to construct a basis with the correct hermiticity preserving structure. We solve eq. (49) by constructing the matrix: and finding its kernel. If Nullity(A) = dim(ker(A)) is equal to |degenerate_set| then there is a hermiticity-preserving basis which spans the degenerate eigenspace (in fact there are uncontably infinitely many choices of hermiticypreserving bases). The algorithm returns the basis vectors: j is the value of α j in the i th solution to eq. (49).
The pseudo-code to construct a basis for a degenerate, positive real eigenvalue is given in Algorithm 5. It follows a similar idea to Algorithm 4, with some additional processing to account for the possibility of self-adjoint eigenvectors.
In the case of a real, positive, degenerate eigenvalue the equation to solve in order to construct a basis with the correct hermiticity-preserving structure is: for α, β,...γ. As in the case of complex or negative real eigenvalues, if Nullity(A) is equal to |degenerate_set| then there is a hermiticity preserving basis which spans the degenerate eigenspace given by: However, now the basis may include self-adjoint eigenvectors, as well as Hermitian conjugate eigenvectors. These need to be treated separately in the next stage of the pre-processing.
To verify if the i th basis-vector is self-adjoint we check whether α If the total number of non-self-adjoint basis vectors is even, the algorithm returns two sets of basis vectors -the self-adjoint set and the nonself-adjoint set. If the total number of non-self-adjoint basis vectors is odd, the algorithm finds the non-self-adjoint basis vector for which j α j is minimised, and returns this with the set of self-adjoint basis vectors instead of the non-self-adjoint set.
There are an infinite number of possible basis choices which respect the hermiticity preserving structure, and which basis is chosen will affect the distance to the closest Lindbladian (see subsection 4.2). To handle this, we use a randomised construction to generate r bases which respect the hermiticity preserving structure, and we run the convex optimisation algorithm on each basis, keeping the optimal result. The number of random bases, r, is an input to Algorithm 3. Running the algorithm with more random bases will give a better result (numerical results on degenerate 1and 2-qubit channels are given in Section 6).
The algorithm to construct a random basis is given in Algorithm 6. It takes as input M , lists indicating the indices of eigenvectors associated to each degenerate eigenvalue {sl λ | for degenerate λ}, and the bases constructed by Algorithm 4 and Algorithm 5. It returns a random choice of basis, S.
In order to construct S, Algorithm 6 checks whether i is in any of the sl λ for each i ∈ dim(M ). If it isn't, we set the i th column of S equal to the i th eigenvector of M . If it is we check whether the i th column needs to be constructed randomly from self-adjoint basis vectors, or constructed randomly from non-self-adjoint basis vectors or constructed by conjugating another column.
For λ real, which columns are constructed randomly (from self-adjoint or non-self-adjoint basis vecors) and which are constructed by conjugation is arbitrary. For λ complex, either all columns associated to λ are constructed randomly, or all are constructed by conjugating the basis vectors associated to λ * . This is determined by which order sl λ and sl λ * are given as input to Algorithm 4. If the column needs to be constructed randomly then we generate |sl λ | random seeds κ j . The i th column of S is then given by: where v j are the hermiticity-structure preserving basis vectors associated with λ, and the sum is over either all self-adjoint or all non-self-adjoint basis vectors. If the column needs to be constructed by conjugation we find which column it is conjugate to, denoted k, and the i th column of S is then given by: The final step of preprocessing is to compute R = SM S −1 for each random basis. We then run the convex optimisation algorithms (see subsection 5.2) on each R, and output the optimal result.

Core algorithms
At this stage, we are given a the tomographic snapshot M and an ND 2 matrix R from the pre-processing phase, (clearly, if M is ND 2 and has no cluster of eigenvalues, then R = M ).
Algorithm 1 (here below is the restatement of the first algorithm presented in the main results, but this time incorporating the matrix R resulting from the pre-processing) looks for the closest Lindbladian to log R in the ε-neighborhood of M by formulating a convex optimisation task whose constraints are exactly the necessary and sufficient conditions for a Lindbladian generator, as discussed in Section 3.

Algorithm 1 (restatement): Retrieve best-fit Lindbladian
The first remark about the algorithm is that we are not searching for the closest Markovian channel, but for the closest Lindbladian to its matrix logarithm. A natural question is then whether the two objects are precisely related, that is, if the closest Lindbladian generates the closest Markovian channel. In the general case, this is not true. A simple counter-example is provided by the perturbed X-gate discussed in subsection 4.2, where the closest Lindbladian generates the identity map, although the unperturbed X-gate is a closer Markovian channel. However, consider the upper bound for general matrices A and B (Lemma 11 in [9]) If we now ask B to be Lindbladian, then the closest Lindbladian L to A is also the operator minimizing this upper bound. This implies that the distance between the Markovian channel exp L retrieved by our algorithm and the input matrix M = exp L R is upper-bounded by where we assume exp C to be the closest Markovian operator to M . Thus, setting ε as our tolerance parameter for the algorithms, we are guaranteed to find a compatible Markovian evolution for the input M . Note that the convex optimisation problem contains a second-order cone constraint (when formulating the Frobenius norm minimization in epigraph form), a linear matrix inequality (LMI) constraint to ensure the conditionally completely positive condition and a first-order cone constraint representing the trace preserving condition. As in any convex optimisation program, every minimum of this program will be a global minimum. Moreover we have A F = A Γ F for A in the elementary basis representation. Thus, L − L F = (L − L) Γ F = X − L Γ F . This is useful in order to run the convex optimisation program with the variable X without involutions.
We iterate over different branches of the matrix logarithm and pick the resulting Lindbladian from the convex optimisation task whose generated channel is the closest to M .
If no Lindbladian generator is found in ε-ball of M , Algorithm 2 (already presented in the main results) will be called to find the channel in the εneighborhood with the smallest non-Markovianity parameter according to the definition in eq. (27). This can again be retrieved by formulating a convex optimisation programme.
The variable δ characterizes the neighborhood around L that maps under the exponential into the ε-neighborhood of M ; a lower bound is given by Lemma 11 in [9], that is, ε ≤ exp(δ) · δ · L 0 F . An upper bound is again provided in Corollary 15 of [9]. However, the boundaries of this region j=1 m j P j (branches of the matrix log L 0 ) Run convex optimisation program on variables µ and X: if H is not null then Output: hermiticity-and trace-preserving channel generator H , non-Markovianity parameter µ min end around L related to the ε-ball are not precisely characterised. This can cause problems for the convex optimisation programme, which in some cases can retrieve a Lindbladian for δ in the interval defined by the above bounds whose generated map falls outside the ε-ball. We further discuss this matter with a relevant example in subsection 6.1.3 (cfr. also Fig. 6). To solve this issue in a practical way, we run over increasing δ values up to 10 times the value determined by the lower bound.
The loop in both algorithms runs a brute-force search over the m-branches of the matrix logarithm. There are countably infinitely many such branches. However, Khachiyan and Porkolab [42] prove that setting m max = O(2 2 poly(d) ) always suffices to find a solution if one exists. In fact, the same authors prove that there is a polynomial-time algorithm for integer semidefinite programming with any fixed number of variables, so for any fixed Hilbert space dimension in our setting, based on the ellipsoid method (which is more efficient that performing brute-force search up to the upper-bound). This polynomial-time algorithm is important theoretically, but not practical. In practice, integer programming solvers use branch-and-cut methods [44].
However, the increased implementation complexity of these more sophisticated techniques is unlikely to be justified here. Instead, we implement a simple, brute-force search over the branches of the logarithm, ordered by increasing | m|. There are physical reasons why a naive, brute-force search is likely to work well here. Large values of m j correspond to high-energy / frequency components of the noise. It is unlikely that very high energy underlying physical processes play a significant role, and it is also unlikely that our tomographic data will be sensitive enough to resolve very high-frequency components. Therefore, if there is a Lindblad generator consistent with the tomographic snapshot, it is most likely to occur at low values of | m|. Numerical studies on synthesised examples of 1-qubit tomography (see Fig. 2) confirm that setting a small value of m max suffices in practice; indeed, even m max = 1 sufficies in all cases we tested. This is also corroborated by the numerics that was carried out in [8].
It would be straightforward to replace the brute-force search by a call to an integer program solver if this ever proved necessary.

Multiple snapshots
The approach for single snapshot can be extended to multiple tomographic snapshots taken at a sequence of times. Augmenting the number of measurements is a way to make the conditions for a compatible Markovian dynamics, or detection of non-Markovian dynamics, more stringent, since the requirement is that there exist a single time-independent Lindbladian that generates a dynamical trajectory that passes close to every snapshot.
Consider a sequence of tomographic snapshots M 1 , . . . , M q associated with measurement times t 1 , . . . , t q . In Algorithm 7 we formulate a convex (a) Distance between tomographic results and the Markovian channel output by our algorithm for an X-gate, running Algorithm 3 for 10,000 random samples for each value of m.
(b) Distance between tomographic results and the Markovian channel output by our algorithm for a depolarizing channel with p = 0.3. For each m-value, Algorithm 3 was run for 100 random samples.

Figure 2:
Investigating how the optimal Lindbladian found by our algorithm varies with the number of branches of the logarithm searched over for the X-gate and depolarizing channel. There is some fluctuation in distance in the case of the X-gate -this is due to the randomised nature of the algorithm (see Section 6 for full results on the X-gate).
optimisation program that finds the Lindbladian minimising the sum of the distances from the logarithms of the input matrices. (For simplicity, we present the algorithm for the case of ND 2 matrices with no cluster of eigenvalues, and discuss the more general case later in this section.) We then iterate over different branches and pick the Lindbladian L for which c M c − exp t c L F is the smallest. We also require that the distance with respect to any logarithm is not larger than some value δ varying in the same manner as in the single-snapshot case. We certify the evolution as Markovian if the exponential of the retrieved logarithm is within the ε-ball of each input tomographic snapshot.
Once again, dealing with degeneracies requires more work. If the perturbation of an n-degenerate complex eigenvalue is identified, then we should check whether this is consistent with all other snapshots, that is, all measurements show a cluster of n eigenvalues. Because of the hermiticity-preserving condition we also expect a partner cluster of n eigenvalues, corresponding to the complex-conjugate partner subspace. Note that the two clusters will overlap when they approach the negative axis, and conversely that any cluster close to the negative axis is expected to split in two partner clusters representing the perturbation of two eigenvalues related by complex conjugation: this means that when taking successive snapshots, we can avoid the case of perturbed negative eigenvalues by choosing suitable measurement times. Secondly, we should also verify if in all snapshots both the multi-dimensional subspace of a cluster and its partner subspace are consistent with the same unperturbed pair of hermitian-related eigenspaces, up to some approximation. Indeed, recall that eigenspaces are stable with respect to perturbation. Moreover, since here we are checking for time-independent Markovianity from multiple snapshots, we are interested in the case where the Lindbladians for all snapshots are the same.
Given this, and the results of subsection 4.2.3, our approach to handling an n-degenerate complex eigenvalue λ is as follows. We select one of the snapshots, c , at random, and apply the pre-processing steps from Algorithm 3 to M c in order to obtain a random hermitian-related basis of vectors for the λ and λ * subspaces. Denote these by {|v j } and {|v j † } respectively, and denote the projections onto the subspace spanned by these vectors by Π v and Π v † . We then determine whether this choice of basis is compatible with the other snapshots, by checking that where ς 1 is a tolerance parameter to be estimated from eq. (48). A set of basis vectors retrieved with the above approach can then be used as the new eigenvectors of the input matrices {M c } c for the clusters of eigenvalues, in analogous fashion as for the single-snapshot case, and then run Algorithm 7 on these modified operators. This procedure for constructing a random basis should be repeated a number of times, and the optimal result kept, as in the single snapshot case. Note, that which snapshot the pre-processing steps are applied to should be picked at random in each run of the algorithm.
A special case is the perturbation of degenerate real eigenvalues that do not turn into complex ones through the sequence of snapshots, corresponding to real eigenvalues for the Lindbladian generator. In this case, we do not have a partner invariant subspace. As discussed already, an eigenspace of an hermiticity-preserving operator with respect to a real eigenvalue admits self-adjoint eigenvectors in addition to hermitian related pairs. Consider an n-degenerate real eigenvalue κ. As in the complex case, we pick a snapshot c at random, and apply the pre-processing steps from Algorithm 3 to M c to find a basis for the subspace Π c (κ). This basis will be composed of a set of vectors {|v j } p j=1 , a set of hermitian-related vectors {|v j † } p j=1 , and a set of self-adjoint vectors {|s j } n−2p j=1 , for some value p = 1, . . . , n/2. Denoting by Π v , Π v † , Π s the projectors onto the subspaces spanned by these sets of vectors, the constraint eq. (58) turns into: where ς 2 is characterized by the bounds in eqs. (43) and (44). As in the complex case, we run Algorithm 7 on the input matrices resulting by using these vectors as the bases for the degenerate subspace in each M c , and then repeat the randomised process a number of times.
if L is not null then Output: Lindbladian L end

Numerical Examples with Cirq
In this section, we present the results of testing our algorithm numerically on noisy dynamics synthesised in Cirq [11]. The numerics serve as a benchmark of both our convex optimisation and pre-processing algorithms. Since we know the 'ideal' channel in each test case, we can benchmark the outcomes against the true values. The algorithms performed well in all cases. (It is worth emphasisiing that the algorithm itself does not require any information about what the 'ideal' channel is -it simply takes as input the synthesised noisy tomographic data. The ideal channel is only used to benchmark the results against.)

One-qubit numerics
In every one-qubit example each measurement in the simulated process tomography was repeated 10,000 times.

Unitary 1-qubit example: X-gate
In subsection 4.2 we demonstrate how the naive algorithm (with no preprocessing) is not guaranteed to find the closest Markovian channel if the input is a 'noisy' X-gate, as the degenerate eigenbases are not stable with respect to perturbations. 5 We used Cirq's density matrix simulator to simulate process tomography on a 1-qubit X-gate. Denoting the result of the process tomography by M X we found: Clearly we expect the closest Markovian channel to be at least as close to M X as this value. We then used our convex optimisation algorithm to construct the closest Markovian channel, T X . In Fig. 3 we show how the distance M X − T X varies with the number of random samples we allowed the code to run for. As expected, increasing the number of random samples decreases the distance between M X and T X , although there is some fluctuation due to the randomised nature of the algorithm. This shows that the 'direction' which the matrix is perturbed in is crucial -as we saw in subsection 4.2, perturbations of the same magnitude along different directions can have very different effects on whether a compatible Lindbladian is found. : Results for a simulated 1-qubit X-gate. In both cases the algorithm ran with ε = 1. A y-axis value greater than 1 indicates that no Markovian channel was found in that run.

Markovian 1-qubit example: depolarizing channel
The depolarizing channel, implementing the evolution: is an example of a non-unitary, but Markovian, quantum channel. We simulated process tomography on a depolarizing channel with p = 0.3. The transfer matrix for this channel has a non-degenerate +1 eigenvalue, and a three-fold degenerate eigenspace with eigenvalue 0.6. We used our convex optimisation algorithm to construct the closest Markovian channel. Denoting the result of the process tomography by M depol and the transfer matrix of the depolarizing channel by E depol , the algorithm found: In Fig. 4 we show how the distance between the tomography result and the closest Markovian channel varies with the number of random samples we allowed the code to run for. Comparing Fig. 4 with Fig. 3 (and noticing the different Y-axis scales in the two graphs) we see that unlike with the X-gate, Algorithm 3 always finds a good approximation to the tomographic result, even for very low random samples. Note also that the direction of perturbation is less important for the depolarizing channel than for the X-gate.

Non-Markovian and Markovian 1-qubit examples: unital quantum channel
In [45] Kraus operators are constructed for a unital quantum channel, and conditions for the channel to be Markovian are found. For a master equation: The Kraus operators for the evolution are: for: where {i, j, k} is a permutation of {1, 2, 3}. The inequality for the channel to be completely positive is given by: where {i, j, k} is a permutation of {1, 2, 3}. The condition for the channel to be Markovian is given by γ i (t) > 0, ∀t, i. We simulated process tomography in Cirq on a non-Markovian unital quantum channel with γ 1 = −200, γ 2 = 201, γ 3 = 200.5. The transfer matrix for this channel is non-degenerate. Denoting the tomography result by M unital , and the transfer matrix of the unital quantum channel by E unital we obtained: Applying Algorithm 3 with a large error tolerance parameter (ε = 1) found a Markovian channel, T , satisfying T − M unital = 0.609.
In Fig. 5 we show the results of comparing the non-Markovianity parameter µ against the accuracy parameter ε. We first ran Algorithm 3 with varying ε using a step size of 0.5 in the Algorithm 2 subroutine. The results are shown in Fig. 5a. For ε ∈ (0, 0.186) no non-Markovianity parameter was found by the algorithm. However in Fig. 5b, where the same analysis is run with a step size of 0.01 a non-Markovianity parameter was found for all ε ≥ 0.025. This can be understood by recalling that Algorithm 2 doesn't search for the channel T that minimises µ and is within an ε-ball of M . It searches for the Lindbladian L that minimises µ and is within a δ-ball of log(M ), then checks whether e L − M ≤ ε. Since the bounds on δ aren't tight, it repeats this process for δ ∈ {δ i = δ min + δ step | δ i < δ max }, and keeps the best result. If δ step is too large, this can lead to the algorithm failing to find a µ within the ε-ball of M , even when one exists. An illustration of this is given in Fig. 6. This is what has happened in Fig. 5a This problem can be eliminated by running Algorithm 2 with a smaller step size. In Figs. 5b and 5c we show the results of running Algorithm 3 with a varying ε parameter, and where the Algorithm 2 subroutine has a step size 0.01. For ε = 0 no µ is found by the convex optimisation algorithm, indicat- δ min δ max Figure 6: A schematic illustration of the distortion of the ε-ball (gray area) under the matrix logarithm, whose boundaries are included between two balls of radius δ min and δ max . Let T i (for i = 1, 2) be quantum channels with non Markovianityparameter µ i such that log(T i ) − log(M ) = δ i and T i − M = ε i where µ 1 > µ 2 , ε 2 > ε > ε 1 and δ 2 > δ 1 . If Algorithm 2 is run on the channel M with δ step > δ 2 −δ 1 the algorithm will not find a non-Markovianity parameter because the convex optimisation finds T 2 , but this is rejected by the step which checks if e L −M ≤ ε. Running Algorithm 2 with a smaller δ step solves this issue.
ing that there is no compatible Markovian channel in the ε-ball around M , even with white noise addition. Above ε ≈ 0.025 the convex optimsation algorithm consistently finds a µ. As ε approaches 0.609 µ decreases to zero, indicating that increasing the size of the ε-ball decreases the white noise needed to render the channel Markovian, as expected We also simulated process tomography on a Markovian version of the unital channel with γ 1 = 200. Running Algorithm 3 on the result of the tomography, we found a Markovian channel T satisfying T − M = 0.02.

Analytical derivation of the non-Markovianity parameter
In order to benchmark the output of Algorithm 2, we calculate by hand the value of µ for the unital quantum channel numerically investigated in subsection 6.1.3. This particular example is convenient because there is no cluster of eigenvalues; indeed, all eigenvalues of the simulated perturbed channel are positive and far apart from each other. This means that all subspaces are 1-dimension and thus we don't have to undergo a search over (infinitely many) feasible basis vectors for the unperturbed channel.
Noting that the white noise addition cannot influence either the hermiticity-preserving or the trace-preserving condition, our strategy will be to construct a matrix from the tomographic data by "filtering" these two properties, and then calculate the minimal value of µ to increase the eigenvalues of its Choi representation in order to satisfy the conditionally complete positivity condition. This approach does not guarantee that we will obtain the closest map compatible with Markovian dynamics through white noise ad-dition, but we can reasonably expect to retrieve an operator that is very close to the optimal one. We will then compare the value for µ from this calculation with the one returned by algorithm when setting ε equal to the distance between the input matrix and the matrix exponential of the map calculated with this analytical method.
Consider the constraints on the eigenvalues and eigenvectors of an hermiticity-and trace-preserving channel having non-degenerate and non-negative spectrum. Hermiticity implies that all eigenvectors must be self-adjoint; the trace-preservaing condition means that one eigenvalue must be equal to 0 with respect to a left eigenvector being proportional to the maximally entangled state ω| = (1, 0, 0, 1)|. By the biorthogonality conditions between left and right eigenvectors, this forces the three right eigenvectors related to the complement of the kernel to have first and fourth components of opposite sign. Hence, we will go through the following steps to turn the simulated perturbed unital channel U into an hermiticity-and trace-preserving matrix H and calculate its non-Markovianity measure by hand: (1) From the eigenvalues λ 1 > λ 2 > λ 3 > λ 4 , set λ 4 = 0 and define the matrix D := diag(λ 1 , λ 2 , λ 3 , 0).
(3) construct three vectors orthogonal to ω| whose first and fourth components are respectively v For the simulated channel, under this procedure we obtain a value of µ = 5.76983 and ε = 0.01788. As a comparison, by setting ε = 0.01788 in our algorithm (with δ step = 10 −5 ) this returns µ min = 5.76539819; also, as we can see in Fig. 5c, (0.012 -0.025) is indeed the first segment for ε where we retrieve a valid value for µ min . We refer to the accompanying Mathematica notebook in the Supplementary Material for the explicit calculation following this procedure.

Two-qubit examples
In every two-qubit example each measurement in the simulated process tomography was repeated 100,000 times.

Unitary 2-qubit example with degenerate eigenvalues: ISWAP gate
The action of the ISWAP gate is to swap two qubits, and introduce a phase of i to the |01 and |10 amplitudes. It has a six-fold degenerate eigenspace with eigenvalue +1, a two-fold degenerate eigenspace with eigenvalue −1, and two four-fold degenerate eigenspaces with eigenvalues ±i. We used Cirq's density matrix simulator to simulate process tomography on the ISWAP-gate. Denoting the result of the process tomography by M I , the algorithm found: We then used our convex optimisation algorithm to construct the closest Markovian channel, T I . We ran Algorithm 3 for 100 random samples, 85 times. The results are shown in Fig. 7.
It is clear that a higher number of random samples is required than in the 1-qubit case with the X-gate. This is to be expected, since the degenerate eigenspaces are higher-dimensional, and there are more of them. So the probability of randomly choosing a 'good' basis is lower.

Markovian 2-qubit example: depolarizing CZ channel
The CZ is a two-qubit quantum gate. Its action is to apply a Z-gate to the second qubit if the first qubit is in the |1 state. Otherwise it acts as the identity on both qubits.
We simulated process tomography on a depolarizing CZ channel, which applied the CZ gate with probability 0.1, or an XX / Y Y / ZZ -gate with probabilities 0.07,0.08 and 0.09 respectively. This is an example of a non-unitary, but Markovian, quantum channel. The transfer matrix of the channel has six two-fold degenerate eigenspace with eigenvalue +1, 0.7, 0.93, 0.57, 0.67, and 0.47, and non-degenerate eigenvalues of 0.68, 0.66, 0.48 and 0.46.
Denoting the result of the process tomography by M CZdepol and the transfer matrix of the depolarizing channel by E CZdepol the algorithm found: We used our convex optimisation algorithm to find the closest Markovian channel, T CZdepol . We ran Algorithm 3 for 100 random samples, 10 times. The results are shown in Fig. 8. Comparing Fig. 8 with Fig. 7 (and noticing the different Y -axis scales in the two graphs) we see that unlike with the ISWAP-gate, Algorithm 3 always finds a good approximation to the tomographic result, even for very low random samples. As in the 1-qubit example, this shows that the direction of the perturbation is less important for the non-unitary channel than it is for the unitary channel.
(a) The distance between the tomography result and the Markovian channel constructed by Algorithm 3 in each 100random sample run.
(b) Results from the same simulation as Fig. 8a, where now we show minimum distance achieved against total number of random samples.

Conclusions
We have developed novel methods and algorithms, based on previous work [8], to retrieve the best-fit Lindbladian to a quantum channel. We have implemented these algorithms in Python, and benchmarked them on synthetic tomography data generated in Cirq.
The key strengths of our method is that it can be applied to a single tomographic snapshot, is completely assumption-free regarding the structure of the analysed operator, and does not rely on any prior knowledge of the environment or the noise model. At the core of the method is a convex optimisation programme which searches for the closest Lindbladian generator within a given distance from the matrix logarithm of the input. This approach is successful in dealing with imprecise tomographic data, extracting Markovian dynamics within any desired regime of tolerance. If no Markovian channel is found, the scheme provides a well-defined quantitative measure of non-Markovianity in terms of the minimal addition of white noise required to "wash-out" memory effects, and render the evolution Markovian.
A significant part of the work is focused on the treatment of input matrices that are perturbations of some unknown process with a degenerate spectrum. This situation commonly arises when analysing noisy unitary gates in quantum computation. In order to address the susceptibility of the convex optimisation programme with respect to perturbation of multidimensional eigenspaces, we have designed a series of pre-processing algorithms, rigorously rooted in the theory of matrix perturbation, which identify and re-contruct the unperturbed hermiticity-preserving structure of the original channel. To benchmark our theoretical formulation, we have implemented and numerically benchmarked our algorithm on simulated data from the Cirq platform, demonstrating that our algorithms are able both to successfully identify channels consistent with an underlying Markovian dynamics, as well detect and quantitify non-Markovianity.
One drawback of the algorithms developed here is that they require a full tomographic snapshot, thus become infeasible for large numbers of qubits. Extending these algorithms to the case where we only have access to measurement data that gives incomplete information about the full dynamics, or incorporating partial prior information about the dynamics, is feasible and will be explored in future work. Finally, as stated in the introduction, another key avenue for future work is to extend the analysis and methods of this paper to the case of time-dependent noise models, allowing the techniques to be applied to quantum dynamics at longer timescales over which the noise processes are likely to vary.

Code availability
The python code implementing the algorithms presented in this work and used to produce the numerical results is available at [46]. The 1-qubit X-gate analysis (with 10,000 random samples) took 2 hours on a standard Intel x86 2.0GHz laptop, and the 2-qubit ISWAP analysis (8,500 random samples) took 2 weeks on a standard Intel x86 3.40Ghz desktop machine. We made no effort to optimise the algorithm implementation. In particular, the most costly part of the algorithms, namely the random sampling in the pre-processing, is trivially parallelisable, but we didn't do this here. The run-time can certainly be reduced significantly if desired.