A simplex algorithm for rational cp-factorization

In this paper we provide an algorithm, similar to the simplex algorithm, which determines a rational cp-factorization of a given matrix, whenever the matrix allows such a factorization. This algorithm can be used to show that every integral completely positive $2 \times 2$ matrix has an integral cp-factorization.


Introduction
Copositive programming gives a common framework to formulate many difficult optimization problems as convex conic ones. In fact, many NP-hard problems are known to have such reformulations (see for example the surveys [3,10]). All the difficulty of these problems appears to be "converted" into the difficulty of understanding the cone of copositive matrices COP n which consists of all symmetric n × n matrices B ∈ S n with x T Bx ≥ 0 for all x ∈ R n ≥0 . Its dual cone is the cone CP n = cone{xx T : x ∈ R n ≥0 } of completely positive n × n matrices. Therefore, it seems no surprise that many basic questions about this cone are still open and appear to be very difficult.
One important problem is to find an algorithmic test deciding whether or not a given symmetric matrix A is completely positive. If possible one would like to obtain a certificate for either A ∈ CP n or A ∈ CP n . Dickinson and Gijben [7] showed that this (strong) membership problem is NP-hard.
In terms of the definitions the most natural certificate for A ∈ CP n is giving a cp-factorization For A ∈ CP n it is natural to give a separating hyperplane defined by a matrix B ∈ COP n so that the inner product of A and B satisfies B, A < 0. From the algorithmic side, different ideas have been proposed, in particular by Jarre and Schmallowsky [19], Nie [24], Sponsel and Dür [28], Elser [12], Groetzner and Dür [15], and Anstreicher, Burer, and Dickinson [6,Section 3.3]. Latter approach is is based on the ellipsoid method, while all others are non-exact numerical heuristics.
In this paper, in Section 3, we describe a new procedure that is based on pivoting like the simplex algorithm. For the defining the pivoting we apply the notion of the copositive minimum which we introduce in Section 2. Our algorithm works for all matrices in the coneC P n = cone{xx T : x ∈ Q n ≥0 }. Moreover, we conjecture that it always computes certificates, if the input matrix is not completely positive. In contrast to the non-exact approaches mentioned above, (the exception being the approach by Anstreicher, Burer and Dickinson whose factorization method is based on the ellipsoid method), our algorithm uses rational numbers only if the input matrix is rational and so does not face numerical instabilities. As a consequence, to the best of our knowledge, our algorithm is currently the only one that can find a rational cp-factorization whenever it exists. In [6] a similar result was obtained, but restricted to matrices in the interior of CP n .
If the input matrix A integral, one can also ask if it admits an integral cpfactorization, i.e. a cp-factorization of the form A = m i=1 x i x T i with x i ∈ Z n ≥0 and i = 1, . . . , m. For n ≥ 3 it is known that there are integral matrices A ∈ CP n which do not have an integral cp-factorization, see [1,Theorem 6.4]. For n = 2 it was conjectured by Berman and Shaked-Monderer [1, Conjecture 6.13] that every integral matrix A ∈ CP 2 possesses an integral cp-factorization. This conjecture was recently proved by Laffey andŠimgoc [21]. In Section 4 we show that our simplex algorithm can be used to give a short, alternative proof of this result.
In Section 5 we describe how an implementation of our algorithm performs on some examples.

