A Practical Approach to SOS Relaxations for Detecting Quantum Entanglement

This article introduces PnCP, a MATLAB toolbox for constructing positive maps which are not completely positive. We survey optimization and sum of squares relaxation techniques to find the most numerically efficient methods and establish some benchmarks for this construction. We also show how this package can be applied to the problem of classifying entanglement in quantum states.


Introduction
For n ∈ N, let M n (R) be the vector space of n × n matrices over R and let S n be the subspace of symmetric matrices {A ∈ M n (R) : A T = A}. A matrix A ∈ S n is called positive semidefinite if all eigenvalues of A are nonnegative; in this case, we write A 0. Given two matrix spaces A and B, a linear map Φ : A → B with Φ(A T ) = Φ(A) T for all A ∈ A, is called positive if for all (symmetric) A 0, one has Φ(A) 0. Given k ∈ N, such linear maps induce the ampliation where ⊗ is the standard Kronecker tensor product of matrices. If Φ (k) is positive, 1 then we call Φ k-positive. If Φ (k) is positive for all k ∈ N, then Φ is called completely positive. Positive and completely positive maps arise naturally in matrix theory and operator algebras (e.g., positive linear functionals) [34,40], frequently in quantum information theory [18,32,38], and have recently even been used in semi-definite programming [21].
Our interest in these maps comes from quantum information, where positive maps that are not completely positive (or pncp maps for short) can detect quantum entanglement. Quantum entanglement is an interesting feature arising from quantum mechanics, which allows physical systems to be correlated non-locally. There are many applications of quantum entanglement in the field of quantum information theory, such as quantum teleportation [6], quantum computation [20] and quantum cryptography [12]. For any such application, it is essential to first identify quantum states which are entangled. The separability problem asks precisely this question; given a state (density matrix) ρ in a composite system, determine whether it is entangled. This classification problem can be reformulated in terms of pncp maps, which preserve their positivity on separable states, and however, they may fail to preserve positivity on entangled states, providing the following classification criterion.

Criterion 1.1 (The general criteria, [2] section 8.4) A quantum state ρ ∈ M t (R) is entangled if there is a pncp map Φ such that the ampliation (I ⊗ Φ)(ρ) 0, where I is the identity matrix.
As an example, consider the Bell State, which has density matrix (see Sect. 5) which has a negative eigenvalue of −1/2 and serves as evidence of entanglement in the Bell State. While the transposition map was sufficient in this simple example, in general finding a suitable map is difficult.
In their recent work [22], the authors presented an algorithmic construction of pncp maps via their correspondence to positive and nonnegative polynomials. Restricting these involution-preserving maps to the space of symmetric matrices, each Φ : S n → S m gives rise to a biquadratic, bihomogeneous polynomial p Φ ∈ R[x, y], with It is a classical result [9] that Φ is positive if and only if p Φ is nonnegative on R n+m , and Φ is completely positive if and only if p Φ is a sum of squares (SOS) on R n+m (a nice proof is given in [22]).
The algorithm of [22] (which is a specialization of the work of [4]) builds generic biquadratic, bihomogeneous polynomials over the Segre variety which are positive, but not SOS, i.e., pncp maps. Letting n, m > 2, d = n + m − 2 and N = n + m, the algorithm of [22] ( [4]) can be summarized as: The implementation of Steps 1-2 is done through computational linear algebra. While there are of course subtleties and complication in Steps 1-2 as well, these can be overcome in principle by using randomly generated data. For a more detailed discussion of this, we refer the interested reader to the Appendix, where we present Algorithm 1 in full detail.
Step 3 on the other hand is a difficult problem, and SOS relaxations are now a common solution technique. In fact, SOS relaxations have even been previously employed in the context of entanglement detection [31].
Our contribution in this work is to find the most practical SOS relaxation for Step 3 and to establish benchmarks for this type of construction. In doing so, we also introduce the MATLAB package PnCP, currently the only implementation of Algorithm 1. We survey recent optimization techniques for verifying Step 3 and specify relaxations theoretically superior to those presented in [22]. We implement and test these methods in PnCP. Our package and test data are made available at https://github.com/Abhishek-B/PnCP. We also consider rationalizations of the forms obtained with Algorithm 1 to obtain exact certificates of nonnegativity (PnCP is able to construct pncp maps with rational coefficients).
The article is organized as follows: Section 2 reviews some notation and algebraic geometry background for the optimization involved in Step 3. In Sect. 3, we present some of the relaxations we surveyed and thought to be promising for use in Step 3. We also present our implementation of these methods using MATLAB and show their performance via computational efficiency (w.r.t. time) and success rate. Section 4 details issues in generating pncp maps with rational coefficients using Algorithm 1. We also show the difference in computational requirements for constructing maps with floating point coefficients and those with rational coefficients. Section 5 explains how we use PnCP to identify entanglement in quantum states. We demonstrate this usefulness through illustrative examples.

