Department of Social Systems and Management
Discussion Paper Series

No.1284

A Primal Barrier Function Phase I Algorithm for
Nonsymmetric Conic Optimization ProblemsbyYasuaki MATSUKAWA and Akiko YOSHISE

November 2011

UNIVERSITY OF TSUKUBA
Tsukuba, Ibaraki 305-8573
JAPANA Primal Barrier Function Phase I Algorithm for Nonsymmetric
Conic Optimization Problems
Yasuaki Matsukawa∗ and Akiko Yoshise†
November 2011
Revised March 2012

Abstract
We call a positive semidefinite matrix whose elements are nonnegative a doubly nonnegative
matrix, and the set of those matrices the doubly nonnegative cone (DNN cone). The DNN cone
is not symmetric but can be represented as the projection of a symmetric cone embedded in ahigher dimension. In [16], the authors demonstrated the eﬃciency of the DNN relaxation using
the symmetric cone representation of the DNN cone. They showed that the DNN relaxation gives
significantly tight bounds for a class of quadratic assignment problems, but the computational time
is not aﬀordable as long as we employ the symmetric cone representation. They then suggested a
primal barrier function approach for solving the DNN optimization problem directly, instead of using
the symmetric cone representation. However, most of existing studies on the primal barrier function
approach have assumed the availability of a feasible interior point. This fact means that thosestudies are not inextricably tied to the practical usage. Motivated by these observations, we propose
a primal barrier function Phase I algorithm for solving conic optimization problems over the closed
convex cone K having the following properties: (a) its interior int K is not necessarily symmetric,(b) a self-concordant function f is defined over int K, and (c) its dual cone K ∗ is not explicit oris intractable, all of which are observed when K is the DNN cone. We analyze the algorithm andprovide a suﬃcient condition for finite termination.

1

Introduction

This paper deals with the conic optimization problem given by
Minimizesubject to�c, x�
Ax = b,x∈K(1)where K ⊆ Rp is a closed convex cone, A : Rp → Rq is a linear operator, b ∈ Rq , c ∈ Rp and �·, ·� isan inner product on Rp . In the last two decades, many studies have been done on the symmetric conic∗ Graduate School of Systems and Information Engineering, University of Tsukuba, Tsukuba, Ibaraki 305-8573, Japan.
e-mail: matsuk50@sk.tsukuba.ac.jp† Corresponding author. Graduate School of Systems and Information Engineering, 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optimization where the closed convex cone K is self-dual, i.e., K satisfies
K = K ∗ := {y | �x, y� ≥ 0 for all x ∈ K}.
A typical example of such problems is the semidefinite optimization problem with the semidefinite cone
(SDP cone) given by��n
K = S+:= X ∈ Rn×n | X = X T , ∀d ∈ Rn , dT Xd ≥ 0 .

The eﬀectiveness of semidefinite relaxation (SDP relaxation) approach for combinatorial optimization
problems has been supported in many papers (see, e.g., [2, 4, 15, 5, 3, 10, 7]).
n
In the SDP relaxation, the solution matrix X is restricted to X ∈ S+
and often meant to be nonnegative,
i.e.,
�
�
X ∈ Rn×n:= X ∈ Rn×n | xij ≥ 0 (i, j = 1, 2, . . . , n) .
+We call a positive semidefinite matrix whose elements are nonnegative a doubly nonnegative matrix, andthe setn
Dn := S+
∩ Rn×n+of those matrices the doubly nonnegative cone (DNN cone). The DNN cone can be represented as the
projection of the direct sum of the SDP cone and the nonnegative orthant as follows:n
Dn = {X | (X, Y ) ∈ S+
× Rn×n, X = Y }.
+

(2)

n
Note that the set S+
×Rn×n
is a symmetric cone, and we can adopt existing symmetric conic optimization+
solvers to solve the DNN relaxation. Matsukawa and Yoshise [16] demonstrated the eﬃciency of thisapproach and reported that- the DNN relaxation gives significantly tighter bounds than the SDP relaxation for a class of quadratic
assignment problems: more than 95% accuracy can be achieved for many instances from the
QAPLIB library [1], but- the size of the problem grows too much large due to the symmetric cone representation (2): even a
mid-size problem with n = 15 requires more than 7 hours to solve the relaxation problem with
SDPA Online Solver [13] in 2010.
Motivated by the above observations, in [16], the authors provided basic properties of the DNN cone
aiming to develop another approach, and showed that(i) the dual cone (Dn )∗ of the DNN cone isnn(Dn )∗ = S+
+ Rn×n:= {x + y | x ∈ S+
, y ∈ Rn×n
}
+
+

