Purifications of multipartite states: limitations and constructive methods

We analyze the description of quantum many-body mixed states using matrix product states and operators. We consider two such descriptions: (i) as a matrix product density operator of bond dimension D, and (ii) as a purification that is written as a matrix product state of bond dimension D'. We show that these descriptions are inequivalent in the sense that D' cannot be upper bounded by D only. Then we provide two constructive methods to obtain (ii) out of (i). The sum of squares (sos) polynomial method scales exponentially in the number of different eigenvalues, and its approximate version is formulated as a Semidefinite Program, which gives efficient approximate purifications whose D' only depends on D. The eigenbasis method scales quadratically in the number of eigenvalues, and its approximate version is very efficient for rapidly decaying distributions of eigenvalues. Our results imply that a description of mixed states which is both efficient and locally positive semidefinite does not exist, but that good approximations do.


INTRODUCTION
Quantum many-body systems appear in a variety of fields in physics, such as condensed matter, quantum chemistry, or high-energy physics. Since their Hilbert space description is intractable (as it scales exponentially with the system size), a number of methods have been proposed to describe them efficiently. One of them are tensor networks [1], which have been particularly successful in describing one-dimensional pure states with Matrix Product States (MPS) [2]. Its canonical form has facilitated the distinction between injective and non-injective MPS [3], (which determines the ground state degeneracy of the parent Hamiltonian, and is linked to many other physical properties, see e.g. [4][5][6]), which has led to the classification of gapped phases of one-dimensional systems [7]. This mathematical understanding has allowed to characterize global properties of the state, such as topological order or symmetries, in a local way [8,9]. This is in contrast with one-dimensional mixed states, whose description with tensor networks is much more scarce. While the class of Matrix Product Density Operators (MPDOs) has been defined [10,11], a canonical form has not been found. This is despite the fact that mixed states are used to describe quantum many-body systems at finite temperature, systems out of equilibrium, dissipative dynamics, or lack of knowledge of the state of the system.
One of the difficulties in defining a canonical form for mixed states is that with MPDOs one cannot verify locally that the global tensor is positive semidefinite. This implies that local truncations of the tensors generally destroy the global positivity, causing instability in numerical algorithms. An alternative is to use the local purification form, where the mixed state is purified and the purification is expressed as an MPS. While its local tensors are positive semidefinite, it is not known if an efficient MPDO form guarantees an efficient local purifi-cation. That is, are the two descriptions equivalent, or can the local purification be much more costly than the MPDO form (see figure 1)?
... In this paper we address this question and show that these two descriptions are inequivalent. Specifically, we provide classical multipartite states whose MPDO form has a fixed cost, but whose local purification form has an unbounded cost. Then, we provide two constructive purification methods applicable to all multipartite density matrices, which relate the two forms but also involve the number of (different) eigenvalues. The sum of squares (sos) polynomial method has an exact version which scales exponentially in the number of different eigenvalues. Its approximate version is formulated as a Semidefinite Program (SDP), which shows an efficient and robust behavior for all the tested distributions of eigenvalues. The eigenbasis method has an exact version which scales multiplicatively with the number of eigenvalues, and its approximate version gives very efficient purifications for rapidly decaying distributions of eigenvalues.
This paper is organized as follows. First, we present the problem in Sec. 2. Then we show the inequivalence of the MPDO form and the local purification form in Sec. 3. In Sec. 4 we present the two purification methods: the sos polynomial method (Sec. 4.1), and the eigenbasis method (Sec. 4.2), both with its main idea, its exact and its approximate version, and in Sec. 4.3 we compare the two approximate methods. Finally we conclude and mention further directions in Sec. 5.

THE SETTING
We now present the question that concerns us, which is whether the MPDO form and the local purification form of a mixed state can be related. We will first present some notation and definitions (Sec. 2.1), and then introduce the problem (Sec. 2.2).