Background
We write R[x] for the set of polynomials in x = (x 1 , . . . , x n ), R[x, y] for the set of polynomials in the variables x and y = (y 1 , . . . , y m ), and R[x] ≤k denotes the polynomials in x of degree at most k. We use the standard notions of ideal and varieties as in [10]. Recall that for an ideal I ⊆ R[x], its radical is the ideal and I is said to be radical if I = √ I . We denote the set of polynomials which are sums of squares as

Minimization
Consider a general constrained minimization problem Step 3 of Algorithm 1 involves solving a maximization problem similar to (1). Specifically, Step 3 requires a solution to the following problem max δ>0 δ, However, as is well known, testing the nonnegativity of a polynomial p is an NPhard problem [29,33]. Thanks to recent advances in computational mathematics, SOS relaxations, which are computationally tractable, have become an efficient approximation tool (the idea is to decompose p as p = i q 2 i modulo the constraints). For the generic minimization problem (1), the standard relaxation is given by: Letting p min and γ (k) be the solutions to (1) and (3), respectively, Lasserre [24] has shown that γ (k) → p min as k → ∞ under some natural conditions. The recommended SOS relaxation in [22] for F δ in Step 3 of Algorithm 1 is: These relaxation problems can be stated and solved as an appropriate optimization program (semi-definite, second-order cone, quadratically constrained, etc.), and the majority of solvers can handle the broad class of these problems.

Alternative Relaxations and Performance
We now present alternate SOS relaxations to solving problem (1). To test the success rate of each relaxation, we generate 50 random forms using Algorithm 1: Step 1-Step 2 (with the standard rand function from MATLAB). Then, we employ Step 3 with each relaxation, and note how many forms are identified as positive and not completely positive.