(3)

(see Proposition 4.1 of [16]) and(ii) the interior of the DNN cone int Dn is a hyperbolic cone characterized by the self-concordant function
f : int Dn → Rn �n�f (X) := − log det (X) −log xij(4)
i=1 i=1(see Proposition 4.2 of [16] and cf. also [6, 12]).

2

The property (ii) above ensures that we can apply the primal barrier function approach provided in,e.g., [8, 11, 9] to solve the DNN optimization problem. In particular, the paper [9] dealt with thenonsymmetric conic optimization where only the primal cone K has a self-concordant function on its
interior int K. However, all of these studies assume that a feasible interior point x0 satisfying
Ax0 = b, x0 ∈ int K(5)

is obtained a priori, which means that those studies are not inextricably tied to the practical usage. In
fact, with the appearance of the primal-dual interior point algorithm in 1989 for linear optimization, thenumber of studies on the primal barrier function algorithm diminished and there are few papers on the
primal infeasible interior point algorithm which may start from an infeasible interior point.
Recently, Skajaa, Jørgensen and Hansen [14] proposed a homogeneous interior-point algorithm for nonsymmetric conic optimization which does not need to assume the existence of a feasible interior point.
Note that, in either [14] or [9], it has been assumed that we can easily check whether an element s lies inthe dual cone K ∗ of K. However, as described in (ii) above, the dual cone (Dn )∗ of the DNN cone Dn
is given by (3) and it is not necessarily easy to find whether s ∈ (Dn )∗ or not. A new and purely primalalgorithm would be expected to find a feasible interior point satisfying (5) for the DNN optimization
problem.
In this paper, we propose an algorithm to find a feasible interior point of the conic optimization (1)where the closed convex cone has the following properties:(a) int K is not necessarily symmetric, i.e., it can be not self-dual or not homogeneous,(b) a self-concordant function f is defined over int K, and(c) its dual cone K ∗ is not explicit or is intractable.
Our algorithm is a purely primal barrier function algorithm based on the self-concordant function f ,and does not require any knowledge of the dual cone K ∗ .
The paper is organized as follows. In Section 2, we summarize basic properties of the self-concordant
function according to [11]. The primal barrier function Phase I algorithm is provided in Section 3.
The algorithm consists of the feasibility step and the well-known centering steps. After analyzing thefeasibility step in Section 4, we provide a suﬃcient condition for finite termination of the algorithm andestimate the required number of iterations in Section 5. Concluding remarks are given in Section 6.

2

Basic properties of the self-concordant function

Our analysis crucially relies on Section 2 of [11] where many properties of the self-concordant function
have been investigated. In this section, we provide some basic properties of the self-concordant functionused in the paper.
Let Df be an open convex set and f : Df → R be a functional f ∈ C 2 . We denote by g(x) and H(x)the gradient and the Hessian of f (x) at x, respectively. We obtain the following proposition.
Proposition 2.1 (Proposition 1.5.7 of [11]). If x, y ∈ Df then� 1�
�g(y) = g(x) + H(x)(y − x) +H(x + t(y − x)) − H(x) (y − x)dt.
03

In what follows, we suppose that H(x) is positive definite for all x ∈ Df which implies that f is strictly
convex. The local inner (or the intrinsic) product and the induced norm are defined as follows:��u, v�x := �u, H(x)v�, �u�x := �u, u�x .
(6)
Let Bx (x, λ) be the open ball given by

Bx (x, λ) := {y | �y − x�x < λ}.

In [11], the self-concordant function is defined as follows.

Definition 2.2 (Section 2.2.1 of [11]). A functional f : Df → R is said to be (strongly nondegenerate)
self-concordant if for all x ∈ Df we haveBx (x, 1) ⊆ Df

and if whenever y ∈ Bx (x, 1) we have
1 − �y − x�x ≤

(7)

�v�y
1
≤
for all v �= 0.
�v�x
1 − �y − x�x

Since H(x) is positive definite, H(x)−1 exists for any x ∈ Df . By introducing the notation
Hx (y) := H(x)−1 H(y),

(8)

(9)

we see the following properties.
Theorem 2.3 (Theorem 2.2.1 of [11]). Assume that the functional f has the property that Bx (x, 1) ⊆ Df
for all x ∈ Df .
(i) f is self-concordant iﬀ for all x ∈ Df and y ∈ Bx (x, 1)
�Hx (y)�x , �Hx (y)−1 �x ≤1.
(1 − �y − x�x )2(i) f is self-concordant iﬀ for all x ∈ Df and y ∈ Bx (x, 1)
�I − Hx (y)�x , �I − Hx (y)−1 �x ≤1
− 1.
(1 − �y − x�x )2

