A sequential semidefinite programming method and an application in passive reduced-order modeling, ,

A sequential semidefinite programming method and an application in passive reduced-order modeling

Roland W. Freund Department of Mathematics University of California at Davis One Shields Avenue Davis, California 95616, U.S.A. e-mail : freund@math.ucdavis.edu Florian Jarre and Christoph Vogelbusch Institut für Mathematik Universität Düsseldorf Universitätsstraße 1 D–40225 Düsseldorf, Germany e-mail : {jarre,vogelbuc}@opt.uni-duesseldorf.de

Abstract. We consider the solution of nonlinear programs with nonlinear semidefiniteness constraints. The need for an efficient exploitation of the cone of positive semidefinite matrices makes the solution of such nonlinear semidefinite programs more complicated than the solution of standard nonlinear programs. In particular, a suitable symmetrization procedure needs to be chosen for the linearization of the complementarity condition. The choice of the symmetrization procedure can be shifted in a very natural way to certain linear semidefinite subproblems, and can thus be reduced to a well-studied problem. The resulting sequential semidefinite programming (SSP) method is a generalization of the well-known SQP method for standard nonlinear programs. We present a sensitivity result for nonlinear semidefinite programs, and then based on this result, we give a self-contained proof of local quadratic convergence of the SSP method. We also describe a class of nonlinear semidefinite programs that arise in passive reduced-order modeling, and we report results of some numerical experiments with the SSP method applied to problems in that class.

Key words. semidefinite programming, nonlinear programming, sequential quadratic programming, quadratic semidefinite program, sensitivity, convergence, reduced-order modeling, passivity, positive realness

1 Introduction

In recent years, interior-point methods for solving linear semidefinite programs (SDPs) have received a lot of attention, and as a result, these methods are now very well developed; see, e.g., [30, 32] , the papers in [35] , and the references given there. At each iteration of an interior-point method, the complementarity condition is relaxed, symmetrized, and linearized.

Various symmetrization operators are known. The choice of the symmetrization operator and of the relaxation parameter determine the step length at each iteration, and thus the efficiency of the overall method.

In this paper, we are concerned with the solution of nonlinear semidefinite programs (NLSDPs). Interior-point methods for linear SDPs can be extended to NLSDPs. However, some additional difficulties arise. First, the step length now also depends on the quality of the linearization of the nonlinear functions. Second, the choice of the symmetrization procedure is considerably more complicated than in the linear case since the system matrix is no longer positive semidefinite. To address these difficulties, in this paper, we consider an approach that separates the linearization and the symmetrization in a natural way, namely a generalization of the sequential quadratic programming (SQP) method for standard nonlinear programs. Such a generalization has already been mentioned by Robinson [26] within the more general framework of nonlinear programs over convex cones. This framework includes NLSDPs as a special case.

While Robinson did not discuss implementational issues of such a generalized SQP approach, the recent progress in the solution of linear SDPs makes this approach especially interesting for the solution of NLSDPs. We first present a derivation of a generalized SQP method, namely the sequential semidefinite programming (SSP) method, for solving NLSDPs. In order to analyze the convergence of the SSP method, we present a sensitivity result for certain local optimal solutions of general, possibly nonconvex, quadratic semidefinite programs. We then use this result to derive a self-contained proof of local quadratic convergence of the SSP method under the assumptions that the optimal solution is locally unique, strictly complementary, and satisfies a second-order sufficient condition.

One of the first numerical approaches for solving a class of NLSDPs was given in [23, 24] .

Other recent approaches for solving NLSDPs are the program package LOQO [33] based on a primal-dual method; see also [34] . Another promising approach for solving large-scale SDPs is the modified-barrier method proposed in [17] . This modified-barrier approach does not require the barrier parameter to converge to zero, and may thus overcome some of the problems related to ill-conditioning in traditional interior-point methods. Further approaches to solving NLSDPs have been presented in [11, 12, 31] . In [11] , the augmented Lagrangian method is applied to NLSDPs, while the approach proposed in [12] is based on an SQP method generalized to NLSDPs. The paper [12] also contains a proof of local quadratic convergence. However, in contrast to this paper, the algorithm [12] is not derived from a comparison with interior-point algorithms, and the proof of convergence does not use any differentiability properties of the optimal solutions. In [10] , Correa and Ramírez present a proof of global convergence of a modification of the method proposed in [12] . The modification employs certain merit functions to control the step lengths of the SQP algorithm.

The remainder of this paper is organized as follows. In Section 2 , we introduce some notation. In Section 3 , we describe a class of nonlinear semidefinite programs that arise in passive reduced-order modeling. In Section 4 , we recall known results for linear SDPs in a form that can be easily transferred to NLSDPs. In Section 5 , we discuss primal-dual systems for NLSDPs, and in Section 6 , the SSP method is introduced as a generalized SQP method. In Section 7 , we present sensitivity results, first for a certain class of quadratic SDPs, and then for general NLSDPs. Based on these sensitivity results, in Section 8 , we give a self-contained proof of local quadratic convergence of the SSP method. In Section 9 , we present results of some numerical experiments. Finally, in Section 10 , we make some concluding remarks.

2 Notation

Throughout this article, all vectors and matrices are assumed to have real entries. As usual,

Y^{T} = [\begin{matrix} y_{j i} \end{matrix}]

denotes the transpose of the matrix

Y = [\begin{matrix} y_{i j} \end{matrix}]

. The vector norm

∥ x ∥ : = \sqrt{x^{T} x}

is always the Euclidean norm and

∥ Y ∥ : = {max}_{∥ x ∥ = 1} ∥ Y x ∥

is the corresponding matrix norm. For vectors

x \in R^{n}

x \geq 0

means that all entries of

x

are nonnegative, and

D i a g (x)

denotes the

n \times n

diagonal matrix the diagonal entries of which are the entries of

x

. The

n \times n

identity matrix is denoted by

I_{n}

The trace inner product on the space of real

n \times m

matrices is given by

〈 Z, Y 〉 : = Z ∙ Y : = t r a c e (Z^{T} Y) = \sum_{i = 1}^{n} \sum_{j = 1}^{m} z_{i j} y_{i j}

for any pair

Y = [\begin{matrix} y_{i j} \end{matrix}]

and

Z = [\begin{matrix} z_{i j} \end{matrix}]

n \times m

matrices. The space of real symmetric

m \times m

matrices is denoted by

S^{m}

. The notation

Y ≽ 0

(

Y ≻ 0

) is used to indicate that

Y \in S^{m}

is symmetric positive semidefinite (positive definite).

Semidefiniteness constraints are expressed by means of matrix-valued functions from

R^{n}

S^{m}

. We use the symbol

A : R^{n} \to S^{m}

if such a function is linear, and the symbol

ℬ : R^{n} \to S^{m}

if such a function is nonlinear.

Note that any linear function

A : R^{n} \to S^{m}

can be expressed in the form

\begin{matrix} A (x) = \sum_{i = 1}^{n} x_{i} A^{(i)} for all x \in R^{n}, \end{matrix}

(1)

with symmetric matrices

A^{(i)} \in S^{m}

i = 1, 2, . . ., n

. Based on the representation ( 1 ) we introduce the norm

\begin{matrix} ∥ A ∥ : = {(\sum_{i = 1}^{n} {∥ A^{(i)} ∥}^{2})}^{\frac{1}{2}} \end{matrix}

(2)

A

. The adjoint operator

A^{*} : S^{m} \to R^{n}

with respect to the trace inner product is defined by

〈 A (x), Y 〉 = 〈 x, A^{*} (Y) 〉 = x^{T} A^{*} (Y) for all x \in R^{n} and Y \in S^{m} .

It turns out that

\begin{matrix} A^{*} (Y) = [\begin{matrix} A^{(1)} ∙ Y & . . . & A^{(n)} ∙ Y \end{matrix}] for all Y \in S^{m} . \end{matrix}

(3)

We always assume that nonlinear functions

ℬ : R^{n} \to S^{m}

are at least

C^{2}

-differentiable. We denote by

B^{(i)} (x) : = \frac{\partial}{\partial x_{i}} ℬ (x) and B^{(i, j)} (x) : = \frac{\partial^{2}}{\partial x_{i} \partial x_{j}} ℬ (x), i, j = 1, 2, . . ., n,

the first and second partial derivatives of

ℬ

, respectively. For each

x \in R^{n}

, the derivative

D_{x} ℬ

x

induces a linear function

D_{x} ℬ (x) : R^{n} \to S^{m}

, which is given by

D_{x} ℬ (x) [Δ x] : = \sum_{i = 1}^{n} (Δ x)_{i} B^{(i)} (x) \in S^{m} for all Δ x \in R^{n} .

In particular,

ℬ (x + Δ x) \approx ℬ (x) + D_{x} ℬ (x) [Δ x], Δ x \in R^{n},

is the linearization of

ℬ

at the point

x

. For any linear function

A : R^{n} \to S^{m}

, we have

\begin{matrix} D_{x} A (x) [Δ x] = A (Δ x) for all x, Δ x \in R^{n} . \end{matrix}

(4)

We always use the expression on the right-hand side of ( 4 ) to describe derivatives of linear functions.

We remark that for any fixed matrix

Y \in S^{m}

, the map

x \mapsto ℬ (x) ∙ Y

is a scalar-valued function of

x \in R^{n}

. Its gradient at

x

is given by

\begin{matrix} \nabla_{x} (ℬ (x) ∙ Y) = {(D_{x} (ℬ (x) ∙ Y))}^{T} = [\begin{matrix} B^{(1)} (x) ∙ Y & . . . & B^{(n)} (x) ∙ Y \end{matrix}] \in R^{n} \end{matrix}

(5)

and its Hessian by

\nabla_{x}^{2} (ℬ (x) ∙ Y) = [\begin{matrix} B^{(1, 1)} (x) ∙ Y & \cdot \cdot \cdot & B^{(1, n)} (x) ∙ Y & . . . & . . . & B^{(n, 1)} (x) ∙ Y & \cdot \cdot \cdot & B^{(n, n)} (x) ∙ Y \end{matrix}] \in S^{n} .

In particular, for any linear function

A : R^{n} \to S^{m}

, in view of ( 1 ), ( 3 ), and ( 5 ), we have

\begin{matrix} \nabla_{x} (A (x) ∙ Y) = A^{*} (Y) . \end{matrix}

(6)

3 An application in passive reduced-order modeling

We remark that applications of linear SDPs include relaxations of combinatorial optimization problems and problems related to Lyapunov functions or the positive real lemma in control theory; we refer the reader to [1, 3, 8, 14, 30, 32] and the references given there. In this section, we describe an application in passive reduced-order modeling that leads to a class of NLSDPs. Roughly speaking, a system is called passive if it does not generate energy. For the special case of time-invariant linear dynamical systems, passivity is equivalent to positive realness of the frequency-domain transfer function associated with the system. More precisely, consider transfer functions of the form

\begin{matrix} Z_{n} (s) = B_{2}^{T} {(G + s C)}^{- 1} B_{1}, s \in C, \end{matrix}

(7)

where

G, C \in R^{n \times n}

and

B_{1}, B_{2} \in R^{n \times m}

are given data matrices. The integer

n

is the state-space dimension of the time-invariant linear dynamical system, and

m

is the number of inputs and outputs of the system. In ( 7 ), the matrix pencil

G + s C

is assumed to be regular, i.e., the matrix

G + s C

is singular for only finitely many values of

s \in C

. Note that

Z_{n}

is an

m \times m

-matrix-valued rational function of the complex variable

s \in C

In reduced-order modeling, one is given a large-scale time-invariant linear dynamical system of state-space dimension

N

, and the problem is to construct a “good” approximation of that system of state-space dimension

n ≪ N