Rational Functions
Let us begin by considering Artin's solution to Hilbert's 17 th problem [5].
This result provides the most fundamental SOS relaxation. For Step 3, we search for the maximal δ * such that F δ * is a sum of rational squares; that is If for some δ, F δ is nonnegative, then by Theorem 3.1 the SOS decomposition in (5) always exists. Of course to solve the general problem (5) using semi-definite programming, we must first bound the degree of σ (x, y) = 0.
Note that (5) is a quadratically constrained optimization program (nonlinear in the decision variables, δ and the coefficients of σ ), which can be solved with solvers such as PENLAB [13], but our early tests indicated that this approach is not ideal. So we instead implement (5) with a "bisection" approach. This is already the suggested method in [22], which tries to solve (4), and increases if a solution is not found. While bisecting may seem like a simple idea, given some tolerance ε, bisection achieves an ε-optimal For the Hilbert method (5), let G be the Gram matrix of σ . We fix δ = 2 0 , k = 1, and solve the following where tr(G) is the trace of the matrix G. If a solution is not found, we first bisect over δ, and if still there is no solution we increase k and repeat. We set the limits of δ to be 2 −6 and k to be 2.
The SOS decomposition and related optimization problems are generated using the symbolic computation package YALMIP [26,27]. Our MATLAB code and data for the experiments are available on https://github.com/Abhishek-B/PnCP, so that the reader may verify the results of our experiments.
To solve the required semi-definite programs, we use the MOSEK solver [28] with our implementations. Verification of the SOS decomposition is done with the YALMIP command sol.problem==0 (where sol is what we name our solution), as well as requiring the residual of the problem to be small (≤ O(10 −6 )).
All of the experiments were carried out on a standard Dell Optiplex 9020, with 12GB of memory, an Intel ®Core TM i5-4590 CPU @ 3.30GHz×4 processor, 500GB of storage and running Ubuntu 18.04 LTS.
The success rate of this relaxation for problems of small size is remarkable, as seen in Fig. 1 and Table 1. Moreover, we observe from the average residual (which includes the failed examples as well) in Table 1, that if we were to allow the residual to be slightly larger (say ≤ O(10 −4 )), we would see a higher success rate. This would also reduce computation times, increasing the appeal of this relaxation.
Remark 3. 1 We add tr(G) = 1 in our constraints to avoid the trivial solution of σ ≡ 0.
After running a few experiments, it becomes apparent that in the Hilbert method, we should initialize k = 2. While there are instances where k = 1 has a solution, it works with very small δ and hence requires a long runtime due to the number of bisections.

Remark 3.2
Under the experimental conditions described above, the Hilbert relaxation (as well as others considered in this work) performs very poorly on problems of larger size. For the Hilbert relaxation on problems of size (m, n) = (4, 4), we observe a success rate of 22%, and computation times exceeding 40 minutes. However, if we change the experiment parameters and aims, by fixing δ = 2 −5 , k = 2, and increase our residual threshold to O(10 −4 ), then we have observed the Hilbert relaxation achieve 100% success rate (even on problems of size (n,m)=(4,5)), with an average computation time of ∼ 180 seconds. For the practitioner working with a specific problem, this kind of heuristic approach is more appropriate.
The relaxation (4) is a simplified version of (5), which fixes the denominator We refer to this simplification as the Coordinate Norm Relaxation (CNR) and implement it similar to the Hilbert method. Since σ is known, we maximize δ and "bisect" over ≤ 2. The verification of a solution is also similar, with the additional requirement δ > O(10 −4 ) as otherwise δ becomes indistinguishable from numerical error. As we can see (Fig. 2 or Table 2), this relaxation is incredibly fast (it is in fact the fastest relaxation). On problems of smaller size, it is not as successful compared to the Hilbert method, but we can see from the residuals, that if we relax our verification criteria, we might improve the success rate of the CNR quite dramatically.
If we consider the variables z i j = x i ⊗ y j over the Segre variety, then the CNR can be written as: For strictly positive polynomials p, there always exists an such that the denominator (x 2 1 + · · · + x 2 n ) allows an SOS decomposition (see the second and third Theorems of   [37]). However the appropriate choice of depends on the minimum of p, and in fact For polynomials with zeros, this denominator has been used in practice (see [25] for instance), but there is little theoretical justification for its use. Algorithm 1 works by fixing some zeros of F δ in Step 1; hence, the relaxation (4), while practically efficient, is not guaranteed to work, jeopardizing the entire construction.

KKT Ideal
Because the constructed F δ is homogeneous, we can consider minimization over the sphere S N −1 (or any suitable compact set). In this case, a more modern SOS relaxation comes from considering the KKT conditions Letting I KKT be the ideal generated from the KKT conditions (stationarity, feasibility, and complementary slackness) and V R KKT its associated real variety. Then, [11,Theorem 3.2] provides us the following result, where we write w for the variables (x, y).
Note that since F δ is randomly generated, the assumption in Theorem 3.2 that I KKT be radical is generically true. For our experiments, we implement the same bisection strategy as Sect. 3.1 to deal with the nonlinearities, and restrict γ in (6) to γ = (0, 0), giving the program To our surprise, this method fails completely on the larger problems (Fig. 3) and has quite poor performance even on the smaller ones of size (3,3). This suggests that the random construction alone is not enough to guarantee the radicalness of I KKT .
Unlike the previous two relaxations, the residuals here do not indicate any room for improvement (Table 3). In our tests, increasing the relaxation degree k offers some success, but this also greatly increases the computation time (to the order of hours), making this relaxation impractical for the problem at hand.