The copositive minimum and copositive perfect matrices
By S n we denote the Euclidean vector space of symmetric n × n matrices with inner product A, B = Trace(AB) = n i,j=1 A ij B ij . With respect to this inner product we have the following duality relations between the cone of copositive matrices and the cone of completely positive matrices COP n = (CP n ) * = {B ∈ S n : A, B ≥ 0 for all A ∈ CP n }, and CP n = (COP n ) * .
So, in order to show that a given symmetric matrix A is not completely positive, it suffices to find a copositive matrix B ∈ COP n with B, A < 0. We call B a separating witness for A ∈ CP n in this case, because the linear hyperplane orthogonal to B separates A and CP n .
For a symmetric matrix B ∈ S n we define the copositive minimum as , and we denote the set of vectors attaining it by . The following proposition shows that matrices in the interior of the cone of copositive matrices attain their copositive minimum. Lemma 2.2. Let B be a matrix in the interior of the cone of copositive matrices. Then, the copositive minimum of B is strictly positive and it is attained by only finitely many vectors.
Proof. Since B is copositive, we have the inequality min COP (B) ≥ 0. Suppose that min COP (B) = 0. Then there is a sequence v i ∈ Z n ≥0 \ {0} of pairwise distinct lattice vectors such that B[v i ] tends to zero when i tends to infinity. From the sequence v i we construct a new sequence u i of vectors on the unit sphere S n−1 by setting v i = v i u i . The sequence u i belongs to the compact set R n ≥0 ∩ S n−1 . Thus, by taking a subsequence if necessary, we may assume that u i converges to a point u ∈ R n ≥0 ∩ S n−1 . The sequence of norms v i tends to infinity since the set of lattice vectors of bounded norm is finite. Thus we get which implies that B[u] = 0, contradicting our assumption B ∈ int(COP n ). Hence, min COP (B) > 0. By the same argument one can show that Min COP (B) only contains finitely many vectors.
In our previous paper [9] and in this paper the set This is because the minimum is attained at a vertex B * of P(A, λ) if and only if A is contained in the (inner) normal cone For a rational matrix A ∈ int(CP n ) we obtain a rational cp-factorization in this way, that is, a decomposition of the form However, for solving the linear program (2) one needs an explicit finite algorithmic description of the set P(A, λ), for example by a finite list of linear inequalities. The proof of the polyhedrality of P(A, λ) in [9, Lemma 2.4] relies on an indirect compactness argument (similar to the one in the proof of Lemma 2.2) which does not yield such an explicit algorithmic description. In the remainder of this paper we are therefore concerned with finding a finite list of linear inequalities.
The following definitions and the algorithm in the following section are inspired by Voronoi's classical algorithm for the classification of perfect positive definite quadratic forms. These can for instance be used to classify all locally densest lattice sphere packings (see for example [22] or [27]). In (3) we use the letter V to denote the normal cone of a vertex, as it is a generalization of the Voronoi cone used in the classical setting. In fact, our generalization of Voronoi's work can be viewed as an example of a broader framework described by Opgenorth [25]. In analogy with Voronoi's theory for positive definite quadratic forms we define the notion of perfectness for copositive matrices: n is a sum of squares. Thus Q An lies in the interior of the copositive cone. Furthermore, where e i is the i-th standard unit basis vector of R n . Thus, the n+1 2 vectors attaining the copositive minimum have a continued sequence of 1s in their coordinates and otherwise 0s. Now it is easy to see that the rank-1-matrices are linearly independent and span the space of symmetric matrices which shows that Q An is COP-perfect.
The matrix Q An is also known as a Gram matrix of the root lattice A n , a very important lattice, for instance in the theory of sphere packings (see for example [4]).

The algorithm
In this section we show how one can solve the linear program (2). Our algorithm is similar to the simplex algorithm for linear programming. It walks along a path of subsequently constructed COP-perfect matrices, which are vertices of the polyhedral set R that are connected by edges of R.
6. Use Algorithm 2 to determine the contiguous COP-perfect matrix B N := B P + λR with λ > 0 and min COP (B N ) = 1. 7. Set B P := B N and goto 2.