; see, e.g., [13] and the references given there. If the large-scale system is passive, then for certain applications, it is crucial that the reduced-order model of state-space dimension

n

preserves the passivity of the original system. Unfortunately, some of the most efficient reduced-order modeling techniques do not preserve passivity. However, the reduced-order models are often “almost” passive, and passivity of the models can be enforced by perturbing the data matrices of the models. Next, we describe how the problem of constructing such perturbations leads to a class of NLSDPs. An

m \times m

-matrix-valued rational function

Z

is called positive real if the following three conditions are satisfied

:

(i) $Z$ is analytic in $C_{+} : = {s \in C | R e (s) > 0}$ ;
(ii) $Z (\bar{s}) = \bar{Z (s)}$ for all $s \in C$ ;
(iii) $Z (s) + {(Z (\bar{s}))}^{T} ≽ 0$ for all $s \in C_{+}$ .

For functions

Z_{n}

of the form ( 7 ) positive realness (and thus passivity of the system associated with

Z_{n}

) can be characterized via linear SDPs; see, e.g., [3, 14] and the references given there. More precisely, if the linear SDP

\begin{matrix} \begin{matrix} P^{T} G + G^{T} P & ≽ 0, \\ P^{T} C = C^{T} P & ≽ 0, \\ P^{T} B_{1} & = B_{2}, \end{matrix} \end{matrix}

(8)

has a solution

P \in R^{n \times n}

, then the transfer function ( 7 ),

Z_{n}

, is positive real. Conversely, under certain additional assumptions (see [14] ), positive realness of

Z_{n}

implies the solvability of the linear SDP ( 8 ).

Now assume that

Z_{n}

in ( 7 ) is the transfer function of a non-passive reduced-oder model of a passive large-scale system. Our goal is to perturb some of data matrices in ( 7 ) so that the perturbed transfer function is positive real. For the special case

m = 1

, such an approach is discussed in [5] . In this case, there is a simple eigenvalue-based characterization [4] of positive realness. However, this characterization cannot be extended to the general case

m \geq 1

Another special case, which leads to linear SDPs, is described in [9] .

In the general case

m \geq 1

, we employ perturbations

X_{G}

and

X_{C}

of the matrices

G

and

C

in ( 7 ). The resulting perturbed transfer function is then of the form

\begin{matrix} {\tilde{Z}}_{n} (s) = B_{2}^{T} {(G + X_{G} + s (C + X_{C}))}^{- 1} B_{1}, \end{matrix}

(9)

and the problem is to construct the perturbations

X_{G}

and

X_{C}

such that

{\tilde{Z}}_{n}

is positive real.

Applying the characterization ( 8 ) of positive realness to ( 9 ), we obtain the following nonconvex NLDSP:

\begin{matrix} \begin{matrix} P^{T} (G + X_{G}) + (G + X_{G})^{T} P & ≽ 0, \\ P^{T} (C + X_{C}) = (C + X_{C})^{T} P & ≽ 0, \\ P^{T} B_{1} & = B_{2} . \end{matrix} \end{matrix}

(10)

Here, the unknowns are the matrices

P, X_{G}, X_{C} \in R^{n \times n}

. If ( 10 ) has a solution

P, X_{G}, X_{C}

, then choosing the matrices

X_{G}

and

X_{C}

as the perturbations in ( 9 ) guarantees passivity of the reduced-order model given by the transfer function

{\tilde{Z}}_{n}

4 Linear semidefinite programs

In this section, we briefly review the case of linear semidefinite programs.

Given a linear function

A : R^{n} \to S^{m}

, a vector

b \in R^{n}

, and a matrix

C \in S^{m}

, a pair of primal and dual linear semidefinite programs is as follows:

\begin{matrix} \begin{matrix} maximize C ∙ Y subject to & Y \in S^{m}, Y ≽ 0, \\ A^{*} (Y) + b = 0, \end{matrix} \end{matrix}

(11)

and

\begin{matrix} \begin{matrix} minimize b^{T} x subject to & x \in R^{n}, \\ A (x) + C ≼ 0 . \end{matrix} \end{matrix}

(12)

We remark that this formulation is a slight variation of the standard pair of primal-dual programs. We chose the above version in order to facilitate the generalization of problems of the form ( 12 ) to nonlinear semidefinite programs in standard form.

If there exists a matrix

Y ≻ 0

that is feasible for ( 11 ), then we call

Y

strictly feasible for ( 11 ) and say that ( 11 ) satisfies Slater's condition. Likewise, if there exists a vector

x

such that

A (x) + C ≺ 0

, then we call ( 12 ) strictly feasible and say that ( 12 ) satisfies Slater's condition.

The following optimality conditions for linear semidefinite programs are well known; see, e.g., [27] . If problem ( 11 ) or ( 12 ) satisfies Slater's condition, then the optimal values of ( 11 ) and ( 12 ) coincide. Furthermore, if in addition both problems are feasible, then optimal solutions

Y^{o p t}

and

x^{o p t}

of both problems exist and

Y : = Y^{o p t}

and

x : = x^{o p t}

satisfy the complementarity condition

\begin{matrix} Y S = 0, where S : = - C - A (x) . \end{matrix}

(13)

Conversely, if

Y

and

x

are feasible points for ( 11 ) and ( 12 ), respectively, and satisfy the complementarity condition ( 13 ), then

Y^{o p t} : = Y

is an optimal solution of ( 11 ) and

x^{o p t} : = x

is an optimal solution of ( 12 ).

These optimality conditions can be summarized as follows. If problem ( 12 ) satisfies Slater's condition, then for a point

x \in R^{n}

to be an optimal solution of ( 12 ) it is necessary and sufficient that there exist matrices

Y ≽ 0

and

S ≽ 0

such that

\begin{matrix} \begin{matrix} A (x) + C + S & = 0, \\ A^{*} (Y) + b & = 0, \\ Y S & = 0 . \end{matrix} \end{matrix}

(14)

Note that, in view of ( 6 ), the second equation in ( 14 ) can also be written in the form

\begin{matrix} \nabla_{x} ((A (x) + C) ∙ Y) + b = A^{*} (Y) + b = 0 . \end{matrix}

(15)

Furthermore, the last equation in ( 14 ) is equivalent to its symmetric form,

Y S + S Y = 0

; see, e.g., [2] . In the case of strict complementarity, the derivatives of

Y S = 0

and

Y S + S Y = 0

are also equivalent. For later use, we state these facts in the following lemma.

Lemma 1 Let

Y, S \in S^{m}

a) If $Y ≽ 0$ or $S ≽ 0$ , then
$\begin{matrix} Y S = 0 ⟺ Y S + S Y = 0 . \end{matrix}$ (16)

b) If

Y

and

S

are strictly complementary, i.e.,

Y, S ≽ 0

Y S = 0

, and

Y + S ≻ 0

, then for any

\dot{Y}, \dot{S} \in S^{m}

\begin{matrix} Y \dot{S} + \dot{Y} S = 0 ⟺ Y \dot{S} + \dot{Y} S + \dot{S} Y + S \dot{Y} = 0 . \end{matrix}

(17)

Moreover,

Y, S

have representations of the form

\begin{matrix} Y = U [\begin{matrix} Y_{1} & 0 & 0 & 0 \end{matrix}] U^{T}, S = U [\begin{matrix} 0 & 0 & 0 & S_{2} \end{matrix}] U^{T}, \end{matrix}

(18)

where

U

is an

m \times m

orthogonal matrix,

Y_{1} ≻ 0

is a

k \times k

diagonal matrix, and

S_{2} ≻ 0

is an

(m - k) \times (m - k)

diagonal matrix, and any matrices

\dot{Y}, \dot{S} \in S^{m}

satisfying 17

()

are of the form

\begin{matrix} \dot{Y} = U [\begin{matrix} {\dot{Y}}_{1} & {\dot{Y}}_{3} & {\dot{Y}}_{3}^{T} & 0 \end{matrix}] U^{T}, \dot{S} = U [\begin{matrix} 0 & {\dot{S}}_{3} & {\dot{S}}_{3}^{T} & {\dot{S}}_{2} \end{matrix}] U^{T}, where Y_{1} {\dot{S}}_{3} + {\dot{Y}}_{3} S_{2} = 0 . \end{matrix}

(19)

Proof. The equivalence ( 16 ) is well known; see, e.g., [2,Page749] .

We now turn to the proof of part b). The strict complementarity of

Y

and

S

readily implies that

Y

and

S

have representations of the form ( 18 ); see, e.g., [18,Page62] . Any matrices

\dot{Y}, \dot{S} \in S^{m}

can be written in the form

\begin{matrix} \dot{Y} = U [\begin{matrix} {\dot{Y}}_{1} & {\dot{Y}}_{3} & {\dot{Y}}_{3}^{T} & {\dot{Y}}_{2} \end{matrix}] U^{T}, \dot{S} = U [\begin{matrix} {\dot{S}}_{1} & {\dot{S}}_{3} & {\dot{S}}_{3}^{T} & {\dot{S}}_{2} \end{matrix}] U^{T}, \end{matrix}

(20)

where

U

is the matrix from ( 18 ) and the block sizes in ( 20 ) are the same as in ( 18 ). Using ( 18 ) and ( 20 ), it follows that the equation on the left-hand side of ( 17 ) is satisfied if, and only if,

Y_{1} {\dot{S}}_{1} = 0, {\dot{Y}}_{2} S_{2} = 0, Y_{1} {\dot{S}}_{3} + {\dot{Y}}_{3} S_{2} = 0 .

Since

Y_{1}

and

S_{2}

are in particular nonsingular, the first two relations imply

{\dot{S}}_{1} = 0

and

{\dot{Y}}_{2} = 0

Thus, any matrices

\dot{Y}, \dot{S} \in S^{m}

satisfying the equation on the left-hand side of ( 17 ) are of the form ( 19 ). Similarly, using ( 18 ) and ( 20 ), it follows that the equation on the right-hand side of ( 17 ) is satisfied if, and only if,

Y_{1} {\dot{S}}_{1} + {\dot{S}}_{1} Y_{1} = 0, {\dot{Y}}_{2} S_{2} + S_{2} {\dot{Y}}_{2} = 0, Y_{1} {\dot{S}}_{3} + {\dot{Y}}_{3} S_{2} = 0 .

Since

Y_{1} ≻ 0

and

S_{2} ≻ 0

, the first two relations imply

{\dot{S}}_{1} = 0

and

{\dot{Y}}_{2} = 0

, and so

\dot{Y}, \dot{S}

are again of the form ( 19 ).

5 Nonlinear semidefinite programs

In this section, we consider nonlinear semidefinite programs, which are extensions of the dual linear semidefinite programs ( 12 ).

Given a vector

b \in R^{n}

and a matrix-valued function

ℬ : R^{n} \to S^{m}

, we consider problems of the following form:

\begin{matrix} \begin{matrix} minimize b^{T} x subject to & x \in R^{n}, \\ ℬ (x) ≼ 0 . \end{matrix} \end{matrix}

(21)

Here, the function

ℬ

is nonlinear in general, and thus ( 21 ) represents a class of nonlinear semidefinite programs. We assume that the function

ℬ

is at least

C^{2}

-differentiable.

For simplicity of presentation, we have chosen a simple form of problem ( 21 ). We stress that problem ( 21 ) may also include additional nonlinear equality and inequality constraints.

The corresponding modifications are detailed at the end of this paper. Furthermore, the choice of the linear objective function

b^{T} x

in ( 21 ) was made only to simplify notation. A nonlinear objective function can always be transformed into a linear one by adding one artificial variable and one more constraint. In particular, all statements about ( 21 ) in this paper can be modified so that they apply to additional nonlinear equality and inequality constraints and to nonlinear objective functions.