We denote by n(x) the Newton direction for f at x which is given byn(x) = −H(x)−1 g(x)(10)

H(x)n̄(x) + g(x) = A∗ y,An̄(x) = 0,

(11)
(12)

and by n̄(x) the projected Newton direction which is the unique solution of the system

where A∗ denotes the adjoint of A satisfying �x, A∗ y� = �Ax, y� for any x ∈ Rp and y ∈ Rq . The
projected Newton direction n̄(x) and the Newton direction n(x) have the following relation:�n̄(x), n(x)�x

=
=
=
=

�n̄(x), n̄(x) − H(x)−1 A∗ y�x(by (10) and (11))

�n̄(x), n̄(x)�x − �n̄(x), A y� (by (6) )�n̄(x)�2x − �An̄(x), y��n̄(x)�2x . (by (12))∗(13)

The theorem below gives an upper bound of the diﬀerence between the function value f (y) and its
second-order approximation at x of the self-concordant function f .
4

Theorem 2.4 (Theorem 2.2.2 of [11]). If f is self-concordant, x ∈ Df and y ∈ Bx (x, 1), then|f (y) − qx (y)| ≤�y − x�3x3(1 − �y − x�x )where qx (y) is the second-order approximation of f (y) at x given by1
qx (y) := f (x) + �g(x), (y − x)� + �(y − x), H(x)(y − x)�
2
1= f (x) − �n(x), (y − x)�x + �y − x�2x2Here the second equality above follows from the definitions of (10) and (6). At the last of this section,
we give two propositions on the complexity value of f .
Proposition 2.5 (Section 2.3.1 of [11]). Let f : Df → R be a strongly nondegenerate self-concordantbarrier. Let ϑf be the complexity value of f given byϑf := sup �n(x)�2x .
x∈Df

Then the following inequalities hold.
�n̄(x)�x ≤ �n(x)�x ≤

�

ϑf .

For the self-concordant function f over the interior int Dn of the DNN cone Dn in (2), we obtain thefollowing lemma by simple calculations.
Proposition 2.6. Let f be a self-concordant barrier on int Dn given by (4). The complexity value
ϑf :=

sup
x∈int D n

�n(x)�2xof the function f is n2 + n.

3

A primal barrier function Phase I algorithm for the nonsymmetric conic optimization problem

Let K be a closed convex cone which is not necessarily symmetric, and f be a self-concordant functional
defined over Df := int K. In this section, we propose an algorithm for finding a feasible interior pointof the primal (nonsymmetric) optimization problem of (1).
Recall that the Newton direction n(x) and the projected Newton direction n̄(x) for f at x are given by(10) and (11) – (12), respectively. Here we outline the algorithm.
Let τ ∈ (0, 1) and choose an initial point x0 ∈ int K. Define
r0 := Ax0 − b, Sγ := {x ∈ int K | Ax − b = γr0 }.

(14)

Note that r0 is not necessarily zero and x0 may be an infeasible interior point. We impose the following
assumption which is often satisfied by the DNN relaxation problems of combinatorial optimization
problems (cf., e.g., [2, 4, 3, 10, 7]).
5

Assumption 3.1. The set
is bounded.

S1 = {x ∈ int K | Ax − b = r0 }If the original problem has a feasible interior point then S0 �= ∅ and otherwise S0 = ∅. In any cases, wehave the following lemma.
Lemma 3.2. Suppose that Assumption 3.1 holds.
(i) The set∪γ∈[0,1] Sγ

is a bounded convex set.

(ii) If Sγ �= ∅ then the analytical center z(γ) for f in Sγ , which is the minimizer of f over Sγ , uniquely
exists and satisfiesn̄(z(γ)) = 0
for any γ ∈ [0, 1].
(iii)

inf{f (x) | x ∈ ∪γ∈[0,1] Sγ } > −∞Proof. (i): The convexity of the set ∪γ∈[0,1] Sγ directly follows from the definition and the convexity of
K. Suppose that ∪γ∈[0,1] Sγ is unbounded. Then there exists a sequence {(γ k , xk )} satisfying
γ k ∈ [0, 1], Axk − b = γ k r0 , xk ∈ int K, �xk � → +∞.
Let us consider the sequencey k :=

