A modified simplex partition algorithm to test copositivity

A real symmetric matrix A is copositive if x⊤Ax≥0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x^\top Ax\ge 0$$\end{document} for all x≥0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x\ge 0$$\end{document}. As A is copositive if and only if it is copositive on the standard simplex, algorithms to determine copositivity, such as those in Sponsel et al. (J Glob Optim 52:537–551, 2012) and Tanaka and Yoshise (Pac J Optim 11:101–120, 2015), are based upon the creation of increasingly fine simplicial partitions of simplices, testing for copositivity on each. We present a variant that decomposes a simplex △\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bigtriangleup $$\end{document}, say with n vertices, into a simplex △1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bigtriangleup _1$$\end{document} and a polyhedron Ω1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega _1$$\end{document}; and then partitions Ω1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega _1$$\end{document} into a set of at most (n-1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(n-1)$$\end{document} simplices. We show that if A is copositive on Ω1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega _1$$\end{document} then A is copositive on △1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bigtriangleup _1$$\end{document}, allowing us to remove △1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bigtriangleup _1$$\end{document} from further consideration. Numerical results from examples that arise from the maximum clique problem show a significant reduction in the time needed to establish copositivity of matrices.


Introduction
Let A be a real symmetric matrix of order n. If x Ax ≥ 0 for all x ∈ R n + = {x ∈ R n | x ≥ 0} then A is called copositive [7]. If x Ax > 0 for all x ∈ R n + \ {0} then A is strictly copositive. The set of all order n copositive matrices COP is a closed, convex, and full-dimensional pointed cone whose interior is the set of strictly copositive matrices. Methods to determine copositivity are important because of their role in algorithms for the solution of copositive optimization problems [2,5], a class of problems that include, for example, the maximum clique problem [1]. The challenge is that Murty and Kabadi [8] have shown that given a square integer matrix A the decision problem "is A not copositive?" is NP-complete.
This algorithm has the same branch-and-bound design as Bundfuss and Dür's [4] BD algorithm, which was improved by Sponsel et al. [9], and then generalized by Tanaka and Yoshise's [10] TY algorithm. All four methods depend on partitions of the standard simplex. Ours is uniquely different in that we decompose a simplex into a simplex 1 and a polyhedron Ω 1 , and then partition Ω 1 into a set of at most (n − 1) simplices, each of which is subsequently decomposed into a simplex and a polyhedron, with the process continuing until copositivity is determined. We will show that if A is copositive on Ω 1 , then A is copositive on 1 , allowing us to remove 1 from further consideration. As our selection of 1 is done in a maximal way, the number of required partitions is, hopefully, significantly reduced.
Similar to the algorithms in [4,9,10], our new algorithm determines copositivity of A by checking whether M AM ∈ M, where M is a matrix whose columns are the vertices of , and where M ⊂ COP. In [4], M = N , the set of all nonnegative n × n matrices. Sponsel et al. [9] comment that while checking M AM ∈ N is relatively easy, the set N is too small and they suggest sets M ⊂ COP with M ⊃ N . While this results in fewer iterations (simplices) the iterations are more time consuming as they require either the calculation of eigenvectors or the solution of semidefinite programs. They reported that in all instances, the BD algorithm with M = N was fastest. It is to the BD algorithm that we compare the algorithm introduced in this paper where we also use M = N .
Tanaka and Yoshise's [10] generalization is based on more sophisticated selections for M. They consider M = N + S + , where S + is the set of all positive semidefinite matrices. The consequence is that a linear program, rather than a semidefinite program, is to be solved to determine whether M AM ∈ M. However, the formulation of the linear program requires a singular value decomposition, and the decomposition is not unique even though the result depends on that decomposition. Tanaka and Yoshise suggest that the TY algorithm is promising for the determination of upper bounds for the maximum clique problem.
Our limited numerical experiments show that for the instances examined there is a significant reduction in the time needed to establish copositivity of matrices that arise in the determination of the clique number of an undirected graph. Further, the experimental results show that for the instances examined there is an advantage to our algorithm when A is near the copositive cone boundary.

Simplicial partition algorithms
We begin with the standard simplex where · 1 is the 1-norm. The matrix A is copositive if and only if it is copositive on s , that is, x Ax ≥ 0 for all x ∈ s . Let P be a simplicial partition of s as defined in [4]. That is, s is the union of all simplices in P; and the intersection of the interior of any two simplices in P is empty. For any ∈ P denote by V the vertices v 1 , . . . , v n of . Thus, the columns of M are the elements of V . Lemma 1, which follows from the work in [9], gives conditions for copositivity. In fact, it shows that if M AM ∈ N , then A is copositive on ∈ P. Lemma 2, which follows from the work in [4], gives a condition for non-copositivity; and, if non-copositivity is not determined, it establishes the existence of a vertex v ∈ V with v Av > 0, needed for the new partitioning process presented in Sect. 3.

Lemma 1
If M AM ≥ 0, i.e., M AM ∈ N , then A is copositive on ∈ P. If A is copositive on all ∈ P then A is copositive on s , and is therefore copositive.
Theorem 1, which is Lemma 2.3 in [9], gives a criterion for termination. Suppose that it has been determined that A is copositive on all ∈P ⊂ P. LetP = P\P and define the diameter ofP by

Theorem 1 If the real symmetric matrix A is not copositive, then there exists an
The BD and TY simplicial partitioning algorithms are branch and bound algorithms that proceed as follows. Start with a partition P of s . This is the first branch step. If Lemma 2 determines that A is not copositive then the algorithm terminates.
If Lemma 1 determines that A is copositive on ∈ P, then can be removed from further consideration. That is, we do not need to create a partition of . This is the bound step.
The branch step is to find simplicial partitions for each of the simplices in the original partition that have not been removed in the bound step. That is, the partitions are further refined. The process continues until it is either determined that A is not copositive (v Av < 0 for some simplex vertex v), that A is copositive on all remaining simplices, in which case A is copositive, or until the remaining simplices are sufficiently small so that the copositivity of A can be determined by the vertices of the remaining simplices.

A new method of partitions
Our new method differs in how the simplicial partitions are created. Consider a simplex ∈ P. If M ˆ AMˆ ≥ 0 then, by Lemma 1, A is copositive onˆ and we can then consider copositivity on the next ∈ P. Let D(M ) be the vector of diagonal elements of M AM . If D(M ) has at least one negative element, then A is not copositive and we are done. If M AM 0 and if D(M ) ≥ 0, then M AM has a negative off-diagonal, say the entry in row i and column j. If both the i-th and j-th diagonal entries are zero, then, by Lemma 2, A is not copositive and we are done. Otherwise, there is at least one strictly positive element in D(M ), say, v Av , that we will use to create a simplex 1 which can be sliced off and eliminated from further consideration.
Let 1 be the simplex with vertices v andv j , 1 ≤ j ≤ n, j = , where thev j are as described in Theorem 2 using the largest value possible for λ j . We have the partition where the convex hull of the points (V 1 ∪ V ) \ {v }. In Fig. 1 we see a slice of a simplex in R 4 and the partition into 1 and Ω 1 that it creates. In this example, = 1. Theorem 3 shows that we can remove 1 from further consideration as copositivity on Ω 1 implies copositivity on 1 .
where v is as described in Theorem 2, we have x Ax > 0. For every x ∈ 1 \v there exist scalars α j ≥ 0, j = 1, . . . , n, summing to unity, such that In this case, x = v so that α = 1. Since the scalars α j /(1 − α ) are nonnegative and sum to unity it follows thatv ∈ Ω 1 . By assumption, A is copositive on Ω 1 so thatv Av ≥ 0. Since and A is copositive on 1 .
If Ω 1 has n vertices then it is a simplex from which we either determine copositivity or partition Ω 1 into a simplex and a polyhedron as above. If Ω 1 has more than n vertices we create the partition Ω 1 = 2 ∪ Ω 2 as follows. Let We then choose any index μ ∈ J 2 . In Fig. 2 we begin with the Ω 1 from Fig. 1. We have J 1 = ∅, J 2 = {2, 3, 4} and we selected μ = 2 giving the simplex 2 with vertices v 2 ,v 2 ,v 3 andv 4 . Lemma 3 shows that these vectors are linearly independent.

Lemma 3
The vectorsv j , 1 ≤ j ≤ n, j = together with v μ are linearly independent.

Proof Consider
Since the vectors in V are linearly independent, and since 0 < λ j < 1, for all j ∈ J 2 , we can conclude that α j = 0, for all j = and that α μ = 0.
By Lemma 3 we conclude that the convex hull of the vectorsv j , 1 ≤ j ≤ n, j = together with v μ is a simplex which we denote by 2 . We denote by Ω 2 the convex hull of the set of vectors V Ω 1 \ {v μ }, which is a polyhedron with one vertex less than Ω 1 . The proposed algorithm will then slice a simplex off Ω 2 creating 3 and Ω 3 , with Ω 3 having one less vertex than Ω 2 , continuing in this way until the Ω polyhedron is itself a simplex. The set of simplices created by this process is denoted by SL( ), where SL represents the slicing process. Since 1 need not be considered further, it is not included in SL( ). The algorithm would continue by restarting the process for each simplex in SL( ). In Fig. 3 we see a graphic of this process with the p ≤ n simplices, labeled 2 to p+1 , in SL( ) down the right-hand side. This, of course, depends on Ω 2 and 2 being a partition of Ω 1 . This is shown in Lemma 4 by showing that the hyperplane passing through v μ ,v j , ∀ j = μ, is a separating hyperplane for 2 and Ω 2 . Denote the hyperplane by H and represent it by w x = b.
Without loss of generality, assume that w v μ < b. Since v μ ,v j , ∀ j = μ, are on H , and sincev μ , v μ ,v j , ∀ j = μ, are linearly independent, it follows that w

Lemma 4
For all x ∈ 2 \ H , we have w x < b and for all x ∈ Ω 2 \ H , we have w x > b.
Proof By assumption we have w x < b when x =v μ . For x ∈ 2 \ H , x =v μ , we can write x as the convex combination Since w v j = b for all j = , which includes j = μ, we have As we have a convex combination it follows that since t ≥ 0 and w v μ < b (by assumption). That w x > b for all x ∈ Ω 2 \ H is proved in a similar manner.

The algorithm
We now present the algorithm. We use P to represent the set of simplices that have not yet been fathomed. For any ∈ P the center and width of are given bȳ respectively. Letδ where ||A|| 1 is the element-wise norm that is the sum of the absolute values of the entries in A.

Algorithm 1: A Modified Simplex Partition Algorithm
Data: An n × n symmetric matrix A. Result: The message that "A is copositive" or that "A is not copositive". Intuitively, we expect δ(P), the diameter of P, to converge to zero, that is, that the sequence of partitions is exhaustive. In [6], Dickinson tells us that whether a sequence of partitions is exhaustive can "often defy our intuition". Fortunately, he provides two methods, that can be adapted to our method, to ensure exhaustion.
Under the assumption that δ(P) converges to zero, if A is strictly copositive or A is not copositive, then the algorithm terminates, but if A is copositive but not strictly copositive, it may or may not terminate [4]. As in [4], we could modify the algorithm by choosing a small tolerance > 0 and terminate when M AM > − for all . In this case, we say that A is -copositive and Theorem 6 in [4] guarantees the termination of the algorithm.

Examples
The first three examples come from graph theory. Let E be the matrix of ones and let A G be the adjacency matrix of a given undirected graph G and define B γ = γ (E − A G ) − E. From [1] we know that the clique number of G is Consider the graphs G 8 , G 10 and G 12 shown in Fig. 4 by way of definition. The figures for G 8 and G 12 are from [9]. We have w(G 8 ) = 3, w(G 10 ) = 3 and w(G 12 ) = 4.
We will determine the copositivity of A = B γ for various values of γ . For the remaining nine examples we check copositivity of A, where A is created as follows. We choose the order n and a γ > 2. ChooseÂ to be symmetric with zeros on the diagonal and Fig. 4 The graphs G 8 , G 10 and G 12 with the upper off-diagonal elements to be zero or one with equal probability. We can think of A as being the adjacency matrix of a randomly generated graph.
We test the proposed algorithm, called SNC, against the BD algorithm. Both algorithms were programmed by the second author and were run on a 2.70 GHz Dual Core Processor machine with 4 GB of RAM. Table 1 shows the results for the graph G 8 . We see that both algorithms correctly determine that w(G 8 ) = 3. For every value of γ , the time taken by SNC is at least an order of magnitude less than for BD. The number of simplices considered by SNC is at most half the number considered by BD. Also, when γ = w(G 8 ), BD failed while SNC correctly determined copositivity. The time taken to determine non-copositivity is two orders of magnitude less for SNC and the number simplices considered an order of magnitude less. Table 2 shows the results for the graph G 10 . We see that both algorithms correctly determine that w(G 10 ) = 3. Both algorithms are able to determine copositivity when γ = w(G 10 ). The number of simplices considered by SNC is an order of magnitude less than for BD when γ > 3.2, and considerably less than BD for γ ≤ 3. Table 3 shows the results for the graph G 12 . For γ ≥ 4 = w(G 12 ) BD is unable to determine copositivity. While SNC is successful, the number of simplices considered is at least three orders of magnitude larger than for G 8 and G 10 . For γ < 4, SNC outperforms BD.
For the three graph theory (clique number) examples, SNC outperforms BD, especially in the determination of non-copositivity.
In Table 4 we see the results for the remaining nine examples. In all cases SNC takes less time and considers fewer simplices than BD.

Conclusion
We have presented a novel method to decompose the standard simplex in order to determine copositivity of real symmetric matrices. Thus, every algorithm that is based on simplicial partitioning in order to test for copositivity can be modified using our method. This includes the algorithms in, for example, [3,4,[9][10][11], and others.
At each major iteration of our algorithm we are able to remove a maximal simplex from future consideration. After removal of the simplex, we are left with a polyhedron, which is then decomposed into at most n simplices. This gives a significant reduction to the number of simplices that need to be considered, and possibly decomposed, making the algorithm run more quickly.  Table 3 Results for the graph G            Table 4 Results for random graphs Our algorithm works well with the set N so that we need not consider more sophisticated choices for M ⊂ COP and hence avoid the need to solve a semidefinite programming problem or linear program. These claims, of fewer simplices and less time, are evidenced by the test results.
The purpose of our limited numerical testing is to validate the new method. As max-clique instances have a very particular structure they are perhaps not sufficient to fully establish a copositivity detection algorithm. Full numerical testing should include instances from the DIMACS challenge as in [4,9]. Rigorous numerical testing and application to other algorithms is the next step.