Note that the class ( 21 ) reduces to linear semidefinite programs of the form ( 12 ) if

ℬ

is an affine function.

The Lagrangian

ℒ : R^{n} \times S^{m} \to R

of ( 21 ) is defined as follows:

\begin{matrix} ℒ (x, Y) : = b^{T} x + ℬ (x) ∙ Y . \end{matrix}

(22)

Its gradient with respect to

x

is given by

\begin{matrix} g (x, Y) : = \nabla_{x} ℒ (x, Y) = b + \nabla_{x} (ℬ (x) ∙ Y) \end{matrix}

(23)

and its Hessian by

\begin{matrix} H (x, Y) : = \nabla_{x}^{2} ℒ (x, Y) = \nabla_{x}^{2} (ℬ (x) ∙ Y) . \end{matrix}

(24)

If the problem ( 21 ) is convex and satisfies Slater's condition [19] , then for each optimal solution

x

of ( 21 ) there exists an

m \times m

matrix

Y ≽ 0

such that the pair

(x, Y)

is a saddle point of the Lagrangian ( 22 ),

ℒ

More generally, for nonconvex problems ( 21 ), let

x \in R^{n}

be a feasible point of ( 21 ), and assume that the Robinson or Mangasarian-Fromovitz constraint qualification [19, 25, 26] is satisfied at

x

, i.e., there exists a vector

Δ x \neq 0

such that

ℬ (x) + D_{x} ℬ (x) [Δ x] ≺ 0

. Then, if

x

is a local minimizer of ( 21 ), the first-order optimality condition is satisfied, i.e., there exist matrices

Y, S \in S^{m}

such that

\begin{matrix} \begin{matrix} ℬ (x) + S & = 0, \\ g (x, Y) & = 0, \\ Y S & = 0, \\ Y, S & ≽ 0 . \end{matrix} \end{matrix}

(25)

The system ( 25 ) is a straightforward generalization of the optimality conditions ( 14 ) and ( 15 ), with the affine function

A (x) + C

in ( 15 ) replaced by the nonlinear function

ℬ (x)

. Primal-dual interior-point methods for solving ( 21 ) roughly proceed as follows. For some sequence of duality parameters

μ_{k} > 0

μ_{k} \to 0

, the solutions of the perturbed primal-dual system,

\begin{matrix} \begin{matrix} ℬ (x) + S & = 0, \\ g (x, Y) & = 0, \\ Y S & = μ_{k} I_{m}, \\ Y, S & ≻ 0, \end{matrix} \end{matrix}

(26)

are approximated by some variant of Newton's method. Since Newton's method does not preserve any inequalities, the parameters

μ_{k} > 0

are used to maintain strict feasibility, i.e.,

Y, S ≻ 0

for all iterates.

The solutions of ( 26 ) coincide with the solutions of the standard logarithmic-barrier problems for ( 21 ). Moreover, the logarithmic-barrier approach for solving ( 21 ) can be interpreted as a certain choice of the `symmetrization operator' for the equation,

Y S = μ_{k} I_{m}

, in the third row of ( 26 ); see Section 6 below. With this choice, the barrier function yields a very natural criterion for the step-size control in trust-region algorithms. The authors have implemented various versions of predictor-corrector trust-region barrier methods for solving ( 21 ). For a number of examples, the running times of the resulting algorithms were comparable to the behavior of interior-point methods for convex programs. However, the authors also encountered several instances in which the number of iterations for these methods was very high compared to the typical number of iterations needed for solving linear SDPs. For such negative examples it may be more efficient to solve a sequence of linear SDPs in order to obtain an approximate solution of ( 21 ). This observation motivated the SSP method described in the next section.

6 An SSP method for nonlinear semidefinite programs

In this section, we introduce the sequential semidefinite programming (SSP) method, which is a generalization of the SQP method for standard nonlinear programs to nonlinear semidefinite programs of the form ( 21 ). For an overview of SQP methods for standard nonlinear programs, we refer the reader to [6] and the references given there.

In analogy to the SQP method, at each iteration of the SSP method one solves a subproblem that is slightly more difficult than the linearization of ( 25 ) at the current iterate. More precisely, let

(x^{k}, Y^{k})

denote the current point at the beginning of the

k

-th iteration. One then determines corrections

(Δ x, Δ Y)

and a matrix

S

such that

\begin{matrix} \begin{matrix} ℬ (x^{k}) + D_{x} ℬ (x^{k}) [Δ x] + S & = 0, \\ b + H^{k} Δ x + \nabla_{x} (ℬ (x^{k}) ∙ (Y^{k} + Δ Y)) & = 0, \\ (Y^{k} + Δ Y^{k}) S & = 0, \\ Y^{k} + Δ Y, S & ≽ 0 . \end{matrix} \end{matrix}

(27)

Here and in the sequel, we use the notation

\begin{matrix} H^{k} : = H (x^{k}, Y^{k}) . \end{matrix}

(28)

Recall from ( 23 ) and ( 24 ) that

g (x, Y)

and

H (x, Y)

denote the gradient and Hessian with respect to

x

, respectively, of the Lagrangian ( 22 ),

ℒ (x, Y)

, of the nonlinear semidefinite program ( 21 ). Moreover, from ( 23 ) it follows that the linearization of

g (x, Y)

at the point

(x^{k}, Y^{k})

is given by

\begin{matrix} g (x^{k} + Δ x, Y^{k} + Δ Y) & \approx b + H^{k} Δ x + \nabla_{x} (ℬ (x^{k}) ∙ (Y^{k} + Δ Y)) . \end{matrix}

Thus the second equation in ( 27 ) is just the linearization of the second equation in ( 25 ).

Furthermore, the first equation of ( 27 ) is a straightforward linearization of the first equation in ( 25 ). This linearization is used in the same way in primal-dual interior methods.

The last two rows in ( 27 ) and ( 25 ) are identical when

Y

in ( 25 ) is rewritten as

Y = Y^{k} + Δ Y

In analogy to SQP methods for standard nonlinear programs, the problem of how to guarantee the nonnegativity constraints, namely

ℬ (x) ≼ 0

, is thus shifted to the subproblem ( 27 ). If the iterates

x^{k}

generated from ( 27 ) converge, then their limit

x

automatically satisfies

ℬ (x) ≼ 0

In contrast, interior methods use perturbations, symmetrizations, and linearizations for the last two rows in ( 25 ), resulting in cheaper linear subproblems that are typically less `powerful' than the subproblems ( 27 ). An important aspect for interior methods for linear SDPs is the choice of the symmetrization procedure for the bilinear equation

Y S = μ_{k} I_{m}

; see, e.g., [20] .

For convex SDPs, theoretical convergence analyses are well developed and also supported by numerical evidence of rapid convergence. However, the generalization of these convergence results to nonlinear nonconvex subproblems is far from obvious. The proposed SSP method allows to apply the symmetrization to linear SDPs, thus reducing this aspect to a well-studied topic.

In summary, both the problem of choosing a suitable symmetrization scheme and the problem of how to guarantee the nonnegativity constraints are shifted to the subproblem ( 27 ).

Note that the conditions ( 27 ) are the optimality conditions for the problem

\begin{matrix} \begin{matrix} minimize b^{T} Δ x + \frac{1}{2} (Δ x)^{T} H^{k} Δ x subject to & Δ x \in R^{n}, \\ ℬ (x^{k}) + D_{x} ℬ (x^{k}) [Δ x] ≼ 0 . \end{matrix} \end{matrix}

(29)

The conditions ( 27 ) and ( 29 ) have been considered in [26,Equations(2.1)and(2.2)] , with the remark that they have “been found to be an appropriate approximation of ” ( 21 ) “for numerical purposes”.

In order to be able to solve the subproblem ( 27 ) efficiently, in practice, one replaces the matrix

H^{k}

in ( 27 ), respectively ( 29 ), by a positive semidefinite approximation

{\hat{H}}^{k}

H^{k}

. As in the case of standard SQP methods, a BFGS update for the Hessian of the Lagrangian ( 22 ),

ℒ

, can be used to approximate

H^{k}

by some positive semidefinite matrix

{\hat{H}}^{k}

. Given an estimate

{\hat{H}}^{k}

H^{k}

for the current,

k

-th, SSP iteration, the quasi-Newton condition to generate a BFGS update

{\hat{H}}^{k + 1}

approximating the matrix

H^{k + 1}

for the next,

(k + 1)

-st, SSP iteration can be derived as follows:

\begin{matrix} \begin{matrix} {\hat{H}}^{k + 1} Δ x & = \nabla_{x} (ℬ (x^{k + 1}) ∙ Y^{k + 1}) - \nabla_{x} (ℬ (x^{k}) ∙ Y^{k + 1}) \\ = \nabla_{x} ℒ (x^{k + 1}, Y^{k + 1}) - \nabla_{x} ℒ (x^{k}, Y^{k + 1}) \\ \approx \nabla_{x}^{2} ℒ (x^{k + 1}, Y^{k + 1}) (x^{k + 1} - x^{k}) . \end{matrix} \end{matrix}

(30)

{\hat{H}}^{k}

is positive semidefinite, the BFGS update with the above condition can be suitably damped such that

{\hat{H}}^{k + 1}

is also positive semidefinite. At each iteration of the SSP method, problem ( 29 ) is solved with

H^{k}

replaced by the matrix

{\hat{H}}^{k + 1}

that is obtained by the BFGS update of

{\hat{H}}^{k}

from the previous SSP iteration. If

{\hat{H}}^{k + 1}

is positive semidefinite, problem ( 29 ) essentially reduces to a linear SDP, since the convex quadratic term in the objective function can be written as a semidefiniteness constraint or a second-order-cone constraint.

While the formulation as a second-order-cone constraint is more efficient, and for example, can be specified as input for the program package SeDuMi [28] in order to solve ( 29 ), it was pointed out by [22] that it may be most efficient to use a program that is designed for SDPs with linear constraints and a convex quadratic objective function.

It seems that many results for standard SQP methods carry over in a straightforward fashion to the SSP method. For example, the SSP method can be augmented by a penalty term in case that the subproblems ( 29 ) become infeasible. In this case, the right-hand sides, “

0

”, of the first three rows in ( 27 ) are replaced by weaker, penalized right-hand sides. Moreover, the convergence analysis of the method proposed in [10, 12] yields results that are comparable to the ones for standard SQP methods.

The standard analysis of quadratic convergence of SQP methods for nonlinear programs that satisfy strict complementarity conditions proceeds by first showing that the active constraints will be identified correctly in the final stages of the algorithm and then using the equivalence of the SQP iteration and the Newton iteration for the simplified KKT-system in which only the active constraints are used.

For nonlinear semidefinite programs the situation is slightly more complicated since it is difficult to identify active constraints. The paper [12] presents a proof that is based on a new approach by Bonnans et al. [7] and uses some general results due to Robinson [26] . It does not require strict complementarity and allows for approximate Hessian matrices in ( 30 ).

In the next two sections, we present a more elementary and self-contained approach to analyze convergence of the SSP method under a strict complementarity condition.

7 Sensitivity results

In this section, we establish sensitivity results, first for the special case of quadratic semidefinite programs and then for general nonlinear semidefinite programs of the form ( 21 ).

More precisely, we show that strictly complementary solutions of such problems depend smoothly on the problem data.

We start with quadratic semidefinite programs of the form

\begin{matrix} \begin{matrix} minimize f (x) subject to & x \in R^{n}, \\ A (x) + C ≼ 0 . \end{matrix} \end{matrix}

(31)

Here,

A : R^{n} \to S^{m}

is a linear function,

C \in S^{m}

, and

f : R^{n} \to R

is a quadratic function defined by

f (x) = b^{T} x + \frac{1}{2} x^{T} H x

, where

b \in R^{n}

and

H \in S^{n}