xk∈ int K.
�xk �

Since {y k } and {γ k } are bounded and �y k � = 1 for any k, by taking a subsequence if necessary, we findȳ �= 0 satisfying
Aȳ = 0, ȳ ∈ K.
Thus, for any α > 0, x0 + αȳ is in the set S1 which contradicts the assumption that S1 is bounded.

(ii): Since Sγ is bounded for any γ ∈ [0, 1], by Theorem 2.2.8 of [11], the function f (x) has a minimizerz(γ) ∈ Sγ . Since Karush-Kuhn-Tucker optimality conditions imply that z(γ) satisfies
g(z(γ)) = A∗ y,
it follows from the system of (11) and (12) that n̄(z(γ)) = 0.
(iii): Since the function f is convex, it holds thatf (x) ≥≥f (x0 ) + �g(x0 ), x − x0 �
f (x0 ) − �g(x0 )��x − x0 �

for any x in the bounded convex set ∪γ∈[0,1] Sγ . Thus the function f (x) is bounded from below on∪γ∈[0,1] Sγ .
6

The assertion (ii) of Lemma 3.2 ensures that the value �n̄(x)�x can be used to measure proximity to theanalytical center. Moreover, it follows from (iii) of Lemma 3.2 that the value∆ := f (x0 ) − inf{f (x) | x ∈ ∪γ∈[0,1] Sγ } ≥ 0(15)

is well defined. It is not required to compute the value ∆ in our algorithm, but plays an important roleto estimate the iteration number of centering steps described below.
Let x be the current iterate of our algorithm and define r := Ax − b. In our algorithm, the residual
vector r is decreasing on a line segment (0, r0 ] whenever r0 = Ax0 − b �= 0 by using the shifted projected
Newton direction given below (see (20)). So, we can take a γ ∈ (0, 1] so that r = γr0 .
If x is close to the analytical center z(γ), e.g., �n̄(x)�x ≤ 1 and γ is suﬃciently small, i.e., γ ≤ �, ouralgorithm will stop with a suﬃciently approximated feasible interior point x. We may start several
primal interior point algorithm with x as an initial point.
If �n̄(x)�x ≤ 1 but γ > �, our algorithm will proceed to the feasibility step where we compute the shiftedprojected Newton direction n̄r (x) which is the solution of�
H(x)n̄r (x) + g(x) = A∗ y,(16)
An̄r (x) = −r.
DefineThen it holds thatdr (x) := −H(x)−1 A∗ (AH(x)−1 A∗ )−1 r.

(17)n̄r (x) = n̄(x) + dr (x).

(18)

We calculate the new iteration x+ as
x+ := x + αn̄r (x), α :=1−τ∈ (0, 1 − τ ].
1 + �dr (x)�x(19)

Note that by (16), x+ satisfies
Ax+ − b = A(x + αn̄r (x)) − b = (Ax − b) + αAn̄r (x) = (1 − α)r.

(20)

The above relation and the fact α ∈ (0, 1) ensure that the residual vector r is decreasing on the line
segment (0, r0 ].
If �n̄(x)�x > 1 our algorithm will proceed to the centering steps described as in Theorem 2.2.3 of [8]until we find a point y such that �n̄(y)�x ≤ 1.
We summarize our algorithm below.
A primal barrier function Phase I algorithm
Step 0: Choose an accuracy � > 0, a step size parameter τ ∈ (0, 1), and an initial point x0 ∈ int K. Set
x := x0 .
Step 1: Define r := Ax − b and let γ ∈ (0, 1] be the value satisfying r = γr0 . If�n̄(x)�x ≤ 1, γ ≤ �then stop. We have obtained a suﬃciently approximated feasible interior point x. Here �n̄(x)�x
can be calculated by solving the system (11) – (12) and by (6).
7

Step 2: If�n̄(x)�x ≤ 1, γ > �then proceed to the feasibility step. Compute the shifted projected Newton direction n̄r (x) bysolving (16). Set the new iterate x+ as in (19). Set x := x+ . Go to Step 1.
Step 3: If�n̄(x)�x > 1

then proceed to the centering steps described as in Theorem 2.2.3 of [8] until we find a point y ∈ Sγ
such that �n̄(y)�x ≤ 1. Set x := y. Go to Step 1.
Note that at the beginning of the centering steps, we have an interior initial point x ∈ Sγ to find apoint y ∈ Sγ which is suﬃciently close to the analytical center z(γ) of Sγ . This fact implies that thecentering steps in Step 3 exactly corresponds to the ones in existing primal feasible algorithms as in[11] or in [8] and has been well analyzed. Thus we only have to know the behavior of the feasibility step.
In the succeeding sections, we will analyze the feasibility step and give a suﬃcient condition for finding
an approximate feasible interior point in a finite number of iterations.

