

Department of Social Systems and Management
Discussion Paper Series

No.1320

An LP-based Algorithm to Test Copositivity
byAkihiro TANAKA and Akiko YOSHISE
May 2014

UNIVERSITY OF TSUKUBA
Tsukuba, Ibaraki 305-8573
JAPANAn LP-based Algorithm to Test Copositivity
Akihiro Tanaka∗ and Akiko Yoshise†
April 2014
Revised October 2014

Abstract
A symmetric matrix is called copositive if it generates a quadratic form taking no negativevalues over the nonnegative orthant, and the linear optimization problem over the set of copositive
matrices is called the copositive programming problem. Recently, many studies have been doneon the copositive programming problem (see, for example, [14, 5]). Among others, several branch
and bound type algorithms have been provided to test copositivity in the context of the fact that
deciding whether a given matrix is copositive is co-NP-complete [23, 13]. In this paper, we propose
a new branch and bound type algorithm for this testing problem based on Sponsel, Bundfuss andDür’s algorithm [27]. Two features of our algorithm are: (1) we introduce new classes of matrices!ns which are relatively large subsets of the set of copositive matrices and work well to
Gns and G!nscheck copositivity of a given n × n symmetric matrix, and (2) for incorporating the sets Gns or G
in checking copositivity, we only have to solve a linear optimization problem with n + 1 variables
and O(n2 ) constraints after computing a singular value matrix decomposition, which implies that
our algorithm is not so time-consuming. Our preliminary numerical experiments suggest that ouralgorithm is promising for determining upper bounds of the maximum clique problem.

Key words. Copositive programming, Matrix decomposition, Linear programming, Branch and boundalgorithm, Maximum clique problem

1

Introduction

Recently, many studies have been done on the conic programming problem of the form(P)!c, x"
Ax = b, x ∈ K

Minimizesubject towhere K ⊂ Rn is a proper cone, i.e., K is a pointed closed convex cone with nonempty interior; !·, ·" :Rn × Rn → R is an inner product over Rn ; A : Rn → Rm is a linear operator and b ∈ Rm . The dual
problem (D) of (P) is given by(D) Maximize
subject to!b, y"
c − A∗ y ∈ K∗

∗ Graduate School of Systems and Information Engineering, University of Tsukuba, Tsukuba, Ibaraki 305-8573, Japan.
e-mail: tanaka.akihiro@sk.tsukuba.ac.jp† Corresponding author. Faculty of Engineering, Information and Systems, University of Tsukuba, Tsukuba, Ibaraki
305-8573, Japan. e-mail: yoshise@sk.tsukuba.ac.jp Research supported in part by Grants-in-Aid for Scientific Research((B)23310099) of the Ministry of Education, Culture, Sports, Science and Technology of Japan.