. Note that we make no further assumptions on the matrix

H

. Thus, problem ( 31 ) is a general, possibly nonconvex, quadratic semidefinite program.

The problem ( 31 ) is described by the data

\begin{matrix} D : = [A, b, C, H] . \end{matrix}

(32)

In Theorem 1 below, we present a sensitivity result for the solutions 31

()

when the data

D

is changed to

D + Δ D

where

\begin{matrix} Δ D : = [Δ A, Δ b, Δ C, Δ H] \end{matrix}

(33)

is a sufficiently small perturbation. We use the norm

∥ D ∥ : = {(∥ A ∥^{2} + ∥ b ∥^{2} + ∥ C ∥^{2} + ∥ H ∥^{2})}^{\frac{1}{2}}

for data ( 32 ) and perturbations ( 33 ). Recall that

∥ A ∥

is defined in ( 2 ).

We denote by

ℒ^{(q)} (x, Y) : = f (x) + (A (x) + C) ∙ Y

the Lagrangian of problem ( 31 ). Note that

\nabla_{x} f (x) = b + H x

. Together with ( 6 ), it follows that

\nabla_{x} ℒ^{(q)} (x, Y) = b + H x + A^{*} (Y) and \nabla_{x}^{2} ℒ^{(q)} (x, Y) = H .

Recall that problem ( 31 ) is said to satisfy Slater's condition if there exists a vector

x \in R^{n}

with

A (x) + C ≺ 0

. Moreover, the triple

(\bar{x}, \bar{Y}, \bar{S})

, where

\bar{x} \in R^{n}

and

\bar{Y}, \bar{S} \in S^{m}

, is called a stationary point of ( 31 ) if

\begin{matrix} \begin{matrix} A (\bar{x}) + C + \bar{S} & = 0, \\ b + H \bar{x} + A^{*} (\bar{Y}) & = 0, \\ \bar{Y} \bar{S} + \bar{S} \bar{Y} & = 0, \\ \bar{Y}, \bar{S} & ≽ 0 . \end{matrix} \end{matrix}

(34)

Here, we have used equivalence ( 16 ) of Lemma 1 and replaced

\bar{Y} \bar{S} = 0

by its symmetric version, which is stated as the third equation of ( 34 ). If in addition to ( 34 ), one has

\begin{matrix} \bar{Y} + \bar{S} ≻ 0, \end{matrix}

(35)

then

(\bar{x}, \bar{Y}, \bar{S})

is said to be a strictly complementary stationary point of 31

()

. Let

\bar{x} \in R^{n}

be a feasible point of ( 31 ). We say that

h \in R^{n}

h \neq 0

, is a feasible direction at

\bar{x}

x = \bar{x} + ε h

is a feasible point of ( 31 ) for all sufficiently small

ε > 0

. Following [26,Definition2.1] , we say that the second-order sufficient condition holds at

\bar{x}

with modulus

μ > 0

if for all feasible directions

h \in R^{n}

\bar{x}

with

h^{T} (b + H \bar{x}) = h^{T} \nabla_{x} f (\bar{x}) = 0

one has

\begin{matrix} h^{T} H h = h^{T} (\nabla_{x}^{2} ℒ^{(q)} (\bar{x}, \bar{Y})) h \geq μ ∥ h ∥^{2} . \end{matrix}

(36)

After these preliminaries, our main result of this section can be stated as follows.

Theorem 1 Assume that problem 31

()

satisfies Slater's condition. Let

(\bar{x}, \bar{Y}, \bar{S})

be a locally unique and strictly complementary stationary point of 31

()

with data 32

()

D

, and assume that the second-order sufficient condition holds at

\bar{x}

with modulus

μ > 0

. Then, for all sufficiently small perturbations 33

()

Δ D

, there exists a locally unique stationary point

(\bar{x} (Δ D), \bar{Y} (Δ D), \bar{S} (Δ D))

of the perturbed program 31

()

with data

D + Δ D

. Moreover, the point

(\bar{x} (Δ D), \bar{Y} (Δ D), \bar{S} (Δ D))

is a differentiable function of the perturbation 33

()

, and for

Δ D = 0

(\bar{x} (0), \bar{Y} (0), \bar{S} (0)) = (\bar{x}, \bar{Y}, \bar{S})

. The derivative

D_{D} (\bar{x} (0), \bar{Y} (0), \bar{S} (0))

(\bar{x}, \bar{Y}, \bar{S})

is characterized by the directional derivatives

(\dot{x}, \dot{Y}, \dot{S}) : = D_{D} (\bar{x} (0), \bar{Y} (0), \bar{S} (0)) [Δ D]

for any

Δ D

. Here

(\dot{x}, \dot{Y}, \dot{S})

is the unique solution of the system of linear equations,

\begin{matrix} \begin{matrix} A (\dot{x}) + \dot{S} & = - Δ C - Δ A (\bar{x}), \\ H \dot{x} + A^{*} (\dot{Y}) & = - Δ b - Δ H \bar{x} - Δ A^{*} (\bar{Y}), \\ \bar{Y} \dot{S} + \dot{Y} \bar{S} + \dot{S} \bar{Y} + \bar{S} \dot{Y} & = 0, \end{matrix} \end{matrix}

(37)

for the unknowns

\dot{x} \in R^{n}

\dot{Y}, \dot{S} \in S^{m}

. Finally, the second-order sufficient condition holds at

\bar{x} (Δ D)

whenever

Δ D

is sufficiently small.

Remark 1 Theorem 1 is an extension of the sensitivity result for linear semidefinite programs presented in [15] . A related sensitivity result for linear semidefinite programs for a more restricted class of perturbations, but also under weaker assumptions, is given in [29] .

A local Lipschitz continuity property of unique and strictly complementary solutions of linear semidefinite programs is derived in [21] .

Remark 2 While we did not explicitly state a linear independence constraint qualification, commonly referred to as LICQ, it is implied by our condition of uniqueness of the stationary point; see, e.g., [15] . Moreover, our assumptions on the stationary point

(\bar{x}, \bar{Y}, \bar{S})

imply that

\bar{x}

is a strict local minimizer of ( 31 ).

Remark 3 The first and third equations in ( 37 ) are symmetric

m \times m

matrix equations, and so only their upper triangular parts have to be considered. Thus the total number of scalar equations in ( 37 ) is

m^{2} + m + n

. On the other hand, there are

m^{2} + m + n

unknowns, namely the entries of

\dot{x} \in R e^{n}

and of the upper triangular parts of

\dot{Y}, \dot{S} \in S^{m}

. Hence, ( 37 ) is a square system.

Remark 4 In view of part b) of Lemma 1 , the last equation of ( 37 ) is equivalent to

\begin{matrix} \bar{Y} \dot{S} + \dot{Y} \bar{S} = 0 . \end{matrix}

(38)

Thus, Theorem 1 can be stated equivalently with ( 38 ) in ( 37 ). However, the resulting system of equations ( 37 ) would then be overdetermined.

Proof of Theorem 1 . The proof is divided into four steps.

Step 1. In this step, we establish the following result. If the perturbed program has a local solution that is a differentiable function of the perturbation, then the derivative is indeed a solution of ( 37 ).

Slater's condition is invariant under small perturbations of the problem data. Hence, if there exists a local solution

\bar{x} + Δ x

\bar{S} + Δ S

of the perturbed problem near

\bar{x}

\bar{S}

, then the necessary first-order conditions of the perturbed problem apply at

\bar{x} + Δ x

\bar{S} + Δ S

, and state that there exists a matrix

Δ Y

such that

\bar{Y} + Δ Y ≽ 0

\bar{S} + Δ S ≽ 0

, and

\begin{matrix} \begin{matrix} (A + Δ A) (\bar{x} + Δ x) + C + Δ C + \bar{S} + Δ S & = 0, \\ b + Δ b + (H + Δ H) (\bar{x} + Δ x) + (A^{*} + Δ A^{*}) (\bar{Y} + Δ Y) & = 0, \\ (\bar{Y} + Δ Y) (\bar{S} + Δ S) + (\bar{S} + Δ S) (\bar{Y} + Δ Y) & = 0 . \end{matrix} \end{matrix}

(39)

Subtracting from these equations the first three equations of ( 34 ) yields

\begin{matrix} \begin{matrix} (A + Δ A) (Δ x) + Δ S & = - Δ C - Δ A (\bar{x}), \\ (H + Δ H) Δ x + (A^{*} + Δ A^{*}) (Δ Y) & = - Δ b - Δ H \bar{x} - Δ A^{*} (\bar{Y}), \\ \bar{Y} Δ S + Δ Y \bar{S} + Δ S \bar{Y} + \bar{S} Δ Y & = - Δ Y Δ S - Δ S Δ Y . \end{matrix} \end{matrix}

(40)

Neglecting the second-order terms in ( 40 ), and using ( 38 ), we obtain the result claimed in ( 37 ). It still remains to verify the existence and differentiability of

Δ x

Δ Y

Δ S

Step 2. In this step, we prove that the system of linear equations ( 37 ) has a unique solution.

To this end, we show that the homogeneous version of ( 37 ), i.e., the system

\begin{matrix} \begin{matrix} A (\dot{x}) + \dot{S} & = 0, \\ H \dot{x} + A^{*} (\dot{Y}) & = 0, \\ \bar{Y} \dot{S} + \dot{Y} \bar{S} + \dot{S} \bar{Y} + \bar{S} \dot{Y} & = 0, \end{matrix} \end{matrix}

(41)

only has the trivial solution

\dot{x} = 0

\dot{Y} = \dot{S} = 0

Let

\dot{x} \in R^{n}

\dot{Y}, \dot{S} \in S^{m}

be any solution of ( 41 ). Recall that, in view of part b) by Lemma 1 we may assume that

\bar{Y}

and

\bar{S}

are given in diagonal form:

\begin{matrix} \bar{Y} = [\begin{matrix} {\bar{Y}}_{1} & 0 & 0 & 0 \end{matrix}], \bar{S} = [\begin{matrix} 0 & 0 & 0 & {\bar{S}}_{2} \end{matrix}], where {\bar{Y}}_{1}, {\bar{S}}_{2} ≻ 0 and {\bar{Y}}_{1}, {\bar{S}}_{2} are diagonal . \end{matrix}

(42)

Indeed, this can be done by replacing, in ( 41 ),

\bar{Y}

\bar{S}

\dot{Y}

\dot{S}

A (x)

U^{T} \bar{Y} U

U^{T} \bar{S} U

U^{T} \dot{Y} U

U^{T} \dot{S} U

U^{T} A (x) U

, respectively, where

U

is the matrix in ( 18 ), and then multiplying the first and third rows from the left by

U

and from the right by

U^{T}

. Furthermore, in view of part b) by Lemma 1 , any matrices

\dot{Y} \dot{S} \in S^{m}

satisfying the last equation of ( 41 ) are then of the form

\begin{matrix} \dot{Y} = [\begin{matrix} {\dot{Y}}_{1} & {\dot{Y}}_{3} & {\dot{Y}}_{3}^{T} & 0 \end{matrix}], \dot{S} = [\begin{matrix} 0 & {\dot{S}}_{3} & {\dot{S}}_{3}^{T} & {\dot{S}}_{2} \end{matrix}], where {\dot{Y}}_{3} {\bar{S}}_{2} + {\bar{Y}}_{1} {\dot{S}}_{3} = 0 . \end{matrix}

(43)

Next, we establish the inequality

\begin{matrix} {\dot{x}}^{T} H \dot{x} \geq μ ∥ \dot{x} ∥^{2}, \end{matrix}

(44)

where

μ > 0

is the modulus of the second-order sufficient condition ( 36 ). Assume that

\dot{x} \neq 0

Let

\overset{˘}{x} \in R^{n}