4

Analysis of the feasibility stepIn this section, we analyze the feasibility step in Section 3. Recall that K is a closed convex cone and f
is a self-concordant function over Df := int K with the complexity value ϑf . At the feasibility step, wechoose a point x+ as in (19). As we will see in Theorem 4.6, the value of the barrier function f will be
increased at most O(ϑf ) at x+ .
Lemma 4.1. At the feasibility step, we have(i)α�n̄r (x)�x ≤ α (�n̄(x)�x + �dr (x)�x ) ≤ 1 − τ.

(ii)
x+

:= x + α (n̄(x) + dr (x)) ∈ Bx (x, 1) ⊆ int K,
x + αn̄(x) ∈ Bx (x, 1) ⊆ int K.

(iii)
0 < α�dr (x)�x+αn̄(x) ≤(iv)

α�dr (x)�x≤ 1 − τ.
1 − α�n̄(x)�xx+ ∈ Bx+αn̄(x) (x + αn̄(x), 1) ⊆ int K.

Proof. (i): The first inequality follows from (18). Since �n̄(x)�x ≤ 1 holds at Step 2, by the definition(19) of α, we haveα (�n̄(x)�x + �dr (x)�x ) ≤ α (1 + �dr (x)�x ) ≤ 1 − τ.
(ii): By the definitions (19) and (18), we see that�x+ − x�x = α�n̄r (x)�x = α�n̄(x) + dr (x)� ≤ α (�n̄(x)�x + �dr (x)�x ) .
8Thus, the assertion (i) and the property (7) imply that x+ ∈ Bx (x, 1) ⊆ Df = int K. By a similar
discussion, we also see thatα�n̄(x)�x ≤ 1 − τ and x + αn̄(x) ∈ Bx (x, 1) ⊆ Df = int K.
(iii): Note that �dr (x)� > 0 whenever the algorithm enters the feasibility step. Since x + αn̄(x) ∈ Df =int K by (ii) above, by substituting x + αn̄(x) and dr (x) �= 0 into y and v in (8), respectively, we seethat�dr (x)�x+αn̄(x)1≤r�d (x)�x
1 − α�n̄(x)�x
or equivalently

0 < �dr (x)�x+αn̄(x) ≤holds. Since the assertion (i) ensures that

�dr (x)�x
1 − α�n̄(x)�xα�dr (x)�x≤ 1 − τ,1 − �αn̄(x)�x
we obtain the assertion (iii).
(iv): We see that�x+ − (x + αn̄(x))�x+αn̄(x) = �αdr (x)�x+αn̄(x) .

Thus the assertion follows from (iii) and the property (7).

We will derive an upper bound of f (x+ ) using the above and the results in Section 2.
First we observe the relation between f (x+ ) and f (x + αn̄(x)). Since x+ ∈ Bx+αn̄(x) (x + αn̄(x), 1) by(iv) of Lemma 4.1, by substituting x+ := x + αn̄(x) + αdr (x) and x + αn̄(x) into y and x in Theorem
2.4, we obtain the following inequality:f (x+ ) ≤f (x + αn̄(x))+�g(x + αn̄(x)), αdr (x)�1
+ �αdr (x), H(x + αn̄(x))αdr (x)�2
α3 �dr (x)�3x+αn̄(x)+.
3(1 − α�dr (x)�x+αn̄(x) )We derive the following bounds for each term in (21) - (24).
Lemma 4.2 (An upperbound of the term in (21)).
f (x + αn̄(x)) ≤ f (x) +9(1 − τ )2
τ

(21)
(22)
(23)
(24)

Proof. By (ii) of Lemma 4.1, we know that x + αn̄(x) ∈ Bx (x, 1). Thus by Theorem 2.4, it holds thatf (x + αn̄(x))