Definitions
Let us first fix some notation. We write A 0 to denote that the matrix A is positive semidefinite (i.e. Hermitian with nonnegative eigenvalues), and A ≥ 0 to denote that it is nonnegative (i.e. with nonnegative entries). We also write rank(ρ) (rank(A)) to denote the number of non-zero eigenvalues (non-zero singular values) of the matrix ρ (A). The trace norm of a matrix A is defined as ||A|| 1 = r i=1 s i , where r = rank(A) and s i are its singular values. Given a multipartite state with N subsystems arranged on a one-dimensional line (henceforth called a 1D state), we call a linear bipartition a splitting of the form 1 . . . k vs. k + 1, . . . , N for any k. Finally, given a multipartite state |ψ = ψ i1...iN |i 1 . . . i N with Schmidt rank D ′ k across the bipartition i 1 . . . i k vs. the rest, we define the Schmidt rank of |ψ as SR(ψ) = max k D ′ k . We focus on 1D mixed states with N d-level systems and open boundary conditions, described by a density matrix ρ, The are two natural ways of describing ρ locally. The first one is MPDO form [10,11], which is defined as where M ) is a row (column) vector of size d. Note that the subindex only indicates the subsystem that the matrix is describing. Here D k (for all k) is the minimal dimension such that (2) holds [26]. The operator Schmidt rank of ρ is defined as [27] OSR(ρ) := max The second form is the local purification form, which is obtained by purifying the mixed state ρ (living in system S) into a pure state |Ψ (living in S, S ′ ), and expressing |Ψ as an MPS [28], Here D ′ k is the Schmidt rank of |Ψ of the bipartition Tr k+1...N |Ψ Ψ| vs. the rest. We define the purification rank of ρ as

The problem
We want to find out if the MPDO and the local purification form are equivalent, or if the latter can be arbitrarily more costly than the former. The advantage of the local purification form is that the local tensors are positive semidefinite, since they are of the form A k A † k . In contrast, in the MPDO form this is generally not true, i.e. M α k−1 ,α k k 0, thus it is locally invisible that ρ 0. This is a problem theoretically as well as numerically, as local truncations destroy the global positivity.
On the other hand, the problem of the local purification form is that it is not known how much larger the purification rank D ′ may be compared to the operator Schmidt rank D (see figure 2). Thus, we want to see if there is a transformation from the MPDO to the local purification form in which the bond dimension increases in a controlled way. That is, we want to find out if D ′ can be upper bounded by a function of D.
The Hilbert space H of mixed states. How much larger is the family of states with fixed operator Schmidt rank D compared to that with fixed purification rank D ′ ? Are they comparable in size, or are there states requiring very small D and very large D ′ (such as the one marked with a cross)?
Observe that it is very easy to obtain a bound in the other direction, namely D ≤ D ′2 , simply by choosing where i k , j k , z k are physical indices. On the other hand, both ranks can be upper bounded by the physical dimension,

INEQUIVALENCE OF THE TWO FORMS
We now show that the MPDO form and the local purification form are inequivalent. More precisely, we provide a family of multipartite classical states with a constant operator Schmidt rank across every linear bipartition, and an unbounded purification rank (Result 2). We will first show the separation for bipartite states (Sec. 3.1), and then for multipartite states (Sec. 3.2).

Inequivalence for bipartite classical states
We now show the inequivalence of the two forms for bipartite classical states, which follows from recent results in convex analysis shown by Gouveia, Parrilo and Thomas [12]. Result 1 Let ρ denote a bipartite classical mixed state. Let D = OSR(ρ) and D ′ = rank puri (ρ). Then D ′ cannot be upper bounded by a function f that only depends on D (in particular, f does not depend on ρ), i.e.
Before proving the result, we establish a relation between the MPDO and the local purification form of classical states and two other known decompositions. Consider the family of classical bipartite states (i.e. diagonal in the computational basis {|x }) of physical dimension t where S t is a nonnegative matrix of size t × t and S t (x, y) denotes its x, y component. Now, the MPDO form of ρ t corresponds to the singular value decomposition of S t . Hence, the operator Schmidt rank of ρ t corresponds to the rank of S t , rank(S t ). On the other hand, the local purification form of ρ t corresponds to the positive semidefinite factorization of S t , in which S t is expressed as for all x, y (the local purification form is obtained by expressing E x = A x A † x and F y = B y B † y ). Thus, the purification rank of ρ t corresponds to the positive semidefinite rank of S t , rank psd (S t ), which is defined as the minimal r such that there exist matrices E x , F y 0 of size r×r such that (11) holds (see figure 3). The positive semidefinite factorization was very recently introduced in Ref. [13] and shown to be related to the quantum communication complexity [13] and the quantum correlation complexity of S t [14]. In summary, for classical states such as ρ t , it holds that [29] OSR(ρ t ) = rank(S t ) rank puri (ρ t ) = rank psd (S t ) . Now we are ready to prove Result 1. Proof We consider classical bipartite states of the form (10), and we focus on a class of nonnegative matrices S t called slack matrices of polytopes, defined as follows (see e.g. [12]). A convex polytope is defined as the intersection of a finite set of halfspaces {h j (x) ≤ b j }, or as the convex hull of a set of vertices {v i } [15]. Its slack matrix S is defined so that its (i, j) entry contains the distance from hyperplane j to vertex i, i.e. S(i, j) = b j − h j (v i ) (see figure 4). Now we let S t be the slack matrix of the twodimensional regular polytope with t vertices (and thus also t faces), called the regular t-gon [30]. Gouveia, Parrilo and Thomas show that [12] [31] Using the equivalences (12), this implies that That is, the operator Schmidt rank of ρ t is constant for all t, whereas its purification rank grows unboundedly with t. It follows that there does not exist an upper bound of rank puri (ρ t ) that depends only on OSR(ρ t ). The slack matrix of the regular hexagon S6 contains in its (i, j) entry the distance from vertex i (here labeled with lowercase letters) to hyperplane j (labeled with uppercase letters). Note that the matrix is circulant.