Jacobian Relaxation
We now present an exact relaxation which (in theory) always works for our problem of interest. This approach is similar to the KKT relaxation, only now to establish the dependence between derivatives of the constraints (S N −1 ) and the function (F δ ), we consider determinants of an associated Jacobian matrix. We let S denote the equation of sphere as before and once again write w for the variables (x, y). Consider the matrix B(w) = (∇ F δ (w) ∇ S(w)). Denote by ϕ 1 , . . . , ϕ the minors of B(w); define J to be the ideal generated by S and ϕ 1 , . . . , ϕ . Then, it is shown in [30] that minimizing F δ on the sphere is equivalent to the following jacobian system Moreover, letting F * be the solution of (7) and F min the minimum of F δ , the following holds.
Due to nonlinearity in the constraints of (8), we employ the bisection approach similar to the other methods and solve find q(w) ∈ J , again with the lower limits of δ being 2 −6 and k being 2.
Unsurprisingly, this relaxation is quite slow. The solve time on test cases of size (3, 5) was close to one hour, and so we do not test the Jacobian relaxation on this set. We can also see (Fig. 4 or Table 4) that this relaxation exhibits low success rates and high residuals. Similar to KKT, the Jacobian relaxation is somewhat impractical in our context.

Remark 3.3
It should be noted again that these tests were conducted with limited freedom on the degrees of the relaxations. Based on our experience, we recommend using the Hilbert method with a high relaxation degree (k = 3) if memory is not a concern and the user wants more successful constructions. When memory becomes an issue, the CNR seems to be a better choice; although its success rate is lower, the speed of computation makes generating random examples more practical.

Rationalization
Constructing PnCP maps over floating point numbers provides quick numerical tests, which can indicate nonnegativity, but ideally we would like to have rational PnCP maps with exact certificates of nonnegativity. The semi-definite programs arising from our SOS relaxations are feasibility problems of the form max R 1, where A i and b i are obtained from the problem data (see [33] for a nice presentation of this). The following theorem, first proved in [35], provides a means to obtain rational solutions of (9) from numerical ones. (1) Compute a rational approximationG with τ := ||G −G|| satisfying τ 2 + 2 ≤ μ 2 .
(2) ProjectG onto the affine subspace L defined by the equations A i , G = b i , to obtainĜ.
For our problems, there are two key issues with using this rationalization. Firstly, our semi-definite programs will never satisfy the strict feasibility requirements of G being positive definite. This is because, by construction, the form F δ will always have non-trivial zeros chosen in Step 1 of Algorithm 1. We note that even the facial reduction methods, like those in [22] or [23], will in general not work in our setting. This is because, even though our initial choice of zeros for F δ may be rational, it is still possible for F δ to have some irrational zeros, which then cannot be removed via facial reduction. Second, and more importantly, the numbers b i are obtained from the coefficients of the polynomial being tested, in our case F δ . This means that the affine subspace L is being defined by floating point numbers, and any sort of rationalization of G will perturb this subspace.
In our package PnCP, we combat this by restricting the randomization in the linear algebra steps of Algorithm 1; we restrict the choices of the initial points x i , y j , so that the generated linear/quadratic forms have rational entries with small (single digit) denominators. Specifically, for constructing the initial points x i , y i , instead of using rand in MATLAB to generate uniform random numbers, we use randi to generate pseudorandom integers, with the integers belonging to the interval [−1, 1]. This ensures that the linear forms h i constructed using Step 1 of Algorithm 1 have rational coefficients with numerators and denominators having a few number of digits (two or three). We similarly restrict the randomization for the construction of the quadratic form f in Step 2.
As expected this reduces the base success rate of Algorithm 1, but it successfully constructs F δ with rational coefficients. We also observe a significant increase in computation time to construct forms with rational coefficients; we test this by constructing 50 random forms with rational coefficients and comparing the timing costs to constructing forms with floating point coefficients.
As we can see from Fig. 5, constructing rational forms is far more expensive than floating point forms. In fact, the average time taken to construct forms with floating point coefficients remains almost constant (∼2 s). In contrast, the (average) construction time for forms with rational coefficients can take close to 10 min. This rational construction can be used in PnCP with the command Gen_PnCP and setting the 'rationalize' argument to 1. Currently, PnCP provides numerical verification of the constructed rational F δ , via the techniques of Sect. 3. This construction can be used in conjunction with the many rational SOS packages (such as RationalSOS, RealCertify, multivsos, etc.) to obtain exact certificates of nonnegativity.