be a Slater point for problem ( 31 ). This guarantees that

\begin{matrix} M = [\begin{matrix} M_{1} & M_{3} & M_{3}^{T} & M_{2} \end{matrix}] : = - (A (\overset{˘}{x}) + C) ≻ 0, \end{matrix}

(45)

where the block partitioning

M

is conforming with ( 43 ). For

η > 0

, set

\begin{matrix} h_{η}^{+} : = \dot{x} + η (\overset{˘}{x} - \bar{x}) and h_{η}^{-} : = - \dot{x} + η (\overset{˘}{x} - \bar{x}) . \end{matrix}

(46)

Since

\dot{x} \neq 0

, there exists an

η_{0} > 0

such that

h_{η}^{\pm} \neq 0

for all

0 < η \leq η_{0}

. Next, we prove that for all such

η

, both vectors

h_{η}^{+}

and

h_{η}^{-}

are feasible directions for ( 31 ) at

\bar{x}

. Let

0 < η \leq η_{0}

be arbitrary, but fixed. We then need to verify that

A (\bar{x} + ε h_{η}^{\pm}) + C ≼ 0

for all sufficiently small

ε > 0

. Recall that

A

is a linear function. Using ( 46 ), ( 42 ), ( 43 ), ( 45 ), the first equation of ( 34 ), and the first equation of ( 41 ), one readily verifies that

\begin{matrix} \begin{matrix} A (\bar{x} + ε h_{η}^{\pm}) + C & = (1 - ε η) (A (\bar{x}) + C) + ε η (A (\overset{˘}{x}) + C) \pm ε A (\dot{x}) \\ = - ([\begin{matrix} 0 & 0 & 0 & {\bar{S}}_{2} \end{matrix}] + ε [\begin{matrix} η M_{1} & η M_{3} \pm {\dot{S}}_{3} & (η M_{3} \pm {\dot{S}}_{3})^{T} & η (M_{2} - {\bar{S}}_{2}) \pm {\dot{S}}_{2} \end{matrix}]) . \end{matrix} \end{matrix}

(47)

Recall that

η > 0

is fixed. Since, by ( 42 ) and ( 45 ),

{\bar{S}}_{2} ≻ 0

and

M_{1} ≻ 0

, a standard Schur-complement argument shows that the matrix on the right-hand side of ( 47 ) is negative definite for all sufficiently small

ε > 0

. Thus the vectors ( 46 ) are feasible directions for ( 31 ) at

\bar{x}

for any

η > 0

. This in turn implies

\begin{matrix} {\dot{x}}^{T} (b + H \bar{x}) = {\dot{x}}^{T} \nabla_{x} f (\bar{x}) = 0 . \end{matrix}

(48)

Indeed, suppose that

{\dot{x}}^{T} \nabla_{x} f (\bar{x}) < 0

. Then, for sufficiently small

η > 0

, the feasible direction

h_{η}^{+}

also satisfies

{(h_{η}^{+})}^{T} \nabla_{x} f (\bar{x}) < 0

, and thus

h_{η}^{+}

is a descent direction for the objective function

f

of ( 31 ) at the point

\bar{x}

. This contradicts the local optimality of

\bar{x}

. Likewise, if

{\dot{x}}^{T} \nabla_{x} f (\bar{x}) > 0

, then, for sufficiently small

η > 0

h_{η}^{-}

is a descent direction for the objective function

f

of ( 31 ) at the point

\bar{x}

, leading to the same contradiction. The second-order sufficient condition ( 36 ) also holds true on the closure of the feasible directions. Since

\dot{x}

is the limit of the feasible directions ( 46 ) for

η \to 0

and

\dot{x}

satisfies ( 48 ), the inequality ( 44 ) follows from ( 36 ).

Next recall from ( 42 ) that

{\bar{Y}}_{1}

and

{\bar{S}}_{2}

are positive definite diagonal matrices. The last relation in ( 43 ) thus implies that corresponding entries of the matrices

{\dot{Y}}_{3}

and

{\dot{S}}_{3}

are either zero or of opposite sign. It follows that

〈 {\dot{Y}}_{3}, {\dot{S}}_{3} 〉 \leq 0

, and equality holds if, and only if,

{\dot{Y}}_{3} = {\dot{S}}_{3} = 0

. Using this inequality, together with the first two relations in ( 43 ), the first two equations of ( 41 ), and ( 44 ), one readily verifies that

0 \geq 2 〈 {\dot{Y}}_{3}, {\dot{S}}_{3} 〉 = 〈 \dot{Y}, \dot{S} 〉 = - 〈 \dot{Y}, A (\dot{x}) 〉 = - 〈 A^{*} (\dot{Y}), \dot{x} 〉 = 〈 H \dot{x}, \dot{x} 〉 = {\dot{x}}^{T} H \dot{x} \geq μ ∥ \dot{x} ∥^{2} .

Since

μ > 0

, this implies

\begin{matrix} \dot{x} = 0 and {\dot{Y}}_{3} = {\dot{S}}_{3} = 0 . \end{matrix}

(49)

By the first row of ( 41 ), it further follows that

\begin{matrix} \dot{S} = - A (\dot{x}) = - A (0) = 0 . \end{matrix}

(50)

Thus it only remains to show that

\dot{Y} = 0

. In view of ( 43 ) and ( 50 ), we have

\begin{matrix} \dot{Y} = [\begin{matrix} {\dot{Y}}_{1} & 0 & 0 & 0 \end{matrix}] . \end{matrix}

(51)

Now suppose that

{\dot{Y}}_{1} \neq 0

. Then, by ( 42 ) and ( 51 ), we have

{\bar{Y}}_{ε} : = \bar{Y} + ε \dot{Y} ≽ 0 and {\bar{Y}}_{ε} \neq \bar{Y}

for all sufficiently small

| ε |

. Moreover, using ( 34 ), ( 41 ), and ( 50 ), one readily verifies that the point

(\bar{x}, {\bar{Y}}_{ε}, \bar{S})

also satisfies ( 34 ) for all sufficiently small

| ε |

. This contradicts the assumption that

(\bar{x}, \bar{Y}, \bar{S})

is a locally unique stationary point. Hence

{\dot{Y}}_{1} = 0

and, by ( 51 ),

\dot{Y} = 0

This concludes the proof that the square system ( 41 ) is nonsingular.

Step 3. In this step, we show that the nonlinear system ( 34 ) has a local solution that depends smoothly on the perturbation

Δ D

. To this end, we apply the implicit function theorem to the system

\begin{matrix} A (x) + C + S = 0, H x + b + A^{*} (Y) = 0, Y S + S Y = 0 . \end{matrix}

(52)

As we have just seen, the linearization of ( 52 ) at the point

(\bar{x}, \bar{Y}, \bar{S})

is nonsingular, and hence ( 52 ) has a differentiable and locally unique solution

(\bar{x} (Δ D), \bar{Y} (Δ D), \bar{S} (Δ D))

. Furthermore, we have

\bar{Y} (Δ D), \bar{S} (Δ D) ≽ 0

. This semidefiniteness follows with a standard continuity argument:

The optimality conditions of the nonlinear SDP coincide with the optimality conditions of the linearized SDP. Under our assumptions, the latter one has a unique optimal solution that depends continuously on small perturbations of the data; see, e.g., [15] . Hence the linearized problem at the data point

D + Δ D

has an optimal solution

(\tilde{x}, \tilde{Y}, \tilde{S})

that satisfies the same optimality conditions as

(\bar{x} (Δ D), \bar{Y} (Δ D), \bar{S} (Δ D))

. The solution of the linearized problem also satisfies

\tilde{Y} ≽ 0

\tilde{S} ≽ 0

. Since

(\bar{x} (Δ D), \bar{Y} (Δ D), \bar{S} (Δ D))

is locally unique, it must coincide with

(\tilde{x}, \tilde{Y}, \tilde{S})

, i.e.,

\bar{Y} (Δ D), \bar{S} (Δ D)

satisfy the semidefiniteness conditions.

Step 4. In this step, we prove that the second-order sufficient condition is satisfied at the perturbed solution. Since feasible directions

h

are defined only up to a positive scalar factor, without loss of generality, one may require that

∥ h ∥ = 1

. For the unperturbed problem ( 31 ), the second-order sufficient condition at

\bar{x}

then states that

h^{T} H h \geq μ

for all

h \in R^{n}

with

A (x + ε h) + C ≼ 0

ε = ε (h) > 0

∥ h ∥ = 1

, and

h^{T} (\bar{b} + H \bar{x}) = 0

. To prove that the second-order sufficient condition is invariant under small perturbations

Δ D

of the problem data

D

, we thus need to show that for some fixed

\tilde{μ} > 0

, we have

\begin{matrix} h^{T} (H + Δ H) h \geq \tilde{μ} \end{matrix}

(53)

for all solutions

h \in R^{n}

of the inequalities

\begin{matrix} (A + Δ A) (\bar{x} (Δ D) + ε h) + C + Δ C ≼ 0, & ε = ε (h) > 0, \end{matrix}

\begin{matrix} ∥ h ∥ = 1, & h^{T} (b + Δ b + (H + Δ H) \bar{x} (Δ D)) = 0 . \end{matrix}

In view of the first two relations in ( 39 ), the above is equivalent to

\begin{matrix} ε (A + Δ A) (h) ≼ \bar{S} (Δ D), ε = ε (h) > 0, ∥ h ∥ = 1, h^{T} (A^{*} + Δ A^{*}) (\bar{Y} (Δ D)) = 0 . \end{matrix}

(54)

It remains to show that the set of solutions

h

of ( 54 ) varies continuously with

Δ D

. Indeed, for any fixed

\tilde{μ}

with

0 < \tilde{μ} < μ

, the second-order condition at

\bar{x}

then readily implies that ( 53 ) is satisfied for all solutions

h

of ( 54 ), provided

∥ Δ D ∥

is sufficiently small. In Step 2, we have shown that both

\bar{S} (Δ D)

and

\bar{Y} (Δ D)

are continuous functions of

Δ D

and that the dimension of the null space of

\bar{S} (Δ D)

is constant, namely equal to

k

, for all sufficiently small

∥ Δ D ∥

. Moreover, the null space of

\bar{S} (Δ D)

varies continuously with

Δ D

Let

Δ D_{k}

be a sequence of perturbations with

Δ D_{k} \to 0

. Let

h_{k}

be a sequence of associated solutions of ( 54 ). It suffices to show that any accumulation point

\bar{h}

of the sequence

h_{k}

satisfies ( 54 ) for

Δ D = 0

and the associated values

Δ A = 0

\bar{Y} (0) = \bar{Y}

\bar{S} (0) = \bar{S}

. Since

\bar{Y} (Δ D)

and

Δ A^{*}

vary continuously with

Δ D

, it follows that

\bar{h}

satisfies the last two relations of ( 54 ) for

Δ D = 0

We now assume by contradiction that

ε A (\bar{h}) ⋠ \bar{S}

for any

ε > 0

. Since

\bar{S} ≽ 0

this implies that there exists a vector

z \in R^{m}

with

∥ z ∥ = 1

z^{T} A (\bar{h}) z = \tilde{ε} > 0

, and

z^{T} \bar{S} z = 0

. It follows that

z^{T} (A + Δ A_{k}) (h_{k}) z \geq \frac{\tilde{ε}}{2}

k

is sufficiently large. Since the null space of

\bar{S} (Δ D)

varies continuously with

Δ D

, we have

(z + Δ z_{k})^{T} \bar{S} (Δ D_{k}) (z + Δ z_{k}) = 0

for some small

Δ z_{k} \in R^{m}

whenever

∥ Δ D_{k} ∥

is sufficiently small. We now choose

∥ Δ D_{k} ∥

so small, i.e.,

k

