365体育网站

Advertisement

Springer Nature is making SARS-CoV-2 and COVID-19 research free. View research | View latest news | Sign up for updates

A simplex algorithm for rational cp-factorization

Abstract

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\)365体育网站 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 \(\text {NP}\)-hard problems are known to have such reformulations (see for example the surveys [4, 11]). All the difficulty of these problems appears to be “converted” into the difficulty of understanding the cone of copositive matrices \({\mathcal {COP}}^{n}\) which consists of all symmetric \(n \times n\) matrices \(B \in {\mathcal {S}}^n\) with \(x^\mathsf{T} B x \ge 0\) for all \(x \in {\mathbb {R}}^n_{\ge 0}\). Its dual cone is the cone

$$\begin{aligned} \begin{aligned} {\mathcal {CP}}^{n}\; = \;&{{\,\mathrm{cone}\,}}\{xx^\mathsf{T} : x \in {\mathbb {R}}^n_{\ge 0}\}\\ \; = \;&\left\{ \sum _{i=1}^m \alpha _i x_i x_i^\mathsf{T} : m \in {\mathbb {N}}, \alpha _i \in {\mathbb {R}}_{\ge 0}, x_i \in {\mathbb {R}}^n_{\ge 0}, i = 1, \ldots , m\right\} \end{aligned} \end{aligned}$$

of completely positive \(n\times 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 \in {\mathcal {CP}}^{n}\) or \(A \not \in {\mathcal {CP}}^{n}\). Dickinson and Gijben [8] showed that this (strong) membership problem is \(\text {NP}\)-hard.

In terms of the definitions the most natural certificate for \(A \in {\mathcal {CP}}^{n}\) is giving a cp-factorization

$$\begin{aligned} A = \sum _{i=1}^m \alpha _i x_i x_i^\mathsf{T} \quad \text {with} \quad m \in {\mathbb {N}}, \alpha _i \in {\mathbb {R}}_{\ge 0}, \; x_i \in {\mathbb {R}}^n_{\ge 0}, \; i = 1, \ldots , m. \end{aligned}$$
(1)

For \(A \not \in {\mathcal {CP}}^{n}\) it is natural to give a separating hyperplane defined by a matrix \(B \in {\mathcal {COP}}^{n}\) so that the inner product of A and B satisfies \(\langle B, A \rangle < 0\).

365体育网站From the algorithmic side, different ideas have been proposed. One can divide the relevant literature according to two complementary approaches:

  1. (1)

    Numerical methods which are practical but “only” can find approximate cp-factorizations. The papers by Jarre and Schmallowsky [20], Nie [25], Sponsel and Dür [30], Elser [13], Groetzner and Dür [16] fall into this category.

  2. (2)

    Theoretical methods which can compute exact cp-factorizations in finitely many algorithmic steps. The factorization method of Anstreicher, Burer, and Dickinson [7, Section 3.3] uses the ellipsoid method and works for all matrices which have a rational cp-factorization and lie in the interior of the cone \({\mathcal {CP}}^{n}\). Berman and Rothblum [1] use quantifier elimination for first order formulae over the reals to compute the \(\mathcal {CP}\)-rank of a given matrix, that is, the minimum number m of vectors used in a cp-factorization (1).

In this paper, in Sect. 3, we describe a new procedure that is based on pivoting like the simplex algorithm. To define the pivoting we apply the notion of the copositive minimum which we introduce in Sect. 2. Our algorithm (Algorithm 1) works for all matrices in the rational cone

$$\begin{aligned} \begin{aligned} \widetilde{\mathcal {CP}}^{n}\; = \;&{{\,\mathrm{cone}\,}}_{{\mathbb {Q}}}\{xx^\mathsf{T} : x \in {\mathbb {Q}}^n_{\ge 0}\} \\ \; = \;&\left\{ \sum _{i=1}^m \alpha _i x_i x_i^\mathsf{T} : m \in {\mathbb {N}}, \alpha _i \in {\mathbb {Q}}_{\ge 0}, x_i \in {\mathbb {Q}}^n_{\ge 0}, i = 1, \ldots , m\right\} . \end{aligned} \end{aligned}$$

Moreover, we conjecture that a variant of our algorithm (Procedure 3) always computes separating hyperplanes, if the input matrix is not completely positive. Overall, our procedure works for matrices with coefficients in any computable subfield F of the real numbers, in that case the coefficients \(\alpha _i\) of the formula above belong to \(F_{\ge 0}\) and the whole algorithmic procedure works similarly as the rational case that we consider in this paper.

Our algorithm uses rational numbers only if the input matrix is rational and so allows in principle exact computations. 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 [7] a similar result was obtained, but restricted to matrices in the interior of \({\mathcal {CP}}^{n}\). A related question is if every rational completely positive matrix has a rational cp-factorization, see the survey [29]. Generally we do not know but from the results in [7] and [10] it follows that this is true for matrices in the interior of \({\mathcal {CP}}^{n}\).

If the input matrix A is integral, one can also ask if it admits an integral cp-factorization, i.e. a cp-factorization of the form \(A = \sum _{i = 1}^m x_i x_i^\mathsf{T}\) with \(x_i \in {\mathbb {Z}}_{\ge 0}^n\) for all \(i = 1, \ldots , m\). For \(n \ge 3\) it is known that there are integral matrices \(A \in {\mathcal {CP}}^{n}\) which do not have an integral cp-factorization, see [2, Theorem 6.4]. For \(n = 2\) it was conjectured by Berman and Shaked-Monderer [2, Conjecture 6.13] that every integral matrix \(A \in {\mathcal {CP}}^{2}\) possesses an integral cp-factorization. This conjecture was recently proved by Laffey and Šimgoc [22]. In Sect. 4365体育网站 we show that our simplex algorithm can be used to give a short, alternative proof of this result.

In Sect. 5365体育网站 we describe how an implementation of our algorithm performs on some examples.

The copositive minimum and copositive perfect matrices

Copositive minimum