1where K∗ is the dual cone of K defined byK∗ := {x ∈ Rn | !x, y" ≥ 0 for all y ∈ K}
and !·, ·" : Rm × Rm → R is an inner product over Rm and A∗ is the adjoint operator of A having thefollowing relationship!Ax, y" = !x, A∗ y".
Typical examples of the proper cone K are the n-dimensional nonnegative orthantRn+ := {x ∈ Rn | x ≥ 0}
in linear programming and the positive semidefinite cone
Sn+ := {X ∈ Sn | dT Xd ≥ 0 for all x ∈ Rn },with the set Sn of n × n symmetric matrices in semidefinite programming.
More recently, the following cones are attracting a lot of attention in a context of the relationship between
combinatorial optimization and conic optimization (see, for example, [14, 5]).
- the nonnegative cone Nn := {X ∈ Sn | xij ≥ 0 for all i, j ∈ {1, 2, . . . , n}},!
"- the copositive cone COP n := X ∈ Sn | dT Xd ≥ 0 for all d ∈ Rn+ ,- the Minkowski sum Sn+ + Nn of Sn+ and Nn ,- the doubly nonnegative cone Sn+ ∩ Nn , i.e., the set of positive semidefinite and componentwise nonnegative matrices,#!
"$- the completely positive cone CP n := conv xxT | x ∈ Rn+ where conv (S) denotes the convex hullof the set S.
All of the above cones are proper (see Section 1.6 of [2] where the proper cone is called a full cone), andwe can easily see from the definitions that the following inclusions hold:
COP n ⊇ Sn+ ⊇ Sn+ ∩ Nn ⊇ CP n .

(1)

It is known that the following proposition holds by defining an inner product between X and Y as
!X, Y " := Tr (Y T X).

(2)

Proposition 1.1 (Properties of the copositive cone). (i) The dual cone of the copositive cone COP n
with respect to the inner product (2) is the completely positive cone CP n and vice versa (see p.57
of [1] and Theorem 2.3 of [2]).
(ii) If n ≤ 4 then COP n = Sn+ + Nn (see [11] and Proposition 1.23 of [2]).
(iii) The dual cone of the doubly nonnegative cone Sn+ ∩ Nn with respect to the inner product (2) is the
Minkowski sum Sn+ + Nn of the positive semidefinite cone Sn+ and the nonnegative cone Nn andvice versa (see Remark 1.2).
Remark 1.2. Proposition 1.1, (iii): The equality (Sn+ ∩ Nn )∗ = cl (Sn+ + Nn ) follows from a well-known
result that (K1 ∩ K2 )∗ = cl (K1 + K2 ) holds for any closed convex cones K1 and K2 (see, e.g., p.11 of[15] or Corollary 2.2 of [1]. The closedness of the set Sn+ + Nn follows from a result in [26]. See also
Proposition 4.1 of [29] where the authors showed the property in a little more general framework.
2

The following inclusions follow from (1) and the above proposition
COP n ⊇ Sn+ + Nn ⊇ Sn+ ⊇ Sn+ ∩ Nn ⊇ CP n(3)

and specially, if n ≤ 4 then we have
COP n = Sn+ + Nn ⊇ Sn+ ⊇ Sn+ ∩ Nn = CP n .

(4)Note that the four cones, COP n , CP n , Sn+ ∩ Nn and Sn+ + Nn lack the self-duality and hence are notsymmetric. Since about 2000, there have been many studies conducted on the above four cones as a newresearch direction in the field of conic optimization [3, 4, 10, 28, 24, 7, 25, 9, 20, 8, 29, 27, 21], and they
are called studies on copositive programming [3].
A growing research interest in the field is to provide efficient algorithms to determine whether a given
matrix belongs to COP n (or CP n , or Sn+ + Nn ). It is known that the problem of testing copositivity,
i.e., deciding A ∈ COP n or not, is co-NP-complete [23, 13]. Bomze and Eichfelder [6] have pointed out
what are desirable algorithms for copositivity detection as follows:
However, there are but a few implemented numerical algorithms which (a) apply to generalsymmetric matrices without any structural assumptions or dimensional restrictions; (b) are
not merely recursive, i.e., do not rely on information taken from all principal submatrices,
but rather focus on generating subproblems in a somehow data-driven way.
After citing the paper [7] as an example satisfying criteria (a) and (b), they have presented their new tests
based upon difference-of-convex (d.c.) decompositions, and have combined them to a branch-and-boundalgorithm of ω-subdivision type employing LP or convex QP techniques.
In this paper, we propose a new branch and bound type algorithm based on Sponsel, Bundfuss andDür’s algorithm [27] which is a generalization of the algorithm in [7]. Two features of our algorithm are%ns which are relatively large subsets of the set of
1. we introduce new classes of matrices Gns and G
copositive matrices and work well to check copositivity of a given n × n symmetric matrix, and%ns in checking copositivity, we only have to solve a linear optimiza2. for incorporating the sets Gns or Gtion problem with n + 1 variables and O(n2 ) constraints after computing a singular value matrix
decomposition, which implies that our algorithm is not so time-consuming.
It should be noted that our algorithm is close to the one in [6] in the sense that its key elements are adecomposition of the quadratic term dT Xd, simplicial partitions and LP techniques, while the numberof constraints of the LP in [6] is smaller and O(n). As we will describe in Section 4, the maximum clique
problem can be solved by checking the copositivity of certain matrices. Our preliminary numericalexperiments suggest that our algorithm and its improved versions are promising for determining upperbounds of the maximum clique problem while Bomze and Eichfelder reported effectiveness of theiralgorithm for determining its lower bounds[6].
This paper is organized as follows. In Section 2, we introduce Sponsel, Bundfuss and Dür’s algorithm[27]for checking copositivity of a given matrix and summarize its theoretical results. Based on their algorithm, we propose our algorithm by employing new classes of matrices which are relatively large subsets
of COP n . We also derive some properties of the classes in Section 3 and derive an improvement of ouralgorithm using them. Numerical results are shown in Section 4 and some refinement strategies are
discussed in Section 5. Section 6 gives concluding remarks.
3

2

Sponsel, Bundfuss and Dür’s algorithm to test copositivity

Our algorithm is based on Sponsel, Bundfuss and Dür’s algorithm to test copositivity[27]. In what
follows, we introduce their arguments.
Defining the standard simplex ∆S by ∆S = {x ∈ Rn+ | ,x,1 = 1}, it can be seen that a given n × nsymmetric matrix A is copositive if and only if
xT Ax ≥ 0 for all x ∈ ∆S(see Lemma 1 of [7]). A family of simplices P = {∆1 , . . . , ∆m } is called a simplicial partition of ∆ if itsatisfies
∆=

m
&

i=1

∆i and int(∆i ) ∩ int(∆j ) -= ∅ for all i -= j.

Such a partition can be generated by successively bisecting simplices in the partition. For a given
simplex ∆ = conv{v1 , . . . , vn }, consider the midpoint vn+1 = 12 (vi + vj ) of the edge [vi , vj ]. Thenthe subdivision ∆1 = {v1 , . . . , vi−1 , vn+1 , vi+1 , . . . , vn } and ∆2 = {v1 , . . . , vj−1 , vn+1 , vj+1 , . . . , vn } of∆ satisfies the above conditions for simplicial partitions. See [18] for a more detailed description of
simplicial partitions.
Denote the set of vertices of partition P by
V (P) = {v | v is a vertex of some ∆ ∈ P}.
Each simplex ∆ is determined by its vertices and can be represented by a matrix V∆ whose columns
are these vertices. Note that V∆ is nonsingular and unique up to a permutation of its columns which isirrelevant in the arguments [27]. Define the set of all matrices corresponding to simplices in partition Pas
M (P) = {V∆ : ∆ ∈ P}.
The “fineness” of a partition P is quantified by the maximum diameter of a simplex in P denoted by
δ(P) = max max ||u − v||.
∆∈P u,v∈∆(5)

Using the above notation, the following results on necessary and sufficient conditions for copositivity
have been shown in [27]. The first theorem gives a sufficient condition for copositivity.
Theorem 2.1 (Theorem 2.1 of [27]). If A ∈ Sn satisfies
V T AV ∈ COP n for all V ∈ M (P)then A is copositive. Hence, for any Mn ⊆ COP n , if A ∈ S n satisfies
V T AV ∈ Mn for all V ∈ M (P)then A is also copositive.
The above theorem implies that by choosing Mn = Nn (see (3)), if V∆T AV∆ ∈ Nn holds for any ∆ ∈ Pthen we find that A is copositive.

4

Before describing further conditions for copositivity, we introduce the definition of strict copositivity.
We say that A ∈ Sn is strictly copositive if it satisfies
xT Ax > 0 for all x ∈ Rn+ \ {0}.
It is well-known (and follows from Proposition 1.24 of [2]) that A ∈ Sn is strictly copositive if and only if
A ∈ int (COP n ). Combining this with Theorem 2.2 of [27], we obtain the following necessary condition
for strict copositivity.
Theorem 2.2 (Theorem 2.2 of [27]). Let A ∈ Sn be strictly copositive, i.e., A ∈ int (COP n ). Then
there exists ε > 0 such that for all partitions P of ∆S with δ(P) < ε we have
V T AV ∈ Nn for all V ∈ M (P).
The above theorem ensures that if A is strictly copositive (i.e., A ∈ int (COP n )) then the copositivityof A (i.e., A ∈ COP n ) can be detected in finitely many steps by an algorithm employing a subdivision
rule with δ(P) → 0. A similar result can be obtained for the case A -∈ COP n by the following lemma.
Lemma 2.3 (Lemma 2.3 of [27]). The following two statements are equivalent.
1. A ∈/ COP n2. There exists an ε > 0 such that for any partition P with δ(P) < ε there exists a vertex v ∈ V (P)such that v T Av < 0ɽ
Based on the above three results, the following algorithm has been provided by Sponsel, Bundfuss andDür [27].
Algorithm 1 Sponsel, Bundfuss and Dür’s algorithm to test copositivity
Input: A ∈ Sn , Mn ⊆ COP n
Output: “A is copositive” or “A is not copositive”
1: P ← {∆S };2: while P -= ∅ do3:
Choose ∆ ∈ P;
4:
if v T Av < 0 for some v ∈ V ({∆}): then5:
return “A is not copositive”;
6:
end if7:
if V∆T AV∆ ∈ Mn then8:P ← P \ {∆};
9:
else
10:
partition ∆ into ∆ = ∆1 ∪ ∆2 ;11:
P ← P \ {∆} ∪ {∆1 , ∆2 };12:
end if
13: end while
14: return “A is copositive”;
As we have already observed, Theorem 2.2 and Lemma 2.3 imply the following corollary.
Corollary 2.4.
1. If A is strictly copositive, i.e., A ∈ int (COP n ) then Algorithm 1 terminates
finitely returning “A is copositive.”
2. If A is not copositive, i.e., A -∈ COP n then Algorithm 1 terminates finitely returning “A is notcopositive.”
5

At Line 8 of Algorithm 1, the algorithm removes the simplex which is determined at Line 7 to be in
no need of further exploration by Theorem 2.1. The accuracy and speed of the determination influence
the total computational time and depend on the choice of the set Mn ⊆ COP n . In the next section, we
introduce three examples of Mn , the set used in [27], our alternative suggestion and its generalization.

3

How to choose Mn to efficiently remove unnecessary simplices

In view of Algorithm 1, a desirable set Mn used at Line 7 would have the following properties.
P1 For any given n × n symmetric matrix A ∈ Sn , we can easily check whether A ∈ Mn , andP2 Mn is a subset of the copositive cone COP n as large as possible.
If we choose the nonnegative cone Nn as the set Mn , we can easily check whether A ∈ Mn or not,but the set Nn is too small a subset of COP n and it may take a long time to check the copositivity by
Algorithm 1. In fact, Mn = Nn was used in [7]. On the other hand, the set Sn+ + Nn is a rather large
subset of COP n , but it is not so easy to check whether A ∈ Mn = Sn+ + Nn or not; a well-known way is
to solve the following doubly nonnegative program (which can be expressed as a semidefinite program)
Minimizesubject to!A, X"!I, X" = 1, X ∈ Sn+ ∩ Nn

but solving the problem takes an awful lot of time [27, 29].
Observing these facts, a new alternate of Mn has been provided in [27]. Before stating its definition, we
need to introduce some additional notation. For any given matrix A ∈ Sn , we denote'
Aij Aij > 0 and i -= j
N (A)ij :=and S(A) := A − N (A).
(6)
0otherwise
In [27], the authors defined the following set
Hn := {A ∈ Sn | S(A) ∈ Sn+ }.

(7)

Note that A = S(A) + N (A) ∈ Sn+ + Nn if A ∈ Hn . Also, for any A ∈ Nn , S(A) becomes a nonnegative
diagonal matrix and hence Nn ⊆ Hn . The detection whether A ∈ Hn is easy and can be done bychecking positivity of Aij (i -= j) and by a Cholesky factorization of S(A) (cf. Algorithm 4.2.4 in [17]).
Thus, by the inclusion relation (3), we see that the set Hn satisfies the desirable properties P1 and P2of Mn . However, S(A) is not necessarily positive semidefinite even if A ∈ Sn+ + Nn or A ∈ Sn+ . The
following theorem summarizes several properties of the set Hn .
Theorem 3.1 ([16] and Theorem 4.2 of [27]). Hn is a convex cone and Nn ⊆ Hn ⊆ Sn+ + Nn . If n ≥ 3,
these inclusions are strict and Sn+ -⊆ Hn . For n = 2, we have Hn = Sn+ ∪ Nn = Sn+ + Nn = COP n .
The construction of the set Hn is based on the idea of “nonnegativity-checking first and positive
semidefiniteness-checking second.” Now, we provide an alternative choice of Mn based on the ideaof “positive semidefiniteness-checking first and nonnegativity-checking second.”
For a given symmetric matrix A ∈ Sn , let P be an orthonormal matrix and Λ = Diag (λ1 , λ2 , . . . , λn ) bea diagonal matrix satisfying
A = P ΛP T .
(8)6

We are interested in decomposing A into a semidefinite matrix and a nonnegative matrix according to
the form A = P ΛP T . By introducing another diagonal matrix Ω = Diag (ω1 , ω2 , . . . , ωn ), consider thefollowing decomposition:
A = P (Λ − Ω)P T + P ΩP T
(9)
If Λ − Ω ∈ Nn , i.e., λi ≥ ωi (i = 1, 2, . . . , n) hold, then the matrix P (Λ − Ω)P T is positive semidefinite.
Thus, if we can find a suitable diagonal matrix Ω satisfying
λi ≥ ωi (i = 1, 2, . . . , n), [P ΩP T ]ij ≥ 0 (i, j = 1, 2, . . . , n, i ≤ j)(10)then (9) and (3) imply
A = P (Λ − Ω)P T + P ΩP T ∈ Sn+ + Nn ⊆ COP n .

(11)

We can determine whether such a matrix exists or not by solving the following linear optimization
problem with variables ωi (i = 1, 2, . . . , n) and α:(LP)P,Λ

Maximizesubject toαω i ≤ λi(n
[P ΩP T ]i,j = k=1 ωk pik pjk ≥ α(i = 1, 2, . . . , n)(i, j = 1, 2, . . . , n, i ≤ j)

(12)

(n
Note that (LP)P,Λ has the feasible solution at which ωi = λi (i = 1, 2, . . . , n) and α = minij k=1 λk pik pjk
and hence has an optimal solution with optimal value α∗ (P, Λ). If α∗ (P, Λ) ≥ 0 then there exists a matrixΩ for which the decomposition (10) holds. Based on these observations, we provide another alternate
Gns of Mn as follows:
Gns := {A ∈ Sn | α∗ (P, Λ) ≥ 0 for some orthonormal matrix P satisfying (8) }.

(13)

As stated above, if α∗ (P, Λ) ≥ 0 for a given decomposition A = P ΛP T then we can determine A ∈ Gns .
In this case, we just need to compute a matrix decomposition and to solve a linear optimization problem
with n + 1 variables and O(n2 ) constraints which implies that it is rather practical to use the set Gns as
an alternate of Mn Suppose that A ∈ Sn has n different eigenvalues. Then the possible orthonormal
matrices P = [p1 , p2 , · · · , pn ] are identifiable except for permutation and sign inversion of {p1 , p2 , · · · , pn }and by the representationn)A=
λi pi pTi
i=1

of (8), we see that the problem (LP)P,Λ is unique for any possible P . In this case, α∗ (P, Λ) < 0 witha specific P implies A -∈ Gns . However, otherwise (i.e., an eigenspace of A has at least dimension 2),α∗ (P, Λ) < 0 with a specific P does not necessarily guarantee that A -∈ Gns . So we cannot say that the
set Gns satisfies the desirable property P1 of Mn . However, as we see in Theorem 3.2 below, Gns may
satisfy the other desirable property P2.
%ns which are closely related to the set Gns and they might be
Let us introduce other new sets Gna and G
useful to clarify some theoretical properties or to improve our algorithm:
Gna
%ns
G

:=
:=

{A ∈ Sn | α∗ (P, Λ) ≥ 0 for any orthonormal matrix P satisfying (8) },{A ∈ Sn | α∗ (P, Λ) ≥ 0 for some arbitrary matrix P satisfying (8) }.

(14)
(15)

Note that if (10) holds for any arbitrary (not necessarily orthonormal) matrix P then (11) also holds,
which implies the following inclusions:%ns ⊆ Sn+ + Nn .
Gna ⊆ Gns ⊆ GMore precisely, the sets Gns , Gna and G*ns have the following properties.
7(16)

Theorem 3.2. The sets Gns , Gna and G*ns are cones and%s ⊆ S + + Nn ⊆ COP n
Sn+ ∪ Nn ⊆ Gna ⊆ Gns = com(Sn+ + Nn ) ⊆ Gnnwhere the set com(Sn+ + Nn ) is defined bycom(Sn+ + Nn ) := {S + N | S ∈ Sn+ , N ∈ Nn , S and N commute}.
For n = 2, we have%ns = Sn+ + Nn = COP n .
Sn+ ∪ Nn = Gna = Gns = com(Sn+ + Nn ) = G

Proof. We assume that A ∈ Sn is diagonalized as in (8) throughout the proof.

Suppose that the associated linear optimization problem (LP)P,Λ has an optimal solution (ω ∗ , α∗ ) :=(ω1∗ , . . . , ωn∗ , α∗ ). Then for any β ≥ 0, βA is diagonalized as in βA = P (βΛ)P T and (βω ∗ , βα∗ ) is an
optimal solution of the associated linear optimization problem (LP)P,βΛ . This implies that βA ∈ Gns%ns ) if A ∈ Gns (respectively A ∈ Gna , respectively A ∈ G%ns ) and
(respectively βA ∈ Gna , respectively βA ∈ G%ns are cones.
hence Gns , Gna and G
We have already seen that (16) holds. So it is sufficient to show that (i) Sn+ ∪ Nn ⊆ Gna and (ii)Gns = com(Sn+ + Nn ).

(i) Sn+ ∪ Nn ⊆ Gna : Let us show that Nn ⊆ Gna and Sn+ ⊆ Gna , respectively. Suppose that A ∈ Nn .
Then for all P the problem (LP)P,Λ has a feasible solution where (ω, α) = (λ1 , . . . , λn , 0) which implies
that A ∈ Gna . Suppose that A ∈ Sn+ , i.e., λi ≥ 0 (i = 1, 2, . . . , n). Then for all P the problem (LP)P,Λ
has a feasible solution where (ω, α) = (0, . . . , 0, 0) which implies that A ∈ Gna . Thus we have shownSn+ ∪ Nn ⊆ Gna .
(ii) Gns = com(Sn+ + Nn ): The inclusion Gns ⊆ com(Sn+ + Nn ) follows from the construction of the set Gns
as in (13) and (12). The converse inclusion Gns ⊇ com(Sn+ + Nn ) is also true since if A ∈ com(Sn+ + Nn )then there exist an orthonormal matrix P and diagonal matrices Θ = Diag (θ1 , θ2 , . . . , θn ) and Ω =Diag (ω1 , ω2 , . . . , ωn ) such thatA = P ΘP T + P ΩP T , P ΘP T ∈ Sn+ , P ΩP T ∈ Nn(see Theorem 1.3.12 of [19]) which implies that θi ≥ 0 (i = 1, 2, . . . , n) and that the problem (LP)P,Λ
with Λ = Θ + Ω has a nonnegative objective value at a solution (ω, α) where α = mini,j {[P ΩP T ]ij } ≥ 0.
The results for n = 2 follow from Theorem 3.1.
As we have seen in Theorem 3.1, Nn ⊆ Hn but Sn+ -⊆ Hn for n ≥ 3. Theorem 3.2 suggests that the set
Gns might be better than the set Hn in the sense of the desirable property (P2) of Mn . The following
examples show some contrasts between Hn , Gns and Gna .
Example 3.3. ConsiderThen, by the definition (6),

1
S(A) = A − N (A) =  1
1

1
1
2
A= 1
1 −1


1
−1  .
2 1
1
0
2 −1  −  1
−1
2
1

1
0
0

8

 
1
1
0
0 = 0
2
0
0 −1


0
−1  ∈ S3+
2

which implies that A ∈ H3 . Moreover,0
N (A)S(A) = S(A)N (A) =  1
1

1
0
0


1
0 
0

which implies that A = S(A) + N (A) ∈ com(S3+ + N3 ), and by Theorem 3.2, A ∈ G3s holds. Thus
H3 ∩ G3s -= ∅.
Example 3.4 (cf. Proof of Theorem 4.2 in [27]). Consider


1 −1
1
1 −1  .
A =  −1
1 −1
1Then A ∈ S3+ and by Theorem 3.2, we see that A ∈ G3s . However,

 
 
1 −1
1
0 0 1
1 −1
1 −1  −  0 0 0  =  −1
1
S(A) = A − N (A) =  −1
1 −1
1
1 0 0
0 −1
which implies that A -∈ H3 . Thus G3s \ H3 -= ∅.
Example 3.5. Consider

and let1A =  −1
1


−1 1
1 1 
1 1



1 −1 00
1 0  and N = A − S =  0
S =  −100 0
1Then S ∈ S3+ , N ∈ N3 and



0 0

SN = N S = 0 0
0 0

0
0
1


0
0 .
0

holds which implies that A ∈ com(S3+ + N3 ) ⊆ G3s . Moreover, if we set



√1
√1
√5
−1 0
3
14
42
 √1 − √3

1
√

−
0 2
P := 
, Λ :=
3
14
42 
1
2
4
√
√
√
0 0
− 3 − 1442

0−1  -∈ S3+
1


1
1 .
1


0
0 
2

then P and Λ satisfy (8) and the corresponding problem (LP)P,Λ is given as follows:
Maximizesubject toαω1 ≤ −1, ω2 ≤ 2, ω3 ≤ 2

1
1
1
− 13
3
3
14
1
1
1 
3


−3
ω1
+ ω2 − 14
3
3
1
1
1
−3 −3
− 17
3

3
− 14
9
14
3
7

− 17
3
7
2
7





 + ω3 

25
42
5
− 42
10
21

5
− 42
1
42
2
− 21

10
21
2
− 21
8
21



 ≥ αE.

By solving this problem, we know that α∗ (P, Λ) < 0. Thus the matrix A lies on G3s but not on G3a . Thus
G3s \ G3a -= ∅.
9In the next section, we will show numerical results of the following three algorithms:
Algorithm 1.1: The set Hn is used for Mn at Line 7 of Algorithm 1, i.e., the original algorithm
proposed in [27].
Algorithm 1.2: The set Gns is used for Mn at Line 7 of Algorithm 1.
Algorithm 2: An improved version of Algorithm 1.2. The set Gns is used for Mn at Line 7 of Algorithm%ns ⊆ COP n , some additional tests to remove simplices have
1. Moreover, based on the fact that G
been incorporated at Lines 7 and 10 of Algorithm 1.
Algorithm 2 An improved version of Algorithm 1.2
Input: A ∈ Sn , Mn ⊆ COP n
Output: “A is copositive” or “A is not copositive”
1: P ← {∆S };2: while P -= ∅ do3:
Choose ∆ ∈ P;
4:
if v T Av < 0 for some v ∈ V ({∆}): then5:
return “A is not copositive”;
6:
end if%ns then7:
if V∆T AV∆ ∈ G
8:
P ← P \ {∆};
9:
else
10:
if V∆T AV∆ ∈ Gns then11:
P ← P \ {∆};12:else
* ← {∆1 , ∆2 };13:
partition ∆ into ∆ = ∆1 ∪ ∆2 and set ∆14:
for p = 1, 2 do
%ns then15:
if V∆Tp AV∆p ∈ G
p
*
*
16:
∆ ← ∆ \ {∆ };17:
end if18:
end for*19:
P ← P \ {∆} ∪ ∆;20:
end if21:
end if
22: end while
23: return “A is copositive”;
The details of the added steps in Algorithm 2 are as follows. Suppose that we have a diagonalization of
the form (8) in advance.
At Line 7, we will solve an additional LP but need not diagonalize V∆T AV∆ . Let P and Λ be matrices
satisfying (8). Then the matrix V∆T P gives the diagonalization of V∆T AV∆ , i.e.,V∆T AV∆ = V∆T (P ΛP T )V∆T = (V∆T P )Λ(V∆T P )T%s by solving (LP) T
while V∆T P is not necessarily orthonormal. Thus we can test V∆T AV∆ ∈ GnV P,Λ .
∆

%ns does not detected at Line 7, we will check whether V T AV∆ ∈ Gns at Line 10. Similarly
If V∆T AV∆ ∈ G∆
to Algorithm 1.2 (where the set Gns is used for Mn at Line 7 of Algorithm 1), we will diagonalize V∆T AV∆
as V∆T AV∆ = P ΛP T with an orthonomal matrix P and a diagonal matrix Λ, and solve (LP)P,Λ .
10

At Line 15, we neither need diagonalize V∆Tp AV∆p nor to solve any further LPs. Let ω ∗ ∈ Rn be an
optimal solution of (LP)V T P,Λ obtained at Line 7 and let Ω∗ := Diag (ω ∗ ). Then the feasibility of ω ∗∆implies the positive semidefiniteness of the matrix V∆Tp P (Λ − Ω∗ )P T V∆p . Thus, if V∆Tp P Ω∗ P T V∆p ∈ Nn
then we see thatV∆Tp AV∆p = V∆Tp P (Λ − Ω∗ )P T V∆p + V∆Tp P Ω∗ P T V∆p ∈ Sn+ + Nn%ns .
and that V∆Tp AV∆p ∈ G

4Numerical results

We implemented Algorithms 1.1, 1.2 and 2 in MATLAB to compare the performance of those algorithms.
As the test-instances, we used the following matrix
Bγ := γ(E − AG ) − E(17)where γ ≥ 1, E ∈ Sn is the matrix whose elements are all one and the matrix AG ∈ Sn is the adjacency
matrix of a given undirected graph G with n nodes. The matrix Bγ comes from the maximum clique
problem. The maximum clique problem is to find a clique (complete subgraph) of maximum cardinality
in G. It has been shown in [10] that the maximum cardinality, the so-called clique number ω(G), isequal to the optimal value of
ω(G) = min{γ ∈ N | Bγ ∈ COP n }.
Thus, the clique number can be found by checking the copositivity of Bγ at most for γ = n, n − 1, . . . , 1.
Note that in [27], the authors showed that Bγ -∈ COP n if γ < ω(G), Bγ ∈ COP n \int (COP n ) if γ = ω(G)and Bγ ∈ int (COP n ) if γ > ω(G) (see Proposition 3.2 of [27]). The results and Corollary 2.4 imply
that the algorithms may fail to terminate if γ = ω(G). In [27], the authors have provided a modified
copositive program to avoid such difficulties as the following theorem.
Theorem 4.1 (Theorem 3.3 in [27]). The clique number ω(G) of a given graph G can be obtained fromthe following modified copositive program:
ω(G) = min{γ ∈ N | Bγ + ρE ∈ COP n }if 0 ≤ ρ < 1/ω(G). Moreover, Bγ + ρE is strictly copositive for any γ ≥ ω(G) and ρ > 0.
An aim of the implementation is to explore the difference of behaviors between the choices Mn = Hn
and Mn = Gns rather than to compute the clique number efficiently. So we conducted our experiment
to examine Bγ for various values of γ at intervals of 0.1 around the value ω(G) (Tables 1 and 2 on page17). We also solved the modified problem in Theorem 4.1 to confirm the efficacy of the modification(Table 5 on page 19).
Figure 1 on page 17 shows our instances for G that have been used in [27]. We know the clique numbersof G8 and G12 are ω(G8 ) = 3 and ω(G12 ) = 4, respectively.
For a given A ∈ Sn , we used the MATLAB command “[P, Λ] = eig(A)” to obtain the diagonalized
form (8). As already mentioned above, α∗ (P, Λ) < 0 with a specific P does not necessarily guarantee%ns ). Thus, it not strictly accurate to say that we have used Gns (G%ns ) for Mn andthat A -∈ Gns (A -∈ G
the algorithms may miss some removable ∆’s. Note that this may have some effect on speed but not
11

on termination of the algorithm since the termination is guaranteed by the subdivision rule satisfying
δ(P) → 0 where δ(P) is defined by (5).
The performance of algorithms is also influenced by strategies to refine the simplex ∆ used at Line 10of Algorithm 1 or Line 13 of Algorithm 2. We employed the most classical longest-edge bisection rule
as a common strategy among the experiments, i.e., we choose the longest edge of a given simplex ∆and bisect the edge with the fixed bisection ratio 1 : 1. It is well known that this type of longest-edge
bisection rule generates a sequence such that δ(P) → 0 ([18], see also [12]).
There have been different strategies on refinement of simplices [7, 8, 27]. Among others, an efficient
strategy for the set Mn = Hn has been provided in [27]. We will discuss our improvement on refinement
for the set Mn = Gns in the next section.
We tested our implementation on a 3.07GHz Core i7 machine with 12 GB of RAM. Tables 1 and 2
represent the numerical results for the graphs G8 and G12 , respectively. In both tables, the symbol “−”means that the algorithm did not terminate within 6 hours. These results may come from the fact thatfor each graph G, the matrix Bγ lies on the boundary of the copositive cone COP n when γ = ω(G)(ω(G8 ) = 3 and ω(G12 ) = 4).
We observe similar trends in Tables 1 and 2, and explain the implications of our results using Table 2
on page 18 for the larger graph G12 . The results in Table 2 imply that- at any γ ≥ 5.2, Algorithm 1.2 terminates in one iteration and its execution time is faster than theone of Algorithm 1.1. The reason may be that, as shown in Theorem 3.2, the set Gns is a relatively
large subset of COP n and useful to check copositivity of the matrix A if A is (strictly) copositive.
- at any γ ∈ {5.1, 5.0, . . . , 4.7}, Algorithm 1.1 terminates within 4.2 hours while Algorithm 1.2 does
not terminate within 6 hours. The results imply that the execution time of Algorithm 1.2 is slower
than the one of Algorithm 1.1 since Algorithm 1.2 requires additional computation for solving an
n + 1 variables and O(n2 ) constraints linear optimization problem at each iteration. The execution
time of the improved version, Algorithm 2, is substantially better than the one of Algorithm 1.2,
but is not better than one of Algorithm 1.1.
- at each γ < 4, Algorithms 1.1, 1.2 and 2 have no significant differences in terms of the number of
iterations since both of algorithms work to find a v ∈ V ({∆}) such that v T (γ(E − AG ) − E)v < 0
while its computational time depends on our choice of simplex refinement strategy but not on the
choice of Mn .
In view of the above observations, we conclude that Algorithm 1.2 with the choice Mn = Gns might be
a promising way to check copositivity of a given matrix A when A is strictly copositive. In addition,
the improved technique used in Algorithm 2 has a pronounced effect on the number of iterations andhence the execution time of Algorithm 1.2. This can be seen in Tables 3 and 4 where the columns “Line8,” “Line 11” and “Line 19” show the number of simplices removed at each line of Algorithm 2. The
elimination of these simplices contributes to improve Algorithm 1.2.
Next we implemented the three algorithms to solve the modified programs provided in Theorem 4.1. The
obtained results are shown in Table 5 on page 19. Note that Theorem 4.1 involves only integer valuesof γ. So we tested the matrices B4 + 0.199E (0.199 < 15 ) and B3 + 0.249E (0.249 < 14 ) of the graph G8with ω(G8 ) = 3, and the matrices B5 + 0.166E (0.166 < 16 ) and B4 + 0.199E (0.199 < 15 ) of the graph
G12 with ω(G12 ) = 4. These results may be compared to the results for γ = 3.9 in Table 1 for the graphG8 and the ones for γ = 4.9 in Table 2 for the graph G12 since the copositivity of B3.9 implies ω(G8 ) ≤ 3
for G8 and the copositivity of B4.9 implies ω(G12 ) ≤ 4 for G12 , respectively. These comparisons suggest
12that for all algorithms, Algorithms 1.1, 1.2 and 2, the modified programs have positive effects especiallyfor detecting ω(G12 ).

5

Improved strategies for refinement of simplices

In this section, we discuss our strategies for refinement of simplices for the improved version, Algorithm
2.
We first introduce the strategy for Algorithm 1.1 provided in [27]. Suppose that we choose Mn = Nnat Line 10 of Algorithm 1. If V∆T AV∆ -∈ Mn then there exist i and j such that viT Avj < 0. If i = j thensince vi ∈ ∆, we find that A -∈ COP n and the algorithm terminates and otherwise it would be natural
to partition an edge {vi , vj } which attains the optimal value of mini,j∈{1,2,...,n},i%=j viT Avj < 0.
Adopting the idea for the case Mn = Hn defined in (7), Sponsel, Bundfuss and Dür [27] suggested astrategy to partition an edge {vi , vj } which gives the optimal value ofmini,j∈{1,2,...,n},i%=j

S(V∆T AV∆ )ij xi xj

and have shown the numerical results using the strategy.
%ns where the definitions are given in (13)
We adopt the same idea for the cases Mn = Gns and Mn = G
T
s
T
s
%
and (15), respectively. If V∆ AV∆ -∈ Gn ( V∆ AV∆ -∈ Gn ) then for any orthonormal matrix (for any
arbitrary) P which gives a diagonalization V∆T AV∆ = P ΛP T , the feasible linear optimization (LP)P,Λ
has the negative optimal value α∗ (P, Λ) < 0.
Our strategy, we call it the “negative-edge bisection rule,” is to partition an edge {vi , vj } which attains
the optimal valuen)ωk∗ pik pjk = α∗ (P, Λ) < 0
k=1

∗

∗

at the optimal solution (ω , α ). Unfortunately, our negative-edge bisection rule does not guarantee that
δ(P) → 0, and we found that Algorithm 2 with the rule fails to terminate for some instances. To improve
the termination behavior, we insert the longest-edge bisection steps periodically during performing thenegative-edge bisection refinement. Note that recently Dickinson [12] showed that this strategy is not
sufficient to ensure δ(P) → 0.
We implemented Algorithm 2 with the negative-edge bisection steps and tested it for checking copositivityof the matrix B3.4 of the graph G8 and the matrix B5.1 of the graph G12 for each of which the numberof iteration of Algorithm 2 has a pronounced jump (see Tables 1 and 2).
Table 6 shows the performance of adding the negative-edge bisection steps for these two instances. The
first row of the table represents the number of inserted longest-edge bisection steps (LEB steps) per
twelve negative-edge bisection steps (NEB steps) where ∞ means that only longest-edge bisection steps
and no negative-edge bisection step have been taken.
Note that at each iteration of Algorithm 2, we may have two optimal solutions of (LP)P,Λ , i.e., the one%ns and the one obtained at Line 10 by detecting V T AV∆ ∈ Gns .
obtained at Line 7 by detecting V∆T AV∆ ∈ G∆
The third column of Table 6 shows which optimal solution is used for the negative-edge bisection steps.

For each case, the number of iterations shows the average number of iterations required after three-times
13execution since we have randomly chosen the negative-edge {vi , vj } to be partitioned if there are multiple
candidates.
We observe from Table 6 that the negative-edge bisection strategy using the solution obtained at Line
10 (by detecting V∆T AV∆ ∈ Gns ) has a positive effect to reduce the number of iterations for checking B5.1
of G12 and that the effect is monotonically increasing with NEB frequency.

6

Concluding remarks

In this paper, we proposed a new branch and bound type algorithm for testing copositivity of a given
symmetric matrix based on the algorithm proposed in [27]. Two features of our algorithm are%s which are relatively large subsets of COP n
1. we have introduced new classes of matrices Gns and Gn
and work well to check copositivity of a given matrix A ∈ Sn (see Theorem 3.2) , and%ns in checking copositivity, we only have to solve a linear optimiza2. for incorporating the sets Gns or Gtion problem with n + 1 variables and O(n2 ) constraints after computing a singular value matrix
decomposition, which implies that our algorithm is not so time-consuming.
Our algorithm determined the copositivity of some instances within a small number of iterations if they
are strictly copositive. We also provided the negative-edge bisection strategy which aims to improverefinement of simplices for checking copositivity of the matrix.
Further research will include- further improvement in checking A ∈ Gns : We only solve the problem (LP)P,Λ for a specific P
to check A ∈ Gns . This is sufficient if A has n different eigenvalues. However, otherwise (i.e., aneigenspace of A has at least dimension 2), we may miss the fact A ∈ Gns . Solving the problem(LP)P,Λ with other possible P s might be effective for further improvement in checking A ∈ Gns , ormore specifically the results at around γ = 5 in Table 2.
%ns . Theorem 3.2 and Examples 3.3, 3.4 and 3.5 show the- more observations on the sets Gns , Gna and G+
relationships among these sets and Sn ∪ Nn , com(Sn+ + Nn ), Sn+ + Nn and COP n . To observe how
those sets are different and what properties they have will be of research interest in the future.
Note that Table 1 shows an interesting result concerning the second point. Let us see the result at
γ = 4.0 of Algorithm 1.2. The multiple number of iterations at γ = 4.0 implies that we could not findB4.0 ∈ Gns at the first iteration for a certain orthonormal matrix P satisfying (8). Recall that the matrix
Bγ is given by (17). It follows from the fact E − AG ∈ Nn ⊆ Gns and from the result at γ = 3.5 in Table
1 that0.5(E − AG ) ∈ Gns and B3,5 = 3.5(E − AG ) − E ∈ Gns .
Thus the fact that we could not find whether the matrix
B4.0 = 4.0(E − AG ) − E = 0.5(E − AG ) + B3.5
lies on the set Gns might suggest that the set Gns = com(Sn + Nn ) is not convex.

14

Acknowledgment
The authors would like to sincerely thank the anonymous reviewers for their thoughtful and valuablecomments which have significantly improved the paper. Among others, one of the reviewers brought us%ns . The authors greatly appreciate his/her constructive
the idea of improving Algorithm 1.2 using the set G
criticism and suggestions on improving this work. They are also grateful to Yasuaki Matsukawa for his
careful reading of an earlier version of the manuscript.

References[1] A. Berman, Cones, Matrices and Mathematical Programming, Lecture Notes in Economics andMathematical Systems 79, Springer Verlag 1973.
[2] A. Berman and N. S. Monderer, Completely Positive Matrices, World Scientific Publishing, 2003.
[3] I. M. Bomze, M. Dür, E. De Klerk, C. Roos, A. J. Quist and T. Terlaky, On copositive programmingand standard quadratic optimization problems, Journal of Global Optimization 18(2000) 301-320.
[4] I. M. Bomze and E. De Klerk, Solving standard quadratic optimization problems via linear, semidefinite and copositive programming, Journal of Global Optimization 24(2002) 163-185.
[5] I. M. Bomze, Copositive optimization - recent developments and applications, European Journal ofOperational Research 216(2012) 509-520.
[6] I. M. Bomze and G. Eichfelder Copositivity detection by difference-of-convex decomposition andω-subdivision, Mathematical Programming, Series A 138(2013) 365-400.
[7] S. Bundfuss and M. Dür, Algorithmic copositivity detection by simplicial partition, Linear Algebra
and its Applications 428(2008) 1511-1523.
[8] S. Bundfuss, Copositive matrices, copositive programming, and applications, Ph.D. Dissertation,
TU Darmstadt, 2009.
Online at http://www3.mathematik.tu-darmstadt.de/index.php?id=483[9] S. Burer, On the copositive representation of binary and continuous nonconvex quadratic programs,
Mathematical Programming 120(2009) 479-495.
[10] E. de Klerk, D. V. Pasechnik, Approximation of the stability number of a graph via copositive
programming, SIAM Journal on Optimization 12(2002) 875-892.
[11] P. H. Diananda, On non-negative forms in real variables some or all of which are non-negative,
Mathematical Proceedings of the Cambridge Philosophical Society 58(1962) 17-25.
[12] P. J. C. Dickinson, On the exhaustivity of simplicial partitioning, Journal of Global Optimization
58(2014) 189-203.
[13] P. J. C. Dickinson and L. Gijben, On the computational complexity of membership problems forthe completely positive cone and its dual, Computational Optimization and Applications 57(2014)403-415.
[14] M. Dür, Copositive programming - a survey, Recent Advances in Optimization and Its Applications
in Engineering Springer-Verlag, pp.3-20, 2010.

15

[15] W. Fenchel, Convex cones, sets and functions, mimeographed notes by D. W. Blackett, Princeton
Univ. Press, Princeton, N. J., 1953.
[16] M. Fiedler and V. Pták, On matrices with non-positive off-diagonal elements and positive principal
minors, Czechoslovak Mathematical Journal 12(1962) 382-400.
[17] G. H. Golub and C. F. Van Loan. Matrix Computations, Johns Hopkins University Press; Third
edition, 1996.
[18] R. Horst, On generalized bisection of n-simplices, Mathematics of Computation 66(1997) 691-698.
[19] R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press, 1985.
[20] F. Jarre and K. Schmallowsky, On the computation of C ∗ certificates, Journal of Global Optimization
45(2009) 281-296.
[21] Y. Matsukawa and A. Yoshise, A primal barrier function Phase I algorithm for nonsymmetric conic
optimization problems, Japan Journal of Industrial and Applied Mathematics 29(2012) 499-517.
[22] T. S. Motzkin, E. G. Straus, Maxima for graphs and a new proof of a theorem of Turan, Canadian
Journal of Mathematics 17(1965) 533-540.
[23] K. G. Murty, S. N. Kabadi, Some NP-complete problems in quadratic and nonlinear programming,
Mathematical Programming 39(1987) 117-129.
[24] J. Povh and F. Rendl, A copositive programming approach to graph partitioning, SIAM Journal
on Optimization 18(2007) 223-241.
[25] J. Povh and F. Rendl, Copositive and semidefinite relaxations of the quadratic assignment problem,
Discrete Optimization 16(2009) 231-241. .
[26] R.T. Rockafellar and R.J-B. Wets. Variational Analysis. Springer, 1998.
[27] J. Sponsel, S. Bundfuss and M. Dür, An improved algorithm to test copositivity, Journal of Global
Optimization 52(2012) 537-551.
[28] J. F. Sturm and S. Zhang. On cones of nonnegative quadratic functions. Mathematics of Operations
Research, 28(2003) 246-267.
[29] Y. Yoshise, Y. Matsukawa, On optimization over the doubly nonnegative cone, Proceedings of 2010
IEEE Multi-conference on Systems and Control (2010) 13-19.

16

Figure 1: The graphs G8 with ω(G8 ) = 3 (left) and G12 with ω(G12 ) = 4 (right).

Table 1: Results for the graph G8

Alg. 1.1. (Hn )
γ
2.8
2.9
3.0
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
4.0
4.1
4.2
4.3
4.4
4.5

# of iterations
2247
1609
3003
1509
469
395
369
209
115
79
63
39
23
17
17
7
7

CPU time (s)
0.368
0.243
0.324
0.152
0.047
0.039
0.037
0.021
0.012
0.009
0.007
0.021
0.002
0.002
0.002
0.001
0.001

Alg. 1.2. (Gns )# of iterations
2599
2273
11227
6243
4365
3295
1
1
1
1
1
385
1
1
1
1
1

17

CPU time (s)35.74131.714137.45378.62755.30842.249
0.012
0.011
0.012
0.011
0.011
4.769
0.011
0.011
0.011
0.011
0.012

%ns )
Alg. 2. (Gns and G# of iterations
2478
2190
6326
3273
1585
1491
1
1
1
1
1
339
1
1
1
1
1

CPU time (s)65.16055.401
129.108
69.27933.20433.043
0.013
0.012
0.012
0.012
0.011
7.366
0.014
0.011
0.011
0.011
0.011

Table 2: Results for G12

Alg. 1.2. (Gns )Alg. 1.1. (Hn )
γ
3.8
3.9
4.0
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
5.0
5.1
5.2
5.3
5.4
5.5
5.6
5.7
5.8
5.9
6

# of iterations
4084
4080
1125035
762931610071
569661407201305521206949141383110641
90877
44731
26171
15045
10239
6977
4717

CPU time (s)
1.753
1.755
14789.042
8166.234
6121.059
5375.867
3225.592
1693.459
611.519262.289154.699
102.006
22.361
8.353
3.593
2.167
1.325
0.839

# of iterations
4088
4088
1
1
1
1
1
1
1
1
1

CPU time (s)134.426
136.621
0.052
0.043
0.056
0.043
0.052
0.055
0.058
0.063
0.053

%s )
Alg. 2. (Gns and Gn# of iterations
4088
4088
219156141907
59063
1
1
1
1
1
1
1
1
1

Table 3: The number of simplices removed at each line of Algorithm 2 for G8
γ
2.8
2.9
3.0
3.1
3.2
3.3
3.4

Line 8
169
179
2150
1003
515
359

Line 11
175
220
1000
586
235
325

18

Line 19
5
3
27
96
86
124

CPU time (s)253.307269.66811135.302
6882.747
2848.623
0.053
0.061
0.052
0.059
0.056
0.058
0.061
0.057
0.064

Table 4: The number of simplices removed at each line of Algorithm 2 for G12
γ
3.8
3.9
4.0
4.8
4.9
5.0
5.1

Line 8
0
0
89561
56157
24983

Line 11
19
21
18546
12744
4549

Line 19
0
0
2943
4106
0

Table 5: Results for modified programs provided in Theorem 4.1

Graph

γ

ρ

G8
G8
G12
G12

4
3
5
4

0.199 < 1/5
0.249 < 1/4
0.166 < 1/6
0.199 < 1/5

Alg. 1.1. (Hn )# ofCPU
iterations time (s)
3
0.015
47
0.021
4771
0.166
337997 1829.384

Alg. 1.2. (Gns )# ofCPU
iterations
time (s)
1
0.055
1
0.056
1
0.068
997653 41301.295%s )Alg. 2. (Gns + G
n
# of
CPU
iterations time (s)
1
0.084
1
0.087
1
0.109
168080 8580.701

Table 6: Results of Algorithm 2 with improved strategies for refinement of simplices

Graph

γ

(LP)P,Λ

G8
G8
G12
G12

3.4
3.4
5.1
5.1

Line 7
Line 10
Line 7
Line 10

∞/12
# of
CPU
iterations time (s)
1491
33
1491
33
59063
2849
59063
2849

6/12
# of
CPU
iterations time (s)
3185
66
1545
33
44583
2111
43499
2063

19

4/12
# of
CPU
iterations time (s)
6850
137
1284
27
125443
6035
36973
1755

3/12
# of
CPU
iterations time (s)
7367
158
1595
34
1304566403
33945
1612