so large, that

(z + Δ z_{k})^{T} (A + Δ A_{k}) (h_{k}) (z + Δ z_{k}) \geq \frac{\tilde{ε}}{4} .

This implies that

h_{k}

does not satisfy ( 54 ), and thus yields the desired contradiction. Hence

\bar{h}

satisfies ( 54 ) for

Δ D = 0

Theorem 1 can be sharpened slightly.

Corollary 1 Under the assumptions of Theorem 1 there exists a small neighborhood

N

of zero in the data space of 31

()

such that for all perturbations

Δ D \in N

of the problem data 32

()

D

, of 31

()

, there exists a local solution

(x_{Δ}, Y_{Δ}, S_{Δ})

of 39

()

near

(\bar{x}, \bar{Y}, \bar{S})

, at which the assumptions of Theorem 2 are also satisfied. Moreover, the second derivatives

\nabla_{D}^{2} (x_{Δ}, Y_{Δ}, S_{Δ}) [Δ D]

of such local solutions

(x_{Δ}, Y_{Δ}, S_{Δ})

are uniformly bounded for all

Δ D \in N

Proof. The first part of the corollary is an immediate consequence of Theorem 1 . For the second part observe that the second derivative is obtained by differentiating the system ( 37 ). For sufficiently small perturbations $Δ D$ , the singular values of this system are uniformly bounded away from zero, and hence the second derivatives are uniformly bounded.

Theorem 1 can be generalized to the class of NLSDPs of the form ( 21 ). Recall that, by ( 22 ) and ( 24 ), the Lagrangian of ( 21 ) and its Hessian are given by

\begin{matrix} ℒ (x, Y) = b^{T} x + ℬ (x) ∙ Y and H (x, Y) = \nabla_{x}^{2} (ℬ (x) ∙ Y), \end{matrix}

(55)

respectively. The generalization of Theorem 1 to problems ( 21 ) is then as follows.

Theorem 2 Let

x^{*}

be a local solution of 21

()

, and let

Y^{*}

be an associated Lagrange multiplier. Assume that the Robinson constraint qualification is satisfied at

x^{*}

and that the point

(x^{*}, Y^{*})

is strictly complementary and locally unique. Finally, assume that the second-order sufficient condition holds at

x^{*}

with modulus

μ > 0

. Then, 21

()

has a locally unique solution for small perturbations of the data

(ℬ, b)

, and the solution depends smoothly on the perturbations.

Proof. First, we define the linear function

A : = D_{x} ℬ (x^{*}) : R^{n} \to S^{m}

, and the matrices

C : = ℬ (x^{*})

and

H : = H (x^{*}, Y^{*})

. Then, the SSP approximation ( 29 ) (with

p = q = 0

) of ( 21 ) at the point

(x^{*}, Y^{*})

is just the quadratic semidefinite problem

\begin{matrix} \begin{matrix} minimize b^{T} Δ x + \frac{1}{2} (Δ x)^{T} H Δ x subject to & Δ x \in R^{n}, \\ A (Δ x) + C ≼ 0 . \end{matrix} \end{matrix}

(56)

Note that ( 56 ) is a problem of the form ( 31 ) with data ( 32 ),

D

. Moreover,

Δ \bar{x} : = 0

\bar{Y} : = Y^{*}

, and

\bar{S} : = - A (0) - C

satisfy the conditions ( 34 ). These conditions coincide with the first-order conditions of ( 21 ), and thus the point

(Δ \bar{x}, \bar{Y}, \bar{S})

is also the unique solution of ( 34 ).

Furthermore, the second-order sufficient condition for ( 56 ) and ( 21 ) coincide. This condition guarantees that

Δ \bar{x}

is a locally unique solution of ( 56 ). Finally, the Robinson constraint qualification for problem ( 21 ) at

x^{*}

implies that problem ( 56 ) satisfies Slater's condition. In particular, all assumptions of Theorem 1 are satisfied. Small perturbations

Δ D

of the data of ( 21 ) result in small changes of the corresponding SSP problem ( 56 ). Since Theorem 1 allows for arbitrary changes in all of the data of ( 56 ), the claims follow.

8 Convergence of the SSP method

In this section, we prove that the plain SSP method with step size 1 is locally quadratically convergent.

For pairs

(x, Y)

, where

x \in R^{n}

Y \in S^{m}

, we use the norm

∥ (x, Y) ∥ : = {(∥ x ∥^{2} + ∥ Y ∥^{2})}^{\frac{1}{2}} .

The main result of this section can then be stated as follows.

Theorem 3 Assume that the function

ℬ

in 21

()

C^{3}

-differentiable and that problem 21

()

has a locally unique and strictly complementary solution

(\bar{x}, \bar{Y})

that satisfies the Robinson constraint qualification and the second-order sufficient condition with modulus

μ > 0

, cf. 36

()

Let some iterate

(x^{k}, Y^{k})

be given and let the next iterate

(x^{k + 1}, Y^{k + 1}) : = (x^{k}, Y^{k}) + (Δ x, Δ Y)

be defined as the local solution of 27

()

, or, equivalently, 29

()

, that is closest to

(x^{k}, Y^{k})

. Then there exist

ε > 0

and

γ < 1 / ε

such that

∥ (x^{k + 1}, Y^{k + 1}) - (\bar{x}, \bar{Y}) ∥ \leq γ {∥ (x^{k}, Y^{k}) - (\bar{x}, \bar{Y}) ∥}^{2}

whenever

∥ (x^{k}, Y^{k}) - (\bar{x}, \bar{Y}) ∥ < ε

Proof. The proof is divided into three steps. In the first step, we establish the exact relation of problems ( 29 ) and ( 31 ). In a second step, we consider a point

x^{k}

near

\bar{x}

. We show that

x^{k}

is the optimal solution of an SSP subproblem the data of which is at most

O (∥ x^{k} - \bar{x} ∥)

away from the data of the SSP subproblem at

(\bar{x}, \bar{Y})

. We remark that

x^{k}

is always the optimal solution of the

(k - 1)

-st subproblem, but the data of this subproblem lies

O (∥ x^{k - 1} - \bar{x} ∥)

away from the SSP subproblem at

(\bar{x}, \bar{Y})

. In a third step, we then show by a perturbation analysis that the correction

Δ x = x^{k + 1} - x^{k}

is of size

O (∥ x^{k} - \bar{x} ∥ + ∥ Y^{k} - \bar{Y} ∥)

and that the residual for the SSP subproblem in the

(k + 1)

-st step is of size

O ((∥ x^{k} - \bar{x} ∥ + ∥ Y^{k} - \bar{Y} ∥)^{2})

Step 1. We first show how the SSP subproblem ( 29 ) can be written in the form ( 31 ). To this end, we define the linear function

A : = D_{x} ℬ (x^{k}) : R^{n} \to S^{m}

, and the matrices

C : = ℬ (x^{k})

and

H : = H (x^{k}, Y^{k})

. Note that the linear constraint

A (Δ x) + C ≼ 0

is just the linearization of the nonlinear constraint

ℬ (x^{k} + Δ x) ≼ 0

about the point

x^{k}

. Finally, let

b

be as in ( 21 ). The SSP subproblem ( 29 ) then takes the simple form ( 56 ), and in particular, it conforms with the format ( 31 ) of Theorem 2 .

Step 2. Let any point

x^{k}

close to

\bar{x}

be given. We show that

Δ x = 0