Detecting Quantum Entanglement
We will now show how we can use PnCP for detecting quantum entanglement. We start with a brief (and simplified) exposition into quantum states, the core object of interest for us, presenting some terminology and commonly known facts (for a more detailed introduction, we refer the reader to [2,19,39], or any graduate text on Quantum Information Theory). We then state two entanglement criteria and then give examples demonstrating how PnCP is used to implement the most general one.
A quantum state is a vector φ ∈ R n , and with any quantum state there is an associated density matrix φφ T =: ρ ∈ S n (normalized to have unit trace). A density matrix with {φ i } an orthonormal system, p i ≥ 0 and i p i = 1, represents a quantum system in one of several states φ i with associated probabilities p i . We use the following terminology; ρ is a pure state if ρ = φφ T , otherwise if ρ is of the form (10), then it is a mixed state. It should be noted that any positive semi-definite matrix ρ with tr(ρ) = 1 is a density matrix. It is known that pure states satisfy tr(ρ 2 ) = 1 while for mixed states tr(ρ 2 ) < 1. Given a composite quantum system S nm = S n ⊗ S m and a state ρ nm ∈ S nm , we call ρ nm simply separable if and entangled if it is not separable. A problem of interest in quantum information theory is the so-called separability problem; given a state (density matrix) ρ in a composite system, determine whether it is entangled.
There are many different criteria and measures of entanglement throughout the literature. For pure states, things are relatively simple and separability can be determined by checking whether the state is in the image of the Segre embedding. For mixed states however, the situation is more complicated.
In low-dimensional composite systems, we have the Peres-Horodecki criterion, also known as the positive partial transpose (PPT) criterion; for ρ nm = i p i ρ n i ⊗ ρ m i define the partial transpose map ( For systems of size (n, m) = (2, 2) or (2, 3), this criterion is both necessary and sufficient. In higher-dimensional systems, we lose the sufficiency of this test, i.e., there are entangled states ρ ent with (I ⊗ T )(ρ ent ) 0 (see [16] for the first such example). In this situation, we instead have the more general entanglement criteria.

Criterion 5.2 (The general criterion, [2] section 8.4) A quantum state
The PPT entanglement criterion is a special case of Criterion 5.2, with Φ being the transposition map. With PnCP, we can apply this test with many different random Φ as in Algorithm 2.

Algorithm 2: Entanglement Detection
Input: Quantum state ρ, maximum iteration number S. Output: Status. i = 0; Status = "Unknown"; while i < S do Generate random Φ, and compute I ⊗ Φ(ρ); if I ⊗ Φ(ρ) 0 then Status = "Entangled"; Stop; end end Example 5.1 As an example, consider the following generalized Bell state where each E i, j ∈ M 3 (R) is the matrix unit with 1 in row i, column j and zeros everywhere else. The Bell state is known to be entangled. We use PnCP to generate the following pncp map Φ Φ(E 1,1 ) = Since we construct Φ on S 3 , we make the canonical extension to M 3 (R) by setting With this extension, we find that (I ⊗ Φ)(Δ) has lowest (numerical) eigenvalue −8.45.
Example 5. 2 We consider now an example of a Bound Entangled State, which are known to be entangled whilst having a positive partial transpose (see [17] or [19,Section 6.11]). We take the example from [1, Section 4.6] (with parameter α = 4) It is known that this state is entangled and that the (pncp) Choi Map [9] demonstrates its entanglement. Using PnCP, we can find the alternative pncp map which shows that (I ⊗Φ)(Γ ) has the lowest (numerical) eigenvalue of −1.17. Searching for this map Φ took PnCP ∼ 9.4 seconds, and it took two iterations of the KMSZ construction. The point of this example is that with our package PnCP, there is now a viable alternative (which can be run in parallel) to searching the literature for an appropriate pncp map. Example 5. 3 We take another bound entangled example from [14], with Note that T r(Θ 2 ) = 0.2 < 1, and so Θ is a mixed state (meaning we cannot simply check whether it is in the image of the Segre embedding). PnCP generates the following We find the ampliation (I ⊗ Φ)(Θ) has lowest (numerical) eigenvalue −0.14. For this example, PnCP took ∼10 s to numerically check the entanglement status of the state, with majority of the time spent constructing the rational Φ. If we desired only an indication of entanglement, we could repeat this with Φ having floating point entries, for which the process takes ∼5 s.  [3]. Distillation of quantum states is beyond the scope of this article.
There are many other entanglement criteria that rely on testing a condition with some PnCP map. As we can see from the examples, PnCP provides a means to implement these criteria by being able to generate random (rational) pncp maps.

Conclusion
In this article, we present PnCP, a MATLAB package for constructing positive maps which are not completely positive, with a focus on the practicality of this construction, and we demonstrate its applicability to testing entanglement of quantum states.
PnCP is an open-source package available from https://github.com/Abhishek-B/ PnCP. The package implements state-of-the-art optimization techniques to numerically ensure positivity of the constructed maps. PnCP is even able to construct pncp maps with rational coefficients, which can be used in conjunction with existing software to obtain not only numerical, but exact certificates of positivity.
We use the KMSZ construction, which additionally provides a priori knowledge of some of the zeros of the constructed polynomial. While there is work on optimizing polynomials with zeros [8,36], there are restrictions on the zeros in these methods. Whether it is possible to adapt the zeros of the KMSZ construction to suit these methods, is something we wish to study in the future.
As the only package for this kind of construction, we intend to maintain and improve PnCP in various means; implementing better nonnegativity tests as they become available, optimizing the existing code (perhaps even pursuing parallel computing where possible), and including more entanglement criteria to improve the classification of quantum states.
Our primary focus moving forward will be to strengthen PnCP as a classification tool for quantum states; primarily by implementing a rational SOS decomposition method, which will automatically provide exact certificates of positivity. We will also be cataloging entangled quantum states and pncp maps which classify them on https:// github.com/Abhishek-B/PnCP. [15]. As in [22], we will write V (I n,m ) = [z 11 : . . . : z nm ] ∈ P nm−1 : f (z) = 0 for every f ∈ I n,m , for the Segre variety, where z = (z 11 , z 12 , . . . , z 1m , . . . , z nm ), and V R (I n,m ) for the subset of its real points. Finally, biquadratic forms in R[x, y] 2,2 are in a bijective correspondence with quadratic forms in R[z]/I n,m (see [22]). Let us write To obtain a quadratic form in P(V R (I n,m ))\Σ(V R (I n,m )) proceed as follows.
Step 1 Construction of linear forms h 0 , . . . , h d . In this step, we wish to find general linear forms h 0 , . . . , h d which intersect in n+m−2 n−1 distinct linearly general positions with at least (e+1) real, smooth intersection points; the existence of such forms is guaranteed (see the proof of [4, Proposition 3.2]). To construct these forms, we "work backwards", by choosing (e + 1) random points (Step 1.1) from which we build h 0 , . . . , h d , and then verify that they intersect in n+m−2 n−1 distinct points (Step 1.2). In PnCP we do not currently implement this verification (and this is one of the reasons why PnCP sometimes outputs a form which is a sum of squares). For small values of m, n, this verification can be done using Gröbner basis, but even for moderate sizes of m, n we found this verification to be computationally more expensive than the rest of the algorithm. However, since the points in Step 1.1 are chosen randomly, the number of points in the intersection "always" agrees with deg(V (I n,m )).
Step 1.1 Choose e + 1 random points x (i) ∈ R n and y (i) ∈ R m and calculate the Kronecker tensor products z (i) = x (i) ⊗ y (i) ∈ R nm . Step 1.2 Choose d random vectors v 1 , . . . v d ∈ R nm from the kernel of the matrix z (1) . . . z (e+1) * .
The corresponding linear forms h 1 , . . . , h d are If the number of points in the intersection is not equal to deg(V (I n,m )) = n+m−2 n−1 or if the points in the intersection are not in linearly general position, then repeat Step 1.1.

In
Step 2.1, we are computing a basis of each of these kernels, so that we can choose a quadratic form which satisfies (12). Note that once (12) is satisfied, then ∇ ⎛ ⎝ f − j λ j g j ⎞ ⎠ (z (i) ) = 0, together with the identity 2q(z) = (∇q(z)) * z, (13) for any quadratic form q, will yield f (z (i) ) = 0. So, once we have a feasible set for f in Step 2.1, we choose f randomly in Step 2.2 by utilizing the identity (13) (the first set of matrices in Step 2.2 is a tensor representation of the right side of (13)). Step 2.1 Let g 1 (z), . . . , g ( n 2 )( m 2 ) (z) be the generators of the ideal I n,m , i.e., the 2 × 2 minors z i j z kl − z il z k j for 1 ≤ i < k ≤ n, 1 ≤ j < l ≤ m. For each i = 1, . . . , e, compute a basis {w (i) 1 , . . . , w (i) d+1 } ⊆ R nm of the kernel of the matrix Note that this kernel is always (d +1)-dimensional, since the variety V (I n,m ) is d-dimensional (in P nm−1 ) and smooth.
Step 2.2 Let e i denote the i-th standard basis vector of the corresponding vector space, i.e., the vector with 1 on the i-th component and 0 elsewhere. Choose a random vector v ∈ R n 2 m 2 from the intersection of the kernels of the matrices with the kernels of the matrices e i ⊗ e j − e j ⊗ e i * , for 1 ≤ i < j ≤ nm.
The latter condition ensures v is a symmetric tensor in R nm ⊗ R nm . Note also that we have omitted the point z (e+1) . Define the quadratic form For 1 ≤ i, k ≤ n and 1 ≤ j, l ≤ m denote E i jkl = (e i ⊗ e j ) ⊗ (e k ⊗ e l ) + (e k ⊗ e l ) ⊗ (e i ⊗ e j ) ∈ R n 2 m 2 . Let and note that A is a representation of (symmetrized) a 2 , while B is a representation of the (symmetrized) generators of the ideal I n,m , namely g 1 , . . . , g ( n 2 )( m 2 ) . Since we wish to have f ∈ (R[z]/I n,m )\a 2 , we need f / ∈ a 2 as well as f / ∈ I n,m . In other words, if v is in does not belong to a 2 .
Step 3 Construction of a quadratic form in R[z]/I n,m that is positive but not a sum of squares. Calculate the greatest δ 0 > 0 such that δ 0 f + d i=0 h 2 i is nonnegative on V R (I n,m ). Then for every 0 < δ < δ 0 the quadratic form is nonnegative on V R (I n,m ) but is not a sum of squares.
As is explained in [22], with random data this algorithm works "well" without implementing verifications (for Step 1.2, etc.).