Algorithm 1. Algorithm to find a cp-factorization [ or a separating witness ]
3.1. Description of the algorithm. In Algorithm 1 we distinguish between two different kinds of inputs, as indicated by the brackets [ . . . ]. If A ∈CP n , then the algorithm always terminates when choosing the right pivots R in Step 4 (cf. Theorem 3.2). If A is in the cone S n >0 of positive definite matrices, then we can use Step 2 and Step 5 to certify A ∈ CP n with a witness matrix W . However, we do not know yet if Algorithm 1 always terminates in this case (cf. Conjecture 3.3).
In the following we describe the steps of Algorithm 1 in more detail: In Step 1, we can choose for instance the initial vertex B P = 1 2 Q An of R with Q An as in (4). Then the algorithm subsequently constructs vertices of R.
In Step 2 we check whether or not the current COP-perfect matrix is already a separating witness. By this, the algorithm subsequently constructs an outer approximation of the CP n cone: We shall show in Theorem 3.1 that this procedure gives a tighter and tighter outer approximation of the completely positive cone which eventually converges.
In Step 3 we determine whether A lies in the polyhedral cone V(B P ). This can be done by solving an auxiliary linear program The minimum equals 0 if and only if A lies in V(B P ). If A ∈ V(B P ), then we can find non-negative coefficients λ v , with v ∈ Min COP (B P ), to get a cp-factorization Using an algorithmic version of Carathéodory's theorem we can choose a sub- If the minimum of the auxiliary linear program (5) is negative we can find in Step 4 a generator R of an extreme ray of (V(B P )) * with A, R < 0. Here, several choices of R with A, R < 0 may be possible and the performance depends on the choices made in this "pivot step". A good heuristic for a "pivot rule" seems to be the choice of R with A, R/ R minimal, where R 2 = R, R . Also a random choice of R among the extreme rays of (V(B P )) * with A, R < 0 seems to perform quite well. In Step 5 it is checked whether or not R is a separating witness for A, that is, if not only A, R < 0 but also R ∈ COP n holds. We will explain the copositivity test of R in Section 3.3. In Step 6 Algorithm 2 is used to determine in direction of R ∈ COP n a new contiguous COP-perfect matrix B N of B P , that is, a neighboring vertex of B P on R, connected via an edge in direction R. Note that such a vertex exists (and R is not unbounded in the direction of R) under the assumption R ∈ COP n , because R ⊆ int(COP n ), see [9, Lemma 2.3].
Finally, we observe that since A, R < 0, we have A, B N < A, B P in each iteration (Steps 2 through 6) of the algorithm.

3.2.
Computing contiguous COP-perfect matrices. Our algorithm for computing contiguous COP-perfect matrices is inspired by a corresponding algorithm for computing contiguous perfect positive definite quadratic forms which is a subroutine in Voronoi's classical algorithm. In fact, the following algorithm is, after performing obvious changes, essentially identical with the one in [8, Section 6, Erratum to algorithm of Section 2.3].
Computationally the most involved parts of Algorithm 2 are checking if a matrix lies in the interior of the cone of copositive matrices, and if so, computing its copositive minimum min COP and all vectors Min COP attaining it. We discuss these tasks in Sections 3.3 and 3.4.
In the first while loop of Algorithm 2, lower and upper bounds l and u for the desired value λ are computed, such that B P + lR and B P + uR are lying in int(COP n ) satisfying min COP (B P + lR) = min COP (B P ) > min COP (B P + uR).
In the second while loop, the exact value of λ is determined. Note that replacing the assignment of u by the simpler assignment u := γ corresponds to a binary search coming at least arbitrarily close to λ. However, it may never reach the exact value. The assignment of u ensures that Input: COP-perfect matrix B P ∈ R, generator R ∈ COP n of extreme ray of the polyhedral cone (V(B P )) * Output: Neighboring vertex B N of B P on R, connected via an edge in direction R, i.e. 3.3. Checking copositivity. From a complexity point of view, checking whether or not a given symmetric matrix is copositive is known to be co-NP-complete by a result of Murty and Kabadi [23].
Nevertheless, in our algorithms we need to check whether or not a given symmetric matrix lies in the cone of copositive matrices (Step 5 in Algorithm 1) or in its interior (Step 2 of Algorithm 2). From a practical point of view, this can be checked by the following recursive characterization of Gaddum [14, Theorem 3.1 and 3.2]: By we denote the (n − 1)-dimensional standard simplex in dimension n. A matrix B ∈ S n lies in COP n (in int(COP n )) if and only if every of its principal minors of size (n − 1) × (n − 1) lies in COP n (in int(COP n−1 )) and the value of the two-player game with payoff matrix B is non-negative (strictly positive). One can compute the value of v in (6) by a linear program: v = max{λ : λ ∈ R, y ∈ ∆, By ≥ λe}, where e = (1, . . . , 1) T is the all-ones vector.

