c The Operations Research Society of Japan⃝

Journal of the Operations Research Society of Japan
Vol. 58, No. 1, January 2015, pp. 24–60

A SURVEY OF NUMERICAL METHODS FOR NONLINEARSEMIDEFINITE PROGRAMMING

Hiroshi Yamashita
NTT DATA Mathematical Systems Inc.

Hiroshi Yabe
Tokyo University of Science(Received September 16, 2014; Revised December 22, 2014)
Abstract Nonlinear semideﬁnite programming (SDP) problems have received a lot of attentions because oflarge variety of applications. In this paper, we survey numerical methods for solving nonlinear SDP problems.
Three kinds of typical numerical methods are described; augmented Lagrangian methods, sequential SDP
methods and primal-dual interior point methods. We describe their typical algorithmic forms and discuss
their global and local convergence properties which include rate of convergence.

Keywords: Optimization, nonlinear semideﬁnite programming, augmented Lagrangian
method, sequential SDP method, primal-dual interior point method

1. Introduction
This paper is concerned with the nonlinear SDP problem:
minimizesubject tof (x),x ∈ Rn ,g(x) = 0, X(x) ≽ 0(1.1)

where the functions f : Rn → R, g : Rn → Rm and X : Rn → Sp are suﬃciently smooth,and Sp denotes the set of pth-order real symmetric matrices. We also deﬁne Sp+ to denotethe set of pth-order symmetric positive semideﬁnite matrices. By X(x) ≽ 0 and X(x) ≻ 0,
we mean that the matrix X(x) is positive semideﬁnite and positive deﬁnite, respectively.
When f is a convex function, X is a concave function and g is aﬃne, this problem is aconvex optimization problem. Here the concavity of X(x) means thatX(λu + (1 − λ)v) − λX(u) − (1 − λ)X(v) ≽ 0
holds for any u, v ∈ Rn and any λ satisfying 0 ≤ λ ≤ 1. We note that problem (1.1) isan extension of a linear SDP problem. When all the functions f and g are linear and thematrix X(x) is deﬁned byn∑X(x) =
xi A i − B
i=1

with given matrices Ai ∈ Sp , i = 1, . . . , n, and B ∈ Sp , the problem reduces to a linear SDP
problem. The linear SDP problems include linear programming problems, convex quadratic
programming problems and second-order cone programming problems [1] as special cases.
Linear SDP model has been one of the most active research ﬁeld for several decades. There
are many researches on theories and applications, and polynomial-time algorithms based oninterior point methods for linear SDP problems. These results can be found in survey papers
by Vandenberghe and Boyd [68] and Todd [66], and in the books by Boyd, Ghaoui, Feron

24

A Survey of Numerical Methods for Nonlinear SDP

25

and Balakrishnan [10], Wolkowicz, Saigal and Vandenberghe [70], Ben-Tal and Nemirovski[4], and Anjos and Lasserre [3], for example.
Though the linear SDP model is very useful in practical applications, it is insuﬃcient
if one wants to deal with more general problems. The nonlinear SDP problems arise fromseveral application ﬁelds, for example, control theory (especially LMI-constrained problems
and BMI-constrained problems), structural optimization, material optimization, eigenvalue
problems, ﬁnance and so forth. See [17, 18, 23, 34, 40, 43, 57, 63, 69] and references therein.
Thus, it is desired to develop numerical methods for solving nonlinear SDP problems.
Though the interior point methods are main tools for linear SDP problems, nonlinear SDP
problems (1.1) can have various algorithms. Typical studies on numerical methods include
the three categories; (1) augmented Lagrangian method, (2) the sequential linear programming (SLP) method, or the sequential quadratic programming (SQP) method, and (3)interior point methods.
The nonlinear SDP contains the nonlinear second-order cone programming (SOCP).
There are several researches on numerical methods for nonlinear SOCP problems, but these
topics are not included in this paper. We list Kanzow, Ferenczi and Fukushima [35], Yamashita and Yabe [74] and Fukuda, Silva and Fukushima [24] for references.
The present paper is organized as follows. In Section 2, we introduce optimality conditions for problem (1.1) and fundamental notions that are used in the subsequent sections.
In Section 3, we review the augmented Lagrangian method and its convergence properties.
Section 4 describes sequential SDP methods from the view point of the rate of local convergence and the global convergence properties within a framework of a line search strategy, thetrust region strategy and the ﬁlter method. In Section 5, we focus on primal-dual interior
point methods and their local and global convergence properties. Since researches by thecurrent authors in nonlinear SDP area have been on primal-dual interior point methods, thedescription of this section may be more detailed than the other methods. Finally, we givesome concluding remarks in Section 6.
Notations: Throughout this paper, we deﬁne the inner product ⟨U, V ⟩ by ⟨U, V ⟩ =tr(U V ) for any matrices U and V in Sp , where tr(M ) denotes the trace of the matrix M .
The superscript T denotes the transpose of a vector or a matrix. For U, V ∈ Sp , we deﬁne
the multiplication U ◦ V by
UV + V U
U ◦V =.
(1.2)
2
We will implicitly make use of various useful relations described in [2] and Appendix of [67].
For P ∈ Rp1 ×p2 and Q ∈ Rq1 ×q2 , the Kronecker product is deﬁned by
P ⊗ Q = [Pij Q] ∈ Rp1 q1 ×p2 q2 .
For U ∈ Rp2 ×q2 , we have(P ⊗ Q)vec(U ) = vec(P U QT ),where the notation vec(U ) is deﬁned byvec(U ) = (U11 , U21 , . . . , Up2 1 , U12 , U22 , . . . , Up2 2 , U13 , . . . , Up2 q2 )T ∈ Rp2 q2 .
For U ∈ Sp , P ∈ Rp×p and Q ∈ Rp×p , we deﬁne the operator
1(P ⊙ Q)U = (P U QT + QU P T )2c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝26

H. Yamashita & H. Yabeand the symmetrized Kronecker product(P ⊗S Q)svec(U ) = svec((P ⊙ Q)U ),where the operator svec is deﬁned by√√√√svec(U ) = (U11 , 2U21 , . . . , 2Up1 , U22 , 2U32 , . . . , 2Up2 , U33 , . . . , Upp )T ∈ Rp(p+1)/2 .
We note that, for any U, V ∈ Sp ,⟨U, V ⟩ = tr(U V ) = svec(U )T svec(V )(1.3)

holds, and that for any U, V ∈ Rp1 ×p2 ,tr(U T V ) = vec(U )T vec(V ).

(1.4)

In what follows, ∥ • ∥, ∥ • ∥1 and ∥ • ∥∞ denote the l2 , l1 and l∞ norms for vectors, and∥ • ∥F denotes the Frobenius norm for matrices. λmin (M ) denotes the minimum eigenvalueof a matrix M ∈ Sp .
2. Optimality Conditions and Preliminaries for Analysis of Local Behavior
This section introduces optimality conditions for problem (1.1) and related quantities. We
ﬁrst deﬁne the Lagrangian function of problem (1.1) by
L(w) = f (x) − y T g(x) − ⟨X(x), Z⟩ ,where w = (x, y, Z), and y ∈ Rm and Z ∈ Sp are the Lagrange multiplier vector andmatrix for the equality and positive semideﬁniteness constraints, respectively. We alsodeﬁne matrices∂X(x)Ai (x) =∂xi
for i = 1, . . . , n. Then the Karush-Kuhn-Tucker (KKT) conditions for optimality of problem(1.1) are given by the following (see [11]):  ∇x L(w)
0



0 
g(x)
=
(2.1)
r0 (w) ≡
0
X(x)Z
and
X(x) ≽ 0,

Z ≽ 0.

Here ∇x L(w) is the gradient vector of the Lagrangian function given by∇x L(w) = ∇f (x) − ∇g(x)y − A∗ (x)Z,
where ∇g(x) is deﬁned by∇g(x) = (∇g1 (x), . . . , ∇gm (x)) ∈ Rn×m
and A∗ (x) is the operator such that for Z,⟨A1 (x), Z⟩..
A∗ (x)Z = .
.
⟨An (x), Z⟩

c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝(2.2)

A Survey of Numerical Methods for Nonlinear SDP

27

We will use the norm ∥r0 (w)∥ deﬁned by√°()°° ∇x L(w) °2°° + ∥X(x)Z∥2∥r0 (w)∥ = °F°g(x)in this paper.
The complementarity condition X(x)Z = 0 will appear in various forms in the following.
We will occasionally deal with the multiplication X(x) ◦ Z instead of X(x)Z. It is known
that X(x) ◦ Z = 0 is equivalent to the relation X(x)Z = ZX(x) = 0 for symmetric positive
semideﬁnite matrices X(x) and Z. We also note that for symmetric positive semideﬁnite
matrices X(x) and Z, X(x)Z = 0 ⇐⇒ ⟨X(x), Z⟩ = 0.
In the rest of this section, we brieﬂy present some deﬁnitions that are necessary forthe analysis of local behavior of methods surveyed below. We also describe the deﬁnitionsof a stationary point, the Mangasarian-Fromovitz constraint qualiﬁcation condition, thequadratic growth condition, the strict complementarity condition and the nondegeneracy
condition, and the second order necessary / suﬃcient conditions for optimality. More comprehensive description can be found in [7–9, 56, 59, 60]. We recommend the paper by Shapiro
and Scheinberg [60] for a good introduction to these subjects. Alternative derivation of the
optimality conditions is described in Forsgren [21].
Deﬁnition 2.1. A point x∗ ∈ Rn is said to be a stationary point of problem (1.1) if there
exist Lagrange multipliers (y, Z) such that (x∗ , y, Z) satisﬁes the KKT conditions (2.1) and(2.2).
Let Λ(x∗ ) denote the set of Lagrange multipliers (y, Z) such that (x∗ , y, Z) satisﬁes the
KKT conditions. In this paper, when we refer to a point w∗ = (x∗ , y ∗ , Z ∗ ), then it meansthat w∗ is a KKT point.
Deﬁnition 2.2. We say that the Mangasarian-Fromovitz constraint qualiﬁcation (MFCQ)
condition holds at a feasible point x if the matrix ∇g(x) is of full column rank and there
exists a nonzero vector v ∈ Rn such that∇g(x)T v = 0

and X(x) +n∑vi Ai (x) ≻ 0.

i=1

It can be shown that, if the Mangasarian-Fromovitz condition holds at a stationary point
x , then the set Λ(x∗ ) is bounded.
The set C(x∗ ) denote the critical cone of (1.1) at x∗ that is deﬁned as follows. Let
TSp+ (X(x∗ )) denote the tangent cone of Sp+ at X(x∗ ), which is deﬁned by∗

TSp+ (X(x∗ )) = {D | dist(X(x∗ ) + tD, Sp+ ) = o(t), t ≥ 0},
where dist(P, Sp+ ) = inf{∥P − Q∥F | Q ∈ Sp+ }. The set C(x∗ ), the critical cone at x∗ , isdeﬁned by{}n∑C(x∗ ) = h ∈ Rn | ∇g(x∗ )T h = 0,hi Ai (x∗ ) ∈ TSp+ (X(x∗ )), ∇f (x∗ )T h = 0 .
i=1Deﬁnition 2.3. The second order necessary condition for local optimality of x∗ under the
MFCQ condition is given by
sup(y,Z)∈Λ(x∗ )hT (∇2x L(x∗ , y, Z) + Ĥ(x∗ , Z))h ≥ 0, for all h ∈ C(x∗ ).

c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝

(2.3)

28

H. Yamashita & H. Yabe

Here Ĥ(x, Z) is a matrix whose (i, j)th element is(Ĥ(x, Z))ij = 2tr(Ai (x)X(x)† Aj (x)Z)(2.4)

and † denotes the Moore-Penrose generalized inverse.
The matrix Ĥ(x, Z) contains a curvature information of the domain deﬁned by theconstraint X(x) ≽ 0.
Corresponding to the above necessary condition, we have the following suﬃcient condition.
Deﬁnition 2.4. The second order suﬃcient condition for local optimality of x∗ under the
MFCQ condition is given by
sup(y,Z)∈Λ(x∗ )hT (∇2x L(x∗ , y, Z) + Ĥ(x∗ , Z))h > 0, for all h ∈ C(x∗ )\{0}.

(2.5)

It is said that a quadratic growth condition holds at a feasible point x∗ of problem (1.1)
if there exists c > 0 such that the following inequality holds
f (x) ≥ f (x∗ ) + c∥x − x∗ ∥2(2.6)

for any feasible point x in a neighborhood of x∗ . The quadratic growth condition implies thatx∗ is a strict local optimal solution of problem (1.1). Suppose that the MFCQ condition
holds. Then the quadratic growth condition holds if and only if the above second order
suﬃcient conditions (2.5) are satisﬁed. By comparing conditions (2.3) and (2.5), we notethat these conditions give a pair of “no gap” second order optimality conditions. By “nogap”, we mean that the weak inequality sign in (2.3) is changed to the strict inequality signin (2.5).
The matrix Ĥ(x∗ , Z) can be written in the following form()
Ĥ(x∗ , Z) = 2A(x∗ )T Z ⊗ X(x∗ )† A(x∗ ),whereA(x) = [vec(A1 (x)), . . . , vec(An (x))] ∈ Rp ×n .
2

(2.7)

Since X(x∗ ) ≽ 0 and Z ≽ 0, the matrix Z ⊗ X(x∗ )† is positive semideﬁnite (see Appendix of[67]), Ĥ(x∗ , Z) is positive semideﬁnite. Then the following second order suﬃcient condition
sup(y,Z)∈Λ(x∗ )hT ∇2x L(x∗ , y, Z)h > 0, for all h ∈ C(x∗ )\{0},(2.8)

which will also be used in the following, is a stronger condition than (2.5). If condition (2.8)
holds, condition (2.5) holds. Then the quadratic growth condition is also satisﬁed under the
MFCQ condition.
Since X(x∗ )Z ∗ = 0, the matrices X(x∗ ) and Z ∗ commute. Therefore they can be simultaneously diagonalized. Thus because of the complementarity condition, we can assume
without loss of generality that the matrix X(x∗ ) and Z ∗ are represented by
 ∗