By \({{\mathcal {S}}}^{n}\) we denote the Euclidean vector space of symmetric \(n\times n\) matrices with inner product \(\langle A, B \rangle = {{\,\mathrm{Trace}\,}}(AB) = \sum _{i,j=1}^n 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

$$\begin{aligned} {\mathcal {COP}}^{n}= ({\mathcal {CP}}^{n})^* = \{B \in {\mathcal {S}}^n : \langle A, B \rangle \ge 0 \text { for all } A \in {\mathcal {CP}}^{n}\}, \end{aligned}$$

and

$$\begin{aligned} {\mathcal {CP}}^{n}= ({\mathcal {COP}}^{n})^*. \end{aligned}$$

So, in order to show that a given symmetric matrix A is not completely positive, it suffices to find a copositive matrix \(B \in {\mathcal {COP}}^{n}\) with \(\langle B, A\rangle < 0\). We call B a separating witness for \(A \not \in {\mathcal {CP}}^{n}\) in this case, because the linear hyperplane orthogonal to B separates A and \({\mathcal {CP}}^{n}\).

Using the notation B[x] for \(x^\mathsf{T} B x = \langle B, xx^\mathsf{T} \rangle \), we obtain

$$\begin{aligned} {\mathcal {COP}}^{n}= \{B \in {\mathcal {S}}^n : B[x] \ge 0 \text { for all } x\in {\mathbb {R}}^n_{\ge 0}\}. \end{aligned}$$

Obviously, the cone of positive semidefinite matrices \({{\mathcal {S}}}^{n}_{\ge 0}\), whose interior is the open cone of positive definite matrices \({{\mathcal {S}}}^{n}_{>0}\), lies between the completely positive cone and the copositive cone: \({\mathcal {CP}}^{n}\subseteq {{\mathcal {S}}}^{n}_{\ge 0}\subseteq {\mathcal {COP}}^n\).

Definition 2.1

For a symmetric matrix \(B\in {{\mathcal {S}}}^{n}\) we define the copositive minimum as

$$\begin{aligned} {{\,\mathrm{min_{\mathcal {COP}}}\,}}(B) = \inf \left\{ B[v] : v \in {{\mathbb {Z}}}^n_{\ge 0}{\setminus }\{0\} \right\} , \end{aligned}$$

and we denote the set of vectors attaining it by

$$\begin{aligned} {{\,\mathrm{Min_{\mathcal {COP}}}\,}}(B) = \left\{ v \in {{\mathbb {Z}}}^n_{\ge 0} : B[v] = {{\,\mathrm{min_{\mathcal {COP}}}\,}}(B) \right\} . \end{aligned}$$

365体育网站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 \({{\,\mathrm{min_{\mathcal {COP}}}\,}}(B) \ge 0\). Suppose that \({{\,\mathrm{min_{\mathcal {COP}}}\,}}(B) = 0\). Then there is a sequence \(v_i \in {{\mathbb {Z}}}_{\ge 0}^n {\setminus } \{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 = \Vert v_i\Vert u_i\). The sequence \(u_i\) belongs to the compact set \({{\mathbb {R}}}^n_{\ge 0} \cap S^{n-1}\). Thus, by taking a subsequence if necessary, we may assume that \(u_i\) converges to a point \(u\in {{\mathbb {R}}}^n_{\ge 0} \cap S^{n-1}\). The sequence of norms \(\Vert v_i\Vert \) tends to infinity since the set of lattice vectors of bounded norm is finite. Thus we get

$$\begin{aligned} 0 = \lim _{i \rightarrow \infty } B[v_i] = \lim _{i \rightarrow \infty } \Vert v_i\Vert ^2 B[u_i], \end{aligned}$$

which implies that \(B[u] = 0\), contradicting our assumption \(B \in {{\,\mathrm{int}\,}}({\mathcal {COP}}^{n})\). Hence, \({{\,\mathrm{min_{\mathcal {COP}}}\,}}(B) > 0\).

By the same argument one can show that \({{\,\mathrm{Min_{\mathcal {COP}}}\,}}(B)\) only contains finitely many vectors. \(\square \)

A locally finite polyhedron

In our previous paper [10365体育网站] and in this paper the set

$$\begin{aligned} {\mathcal {R}} = \{B \in {{\mathcal {S}}}^{n}: B[v] \ge 1 \text { for all } v \in {\mathbb {Z}}^n_{\ge 0} {\setminus } \{0\} \} \end{aligned}$$

plays a central role.Footnote 1 The set \({\mathcal {R}}\) is a locally finite polyhedron, meaning that every intersection of \({\mathcal {R}}\) with a polytope is a polytope itself. In [10, Lemma 2.3] we showed that \({\mathcal {R}}\) is contained in the interior of the cone of copositive matrices. Thus, we can rewrite \({\mathcal {R}}\) as

$$\begin{aligned} {\mathcal {R}} = \{B \in {{\mathcal {S}}}^{n}: {{\,\mathrm{min_{\mathcal {COP}}}\,}}(B) \ge 1 \}. \end{aligned}$$
(2)

Note also that [10, Lemma 2.3] together with Lemma 2.2 implies

$$\begin{aligned} {{\,\mathrm{cone}\,}}\mathcal {R} \setminus \{0\} = {{\,\mathrm{int}\,}}{\mathcal {COP}}^n \end{aligned}$$
(3)

The following theorem 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 \({\mathcal {R}}\). Similarly, Yıldırım [32] discusses uniform polyhedral approximations of the the cone of copositive matrices.

Theorem 2.3

We have

$$\begin{aligned} \begin{aligned} {\mathcal {CP}}^{n}= \{ Q\in {{\mathcal {S}}}^{n}: \langle Q, B \rangle \ge 0&\text { for all vertices and for all generators}\\&\hbox { of extreme rays } B \hbox { of } {\mathcal {R}} \}. \end{aligned} \end{aligned}$$
(4)

Proof

We have

$$\begin{aligned} {\mathcal {CP}}^{n}= ({\mathcal {COP}}^{n})^* = ({{\,\mathrm{int}\,}}({\mathcal {COP}}^{n}))^* = ({{\,\mathrm{cone}\,}}{\mathcal {R}})^*, \end{aligned}$$

where the identity \(K^* = ({{\,\mathrm{int}\,}}(K))^*\) is generally true for full dimensional convex cones and the last identity is (3). Since \({\mathcal {R}}\) is a locally finite polyhedron, \(({{\,\mathrm{cone}\,}}{\mathcal {R}})^*\) is equal to the right hand side of (4). \(\square \)

A linear program for finding a rational cp-factorization

In [10, Lemma 2.4] we showed that for \(A\in {{\,\mathrm{int}\,}}({\mathcal {CP}}^{n})\) and all sufficiently large \(\lambda >0\) the set

$$\begin{aligned} {\mathcal {P}}(A,\lambda ) = \{B\in {\mathcal {R}} : \langle A, B \rangle \le \lambda \} \end{aligned}$$

365体育网站is a full-dimensional polytope.

In principle (cf. [10, Proof of Theorem 1.1]), this gives a way to compute a cp-factorization for a given matrix \(A\in {{\,\mathrm{int}\,}}({\mathcal {CP}}^{n})\)365体育网站 by solving the linear program

$$\begin{aligned} \min \left\{ \langle A, B \rangle : B \in {\mathcal {P}}(A,\lambda )\right\} : \end{aligned}$$
(5)

This is because the minimum is attained at a vertex \(B^*\) of \({\mathcal {P}}(A,\lambda )\). Hence, due to the minimality of \(\langle A, B^* \rangle \), the matrix A is contained in the (inner) normal cone

$$\begin{aligned} {{\mathcal {V}}}(B^*) = {{\,\mathrm{cone}\,}}\left\{ vv^\mathsf{T} : v \in {{\,\mathrm{Min_{\mathcal {COP}}}\,}}B^* \right\} \end{aligned}$$
(6)

of \({\mathcal {R}}\) at \(B^*\). For a rational matrix \(A\in {{\,\mathrm{int}\,}}({\mathcal {CP}}^{n})\) we obtain a rational cp-factorization365体育网站 in this way, that is, a decomposition of the form

$$\begin{aligned} A = \sum _{i=1}^m \alpha _i v_i v_i^\mathsf{T} \quad \text {with } \alpha _i \in {{\mathbb {Q}}}_{\ge 0 }\text { and } v_i \in {{\mathbb {Z}}}^n_{\ge 0}, \quad \text {for } i = 1, \ldots , m. \end{aligned}$$
(7)

To find this factorization, we apply Carathéodory’s theorem (see for example [27, Corollary 7.1i]) and choose a subset \(v_1, \ldots , v_m \in {{\,\mathrm{Min_{\mathcal {COP}}}\,}}B^*\) such that \(v_iv_i^\mathsf{T}\) are linearly independent and \(A \in {{\,\mathrm{cone}\,}}\{v_i v_i^\mathsf{T} : i = 1, \ldots , m\}\). So we can find unique non-negative rational coefficients \(\alpha _1, \ldots , \alpha _m\) giving the rational cp-factorization (7).

However, for solving the linear program (5) one needs an explicit finite algorithmic description of the set \({\mathcal {P}}(A, \lambda )\), for example by a finite list of linear inequalities. The proof of the polyhedrality of \({\mathcal {P}}(A, \lambda )\) in [10, Lemma 2.4] relies on an indirect compactness argument (similar to the one in the proof of Lemma 2.2365体育网站) 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.

Copositive perfect matrices

In the next step we characterize the vertices of \({\mathcal {R}}\). 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 [23] or [28]). In (6) we use the letter \({{\mathcal {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 [26]. In analogy with Voronoi’s theory for positive definite quadratic forms we define the notion of perfectness for copositive matrices:

Definition 2.4

A copositive matrix \(B \in {{\,\mathrm{int}\,}}({\mathcal {COP}}^{n})\) is called \({\mathcal {COP}}\)-perfect if it is uniquely determined by its copositive minimum \({{\,\mathrm{min_{\mathcal {COP}}}\,}}B\) and the set \({{\,\mathrm{Min_{\mathcal {COP}}}\,}}B\) attaining it.

In other words, \(B \in {{\,\mathrm{int}\,}}({\mathcal {COP}}^{n})\) is \({\mathcal {COP}}\)-perfect if and only if it is the unique solution X of the system of linear equations

$$\begin{aligned} \langle X, vv^\mathsf{T}\rangle = {{\,\mathrm{min_{\mathcal {COP}}}\,}}B, \text { for all } v \in {{\,\mathrm{Min_{\mathcal {COP}}}\,}}B. \end{aligned}$$

Hence, \({\mathcal {COP}}\)-perfect matrices are, up to scaling, exactly the vertices of \({\mathcal {R}}\).

Lemma 2.5

\({\mathcal {COP}}\)-perfect matrices exist in all dimensions (dimension \(n = 1\) being trivial): For dimension \(n \ge 2\) the following matrix

$$\begin{aligned} Q_{{\mathsf {A}}_n} = \begin{pmatrix} 2 &{}\quad -1 &{}\quad 0 &{}\quad \ldots &{}\quad 0 \\ -1 &{}\quad 2 &{}\quad \ddots &{}\quad \ddots &{}\quad \vdots \\ 0 &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots &{}\quad 0 \\ \vdots &{}\quad \ddots &{}\quad \ddots &{}\quad 2 &{}\quad -1 \\ 0 &{}\quad \ldots &{}\quad 0 &{}\quad -1 &{}\quad 2 \\ \end{pmatrix} \end{aligned}$$
(8)

is \({\mathcal {COP}}\)-perfect; \(\frac{1}{2} Q_{{\mathsf {A}}_n}\) is a vertex of \({\mathcal {R}}\).

The matrix \(Q_{{\mathsf {A}}_n}\) is also known as a Gram matrix of the root lattice \({\mathsf {A}}_n\), a very important lattice, for instance in the theory of sphere packings (see for example [5]).

Proof

The matrix \(Q_{{\mathsf {A}}_n}\) is positive definite since

$$\begin{aligned} Q_{{\mathsf {A}}_n}[x] = x_1^2 + \sum _{i=1}^{n-1} (x_i - x_{i+1})^2 + x_n^2 \end{aligned}$$

is a sum of squares and \(Q_{{\mathsf {A}}_n}[x] = 0\) if and only if \(x = 0\). Thus \(Q_{{\mathsf {A}}_n}\) lies in the interior of the copositive cone. Furthermore,

$$\begin{aligned} {{\,\mathrm{min_{\mathcal {COP}}}\,}}Q_{{\mathsf {A}}_n} = 2 \quad \text {with} \quad {{\,\mathrm{Min_{\mathcal {COP}}}\,}}Q_{{\mathsf {A}}_n} = \left\{ \sum _{i=j}^k e_i : 1 \le j \le k \le n\right\} , \end{aligned}$$

where \(e_i\) is the i-th standard unit basis vector of \({\mathbb {R}}^n\). Thus, the \({n+1 \atopwithdelims ()2}\)365体育网站 vectors attaining the copositive minimum have a continued sequence of 1s in their coordinates and 0s otherwise. Now it is easy to see that the rank-1-matrices

$$\begin{aligned} \left( \sum _{i=j}^k e_i\right) \left( \sum _{i=j}^k e_i\right) ^\mathsf{T}, \text { where } 1 \le j \le k \le n, \end{aligned}$$

are linearly independent and span the space of symmetric matrices which shows that \(Q_{{\mathsf {A}}_n}\) is \({\mathcal {COP}}\)-perfect. \(\square \)

Algorithms

In this section we show how one can solve the linear program (5). Our algorithm is similar to the simplex algorithm for linear programming. It walks along a path of subsequently constructed \({\mathcal {COP}}\)-perfect matrices, which are vertices of the polyhedral set \({\mathcal {R}}\) that are connected by edges of \({\mathcal {R}}\).

We start with a simple version assuming that the input matrix lies in \(\widetilde{\mathcal {CP}}^{n}\). Of course, this assumption can usually not been easily checked beforehand and the rational cp-factorization is only given as the output of the algorithm. In this sense, the algorithm gets the promise that the input matrix possesses a rational cp-factorization. In theoretical computer science promise problems are common; for practical purposes we propose an extended procedure at the end of this section, see Procedure 3.

figurea

Description and analysis of the algorithm

In the following we describe the steps of Algorithm 1 in more detail:

In Step 1, we can choose for instance the initial vertex \(P=\frac{1}{2}Q_{{\mathsf {A}}_n}\) of \({{\mathcal {R}}}\) with \(Q_{{\mathsf {A}}_n}\) as in (8). Then the algorithm subsequently constructs vertices of \({{\mathcal {R}}}\).

In Step 2 we determine whether A lies in the polyhedral cone \({{\mathcal {V}}}(P)\). For this we consider all \(v\in {{\,\mathrm{Min_{\mathcal {COP}}}\,}}(P)\) giving generators \(vv^\mathsf{T}\) of the polyhedral cone \({{\mathcal {V}}}(P)\), respectively defining linear inequalities of the dual cone \(({{\mathcal {V}}}(P))^*\). Testing \(A\in {{\mathcal {V}}}(P)\) can then be done by solving an auxiliary linear program

$$\begin{aligned} \min \left\{ \langle A,Q \rangle : Q \in ({{\mathcal {V}}}(P))^*\right\} . \end{aligned}$$
(9)

The minimum equals 0 if and only if A lies in \({{\mathcal {V}}}(P)\). If \(A \in {{\mathcal {V}}}(P)\), then we can find non-negative coefficients \(\lambda _v\), with \(v \in {{\,\mathrm{Min_{\mathcal {COP}}}\,}}(P)\), to get a cp-factorization

$$\begin{aligned} A = \sum _{v \in {{\,\mathrm{Min_{\mathcal {COP}}}\,}}(P)} \lambda _v vv^\mathsf{T}. \end{aligned}$$

Using (an algorithmic version of) Carathéodory’s theorem we can choose in Step 3 a subset \(\{v_1, \ldots , v_m\} \subseteq {{\,\mathrm{Min_{\mathcal {COP}}}\,}}(P)\) so that we get a rational cp-factorization \(A = \sum _{i = 1}^m \alpha _i v_iv_i^\mathsf{T}\) with non-negative rational numbers \(\alpha _i\); see Sect. 2.3.

If the minimum of the auxiliary linear program (9) is negative we can find in Step 2(a) a generator R of an extreme ray of \(({{\mathcal {V}}}(P))^{*}\) with \(\langle A,R \rangle < 0\). Here, several choices of R with \(\langle A,R \rangle < 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 \(\langle A,R / \Vert R\Vert \rangle \) minimal, where \(\Vert R\Vert ^2 = \langle R, R\rangle \). Also a random choice of R among the extreme rays of \(({{\mathcal {V}}}(P))^{*}\) with \(\langle A,R \rangle < 0\) seems to perform quite well. When choosing the right pivots R in Step 2(a) Algorithm 1 always terminates, as shown by Theorem 3.1 below.

In Step 2(b) Algorithm 2 (see Sect. 3.2 is used to determine a new contiguous\({\mathcal {COP}}\)-perfect matrixN of P in direction of \(R\not \in {\mathcal {COP}}^{n}\), that is, a contiguous vertex of P on \({\mathcal {R}}\), connected via an edge in direction R. Note that such a vertex exists (and \({\mathcal {R}}\) is not unbounded in the direction of R) under the assumption \(R \not \in {\mathcal {COP}}^{n}\), because \({\mathcal {R}} \subseteq {{\,\mathrm{int}\,}}({\mathcal {COP}}^{n})\), see [10, Lemma 2.3]. We can exclude \(R \in {\mathcal {COP}}^{n}\) here, since together with \(\langle A,R \rangle < 0\) it would contradict the promise \(A\in \widetilde{\mathcal {CP}}^{n}\) on the input. Note also that as a byproduct of Algorithm 2 we compute generators of the cone \({\mathcal {V}}(N)\).

Finally, we observe that since \(\langle A,R \rangle < 0\), we have \(\langle A,N \rangle < \langle A,P \rangle \) in each iteration (Step 2) of the algorithm.

The following theorem shows that we can set up an algorithm for the promise problem.

Theorem 3.1

For \(A\in \widetilde{\mathcal {CP}}^{n}\), Algorithm 1 with suitable choices in Step 2(a) ends after finitely many iterations giving a rational cp-factorization of A.

In particular, with breadth-first-search added to Algorithm 1 we can guarantee finite termination (but this of course would be far less efficient).

Proof

For \(A\in {{\,\mathrm{int}\,}}({\mathcal {CP}}^{n}) \cap \widetilde{\mathcal {CP}}^{n}\) the assertion follows from Lemma 2.4 in [10].

So let us assume \(A\in {{\,\mathrm{bd}\,}}{\mathcal {CP}}^{n}\cap \widetilde{\mathcal {CP}}^{n}\). Note that \(\widetilde{\mathcal {CP}}^{n}\) is tessellated into cones \({{\mathcal {V}}}(P)\) of \({\mathcal {COP}}\)-perfect matrices P. In fact, the convex hull of \(D = \{xx^\mathsf{T} : x \in {{\mathbb {Z}}}^n_{\ge 0}, x \ne 0\}\) is a locally finite polyhedral set whose facets are in one-to-one correspondence with the \({\mathcal {COP}}\)-perfect matrices (see [26]). For any \(A\in \widetilde{\mathcal {CP}}^{n}\), the ray \(\{ \lambda A : \lambda \ge 0\}\) meets (at least) one facet of \({{\,\mathrm{conv}\,}}D\) and A is in \({{\mathcal {V}}}(P)\) of the corresponding \({\mathcal {COP}}\)-perfect matrix P.

Let \(\{ R_1, R_2, \ldots \}\) be a possible sequence of generators of rays constructed in Step 2(a) of Algorithm 1. For all of these generators, the inequality \(\langle A, R_i\rangle < 0\) holds. For k such generators, the conditions \(\langle Q, R_i\rangle < 0\) for \(i=1,\ldots , k\) are not only satisfied for \(Q=A\), but also for all Q in an \(\varepsilon \)-neighborhood of A (with a suitable \(\varepsilon \) depending on k). For any k, this neighborhood also contains points of \({{\,\mathrm{int}\,}}({{\mathcal {V}}}(P)) \subseteq {{\,\mathrm{int}\,}}({\mathcal {CP}}^{n})\). For these interior points Q, however, Algorithm 1 finishes after at most finitely many steps (when checking for \(Q \in {{\mathcal {V}}}(P)\) in Step 2). Thus, for some finite number of suitable choices in Step 2(a), the algorithm also ends for A. \(\square \)

Computing contiguous \({\mathcal {COP}}\)-perfect matrices

Our algorithm for computing contiguous \({\mathcal {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. The following algorithm is similar to [9, Section 6, Erratum to algorithm of Section 2.3] and [31].

figureb

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 \({{\,\mathrm{min_{\mathcal {COP}}}\,}}\) and all vectors \({{\,\mathrm{Min_{\mathcal {COP}}}\,}}\) attaining it. We discuss these tasks in Sects. 3.3 and 3.4.

In the while loop (Step 2) of Algorithm 2, lower and upper bounds l and u for the desired value \(\lambda \) are computed, such that \(P+ l R\) and \(P + u R\) are lying in \({{\,\mathrm{int}\,}}({\mathcal {COP}}^{n})\) satisfying

$$\begin{aligned} {{\,\mathrm{min_{\mathcal {COP}}}\,}}(P+ l R)={{\,\mathrm{min_{\mathcal {COP}}}\,}}(P) > {{\,\mathrm{min_{\mathcal {COP}}}\,}}(P + u R). \end{aligned}$$

In other words, \(P + lR\) lies on the edge \([P,N] \subseteq {\mathcal {R}}\), but \(P + uR\) lies outside of \({\mathcal {R}}\).

The set S in Step 3 contains all vectors \(v \in {\mathbb {Z}}^n_{\ge 0}\) defining a separating hyperplane \(\{X \in {\mathcal {S}}^n : \langle X, vv^\mathsf{T} \rangle = 1\}\), separating \({\mathcal {R}}\) and \(P + uR\).

If \(v \in S\), then \(R[v] < 0\) and

$$\begin{aligned} (P+\lambda R)[v] = P[v] + \min _{w \in S} \left( \frac{1-P[w]}{R[w]}\right) R[v] \ge 1. \end{aligned}$$

If \(v \not \in S\) and \(R[v] \ge 0\), then clearly \((P + \lambda R)[v] \ge 1\), since \(\lambda \ge 0\). Finally, if \(v \not \in S\) and \(R[v] < 0\), then since \(\lambda \le u\), we have

$$\begin{aligned} (P+\lambda R)[v] \ge (P+u R)[v] \ge 1. \end{aligned}$$

Therefore, the choice of \(\lambda \) in Step 4 guarantees that \(P + \lambda R\) is the contiguous \(\mathcal {COP}\)-perfect matrix of P. We have \({{\,\mathrm{min_{\mathcal {COP}}}\,}}(P + \lambda R) = 1\) but \({{\,\mathrm{Min_{\mathcal {COP}}}\,}}(P + \lambda R) \not \subseteq {{\,\mathrm{Min_{\mathcal {COP}}}\,}}(P)\).

In practice the set S in Step 3 is maybe too big for a complete enumeration. In this case partial enumerations may help to pick successively smaller u’s first, which are not necessarily equal to the desired \(\lambda \); see [31].

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 [24].

Nevertheless, in our algorithms we need to check whether or not a given symmetric matrix lies in the cone of copositive matrices (Step 2(c) of Procedure 3) or in its interior (Step 2 of Algorithm 2). This can be checked by the following recursive characterization of Gaddum [15, Theorem 3.1 and 3.2], which of course is not computable in polynomial time: By

$$\begin{aligned} \Delta = \left\{ x \in {\mathbb {R}}^n : x \ge 0, \sum _{i=1}^n x_i = 1\right\} \end{aligned}$$

we denote the \((n-1)\)-dimensional standard simplex in dimension n. A matrix \(B \in {\mathcal {S}}_n\) lies in \({\mathcal {COP}}^{n}\) (in \({{\,\mathrm{int}\,}}({\mathcal {COP}}^{n})\)) if and only if every of its principal minors of size \((n-1)\times (n-1)\) lies in \({\mathcal {COP}}^{n}\) (in \({{\,\mathrm{int}\,}}(\mathcal {COP}_{n-1})\)) and the value

$$\begin{aligned} v = \max _{x \in \Delta } \min _{y \in \Delta } x^\mathsf{T} B y = \min _{y \in \Delta } \max _{x \in \Delta } x^\mathsf{T} B y. \end{aligned}$$
(10)

of the two-player game with payoff matrix B365体育网站 is non-negative (strictly positive).

One can compute the value of v in (10365体育网站) by a linear program:

$$\begin{aligned} v = \max \{\lambda : \lambda \in {\mathbb {R}}, y \in \Delta , By \ge \lambda e\}, \end{aligned}$$

where \(e = (1, \ldots , 1)^\mathsf{T}\) is the all-ones vector.

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 [3] to compute its copositive minimum \({{\,\mathrm{min_{\mathcal {COP}}}\,}}(B)\) and all vectors \({{\,\mathrm{Min_{\mathcal {COP}}}\,}}(B)\) attaining it. Again we note that this is not a polynomial time algorithm.

First we recall some facts and results from [3]. A family \({\mathcal {P}} = \{\Delta ^1, \ldots , \Delta ^m\}\) of simplices is called a simplicial partitioning of the standard simplex \(\Delta \) if

$$\begin{aligned} \Delta = \bigcup _{i=1}^m \Delta ^i \quad \text {with} \quad {{\,\mathrm{int}\,}}(\Delta ^i) \cap {{\,\mathrm{int}\,}}(\Delta ^j) = \emptyset \quad \text {whenever } i \ne j. \end{aligned}$$

Let \(v^k_1, \ldots , v^k_n\) be the vertices of simplex \(\Delta ^k\). It is easy to verify that if a symmetric matrix \(B \in {\mathcal {S}}^n\) satisfies the strict inequalities

$$\begin{aligned} (v^k_i)^\mathsf{T} B v^k_j > 0 \text { for all } i, j = 1, \ldots , n, \text{ and } k= 1, \ldots , m, \end{aligned}$$
(11)

then it lies in \({{\,\mathrm{int}\,}}({\mathcal {COP}}^{n})\). Bundfuss and Dür [3, Theorem 2] proved the following converse: Suppose \(B \in {{\,\mathrm{int}\,}}({\mathcal {COP}}^{n})\), then there exists an \(\varepsilon > 0\) so that for all finite simplex partitions \({\mathcal {P}} = \{\Delta ^1, \ldots , \Delta ^m\}\) of \(\Delta \), where the diameter of every simplex \(\Delta ^k\) is at most \(\varepsilon \), strict inequalities (11) hold. Here, the diameter of \(\Delta ^k\) is defined as \(\max \{\Vert v^k_i - v^k_j\Vert : i, j = 1, \ldots , n\}\).

We assume now that \(B \in {{\,\mathrm{int}\,}}({\mathcal {COP}}^{n})\) and that we have a finite simplex partition \({\mathcal {P}}\) so that (11) 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 [3365体育网站, Algorithm 2].

Each simplex \(\Delta ^k = {{\,\mathrm{conv}\,}}\{v^k_1, \ldots , v^k_n\}\) defines a simplicial cone by \({{\,\mathrm{cone}\,}}\{v^k_1, \ldots , 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\)365体育网站’s to have integral coordinates.

The goal is now to find all integer vectors v in \(\Delta ^k\) which minimize B[v]. To do this we adapt the algorithm of Fincke and Pohst [14], 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 \in {{\,\mathrm{int}\,}}({\mathcal {CP}}^{n})\) and a simplicial cone, which is generated by integer vectors \(v_1, \ldots , v_n\) so that \(v_i^\mathsf{T} B v_j \ge 0\) holds, and given a positive constant M, find all integer vectors v in the cone so that \(B[v] \le 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 B365体育网站 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, \ldots , v_n\) as it columns. (see for example Kannan and Bachem [21] or Schrijver [27], where it is shown that computing the Hermite normal form can be done in polynomial time). We find a unimodular matrix \(U \in \mathsf {GL}_n({\mathbb {Z}})\) such that \(UV = W\) holds, where W is an upper triangular matrix with columns \(w_1, \ldots , w_n\) and coefficients \(W_{i,j}\). Note that the diagonal coefficients of W are not zero since W has full rank. Moreover, denoting the matrix \((U^{-1})^\mathsf{T} B U^{-1}\) by \(B'\) we have for all ij

$$\begin{aligned} 0\le v_i^\mathsf{T} B v_j = w_i^\mathsf{T} (U^{-1})^\mathsf{T} B U^{-1} w_j = w_i B' w_j , \end{aligned}$$
(12)

where the inequality is strict for whenever \(i = j\).

We want to find all vectors \(v \in {{\,\mathrm{cone}\,}}\{v_1, \ldots , v_n\} \cap {\mathbb {Z}}^n\) so that \(B[v] \le M\). In other words, the goal is to find all rational coefficients \(\alpha _1, \ldots , \alpha _n\)365体育网站 satisfying the following three properties:

  1. (i)

    \(\alpha _1, \ldots , \alpha _n \ge 0\),

  2. (ii)

    \(\sum _{i=1}^n \alpha _i v_i \in {\mathbb {Z}}^n\),

  3. (iii)

    \(B\left[ \sum _{i=1}^n \alpha _i v_i\right] \le M\).

Since matrix U lies in \(\mathsf {GL}_n({\mathbb {Z}})\), a vector \(\sum _{i=1}^n \alpha _i v_i\) is integral if and only if \(\sum _{i=1}^n \alpha _i w_i\) is integral. Looking at the last vector componentwise we have

$$\begin{aligned} \sum _{i=1}^n \alpha _i w_i = \left( \sum _{j=1}^n \alpha _j W_{1,j}, \sum _{j=2}^n \alpha _j W_{2,j}, \ldots , \alpha _{n-1} W_{n-1,n-1} + \alpha _n W_{n-1,n}, \alpha _n W_{n,n} \right) . \end{aligned}$$

We first consider the possible values of the last coefficient \(\alpha _n\), and then continue to other coefficients \(\alpha _{n-1}, \ldots , \alpha _1\)365体育网站, one by one via a backtracking search. Conditions (i) and (ii) imply that

$$\begin{aligned} \alpha _n \in \{k/W_{n,n} : k = 0, 1, 2, \ldots \}. \end{aligned}$$

Condition (iii) gives an upper bound for \(\alpha _n\): Write \(\alpha = (\alpha _1, \ldots , \alpha _n)^\mathsf{T}\), then

$$\begin{aligned} M \ge (V\alpha )^\mathsf{T} B V\alpha = \alpha ^\mathsf{T} W^\mathsf{T} B' W \alpha = B'\left[ \sum _{i=1}^n \alpha _i w_i\right] \ge B'[\alpha _n w_n] = \alpha ^2_n B'[w_n], \end{aligned}$$

where the last inequality follows from (12). Hence, \(\alpha _n \le \sqrt{M/B'[w_n]}\) and so

$$\begin{aligned} \alpha _n \in \left\{ k/W_{n,n} : k = 0, 1, \ldots , \left\lfloor \sqrt{M/B'[w_n]} \right\rfloor W_{n,n}\right\} . \end{aligned}$$

Now suppose \(\alpha _n\) is fixed. We want to compute all possible values of the coefficient \(\alpha _{n-1}\). Then the second but last coefficient \(\alpha _{n-1} W_{n-1,n-1} + \alpha _n W_{n-1,n}\) should be integral and \(\alpha _{n-1}\)365体育网站 should be non-negative. Thus,

$$\begin{aligned} \alpha _{n-1} \in \left\{ (k - \alpha _n W_{n-1,n})/W_{n-1,n-1} : k = \lceil \alpha _n W_{n-1,n} \rceil , \lceil \alpha _n W_{n-1,n} \rceil + 1, \ldots \right\} . \end{aligned}$$

Again we use condition (iii) to get an upper bound for \(\alpha _{n-1}\):

$$\begin{aligned} \begin{aligned}&M \ge B'\left[ \sum _{i=1}^n \alpha _i w_i\right] \ge B'[\alpha _{n-1} w_{n-1} + \alpha _n w_n]\\&\qquad = \alpha _{n-1}^2 B'[w_{n-1}] + 2\alpha _{n-1}\alpha _n w_{n-1}^\mathsf{T} B' w_{n} + \alpha _n^2 B'[w_n], \end{aligned} \end{aligned}$$

and solving the corresponding quadratic equation gives the desired upper bound.

Now suppose \(\alpha _n\) and \(\alpha _{n-1}\) are fixed. We want to compute all possible values of the coefficient \(\alpha _{n-2}\)365体育网站 and we can proceed inductively.

Modifying the algorithm for general input

In this section we discuss an adaption of Algorithm 1 for general symmetric matrices A as input. If A is not in \({\mathcal {CP}}^{n}\) then the procedure ends with a separating witness matrix W if it terminates. However, we currently do not know if our Procedure 3 always terminates in this case (cf. Conjecture 3.2).

figurec

The difference between Algorithm 1 and Procedure 3 is in the new steps 2(a) and 2(c). Here it is tested, whether or not we can already certify that the input matrix A is not in \({\mathcal {CP}}^{n}\).

In Step 2(a) we check whether or not the current \({\mathcal {COP}}\)-perfect matrix P is already a separating witness. By this, the algorithm subsequently constructs an outer approximation of the \({\mathcal {CP}}^{n}\) cone:

$$\begin{aligned} {\mathcal {CP}}^{n}\subseteq \{ Q\in {{\mathcal {S}}}^{n}: \langle Q, B \rangle \ge 0 \hbox { for all constructed vertices}\ B \hbox { of } {\mathcal {R}} \} \end{aligned}$$

This procedure gives a tighter and tighter outer approximation of the completely positive cone (cf. Theorem 2.3).

In Step 2(c) it is checked whether or not R is a separating witness for A, that is, if not only \(\langle A,R \rangle < 0\) but also \(R\in {\mathcal {COP}}^{n}\) holds. The copositivity test of R can be realized as explained in Sect. 3.3.

For the case of \(A\not \in {\mathcal {CP}}^{n}\), we do not know if it is possible that Procedure 3 does not provide a separating witness W365体育网站 after finitely many iterations. With a suitably chosen rule in Step 2(b), however, we conjecture that the computation finishes with a certificate:

Conjecture 3.2

For \(A\not \in {\mathcal {CP}}^{n}\), Procedure 3 with a suitable “pivot rule” in Step 2(b) 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 \in {{\,\mathrm{bd}\,}}{\mathcal {CP}}^{n}{\setminus } \widetilde{\mathcal {CP}}^{n}\). In this case, Procedure 3 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 \(P^{(i)}\) of \({{\mathcal {R}}}\) constructed in Procedure 3365体育网站, we know however at least the following:

  1. (i)

    The \({\mathcal {COP}}\)-perfect matrix \(P^{(i)}\) is in \( \{ B\in {\mathcal {COP}}^{n}: \langle P^{(i-1)}, A\rangle > \langle B, A\rangle \ge 0 \}. \)

  2. (ii)

    The norms \(\Vert P^{(i)} \Vert \) are unbounded. Otherwise—following the arguments in the proof of Lemma 2.4 in [10]—we could construct a convergent subsequence with limit \(P \in {{\mathcal {R}}}\), for which we could then find a \(u\in {{\mathbb {R}}}^n_{\ge 0}\) of norm \(\Vert u\Vert =1\) with \(P[u]=0\) (contradicting \(P\in {{\,\mathrm{int}\,}}({\mathcal {COP}}^{n})\)).

  3. (iii)

    \(P^{(i)} / \Vert P^{(i)}\Vert \) contains a convergent subsequence with limit \(P \in \{ X \in {{\mathcal {S}}}^{n}: \langle X, A \rangle = 0\} \). It can be shown that this P is in \({{\,\mathrm{bd}\,}}{\mathcal {COP}}^{n}\). Infinite sequences of vertices \(P^{(i)}\) of \({{\mathcal {R}}}\) with such a limit P exist. For \(n=2\) we give an example in Sect. 4, in which A is from the “irrational boundary part” \(({{\,\mathrm{bd}\,}}{\mathcal {CP}}^{n}) {\setminus } \widetilde{\mathcal {CP}}^{n}\).

A 2-dimensional example

In this section we demonstrate how Algorithm 1 respectively Procedure 3 works for \(n = 2\). Thereby, we discover a relation to beautiful classical results in elementary number theory. In particular, we consider the case when the input matrix A lies on the boundary of \({\mathcal {CP}}^{2}\), see Fig. 1.

Fig. 1
figure1

Subdivision of \({\mathcal {CP}}^{2}\) by Voronoi cones \({\mathcal {V}}(P)\). Matrices \(A = (a_{ij})\) are drawn with 2-dimensional coordinates \((x,y) = \frac{1}{a_{11}+a_{22}} (a_{11}-a_{22}, a_{12}).\) Integers \(\alpha ,\beta \) indicate that the shown point is on a ray spanned by the rank-1 matrix \(A=v v^\mathsf{T}\) with \(v=( \alpha , \beta )^\mathsf{T}\)

Input on the boundary

The boundary of \({\mathcal {CP}}^{2}\)365体育网站 splits into a part of diagonal matrices

$$\begin{aligned} A = \begin{pmatrix} \alpha &{}\quad 0 \\ 0 &{}\quad \beta \end{pmatrix} \quad \text {with} \quad \alpha , \beta \ge 0 \end{aligned}$$

and into rank-1 matrices \(A=xx^\mathsf{T}\). In the first case, Procedure 3 finishes already in its first iteration, if we use \(Q_{{\mathsf {A}}_2}\) as a starting perfect matrix,Footnote 2 where

$$\begin{aligned} Q_{{\mathsf {A}}_2} = \begin{pmatrix} 2 &{}\quad -1 \\ -1 &{}\quad 2 \end{pmatrix} \quad \text {and} \quad {{\,\mathrm{Min_{\mathcal {COP}}}\,}}(Q_{{\mathsf {A}}_2}) = \left\{ \begin{pmatrix} 1 \\ 0 \end{pmatrix} , \begin{pmatrix} 0 \\ 1 \end{pmatrix} , \begin{pmatrix} 1 \\ 1 \end{pmatrix} \right\} . \end{aligned}$$
(13)

Let us consider the other boundary cases for \(n=2\), where \(A = xx^\mathsf{T}\) is a rank-1 matrix. Without loss of generality we can assume that \(x = (\alpha , 1)^\mathsf{T}\). As we explain in the following, Procedure 3 will terminate after finitely many iterations with a \({\mathcal {COP}}\)-perfect matrix P satisfying \(x \in {{\,\mathrm{Min_{\mathcal {COP}}}\,}}P\) when \(\alpha \) is rational. For irrational \(\alpha \)365体育网站 the procedure will not terminate.

The first observation is that Procedure 3 subsequently replaces a \({\mathcal {COP}}\)-perfect matrix P by a contiguous \({\mathcal {COP}}\)-perfect matrix N in a way that one of the three vectors in \({{\,\mathrm{Min_{\mathcal {COP}}}\,}}(P)\) is replaced by the sum of the remaining two. Let P365体育网站 be a copositive matrix with

$$\begin{aligned} {{\,\mathrm{Min_{\mathcal {COP}}}\,}}P = \left\{ \begin{pmatrix} a \\ b \end{pmatrix} , \begin{pmatrix} c \\ d \end{pmatrix} , \begin{pmatrix} e \\ f \end{pmatrix} \right\} . \end{aligned}$$
(14)

It is known (see for example [17, Section “Determinants Determine Edges”]) that \(\det \left( {\begin{matrix} a &{} c\\ b &{} d\end{matrix}}\right) = \pm 1\) and we get a contiguous \({\mathcal {COP}}\)-perfect matrix N with

$$\begin{aligned} {{\,\mathrm{Min_{\mathcal {COP}}}\,}}N = \left\{ \begin{pmatrix} a \\ b \end{pmatrix} , \begin{pmatrix} c \\ d \end{pmatrix} , \begin{pmatrix} a+c \\ b+d \end{pmatrix} \right\} \ne {{\,\mathrm{Min_{\mathcal {COP}}}\,}}P \end{aligned}$$

by

$$\begin{aligned} N = P + 4 \begin{pmatrix} bd &{}\quad -\frac{1}{2}(ad + bc)\\ -\frac{1}{2}(ad + bc) &{}\quad ac \end{pmatrix}. \end{aligned}$$

For instance, starting with \(P = Q_{{\mathsf {A}}_2}\) as in (13)

$$\begin{aligned} \begin{pmatrix} 1 \\ 0 \end{pmatrix} \text { is replaced by } \begin{pmatrix} 1 \\ 2 \end{pmatrix} \text { if } \alpha < 1 \text { yielding } N = \begin{pmatrix} 6 &{}\quad -3\\ -3 &{}\quad 2 \end{pmatrix}, \end{aligned}$$

or

$$\begin{aligned} \begin{pmatrix} 0 \\ 1 \end{pmatrix} \text { is replaced by } \begin{pmatrix} 2 \\ 1 \end{pmatrix} \text { if } \alpha > 1 \text { yielding } N = \begin{pmatrix} 2 &{}\quad -3\\ -3 &{}\quad 6 \end{pmatrix}. \end{aligned}$$

Note that for \(\alpha =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 [19] in 1894 inspired by a lecture of Felix Klein; see also the book by Hatcher [17365体育网站], which contains the proofs.

For concreteness, let us choose \(\alpha = \sqrt{2}\). Then \({{\,\mathrm{Min_{\mathcal {COP}}}\,}}(P)\) is changed by replacing a suitable vector subsequently with

$$\begin{aligned} \begin{pmatrix} 2 \\ 1 \end{pmatrix} , \begin{pmatrix} 3 \\ 2 \end{pmatrix} , \begin{pmatrix} 4 \\ 3 \end{pmatrix} , \begin{pmatrix} 7 \\ 5 \end{pmatrix} , \begin{pmatrix} 10 \\ 7 \end{pmatrix} , \begin{pmatrix} 17 \\ 12 \end{pmatrix} , \begin{pmatrix} 24 \\ 17 \end{pmatrix} , \begin{pmatrix} 41 \\ 29 \end{pmatrix} , \begin{pmatrix} 58 \\ 41 \end{pmatrix} , \begin{pmatrix} 99 \\ 70 \end{pmatrix} , \ldots \end{aligned}$$

Note that there is always a unique choice in Step 2(b) of Procedure 3 in case A is a \(2\times 2\) rank-1 matrix. Note also that the vectors represent fractions that converge to \(\sqrt{2}\). Every second vector corresponds to a convergent of the continued fraction expansion of \(\sqrt{2}\): We have

$$\begin{aligned} \sqrt{2} = 1 + \frac{1}{2 +\frac{1}{2+\frac{1}{2 + \frac{1}{2 + \frac{1}{2 + \ddots }}}}} \end{aligned}$$

and

$$\begin{aligned} 3/2 = 1 + \frac{1}{2}, \; 7/5 = 1 + \frac{1}{2 + \frac{1}{2}},\ldots , \; 99/70 = 1 + \frac{1}{2 +\frac{1}{2+\frac{1}{2 + \frac{1}{2 + \frac{1}{2}}}}},\ldots \end{aligned}$$

The \({\mathcal {COP}}\)-perfect matrix after ten iterations of the algorithm is

$$\begin{aligned} P^{(10)}= \begin{pmatrix} 4756 &{}\quad -6726 \\ -6726 &{}\quad 9512 \end{pmatrix}. \end{aligned}$$

It can be shown that the matrices \(P^{(i)}\)365体育网站 converge to a multiple of

$$\begin{aligned} B = \begin{pmatrix} 1 &{}\quad -\sqrt{2} \\ -\sqrt{2} &{}\quad 2 \end{pmatrix} \text { satisfying } \langle A,B\rangle = 0 \text { and } \langle X,B \rangle \ge 0 \text { for all } X \in {\mathcal {CP}}^{2}. \end{aligned}$$

However, every one of the infinitely many perfect matrices \(P^{(i)}\) satisfies

$$\begin{aligned} \langle X,P^{(i)}\rangle > 0 \text { for all } X \in {\mathcal {CP}}^{2}. \end{aligned}$$

Input outside

In case the input matrix \(A=(a_{ij})\) is outside of \({\mathcal {CP}}^{2}\) we distinguish two cases using the starting \({\mathcal {COP}}\)-perfect matrix \(Q_{{\mathsf {A}}_2}\): If \(a_{12} = a_{21} < 0\) then Procedure 3 finishes already in its first iteration (in Step 2(c)) with a separating witness

$$\begin{aligned} W = R = \begin{pmatrix} 0 &{}\quad 1 \\ 1 &{}\quad 0 \end{pmatrix} . \end{aligned}$$

If \(a_{12} = a_{21} \ge 0\), Procedure 3 terminates after finitely many iterations (in Step 2(a)) with a separating \({\mathcal {COP}}\)-perfect witness matrix \(W = P\).

We additionally note that it is a special feature of the \(n=2\) case that we can conclude that the input matrix A is outside of \({\mathcal {CP}}^{2}\) if we have a choice between two possible R with \(\langle A,R \rangle < 0\) in Step 2(b) of Procedure 3.

Integral input

Laffey and Šimgoc [22] showed that every integral matrix \(A \in {\mathcal {CP}}^{2}\) possesses an integral cp-factorization. This can also be seen as follows: If P is a copositive matrix with \({{\,\mathrm{Min_{\mathcal {COP}}}\,}}P\) as in (14) then the matrices

$$\begin{aligned} \begin{pmatrix} a\\ b \end{pmatrix} \begin{pmatrix} a\\ b \end{pmatrix}^\mathsf{T}, \begin{pmatrix} c\\ d \end{pmatrix} \begin{pmatrix} c\\ d \end{pmatrix}^\mathsf{T}, \begin{pmatrix} e\\ f \end{pmatrix} \begin{pmatrix} e\\ f \end{pmatrix}^\mathsf{T} \end{aligned}$$

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 \(P = Q_{\mathsf{A}_2}\). Then all the other cones are equivalent by conjugating with a matrix in \(\mathsf {GL}_2({\mathbb {Z}})\).

Computational experiments

We implemented our algorithm. The source code, written in C++, is available on GitHub [18]. 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 respectively Procedure 3 drastically varies in the considered examples. Most of the computational time is taken by the computation of the copositive minimum as described in Sect. 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 [12] and in [6] characterizations of matrices in the interior of the completely positive cone are given. For example, we have that \(A \in {{\,\mathrm{int}\,}}({\mathcal {CP}}^{n})\) if and only if A has a factorization \(A = BB^\mathsf{T}\) with \(B > 0\) and \({{\,\mathrm{rank}\,}}B = n\).

The matrix

$$\begin{aligned} \begin{pmatrix} 6 &{}\quad 7 &{}\quad 8 &{}\quad 9 &{}\quad 10 &{}\quad 11\\ 7 &{}\quad 9 &{}\quad 10 &{}\quad 11 &{}\quad 12 &{}\quad 13\\ 8 &{}\quad 10 &{}\quad 12 &{}\quad 13 &{}\quad 14 &{}\quad 15\\ 9 &{}\quad 11 &{}\quad 13 &{}\quad 15 &{}\quad 16 &{}\quad 17\\ 10 &{}\quad 12 &{}\quad 14 &{}\quad 16 &{}\quad 18 &{}\quad 19\\ 11 &{}\quad 13 &{}\quad 15 &{}\quad 17 &{}\quad 19 &{}\quad 21 \end{pmatrix} \end{aligned}$$

for example lies in the interior of \(\mathcal {CP}_6\), as it has a cp-factorization with vectors (1, 1, 1, 1, 1, 1), (1, 1, 1, 1, 1, 2), (1, 1, 1, 1, 2, 2), (1, 1, 1, 2, 2, 2), (1, 1, 2, 2, 2, 2) and (1, 2, 2, 2, 2, 2). It is found after 8 iterations of our algorithm.

Matrices on the boundary

For matrices in \(\widetilde{\mathcal {CP}}^{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 [16] and lies in the boundary of \(\widetilde{\mathcal {CP}}^5\):

$$\begin{aligned} \begin{pmatrix} 8 &{}\quad 5 &{}\quad 1 &{}\quad 1 &{}\quad 5\\ 5 &{}\quad 8 &{}\quad 5 &{}\quad 1 &{}\quad 1\\ 1 &{}\quad 5 &{}\quad 8 &{}\quad 5 &{}\quad 1\\ 1 &{}\quad 1 &{}\quad 5 &{}\quad 8 &{}\quad 5\\ 5 &{}\quad 1 &{}\quad 1 &{}\quad 5 &{}\quad 8\\ \end{pmatrix} \end{aligned}$$

Starting from \(Q_{{{\mathsf {A}}_5}}\)365体育网站 our algorithm needs 5 iterations to find the cp-factorization with the ten vectors (0, 0, 0, 1, 1), (0, 0, 1, 1, 0), (0, 0, 1, 2, 1), (0, 1, 1, 0, 0), (0, 1, 2, 1, 0), (1, 0, 0, 0, 1), (1, 0, 0, 1, 2), (1, 1, 0, 0, 0), (1, 2, 1, 0, 0) and (2, 1, 0, 0, 1).

365体育网站While the above example can be solved within seconds on a standard computer, the matrix

$$\begin{aligned} A = \begin{pmatrix} 41 &{}\quad 43 &{}\quad 80 &{}\quad 56 &{}\quad 50\\ 43 &{}\quad 62 &{}\quad 89 &{}\quad 78 &{}\quad 51\\ 80 &{}\quad 89 &{}\quad 162 &{}\quad 120 &{}\quad 93\\ 56 &{}\quad 78 &{}\quad 120 &{}\quad 104 &{}\quad 62\\ 50 &{}\quad 51 &{}\quad 93 &{}\quad 62 &{}\quad 65 \end{pmatrix} \end{aligned}$$

from Example 7.2 in [16] 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). The second algorithm suggested in [16365体育网站] found the following approximate cp-factorization in 0.018 seconds

$$\begin{aligned} A = {\tilde{B}}{\tilde{B}}^\mathsf{T}, \; \text { with } \; {\tilde{B}} = \begin{pmatrix} 0.0000 &{}\quad 3.3148 &{}\quad 4.3615 &{}\quad 3.3150 &{}\quad 0.0000\\ 0.0000 &{}\quad 0.7261 &{}\quad 4.3485 &{}\quad 6.5241 &{}\quad 0.0000\\ 0.0000 &{}\quad 4.5242 &{}\quad 9.9675 &{}\quad 6.4947 &{}\quad 0.0000\\ 0.0000 &{}\quad 0.1361 &{}\quad 7.4192 &{}\quad 6.9955 &{}\quad 0.0000\\ 0.0000 &{}\quad 5.3301 &{}\quad 3.8960 &{}\quad 4.6272 &{}\quad 0.0000 \end{pmatrix}. \end{aligned}$$

We also considered the following family of completely positive \((n+m)\times (n+m)\) matrices, generalizing the family of examples considered in [20]: The matrices

$$\begin{aligned} \begin{pmatrix} n {{\,\mathrm{Id}\,}}_m &{}\quad J_{m,n}\\ J_{n,m} &{}\quad m {{\,\mathrm{Id}\,}}_n \end{pmatrix}, \end{aligned}$$

with \(J_{\cdot ,\cdot }\) 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 \({\mathcal {COP}}\)-perfect matrix \(Q_{{{\mathsf {A}}_{m+n}}}\) for all \(n,m\le 3\)365体育网站 in less than 6 iterations.

Matrices that are not completely positive

365体育网站For matrices that are not completely positive, our algorithm can find a certificate in form of a witness matrix that is copositive.

The following example is taken from [25, Example 6.2].

$$\begin{aligned} A= \begin{pmatrix} 1 &{}\quad 1 &{}\quad 0 &{}\quad 0 &{}\quad 1\\ 1 &{}\quad 2 &{}\quad 1 &{}\quad 0 &{}\quad 0\\ 0 &{}\quad 1 &{}\quad 2 &{}\quad 1 &{}\quad 0\\ 0 &{}\quad 0 &{}\quad 1 &{}\quad 2 &{}\quad 1\\ 1 &{}\quad 0 &{}\quad 0 &{}\quad 1 &{}\quad 6 \end{pmatrix} \end{aligned}$$

is positive semidefinite, but not completely positive. Starting from \(Q_{{{\mathsf {A}}_5}}\)365体育网站 our algorithm needs 18 iterations to find the copositive witness matrix

$$\begin{aligned} B = \begin{pmatrix} 363/5 &{}\quad -2126/35 &{}\quad 2879/70 &{}\quad 608/21 &{}\quad -4519/210\\ -2126/35 &{}\quad 1787/35 &{}\quad -347/10 &{}\quad 1025/42 &{}\quad 253/14\\ 2879/70 &{}\quad -347/10 &{}\quad 829/35 &{}\quad -1748/105 &{}\quad 371/30\\ 608/21 &{}\quad 1025/42 &{}\quad -1748/105 &{}\quad 1237/105 &{}\quad -601/70\\ -4519/210 &{}\quad 253/14 &{}\quad 371/30 &{}\quad -601/70 &{}\quad 671/105 \end{pmatrix} \end{aligned}$$

with \(\langle A, B \rangle = -2/5\), verifying \(A \not \in \mathcal {CP}_5\).

Notes

  1. 1.

    We use the letter \({\mathcal {R}}\) here because Ryshkov used a similar construction in the study of lattice sphere packings, see for example [28, Chapter 3].

  2. 2.

    Strictly speaking we should use \(\frac{1}{2}Q_{{\mathsf {A}}_2}\) here. If we use \(Q_{{\mathsf {A}}_2}\) instead, then the algorithm produces integral matrices and vertices of \(2{\mathcal {R}}\).

References

  1. 1.

    Berman, A., Rothblum, U.G.: A note on the computation of the CP-rank. Linear Algebra Appl. 419, 1–7 (2006)

  2. 2.

    Berman, A., Shaked-Monderer, N.: Completely positive matrices: real, rational, and integral. Acta Math. Vietnam. 43365体育网站, 629–639 (2018)

  3. 3.

    Bundfuss, S., Dür, M.: Algorithmic copositivity detection by simplicial partition. Linear Algebra Appl. 428, 1511–1523 (2008)

  4. 4.

    Bomze, I.M., Schachinger, W., Uchida, G.: Think co(mpletely)positive! Matrix properties, examples and a clustered bibliography on copositive optimization. J. Global Optim. 5365体育网站(2), 423–445 (2012)

  5. 5.

    365体育网站Conway, J.H., Sloane, N.J.A.: Sphere Packings, Lattices and Groups. Springer, New York (1988)

  6. 6.

    Dickinson, P.J.C.: An improved characterisation of the interior of the completely positive cone. Electron. J. Linear Algebra 20, 723–729 (2010)

  7. 7.

    Dickinson, P.J.C.: The copositive cone, the completely positive cone and their generalisations. Ph.D. Thesis, University of Groningen (2013)

  8. 8.

    Dickinson, P.J.C., Gijben, L.: On the computational complexity of membership problems for the completely positive cone and its dual. Comput. Optim. Appl. 57, 403–415 (2014)

  9. 9.

    Dutour Sikirić, M., Schürmann, A., Vallentin, F.: Classification of eight dimensional perfect forms. Electron. Res. Announ. Am. Math. Soc. 13365体育网站 (2007), 21–32, [math.NT]

  10. 10.

    Dutour Sikirić, M., Schürmann, A., Vallentin, F.: Rational factorizations of completely positive matrices. Linear Algebra Appl. 523, 46–51 (2017)

  11. 11.

    Dür, M.: Copositive programming–a survey. In: Diehl, M., Glineur, F., Jarlebring, E., Michiels, W. (eds.) Recent Advances in Optimization and its Applications in Engineering, pp. 3–20. Springer, Berlin (2010)

  12. 12.

    Dür, M., Still, G.: Interior points of the completely positive cone. Electron. J. Linear Algebra 17365体育网站, 48–53 (2008)

  13. 13.

    Elser, V.: Matrix product constraints by projection methods. J. Global Optim. 68, 329–355 (2017)

  14. 14.

    Fincke, U., Pohst, M.: Improved methods for calculating vectors of short length in a lattice, including a complexity analysis. Math. Comput. 44, 463–471 (1985)

  15. 15.

    Gaddum, J.W.: Linear inequalities and quadratic forms. Pacific J. Math. 8365体育网站, 411–414 (1958)

  16. 16.

    Groetzner, P., Dür, M.: A factorization method for completely positive matrices. Preprint (2018)

  17. 17.

    365体育网站Hatcher, A.: Topology of Numbers. Book in preparation (2017). . Accessed 24 Jan 2020

  18. 18.

    Dutour Sikirić, M.: (2018). Copositive. . Accessed 24 Jan 2020

  19. 19.

    Hurwitz, A.: Über die Reduktion der binären quadratischen Formen. Math. Annalen 45, 85–117 (1894)

  20. 20.

    Jarre, F., Schmallowsky, K.: On the computation of \(C^\ast \) certificates. J. Global Optim. 45, 281–296 (2009)

  21. 21.

    Kannan, R., Bachem, A.: Polynomial algorithms for computing the Smith and Hermite normal forms of an integer matrix. SIAM J. Comput. 8365体育网站, 499–507 (1979)

  22. 22.

    Laffey, T., Šmigoc, H.: Integer completely positive matrices of order two. [math.OC]

  23. 23.

    365体育网站Martinet, J.: Perfect Lattices in Euclidean spaces. Springer, Berlin (2003)

  24. 24.

    Murty, K.G., Kabadi, S.N.: Some NP-complete problems in quadratic and nonlinear programming. Math. Program. 39, 117–129 (1987)

  25. 25.

    Nie, J.: The \({\cal{A}}\)-truncated \(K\)-moment problem. Found. Comput. Math. 14, 1243–1276 (2014)

  26. 26.

    Opgenorth, J.: Dual cones and the Voronoi algorithm. Exp. Math. 10365体育网站, 599–608 (2001)

  27. 27.

    Schrijver, A.: Theory of Linear and Integer Programming. Wiley, New York (1986)

  28. 28.

    365体育网站Schürmann, A.: Computational Geometry of Positive Definite Quadratic Forms. American Mathematical Society, Providence (2009)

  29. 29.

    Shaked-Monderer, N., Dür, M., Berman, A.: Complete positivity over the rationals. Pure Appl. Funct. Anal. 3365体育网站, 681–691 (2018)

  30. 30.

    Sponsel, J., Dür, M.: Factorization and cutting planes for completely positive matrices by copositive projection. Math. Program. Ser. A 143, 211–229 (2014)

  31. 31.

    van Woerden, W.P.J.: Perfect quadratic forms: an upper bound and challenges in enumeration. Master thesis, Leiden University (2018)

  32. 32.

    Yıldırım, E.A.: On the accuracy of uniform polyhedral approximations of the copositive cone. Optim. Methods Softw. 27, 155–173 (2012)

Download references

Acknowledgements

365体育网站Open Access funding provided by Projekt DEAL. We like to thank Jeff Lagarias for a helpful discussion about Farey sequences and Renaud Coulangeon for pointing out the link to the work of Opgenorth. We also like to thank Peter Dickinson and Veit Elser for helpful suggestions. Finally, we are very grateful to the referees for carefully reading our paper and for their detailed comments. This helped us to improve the presentation of our paper substantially. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie agreement No 764759 and it is based upon work supported by the National Science Foundation under Grant No. DMS-1439786 while the authors were in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the Spring 2018 semester. Moreover, the first author was supported by the Humboldt foundation and the third author is partially supported by the SFB/TRR 191 “Symplectic Structures in Geometry, Algebra and Dynamics”, funded by the DFG.

Author information

Correspondence to Frank Vallentin.

Additional information

Publisher's Note

365体育网站Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

About this article

Cite this article

Dutour Sikirić, M., Schürmann, A. & Vallentin, F. A simplex algorithm for rational cp-factorization. Math. Program.365体育网站 (2020). http://doi.org/10.1007/s10107-020-01467-4

Download citation

  • Received:

  • Accepted:

  • Published:

Keywords

  • Copositive programming
  • Complete positivity
  • Matrix factorization
  • Copositive minimum

Mathematics Subject Classification

  • 90C20
  • 11H50
  • 11H55