3.4.
Computing the copositive minimum. Once we know that a given symmetric matrix B lies in the interior of the copositive cone (i.e. after Step 2 of Algorithm 2) we apply the idea of simplex partitioning initially developed by Bundfuss and Dür [2] to compute its copositive minimum min COP (B) and all vectors Min COP (B) attaining it.
First we recall some facts and results from [2]. A family P = {∆ 1 , . . . , ∆ m } of simplices is called a simplicial partitioning of the standard simplex ∆ if Let v k 1 , . . . , v k n be the vertices of simplex ∆ k . It is easy to verify that if a symmetric matrix B ∈ S n satisfies the strict inequalities (7) (v k i ) T Bv k j > 0 for all i, j = 1, . . . , n, and k = 1, . . . , m, then it lies in int(COP n ). Bundfuss and Dür [2, Theorem 2] proved the following converse: Suppose B ∈ int(COP n ), then there exists an ε > 0 so that for all finite simplex partitions P = {∆ 1 , . . . , ∆ m } of ∆, where the diameter of every simplex ∆ k is at most ε, strict inequalities (7) hold. Here, the diameter of ∆ k is defined as max{ v k i − v k j : i, j = 1, . . . , n}. We assume now that B ∈ int(COP n ) and that we have a finite simplex partition P so that (7) holds. We furthermore assume that all the vertices v k i have rational coordinates. Such a simplex partition exists as shown by Bundfuss and Dür [2, Algorithm 2].
Each simplex ∆ k = conv{v k 1 , . . . , v k n } defines a simplicial cone by cone{v k 1 , . . . , v k n }. From now on we only work with the simplicial cones and not with the simplices any more, so we may scale the rational v k i 's to have integral coordinates. The goal is now to find all integer vectors v in ∆ k which minimize B [v]. To do this we adapt the algorithm of Fincke and Pohst [13], which solves the shortest lattice vector problem. It is the corresponding problem for positive semidefinite matrices. The adapted algorithm will solve the following problem: Given a matrix B ∈ int(CP n ) and a simplicial cone, which is generated by integer vectors v 1 , . . . , v n so that v T i Bv j ≥ 0 holds, and given a positive constant M , find all integer vectors v in the cone so that B[v] ≤ M holds. Then by reducing M successively to B[v], whenever such a non-trivial integer vector v is found, we can find the copositive minimum of B in the simplicial cone, as well as all integer vectors attaining it.
The first step of the algorithm is to compute the Hermite normal form of the matrix V which contains the the vectors v 1 , . . . , v n as it columns. (see for example Kannan and Bachem [20] or Schrijver [26], where it is shown that computing the Hermite normal form can be done in polynomial time). We find a unimodular matrix U ∈ GL n (Z) such that U V = W holds, where W is an upper triangular matrix with columns w 1 , . . . , w n and coefficients W i,j . Then, By B we denote the matrix (U −1 ) T BU −1 .
Again we use condition (iii) to get an upper bound for α n−1 : , and solving the corresponding quadratic equation gives the desired upper bound. Now suppose α n and α n−1 are fixed. We want to compute all possible values of the coefficient α n−2 and we can proceed inductively.
3.5. Analysis of the algorithm. In the following we prove some features of the algorithm.
We start with a general result which gives a tight outer approximation of the cone of completely positive matrices in terms of the boundary structure (its 1-skeleton to be precise) of the convex set R. Proof. Clearly, by [9,Lemma 2.3], the left hand side is contained in the right hand side because all vertices and all generators of extreme rays of R are elements of the (interior of the) copositive cone.
In order to show the reverse inclusion we use the fact where the identity K * = (int(K)) * is generally true for full dimensional convex cones.
Suppose for contradiction, we may assume that there is a Q in the right hand side of (9), for which there is a matrix B ∈ int(COP n ) so that Q, B < 0 holds. By Lemma 2.2 and appropriate positive scaling we may assume without loss of generality that B lies in R. Therefore, we find vertices B 1 , . . . , B k of R (COPperfect matrices) and generators R 1 , . . . , R l of extreme rays of R such that with suitable positive coefficients λ i and µ i , satisfying λ 1 + · · · + λ k = 1.
Because Q lies in the right hand side of (9) Q, B i ≥ 0 and Q, R i ≥ 0 for all i, and hence Q, B ≥ 0 contradicting our initial assumption. Proof. For A ∈ int(CP n ) ⊆CP n the assertion follows from Lemma 2.4 in [9].
So let us assume A ∈ bd CP n ∩CP n . Then A is in the relative interior of a proper face F ofCP n , spanned by finitely many rank-1 matrices xx T with x ∈ Z n ≥0 . This face F is contained in at least one cone V(B P ) of a perfect matrix B P (being w.l.o.g. a vertex of R). In fact, in any neighborhood of A there are interior points of CP n which are contained in one of the cones V(B P ), having F as one of its faces.
Let {R 1 , R 2 , . . .} be a possible sequence of generators of rays constructed in Step 4 of Algorithm 1. For all of these generators, the inequality A, R i < 0 holds. For k such generators, the conditions Q, R i < 0 for i = 1, . . . , k are not only satisfied for Q = A, but also for all Q in an ε-neighborhood of A (with a suitable ε depending on k). For any k, this neighborhood also contains points of int(V(B P )) ⊆ int(CP n ). For these interior points Q, however, Algorithm 1 finishes after at most finitely many steps (when checking for Q ∈ V(B P ) in Step 3). Thus, for some finite number of suitable choices in Step 4, the algorithm also ends for A (when checking for A ∈ V(B P ), and then returning X = Min COP (B P )).
In principle we could use breadth-first-search in Step 4 of Algorithm 1, to ensure finite termination, but this of course would be far less efficient.
For the case of A ∈ CP n , we do not yet know if it is possible that Algorithm 1 does not provide a separating witness W after finitely many iterations. With a suitable chosen rule in Step 4, however, we conjecture that the computation finishes with a certificate: Conjecture 3.3. For A ∈ CP n , Algorithm 1 with a suitable "pivot rule" in Step 4 ends after finitely many iterations with a separating witness W .
We close this subsection with a few observations that can be made in the remaining "non-rational boundary cases", that is, for A ∈ bd CP n \CP n . In this case, Algorithm 1 may not terminate after finitely many steps, as shown in a 2-dimensional example in the following section. Assuming there is an infinite sequence of vertices B (i) P of R constructed in Algorithm 1, we know however at least the following: (i) The perfect matrix B Integers α, β indicate that the shown point is on a ray spanned by the rank- contains a convergent subsequence with limit B ∈ {X ∈ S n : X, A = 0}. It can be shown that this B is in bd COP n . Infinite sequences of vertices B (i) P of R with such a limit B exist. For n = 2 we give an example in Section 4, in which A is from the "irrational boundary part" (bd CP n ) \CP n .