Inequivalence for multipartite classical states
Now we show a more general form of separation, namely we provide classical multipartite states with a constant operator Schmidt rank across every linear bipartition, and an unbounded purification rank. Result 2 Let ρ denote a mutipartite classical mixed state. Let D = OSR(ρ) and D ′ = rank puri (ρ). Then D ′ cannot be upper bounded by a function f that depends Proof Consider the family of states ρ t (equation (10)) with t = 2 m (with natural m). That is, S t is the slack matrix of the square, octagon, etc. Let the binary representations of the row index x and column index y be (x 1 , . . . , x m ) and (y 1 , . . . , y m ), respectively. Since the polygon is regular, S t is a circulant matrix. Hence, it is diagonalized by the Fourier transform F t , with components F t (x, y) = ω xy , where ω = exp(i2π/2 m ). That is, where and {λ α } are the eigenvalues of S t . Now we decompose the tensors L and R into smaller tensors, where each M k depends only on one bit (x k or y k ), That is, each M k is a 3 × 3 diagonal matrix, except for M 1 and M 2m , which are a row and a column vector, respectively, defined analogously. This shows that ρ t has operator Schmidt rank 3 across every linear bipartition.
The classical state ρt (equation (10)), where St is the slack matrix of the regular t-gon, shows the separation of the MPDO form and the local purification form. (Middle) In the bipartite case, ρt has operator Schmidt rank 3 for all t, and purification rank ∼ log t. (Right) In the multipartite case, for t = 2 m , ρt has operator Schmidt rank 3 across every linear bipartition, and purification rank ∼ log t at least across one bipartition.
We know from equation (14) that the purification rank of ρ t along the x vs. y bipartition grows unboundedly like log t. Clearly, a small purification rank across any bipartition would imply a small purification rank across the x vs. y bipartition. Thus, the purification rank is unbounded at least across one bipartition (see figure 5).
This shows that a small operator Schmidt rank across all linear bipartitions does not imply a small purification rank.
Let us make a final remark. Note that a possible purification of ρ t is where the second and fourth index refer to the ancillary states. Thus, rank puri (ρ t ) ≤ SR(ϕ t ). Above we have seen that rank puri (ρ t ) ∼ log t, thus the Schmidt rank of |ϕ t grows with t, and so does its preparation cost [32]. Now consider the state obtained by taking the square of the coefficients of |ϕ t , From Result 2 it follows that the Schmidt rank of this state across any linear bipartition is 3. Hence its preparation cost is constant with t [16]. Thus, we see that transforming (21) to (20), i.e. taking the Hadamard square root of the coefficients, may have a high cost in the Schmidt rank, and thus in the preparation cost of a pure state.

PURIFICATION METHODS
We will now present two constructive purification methods: the sos polynomial method (Sec. 4.1) and the eigenbasis method (Sec. 4.2). Both are applicable to all multipartite mixed states, and can be used to construct exact and approximate purifications. We will compare both approximation methods and see that they are complementary for various eigenvalue distributions in Sec. 4.3.