1α3 �n̄(x)�3x≤ f (x) + �g(x), αn̄(x)� + �αn̄(x), H(x)αn̄(x)� +2
3(1 − α�n̄(x)�x )331α�n̄(x)�
x= f (x) − �n(x), αn̄(x)�x + �αn̄(x)�2x +23(1 − α�n̄(x)�x )(by (6) and (10))
1α3 �n̄(x)�3x≤ f (x) − α�n̄(x)�2x + α2 �n̄(x)�2x +23(1 − α�n̄(x)�x )(by (13))1
α3 �n̄(x)�3x≤ f (x) + α2 �n̄(x)�2x +.
23(1 − α�n̄(x)�x )Note that at the feasibility step, we have �n̄(x)�x ≤ 1, α ≤ 1 − τ by (19), and
1
1
1
≤
≤
r
1 − α�n̄(x)�x1 − α (�n̄(x)�x + �d (x)�x )τ
by (i) of Lemma 4.1. These facts and τ ∈ (0, 1) imply that1 2
α3 �n̄(x)�3x
α �n̄(x)�2x +23(1 − α�n̄(x)�x )≤<<(1 − τ )2(1 − τ )3+23τ
(1 − τ )2(1 − τ )2+2τ
3τ
(1 − τ )2.
τ

Thus we obtain the lemma.
Lemma 4.3 (An upperbound of the term in (22)).
�g(x + αn̄(x)), αdr (x)� ≤ ϑf +(1 − τ )2.
τ

Proof. In what follows, we denote H(x)−1 g(x) by gx (x) for ease of notation. It follows from (10) that10

gx (x) = −n(x) and �gx (x)�x = �n(x)�x . Then, by the definition (6), we have�g(x + αn̄(x)), αdr (x)�= �gx (x + αn̄(x)), αdr (x)�x≤ �αdr (x)�x �gx (x + αn̄(x))�x
� �
��
� 1
�
�
r
�
= �αd (x)�x �gx (x) + αn̄(x) + αn̄(x)(Hx (x + tαn̄(x)) − I)dt��0x(by Proposition 2.1)��� 1≤ �αdr (x)�x �gx (x)�x + �αn̄(x)�x + �αn̄(x)�x
�Hx (x + tαn̄(x)) − I�x dt
0

�
� 1�
r
≤ �αd (x)�x �gx (x)�x + �αn̄(x)�x + �αn̄(x)�x
0

=
=
=
=

(by Theorem 2.3)�� �
1
− 1 dt(1 − t�αn̄(x)�x )2

�1 �
1
�αd (x)�x �n(x)�x + �αn̄(x)�x + �αn̄(x)�x−t�αn̄(x)�x − �αn̄(x)�2x t0��2�αn̄(x)�x�αdr (x)�x �n(x)�x + �αn̄(x)�x +1 − �αn̄(x)�x���αn̄(x)�xr
�αd (x)�x �n(x)�x +1 − �αn̄(x)�x��α�n̄(x)�xrα�d (x)�x �n(x)�x +.
1 − α�n̄(x)�x�

r

Note that at the feasibility step, we have α�n̄(x)�x ≤ 1 − τ by (19), andα�dr (x)�x ≤ 1 − τ,
1
1
1
≤
≤
r
1 − α�n̄(x)�x1 − α (�n̄(x)�x + �d (x)�x )τ
by (i) of Lemma 4.1. The above facts ensure that��α�n̄(x)�xrα�d (x)�x �n(x)�x +≤1 − α�n̄(x)�x
≤

�
�
1−τ
(1 − τ ) ϑf +
τ(1 − τ )2ϑf +.
τ

and we obtain the lemma.
Lemma 4.4 (An upperbound of the term in (23) ).
1(1 − τ )2�αdr (x), H(x + αn̄(x))αdr (x)� ≤
.
2
2τ 2

11

Proof. Recall that Hx (y) is defined by (9). By the definition (6), we see that1�αdr (x), H(x + αn̄(x))αdr (x)�21=
�αdr (x), Hx (x + αn̄(x))αdr (x)�x21≤
�αdr (x)�2x �Hx (x + αn̄(x))�x211≤�αdr (x)�2x
2(1 − �αn̄(x)�x )2��2
1α�dr (x)�x=
2 1 − α�n̄(x)�x
where the last inequality follows from x + αn̄(x) ∈ Bx (x, 1) by (ii) of Lemma 4.1 and from Theorem 2.3.
As we have seen in the proof of Lemma 4.3, it holds thatα�dr (x)�x ≤ 1 − τ and
Therefore, we have

which completes the proof.

�α�dr (x)�x1 − α�n̄(x)�x11≤ .
1 − α�n̄(x)�x
τ

�2

≤

�

1−τ
τ

�2

Lemma 4.5 (An upperbound of the term in (24)).

α3 �dr (x)�3x+αn̄(x)3(1 − α�dr (x)�x+αn̄(x) )≤(1 − τ )3.
3τ