A 2-dimensional example
In this section we demonstrate how Algorithm 1 works for n = 2. Thereby, we discover a relation of our algorithm to beautiful classical results in elementary number theory. In particular, we consider the case when the input matrix A lies on the boundary of CP 2 . For n = 2 the boundary cases can be understood (up to scaling) to lie on a half circle, see Figure 1.
The boundary of CP 2 splits into a part of diagonal matrices and into rank-1 matrices A = xx T . In the first case, Algorithm 1 finishes already in its first iteration, if we use Q A2 as a starting perfect matrix, 1 where Let us consider the other boundary cases for n = 2, where A = xx T is a rank-1 matrix. Without loss of generality we can assume that x = (α, 1) T . As we explain in the following, Algorithm 1 will terminate after finitely many iterations 1 Strictly speaking we should use 1 2 Q A 2 here. If we use Q A 2 instead, then the algorithm produces integral matrices and vertices of 2R.
with a COP-perfect matrix B P satisfying x ∈ Min COP B P when α is rational. For irrational α the algorithm will not terminate.
The first observation is that Algorithm 1 subsequently replaces a COP-perfect matrix B P by a contiguous COP-perfect matrix B N in a way that one of the three vectors in Min COP (B P ) is replaced by the sum of the remaining two. Let B P be a copositive matrix with (11) Min Then, det ( a c b d ) = ±1 and we get a contiguous COP-perfect matrix B N with For instance, starting with B P = Q A2 as in (10) 1 0 is replaced by Note that for α = 1, Algorithm 1 also finishes already in the first iteration. The way these vectors are constructed corresponds to the way the famous Farey sequence is obtained. This relation between the Farey diagram/sequence and quadratic forms was first investigated in a classical paper of Adolf Hurwitz [18] in 1894 inspired by a lecture of Felix Klein; see also the book by Hatcher [16], which contains the proofs. However, every one of the infinitely many perfect matrices B Laffey andŠimgoc [21] showed that every integral matrix A ∈ CP 2 possesses an integral cp-factorization. This can also be seen as follows: If B P is a copositive matrix with (11) then the matrices e f e f T form a Hilbert basis of the convex cone which they generate. This means that every integral matrix in this cone is an integral combination of the three matrices above.
To show this, one immediately verifies this fact in the special case of B P = Q A2 . Then all the other cones are equivalent by conjugating with a matrix in GL 2 (Z).