Sos polynomial method
We will first present the idea of the sos polynomial method (Sec. 4.1.1), and then explain how to use it to construct exact (Sec. 4.1.2) and approximate purifications (Sec. 4.1.3).

The idea
The idea of the sos polynomial method is the following: given a mixed state ρ, we construct a purifying state as a sum of powers of ρ (up to certain degree), where each power is attached to an ancillary state. If the degree is large enough, there exists a choice of the ancillary states such that this purifying state is an exact purification for ρ. If the degree is not large enough, one can find the ancillary states with an ansatz of sos polynomials or with an SDP.
Specifically, consider a multipartite density matrix ρ of size d N × d N , with spectral decomposition where λ i = 0 for i > n. We construct a purifying state |Ψ k as a sum of powers of ρ, from 0 to k − 1, where each power is attached to an ancillary state (see figure 6), where |ρ denotes a vectorized matrix ρ (i.e. given ρ = . We use |Ψ k as the purifying state of a density matrix σ k , Here p k is a polynomial that can be written as where R k is a positive semidefinite matrix which is the Gram matrix of the ancillary states, i.e. its (i, j) component is R k (i, j) = a i |a j . This polynomial is sum of squares (sos), as can be readily seen by writing R k = A T A, where A u denotes the uth row of A, and r = rank(R k ). Note that, since the polynomial is univariate, the set of sos polynomials is identical to the set of nonnegative polynomials (i.e. with the property p(x) ≥ 0 for all x) [17]. Now, by construction, the purification rank of σ k is at most the Schmidt rank of |Ψ k (equation (24)). Denoting OSR(ρ) = D, this is at most 1 + D + D 2 + . . . + D k−1 (equation (23)), and thus We are interested in making σ k as close as possible to ρ. We will show that if k = m, where m is the number of different nonnegative eigenvalues of ρ [33], there exists ancillary states such that p k (λ i ) = λ i for all i, and thus σ m = ρ, i.e. σ m is an exact purification of ρ (Sec. 4.1.2). If, on the contrary, k < m, one can choose k − 1 points and construct the p k that passes through them, or one can find the p k that minimizes the distance between σ k and ρ (in trace norm) with an SDP (Sec

Exact case
We now show how to build exact purifications with the sos polynomial method. Result 3 Let ρ denote a multipartite density matrix with m different nonnegative eigenvalues. Let D = OSR(ρ) and D ′ = rank puri (ρ). The sos polynomial method provides an exact purification of ρ with Proof To ease the notation we consider that the different nonnegative eigenvalues are the first m (i.e. λ i = λ j for i = j and i, j = 1, . . . , m). Consider the construction of the purifying state σ k as explained in Sec. 4.1.1. Define the vector |v i k as so that the sos polynomial evaluated at λ i can be written as p k (λ i ) = v i k |R k |v i k (compare with equation (25)). Now we choose k = m. The important observation is that the set of vectors {|v i m } m i=1 is linearly independent (since only different eigenvalues are considered). Hence, there exists another set {|w j m } m j=1 which is biorthogonal to it, i.e. v i m |w j m = δ ij for all i, j. Then we choose R m as follows: This satisfies that p m (λ i ) = v i m |R m |v i m = λ i for all i. From equation (24) it follows that σ m = ρ and thus σ m is an exact purification of ρ. Finally, using equation (28) with k = m, the claim of the result follows.
Note that Result 3 depends on the number of different nonnegative eigenvalues m because σ k = ρ requires that p(λ i ) = λ i for all i, and these are only m independent conditions.

Approximate case: SDP
Building approximate purifications with the sos polynomial method can be done analytically or numerically, as we show next. Result 4 Let ρ denote a multipartite density matrix of dimension d N × d N , with rank(ρ) = n. Let D = OSR(ρ). The sos polynomial method provides a density matrix σ k with and where p k is a sos polynomial of degree 2(k − 1). p k can be constructed by choosing k − 1 nonnegative points through which it must pass. Alternatively, the optimal p k can be found with an SDP.
Proof Sos polynomial that passes through k − 1 points. We construct the sos polynomial p k by letting it pass through k − 1 chosen points {µ 1 , . . . , µ k−1 }. We use the Lagrange basis "squared" to this end: where we have omitted the dependence of l j on k. This satisfies p k (µ j ) = µ j for j = 1, . . . , k − 1. Note that the degree of p k is 2(k − 1) as required. The distance (33) depends on the points {µ i }, thus the difficulty lies in choosing them.
The optimization problem as an SDP.
We search for the sos polynomial p k that minimizes the distance (33), i.e. for the positive semidefinite matrix R k that minimizes it (see equation (25)). The optimization problem thus reads The objective function can be made linear by introducing the slack variables z, The optimization variables are now R k 0 and z ≥ 0, and the constraints are linear in them. Thus, this is an SDP optimization problem (see Appendix A for the precise formulation).
In words, the SDP searches for the sos polynomial p k (of degree 2(k − 1)) whose distance to the eigenvalue distribution is minimal in trace norm. This formulation is consistent with the fact that optimization over sos polynomials can be done with SDPs [17].
Note that the non-trivial condition is that the polynomial be sos (equivalently, nonnegative), since otherwise one could take p(λ) = λ. Approximating λ for 0 ≤ λ ≤ λ max and another nonnegative function elsewhere has the problem that the function at λ = 0 is non analytical.
We remark that there may exist exact solutions for m/2 < k < m, since p k (λ) − λ has degree 2(k − 1), and thus can have 2(k − 1) real roots. However, we only know how to construct the exact solution for the case k = m (with the idea of the proof of Result 3).
Remark 5 Non-orthogonal ancillary states. Observe that the exact solution for k = m involves nonorthgonal ancillary states, since R m is non-diagonal (equation (31)). This is so because the basis {|v i m } m i=1 is not orthogonal, and neither is the biorthogonal basis {|ω j m } m j=1 . More generally, the solution for k < m also involves non-orthogonal ancillary states, since orthogonal states result in a diagonal R k , which renders a polynomial only even powers with non-negative coefficients (i.e. a monotonously increasing polynomial for positive λ). In contrast, a non-diagonal R yields polynomials with possibly odd powers with negative coefficients, thus with various minima, rendering a better approximation of λ in the desired interval. Remark 6 Real ancillary vectors. The ancillary vectors {|a j } can be taken real without loss of generality. To see this, first write the polynomial (equation (25)) as where c l = s+t=l R k (s, t). That is, each coefficient c l is the sum of the lth antidiagonal of R k . Thus c l depends on diagonal elements of R k , which are nonnegative, or on the sum of an element and its transposed (R k (s, t) + R k (t, s)), which is real. Therefore c l only depends on the real part of R. If R 0, then Re(R) is positive semidefinite for all real vectors |r , since r|Re(R) + iIm(R)|r ≥ 0, which implies i r|Im(R)|r = 0. Since we only consider contractions with real vectors (namely |v i k ), we can restrict R to its real part. Thus, its spectral decomposition reads R = ODO T , where O is orthogonal, and D is diagonal and nonnegative. This readily yields R = O √ D √ DO T = AA T , where the ith row of A contains the coefficients of the ancillary state a i | expressed in the eigenbasis of R, which are real.
While the exact sos polynomial method depends on m (Result 3) and the approximate on d N (Result 4), we will see in Sec. 4.3 that in practice there exists good ansätze of sos polynomials which make it independent of both. We will also compare this approximate method with the eigenbasis method, which we present next.

Eigenbasis method
Now we turn to the eigenbasis method, for which we present the main idea (Sec. 4

The idea
The idea of the eigenbasis method is to consider the standard purification of ρ, which is the spectral decomposition, and upper bound the Schmidt rank of each eigenstate of ρ. This is done by constructing a basis of the range of ρ where each basis element is the image (under the map ρ) of a product state. Thus, each has Schmidt rank at most D, where D is the operator Schmidt rank of ρ. Expressing each eigenstate in terms of this basis, we see that each can have Schmidt rank at most Dn, where n is the number of non-zero eigenvalues of ρ. Thus the purification rank of ρ is at most Dn 2 .
To be more precise, consider the standard purification of ρ, obtained from its spectral decomposition (see equation (22)) as By construction, the purification rank of ρ is at most the Schmidt rank of |Ψ , which is at most n times the maximum Schmidt rank of |φ i , Our goal is to upper bound max i SR(φ i ) as a function of D. To this end, we build a basis of the range of ρ where each basis element |χ α is the image (under the map ρ) of a certain product state |p α (see figure 7), where α = (α 1 , . . . , α N ), α i is a label of the ith product state, and the subindex outside the ket denotes the subsystem that the state is describing. Now, consider the MPDO form of ρ as in equation (2), with operator Schmidt rank D. The Schmidt rank of |p α is one, and applying ρ to |p α increases the Schmidt rank by at most D. It follows that On the other hand, we consider the spectral decomposition of ρ (equation (22)), and use the eigenstates as a basis of the range. In particular, we express |χ α in terms of this basis, with coefficients f iα of the linear combination, (45) Since the range has dimension n, there are at most n linearly independent |χ α .

Exact case
We now show how to use the eigenbasis method to construct an exact purification. Result 7 Let ρ denote a multipartite density matrix with rank(ρ) = n. Let D = OSR(ρ) and D ′ = rank puri (ρ). Then the eigenbasis method constructs an exact purification of ρ with Proof Let ρ be a multipartite density matrix with spectral decomposition as in (22) (thus with rank(ρ) = n). Building upon Sec. 4.2.1, we only need to show that there exist n product states {|p α } n α=1 such that their images under ρ, {|χ α = ρ|p α } n α=1 , form a basis of the range of ρ. Take a product basis {|x } with x = 1, . . . , d N , which spans the whole space. Then {ρ|x } spans the range of ρ. Then, we can select a basis ρ|x s of the range of ρ with s = 1, . . . , n. We call the latter product states |p α with α = 1, . . . , n.
This implies that the coefficient matrix (f iα ) (equation (45)) is full-rank and hence can be inverted. We denote the elements of the inverse by g αj , i.e. n α=1 f iα g αj = δ ij , and we invert equation (45) to express the eigenvector |φ i as a linear combination of the vectors |χ α , Since SR(χ α ) ≤ D, we obtain for all i. Finally, using (41), the claim of the result follows.

Approximate case
We now show how to use the eigenbasis method to build approximate purifications. Result 8 Let ρ denote a multipartite density matrix with spectral decomposition as in (22), and let D = OSR(ρ). The eigenbasis method provides a density matrix σ s with rank(σ s ) = s such that and its distance to ρ is Proof We construct σ s with the largest s eigenvalues of ρ, i.e.
where N = s i=1 λ i . By direct calculation one can see that (50) holds with equality. Applying Result 7 to σ s yields equation (49).
Clearly, this method yields good approximations for rapidly decaying distributions of eigenvalues, for which the distance (50) is small. In the next section we make these statements precise, and compare this method to the sos polynomial method.

Comparison of approximation methods
We now compare the two approximation methods for a state ρ with rank(ρ) = n with the following eigenvalue distributions (where the eigenvalues are ordered in nonincreasing magnitude).
1. Uniform distribution, defined as λ j = 1/n for all j.

Sos polynomial method
We start with the sos polynomial method, for which we first present analytical and then numerical results obtained with the SDP, both based on Result 4.
First, the uniform distribution is the easiest for this method. Using the exact result (Result 3) with m = 1 we obtain D ′ = 1, since this distribution describes the maximally mixed state. The sos polynomial that constructs this exact purification is p 1 (λ) = 1/n. Second, the sos polynomial that best approximates the line λ in the interval [0, λ 1 ] is a good ansatz for the equally spaced and random distribution. Moreover, these polynomials can be chosen so that their distance to the distribution is independent of n, as we show next. Result 9 Consider a distribution of n eigenvalues whose largest eigenvalue λ 1 ∼ 1/n, such as the equally spaced or the random distribution. Then there exists a sos polynomial p k such that Proof Let p k denote the sos polynomial of degree 2(k−1) that best approximates the straight line λ in the interval [0, 1]. This can be rescaled to approximate the straight line in [0, λ 1 ], Now, let ε = max λ∈[0,1] |p k (λ) − λ|.
If λ 1 ∼ 1/n, this upper bound is independent of n.
For the distribution with one fixed eigenvalue and the rest equally spaced, it holds that λ 2 ∼ 1/n, hence the sos polynomial q k (λ) = p k−1 (λ)(λ − λ 1 ) 2 , where p k−1 is defined in the proof of Result 9, also has a distance independent of n. Denoting the distance ||ρ − σ k || 1 = ǫ = g(k), and using Result 4, this implies that one can upper bound the purification rank of σ of these three distributions by a function independent of n, namely rank puri (σ) ≤ O(D g −1 (ǫ)−1 ). Finally, for the exponentially decaying distribution, we present an ansatz of sos polynomials whose distance decreases exponentially with k and is independent of n. Result 10 Consider the exponentially decaying distribution of eigenvalues λ j = ae −bj (with a, b > 0) for j = 1, . . . , n, and the sos polynomial that passes through the largest k − 2 eigenvalues and 0, constructed as (57) The distance from this polynomial to the distribution decays exponentially in k, i.e.
Proof See Appendix B. The exponential decrease in the distance implies that the purification rank of σ scales polynomially in D, For the exponentially decaying distribution, we provide A and show that B = b in Appendix B.
Let us now analyze how well these eigenvalue distributions are approximated with the SDP. For distributions 2 to 5, the SDP gives an exponential decrease of the distance in k and independent of n (equation (59)) with A ≈ 4 and B ≈ 2 for distributions 2, 3, A ≈ 3 and B ≈ 1.3 for 4, and for 5 with b = 1 we find A ≈ 4 and B ≈ 1.3 (see Figs. 8,9). The sos polynomials found by the SDP for the equally spaced distribution are shown in figure 10. We remark that the numerics does not improve beyond k ≈ 4, as very small numbers (such as powers of small eigenvalues) are numerically treated as 0. We have rescaled the eigenvalues λ i → λ i /λ 1 for all i, which allows us to use the SDP until k ≈ 7 for some distributions.

Eigenbasis method
We now use the eigenbasis method (Result 8) to upper bound the purification rank of σ (equation (53)) for the eigenvalue distributions 1 to 5.
1. Uniform distribution. This is the hardest distribution, since the smallest eigenvalues are as large as possible. We obtain 2. Equally spaced distribution. We obtain 3. Random distribution. The distance depends on the particular random distribution; in the worst case, it is the uniform distribution, hence it is upper bounded by equation (61).
4. One fixed eigenvalue and the rest equally spaced. We obtain 5. Exponentially decaying distribution. Assuming that e −n ≪ ǫ, the purification rank of σ grows lin- early in D and quadratically in ln(1/ǫ), i.e.
In summary, for the uniform distribution, the sos polynomial is the best, as it shows exactly that D ′ = 1. For the equally spaced, random, and one large eigenvalue and the rest equally spaced distributions, the sos polynomial method is better as it is independent of n and scales polynomially in D. Finally, for the exponentially decaying distribution, the eigenbasis method is better, as it scales linearly in D (and quadratically in ln(1/ǫ)). We thus see that the sos polynomial method is very robust, as it yields the same scaling for very different eigenvalue distributions, while the eigenbasis method works well only for rapidly decaying distributions of eigenvalues, in which case it works better than the other.

CONCLUSIONS AND OUTLOOK
In this paper we have analyzed the efficiency of representing a mixed state as an MPDO and as a local purification, and we have shown that the latter can be arbitrarily more costly than the former. In particular, we have provided a family of multipartite classical states with a constant operator Schmidt rank D across each linear bipartition and an unbounded purification rank D ′ (Result 2). This shows that, in the exact case, one cannot upper bound D ′ by a function of D only.
Then we have presented two constructive purifications methods which are applicable to any multipartite density matrix. The exact sos polynomial method implies that D ′ ≤ O(D m−1 ), where m is the number of different eigenvalues. Its approximate version consists of finding the sos polynomial which passes through certain points, and the optimal one can be found with an SDP (Result 4). For the four tested eigenvalue distributions, this method upper bounds D ′ by a polynomial function of D which is independent of n, thus showing a robust and efficient behavior.
The exact eigenbasis method implies that D ′ ≤ Dn 2 , where n is the number of eigenvalues. Its approximate version discards the smallest eigenvalues (Result 8), and for the exponentially decaying distribution it yields D ′ scaling linearly in D (and independent of n).
Our inequivalence result (Result 1) implies that a single canonical form which is both efficient and locally positive semidefinite cannot exist in the exact case. Note that this is also a numerical limitation, as contracting two-dimensional (2D) PEPS requires determining the 1D density matrices which are at the boundary [18]. In order to have an efficient algorithm for 2D PEPS it is thus necessary to use an efficient description of such 1D mixed states, hence to work with the MPDO form. On the other hand, for the algorithm to be stable, it is necessary that the positivity of the operator can be checked efficiently, which means for example, that it can be checked locally. While in practice one can work with MPDOs of low dimension, one would hope to reexpress such MPDOs in terms of a purification fulfilling the latter requirement, but without increasing significantly the bond dimension.
Our results show that this may not be possible.
At the same time, the results of section 4.3 show that one can construct efficient approximate purifications of various relevant eigenvalue distributions. One should nonetheless analyze how successive truncations of this approximate purification affect the total error.
Let us now mention some open questions. Concerning the inequivalence of the two forms, we believe that a larger separation could be obtained with states of the form ρ = I ⊗ I − P 1 ⊗ Q 1 − P 2 ⊗ Q 2 where P i , Q i are Hermitian operators constrained by the fact ρ 0, but otherwise with random entries. Concerning the sos polynomial method, it would be interesting to find sos polynomials whose distance to the equally spaced distribution decreases exponentially with the degree k. It would also be appealing to combine both purification methods in a single one with the best of each, but this requires to split the density matrix into different 'eigenspace sectors' (such as in ρ = ρ 1 + ρ 2 where ρ 1,2 = i≤r,i>r λ i |φ i φ i |), and it is not clear how the operator Schmidt rank of ρ relates to that of ρ 1 or ρ 2 . On a more practical level, it would be worth devising a purification method that works sequentially and does not require to know the spectral decomposition of the density matrix.
Our results also have connections to other research areas. In particular, the results on decomposability of mixed states translate one-to-one to divisibility properties of completely positive (CP) maps via the Choi-Jamio lkowski isomorphism. While the divisibility of CP maps with a fixed dimension of the intermediate space has been studied e.g. in [19,20], our approach would allow to extend this study to varying middle dimension. Finally, our investigations are also related to communication complexity, since the positive semidefinite rank determines the quantum communication complexity [13] and the quantum correlation complexity [14] of the associated matrix. depend on the order on which the decompositions are made. To see this, note that D k is the bond dimension of the bipartition i1j1 . . . i k j k vs. the rest of the pure state |ψ = ρ j 1 ,...,j N i 1 ,...,i N |i1, . . . , iN |j1, . . . , jN .
[27] The operator Schmidt rank of ρ must not be confused with its tensor rank, which is the minimal r such that ρ = r α=1 M α 1 ⊗ M α 2 ⊗ . . . ⊗ M α N .
[28] There are infinitely many purifications of a given mixed state, and they are all related by an isometry on the an-cillary system. We always consider the purification which minimizes the Schmidt rank of the purifying state along every linear bipartition and which is local (i.e. an ancillary subsystem S ′ i is attached to each subsystem Si).
[29] The nonnegative factorization [21] of nonnegative matrices such as St has also been defined, and has traditionally received more attention, in particular in its connections to classical communication complexity. In this case, St is expressed as a product of two nonnegative matrices X and Y . The nonnegative rank of St, rank+(St), is the minimal r such that there exist X, Y ≥ 0 of size t × r and r × t, respectively, such that St = XY . Note that this factorization cannot be defined for general quantum states.
[32] In the sense of the dimension of the ancilla required to prepare the state in a sequential scheme [16].
[33] By nonnegative eigenvalues, we mean that the 0 should also be counted as an eigenvalue. That is, if ρ is rank deficient, m is the number of different non-zero eigenvalues plus one.
where I d N is the identity matrix of size d N and 0 k the zero matrix of size k × k. The matrix constraints are given by A j = diag(0, . . . , −1, . . . , 0) ⊕ (−|v j k v j k |) 1 ≤ j ≤ d N (A5) where the −1 is in the jth position. Thus we have s = 2d N constraints. Note that one can always write the constraints of (A2) with equalities by introducing additional slack variables [23]. Finally, b = (−λ 1 , . . . , −λ d N , λ 1 , . . . , λ d N ) T .
We remark that the SDP is strongly feasible for k < m. To see that, note that the vectors {|v i k } k i=1 are linearly independent, hence there exists a biorthogonal basis {|w i k } k i=1 . We use them to define R = k i=1 |w i k w i k |. This R is positive definite (R ≻ 0) and we choose z to be strictly positive (z > 0). This point is in the interior of the feasible region. It follows that there exists no duality gap in the SDP [23].
We have implemented the SDP with SeDuMi 1.3 [24] and the add-on Yalmip [25], with default options.