XB 0 00 0 0
X(x∗ ) =  0 0 0  and Z ∗ =  0 0 0 
0 0 0
0 0 ZN∗
respectively, where XB∗ and ZN∗ are diagonal and positive deﬁnite matrices.
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

29

Deﬁnition 2.5. We say that the strict complementarity condition holds at x∗ if there exists(y ∗ , Z ∗ ) ∈ Λ(x∗ ) such thatrank(X(x∗ )) + rank(Z ∗ ) = p.
If the strict complementarity condition holds, we have
)
( ∗
)
(
XB 0
0 0
∗
∗
X(x ) =and Z =
,
0 0
0 ZN∗

(2.9)

where rank(XB∗ ) + rank(ZN∗ ) = p. Corresponding to (2.9), we partition the matrices X(x)and Z as()()
XB (x) XU (x)ZB (x) ZU (x)X(x) =and Z =XU (x)T XN (x)ZU (x)T ZN (x)in a neighborhood of w∗ = (x∗ , y ∗ , Z ∗ ). Similarly, we partition the matrix Ai (x) as()
ABi (x) AU i (x)Ai (x) =AU i (x)T AN i (x)for i = 1, . . . , n. Then the critical cone at x∗ can be speciﬁcally represented by{}n∑C(x∗ ) = h ∈ Rn | ∇g(x∗ )T h = 0,
hi AN i (x∗ ) = 0 .
i=1Deﬁnition 2.6. We say that the nondegeneracy condition holds at x∗ if the n-dimensional
vectors
eTi A1 (x∗ )ej
..
∇gi (x∗ ), i = 1, . . . , m and  , i, j = 1, . . . , p − r.
T∗ei An (x )ej
are linearly independent, where r = rank(X(x∗ )) and the vectors e1 , . . . , ep−r form a basisof the null space of the matrix X(x∗ ).
If the strict complementarity condition holds, the nondegeneracy condition means thatthe vectors(AN 1 (x∗ ))ij
..
∇gi (x∗ ), i = 1, . . . , m and  , i, j = 1, . . . , |N |.
∗(AN n (x ))ij
are linearly independent, where |N | denotes the size of ZN∗ . If the strict complementarity
condition holds at x∗ , then Λ(x∗ ) is a singleton if and only if the nondegeneracy conditionis satisﬁed. It is known that the nondegeneracy condition is stronger than the MFCQ
condition, i.e., if the nondegeneracy condition holds at x∗ , then the MFCQ condition also
holds at x∗ .
3. Augmented Lagrangian Method
As in ordinary nonlinear programming, it is possible to device a class of algorithms based on
the augmented Lagrangian function. In this section, we survey two kinds of approaches that
use the augmented Lagrangian type merit function. The ﬁrst one is based on a quadratic
penalty type augmentation and the second one is based on Polyak’s modiﬁed barrier function. For other methods which are not reviewed here, refer to [28–31].
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝30

H. Yamashita & H. Yabe

3.1. Augmented Lagrangian method based on quadratic penalty function
To solve the following nonlinear optimization problem:
minimizesubject tof (x),x ∈ Rn ,g(x) = 0, x ≥ 0,

the augmented Lagrangian by the quadratic penalty is of the form (see, for example, [5])
{(
}
{
})2
n
xi
µ∑
1
max 0, zi −Fµ (x, y, z) = f (x) +− zi2 − y T g(x) +∥g(x)∥2 , (3.1)
2 i=1
µ
2µ
where µ > 0 is a penalty parameter. We repeat alternately an unconstrained minimization ofFµ (x, yk , zk ) for computing next primal estimate xk+1 and a multiplier update for computing
next dual estimates yk+1 , zk+1 for k = 0, 1, . . . . Because
0 = ∇x Fµ (xk+1 , yk , zk )}{n∑(xk+1 )i
= ∇f (xk+1 ) −− ∇g(xk+1 )(yk − g(xk+1 )/µ),max 0, (zk )i −µi=1
if xk+1 is at the minimum of Fµ (x, yk , zk ), the usual practice is to update the dual variables
byyk+1 = yk − g(xk+1 )/µ,{}(xk+1 )i(zk+1 )i = max 0, (zk )i −, i = 1, . . . , n,µwhich gives ∇x L(xk+1 , yk+1 , zk+1 ) = 0. The update of dual variables corresponds to the
maximization step of the dual function.
An extension of the above algorithm to problem (1.1) is rather straightforward, and theaugmented Lagrangian is deﬁned by
{ (¯
}
¯2 )
¯
¯
µ
1
1
Fµ (x, y, Z) = f (x) +tr ¯¯Z − X(x)¯¯∥g(x)∥2 , (3.2)− tr(Z 2 ) − y T g(x) +2µ2µ+where the matrix |A|+ for A ∈ Sp is deﬁned by |A|+ = Sdiag(|λ1 |+ , . . . , |λp |+ )S T where
A = Sdiag(λ1 , . . . , λp )S T is an eigenvalue decomposition of A, and |λi |+ = max{0, λi }, i =1, . . . , p. Therefore the matrix |A|+ is the projection of A onto Sp+ . We will also use thenotation |A| for A ∈ Sp as the matrix |A| = Sdiag(|λ1 | , . . . , |λp |)S T later.
Since
¯
¯
¯
¯
1
∗
¯
(3.3)
∇x Fµ (x, y, Z) = ∇f (x) − A (x) ¯Z − X(x)¯¯ − ∇g(x)(y − g(x)/µ),µ+and ∇x Fµ (xk+1 , yk , Zk ) = 0 at the minimum xk+1 of Fµ (x, yk , Zk ) (see, for example, [65]),
similarly to the nonlinear optimization case, the next multiplier estimates are computed byyk+1 = yk − g(xk+1 )/µ,
¯
¯
¯
¯
1
Zk+1 = ¯¯Zk − X(xk+1 )¯¯ ,µ+c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

31

yielding ∇x L(xk+1 , yk+1 , Zk+1 ) = 0. We note that Zk+1 is always positive semideﬁnite.
From (3.3), it is easy to verify that ∇x Fµ (x∗ , y ∗ , Z ∗ ) = ∇x L(x∗ , y ∗ , Z ∗ ) = 0 at a KKTpoint w∗ .
In order to calculate the values of the merit function and its derivatives at a point x withthis type of the augmented Lagrangian function, we have to calculate all the eigenvalues ofthe matrix Z − µ1 X(x). This may cause some diﬃculties in practical computation, especiallywith large scale problems.
We note that Noll, Torki and Apkarian [54] solved nonlinear LMI/BMI problems byhandling the equality constraints with an ordinary Lagrangian augmentation as above, andtreating the matrix inequalities by interior point methods.
3.1.1. Local properties of quadratic penalty augmented Lagrangian
In this subsection, we describe local properties of the above augmented Lagrangian near the
KKT point w∗ . The following results are given by Sun, Zhang and Wu [65]. In [65], only
problems with inequality constraints are considered, therefore we omit the term containing
g(x) and y in this subsection. We note that it is a straightforward matter to include
equality constraints in the following analysis. The local convexiﬁcation and global saddle
point condition are proved in the following theorem.
Theorem 3.1. Assume that the functions f and X are twice continuously diﬀerentiable.
Assume further that the Mangasarian-Fromovitz constraint qualiﬁcation, the strict complementarity condition and second order suﬃcient condition (2.8) hold at w∗ . Then theaugmented Lagrangian function Fµ (x, Z) deﬁned by (3.2) satisﬁes the following properties:(i) ∇x Fµ (x∗ , Z ∗ ) = 0.
(ii) v T ∇2x Fµ (x∗ , Z ∗ )v > 0, ∀v ∈ C(x∗ )\{0} when µ > 0. Moreover there is a positive
constant µ̄ > 0 such that ∇2x Fµ (x∗ , Z ∗ ) is positive deﬁnite when µ ∈ (0, µ̄].
(iii) There exists a positive µ0 ≤ µ̄ such that the global saddle point condition
Fµ (x∗ , Z) ≤ Fµ (x∗ , Z ∗ ) ≤ Fµ (x, Z ∗ )for all x ∈ Rn , Z ∈ Sp

holds for all positive µ ≤ µ0 .
Based on these properties, various types of augmented Lagrangian methods can be derived. The following simple algorithm is presented for the analysis of local convergenceproperties.
Algorithm 3.1.
Step 0. Let x0 ∈ Rn , Z0 ≽ 0, µ̄ > µ > 0 and ε > 0 be given. Set k = 0.
Step 1. (Primal update) Solve the following subproblem, and calculate its minimizer xk+1 ∈Rn :minimize Fµ (x, Zk ), x ∈ Rn .
Step 2. If Zk X(xk+1 ) = 0, then stop.
Step 3. (Dual update) Calculate the next multiplier estimate by
¯
¯
¯
¯
1
Zk+1 = ¯¯Zk − X(xk+1 )¯¯ .
µ+Step 4. Set k = k + 1 and go to Step 1.

¤

The local convergence property of the above algorithm is given by the following theorem.
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝32

H. Yamashita & H. Yabe

Theorem 3.2. Assume that the functions f and X are twice continuously diﬀerentiable.
Assume further that the Mangasarian-Fromovitz constraint qualiﬁcation, the strict complementarity condition and second order suﬃcient condition (2.8) hold at w∗ . Then there exist
δ > 0, ϵ > 0, µ̄ > 0, µ1 > 0,and µ2 > 0 with µ1 < µ2 ≤ µ̄, such that for any µ ∈ [µ1 , µ2 ] andfor Z ∈ {Z |∥Z − Z ∗ ∥F ≤ δ }, the following statements hold:(i) There exists a unique vector x̂ ∈ Rn such that ∇Fµ (x̂, Z) = 0 andx̂ = arg min {Fµ (x, Z) |∥x − x∗ ∥ ≤ ϵ, x ∈ Rn } .
(ii) Denote U = Z − µ1 X(x̂) and Ẑ = |U |+ . Assume that the inverse of the matrix)(−A(x̂)T∇2x Fµ (x̂, Ẑ)− 12 A(x̂)T − 12 V −µp2 I
is bounded for Z in a neighborhood of Z ∗ and µ suﬃciently small, where A(x) is deﬁned in(2.7), andV = [vec(V1 ), . . . , vec(Vn )]with

Vi = (|U | ⊙ I)−1 (U ⊙ I)Ai (x̂), i = 1, . . . , n.

Then, the following estimates hold:∥x̂ − x∗ ∥ ≤ cµ ∥Z − Z ∗ ∥F ,
°
°
°
∗°
Ẑ
−
Z
°
° ≤ cµ ∥Z − Z ∗ ∥F
F

where the constant c is independent of µ.
The above theorem shows the rate of convergence is linear, but its speed can be adjusted
by decreasing the value of µ.
We note that some local properties of the quadratic penalty type augmented Lagrangianis studied by Shapiro and Sun [61]. Sun, Sun and Zhang [64] proved the similar linear rateof convergence without assuming the strict complementarity condition.
3.1.2. Globally convergent algorithm
It is possible to globalize the quadratic penalty augmented Lagrangian method. The essenceof the modiﬁcations/assumptions of possible variations is to maintain the primal and dual
sequences as bounded somehow, and let the penalty parameter µk → 0. We describe thealgorithm proposed by Wu, Luo, Ding and Chen [71] in this subsection. We note that Luo,
Wu and Chen [47] proposed other variants of globally convergent algorithms. Let
}
{ °¯
°
¯
°
¯
°¯
1
°
¯
¯
σk ≡ max µk °°¯Zk − µk X(xk+1 )¯ − Zk ° , ∥g(xk+1 )∥∞ .
+F
It can be shown that, if ∥∇Fµk (xk+1 , yk , Zk )∥ ≤ ϵ and σk = 0, then the point (xk+1 , yk , Zk )is an ϵ-approximate KKT point which satisﬁes∥∇f (xk+1 ) − ∇g(xk+1 )yk − A∗ (xk+1 )Zk ∥g(xk+1 )X(xk+1 )Zk
X(xk+1 )

≤
=
=
≽

ϵ,
0,
0,
0, Zk ≽ 0.

Therefore searching for a point that satisﬁes ∥∇Fµk (x, yk , Zk )∥ ≤ ϵ and σk = 0 along withthe semideﬁnite conditions may give a desired approximate KKT point.
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

33

Algorithm 3.2.
Step 0. Let x0 ∈ Rn , y0 ∈ Rm , Z0 ≽ 0, µ0 > 0, µ′0 > 0, ϵ > 0 and θ ∈ (0, 1) be given. Set
k = 0.
Step 1. (Primal update) Find a point xk+1 ∈ Rn that satisﬁes∥∇Fµk (x, yk , Zk )∥ ≤ µ′k ,by solving the problem
minimize Fµk (x, yk , Zk ), x ∈ Rn .
Step 2. If σk = 0 and µ′k ≤ ϵ, then stop. If σk = 0 and µ′k > ϵ, then set µ′k+1 = θµ′k , k =k + 1, and go to Step 1. Otherwise go to Step 3.
Step 3. (Dual update) Calculate the next multiplier estimate by
¯
¯
¯
¯
1
¯
Zk+1 = ¯Zk − X(xk+1 )¯¯ , yk+1 = yk − g(xk+1 )/µk .
µk+
Step 4. Set 0 < µk+1 ≤ µk , k = k + 1 and go to Step 1.

¤

The following theorem shows the convergence of the above algorithm.
Theorem 3.3. Assume that the functions f, g and X are continuously diﬀerentiable.
(i) If the sequence generated by Algorithm 3.2 stops in Step 2 at the kth iteration, then{xk+1 , yk , Zk } is an ϵ-approximate KKT point.
(ii) If the generated sequence is not ﬁnite, and if {yk , Zk } is bounded and µk → 0 as
k → ∞, then any limit point x̄ of {xk } is either infeasible, or does not satisfy the MFCQ
condition, or a KKT point.
In [71], another more complicated algorithm that converges under more relaxed conditions
is also described.
3.2. Augmented Lagrangian function based on modiﬁed barrier function
Representatives of this line were ﬁrst done by Mosheyev and Zibulevsky [49], and later
developed by Kočvara and Stingl ([38–40, 62]). Kočvara and Stingl developed the software
PENNON [38, 41], and applied it to many practical problems. The augmented Lagrangian
(
{
})2
µ
xi
function (3.1) for ordinary NLP contains the penalty term 2 { max 0, zi − µ− zi2 }
for the constraint xi ≥ 0, that goes to inﬁnity as xi goes to −∞. Polyak proposed themodiﬁed barrier function method [55] that uses the penalty/barrier term −µ log( xµi + 1)zi ,for example, for the constraint xi ≥ 0, that goes to inﬁnity as xi goes to −µ < 0 from theabove.
In the following, the outline of the algorithm by Kočvara and Stingl is described. The
problem to be solved isminimizef (x),x ∈ Rn ,(3.4)
subject to X(x) ≽ 0.
Here the equality constraints are omitted for the sake of simplicity of discussion. Kočvara
and Stingl proposed two approaches for handling equality constraints in their implementation. One is that an equality constraint gi (x) = 0 is converted to two inequality constraints
gi (x) ≤ 0 and gi (x) ≥ 0, and the other is that gi (x) = 0 is converted to the inequality
constraints −ε ≤ gi (x) ≤ ε, where ε is a small positive constant. They claimed that these
two approaches are equally eﬃcient.
In this section, we describe various properties of this type following the contents in [62].
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝34

H. Yamashita & H. Yabe

3.2.1. Properties of modiﬁed barrier augmented Lagrangian function
In this subsection, we assume that the KKT point x∗ satisﬁesx∗ ∈ arg min{f (x) |x ∈ Ω },where Ω = {x ∈ Rn |X(x) ≽ 0}, the strict complementarity condition, the nondegeneracy
condition, and the second order suﬃcient condition (2.5), and that∃π > 0 and τ > 0 such that max{∥X(x)∥F |x ∈ Ωπ } ≤ τ,

(3.5)

with
Ωµ = {x ∈ Rn |X(x) ≽ −bµI },where b is a positive constant and µ is a positive number. Condition (3.5) is called a growth
condition.
We introduce a penalty/barrier matrix function Φµ : Sp → Sp that satisﬁes a number of
constraints described below. Let X(x) = S(x)diag(λ1 (x), . . . , λp (x))S(x)T be an eigenvalue
decomposition of X(x). Then the primary matrix function Φµ is deﬁned byφµ (λ1 (x))
0
···
0
..


0
φµ (λ2 (x)).
Φµ (X(x)) = S(x)  S(x)T ,
..
..


.
.
0
0
···
0 φµ (λp (x))
where φµ : R → R characterizes the above matrix function and speciﬁc choice of φµ given
below assures
X(x) ≽ 0 ⇐⇒ Φµ (X(x)) ≽ 0,
in particular, for any value of penalty parameter µ > 0.
The augmented Lagrangian of (3.4) is deﬁned by
Fµ (x, Z) = f (x) − ⟨Φµ (X(x)), Z⟩ .

(3.6)

In order for the augmented Lagrangian (3.6) to have desired properties, we require thefunction φ to satisfy the following properties, and let the above φµ be deﬁned by φµ (t) =−µφ(−t/µ):Deﬁnition 3.1. Let φ : (−∞, b) → R, where b ∈ (0, ∞] is a given number, be a function
with the following properties.
(φ0 ) φ is strictly convex, strictly monotone increasing and twice continuously diﬀerentiable,(φ1 ) φ(0) = 0,
(φ2 ) φ′ (0) = 1,
(φ3 ) ∃c1 such that φ′ (t) ≥ c1 for any t < 0,
(φ4 ) ∃c2 such that φ′ (σ/µ) ≤ c2 µ for any σ < 0 and µ > 0,
(φ5 ) ∃c3 such that φ′′ (ρσ) ≤ c3 µ2 for any σ < 0 and µ > 0,
(φ6 ) φ′ is convex, and limt→b φ′ = ∞, limt→−∞ φ′ = 0.
Well known examples of the function φ are:
• The logarithmic penalty function
φlog (t) = − log(1 − t).
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝35

A Survey of Numerical Methods for Nonlinear SDP• The hyperbolic penalty function
φhyp (t) =

1
− 1.
1−t

These penalty functions are introduced by Polyak [55], and he called them the modiﬁed
barrier functions.
Let t1 , . . . , tp be real values, and deﬁne{
φµ (ti )−φµ (tj ), for i ̸= j,ti −tj∆φµ (ti , tj ) =′for i = j,φµ (ti ), ∆φ (t ,t )−∆φ (t ,t )
µ i k
µ j k

, for i ̸= j,
ti −tj
2
∆φ
(t
,t
)−∆φ
(t
,t
)
µ
µ
i
j
j
k
∆ φµ (ti , tj , tk ) =, for i = j ̸= k,ti −tk ′′φµ (ti ),for i = j = k,where φ′µ and φ′′µ are the ﬁrst and the second derivatives of φµ , respectively. Further,let λ′1 (x) < λ′2 (x) < . . . < λ′p′ (x) (x) be distinct eigenvalues of X(x), and let X(x) =S(x)diag(λ′1 (x), . . . , λ′1 (x), . . . , λ′p′ (x) (x), . . . , λ′p′ (x) (x))S(x)T be an eigenvalue decompositionof X(x) where each eigenvalue occurs in its multiplicity at x. The Frobenius covariant
matrices Pi (x) ∈ Sp , i = 1, . . . , p′ (x) of X(x) is deﬁned by
Pi (x) = S(x)diag(0, . . . , 0, 1, . . . , 1, 0, . . . , 0)S(x)T , i = 1, . . . , p′ (x),where the nonzeros in the diagonal matrix occur exactly in the positions of λ′i (x) indiag(λ′1 (x), . . . , λ′p′ (x) (x)) ∈ Sp .
With the above deﬁnition of the augmented Lagrangian, we can prove the following
theorem that shows the desired properties at the KKT point w∗ .
Theorem 3.4. The augmented Lagrangian function Fµ (x, Z) deﬁned by (3.6) satisﬁes thefollowing properties:(i) Fµ (x∗ , Z ∗ ) = f (x∗ ).
(ii) ∇x Fµ (x∗ , Z ∗ ) = ∇f (x∗ ) − A∗ (x∗ )Z ∗ = 0.
(iii) ∇2x Fµ (x∗ , Z ∗ ) = ∇2x L(x∗ , Z ∗ ) + Hµ (x∗ , Z ∗ ) + M/µ, where
⟨
′ ∗

⟩n
p (x )−1∑Hµ (x∗ , Z ∗ ) = 2  Z ∗ , Ai (x∗ ) ∆2 φµ (tk , 0, 0)Pk (x∗ ) Aj (x∗ ) 
k=1

[⟨
⟩]n
M = 2 Z ∗ , Ai (x∗ )φ′′µ (0)Pp′ (x∗ ) (x∗ )Aj (x∗ ) i,j=1 ∈ Sn ,∈ Sn ,

i,j=1

and limµ→0 Hµ (x∗ , Z ∗ ) = Ĥ(x∗ , Z ∗ ) (see (2.4)), KerM = C(x∗ ) and v T M v > 0, ∀v ∈/C(x∗ ).
If f is convex and X is concave, then(iv) Fµ (x, Z) is convex in x for all x ∈ Ωµ .
From (iii) of the above theorem, we see that the matrix ∇2x Fµ (x∗ , Z ∗ ) is positive deﬁnite
when µ > 0 is suﬃciently small. The following theorem shows properties of the augmented
Lagrangian in a neighborhood of x∗ .

c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝36

H. Yamashita & H. Yabe

Theorem 3.5. The augmented Lagrangian function Fµ (x, Z) satisﬁes the following properties:(i) There exists µ̄ > 0 such that Fµ (x, Z ∗ ) is strongly convex for all µ ≤ µ̄ in a neighborhoodof x∗ .
(ii) There exist ε > 0 and µ̄ > 0 such that, for any µ < µ̄,x∗ = arg min {Fµ (x, Z ∗ ) |∥x − x∗ ∥ ≤ ε} .
Moreover, if f is convex and X is concave, then for any µ > 0,x∗ = arg min {Fµ (x, Z ∗ ) |x ∈ Rn } .
3.2.2. Algorithm and its local convergence property
As usual, we perform the unconstrained minimization of Fµ (x, Z) with ﬁxed multiplierestimate Z and current penalty parameter µ. Then we update the multiplier matrix Z anddecrease the penalty parameter µ if necessary, and repeat this primal and dual iteration
until some termination criterion is satisﬁed. This procedure will be described more precisely
below.
⟨⟩(X(x))
By using the matrix function chain rule, it is possible to write ∂Φµ∂x,Z
as a trace ofia product matrix of Ai (x) and some matrix. We write this latter matrix as DΦµ (X(x)) [Z],and interpret it as a directional derivative of Φµ (X(x)) with respect to X(x) along thedirection Z. More precisely, we have (see (6.6.25) of [27])
′

p (x)
∑
∂Φµ (X(x))∆φµ (λ′j (x), λ′k (x))Pj (x)Ai (x)Pk (x),
=
∂xi
j,k=1

then we have∂Fµ (x, Z)∂f (x)
=
−
∂xi
∂xi

⟨

∂Φµ (X(x))
,Z
∂xi

⟩

′

p (x)
⟩
∂f (x) ∑ ⟨=∆φµ (λ′j (x), λ′k (x))Pj (x)Ai (x)Pk (x), Z
−
∂xi
j,k=1
′

p (x)
⟩
∂f (x) ∑ ⟨=−Ai (x), ∆φµ (λ′j (x), λ′k (x))Pk (x)ZPj (x)
∂xi
j,k=1

=

∂f (x)− ⟨Ai (x), DΦµ (X(x)) [Z]⟩ ,
∂xi

where

p′ (x)DΦµ (X(x)) [Z] =∑∆φµ (λ′j (x), λ′k (x))Pj (x)ZPk (x).

j,k=1

The formal statement of the algorithm is as follows:
Algorithm 3.3.
Step 0. Let x0 ∈ Rn , Z0 ≻ 0, µ0 > 0 and ε > 0 be given. Set k = 0.
Step 1. If∥r0 (wk )∥ ≤ ε, X(xk ) ≽ 0, Zk ≽ 0,then stop.
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

37

Step 2. (Primal update) Solve the following subproblem, and calculate the minimizer
xk+1 ∈ Rn :minimize Fµk (x, Zk ), x ∈ Rn .
Step 3. (Dual update) Calculate the next multiplier estimate by
Zk+1 = DΦµk (X(xk+1 )) [Zk ] .
Step 4. Set µk+1 ≤ µk , k = k + 1 and go to Step 1.

¤

By using the above update formula for Zk , we have∇x Fµ (xk+1 , Zk ) = ∇x L(xk+1 , Zk+1 ) = 0,
as desired. It is also possible to show that if Zk ≻ 0, then Zk+1 ≻ 0 with the above update
formula (see Lemma 7.9 in [62]).
Local convergence properties of the above algorithm is given below. Let |N | be thesize of the matrix ZN∗ . For ϵ and Θ (0 < ϵ < (ZN∗ )ii < Θ, i = 1, . . . , |N |), deﬁne the set
V = V (Z ∗ , µ̄, δ, ϵ, Θ) by
V= {(Z, µ) ∈ Sp+ × R |∥Z − Z ∗ ∥F ≤ δ/µ, µ < µ̄ } ∩{(Z, µ) ∈ Sp+ × R |∥Z∥F ≤ Θ} ∩{(Z, µ) ∈ Sp+ × R |(ZN )ii ≥ ϵ, i = 1, . . . , |N |} .

Theorem 3.6. There exist a penalty parameter µ̄ > 0 and a small enough number δ > 0such that for any (Z, µ) ∈ V ,(i) There exists a vectorx̂ = arg min {Fµ (x, Z) |x ∈ Rn }such that ∇x Fµ (x̂, Z) = 0.
(ii) For the pair x̂ and Ẑ = DΦµ (X(x̂)) [Z], the estimate
°
° }
{
°
∗
∗°
max ∥x̂ − x ∥ , °Ẑ − Z °≤ µc ∥Z − Z ∗ ∥F
Fholds, where c is a constant independent of µ.
(iii) x∗ = arg min {Fµ (x, Z ∗ ) |x ∈ Rn } and Z ∗ = DΦµ (X(x∗ )) [Z ∗ ].
(iv) The function Fµ (x, Z) is strongly convex with respect to x in a neighborhood of x̂.
From the above theorem, the following local convergence result is obtained.
Theorem 3.7. For µ0 (< µ̄) small enough and (Z0 , µ0 ) ∈ V , Algorithm 3.3 converges tow∗ with a linear rate. If µk+1 < µk for all k ≥ 0, and limk→∞ µk = 0, then the rate ofconvergence is superlinear.
3.2.3. Globally convergent algorithm
In [62], a globally convergent augmented Lagrangian algorithm was proposed and its convergence property was proved.
Algorithm 3.4.
Step 0. Let x0 ∈ Rn , Z0 ≻ 0, µ0 > 0, µ′0 > 0 and ε > 0 be given. Set k = 0.
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝38

H. Yamashita & H. Yabe

Step 1. If∥r0 (wk )∥ ≤ ε, X(xk ) ≽ 0, Zk ≽ 0,then stop.
Step 2. (Primal update) Find a point xk+1 ∈ Rn that satisﬁes∥∇Fµk (x, Zk )∥ ≤ µ′k ,by minimizing Fµk (x, Zk ).
Step 3. (Dual update) Calculate the next multiplier estimate by
Zk+1 = DΦµk (X(xk+1 )) [Zk ] .
Step 4. Set µk+1 < µk , µ′k+1 < µ′ , k = k + 1 and go to Step 1.

¤

The global convergence of the above algorithm is shown by the following theorem. As
noted in the previous subsection, this theorem also assumes the boundedness of the sequenceand µk → 0.
Theorem 3.8. Assume that there exists µ̄ > 0 such that Ωµ is a compact set for all µ ≤ µ̄,
the nondegeneracy condition holds at any limit point of the sequence {xk }, and the sequence{Zk } is bounded. Assume further that µk → 0 and µ′k → 0. Then the limit point of any
convergent subsequence of {xk , Zk } satisﬁes the KKT conditions.
3.2.4. Actual form of penalty function
In order to have computable forms of multiplier update and derivatives, Kočvara and Stingl
used the following penalty function:
Φµ (X(x)) = −µ2 (X(x) + µI)−1 + µI.
Let

U (x) = (X(x) + µI)−1 .

Since
]
]
∂ [
∂
∂ [
U (x)U (x)−1 =[U (x)] U (x)−1 + U (x)U (x)−1
∂xi
∂xi
∂xi
∂
[U (x)] U (x)−1 + U (x)Ai (x),
=
∂xi

0 =

we have∂Φµ (X(x))= µ2 U (x)Ai (x)U (x)∂xifor the ﬁrst derivative of Φµ (X(x)). In order to have⟨⟩∂Φµ (X(xk )), Zk = ⟨Ai (xk ), Zk+1 ⟩ ,∂xi
we setZk+1 = µ2k U (xk+1 )Zk U (xk+1 )for multiplier update. We note that when Zk is positive deﬁnite, then Zk+1 is also positivedeﬁnite.
Noll [53] also gave a local convergence property similar to the result given in Section
3.2.2 using this particular form of penalty function.
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝39

A Survey of Numerical Methods for Nonlinear SDP

4. Sequential SDP Method
The sequential SDP (S-SDP) method was proposed by using the idea of the SQP (sequential
quadratic programming) method [6, 52] for solving nonlinear optimization problems. Suchan idea has already been mentioned by Robinson [58] within the more general framework ofnonlinear programming problems over convex cones. Given a current point xk , the S-SDP
method generates a search direction ∆xk by solving the quadratic subproblem:
minimizesubject to∇f (xk )T ∆x + 12 ∆xT Gk ∆x, ∆x ∈ Rng(xk ) + ∇g(xk )T ∆x = 0,∑X̂k (∆x) ≡ X(xk ) + ni=1 ∆xi Ai (xk ) ≽ 0,

(4.1)

where Gk is the Hessian matrix ∇2x L(xk , yk , Zk ) of the Lagrangian function or its approximation. Letting (∆xk , yk+1 , Zk+1 ) be a KKT point for minimization problem (4.1), we havethe following:∇f (xk ) + Gk ∆xk − ∇g(xk )yk+1 − A∗ (xk )Zk+1g(xk ) + ∇g(xk )T ∆xk
X̂k (∆xk )
Zk+1
X̂k (∆xk )Zk+1

=
=
≽
≽
=

0,
0,
0,
0,
0,

where yk+1 and Zk+1 are multipliers corresponding to the equality and positive semideﬁniteness constraints, respectively. We note that when ∆xk = 0, the above conditions become theKKT conditions (2.1) and (2.2) of problem (1.1). As noted in [23], if the matrix Gk is positive deﬁnite, then problem (4.1) reduces to a linear SDP, because the convex quadratic termin the objective function can be written as a semideﬁniteness constraint or a second-ordercone constraint.
Fares, Noll and Apkarian [18] proposed a local S-SDP method. Later Freund, Jarre andVogelbusch [23] also analyzed the local rate of convergence of the S-SDP method from adiﬀerent viewpoint. The global convergence of the S-SDP method was shown by Correa andRamı́rez [15] under the line search strategy. As another global convergent method, Gómez
and Ramı́rez [25] applied a ﬁlter method to the S-SDP method. These methods will be
brieﬂy reviewed in this section. Furthermore, other related methods were studied by Jarre[33], Kanzow, Nagel, Kato and Fukushima [36], Diehl, Jarre and Vogelbusch [16], Kanno
and Takewaki [34], Li and Sun [46] and Zhu and Zhu [78], for example.
4.1. Local convergence properties of S-SDP methods
Fares, Noll and Apkarian [18] originally dealt with the following minimization problem of alinear objective function subject to linear matrix inequality (LMI) constraints and nonlinear
matrix equality constraints, which arises from robust control design:
minimizesubject to

cT x,
x ∈ Rn ,L(x) ≽ 0, B(x) = 0,

(4.2)

where c is a given vector, L(x) : Rn → Sp is an aﬃne symmetric matrix function, and B(x)is a nonlinear matrix valued function. The matrix function B(x) is bilinear with respect tox in many cases. Fares et. al. [18] solved the problem via the S-SDP method. The following
local convergence theorem was shown for solving problem (1.1) (See Theorem 11 of [15]).

c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝40

H. Yamashita & H. Yabe

Theorem 4.1. Let (x∗ , y ∗ , Z ∗ ) be a point satisfying the KKT conditions of problem (1.1).
Assume that (x∗ , y ∗ , Z ∗ ) satisﬁes
v T ∇2x L(x∗ , y ∗ , Z ∗ )v > 0(4.3)

for all nonzero vectors v ∈ C(x∗ ), where C(x∗ ) denotes the critical cone of (1.1) at x∗ . Suppose that the matrix ((A1 (x∗ ), . . . , An (x∗ ))T , ∇g(x∗ ))T has full rank√and that Gk approaches∇2x L(x∗ , y ∗ , Z ∗ ). Then there exists δ > 0 such that, if ∥x0 −x∗ ∥ < δ, ∥y0 −y ∗ ∥2 +∥Z0 −Z ∗ ∥2F< δ and ∥Gk − ∇2x L(x∗ , y ∗ , Z ∗ )∥ < δ for all k, the sequence (xk , yk , Zk ) generated by the SSDP method with xk+1 = xk +∆xk is well deﬁned and converges superlinearly to (x∗ , y ∗ , Z ∗ ).
Furthermore, the convergence is quadratic if the following holds√
Gk − ∇2x L(x∗ , y ∗ , Z ∗ ) = O(∥xk − x∗ ∥ + ∥yk − y ∗ ∥2 + ∥Zk − Z ∗ ∥2F ).
We note that assumption (4.3) is stronger than (2.3).
Freund, Jarre and Vogelbusch [23] chose a motivation that contrasts the S-SDP methodwith primal-dual interior-point methods and analyzed the local and quadratic convergencebased on a sensitivity result for nonlinear SDP problems. They dealt with the following
nonlinear SDP problem:
minimizesubject to

cT x,
x ∈ Rn ,X(x) ≽ 0(4.4)

where c ∈ Rn is a given vector. In addition, they extended this result to a more generalclass of nonlinear SDP problems with equality and inequality constraints.
4.2. Globally convergent algorithms
In order to establish the global convergence of the S-SDP method, Correa and Ramı́rez [15]introduced the Han penalty function as a merit function of the line search procedure. Fora penalty parameter σ > 0, the merit function is deﬁned byθσ (x) = f (x) − σ(min{0, λmin (X(x))} + ∥g(x)∥).

(4.5)

Then they show that the directional derivative θσ′ (xk ; ∆xk ) of the function θσ (x) along asearch direction ∆xk at xk satisﬁes the following estimate.
Lemma 4.1. Assume that the functions f , g and X are continuously diﬀerentiable and that
their derivatives are Lipschitz continuous in a neighborhood of xk . Using the penalty function
θσ (x) in (4.5), if the point (∆xk , yk+1 , Zk+1 ) satisﬁes the KKT conditions of subproblem(4.1), then the following holds
θσ′ (xk ; ∆xk ) ≤ ∇f (xk )T ∆xk + σ(min{0, λmin (X(xk ))} − ∥g(xk )∥)Tg(xk )= −∆xTk Gk ∆xk − tr(Zk+1 X(xk )) − yk+1
+σ(min{0, λmin (X(xk ))} − ∥g(xk )∥).

(4.6)

Furthermore, if σ ≥ max{tr(Zk+1 ), ∥yk+1 ∥}, we haveθσ′ (xk ; ∆xk ) ≤ −∆xTk Gk ∆xk .
This result shows that the search direction becomes a descent direction for θσ (x), i.e.

θσ′ (xk ; ∆xk ) < 0, when Gk is positive deﬁnite and σ is suﬃciently large.

The algorithm proposed by Correa and Ramı́rez [15] is described for a given current
point xk and a positive deﬁnite matrix Gk as follows:c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

41

Algorithm 4.1.
Step 0. Give x1 ∈ Rn and σ0 > 0. Set k = 1.
Step 1. Compute a point (∆xk , yk+1 , Zk+1 ) satisfying the KKT conditions of subproblem(4.1).
Step 2. Compute σk in a way that the following properties are satisﬁed:(a) σk ≥ max{tr(Zk+1 ), ∥yk+1 ∥} + σ̄.
(b) For all k ≥ k1 , if σk−1 ≥ max{tr(Zk+1 ), ∥yk+1 ∥} + σ̄, then σk = σk−1 .
where a positive integer k1 and σ̄ > 0 are ﬁxed parameters.
(Remark: The parameter σk is updated by σk = max{1.5σk−1 , max{tr(Zk+1 ), ∥yk+1 ∥} +σ̄}, for example.)
Step 3. The step size αk is computed by using the Armijo rule, that is, αk is set to αk = β lk
for the smallest nonnegative integer lk that satisﬁes the suﬃcient decrease condition˜ k,θσk (xk + β lk ∆xk ) ≤ θσk (xk ) + ε0 β lk ∆˜ k is the upper bound of θ′ (xk ; ∆xk )where 0 < ε0 < 1 and 0 < β < 1 are given, and ∆σkgiven in (4.6).
Step 4. Set xk+1 = xk + αk ∆xk .
Step 5. Set k := k + 1 and go to Step 1.
¤
The global convergence property is shown by the following theorem.
Theorem 4.2. Assume that the functions f , g and X are continuously diﬀerentiable andthat their derivatives are Lipschitz continuous. Consider the global algorithm described in
Algorithm 4.1. Suppose that Gk is a positive deﬁnite matrix for all k. If the sequence {Gk }is bounded together with the sequence {G−1k }, then one of the following situations occurs forthe sequence {(xk , yk+1 , Zk+1 )}:1. The sequences {σk } and {(yk+1 , Zk+1 )} are unbounded.
2. There is an index k2 such that σk is constant for all k ≥ k2 . In this case, one of thefollowing situations occurs:(a) θσk (xk ) → −∞,or(b) ∇x L(xk , yk+1 , Zk+1 ) → 0, g(xk ) → 0, min{0, λmin (X(xk ))} → 0,
and tr(Zk+1 X(xk )) → 0.
As another globally convergent method, Gómez and Ramı́rez [25] proposed a ﬁlter
method for solving problem (1.1). The ﬁlter methods were ﬁrst introduced by Fletcher
and Leyﬀer [19] for nonlinear optimization problems. In general, ﬁlter methods deal simultaneously with the optimality and feasibility. In the methods, the trial points are accepted
when they improve either the objective function or a constraint violation measure. These
criteria are compared to previous iterates collected in a ﬁlter. Global convergence results ofseveral kinds of ﬁlter methods have been studied for nonlinear optimization problems (see[20], for example). Given the current iterate xk and matrix Gk , Gómez and Ramı́rez deﬁned
the following trust-region local semideﬁnite approximation of problem (1.1):
minimizesubject to∇f (xk )T ∆x + 21 ∆xT Gk ∆x, ∆x ∈ RnT
g(xk ) + ∇g(x∑n k ) ∆x = 0,
X(xk ) + i=1 ∆xi Ai (xk ) ≽ 0,∥∆x∥∞ ≤ ∆k ,c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝

(4.7)

42

H. Yamashita & H. Yabe

where ∆k is a trust-region radius. In order to investigate the feasibility, let us deﬁne theconstraint violation measure by
θ(x) = ∥g(x)∥ − min{0, λmin (X(x))}.

(4.8)

A ﬁlter, denoted by F, is a ﬁnite collection of two dimensional vectors. In each vectorof the ﬁlter, the two components refer to the values of the feasibility function θ and theobjective function f , respectively. The new candidate (θ̄, f¯) is called acceptable to the ﬁlter
F = {(θj , fj )}N
j=1 , if, for each j = 1, . . . , N , at least one of the following two conditions is
satisﬁed:θ̄ ≤ βθj ,¯
f + γ θ̄ ≤ fj ,where β ∈ (0, 1) and γ ∈ (0, β) are ﬁxed parameters. If the point (θ̄, f¯) is acceptable for thecurrent ﬁlter F, the new ﬁlter Add((θ̄, f¯), F) is deﬁned as
Add((θ̄, f¯), F) = (F ∪ {(θ̄, f¯)})\{(θj , fj ) ∈ F | θ̄ ≤ θj , f¯ ≤ fj }.
The following is a ﬁlter-type algorithm of the S-SDP method.
Algorithm 4.2.
Step 0. Set the parameters β ∈ (0, 1), γ ∈ (0, β), ∆max > ∆min > 0, σ ∈ (0, 1), u > 0 and
k = 1. Deﬁne F 0 = {(u, −∞)}.
Step 1. (Restoration phase) Find a point xk and a corresponding trust-region radius ∆max
≥ ∆k ≥ ∆min such that(a) (θ(xk ), f (xk )) is acceptable to F k−1 .
(b) Subproblem (4.7) is feasible.
Step 2. (Solving subproblem (4.7))
Step 2.1 Solve subproblem (4.7) to obtain a step ∆x.
If subproblem (4.7) is not feasible, then set F k = Add((θ(xk ), f (xk )), F k−1 ), k :=k + 1 and go to Step 1.
If ∆x = 0 is obtained, then xk is a stationary point of problem (1.1) and stop.
Step 2.2 If (θ(xk +∆x), f (xk +∆x)) is not acceptable to the ﬁlter F k−1 ∪{(θ(xk ),f (xk ))},
then set ∆ := ∆/2 and go to Step 2.1.
Step 2.3 If the following conditions are fulﬁlled
1∇f (xk )T ∆x + ∆xT Gk ∆x < 0,
( 2
)
1 T
T
f (xk ) + σ ∇f (xk ) ∆x + ∆x Gk ∆x < f (xk + ∆x),2
then set ∆ := ∆/2 and go to Step 2.1
Step 2.4 Go to Step 3.
Step 3. If ∇f (xk )T ∆x + 12 ∆xT Gk ∆x ≥ 0, then F k = Add((θ(xk ), f (xk )), F k−1 ), otherwise
set F k = F k−1 .
Step 4. Deﬁne ∆k = ∆, ∆xk = ∆x. Set xk+1 = xk + ∆xk , k := k + 1. Finally, reset ∆ suchthat ∆max ≥ ∆ ≥ ∆min is satisﬁed and go to Step 2.
¤
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

43

In the above algorithm, Step 1 is called restoration phase. The restoration phase is executedat the ﬁrst iteration and when subproblem (4.7) is not feasible in Step 2.1. The loop between
Steps 2 and 4 is a main loop, and the loop between Steps 2.1 and 2.3 solves subproblem(4.7) many times for various trust region radii if necessary.
Now we give the global convergence theorem as follows.
Theorem 4.3. Suppose that the points sampled by Algorithm 4.2 lie in a nonempty compact
set Ω ⊂ Rn and that the MFCQ condition is satisﬁed at each feasible point of problem (1.1)
lying in the set Ω. Assume that there exists a positive constant M such that ∥Gk ∥F ≤ M for
all k. Consider the sequence {xk } generated by Algorithm 4.2. Then one of the following
situations occurs:
1. The restoration phase (Step 1) fails to ﬁnd a point xk satisfying (a) and (b).
2. A stationary point of problem (1.1) is found, that is, ∆x = 0 solves subproblem (4.7) forsome iteration k.
3. There exists an accumulation point of {xk } that is a stationary point of problem (1.1).
It is known that the trust-region method [14] is a robust strategy to obtain the global
convergence property for general nonlinear optimization problems. Kanzow, Nagel, Katoand Fukushima [36] incorporated the idea of the trust-region method into the successive linearization method for nonlinear SDP problems, and showed its global convergence property.
Speciﬁcally, they considered the following nonlinear semideﬁnite programming problem:
minimizesubject tof (X),X ∈ Sp ,h(X) ≤ 0, X ≽ 0(4.9)

where X is a variable, and the functions f : Sp → R and h : Sp → Rm are continuously
diﬀerentiable functions. By deﬁning the exact l1 penalty functionpρ (X) = f (X) + ρm∑

max{0, hi (X)}i=1

with a penalty parameter ρ > 0, Kanzow et al. dealt with the penalized problem:
minimize pρ (X) subject to X ≽ 0.

(4.10)

Letting hi denote the ith component of h, and letting Df (X) and Dhi (X) be the Fréchetderivatives of f and hi , respectively, at X, they deﬁned the ﬁrst-order approximation ofpρ (X + ∆X) byΦρ (X, ∆X) = f (X) + ⟨Df (X), ∆X⟩ + ρm∑

max{0, hi (X) + ⟨Dhi (X), ∆X⟩}.

i=1

Based on the above, they calculated a search direction ∆Xk ∈ Sp by solving the following
subproblem:
minimize1
ck ⟨∆X, ∆X⟩ + Φρk (Xk , ∆X) subject to Xk + ∆X ≽ 0,
2

(4.11)

where ck and ρk are suitably chosen parameters. Since the quadratic term of the objectivefunction in (4.11) means that the trust-region idea is implicitly used, this algorithm may bec by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝44

H. Yamashita & H. Yabe

regarded as a successive linearization method for problem (4.9) that employs a trust-regiontype globalization technique. We can also consider that the parameter 1/ck plays the roleof the trust-region radius. Furthermore, the quadratic term guarantees the strict convexityof the subproblem. In their algorithm, they used criteria based on the ratiork =pρk (Xk ) − pρk (Xk + ∆Xk ),pρk (Xk ) − Φρk (Xk , ∆Xk )which is the quotient of the actual and the predicted reductions for the function value ofthe penalty function. Then the point Xk and the parameter ck are updated by the same
method as the trust-region method. In a similar way to this method, Kanno and Takewaki[34] proposed a method for maximum robustness design of structures. Furthermore, Li andSun [46] applied the ﬁlter method to the successive linearization method.
In the algorithms given in this section, the matrix Gk approximates the Hessian matrix
2∇x L(xk , yk , Zk ) of the Lagrangian function by using the quasi-Newton updating formula, if
necessary. In this case, Gk+1 must satisfy the following secant condition
Gk+1 sk = qk ,where sk = xk+1 − xk andqk = ∇x L(xk+1 , yk+1 , Zk+1 ) − ∇x L(xk , yk+1 , Zk+1 )= (∇f (xk+1 ) − ∇g(xk+1 )yk+1 − A∗ (xk+1 )Zk+1 )−(∇f (xk ) − ∇g(xk )yk+1 − A∗ (xk )Zk+1 )= ∇f (xk+1 ) − ∇f (xk ) − (∇g(xk+1 ) − ∇g(xk ))yk+1 − (A∗ (xk+1 ) − A∗ (xk ))Zk+1 .
In order to preserve the positive deﬁniteness of the matrix Gk , we can use the modiﬁed
BFGS update proposed by Powell, which is given by the form
Gk+1 = Gk −Gk sk sTk Gk q̂k q̂kT
+ T ,sTk Gk sk
sk q̂kwhereq̂k = ψk qk + (1 − ψk )Gk sk ,

if sTk qk ≥ 0.2sTk Gk sk
 1T
0.8sk Gk skψk =otherwise.
 T
sk (Gk sk − qk )5. Primal-Dual Interior Point Method
As in ordinary NLP problems and linear SDP problems, the interior point methods can beone of most useful numerical methods for solving nonlinear SDP problems. In this section,
we introduce the primal-dual interior point methods proposed by Yamashita, Yabe andHarada [76], and Yamashita and Yabe [75]. The proposed algorithm is implemented in thesoftware NUOPT (the current name is Numerical Optimizer), and its numerical performanceis described in [76].
We note Jarre [32] gave a primal predictor-corrector type interior point method forsolving nonlinear SDP problem, and Leibfritz and Mostafa [45] gave a primal interior point
method for solving a special class of nonlinear SDP problem. Burer, Monteiro and Zhang [12,c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

4513] converted problem (1.1) to an equivalent NLP problem, and solved the resultant problem
by the interior point method. In [77], Yang and Yu proposed a primal-dual predictorcorrector method based on the homotopy method which is not an interior point type.
A point w = (x, y, Z) satisfying X(x) ≻ 0 and Z ≻ 0 is called an interior point, andinterior point methods generate such points. To construct an interior point algorithm, apositive barrier parameter µ is introduced as in the interior point methods for ordinary
NLP, and the complementarity condition X(x)Z = 0 is replaced by X(x)Z = µI, where
I denotes the identity matrix. The primal-dual interior point methods try to ﬁnd a pointthat satisﬁes the barrier KKT (BKKT) conditions:  ∇x L(w)
0



g(x)
0 
r(w, µ) ≡=X(x)Z − µI
0

(5.1)

and the positivity conditions:X(x) ≻ 0,

Z ≻ 0,

for a given µ. By using the symmetrization (1.2), we also deﬁne the symmetrized residualrS (w, µ) by∇x L(w)
.
g(x)rS (w, µ) = (5.2)
X(x) ◦ Z − µI
This will be used for constructing Newton iteration later, and we denote rS (w, 0) by r0S (w).
We deﬁne the norms ∥r(w, µ)∥ and ∥rS (w, µ)∥ by√°()°° ∇x L(w) °2° + ∥X(x)Z − µI∥2°∥r(w, µ)∥ = °
F
°
g(x)
and

√°(
)°
° ∇x L(w) °2°° + ∥X(x) ◦ Z − µI∥2 ,∥rS (w, µ)∥ = °
F
°
g(x)

respectively. Note that ∥rS (w, µ)∥ ≤ ∥r(w, µ)∥ is satisﬁed, because ∥X(x) ◦ Z − µI∥F ≤∥X(x)Z − µI∥F .
5.1. Algorithm for ﬁnding a KKT point
We ﬁrst describe a procedure for ﬁnding a KKT point by using the BKKT conditions. In
this subsection, the subscript k denotes an iteration count of the outer iterations.
Algorithm 5.1.
Step 0. (Initialize) Set ε > 0, Mc > 0 and k = 0. Let a positive sequence {µk } , µk ↓ 0 be
given.
Step 1. (Approximate BKKT point) Find an interior point wk+1 that satisﬁes∥r(wk+1 , µk )∥ ≤ Mc µk .
Step 2. (Terminate) If ∥r0 (wk+1 )∥ ≤ ε, then stop.
Step 3. (Update) Set k := k + 1 and go to Step 1.

c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝

(5.3)

¤

46

H. Yamashita & H. Yabe

The condition (5.3) is called the approximate BKKT condition, and a point that satisﬁes
this condition the approximate BKKT point. The barrier parameter sequence {µk } in theabove needs not be determined beforehand, and the value of each µk may be set adaptively
as the iteration proceeds. The concrete procedure in Step 1 of Algorithm 5.1 will be given
as Algorithm 5.2 in Section 5.2.
The following theorem shows the convergence property of Algorithm 5.1 under the MFCQ
condition.
Theorem 5.1. Assume that the functions f , g and X are continuously diﬀerentiable. Let{wk } be an inﬁnite sequence generated by Algorithm 5.1. Suppose that the sequence {xk } isbounded and that the MFCQ condition is satisﬁed at any accumulation point of the sequence{xk }. Then the sequences {yk } and {Zk } are bounded, and any accumulation point of {wk }satisﬁes the KKT conditions (2.1) and (2.2).
5.2. Newton method for ﬁnding a barrier KKT point and primal-dual merit
function
In [76], a globally convergent algorithm to ﬁnd a BKKT point for a given ﬁxed barrierparameter µ > 0 was proposed. The algorithm given below is used as an inner iteration of
Algorithm 5.1. This subsection describes a Newton-like method to the system of equations(5.4). We denote the Newton directions ∑for the primal and dual variables by ∆x ∈ Rn , ∆y ∈mp
R and ∆Z ∈ S . We deﬁne ∆X = ni=1 ∆xi Ai (x), and note that ∆X ∈ Sp . In orderto have a Jacobian operator from Rn × Rm × Sp to Rn × Rm × Sp (a square matrix whenrepresented in matrix form), we employ the symmetrized residual representation of the
BKKT conditions rS (w, µ) = 0.
As in the case of linear SDP problems, a scaling of the primal-dual pair (X(x), Z) is
used in order to make the resulting Newton equation solvable. In what follows, we denote
X(x) simply by X when it is not confusing. Throughout this subsection, we assume thatX ≻ 0 and Z ≻ 0 hold. We introduce a nonsingular matrix T ∈ Rp×p and scale X and Z
bye = T XT T and Ze = T −T ZT −1 ,Xe ◦ Ze = µI instead
respectively. Using the scaling matrix T , we try to solve the equation X
of XZ = µI, and consider the scaled symmetrized residual:  ∇x L(w)0 =  0 ,g(x)r̃S (w, µ) ≡ (5.4)
e ◦ Ze − µI0X
to form the Newton directions.
Scalings of ∆X and ∆Z are similarly done by
e = T ∆XT T∆Xand ∆Ze = T −T ∆ZT −1 .

Newton equations for (5.4) are given byG∆x − ∇g(x)∆y − A∗ (x)∆Z = −∇x L(x, y, Z)∇g(x)T ∆x = −g(x)1e Ze + Z∆
e X
e + X∆e Ze + ∆ZeX)e = µI − 1 (Xe Ze + ZeX),
e
(∆X
2
2

(5.5)
(5.6)
(5.7)

where G denotes the Hessian matrix of the Lagrangian function L(w) or its approximation.
By solving these equations, we obtain the Newton direction ∆w = (∆x, ∆y, ∆Z) ∈ Rn ×
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

47

Rm × Sp . Writing the last part of Newton equations as
e + (Xe ⊙ I)∆Ze = µ(Xe ⊙ I)Xe −1 − (Xe ⊙ I)Z,e(Ze ⊙ I)∆Xe
we can formally solve this equation for ∆Z.
The following theorem [76] gives the desired form of the Newton directions.
e ⊙ I is invertible. Then the direction ∆Ze ∈ Sp
Theorem 5.2. Suppose that the operator X
is given by the forme −1 − Ze − (Xe ⊙ I)−1 (Ze ⊙ I)∆X,e∆Ze = µX
(5.8)
or equivalentlye ⊙ I)−1 (Ze ⊙ I)(T ⊙ T )∆X.
∆Z = µX −1 − Z − (T T ⊙ T T )(X
Furthermore, the directions (∆x, ∆y) ∈ Rn × Rm satisfy
(
)(
)
(
)
G + H −∇g(x)∆x∇f (x) − ∇g(x)y − µA∗ (x)X −1=−,−∇g(x)T0∆y−g(x)where the elements of the matrix H ∈ Rn×n are represented by the form
⟨
⟩
−1 e
e
e
e
Hij = Ai (x), (X ⊙ I) (Z ⊙ I)Aj (x)(5.9)(5.10)(5.11)ei (x) = T Ai (x)T T .
with AIn addition, if the matrix G + H is positive deﬁnite and the matrix ∇g(x) is of full rank,
then the Newton equations (5.5) – (5.7) give a unique search direction ∆w = (∆x, ∆y, ∆Z)∈ Rn × Rm × Sp .
e and Ze commute, i.e.,
In the following, we choose a nonsingular matrix T so that X
e
e
e
e
e
e
X Z = Z X. In this case, the matrices X and Z share the same eigensystem, and XZ =e ◦ Ze = µI. In the algorithm described here, this commutativity is necessaryµI ⇐⇒ X
to prove a descent property of the Newton direction with respect to the merit function
described later.
Well known examples of the scaling matrix T that satisfy such a condition are given
below.
(i) HRVW/KSH/M directione = I and Ze = X 1/2 ZX 1/2 , which corresponds to theIf we set T = X −1/2 , then we have X
HRVW/KSH/M direction for linear SDP problems [26, 42, 48]. In this case, the matrices Hand ∆Z can be represented by the form:()
Hij = tr Ai (x)X −1 Aj (x)Z ,1∆Z = µX −1 − Z − (X −1 ∆XZ + Z∆XX −1 ).
2(ii) NT directione = W −1/2 XW −1/2
If we set T = W −1/2 with W = X 1/2 (X 1/2 ZX 1/2 )−1/2 X 1/2 , then we have X
1/2
1/2
e which corresponds to the NT direction for linear SDP problems [50, 51].
= W ZW= Z,
In this case, the matrices H and ∆Z can be represented by the form:{}
Hij = tr Ai (x)W −1 Aj (x)W −1 ,∆Z = µX −1 − Z − W −1 ∆XW −1 .
e ⊙ I)−1 in (5.8) tractable by letting (Xe ⊙ I)−1 = I (HRVW/KSH
Both examples make (Xe ⊙ I)−1 (Ze ⊙ I) = I (NT).
/M) or (XUnder the above scalings, if the matrices Ai (x) (i = 1, . . . , n) are linearly independent,
the matrix H is symmetric positive deﬁnite.
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝48

H. Yamashita & H. Yabe

5.3. Global convergence of primal-dual interior point methods
In [76], the line search strategy was used, and the following merit function in the primal-dual
space was proposed:
F (x, Z) = FBP (x) + νFP D (x, Z).
(5.12)
Here FBP (x) and FP D (x, Z) are the primal barrier penalty function and the primal-dual
barrier function, respectively, and they are given by
FBP (x) = f (x) − µ log(detX) + ρ∥g(x)∥1 ,FP D (x, Z) = ⟨X, Z⟩ − µ log(detXdetZ),(5.13)(5.14)

where ν and ρ are positive parameters. Though the functions FBP (x) and FP D (x, Z) dependon the parameters ν, ρ and ⟨µ, we ⟩use the notation F (x, Z) for simplicity. It follows from thee Ze = T XZT −1 that X,
e Ze = ⟨X, Z⟩ and the value of FP D (x, Z) is invariant underfact X
the choice of T .
We deﬁne the ﬁrst order approximation Fl of the merit function by
Fl (x, Z; ∆x, ∆Z) = F (x, Z) + ∆Fl (x, Z; ∆x, ∆Z).
Here ∆Fl (x, Z; ∆x, ∆Z) corresponds to the directional derivative and it is deﬁned by theform∆Fl (x, Z; ∆x, ∆Z) = ∆FBP l (x; ∆x) + ν∆FP Dl (x, Z; ∆x, ∆Z),
where
∆FBP l (x; ∆x) = ∇f (x)T ∆x − µtr(X −1 ∆X)()+ρ ∥g(x) + ∇g(x)T ∆x∥1 − ∥g(x)∥1 ,∆FP Dl (x, Z; ∆x, ∆Z) = tr(∆XZ + X∆Z − µX −1 ∆X − µZ −1 ∆Z).

(5.15)

With these deﬁnitions, it is possible to prove that the Newton direction ∆w satisﬁes∆Fl (x, Z; ∆x, ∆Z) ≤ −∆xT (G + H) ∆x − (ρ − ∥y + ∆y∥∞ )∥g(x)∥1 .
By using this inequality, the following theorem is obtained, which shows that the Newton
direction ∆w becomes a descent search direction for the proposed primal-dual merit function(5.12).
Theorem 5.3. Assume that ∆w solves (5.5) – (5.7) and that the matrix G+H is symmetricpositive deﬁnite. Suppose that the penalty parameter ρ satisﬁes ρ > ∥y + ∆y∥∞ . Then thefollowing hold:(i) The direction ∆w becomes a descent search direction for the primal-dual merit function
F (x, Z), i.e. ∆Fl (x, Z; ∆x, ∆Z) ≤ 0.
(ii) If ∆x ̸= 0, then ∆Fl (x, Z; ∆x, ∆Z) < 0.
(iii) ∆Fl (x, Z; ∆x, ∆Z) = 0 holds if and only if (x, y + ∆y, Z) is a BKKT point.
In order to construct a globally convergent algorithm to a BKKT point for a ﬁxed µ > 0,
the iterations take the form
xk+1 = xk + αk ∆xk ,Zk+1 = Zk + αk ∆Zk

and yk+1 = yk + ∆yk

where αk is a step size determined by the line search procedure described below. Throughout
this subsection, the index k denotes the inner iteration count for a given µ > 0. We also
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

49

denote X(xk ) by Xk for simplicity. We note that Xk ≻ 0 and Zk ≻ 0 are maintained for all
k in the following.
Since the main iteration is to decrease the value of the merit function (5.12), the step
size is determined by the suﬃcient-decrease rule of the merit function. At the current point(xk , Zk ), the initial step size is calculated by{− λ (Xγ−1 ∆X ) if X(x) is linearminkkᾱxk =(5.16)
1otherwiseandᾱzk = −γ,λmin (Zk−1 ∆Zk )(5.17)

where γ ∈ (0, 1) is a constant. If the minimum eigenvalue in either expression (5.16) or(5.17) is positive, the corresponding step size is set to 1. By trying the values ℓk = 0, 1, . . .
successively, we ﬁnd a step to the next iterate given by
αk = ᾱk β lk ,ᾱk = min {ᾱxk , ᾱzk , 1} ,where β ∈ (0, 1) is a constant, and lk is the smallest nonnegative integer such that thesuﬃcient-decrease condition
F (xk + ᾱk β lk ∆xk , Zk + ᾱk β lk ∆Zk ) ≤ F (xk , Zk ) + ε0 ᾱk β lk ∆Fl (xk , Zk ; ∆xk , ∆Zk )(5.18)and the positive deﬁniteness condition
X(xk + ᾱk β lk ∆xk ) ≻ 0(5.19)

hold, where ε0 ∈ (0, 1) is a constant.
Now we summarize the line search algorithm in the following. Since this algorithm
should be regarded as the inner iteration of Algorithm 5.1 (Step 1 of Algorithm 5.1), ε′given below corresponds to Mc µ and an initial point can be set to the approximate BKKT
point obtained at the previous outer iteration.
Algorithm 5.2.
Step 0. (Initialize) Let w0 ∈ Rn × Rm × Sp (X0 ≻ 0, Z0 ≻ 0), µ > 0, ρ > 0 and ν > 0 be
given. Set ε′ > 0, γ ∈ (0, 1), β ∈ (0, 1) and ε0 ∈ (0, 1). Let k = 0.
Step 1. (Termination) If ∥r(wk , µ)∥ ≤ ε′ , then stop.
Step 2. (Search direction) Calculate the matrix Gk and the scaling matrix Tk . Determine
the direction ∆wk by solving (5.5) – (5.7). If ∆xk = 0, then set wk+1 = (xk , yk +∆yk , Zk +∆Zk ), and stop.
Step 3. (Step size) Find the smallest nonnegative integer lk that satisﬁes the criteria (5.18)
and (5.19), and calculate
αk = ᾱk β lk .
Step 4. (Update) Set
xk+1 = xk + αk ∆xk ,Zk+1 = Zk + αk ∆Zk

and yk+1 = yk + ∆yk .

Step 5. Set k := k + 1 and go to Step 1.

c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝¤50

H. Yamashita & H. Yabe

If the matrix Gk approximates the Hessian matrix ∇2x L(wk ) by the quasi-Newton formula, the method stated in the last part of Section 4 can be used.
If we want to use the Hessian matrix ∇2x L(wk ) itself for the matrix Gk , the LevenbergMarquardt type modiﬁcation of ∇2x L(wk ) may be used to obtain a positive semideﬁnite Gk
for global convergence property shown below. Namely, we compute a parameter σ ≥ 0which gives a positive semideﬁnite matrix ∇2x L(wk ) + σI. This strategy may be useful forsolving large scale problems.
The following global convergence theorem was proved by Yamashita, Yabe and Harada[76].
Theorem 5.4. If ∆xk = 0 in Step 2 of Algorithm 5.2, then (xk , yk + ∆yk , Zk ) is a BKKT
point. Otherwise let an inﬁnite sequence {wk } be generated by Algorithm 5.2. Suppose thefollowings:(i) The functions f , gi , i = 1, . . . , m, and X are twice continuously diﬀerentiable;(ii) The sequence {xk } remains in a compact set Ω of Rn ;(iii) The matrix ∇g(xk ) is of full rank and the matrices A1 (xk ), . . . , An (xk ) are linearly
independent for all xk in Ω;
ek and Zek commute;(iv) The scaling matrix Tk is chosen so that X(v) The matrix Gk is uniformly bounded and positive semideﬁnite;(vi) Both of the sequences {Tk } and {Tk−1 } are bounded;(vii) The penalty parameter ρ is suﬃciently large so that ρ > ∥yk + ∆yk ∥∞ holds for all k.
Then there exists at least one accumulation point of {wk }, and any accumulation pointof the sequence {wk } is a BKKT point.
5.4. The rate of convergence of primal-dual interior point methods
In this section, we survey the local behavior of the primal-dual interior point methods
following the study by Yamashita and Yabe [75]. Consider the scaled Newton equations(5.5)-(5.7) and denote the Newton equations by
JeS (w)∆w = −r̃S (w, µ),(5.20)

where JeS (w) is the linear operator from Rn × Rm × Sp to Rn × Rm × Sp and r̃S (w, µ) is
deﬁned by (5.4). If we choose T = I, we call the above equations the unscaled Newton
equations and denote the operator JeS (w) by JS (w) in this case. This case corresponds to
the AHO direction for linear SDP problems [2].
In order to analyze the local behavior of the methods, we make the following assumptionsin this subsection: The second derivatives of the functions f , gi , i = 1, . . . , m, and X are
Lipschitz continuous at x∗ ; The second order suﬃcient condition (2.5) holds at x∗ ; The strict
complementarity condition holds at x∗ ; The nondegeneracy condition is satisﬁed at x∗ .
We will discuss a local behavior of the unsymmetric residual r0 (w) in (2.1) or r(w, µ) in(5.1). For this purpose, we deﬁne a linear operator J : Rn × Rm × Sp → Rn × Rm × Rp×p
at w by 2
∇x L(w)∆x − ∇g(x)∆y − A∗ (x)∆Z∇g(x)T ∆x
J(w)∆w = ∆XZ + X∆Z
for ∆w = (∆x, ∆y, ∆Z) ∈ Rn × Rm × Sp , which is an estimate of the ﬁrst order change ofr0 (w + ∆w) or r(w + ∆w, µ).
For the unscaled Newton method, the Newton equations (5.20) is
JS (w)∆w = −rS (w, µ).
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝(5.21)

A Survey of Numerical Methods for Nonlinear SDP

51

Under the above assumptions, the regularity of the operators JS (w∗ ) and J(w∗ ) is shown as
follows.
Theorem 5.5. ([75]) The operator JS (w∗ ) is nonsingular, and the operator J(w∗ ) is left
invertible.
With the above properties of JS (w∗ ) and J(w∗ ), an algorithm for the analysis of local
behavior can be deﬁned as follows.
Algorithm 5.3.
Step 0. (Initialize) Set ε > 0 and 0 < τ < 1. Choose w0 = (x0 , y0 , Z0 ) ∈ Rn × Rm × Sp
(X(x0 ) ≻ 0, Z0 ≻ 0). Set k = 0.
Step 1. (Termination) If ∥r0 (wk )∥ ≤ ε, then stop.
Step 2. (Newton step) Choose a barrier parameter µk such thatµk = ξk ∥r0 (wk )∥1+τ

(5.22)

with ξk = Θ(1). Compute the direction ∆wk by solving the Newton equations (5.21).
Set wk+1 = wk + ∆wk .
Step 3. (Update) Set k := k + 1 and go to Step 1.
¤
We note that in the above algorithm, the barrier parameter is decreased faster than∥r0 (wk )∥, and the step sizes are set to 1. By Theorem 5.5, if the iterate wk is suﬃciently close
to w∗ , the Jacobian operator JS (wk ) is nonsingular and its inverse is uniformly bounded.
Thus the Newton equations have a unique solution and the following relations hold∆wk = Θ(∥rS (wk , µk )∥) = O(∥r0S (wk )∥) + O(µk ) = O(∥r0 (wk )∥).

(5.23)

Here the last equality can be obtained by equation (5.22) which shows that the parameterµk is decreased suﬃciently fast. Then it is possible to show that the unit step iteration
gives an interior point and an approximate BKKT point successively. Thus we can prove
the superlinear convergence of Algorithm 5.3 in the following theorem.
Theorem 5.6. Let Mc and τ be given constants satisfying 0 < Mc < 1 and 0 < τ < 1.
Let µ−1 be a suﬃciently small positive number. Assume that an initial interior point w0
is suﬃciently close to w∗ and satisﬁes the approximate BKKT condition ∥r(w0 , µ−1 )∥ ≤Mc µ−1 . Then the sequence {wk } generated by Algorithm 5.3 satisﬁes∥r(wk , µk−1 )∥ ≤ Mc µk−1 ,X(xk ) ≻ 0 and Zk ≻ 0(5.24)

for all k ≥ 0 and converges locally and superlinearly to w∗ .
We note that although no globally convergent algorithm for unscaled Newton method
is proposed yet, this method has favorable local properties as seen in the above theorem.
We also note that, to solve the unscaled Newton equation, we have to represent it by thefollowing matrix-vector form:


 

−∇g(x) −AS (x)T∆x
∇2x L(w)−∇x L(x, y, Z)


=
,
∇g(x)T00∆y−g(x)(Z ⊗S I)AS (x)0(X ⊗S I)
svec(∆Z)
svec(µI − X ◦ Z)
where
AS (x) = [svec(A1 (x)), . . . , svec(An (x))] ∈ Rp(p+1)/2×n .
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝52

H. Yamashita & H. Yabe

Next, we proceed to results on the local behavior of the scaled Newton equations. Specifically we show that local and two-step superlinear convergence properties are proved for the
primal-dual interior point methods which use the HRVW/KSH/M and the NT directions.
In the following, we describe an algorithm which calculates a KKT point by using thescaled Newton method. Roughly speaking, if the current point wk is suﬃciently close to the
BKKT point with suﬃciently small µk , then the Newton equation (5.20) can be written as
JS′ (wk )∆wk = −rS (wk , µk )with an operator JS′ (wk ) which is close to nonsingular JS (wk ), and ∆wk satisﬁes∆wk = O(∥r(wk , µk )∥).
With these properties, the scaled Newton step gives a local and two-step superlinear convergence as the following algorithm and theorem show.
Algorithm 5.4.
Step 0. (Initialize) Set ε > 0 and 0 < τ < 1. Choose w0 = (x0 , y0 , Z0 ) ∈ Rn × Rm × Sp
(X(x0 ) ≻ 0, Z0 ≻ 0). Set k = 0.
Step 1. (Termination) If ∥r0 (wk )∥ ≤ ε, then stop.
Step 2. (Scaled Newton steps)
Step 2.1 Choose µk = ξk ∥r0 (wk )∥1+τ with ξk = Θ(1).
Step 2.2 Calculate the direction ∆wk by solving the scaled Newton equations
JeS (wk )∆wk = −r̃S (wk , µk ) at wk . Set wk+ 1 = wk + ∆wk .
2
Step 2.3 Calculate the direction ∆wk+ 1 by solving the scaled Newton equations
2

JeS (wk+ 1 )∆wk+ 1 = −r̃S (wk+ 1 , µk )22

2

at wk+ 1 . Set wk+1 = wk+ 1 + ∆wk+ 1 .
2

2

2

¤

Step 3. (Update) Set k := k + 1 and go to Step 1.

In the above algorithm, the decreasing speed of the barrier parameter is set to the same as
the unscaled algorithm. But the scaled Newton steps with unit step sizes are repeated twicewith the same barrier parameter value. The following theorem gives a detailed behavior of
the scaled Newton steps near the KKT point.
Theorem 5.7. Let Mc be a positive constant, and let τ and τ ′ be positive constants thatsatisfy2τ
1 > τ ′ > τ and τ ′ >.
1−τ
Let µ−1 be a suﬃciently small positive number and satisfy)1/τ ′(
1≥ µ−1 .
2Mc
Assume that an initial interior point w0 is suﬃciently close to w∗ and satisﬁes the ap′
proximate BKKT condition ∥r(w0 , µ−1 )∥ ≤ Mc µ1+τ
−1 . Then the sequence {wk } generated by
−1/2
−1/2
Algorithm 5.4 with Tk = Xk
or Tk = Wksatisﬁes′∥r(wk , µk−1 )∥ ≤ Mc µ1+τ
k−1 ,X(xk ) ≻ 0 and Zk ≻ 0

for all k ≥ 0 and converges two-step superlinearly to w∗ in the sense that′∥wk + ∆wk + ∆wk+ 1 − w∗ ∥ = O(∥wk − w∗ ∥1+τ ) for all k.
2c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

53

5.5. Primal-dual interior point methods with shifted KKT conditions
Interior point methods usually replace the complementarity condition X(x)Z = 0 by X(x)Z= µI and manage the resulting barrier KKT conditions instead of conditions (2.1) and(2.2). These conditions are intimately related to the necessary condition for optimality of
the following problemmin
X(x)≻0FBP (x, µ) = f (x) + ρ∥g(x)∥1 − µ log(detX(x)), x ∈ Rn .

(5.25)

In connection with (1.1), we consider the following problem instead of (5.25):
min
X(x)≻0FBP 2 (x, µ) = f (x) +1∥g(x)∥2 − µ log(detX(x)), x ∈ Rn .
2µ(5.26)The necessary conditions for optimality of this problem are given by∇FBP 2 (x, µ) = ∇f (x) +1∇g(x)g(x) − µA∗ (x)X −1 (x) = 0,µand X(x) ≻ 0. If we let the variables y and Z satisfy the relations y = −g(x)/µ andZ = µX(x)−1 , respectively, then the above conditions are written as
  

0
∇x L(w)
0 
g(x) + µy(5.27)=
r(w, µ) =0
X(x)Z − µI
andX(x) ≻ 0,

Z ≻ 0.

These conditions are called the shifted barrier KKT (SBKKT) conditions, and the point
which satisﬁes SBKKT conditions are called the SBKKT point. Motivated by the prior workby Forsgren and Gill [22], Yamashita and Yabe [73] proposed a method that uses SBKKT
conditions for ordinary NLP. An extension of the method for solving nonlinear SDP wasdone by Kato, Yabe and Yamashita [37]. We note that similar approach for nonlinear SDP
was also done by Yamakawa and Yamashita [72].
In order to obtain the KKT point, the same prototype algorithm as Algorithm 5.1 isused, where condition (5.3) in Step 1 of Algorithm 5.1 should use the above residual (5.27).
As in the previous section, the Newton method for equations (5.27) is considered to search
the SBKKT point for ﬁxed µ, and we obtain the Newton step ∆w by solving the following
linear system:G∆x − ∇g(x)∆y − A∗ (x)∆Z = −∇x L(x, y, Z)∇g(x)T ∆x + µ∆y = −g(x) − µye Ze + Z∆
e X
e + X∆e Ze + ∆ZeX
e = 2µI − X
e Ze − ZeX.
e
∆X

(5.28)(5.29)(5.30)

Then, the equations (5.28) – (5.30) are rewritten by{G + H +1∇g(x)∇g(x)T }∆x = −∇FBP 2 (x, µ),µ1∆y = − {g(x) + µy + ∇g(x)T ∆x}µc by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝(5.31)(5.32)

54

andH. Yamashita & H. Yabee ⊙ I)−1 (Ze ⊙ I)(T ⊙ T )∆X,
∆Z = µX −1 − Z − (T T ⊙ T T )(X(5.33)

where the elements of the matrix H ∈ Rn×n are given by (5.11).
If the matrix G + H + µ1 ∇g(x)∇g(x)T is nonsingular, then the Newton equations (5.28)– (5.30) have a unique search direction ∆w = (∆x, ∆y, ∆Z) ∈ Rn × Rm × Sp obtained by(5.31) – (5.33). To obtain a unique direction, the interior point method in [76] needs theassumption that ∇g(x) is of full rank, while this method does not need such an assumption.
In contrast to the merit function FBP (x, µ) which contains the non-diﬀerentiable term
and which does not depend on y, the merit function in [37] has the following form:
Fℓ2 (w, µ) = FBP 2 (x, µ) + νFP D2 (w, µ),where FBP 2 (x, µ) is deﬁned by (5.26) and FP D2 (w, µ) is a primal-dual barrier function givenby1
1
1⟨X, Z⟩/p + ∥Z 2 XZ 2 − µI∥2F2
FP D2 (w, µ) = ∥g(x) + µy∥ + log.
2{det(XZ)}1/p
It is possible to prove that the Newton direction is a descent direction with respect tothis merit function. The Newton direction ∆w satisﬁes
1
D(Fℓ2 (w, µ); ∆w) ≤ −∆xT {G + H + ∇g(x)∇g(x)T }∆xµ1
1∥Z 2 XZ 2 − µI∥2F−ν− ν∥g(x) + µy∥2 ,1
1
2⟨X, Z⟩/p + ∥Z 2 XZ 2 − µI∥F
where D(Fℓ2 (w, µ); ∆w) is the directional derivative of the function Fℓ2 (w, µ) along thedirection ∆w. If the matrix G + H + µ1 ∇g(x)∇g(x)T is positive deﬁnite, then the Newton
direction ∆w is a descent direction for the merit function Fℓ2 (w, µ).
The algorithm that ﬁnds an approximate SBKKT point given in [37] is similar to Algorithm 5.2. The only diﬀerence is that ∆Fl (xk , Zk ; ∆xk , ∆Zk ) in Algorithm 5.2 is replaced by the directional derivative D(Fℓ2 (wk , µ); ∆wk ), because the merit function Fℓ2 is
diﬀerentiable.
Algorithm 5.5.
Step 0. (Initialize) Set an initial interior point w0 = (x0 , y0 , Z0 ) ∈ Rn × Rm × Sp (X(x0 ) ≻0, Z0 ≻ 0), the ﬁxed barrier parameter µ > 0 and a parameter Mc > 0. Choose theparameters ν > 0, γ ∈ (0, 1), β ∈ (0, 1) and ε0 ∈ (0, 1). Set k = 0.
Step 1. (Termination) If ∥r(wk , µ)∥ ≤ Mc µ, then stop.
Step 2. (Search direction) Calculate the matrix Gk and the scaling matrix Tk . Determine
a search direction ∆wk by (5.31) – (5.33).
Step 3. (Step size) Calculate an initial step size ᾱk byᾱk = min {ᾱxk , ᾱzk , 1} ,where ᾱxk and ᾱzk are calculated by (5.16) and (5.17). Find the smallest nonnegative
integer ℓk that satisﬁes the conditionskk

Fℓ2 (wk + ᾱk β ℓ ∆wk , µ) ≤ Fℓ2 (wk , µ) + ε0 ᾱk β ℓ D(Fℓ2 (wk , µ); ∆wk )andX(xk + ᾱk β ℓk ∆xk ) ≻ 0.
Set αk = ᾱk β ℓk .
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

55

Step 4. (Update) Set
wk+1 = wk + αk ∆wk .
Step 5. Set k := k + 1 and go to Step 1.
The global convergence of the proposed algorithm is proved in [37]:
Theorem 5.8. Let {wk } be the sequence generated by Algorithm 5.5. Suppose the followings:(i) The functions f , gi , i = 1, . . . , m, and X are twice continuously diﬀerentiable;(ii) The sequence {xk } remains in a compact set Ω of Rn ;(iii) The matrix Gk + Hk + µ1 ∇g(xk )∇g(xk )T is uniformly positive deﬁnite and the matrix
Gk + Hk is uniformly bounded;
ek and Zek commute and both of the sequences(iv) The scaling matrix Tk is chosen so that X−1{Tk } and {Tk } are bounded.
Then there exists at least one accumulation point of {wk }, and any accumulation pointof the sequence {wk } is an SBKKT point.
6. Concluding Remarks
In this paper, current status of numerical methods for solving the nonlinear SDP problems
has been surveyed. It is seen that typical methods for ordinary NLP problems, such as
the augmented Lagrangian method, the SQP method and the primal-dual interior point
method, have their counterparts in this area.
The augmented Lagrangian method, especially the one with the modiﬁed barrier typeaugmentation, has been applied to the wide area of problems using PENNON [38] by
Kočvara and Stingl. The prototype algorithm is rather simple, and the subproblems tobe solved are ordinary nonlinear optimization problems. The globally convergent algorithm
that assures bounded sequences could be somewhat complex though. The primal-dual interior point method proposed by Yamashita, Yabe and Harada [76] has nice global andlocal convergence properties, and implemented in the software Numerical Optimizer (its oldname is NUOPT). The sequential SDP methods seem to have rooms for more research intheory and practice. Also it is not so obvious that the methods that have linear/quadratic
SDP subproblems can be practical or not, because these subproblems should be solved byan interior point method anyway. Notable practical implementation of the sequential SDP
method is not yet known to the current authors.
In order to have a brief comparison of the algorithms described in this paper, we cite
the comparative study of numerical results of a class of problems solved by PENNON and
Numerical Optimizer given in [76]. For detailed explanation about computational aspects ofthe comparison, we refer to the original paper. The following example problem area is theso called static output feedback (SOF) problems from COM P le ib library [44]. The SOF-H2
type problem is deﬁned as follows:
minimize tr(X),subject to Q ≽ 0,T
A(FB1T ≼ 0,( )Q + QA(F ) + B1)
X
C(F )Q≽ 0,
T
QC(F )
Q
where X ∈ Snz ×nz , F ∈ Rnu ×ny and Q ∈ Snx ×nx are variable matrices to be determined.
The matrices A ∈ Rnx ×nx , B ∈ Rnx ×nu , B1 ∈ Rnx ×nw , C ∈ Rny ×nx , C1 ∈ Rnz ×nx , D11 ∈c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝56

H. Yamashita & H. YabeRnz ×nw , D12 ∈ Rnz ×nu and D21 ∈ Rny ×nw are given constant matrices, and form the matrices
A(F ), B(F ), C(F ), D(F ) which appear in the problem deﬁnition as follows:
A(F )
B(F )
C(F )
D(F )

=
=
=
=

A + BF C,
B1 + BF D21 ,C1 + D12 F C,
D11 + D12 F D21 .

Table 8 in [76] shows the required CPU(sec) data for solutions of the SOF-H2 problems
by Numerical Optimizer and PENBMI which is a specialized BMI-version of PENNON.
The number of variables of the test problem ranges from 7 to 160. The required CPU by
Numerical Optimizer varies from 0.03 sec to 12.19 sec, and that by PENBMI varies from
0.18 sec to 223 sec. As the table shows, the performances of these two codes are comparable,
and it is diﬃcult to draw any ﬁrm conclusion about relative merits of these two algorithms
from the table. However we should point out that interior point methods need an interior
point for starting the calculation, and it is not so easy to ﬁnd an interior point for this typeof problems.
In summary, the area of nonlinear SDP problems seems to be still in the early stage.
We believe there appear new proposals of algorithms and their variations in the future as
in ordinary NLP world.
Acknowledgement. The authors are grateful to the referees for their valuable comments that greatly helped to improve the paper. The second author was supported in partby the Grant-in-Aid for Scientiﬁc Research (C) 25330030 of Japan Society for the Promotionof Science.
References[1] F. Alizadeh and D. Goldfarb: Second-order cone programming. Mathematical Programming, 95 (2003), 3–51.
[2] F. Alizadeh, J.A. Haeberly, and M.L. Overton: Primal-dual interior-point methods forsemideﬁnite programming – Convergence rates, stability and numerical results. SIAM
Journal on Optimization, 8 (1998), 746–768.
[3] M.F. Anjos and J.B. Lasserre (eds.): Handbook on Semideﬁnite, Conic and Polynomial
Optimization (Springer, New York, 2012).
[4] A. Ben-Tal and A. Nemirovski: Lectures on Modern Convex Optimization, MPS-SIAM
Series on Optimization (SIAM, Philadelphia, 2001).
[5] D.P. Bertsekas: Nonlinear Programming – Second Edition (Athena Scientiﬁc, Belmont,
MA, 1999).
[6] P.T. Boggs and J.W. Tolle: Sequential quadratic programming. Acta Numerica, 4(1995), 1–51.
[7] J.F. Bonnans, R. Cominetti, and A. Shapiro: Sensitivity analysis of optimization problems under second order regular constraints. Mathematics of Operations Research, 23(1998), 806–831.
[8] J.F. Bonnans, R. Cominetti, and A. Shapiro: Second order optimality conditions basedon parabolic second order tangent sets. SIAM Journal on Optimization, 9 (1999), 466–
492.
[9] J.F. Bonnans and A. Shapiro: Perturbation Analysis of Optimization Problems(Springer Verlag, New York, 2000).
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

57[10] S. Boyd, L.E. Ghaoui, E. Feron, and V. Balakrishnan: Linear Matrix Inequalities
in System and Control Theory, SIAM Studies in Applied Mathematics, 15 (SIAM,
Philadelphia, 1994).
[11] S. Boyd and L. Vandenberghe: Convex Optimization (Cambridge University Press,
2004).
[12] S. Burer, R.D.C. Monteiro, and Y. Zhang: Solving a class of semideﬁnite programs via
nonlinear programming. Mathematical Programming, 93 (2002), 97–122.
[13] S. Burer, R.D.C. Monteiro, and Y. Zhang: Interior-point algorithms for semideﬁnite
programming based on a nonlinear formulation. Computational Optimization and Applications, 22 (2002), 49–79.
[14] A.R. Conn, N.I.M. Gould, and P.L. Toint: Trust-Region Methods, MPS-SIAM Series
on Optimization (SIAM, Philadelphia, 2000).
[15] R. Correa and C.H. Ramı́rez: A global algorithm for nonlinear semideﬁnite programming. SIAM Journal on Optimization, 15 (2004), 303–318.
[16] M. Diehl, F. Jarre, and C.H. Vogelbusch: Loss of superlinear convergence for an SQPtype method with conic constraints. SIAM Journal on Optimization, 16 (2006), 1201–
1210.
[17] B. Fares, P. Apkarian, and D. Noll: An augmented Lagrangian method for a class of
LMI-constrained problems in robust control theory. International Journal of Control,
74 (2001), 348–360.
[18] B. Fares, D. Noll, and P. Apkarian: Robust control via sequential semideﬁnite programming. SIAM Journal on Control and Optimization, 40 (2002), 1791–1820.
[19] R. Fletcher and S. Leyﬀer: Nonlinear programming without a penalty function. Mathematical Programming, 91 (2002), 239–269.
[20] R. Fletcher, S. Leyﬀer, and P. Toint: A brief history of ﬁlter methods. Preprint
ANL/MCS-P1372-0906, September 26, 2006 (revised October 9, 2006), Argonne National Laboratory, Illinois, USA.
[21] A. Forsgren: Optimality conditions for nonconvex semideﬁnite programming. Mathematical Programming, 88 (2000), 105–128.
[22] A. Forsgren and Ph.E. Gill: Primal-dual interior methods for nonconvex nonlinear
programming. SIAM Journal on Optimization, 8 (1998), 1132–1152.
[23] R.W. Freund, F. Jarre, and C.H. Vogelbusch: Nonlinear semideﬁnite programming:
sensitivity, convergence, and an application in passive reduced-order modeling. Mathematical Programming, 109 (2007), 581–611.
[24] E.H. Fukuda, P.J.S. Silva, and M. Fukushima: Diﬀerentiable exact penalty functions
for nonlinear second-order cone programs. SIAM Journal on Optimization, 22 (2012),
1607–1633.
[25] W. Gómez and H. Ramı́rez: A ﬁlter algorithm for nonlinear semideﬁnite programming.
Computational and Applied Mathematics, 29 (2010), 297–328.
[26] C. Helmberg, F. Rendl, R.J. Vanderbei, and H. Wolkowicz: An interior-point method
for semideﬁnite programming. SIAM Journal on Optimization, 6 (1996), 342–361.
[27] R.A. Horn and C.R. Johnson: Topics in Matrix Analysis (Cambridge University Press,
Cambridge, 1991).
[28] X.X. Huang, K.L. Teo, and X.Q. Yang: Approximate augmented Lagrangian functions and nonlinear semideﬁnite programs. Acta Mathematica Sinica, English Series,
22 (2006), 1283–1296.
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝58

H. Yamashita & H. Yabe[29] X.X. Huang, X.Q. Yang, and K.L. Teo: A sequential quadratic penalty method for
nonlinear semideﬁnite programming. Optimization, 52 (2003), 715–738.
[30] X.X. Huang, X.Q. Yang, and K.L. Teo: Augmented Lagrangian and nonlinear semideﬁnite programs. Variational Analysis and Applications. Nonconvex Optimization and Its
Applications, 79 (2005), 513–529 .
[31] X.X. Huang, X.Q. Yang, and K.L. Teo: Lower-order penalization approach to nonlinear semideﬁnite programming. Journal of Optimization Theory and Applications, 132(2007), 1–20.
[32] F. Jarre: An interior method for nonconvex semideﬁnite programs. Optimization andEngineering, 1 (2000), 347–372.
[33] F. Jarre: Some aspects of nonlinear semideﬁnite programming. In E.W. Sachs and R.
Tichatschke (Eds.): System Modeling and Optimization XX (Kluwer Academic Publishers, 2003), 55–69.
[34] Y. Kanno and I. Takewaki: Sequential semideﬁnite program for maximum robustness
design of structures under load uncertainty. Journal of Optimization Theory and Applications, 130 (2006), 265–287.
[35] C. Kanzow, I. Ferenczi, and M. Fukushima: On the local convergence of semismooth
Newton methods for linear and nonlinear second-order cone programs without strict
complementarity. SIAM Journal on Optimization, 20 (2009), 297–320.
[36] C. Kanzow, C. Nagel, H. Kato, and M. Fukushima: Successive linearization methods
for nonlinear semideﬁnite programs. Computational Optimization and Applications, 31(2005), 251–273.
[37] A. Kato, H. Yabe, and H. Yamashita: An interior point method with a primal-dual
quadratic barrier penalty function for nonlinear semideﬁnite programming. Journal ofComputational and Applied Mathematics, 275 (2015), 148–161.
[38] M. Kočvara and M. Stingl: PENNON – A code for convex nonlinear and semideﬁnite
programming. Optimization Methods and Software, 18 (2003), 317–333.
[39] M. Kočvara and M. Stingl: PENNON – A generalized augmented Lagrangian method
for semideﬁnite programming. In G. Di Pillo and A. Murli (eds.): High Performance
Algorithms and Software for Nonlinear Optimization (Kluwer Academic Publishers,
Boston, 2003), 303–321.
[40] M. Kočvara and M. Stingl: Solving nonconvex SDP problems of structured optimization
with stability control. Optimization Methods and Software, 19 (2004), 595–609.
[41] M. Kočvara and M. Stingl: On the solution of large-scale SDP problems by the modiﬁed
barrier method using iterative solvers. Mathematical Programming, 109 (2007), 413–
444.
[42] M. Kojima, S. Shindoh, and S. Hara: Interior-point methods for the monotone semidefinite linear complementarity problem in symmetric matrices. SIAM Journal on Optimization, 7 (1997), 86–125.
[43] H. Konno, N. Kawadai, and D. Wu: Estimation of failure probability using semi-deﬁnite
Logit model. Computational Management Science, 1 (2003), 59–73.
[44] F. Leibfritz and W. Lipinski: Description of the benchmark examples in COM P le ib
1.0. Technical Report, Department of Mathematics, University of Trier, D-54286, Trier,Germany, 2003.
[45] F. Leibfritz and E.M.E. Mostafa: An interior point constrained trust region method
for a special class of nonlinear semideﬁnite programming problems. SIAM Journal on
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝A Survey of Numerical Methods for Nonlinear SDP

59

Optimization, 12 (2002), 1048–1074.
[46] C.J. Li and W.Y. Sun: On ﬁlter-successive linearization methods for nonlinear semidefinite programming. Science in China Series A, Mathematics, 52 (2009), 2341–2361.
[47] H.Z. Luo, H.X. Wu, and G.T. Chen: On the convergence of augmented Lagrangian
methods for nonlinear semideﬁnite programming. Journal of Global Optimization, 54(2012), 599–618.
[48] R.D.C. Monteiro: Primal-dual path-following algorithms for semideﬁnite programming.
SIAM Journal on Optimization, 7 (1997), 663–678.
[49] L. Mosheyev and M. Zibulevsky: Penalty/barrier multiplier algorithm for semideﬁnite
programming. Optimization Methods and Software, 13 (2000), 235–261.
[50] Y.E. Nesterov and M.J. Todd: Self-scaled barriers and interior-point methods for convex programming. Mathematics of Operations Research, 22 (1997), 1–42.
[51] Y.E. Nesterov and M.J. Todd: Primal-dual interior-point methods for self-scaled cones.
SIAM Journal on Optimization, 8 (1998), 324–364.
[52] J. Nocedal and S.J. Wright: Numerical Optimization (Springer, New York, 1999).
[53] D. Noll: Local convergence of an augmented Lagrangian method for matrix inequality
constrained programming. Optimization Methods and Software, 22 (2007), 777–802.
[54] D. Noll, M. Torki, and P. Apkarian: Partially augmented Lagrangian method for matrix
inequality constraints. SIAM Journal on Optimization, 15 (2004), 161–184.
[55] R. Polyak: Modiﬁed barrier functions (theory and methods). Mathematical Programming, 54 (1992), 177–222.
[56] H.D. Qi: Local duality of nonlinear semideﬁnite programming. Mathematics of Operations Research, 34 (2009), 124–141.
[57] H. Qi and D. Sun: A quadratically convergent Newton method for computing thenearest correlation matrix. SIAM Journal on Matrix Analysis and Applications, 28(2006), 360–385.
[58] S.M. Robinson: Generalized equations and their solutions, part II: Applications tononlinear programming. Mathematical Programming Study, 19 (1982), 200–221.
[59] A. Shapiro: First and second order analysis of nonlinear semideﬁnite programs. Mathematical Programming, 77 (1997), 301–320.
[60] A. Shapiro and K. Scheinberg: Duality and optimality conditions. In H. Wolkowicz, R.
Saigal and L. Vandenberghe (eds.): Handbook of Semideﬁnite Programming, Theory,
Algorithms and Applications (International Series in Operations Research and Management Science, Kluwer, Boston, 2000), 67–110.
[61] A. Shapiro and J. Sun: Some properties of the augmented Lagrangian in cone constrained optimization. Mathematics of Operations Research, 29 (2004), 479–491.
[62] M. Stingl: On the Solution of Nonlinear Semideﬁnite Programs by Augmented Lagrangian Methods. Doctoral Thesis, University of Erlangen, 2005.
[63] M. Stingl, M. Kočvara, and G. Leugering: A sequential convex semideﬁnite programming algorithm with an application to multiple-load free material optimization. SIAM
Journal on Optimization, 20 (2009), 130–155.
[64] D. Sun, J. Sun, and L. Zhang: The rate of convergence of the augmented Lagrangian method for nonlinear semideﬁnite programming. Mathematical Programming,
114 (2008), 349–391.
c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝60

H. Yamashita & H. Yabe[65] J. Sun, L.W. Zhang, and Y. Wu: Properties of the augmented Lagrangian in nonlinear semideﬁnite optimization. Journal of Optimization Theory and Applications, 129(2006), 437–456.
[66] M.J. Todd: Semideﬁnite optimization. Acta Numerica, 10 (2001), 515–560.
[67] M.J. Todd, K.C. Toh, and R.H. Tütüncü: On the Nesterov-Todd direction in semidefinite programming. SIAM Journal on Optimization, 8 (1998), 769–796.
[68] L. Vandenberghe and S. Boyd: Semideﬁnite programming. SIAM Review, 38 (1996),
49–95.
[69] L. Vandenberghe, S. Boyd, and S.P. Wu: Determinant maximization with linear matrix
inequality constraints. SIAM Journal on Matrix Analysis and Applications, 19 (1998),
499–533.
[70] H. Wolkowicz, R. Saigal, and L. Vandenberghe (eds.): Handbook of Semideﬁnite Programming – Theory, Algorithms and Applications (International Series in Operations
Research and Management Science, Kluwer, Boston, 2000).
[71] H. Wu, H. Luo, X. Ding, and G. Chen: Global convergence of modiﬁed augmented
Lagrangian methods for nonlinear semideﬁnite programming. Computational Optimization and Applications, 56 (2013), 531–558.
[72] Y. Yamakawa and N. Yamashita: Diﬀerentiable merit function for shifted perturbed
Karush-Kuhn-Tucker conditions of nonlinear semideﬁnite programming. to appear in
Paciﬁc Journal of Optimization.
[73] H. Yamashita and H. Yabe: An interior point method with a primal-dual quadratic
barrier penalty function for nonlinear optimization. SIAM Journal on Optimization, 14(2003), 479–499.
[74] H. Yamashita and H. Yabe: A primal-dual interior point method for nonlinear optimization over second-order cones. Optimization Methods and Software, 24 (2009),
407–426.
[75] H. Yamashita and H. Yabe: Local and superlinear convergence of a primal-dual interior
point method for nonlinear semideﬁnite programming. Mathematical Programming, 132(2012), 1–30.
[76] H. Yamashita, H. Yabe, and K. Harada: A primal-dual interior point method for
nonlinear semideﬁnite programming. Mathematical Programming, 135 (2012), 89–121.
[77] L. Yang and B. Yu: A homotopy method for nonlinear semideﬁnite programming.
Computational Optimization and Applications, 56 (2013), 81–96.
[78] Z.B. Zhu and H.L. Zhu: A ﬁlter method for nonlinear semideﬁnite programming with
global convergence. Acta Mathematica Sinica, English Series, 30 (2014), 1810–1826.
Hiroshi Yabe
Department of Mathematical Information Science
Tokyo University of Science
1-3 Kagurazaka, Shinjuku-ku
Tokyo 162-8601, Japan
E-mail: yabe@rs.kagu.tus.ac.jp

c by ORSJ. Unauthorized reproduction of this article is prohibited.
Copyright ⃝