Proof. The definition of the local norm (6) and the assertion (iii) of Lemma 4.1 ensure thatα�dr (x)�x+αn̄(x) ≤ 1 − τ.
Thus, it immediately follows thatα3 �dr (x)�3x+αn̄(x)3(1 − α�dr (x)�x+αn̄(x) )≤(1 − τ )33τ

which completes the proof.
Combining the results in Lemmas 4.2 to 4.5, we obtain the following theorem which gives an upperboundof the value f (x+ ).
Theorem 4.6. At the feasibility step, the new iterate x+ ∈ Df = int K satisfies
f (x+ ) ≤ f (x) + ϑf + 3
Thus

�

1−τ
τ

�2

f (x+ ) − f (x) = O(1 + ϑf ).
12.

Proof. By replacing the terms in (21) – (24) by the upperbounds in Lemmas 4.2 to 4.5, we havef (x+ ) ≤≤(1 − τ )2τ(1 − τ )2+ϑf +τ
(1 − τ )2+2τ 2(1 − τ )3
+
3τ
�
�2
1−τ
f (x) + ϑf + 3.
τ
f (x) +

Since τ ∈ (0, 1), we obtain the theorem.

5

Finite termination of the algorithm

In this section, we analyze the total number of steps of the primal barrier function Phase I algorithm
proposed in Section 3. We impose the following assumption.
Assumption 5.1. There exists δ < +∞ such that��sup �dr (x)�x | x ∈ ∪γ∈[0,1] Sγ , �n̄(x)�x ≤ 1 ≤ δ.

Here � · �x , dr (x), Sγ and n̄(x) are defined by (6), (17), (14) and (11) – (12), respectively.
Note that by the definitions � · �x and dr (x), �dr (x)�x can be represented equivalently as���dr (x)�x = �dr (x), H(x)dr (x)� = �r, (AH(x)−1 A∗ )−1 r�.
In what follows, we show that Assumption 5.1 is a suﬃcient condition for finite termination of thealgorithm. We obtain the following lemma under Assumption 5.1.
Lemma 5.2. Suppose that Assumptions 3.1 and 5.1 hold. Then the maximum number of executions ofStep 2 (the feasibility step) is bounded by
�
�
1+δ
1
1
log = O (1 + δ) log
.
1−τ
�
�
Proof. We have seen that (20) holds at the feasibility step. It follows from the definition (19) that
Assumption 5.1 ensures that
1−τ
α≥
.
(25)
1+δ
Note that the value of γ is not changed during centering steps. Thus, at the N -th iteration xN of thefeasibility step, we have
�
�N
1−τ
N
N 0
N
Ax − b = γ r , γ ≤ 1 −
.
1+δ
13

By the construction of the algorithm, we will stop if γ N ≤ �. Thus if
�
�N
1−τ
1−
≤�
1+δ

(26)

then the algorithm will be terminated. By taking the logarithm of both sides of (26) and using the factlog(1 − ξ) ≤ −ξ for any ξ ∈ [0, 1)
we obtain
N
Thus we complete the proof.

�

1−τ
1+δ

�

1
≥ log .
�The above lemma ensures that if Assumptions 3.1 and 5.1 hold then the original problem has an approximated feasible interior point. Using the above lemma, we obtain the following iteration bound of
the algorithm.
Theorem 5.3. Suppose that Assumptions 3.1 and 5.1 hold. Then the total number of steps of the primal
barrier function Phase I algorithm is
�
�
�2 �
1
1
O ∆(1 + δ) log + (1 + ϑf ) (1 + δ) log
��
where ∆ is defined by (15).
Proof. It is well known that the number of the centering steps bounded by O(f 0 − f ∗ ) where f 0 is thevalue of f at the beginning of Step 3 and f ∗ is the value of f at the analytical center of f ( see Theorem
2.2.3 of [8] or Section 3 of [9]).
In our algorithm, the value of the function f is increased only at Step 2 of the feasibility step. As wehave shown in Theorem 4.6, the increment is O(1 + ϑf ) at each feasibility step. Moreover, Lemma 5.2
ensures that the total increment throughout the algorithm is��1
O (1 + ϑf )(1 + δ) log.
�
Thus, the function value f 0 at the beginning of Step 3 is
�
�
1
0
0
f = f (x ) + O (1 + ϑf )(1 + δ) log.
�
At each Step 3, the value f (z(γ)) at the analytical center z(γ) ∈ Sγ satisfies
f (z(γ)) ≥ inf{f (x) | x ∈ ∪γ∈[0,1] Sγ }where the set Sγ is defined by (14). Thus, by the definition (15), we have��10
f − f (z(γ)) ≤ ∆ + O (1 + ϑf )(1 + δ) log
�

14which ensures that the number of corrector steps in each Step 3 is��1O ∆ + (1 + ϑf )(1 + δ) log.
�
Let N be the maximum number of executions of Step 2. Here the maximum number of executions ofStep 3 is bounded by N + 1 since Step 3 is never executed consecutively by the construction of thealgorithm.
Thus, by Lemma 5.2 again, the total number of centering steps throughout the algorithm is
�
�
�2 �
1
1
O ∆(1 + δ) log + (1 + ϑf ) (1 + δ) log��
which completes the proof.
By Proposition 2.6, we obtain the following corollary for the special case where K is the DNN cone Dn .
Corollary 5.4. Let K be the DNN cone Dn given by (2). Suppose that Assumptions 3.1 and 5.1 hold.
Then the total number of steps of the primal barrier function Phase I algorithm is
�
�
�2 �
1
1
2
O ∆(1 + δ) log + n (1 + δ) log
��
where ∆ is defined by (15).

6Concluding remarks

In this paper, we proposed a purely primal barrier function Phase I algorithm for solving nonsymmetricconic optimization problems including the DNN optimization problem. In order to develop eﬃcient andpractical algorithms for such problems, many issues remain to be investigated.
The first one is to find much clearer justification for finite termination of the algorithm. Assumption
5.1 gives a suﬃcient condition, but we have to figure out which instances satisfy the assumption. The
second one is to provide more sophisticated algorithms. In our algorithm, we only use the self-concordantbarrier function f and do not make any use of the objective function �c, x�. Using the hybrid functionη�c, x� + f (x)instead of f will be a promising research direction for developing practical primal barrier algorithms. It
is of course important to conduct numerical experiments to observe the behavior of the algorithms.

References[1] R.E. Burkard, S.E. Karisch and F. Rendl, QAPLIB: A Quadratic Assignment Problem Library,
http://www.seas.upenn.edu/qaplib/15[2] 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, vol. 18, pp.301-320,
2000.
[3] S. Burer, On the copositive representation of binary and continuous nonconvex quadratic programs,
Mathematical Programming, vol. 120, pp.479-495, 2009.
[4] E. De Klerk and D. V. Pasechnik, Approximation of the stability number of a graph via copositive
programming, SIAM Journal on Optimization, vol. 12, pp.875-892, 2002.
[5] E. De Klerk, D. V. Pasechnik, R. Sotirov, On semidefinite programming relaxations of the travelingsalesman problem, SIAM Journal on Optimization vol. 19, pp.1559-1573, 2008.
[6] O. Güler, Hyperbolic polynomials and interior point methods for convex programming, Mathematics
of Operations Research, vol. 22, pp.350–377, 1997.
[7] G. Iyengar, D.J. Phillips and C. Stein, Feasible and accurate algorithms for covering semidefinite
programs, H. Kaplan (Ed.): SWAT 2010, LNCS 6139, pp.150–162, 2010.
[8] Yu.E. Neasterov and A.S. Nemirovskii, Interior-Point Polynomial Algorithm in Convex Programming, SIAM, Philadelphia, 1994.
[9] Yu.E. Nesterov, Towards nonsymmetric conic optimization, CORE Discussion Paper No. 2006/28,
2006.
[10] J. Povh and F. Rendl, Copositive and semidefinite relaxations of the quadratic assignment problem,
Discrete Optimization vol. 16, pp.231–241, 2009.
[11] J. Renegar, A Mathematical View of Interior Point Methods for Convex Optimization, MPS-SIAM
Series on Optimization, 2001.
[12] J. Renegar, Hyperbolic programs, and their derivative relaxations, Foundations of Computational
Mathematics vol. 6, pp.59-79, 2005.
[13] SDPA Online Solver. http://sdpa.indsys.chuo-u.ac.jp/portal/[14] A. Skajaa, J. B. Jørgensen and P. C. Hansen, On implementing a homogeneous interior-point
algorithm for nonsymmetric conic optimization, IMM-Technical Report 2011-02, 2011.
[15] J. F. Sturm and S. Zhang. On cones of nonnegative quadratic functions. Mathematics of Operations
Research, vol. 28, 2003, pp.246-267.
[16] A. Yoshise and Y. Matsukawa, On optimization over the doubly nonnegative cone, Proceedings of2010 IEEE Multi-conference on Systems and Control, pp.13-19, 2010.

16