is a local solution of a problem of the form ( 56 ), where the data is `close' to the data of ( 29 ) at

(\bar{x}, \bar{Y})

. Let

Δ C : = ℬ (\bar{x}) - ℬ (x^{k})

. By continuity of

ℬ

∥ Δ C ∥

is of order

O (∥ x^{k} - \bar{x} ∥)

. Let

Δ b : = - \nabla_{x} (ℬ (x^{k}) ∙ \bar{Y}) - b = - A^{*} (\bar{Y}) - b .

From ( 23 ) and the second row of ( 25 ), we obtain the estimate

∥ Δ b ∥ = O (∥ x^{k} - \bar{x} ∥)

, and the point

(0, \bar{Y}, \bar{S})

satisfies the first-order conditions,

A (0) + C + Δ C + \bar{S} = 0, b + Δ b + H \cdot 0 + A^{*} (\bar{Y}) = 0, \bar{Y} \bar{S} = 0,

for the quadratic semidefinite program

\begin{matrix} \begin{matrix} minimize (b + Δ b)^{T} Δ x + \frac{1}{2} (Δ x)^{T} H Δ x subject to & Δ x \in R^{n}, \\ A (Δ x) + C + Δ C ≼ 0 . \end{matrix} \end{matrix}

(57)

Let

\begin{matrix} \bar{A} : = D_{x} ℬ (\bar{x}), \bar{b} : = b, \bar{C} : = ℬ (\bar{x}), \bar{H} : = \nabla_{x}^{2} ℒ (\bar{x}, \bar{Y}), \end{matrix}

(58)

be the data of the SSP subproblem ( 29 ) at the point

(\bar{x}, \bar{Y})

. Then, the data of ( 57 ) differs from the data ( 58 ) by a perturbation of norm

O (∥ x^{k} - \bar{x} ∥ + ∥ Y^{k} - \bar{Y} ∥)

. Here, the term

∥ Y^{k} - \bar{Y} ∥

reflects the fact that

H

and

\bar{H}

also differ by the choice of

Y

. Note that the point

(Δ x, Y, S) = (0, \bar{Y}, \bar{S})

satisfies the first-order optimality conditions,

\begin{matrix} \bar{A} (Δ x) + \bar{C} + S = 0, b + \bar{H} Δ x + {\bar{A}}^{*} (Y) = 0, Y S = 0, \end{matrix}

(59)

of the quadratic problem ( 56 ) with data ( 58 ). As shown in Theorem 2 , the assumptions for the nonlinear SDP ( 21 ) at

(\bar{x}, \bar{Y})

imply that problem ( 56 ) with data ( 58 ) satisfies all conditions of Theorem 2 at

(0, \bar{Y}, \bar{S})

Since the second-order sufficient condition depends continuously on the data of ( 31 ), it follows that for ( 57 ), the second-order condition at

(0, \bar{Y}, \bar{S})

is satisfied, provided that

∥ x^{k} - \bar{x} ∥

and

∥ Y^{k} - \bar{Y} ∥

are sufficiently small. Thus problem ( 57 ) fulfills all assumptions of Theorem 2 .

Step 3. By definition,

(Δ x, Y) = (0, \bar{Y})

is the optimal solution (with associated multiplier) of ( 57 ). Let

(x^{k}, Y^{k})

be close to

(\bar{x}, \bar{Y})

. The SSP subproblem replaces the data

Δ b

and

Δ C

of ( 57 ) by 0 (of the respective dimension). Thus, the data of ( 57 ) is changed by a perturbation of order

O (∥ x^{k} - \bar{x} ∥ + ∥ Y^{k} - \bar{Y} ∥)

. We assume that this perturbation lies in the neighborhood

N

about zero as guaranteed by Corollary 1 . Denote the optimal solution of the SSP subproblem by

(Δ x, \bar{Y} + Δ \hat{Y})

The SSP subproblem is then used to define

(x^{k + 1}, Y^{k + 1})

. Let

A^{+} : = D_{x} ℬ (x^{k + 1}), C^{+} : = ℬ (x^{k + 1}), H^{+} : = \nabla_{x}^{2} ℒ (x^{k + 1}, Y^{k + 1})

be the data of the SSP subproblem at the next,

(k + 1)

-st iteration.

Corollary 1 states that

(Δ x, Δ \hat{Y})

are given by the tangent equations ( 37 ) plus some uniformly bounded second-order terms. Thus,

(Δ x, Δ \hat{Y})

are of the order

O (∥ x^{k} - \bar{x} ∥ + ∥ Y^{k} - \bar{Y} ∥)

. Here,

Δ \hat{Y}

is a correction of the unknown point

\bar{Y}

, while the correction

Δ Y = Y^{k + 1} - Y^{k}

produced by the SSP subproblem has the form

Δ Y = Δ \hat{Y} + Y^{k} - \bar{Y}

. Obviously, also the norm

∥ Δ Y ∥

of this correction is of the order

O (∥ x^{k} - \bar{x} ∥ + ∥ Y^{k} - \bar{Y} ∥)

Next, we compute an upper bound on the size of the residuals of the first and second equations in ( 59 ) at

(x^{k + 1}, Y^{k + 1}, S^{k + 1})

. Note that the residual term of the third equation in ( 59 ) is zero. By definition of

(Δ x, Y^{k + 1}, S^{k + 1})

, it follows that

\begin{matrix} A (Δ x) + C + S^{k + 1} = 0, b + H Δ x + A^{*} (Y^{k + 1}) = 0, Y^{k + 1} S^{k + 1} = 0 . \end{matrix}

(60)

If the data of ( 21 ) is

C^{3}

-smooth, this implies that

\begin{matrix} (A^{+})^{*} (Y^{k + 1}) + b & = A^{*} (Y^{k + 1}) + b + Δ A^{*} (Y^{k + 1}) \end{matrix}

\begin{matrix} = - H Δ x + Δ A^{*} (Y^{k}) + Δ A^{*} (Δ Y) \end{matrix}

\begin{matrix} = - \nabla_{x}^{2} (ℬ (x^{k}) ∙ Y^{k}) Δ x + (\nabla_{x} ℬ (x^{k} + Δ x) - \nabla_{x} ℬ (x^{k})) ∙ Y^{k} + Δ A^{*} (Δ Y) \end{matrix}

\begin{matrix} = O (∥ Δ x ∥^{2} + ∥ Δ Y ∥^{2}), \end{matrix}

where

Δ A : = A^{+} - A

. Likewise, it follows from ( 60 ) that

\begin{matrix} C^{+} + S^{k + 1} & = Δ C + C + S^{k + 1} = Δ C - A (Δ x) \end{matrix}

\begin{matrix} = ℬ (x^{k + 1}) - ℬ (x^{k}) - D_{x} ℬ (x^{k}) [Δ x] = O (∥ Δ x ∥^{2}), \end{matrix}

where

Δ C : = C^{+} - C

. Hence, we can define perturbations

\hat{b}

and

{\hat{C}}^{+}

b

and

C^{+}

of order

∥ \hat{b} - b ∥ + ∥ {\hat{C}}^{+} - C^{+} ∥ = O ((∥ x^{k} - \bar{x} ∥ + ∥ Y^{k} - \bar{Y} ∥)^{2})

such that

(Δ x, Y, S) = (0, Y^{k + 1}, S^{k + 1})

is an optimal solution of the problem ( 56 ) with data

A^{+}

\hat{b}

{\hat{C}}^{+}

H^{+}

. By the same derivation as above, the next SSP step has length

O ((∥ x^{k} - \bar{x} ∥ + ∥ Y^{k} - \bar{Y} ∥)^{2})

and generates residuals of order

O ((∥ x^{k} - \bar{x} ∥ + ∥ Y^{k} - \bar{Y} ∥)^{4})

. Repeating this process, it then follows by a standard argument that

∥ x^{k + 1} - \bar{x} ∥

and

∥ Y^{k + 1} - \bar{Y} ∥

are of order

O ((∥ x^{k} - \bar{x} ∥ + ∥ Y^{k} - \bar{Y} ∥)^{2})

as well.

Remark 5 As mentioned before, one will typically choose to solve SSP subproblems with a positive semidefinite approximation

\hat{H}

to the Hessian of the Lagrangian. A proof of convergence for such modifications is the subject of current research; see, e.g., [10] . Since all the data enters in a continuous fashion in the preceding analysis, it follows that the SSP method with step size one is still locally superlinearly convergent if the matrices

H^{k}

in ( 28 ) are replaced by approximations

{\hat{H}}^{k}

with

∥ H^{k} - {\hat{H}}^{k} ∥ \to 0

Remark 6 The assumption of

C^{3}

-differentiability of the function

ℬ

in Theorem 3 can be weakened to

C^{2}

-differentiability and a Lipschitz condition for the Hessian at

\bar{x}

The result of Theorem 3 can be extended to the following slightly more general class of NLSDPs. Given a vector

b \in R^{n}

, a matrix-valued function

ℬ : R^{n} \to S^{m}

, and two vector-valued functions

c : R^{n} \to R^{p}

and

d : R^{n} \to R^{q}

, we consider problems of the following form:

\begin{matrix} \begin{matrix} minimize b^{T} x subject to & x \in R^{n}, \\ ℬ (x) ≼ 0, \\ c (x) \leq 0, \\ d (x) = 0 . \end{matrix} \end{matrix}

(61)

The Lagrangian of problem ( 21 ) takes the form

ℒ : R^{n} \times S^{m} \times R^{p} \times R^{q} \to R

\begin{matrix} ℒ (x, Y, u, v) : = b^{T} x + ℬ (x) ∙ Y + u^{T} c (x) + v^{T} d (x) . \end{matrix}

(62)

Its gradient with respect to

x

is given by

\begin{matrix} g (x, Y, u, v) : = \nabla_{x} ℒ (x, Y, u, v) = b + \nabla_{x} (ℬ (x) ∙ Y) + \nabla_{x} c (x) u + \nabla_{x} d (x) v \end{matrix}

(63)

and its Hessian by

\begin{matrix} H (x, Y, u, v) : = \nabla_{x}^{2} ℒ (x, Y, u, v) = \nabla_{x}^{2} (ℬ (x) ∙ Y) + \sum_{i = 1}^{p} u_{i} \nabla_{x}^{2} c_{i} (x) + \sum_{j = 1}^{q} v_{j} \nabla_{x}^{2} d_{j} (x) . \end{matrix}

(64)

Note that in ( 63 ), the gradients of the vector-valued functions

c (x)

and

d (x)

are defined as

\nabla_{x} c (x) : = (D_{x} c (x))^{T}

and

\nabla_{x} d (x) : = (D_{x} d (x))^{T}

For NLSDPs ( 61 ), the SSP subproblems are of the form

\begin{matrix} \begin{matrix} minimize b^{T} Δ x + \frac{1}{2} (Δ x)^{T} H^{k} Δ x subject to & Δ x \in R^{n}, \\ ℬ (x^{k}) + D_{x} ℬ (x^{k}) [Δ x] ≼ 0, \\ c (x^{k}) + D_{x} c (x^{k}) Δ x \leq 0, \\ d (x^{k}) + D_{x} d (x^{k}) Δ x = 0 . \end{matrix} \end{matrix}

(65)

The extension of Theorem 3 to problems ( 61 ) is as follows.

Theorem 4 Assume that the functions

ℬ

c

, and

d

in ( 61 ) are

C^{3}

-differentiable, and that problem 61

()

has a locally unique and strictly complementary solution

(\bar{x}, \bar{Y}, \bar{u}, \bar{v})

that satisfies the Robinson constraint qualification and the second-order sufficient condition with modulus

μ > 0

, cf. 36

()

. Let some iterate

(x^{k}, Y^{k}, u^{k}, v^{k})

be given and let the next iterate

(x^{k + 1}, Y^{k + 1}, u^{k + 1}, v^{k + 1}) : = (x^{k}, Y^{k}, u^{k}, v^{k}) + (Δ x, Δ Y, Δ u, Δ v)

be defined as the local solution of 65

()

that is closest to

(x^{k}, Y^{k}, u^{k}, v^{k})

. Then there exist

ε > 0

and

γ < 1 / ε

such that

∥ (x^{k + 1}, Y^{k + 1}, u^{k + 1}, v^{k + 1}) - (\bar{x}, \bar{Y}, \bar{u}, \bar{v}) ∥ \leq γ {∥ (x^{k}, Y^{k}, u^{k}, v^{k}) - (\bar{x}, \bar{Y}, \bar{u}, \bar{v}) ∥}^{2}

whenever

∥ (x^{k}, Y^{k}, u^{k}, v^{k}) - (\bar{x}, \bar{Y}, \bar{u}, \bar{v}) ∥ < ε

Proof. By our assumption on strict complementarity, all entries of the vector $\bar{v}$ of the Lagrange multipliers associated with the equality constraints $d (x) = 0$ of ( 21 ) are different from zero.
Without loss of generality, we assume that $\bar{v} > 0$ . Indeed, for any entry ${\bar{v}}_{j} < 0$ we replace the corresponding constraint $d_{j} (x) = 0$ by the equivalent constraint $- d_{j} (x) = 0$ . These sign changes do not change the iterates generated by ( 29 ). Moreover, for $(x^{k}, Y^{k}, u^{k}, v^{k})$ sufficiently close to $(\bar{x}, \bar{Y}, \bar{u}, \bar{v})$ it follows from $\bar{v} > 0$ that the iterates do not change when the constraints $d (x) = 0$ are replaced by $d (x) \leq 0$ . We can thus assume that $q = 0$ , i.e., there are no equality constraints in ( 61 ).
We further assume that, without loss of generality, the matrix $ℬ$ is augmented to a $2 \times 2$ block diagonal matrix, where the $(2, 2)$ -block is the diagonal matrix $D i a g (c (x))$ . Thus, for the analysis of the SSP method we may assume that $p = q = 0$ in ( 21 ), i.e., we only need to consider problems of the form ( 21 ).

9 Numerical results

In this section, we present results of some numerical experiments with a Matlab implementation of the SSP method. Actually, our Matlab program is for a slightly more general class of nonlinear programs with conic constraints (NLCPs). The numerical experiments with our program illustrate the theoretical results of the preceding sections. In particular, quadratic convergence is observed for problems where the Hessian

H

of the Lagrangian at the optimal solution is positive semidefinite. In cases where

H

is not positive semidefinite, our implementation uses perturbations of the nonconvex SSP subproblems in order to obtain convex conic subproblems.

In these cases, typically, the rate of convergence of the algorithm based on such perturbed problems is only linear.

The Matlab program generates its search directions by solving conic quadratic subproblems using Version 1.05R5 of SeDuMi [28] . SeDuMi allows free and positive variables as well as Lorentz-cone (“ice-cream cone”) constraints, rotated Lorentz-cone constraints, and semidefinite cone constraints. The NLCPs can also be formulated in terms of these cones. In order to simplify the use of SeDuMi for the SSP subproblems, the NLCPs are rewritten in the following standard format:

\begin{matrix} \begin{matrix} minimize c^{T} x subject to & x \in K, \\ F (x) = 0 . \end{matrix} \end{matrix}

(66)

Here,

K

is a Cartesian product of free variables and several cones of the types allowed in SeDuMi. We tested the following techniques for generating positive semidefinite approximations of

H

: a BFGS approach, the Hessian of the augmented Lagrangian, and the orthogonal projection of

H

onto the cone of positive semidefinite matrices. Our experiences with these techniques are as follows.

1. The BFGS approach can result in considerably more SSP iterations compared to the projection of the Hessian of the Lagrangian. Moreover, the BFGS approach strongly depends on the initial matrix $H_{0}$ . A good choice is $H_{0} : = V max (D, ε I) V^{T}$ , where $H = V D V^{T}$ is the eigenvalue decomposition of $H$ .
2. The use of the Hessian of the augmented Lagrangian can be a good choice for some problems, but for most of our test problems the penalty parameter had to be very large to obtain a semidefinite Hessian. This, in turn, significantly reduced the precision of our computations.
3. In spite of not being affinely invariant, the use of the projection of the Hessian of the Lagrangian resulted in the most efficient overall algorithm.

We also tested different step length strategies.

1. The following penalty line search with a quadratic correction gave good results for all test cases. The SSP subproblem provides a search direction $Δ x$ for problem ( 66 ). By solving a least-squares problem, a vector $q$ is computed satisfying $D_{x} F (x) q = - F (x + Δ x)$ . For $λ \in [0, 1]$ , a line search along the points $x (λ) : = x + λ Δ x + λ^{2} q$ is performed based on the penalty function $M ∥ F (x (λ)) ∥ + c^{T} x (λ),$ where $M > 0$ is a penalty parameter.
2. For some examples, the choice of a filter approach was slightly better. In the filter approach used here, a Euclidean trust-region radius was always set to be 1.5 times larger than the previous step, and non-acceptable steps were not discarded, but instead an Armijo-type step-length reduction was used to generate an acceptable step. The motivation for this modified filter strategy lies in the fact that the computation of a solution of a subproblem is very expensive, and therefore discarding the solution of a subproblem is avoided. The above filter approach led to very fast convergence, especially for convex problems.
3. For the examples presented here, the trust-region approach was the best choice. The SSP subproblem was restricted by an additional Euclidean trust-region constraint. For problems of the form ( 67 ) below, it was sufficient to apply the trust-region constraint only to the variables $X_{G}, X_{C}$ , while $P, S$ remain free. For these examples, an additional corrector step significantly accelerated the convergence. For this corrector step, $X_{G}, X_{C}$ is kept fixed, and $P, S$ are updated by solving an additional linear SDP. At each iteration, the ratio between predicted and actual reduction was computed. Depending on that ratio, the step was accepted and the trust region was increased or decreased, or the step was rejected and the trust region was decreased.

For our numerical examples, we use nonlinear nonconvex SDPs of the form ( 10 ), which we rewrite in the form ( 67 ) below. Recall that in ( 10 ),

G, C \in R^{n \times n}

and

B_{1}, B_{2} \in R^{n \times m}

are given data matrices. The nonconvex NLSDP used for the numerical examples is then as follows:

\begin{matrix} \begin{matrix} minimize ∥ S ∥ subject to & P \in R^{n \times n}, S \in R^{n \times m}, \\ X_{G} \in R^{n \times n}, ∥ X_{G} ∥ \leq r_{G}, \\ X_{C} \in R^{n \times n}, ∥ X_{C} ∥ \leq r_{C}, \\ P^{T} B_{1} + S = B_{2}, \\ P^{T} (G + X_{G}) + (G + X_{G})^{T} P ≽ ɛ_{G} I, \\ P^{T} (C + X_{C}) + (C + X_{C})^{T} P ≽ ɛ_{C} I, \\ P^{T} (C + X_{C}) - (C + X_{C})^{T} P = 0 . \end{matrix} \end{matrix}

(67)

Furthermore, in ( 67 ), in addition to the constraints on the norms of the perturbations

X_{G}

and

X_{C}

, we restrict

X_{G}

and

X_{C}

to have possible nonzero entries only in certain positions, which depend on the nonzero structure of the given matrices

G

and

C

, respectively. For our numerical tests, the data matrices

G

and

C

in ( 67 ) are generated as follows. First, two matrices

C_{o r g}

and

G_{o r g}

were constructed such that the associated transfer function is guaranteed to be positive real. Then, certain entries of

G_{o r g}

and

C_{o r g}

were perturbed by random perturbations of norm at most

ɛ_{G}

and

ɛ_{C}

respectively, to define the data matrices

G

and

C

. In all our examples, the transfer functions of the systems given by the resulting matrices

C

and

G

were not positive real.

All our computations were run on a Xeon with a clock rate of 2.8 GHz and 3 GB RAM. All solutions were computed to a precision of 12 decimal digits. In the following table, we list the problem dimension

n

, the total number

M (n)

of equality constraints, the total number

N (n)

of scalar unknowns, the number of iterations, and the cpu time (in seconds).

$n$	$M (n)$	$N (n)$	iterations	cpu time
8	118	285	5	3.71
9	146	348	5	4.43
10	177	417	7	8.06
11	211	492	5	7.18
12	248	573	8	16.05
13	288	660	4	10.88
14	331	753	6	20.77
15	377	852	7	30.12
16	426	957	6	34.38
17	478	1068	5	37.40
18	533	1185	10	91.17
19	591	1308	4	47.61
20	652	1437	5	83.66
21	716	1572	4	289.48
22	783	1713	4	416.31
23	853	1862	5	151.33
24	926	2013	6	683.54
25	1002	2176	3	145.60
26	1081	2337	5	612.22
27	1163	2508	7	518.92
28	1248	2685	5	789.41
29	1336	2868	4	475.52
30	1427	3057	7	4213.50
31	1521	3252	4	784.34
32	1618	3455	6	4659.64
33	1718	3660	5	1130.44
34	1821	3877	2	630.53
35	1927	4092	6	1799.36

Table 9 shows that the number of iterations is nearly independent of the dimension

n

of the problem, while—as expected—the cpu time increases with

n

. The total number of constraints is approximately

M (n) \approx \frac{3}{2} n^{2}

, and the total number of scalar variables is approximately

N (n) \approx 3 n^{2}

. The number of iterations to solve the linear semidefinite subproblems not only depends on the dimension, but also on other properties of the problem as, for example, a comparison of the problems of dimension 32 and 33 shows. In this case, the iteration counts differ only by one, yet the cpu time quadruples, since SeDuMi needs more iterations to solve the subproblems. Some of the linear semidefinite subproblems are nearly infeasible, a situation for which SeDuMi (and other solvers) needs a higher number of interior-point iterations.

10 Concluding remarks

We have discussed the SSP method, which is a generalization of the SQP method for standard nonlinear programs to nonlinear semidefinite programming problems. For the derivation of this generalization, we have chosen a motivation that contrasts the SSP method with primal-dual interior methods. For interior methods that are applied directly to nonlinear semidefinite programs, the choice of the symmetrization procedure is considerably more complicated than in the linear case since the system matrix is no longer positive semidefinite. In the proposed method, the choice of the symmetrization scheme is shifted in a very natural way to the subproblems, and is thus reduced to a well-studied problem. Our convergence analysis differs from the convergence analyses of standard SQP methods in that it is based on a sensitivity result for certain optimal solutions of quadratic semidefinite programs. The derivation of this sensitivity result is also of independent interest.

References

Alizadeh, F. (1995): Interior point methods in semidefinite programming with applications to combinatorial optimization. SIAM J. Optim. 5, 13–51
Alizadeh, F., Haeberly, J.-P.A., Overton, M.L. (1998): Primal-dual interior-point methods for semidefinite programming: convergence rates, stability and numerical results, SIAM J. Optim. 8, 746–768
Anderson, B.D.O., Vongpanitlerd, S. (1973): Network Analysis and Synthesis. Prentice-Hall, Englewood Cliffs, New Jersey
Bai, Z., Freund, R.W. (2000): Eigenvalue-based characterization and test for positive realness of scalar transfer functions, IEEE Trans. Automat. Control 45, 2396–2402
Bai, Z., Freund, R.W. (2001): A partial Padé-via-Lanczos method for reduced-order modeling, Linear Algebra Appl. 332–334, 139–164
Boggs, P.T., Tolle, J.W. (1995): Sequential quadratic programming. Acta Numer. 4, 1–51
Bonnans, J.F., Gilbert, J.C., Lemaréchal, C., Sagastizábal, C. (1997): Optimisation Numérique. Springer-Verlag, Berlin
Boyd, S., El Ghaoui, L., Feron, E., Balakrishnan, V. (1994): Linear Matrix Inequalities in System and Control Theory. SIAM Publications, Philadelphia
Coelho, C.P., Phillips, J.R., Silveira, L.M. (2004): A convex programming approach for generating guaranteed passive approximations to tabulated frequency-data. IEEE Trans. Computer-Aided Design 23, 293–301
Correa, R., Ramírez C., H. (2004): A global algorithm for nonlinear semidefinite programming. SIAM J. Optim. 15, 303–318
Fares, B., Apkarian, P., Noll, D. (2001): An augmented Lagrangian method for a class of LMI-constrained problems in robust control theory. Internat. J. Control 74, 348–360
Fares, B., Noll, D., Apkarian, P. (2002): Robust control via sequential semidefinite programming. SIAM J. Control Optim. 40, 1791–1820
Freund, R.W. (2003): Model reduction methods based on Krylov subspaces. Acta Numer. 12, 267–319
Freund, R.W., Jarre, F. (2004): An extension of the positive real lemma to descriptor systems. Optim. Methods Softw. 19, 69–87
Freund, R.W., Jarre, F. (2004): A sensitivity result for semidefinite programs. Oper. Res. Lett. 32, 126–132
Hörmander, L. (1973): An Introduction to Complex Analysis in Several Variables, Second Edition. North-Holland Publishing Company, Amsterdam
Kocvara, M., Stingl, M. (2003): PENNON — A code for convex nonlinear and semidefinite programming. Optim. Methods Softw. 18, 317–333
Luo, Z.-Q., Sturm, J.F., Zhang, S. (1998): Superlinear convergence of a symmetric primal-dual path following algorithm for semidefinite programming, SIAM J. Optim. 8, 59–81
Mangasarian, O.L. (1969): Nonlinear Programming. MvGraw-Hill, New York
Monteiro, R.D.C., Zhang, Y. (1998): A unified analysis for a class of long-step primal-dual path-following interior-point algorithms for semidefinite programming. Math. Program. 81, 281–299
Nayakkankuppam, M.V., Overton, M.L. (1999): Conditioning of semidefinite programs. Math. Program. 85, 525–540
D. Noll, personal communication, Sept 2004.
Overton, M.L. (1988): On minimizing the maximum eigenvalue of a symmetric matrix. SIAM J. Matrix Anal. Appl. 9, 256–268
Overton, M.L. (1992): Large-scale optimization of eigenvalues. SIAM J. Optim. 2, 88–120
Robinson, S.M. (1976): First order conditions for general nonlinear optimization. SIAM J. Appl. Math. 30, 597–607
Robinson, S.M. (1982): Generalized equations and their solutions, part II: applications to nonlinear programming. Math. Prog. Study 19, 200–221
Shapiro, A., Scheinberg, K. (2000): Duality and optimality conditions. In: Wolkowicz, H., Saigal, R., Vandenberghe, L. (eds.) Handbook of Semidefinite Programming: Theory, Algorithms and Applications. Kluwer Academic Publishers, Boston, pp. 67–110
Sturm, J.F. (1999): Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones. Optim. Methods Softw. 11–12, 625–653
Sturm, J.F., Zhang, S. (2001): On sensitivity of central solutions in semidefinite programming. Math. Program. 90, 205–227
Todd, M.J. (2001): Semidefinite optimization. Acta Numer. 10, 515–560
Tuan, H.D., Apkarian, P., Nakashima, Y. (2000): A new Lagrangian dual global optimization algorithm for solving bilinear matrix inequalities. Internat. J. Robust Nonlinear Control 10, 561–578
Vandenberghe, L., Boyd, S. (1996): Semidefinite programming. SIAM Rev. 38, 49–95
Vanderbei, R.J. (1997): LOQO user's manual—version 3.10. Report SOR 97-08, Princeton University, Princeton, New Jersey
Vanderbei, R.J., Benson, H.Y., Shanno, D.F. (2002): Interior-point methods for nonconvex nonlinear programming: filter methods and merit functions. Comput. Optim. Appl. 23, 257–272
Wolkowicz, H., Saigal, R., Vandenberghe, L., eds. (2000): Handbook of Semidefinite Programming: Theory, Algorithms and Applications. Kluwer Academic Publishers, Boston