Computational Experiments
We implemented our algorithm. The source code, written in C++, is available on GitHub [17]. In this section we report on the performance on several examples, most of them previously discussed in the literature. Generally, the running time of the procedure is hard to predict. The number of necessary iterations in Algorithm 1 drastically varies in the considered examples. Most of the computational time is taken by the computation of the copositive minimum as described in Section 3.4.

Matrices in the interior.
For matrices in the interior of the completely positive cone, our algorithm terminates with a certificate in form of a cp-factorization. Note that in [11] and in [5] characterizations of matrices in the interior of the completely positive cone are given. For example, we have that A ∈ int(CP n ) if and only if A has a factorization A = BB T with B > 0 and rank B = n.

5.2.
Matrices on the boundary. For matrices inCP n there exists a cp-factorization by definition. However, on the boundary of the cone these are often difficult to find.
The following example is from [15] and lies in the boundary ofCP 5 :         [15] took roughly 10 days and 70 iterations to find a factorization with only three vectors (3,5,8,8,2), (4, 1, 7, 2, 5) and (4,6,7,6,6). Both algorithms suggested in [15] were not able to find any cp-factorization. We also considered the following family of completely positive (n + m) × (n + m) matrices, generalizing the family of examples considered in [19]: The matrices n Id m J m,n J n,m m Id n , with J ·,· denoting an all-ones matrix of suitable size, are known to have cp-rank nm, that is, they have a cp-factorization with nm vectors, but not with less. These factorizations are found by our algorithm with starting COP-perfect matrix Q Am+n for all n, m ≤ 3 in less than 6 iterations.