Step by step design of an interior-point solver in self-dual conic optimization with applications Jean Charles Gilbert, INRIA December 18, 2019 These notes present a project in numerical optimization dealing with the implementation of an interior-point method for solving a self-dual conic optimization (SDCO) problem. The cone is the Cartesian product of cones of positive semidefinite matrices of various dimensions (imposing to matrices to be positive semidefinite) and of a positive orthant (imposing to a vector to have nonnegative components). Therefore, the solved problem encompasses semidefinite and linear optimization.

The project was given in a course entitled Advanced Continuous Optimization II at the University Paris-Saclay, in 2016-2020. The solver is designed step by step during a series of 5 sessions of 4 hours each. Each session corresponds to a chapter of these notes (or a part of it). The correctness of the SDCO solver is verified during each session on small academic problems, having diverse properties. During the last session, the developed piece of software is used to minimize a univariate polynomial on an interval and to solve a few small size rank relaxations of QCQO (quadratically constrained quadratic optimization) problems, modeling various instances of the OPF (optimal power flow) problem. The student has to master not ony the implementation of the interior-point solver, but is also asked to understand the underliying theory by solving exercises consisting in proving some properties of the implemented algorithms.

The goal of the project is not to design an SDCO solver that would beat the best existing solver but to help the students to understand and demystify what there is inside such a piece of software. As a side outcome, this course also shows that a rather performent SDCO solver can be realized in a relatively short time. In addition, the student will be able to improve the developed piece of software, in case this is required by professionel needs. For a review of SDCO codes and more details on their development, we refer the reader to [START_REF]Handbook of Semidefinite Programming -Theory, Algorithms, and Applications[END_REF][START_REF] Anjos | Handbook on Semidefinite, Conic and Polynomial Optimization[END_REF], in particular [START_REF] Mittelmann | The state-of-the-art in conic optimization software[END_REF][START_REF] Toh | On the implementation and usage of SDPT3a Matlab software package for semidefinite-quadratic-linear programming, version 4[END_REF], which have guided us in the composition of these notes.

In case you use these notes in a paper, thanks to cite them like in [START_REF] Ch | Step by step design of an interior-point solver in self-dual conic optimization with applications[END_REF]. It also shows how to compute a direction of displacement, which will be used at each iteration by the interior point algorithms described in the other sessions.

1.1 The problem to solve

Space structures

We denote by E the Euclidean vector space of the variables defining the primal form of the convex optimization considered in all this project; see section 1.1.2. It is equipped with a scalar product denoted by •, • E . Conic optimization is formulated thanks to a cone of E, which is denoted by K. By definition a cone K of E is a set verifying R ++ K ⊂ K, where R ++ := {t ∈ R : t > 0}, meaning that tx must belong to K as soon as t > 0 and x ∈ K. The dual cone of K is the closed convex cone defined by K + := {y ∈ E : y, x E 0 for all x ∈ K}.

The cone K is said to be self-dual if K + = K. Hence, a self-dual cone is necessarily closed and convex. The unknowns of the optimization problem we consider below are elements of a vector space E that are constrained to belong to some self-dual cone K and to satisfy some affine constraint.

More specifically, the vector space E considered in this project is the Cartesian product

E := S n 1 × • • • × S ns × R n l ,
where, for j ∈ [1 : s], S n j is the space of symmetric real matrices of order n j ; it has dimension n j (n j +1)/2. This means that the unknowns of the optimization problem below are s symmetric matrices of order n 1 , . . . , n s and a vector of variables of dimension n l . An element x ∈ E will be denoted by

x = (x s,1 , . . . , x s,s , x l ), where x s,j ∈ S n j , for j ∈ [1 : s], and x l ∈ R n l (unusually, we have chosen to represent matrices by lower case letters, because in E they are assembled with a vector denoted by a lower case letter). The scalar product on some S n j is defined and denoted by where tr a := i a ii denotes the trace of the matrix a ∈ R n j ×n j . The associated norm is the so-called Frobenius norm; its value at a ∈ S n j is denoted by

a S n j =   n j k,l=1 a 2 kl   1/2 .
The Euclidean scalar product is supposed to equip R n l . Finally the scalar product on E takes at (x, s) ∈ E × E the value

x, s E := s j=1

x s,j , s s,j S n j + (x l ) T (s l ).

The associated norm is

x E :=   s j=1 x s,j 2 S n j + x l 2 2   1/2 .
A product of two symmetric matrices is not necessary symmetric. Since these products will appear below, it is natural to introduce the vector space

F := R n 1 ×n 1 × • • • × R ns×ns × R n l ,
where R n j ×n j denotes the set of square real matrices of order n j . Of course, E ⊂ F. The dimension of F is given by

n F := s j=1 n 2 j + n l . (1.1)
The space F has its own scalar product, which takes at (x, s) ∈ F × F the value

x, s F := s j=1 tr(x s,j ) T (s s,j ) + (x l ) T (s l ). (1.2)
The associated norm is

x F := x, x 1/2
F . Clearly, when x and s ∈ E, x, s F = x, s E . The vector space F appears for example, when one defines the Hadamard product of two vectors x and s ∈ E, which is the vector x s = (x s,1 s s,1 , . . . , x s,s s s,s , x l • s l ) ∈ F, where x s,j s s,j is the product of the matrices x s,j and s s,j and x l • s l is the standard Hadamard product of two vectors, which is

(x l • s l ) i = x l i s l i for all i ∈ [1 : n l ].
The fact that x s may not be in E when x and s ∈ E will be a source of trouble in the sequel. The double Hadamard product x s v is equivalently defined by (x s) v or x (s v). For x = (x s,1 , . . . , x s,s , x l ) ∈ F, one defines x T = ((x s,1 ) T , . . . , (x s,s ) T , x l ) ∈ F.

(1.3)

Note the calculus rules

(a b) T = b T a T , (1.4a) a b, c F = a, c b T F = b, a T c F . (1.4b)
Note also that, when a ∈ F, a + a T ∈ E.

The self-dual cone K of E considered in this project is a Cartesian product of cones, namely

K := S n 1 + × • • • × S ns + × R n l + , where S n j + is the cone of positive semidefinite matrices of S n j and R n l + := {x ∈ R n l : x
0} is the nonnegative orthant of R n l (vector inequalities have to be understood componentwise: x 0 iff x i 0 for all i). The cone K is self-dual since this is the case for S n j + and R n l + . This choice of K implies that the unknown matrices in S n j of the optimization problem are forced to be positive semidefinite and that the unknown vector in R n l is forced to have nonnegative components. We associate with K the operator min K : E → R defined at x ∈ E by min K (x) = min λ min (x s,1 ), . . . , λ min (x s,s ), min

i∈[1 : n l ] x l i , (1.5)
where λ min (x s,j ) is the smallest eigenvalue of the symmetric matrix x s,j . Therefore x ∈ K if and only if min K (x) 0. The strict feasible sets of the optimization problems (P ) and (D) defined below make use of the strict cone

K s := S n 1 ++ × • • • × S ns ++ × R n l ++ ,
where S

n j
++ is the cone of positive definite matrices of S n j and R n l ++ := {x ∈ R n l : x > 0} is the positive orthant (strict inequalities also act componentwise: x > 0 iff x i > 0 for all i). Clearly, x ∈ K s if and only if min K (x) > 0. For x = (x s,1 , . . . , x s,s , x l ) in K s , one can define

x -1 := ((x s,1 ) -1 , . . . , (x s,s ) -1 , (x l ) -1 ) ∈ K s , x -T := (x T ) -1 = (x -1 ) T ∈ K s , x 1/2 := ((x s,1 ) 1/2 , . . . , (x s,s ) 1/2 , (x l ) 1/2 ) ∈ K s ,
where the exponent -1 (resp. 1/2) refers to the inverse (resp. square root) of a positive definite matrix or the componentwise inverse (resp. square root) of a positive vector.

Directions of displacement

A particular element of K s , which is used continually below, is e := (I n 1 , . . . , I ns , e l ), where I n j is the identity matrix of order n j and e l is the vector of all ones in R n l . Observe that x e = e x = x, x s, e E = x, s E , and

n c := e 2 E = s j=1 n j + n l . (1.6)
The number n c is called the complexity module of the problem (hence the index c).

Note the difference with the dimension n F of the space F, defined in (1.1), in which the square of the dimensions n j appears instead. We will see with (2.9) that the number of iterations to reach optimality at a given precision depends on the square root of n c .

The primal and dual SDCO problems

The primal (P ) and dual (D) self-dual conic optimization (SDCO) problems read

(P )    inf c, x E A(x) = b x ∈ K and (D)    sup b T y A * (y) + s = c s ∈ K, (1.7)
where c ∈ E, A : E → R m is a linear map, b ∈ R m , and A * : R m → E is the adjoint of A when R m is equipped with the Euclidean scalar product. These problems are Lagrangian dual to each other. They are convex in the sense that their objective is convex and their feasible set is also convex. When s = 0 and n l = 0, one recovers the linear optimization (LO) problem; when s = 1 and n l = 0, one recovers the standard semidefinite optimization (SDO) problem. The feasible sets of (P ) and (D) are respectively denoted by

F P := {x ∈ K : A(x) = b} and F D := {(y, s) ∈ R m × K : A * (y) + s = c},
and their strictly feasible sets are denoted by

F s P := {x ∈ K s : A(x) = b} and F s D := {(y, s) ∈ R m × K s : A * (y) + s = c}.
Accordingly, a point x ∈ F s P is said to be strictly feasible for (P ) and a pair (y, s) ∈ F s D is said to be strictly feasible for (D). We also introduce the following Cartesian products

F := F P × F D and F s := F s P × F s D .
The solution sets of these problems are denoted by Sol(P ) and Sol(D)

and their optimal values by val(P ) and val(D).

The duality gap is zero if val(D) = val(P ) (with possible infinite values) and, otherwise, is the positive difference val(P )val(D) > 0.

By the Riesz-Fréchet theorem, there are elements a i ∈ E such that for all x ∈ E, there holds

A(x) =    a 1 , x E . . . a m , x E    . (1.8)
Each a i ∈ E is therefore an assembling of s matrices a s,j i ∈ S n j , for j ∈ [1 : s], and a vector a l i in R n l . It will be also assumed that the Euclidean scalar product is used on R m , in which case A * (y) is a weighted sum of the vectors a i (see question 1.1):

A * (y) = m i=1 y i a i ∈ E, (1.9)
where the y i 's are the components of y ∈ R m .

The central path

It is known that when F s = ∅, z := (x, y, s) is a primal-dual solution to (P ) and (D) if and only if

   A * (y) + s = c, s ∈ K A(x) = b, x ∈ K x, s E = 0.
(1.10)

These may be considered as the necessary and sufficient optimality conditions of the convex problems (P ) or (D), under the "constraint qualification condition" F s = ∅ (it is known, indeed, that "F s = ∅ and the surjectivity of A" is equivalent to Robinson's constraint qualification [START_REF] Robinson | Stability theory for systems of inequalities, part II: differentiable nonlinear systems[END_REF] at one or any point of the feasible sets of the problems (P ) and (D) [START_REF] Ch | Fragments d'Optimisation Différentiable -Théorie et Algorithmes[END_REF]). It can be shown (see question 1.2) that for x and s ∈ K, there holds

x, s E = 0 ⇐⇒ x s = 0. (1.11)
Therefore, (1.10) also reads

   A * (y) + s = c, s ∈ K A(x) = b, x ∈ K x s = 0.
(1.12)

We prefer (1.12) to (1.10), since then the "number of equations" in (1.12) is equal to the "number of unknowns". Indeed, the triple (x, y, s) lies in E × R m × E and the equation values are in the same space when x s ∈ E (this is certainly the case when x s = 0).

When A is surjective, the central path is the (image of the) map µ ∈ R ++ → z(µ) ∈ E×R m ×E, where z = z(µ) is the unique solution to the perturbed optimality conditions (see question 1.3)

   A * (y) + s = c, s ∈ K s A(x) = b, x ∈ K s x s = µe.
(1.13)

When F s = ∅, the central path converges to a primal-dual solution, when µ ↓ 0. The goal of the path-following algorithms that we shall consider is to follow the central path to reach a solution asymptotically.

Two directions

Overview

Given a point z = (x, y, s) ∈ E × R m × E, the goal of one iteration of the primal-dual path-following interior-point (IP) algorithm we consider is to make a displacement d in E × R m × E towards a central point z(µ), for some µ 0, that is "closer" to the solution than the current iterate z.

This displacement d = (d x , d y , d s ) ∈ E × R m × E is
a Newton-like step on the perturbed system (1.13) without its conic conditions x ∈ K s and s ∈ K s (these will be enforced at each iteration), namely

   A * (y) + s = c A(x) = b x s = µe. (1.14)
Computing the standard Newton direction on this system raises a serious difficulty, however. When z is not on the central path, which is usually the case, the point x s may not be in E, since the matrix components of x s may not be symmetric matrices.

In that case, one should work with more equations than unknowns in (1.14), which a situation we wish to avoid, since the linearized system may have no solution. To tell it otherwise, one would like z + d to satisfy the system (1.14), which reads

   A * (d y ) + d s = r d A(d x ) = r p d x s + x d s = µe -(x s + d x d s ), (1.15)
where we have introduced the primal r p and dual r d residuals

r p := b -A(x) ∈ R m and r d := c -A * (y) -s ∈ E. (1.16)
The last equation of (1.15) is nonlinear in d. Its linearization is obtained by dropping the term d x d s , resulting in the system

   A * (d y ) + d s = r d A(d x ) = r p d x s + x d s = µe -x s.
(1.17)

If d s is in E by the first equation of (1.17) (since A * takes its values in E, and c and s ∈ E), there are examples [START_REF] Lin | A predictor-corrector method for semi-definite linear programming[END_REF] showing that d x might not be in E (because x d s s -1 in the third equation is not necessarily in E, by lack of symmetry of its matrix components). For this reason a "symmetrization" operation must be introduced in addition to the linearization of the system, in order to keep x + d x in E. This can be done using a large number of methods, leading to many different IP directions; not all of them have the desired properties. We present below two of the most often implemented techniques: the NT and HKM directions. These directions have the following features: they preserve the primal and dual displacements d x and d s in E, they are uniquely defined at any strictly feasible iterate z ∈ F s , they are less computationally demanding than most other directions.

The NT direction

This section presents the Nesterov-Todd (NT) direction, called that way because of [35, 36, 47; 1997-1998]. This direction has the additional property of being scale invariant.

Derivation of the NT direction

The NT direction d = (d x , d y , d s ) ∈ E × R m × E at z = (x, y, s) ∈ K s × R m × K s
can be shortly presented as the result of three operations [START_REF] De Klerk | Aspects of Semidefinite Programming -Interior Point Algorithms and Selected Applications[END_REF].

1. Scaling. A weight w ∈ K s is introduced, which is the unique element in K s that satisfies the identity (see question 1.4)

w s w = x. (1.18)
It is given by the formulas

w := x 1/2 x 1/2 s x 1/2 -1/2 x 1/2 = s -1/2 s 1/2 x s 1/2 1/2 s -1/2 . (1.19)
Let us also introduce the scaled variable

v := w -1/2 x w -1/2 = w 1/2 s w 1/2 ∈ K s . (1.20)
Note that the vector component of w and v simply read for some µ 0 to be defined. After a left-Hadamard-multiplication by w -1/2 and a right-Hadamard-multiplication by w 1/2 , this equation takes the equivalent form We now symmetrize the left-hand side of this equation as follows

w l = (x l ) 1/2 • (s l ) -1/2 and v l = (x l ) 1/2 • (s l ) 1/2 . ( 1
1 2 (v + dx ) (v + ds ) + (v + ds ) (v + dx ) = µe.
The left-hand side is now an element in E.

3. Pseudo-linearization (the term "pseudo" is used because the weight w, which depends on the current iterate, is not linearized). The linearization consists in dropping the nonlinear terms from the last identity, namely the Hadamard products dx ds and ds dx . The identity then becomes

1 2 v ( dx + ds ) + ( dx + ds ) v = µe -v v.
This is a Lyapunov equation in the unknown dx + ds . Since v ∈ K s , it has a unique solution, which is dx + ds = µv -1 -v.

Left and right-Hadamard-multiplication of the two sides of this last equation by w 1/2 yield the symmetric pseudo-linearization at z of the third equation in (1.14):

d x + w d s w = µs -1 -x.
Therefore, the NT direction is obtained by solving the following linear system in

d = (d x , d y , d s ) ∈ E × R m × E: A * (d y ) + d s = r d A(d x ) = r p d x + w d s w = µs -1 -x, (1.22)
in which the value of µ still needs to be specified. It is a consequence of the calulation made in section 1.2.2, that, when A is surjective, the system (1.22) has a unique solution, whatever the right-hand side is (in particular, whatever µ > 0 is).

Computation of the weight w

The weight w given by (1.19) intervenes in the computation of the NT direction and must be computed explicitly. The vector component w l of w is easily computed by (1.21), so that we only consider a particular matrix component w s,j of w below.

Each matrix components w s,j of w can be computed using two Cholesky factorizations and one singular value decompositions (SVD) [47; § 4.1]. The Cholesky factorizations are those of the matrix components in x ∈ K s and s ∈ K s :

x s,j = L x s,j L T x s,j and s s,j = L s s,j L T s s,j , where L x s,j and L s s,j are lower triangular nonsingular matrices of order n j . The SVD is the one of L T s s,j L x s,j :

L T s s,j L x s,j = U j Σ j V T j .
where U j and V j are orthogonal matrices and Σ j is the diagonal matrix formed of the positive singular values of the nonsingular matrix L T s s,j L x s,j . Define

Q j := L -1 x s,j (x s,j ) 1/2 ,
which is an orthogonal matrix. It needs not be computed. Now, there holds

(x s,j ) 1/2 s s,j (x s,j ) 1/2 = (x s,j ) 1/2 L -T x s,j L T x s,j I n j L s s,j L T s s,j s s,j L x s,j L -1
x s,j

I n j (x s,j ) 1/2 = Q T j L T x s,j L s s,j V j Σ j U T j L T s s,j L x s,j U j Σ j V T j Q j = Q T j V j Σ 2 j V T j Q j .
From this identity and the orthogonality of Q T j V j , it results

((x s,j ) 1/2 s s,j (x s,j ) 1/2 ) -1/2 = Q T j V j Σ -1 j V T j Q j .
From the first identity in (1.18), we get

w s,j = (x s,j ) 1/2 Q T j L x s,j V j Σ -1 j V T j Q j (x s,j ) 1/2 L T x s,j = L x s,j V j Σ -1 j V T j L T x s,j .
Finally

w s,j = G j G T j , where G j := L x s,j V j Σ -1/2 j = L -T s s,j U j Σ 1/2 j . (1.23)
The matrices G j and their transpose intervene below in (1.29), for which it is useful to introduce the vectors in F :

g := (G 1 , . . . , G s , (w l ) 1/2 ) ∈ F, (1.24)
allowing us to write (recall the definition of g T in (1.3))

w = g g T ∈ E. (1.25)
It will also be useful to note that

g T s g = g -1 x g -T = σ and g -1 x s g = σ 2 , (1.26)
where the matrix components of σ ∈ E are the diagonal positive definite matrices Σ j .

Computation of the NT direction

A standard way of solving the system (1.22) defining the NT direction d consists in eliminating first d x using the second and third equation, to get

A(w d s w) = -r p + A(µs -1 -x).
Next, eliminating d s using the first equation, we get

A(w A * (d y ) w) = r p + A(w r d w) + A(x -µs -1
M(d y ) = r p + A(w r d w) + A(x -µs -1 ), (1.28)
where M :

d y ∈ R m → A(w A * (d y ) w) ∈ R m
is a symmetric positive definite (hence nonsingular) linear map and to form the "compact" (when m is small) m × m coefficient matrix M. For i and j ∈ [1 : m], there holds

M kl = [A(w A * (e l ) w)] k = a k , w a l w E [(1.8) and (1.9)] = a k , g g T a l g g T E [(1.25)] = g T a k g, g T a l g E . (1.29)
Then, one takes the Cholesky factorization M = L M L T M and solve (1.28). Forming M mainly requires 2mn 3 operations (for getting the products g T a k g, with i ∈ [1 : m]). When n is large, this is the most computationally expensive part of the solver, limiting its use to problems with n not exceeding a few hundreds. Problem with higher dimension must be sparse to be considered by SCDO, and this one must be able to consider such sparse problems, which is not the case of the solver presented in these notes.

Another way of computing d y is to observe that (1.27) or . . . Most implementations use the first approach.

The HKM direction

Another efficient direction is the so-called HKM direction (obtained independently in [START_REF] Helmberg | An interior-point method for semidefinite programming[END_REF][START_REF] Kojima | Interior-point methods for the monotone semidefinite linear complementarity problem in symmetric matrices[END_REF], see also [START_REF] Rendl | A simplified/improved HKM direction for certain classes of semidefinite programming[END_REF]). It is also implemented in SDPT3, for instance. . . .

Implementation

Data representation

We propose to follow the apparent data structure used in SeDuMi.

An element x ∈ E is represented by a large vector of dimension n F . This large vector is formed of the matrices x s,j ∈ S n j represented by a vector containing its n j × n j elements (not only its symmetric part having n j (n j + 1)/2 elements) and the vector

x l ∈ R n l .
This representation has the advantage of making many computations easy. For instance the scalar product x, s , for x and s ∈ E, is obtained in Matlab by x'*s if x is the n F -vector representing x and s is the n F -vector representing s. This makes the design of the piece of software easier and also makes it more efficient computationally (since there is no loop).

The above representation of the elements of E makes it necessary to often switch between the vector and matrix representations of matrices. This can be done with the CVX functions mat and vec (see below), probably inexpensively if the result is not safeguarded.

The linear operator A : E → R m introduced in the primal and dual problems in (1.7) is first represented by the vectors a i ∈ E, for i ∈ [1 : m], using (1.8). Then each vector a i ∈ E is represented as an n F -vector, as described in the first point, whose transpose forms the ith row of a matrix A.

With these representations of A and x, the value A(x) is obtained in Matlab by A*x and the value of A * (y) is obtained by A'*y, which are elegant expressions, easy to program and debug, and computationally efficient.

Calling statement

We propose to give to the sdco solver, the following Matlab function structure

[x,y,s,info] = sdco (A,b,c,K,x0,y0,options)
where A is an m × n F representation of the linear operator A, as described in section 1.3.1; b is an m real vector, representing the right hand side of the equality constraint in (P ); c is an n F -vector representing the vector c ∈ E, which determines the linear objective in (P ); K is a Matlab structure describing the structure of the vectors x ∈ E:

• K.s is the Matlab vector [n1; n2; ...; ns], providing the orders nj = n j , for j ∈ [1 : s], of the matrices x s,j ; • K.l is the length n l of the vector x l ;

1. Directions of displacement x0 is an n F -vector representing the vector x 0 ∈ E, which is the initial guess of the primal solution x; y0 is an m-vector, which is the initial guess of the dual solution y ∈ R m ; s 0 is then deduced form y 0 by s 0 := c -A * (y 0 ); options is a structure of possible options, which will be described in other session; this one can be used to select some features of the algorithm, like the threshold to decide optimality;

x is an n F -vector representing a vector x ∈ E, which gives the primal solution x if any; y is an m vector, giving the dual solution y if any; s is an n F -vector, giving the dual solution s if any; it is linked to y by s = c -A * (y); info is a structure providing information of the run (success, failure, primal/dual infeasibility, primal/dual unboundedness, duality gap, etc).

Matlab functions

Here are a few useful Matlab functions. L = chol(A,'lower') computes the Cholesky factorization of A = LL T , where L is lower triangular (only the lower triangular part of A is used).

[U,S,V] = svd(A) computes the singular value decomposition (SVD) of a real n F × n F matrix A = U SV T , where U and V are orthogonal matrices and S is diagonal with nonnegative entries (the singular values of A).

With CVX, M = mat(v) converts an n 2 vector v into an n × n matrix M, such that the columns of M are stacked below each other in v. This operation is frequently used, so that it is appropriate to write a personnal version, using reshape, if the function is not available. The reverse operation is realized by vec.

[Q,R] = qr(A,0) computes the QR factorization of a real m × n matrix A = QR, where Q is orthogonal and R is upper triangular. The second zero argument is useful when m > n (our case), in which case, only the first n column of Q and the first n rows of R are computed. With CVX, v = vec(M) converts an m × n matrix M into an mn vector v, stacking the columns of M below each other in v. This operation is frequently used, so that it is appropriate to write a personnal version, using reshape, if the function is not available. The reverse operation is realized by mat.

It is useful to introduce a Matlab function computing the Hadamard product of two vectors x and s in E:

[hdot] = sdco_hdot (x,s,K)

Using this function and introducing a matrix B whose kth column is g T a k g, the expression (1.29) of M simply reads B T B.

Recommendations

In order to get convergence of the algorithm of the next section, it is imperative to have a correctly computed direction. There are very few things that can certify the validity of the realized computation. We can recommend to verify that the vector w verifies (1.18), the matrix M defined by (1.29) is symmetric and positive definite, d x , d s = 0 (this is only true when the current iterate (x, y, s) is feasible, meaning that r p = 0 and r d = 0), the residual of the system (1.22) is close to zero.

Test cases

Test case 1a: an easy SDO problem of dimension 3

Consider the semidefinite optimization (SDO) problem in y ∈ R 3 , written in standard dual form

       sup e T y   1 -y 1 -y 3 -y 2 -y 3 1 -y 2 0 -y 2 0 1 -y 3   0,
where e ∈ R 3 is the vector of all ones. It is easily verified that x 0 = I 3 and (y 0 , s 0 ) = (0, I 3 ) are strictly feasible for the primal and dual problems respectively, so that these problems have a solution without duality gap. There holds b T y 0 = 0 < 3 = c, x 0 , so that z 0 := (x 0 , y 0 , s 0 ) is not a primal-dual solution.

Test case 1b: a simple LO problem of dimension 2

Consider the following linear optimization (LO) problem in only 2 variables and a single affine constraint, which reads

   inf x 1 + x 2 x 1 + 2x 2 = 1 x 0.
The initial primal and dual variables are

x 0 := 1/3 1/3 and y 0 := 0,
which are strictly feasible points.

1. Directions of displacement

Test case 1c: two SDO problems and one LO problem in parallel

The goal of this problem is to test the capacity of the developed solver to deal simultaneously with two positive semidefinite matrices and one nonnegative vector as primal variables. The problem consists in solving in parallel two copies of the semidefinite optimization problems of test case 1a, each one starting from a different matrix, and the additional linear optimization (LO) problem of test case 1c.

Denote by c the gradient of the objective of the primal expression of test case 1a and by A(x) = b its affine constraint. Then the primal expression of test case 1c is the

following problem in x = (x s,1 , x s,2 , x l ) ∈ S 3 × S 3 × R 2            inf c, x s,1 + c, x s,2 + (e 2 ) T x l A(x s,1 ) = b A(x s,2 ) = b x l 1 + 2x l 2 = 1 x s,1 0, x s,2 0, x l 0,
The initial primal and dual variables are

x s,1 0 := I 3 , x s,1 0 :=   1 0 1/4 0 1/2 0 1/4 0 1   , x l 0 := 5- √ 17 2 -3+ √ 17 4
, and

y 0 := (0 R 6 , -1+ √ 17 4
), which are strictly feasible points.

Notes

There are various Matlab SDCO solvers. Having a look at SeDuMi [START_REF] Sturm | Using SeDuMi 1.02, a Matlab toolbox for optimization over symmetric cones[END_REF][START_REF] Sedumi | [END_REF] is certainly a good idea, since sdco may be viewed as a reduced feature version of that solver. Another source of inspiration is SDPT3 [START_REF] Toh | On the implementation and usage of SDPT3a Matlab software package for semidefinite-quadratic-linear programming, version 4[END_REF].

Questions

When a question starts with a number within braces, the number gives an indication on the difficulty of the question, with respect to the chapter (some questions will be easier to answer when subsequent chapters will be known). This indication ranges from 1 to 5: 1 for easy-or-classical and short arguments, 2 for easy-or-classical arguments, 3 for arguments requiring specific knowledge, 4 for more difficult arguments, 5 for very difficult arguments or when a advanced computer program must be written to answer the question.

1.1. {1} Adjoint of A. Show that the adjoint of the linear map A : E → R m defined by (1.8) is given by (1.9) when R m is equipped with the Euclidean scalar product.

1.2. {2} Towards a square optimality system. Show the equivalence (1.11) when x and s are in the cone K.

1.3. {2} Unique central point. Show that (1.13) has a unique solution z = (x, y, s) ∈ E × R m × E, provided A is surjective.
1.4. {2} Weight w. Show that, when x and s ∈ K s , the identity (1.18) has a unique solution w ∈ K s , which is given by (1.19).

1.5. {5} On test cases 1a, 1b, and 1c. Consider test cases 1a, 1b, and 1c. Let z 0 := (x 0 , y 0 , s 0 ) and d 0 be the NT direction computed at z 0 for some µ. 1) {1} Is z 0 on the central path?

2) {5} Compute the NT direction for µ = 0 by a computer program.

3) {1} Is z 1 = z 0 + d 0 a primal-dual solution to the problem if µ = 0?

2 Moving along the central path

Algorithmic concepts

In this chapter, we force the iterates z = (x, y, s) to be strictly feasible, meaning that they will be in F s . For this to work, it will be necessary to have a first iterate that is strictly feasible, and even sufficiently close to the central path.

The closest central point

In a path-following algorithm, the first question that arises deals with the determination of the central point z(µ) to which the current iterate z := (x, y, s) is the closest. Indeed, in such an algorithm, the next iterate aims a point on the central path C and gets close to it by a Newton step, provided this central point is not too far from z. If we want the iteration to make a progress towards the solution it is necessary that this aimed central point be closer to the solution than the central point that is the closest to z, hence the necessity to determine the latter. Finding the closest central point is not a well defined concept, actually, since the central path is not a convex set, making the projection on C not well defined. To overcome the difficulty, we look instead to the closest central point in a transformed space in which the central path becomes a half-line. The transformation is

τ : z = (x, y, s) ∈ F → x s ∈ F.
This transformation is certainly not a bijection, but it turns out to be useful, which is enough to make it appropriate. Since C := {z ∈ F : x s ∈ R ++ e}, it follows that τ (C ) = R ++ e. Then, by the convexity and closedness of τ (C ), it makes sense to look for the point in τ (C ) that is the closest to τ (z). That problem simply reads

min µ 0 x s -µe 2 F . (2.1)
It has a unique solution (see question 2.1), which is

µ(z) := x, s E n c , (2.2a)
where n c is given by (1.6). The closest central point is therefore "considered" to be the one on the central path with the parameter µ = µ(z). This is not irrelevant, since when z ∈ C , it is clear that the closest central point to z in the preceding sense is z itself; reassuring.

Neighborhood of the central path

Experimentation has convinced the algorithmicians that it is not good to let the iterates going too far from the central path C or, more correctly, too close to the boundary of the feasible set F (it is not the same thing actually, since the "central" path may be close to the boundary of F [START_REF] Deza | The central path visits all the vertices of the Klee-Minty cube[END_REF]). When z is close to that boundary, the stepsize along the Newton direction may be very small, so much as to prevent any progress to the solution: the iterates get stuck on the boundary of F. For this reason, the iterates are maintained in some neighborhood of the central path, often the one that we now define.

Recall that the central path is an analytic concept (not a geometric one, depending only on the feasible set), so that this is a second best solution.

Observe first that for a point z ∈ F s :

x s = µe ⇐⇒ v v = µe ⇐⇒ µ 1/2 v -1 = µ -1/2 v, (2.3)
where v := w -1/2 x w -1/2 = w 1/2 s w 1/2 ∈ K s was already defined in (1.20) and w ∈ K s was already defined in (1.19) as the unique solution to (1.18). Then one defines a proximity measure of the central point z(µ) by ( [START_REF] Jiang | A long step primal-dual path following method for semidefinite programming[END_REF] and [8; § 7.1]):

δ µ : z ∈ F s → δ µ (z) := 1 2 µ 1/2 v -1 -µ -1/2 v E ∈ R. (2.4)
Observe now that for a point z ∈ F s :

z ∈ C ⇐⇒ x s = µ(z)e. (2.5)
Therefore, it makes sense to define a proximity measure of the central path by

δ : z ∈ F s → δ(z) := δ µ(z) (z). (2.6)
Accordingly, for θ ∈ [0, 1), one introduces the following neighborhood of the central path V(θ) := {z ∈ F s : δ(z) θ}.

Using (1.4b) and w = g g T from (1.25), we deduce that

δ(z) = 1 2 g T µ(z) 1/2 x -1 -µ(z) -1/2 s g E , (2.7)
so that the distance can be computed efficiently.

Two algorithms

A predictor-corrector algorithm

The predictor-corrector algorithms are the most often implemented interior-point methods. We essentially give the description of the Mizuno-Todd-Ye predictor-corrector method [START_REF] Mizuno | On adaptive-step primal-dual interior-point algorithms for linear programming[END_REF][START_REF] De Klerk | Aspects of Semidefinite Programming -Interior Point Algorithms and Selected Applications[END_REF], whose iteration is composed of two phases: a corrector phase, followed by a predictor phase (in reverse order, in a way).

The corrector phase starts with an iterate z ∈ V(1/3). Its goal is to compute an intermediate iterate z ′ , close enough to the central path, so that the predictor step that follows will compute a subsequent iterate z + again in the neighborhood V(1/3). This goal is achieved by making a displacement along the NT direction with µ = µ(z), say d c (known as the centering direction), and a unit stepsize. The intermediate iterate is then

z ′ := z + d c .
The predictor phase, which follows the corrector step, starts with z ′ ∈ F s and consists in making a displacement along the NT direction with µ = 0, say d a (known as the affine-scaling direction). The stepsize α > 0 along that direction is given by the formula (no need of linesearch, provided a condition given below is fulfilled)

α = 2 1 + [1 + 13 da x da s + da s da x E /(2µ(z ′ ))] 1/2 , ( 2.8a)
where

da x = w -1/2 d a x w -1/2
and da s = w 1/2 d a s w 1/2 . Note that the matrix w is computed at z ′ , not at z (hence by (1.19) with (x, s) changed into (x ′ , s ′ )). It can be shown that

da x da s + da s da x E = g -1 d a x d a s g + g T d a s d a x g -T E , (2.8b)
where g and g T are defined by (1.24), so that the computation of this norm can be done efficiently.

Note that for the vector parts of d a x and d a s , there holds

da x l • da s l = da s l • da x l = d a x l • d a s l ,
so that the weight w does not intervene.

The stepsize (2.8) is small enough to ensure that the next iterate

z + := z ′ + αd a
is in the neighborhood V(1/3) and large enough to ensure the polynomiality of the algorithm (see below).

Algorithm 2.2.1 (predictor-corrector) A tolerance ε > 0 on µ(z) is given to determine when stopping the iterations. The iteration from z to z + starts at some z ∈ V(1/3).

1. Stopping test. If µ(z) ε, stop.

2. Corrector phase. Compute the NT direction at z with µ = µ(z), denote it d c (centering direction). Take as next iterate

z ′ := z + d c .
3. Predictor phase. Compute the NT direction at z ′ with µ = 0, denote it d a (affine-scaling direction), and the stepsize α given by (2.8). Set

z + := z ′ + αd a .
It can be shown [START_REF] De Klerk | Aspects of Semidefinite Programming -Interior Point Algorithms and Selected Applications[END_REF] that α given by (2.8) satisfies

α 2 1 + [1 + 13n c /2] 1/2 ,
which implies the polynomiality of the algorithm, with the iteration complexity

K := 1 + 1 + 13n c 2 log µ(z 0 ) ε . (2.9)
This means that if the algorithm starts with z 0 ∈ V(1/3), it ends up with z K satisfying µ(z K ) ε, where the number K of iterations is given by (2.9).

A large step algorithm

Implementation

Here are a few recommendations.

It is important to check numerically that the following formula is verified after each step αd from a point z ∈ F s , where α > 0 and d is the NT direction computed for a certain µ (see question 2.2):

µ(z + αd) = (1 -α)µ(z) + αµ. (2.10)
In particular, µ(z ′ ) = µ(z) and µ(z

+ ) = (1 -α)µ(z),
showing that the progress to the solution is made by the prediction step. Observe indeed that z * is a solution if and only if z * ∈ F s and µ(z * ) = 0, so that the goal of the algorithm is to force the decrease of µ(z).

Note finally that if µ = µ(z), then (2.10) shows that µ(z + αd) is the constant µ(z), whatever α 0 is.

It is important to check that the iterates z remains in the neighborhood V(1/3).

Even if the distance to the central path does not intervene in the predictor-corrector algorithm, in the development phase, it is useful to have a function that can compute the distance to the central path and to verify that this one remains below 1/3.

The value of α given by (2.8) may be very close to one, or even equal to one, due to rounding errors. This usually yields vectors x and/or s that are no longer in the cone K s . In this case, the algorithm gets stuck in a Cholesky factorization. Limiting this computed value of α to options.max_stepsize = 0.999 or so, may help the algorithm to go further and compute a more precise solution.

Of course, taking a too small value for this parameter will prevent the algorithm from converging rapidly.

Test case

Test case 2: minimum matrix norm

We consider the minimum matrix norm problem which reads

min v∈R p B(v) 2 , (2.11)
where B(v) := B 0 + p i=1 v i B i , the matrices B i ∈ R q×r and • 2 denotes the ℓ 2 -norm. We assume that the B i 's are all different, which is a quite reasonable assumption. This problem is convex, but nonsmooth. It can be expressed as an SDO problem as follows.

Problem (2.11) can be rewritten min t,v {t : B(v) 2 2 t 2 }. Using the Schur complement technique, it can be verified that this last problem can also be written as the following SDO problem in dual form [START_REF] Vandenberghe | Semidefinite programming[END_REF]:

   max (t,v)∈R×R p -t tI q B(v) B(v) T tI r 0. (2.12)
This problem has dimension n = q + r and m = p + 1. One can also easily find a strictly feasible primal-dual point (question 2.4). What happens with the SDCO solver when two B i 's are identical? Can you explain its behavior?

Notes

The proximity measure δ defined by (2.6) was introduced by Jiang [START_REF] Jiang | A long step primal-dual path following method for semidefinite programming[END_REF], extending to semidefinite optimization a similar measure introduced by Jansen et al. [START_REF] Janssen | Primal-dual algorithms for linear programming based on the logarithmic barrier method[END_REF] in linear optimization.

Questions 2.1. {1} Formula of µ(z). Show that the solution to problem (2.1) is given by (2.2).

{2}

Evolution of µ(z) along the NT direction. Show formula (2.10).

{5}

On test cases 1a, 1b, and 1c. Using your solver, find a solution to test cases 1a, 1b, and 1c.

2.4. {5} On test case 2.

1) {1} Show that, in standard notation, z 0 = (x 0 , y 0 , s 0 ), with

x 0 = 1 q+r I q+r , y 0 = (t 0 , v 0 ) = ( B 0 2 + 1, 0), and s 0 =

t 0 I q B 0 B T 0 t 0 I r
is a strictly feasible primal-dual point for problem (2.12). 2) {5} Using your solver, find a solution to problem (2.11), equivalent to problem (2.12), with, say p = 2, q = 2, r = 2, and random data for the matrices

B i ∈ R q×r , with i ∈ [0 : 2].
The predictor-corrector algorithm that was presented in section 2.2.1 assumes that the first iterate is strictly feasible and is in some neighborhood V(θ) of the central path.

In this section, we consider a method to get such a point, provided a strictly feasible primal-dual point z 0 = (x 0 , y 0 , s 0 ) is known. The idea is to minimize a primal-dual merit function that has a central point as unique minimizer, starting the minimization process at z 0 . The NT direction can be used to force the decrease of that merit function.

Getting a feasible point close to the central path

We use on E × R m × E the scalar product inherited from those on E and R m , which, for z = (x, y, s) and z ′ = (x ′ , y ′ , s ′ ), is denoted by and takes the value

z, z ′ = x, x ′ E + y T y ′ + s, s ′ E .
The associated norm is denoted by • .

A primal-dual merit function

We investigate the simple idea that consists in minimizing approximately on F s a function

ϕ µ : E × R m × E → R ∪ {+∞},
whose unique minimizer in F s is the central point z(µ), that is the unique solution to the perturbed optimality system (1.13). There are several possibilities for ϕ µ , including a function defined only on E (see the indication given for solving exercise 1.3). We prefer, however, using a function defined on E × R m × E that can be minimized using the NT direction, which has already been implemented. With that objective in mind, we follow [8; § 7.1] by defining the function

ϕ µ at z ∈ E × R m × E by ϕ µ (z) := x, s E + µ ψ(x s), (3.1)
where µ > 0 is a fixed parameter and ψ :

F → R ∪ {+∞} is defined at v ∈ F by ψ(v) = s j=1 ld v s,j - n l i=1 log v l i .
This expression makes use of the self-concordant barrier [START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF][START_REF] Renegar | A Mathematical View of Interior-Point Methods in Convex Optimization[END_REF] ld : R n×n → R ∪ {+∞} of the cone of positive definite matrices that are not necessarily symmetric, defined at

M ∈ R n×n by ld(M ) = -log det M if det M > 0 +∞ otherwise. (3.2)
It can be shown that, when F s = ∅ and A is surjective, ϕ µ is strictly convex on F s and has the central point z(µ) as unique minimizer on F s (see exercise 3.1). This has for consequence that generating iterates in F s that minimize ϕ µ will eventually provide a point in V(θ) for some prescribed θ > 0.

The gradient and Hessian of ϕ µ at z ∈ E × R m × E can be computed by

∇ϕ µ (z) =   s -µx -1 0 x -µs -1   and ∇ 2 ϕ µ (z) d =   µ x -1 d x x -1 + d s 0 d x + µ s -1 d x s -1   , (3.3)
where

d ∈ E × R m × E.

Use of the NT direction

Now, the question arises to see whether ϕ µ can be minimized using the NT direction.

The next result shows that the NT direction d defined at z ∈ F s , with the parameter µ > 0, is a descent direction of ϕ µ . It also provides with (3.4) the value of the directional derivative of ϕ µ along the NT direction.

Proposition 3.1.1 (decrease of the primal-dual barrier along the NT direction) Let d be the NT direction at a point z ∈ F s for some µ > 0. Then

∇ϕ µ (z), d = -4 µ δ µ (z) 2 , (3.4)
where δ µ is defined by (2.4). Furthermore, the following properties are equivalent:

(i) z = z(µ), (ii) d = 0, (iii) ∇ϕ µ (z), d < 0.
Proof. Using point 1 of exercise 3.1, one gets

∇ϕ µ (z), d = s -µx -1 , d x E + x -µs -1 , d s E .
Using the scaling vector w defined by (1.19) and the scaled directions dx := w -1/2 d x w -1/2 and ds := w 1/2 d s w 1/2 , one gets

∇ϕ µ (z), d = s -µx -1 , w 1/2 dx w 1/2 E + x -µs -1 , w -1/2 ds w -1/2 E .
It is known, from (1.20), that v := w -1/2 x w -1/2 = w 1/2 s w 1/2 and, from (1.22) 3 , that dx + ds = µv -1 -v, so that

∇ϕ µ (z), d = v -µv -1 , dx + ds E = -v -µv -1 2 E = -µ µ -1/2 v -µ 1/2 v -1 2 E = -4µδ µ (z) 2 [(2.4)],
which is (3.4). The next equivalences follow easily. ✷

Note that if µ is set to µ(z), if d is the NT direction computed for that µ, and if a stepsize α > 0 is taken along d to force the decrease of ϕ µ , then µ(z + αd) = µ(z) by (2.10), so that an NT direction at z + αd with µ = µ(z + αd) will also be a descent direction of the same merit function ϕ µ . In other words, if one sets µ = µ(z) at each iterate z, the NT direction is a descent direction for the same merit function ϕ µ . Proposition 3.1.2 (distance to the central path after a NT step) Let d be the NT direction at a point z ∈ F s for µ := µ(z). Then the distance δ to the central path satisfies the inequality

δ(z + d) δ(z) 2 2(1 -δ(z) 2 ) . (3.5)
Proof. Adapt the proof of [8; lemma 7.4] to the present framework. ✷

An algorithm

A consequence of the fact that z(µ) uniquely minimizes ϕ µ (exercise 3.1) and proposition 3.1.1 is that, to get closer to the central path, it makes sense to minimize the function ϕ µ , with µ = µ(z), along the NT direction, using a linesearch technique. Furthermore, from (3.5) in proposition 3.1.2,

δ(z) 2τ 2 1 + 2τ 2 =⇒ δ(z + d) τ δ(z),
so that, choosing τ < 1, the distance decreases linearly without the need to make linesearch, as soon as δ(z) 2τ 2 /(1 + 2τ 2 ), which is < 2/3. Inequality (3.5) also tells us that the distance δ(z) converges to zero quadratically.

This discussion leads to the following algorithm, which computes a point in the neighborhood V(θ) of the central path, starting from some given point z ∈ F s . Algorithm 3.1.3 (getting a point in V(θ)) The algorithm uses three parameters: the scalar θ > 0 specifying the neighborhood to reach, a desired linear speed of convergence τ < 1 close to 1, and a linesearch parameter ω > 0 close to zero. Typically θ = 1/3, τ ≃ 0.99, and ω ≃ 10 -4 . Let z ∈ F s .

Repeat until z ∈ V(θ).

1. NT direction. Compute the NT direction at z with µ := µ(z), denote it d c .

Stepsize α. If δ(z)

2τ 2 /(1 + 2τ 2 ), take α = 1. Otherwise, find the largest α of the form 2 -i with i ∈ N such that

ϕ µ (z + αd c ) ϕ µ (z) -4 ω µ α δ(z) 2 .
(3.6)

3. Next iterate. Set the next iterate z to z + αd c .

The previous algorithm hides a serious difficulty, which may occur during the linesearch in (3.6), which evaluates the quality of the trial stepsize α. By definition, the function ϕ µ (z) = +∞ when z / ∈ K s . This fact may not be seen by simply evaluating x s,j s s,j (resp. x l i s l i ), which may be positive definite (resp. positive) even when x s,j | ≻ 0 or s s,j | ≻ 0 (resp. when both x l i < 0 and s l i < 0). Therefore, the correct evaluation of ϕ µ requires to check that, for all j ∈ [1 : s] and all i ∈ [1 : n l ] : x s,j ≻ 0 and s s,j ≻ 0 (either by the computation of the minimum eigenvalues of x s,j and s s,j or by doing their Cholesky factorizations) and that x l i > 0 and s l i > 0. where p, q, r are the parameters (p, q, r).

Notes

In this chapter, we consider algorithms that do not generate iterates in the feasible set F. It is interesting to have such algorithms since in practice it is usually not obvious to have a strating point in F. Two approaches are considered.

In the big M approach of section 4.1, the starting iterate z 0 = (x 0 , y 0 , s 0 ) is supposed to satisfy the affine constraints A(x 0 ) = b and A * (y 0 ) + s 0 = c but x 0 and s 0 are not supposed to belong to the cone K s . Such an infeasible starting point is easy to find by linear algebra techniques. One then constructs from (P ) and (D), other primal-dual pairs of problems, for which z 0 is strictly feasible. These pairs have parameters that must be driven to zero. Hence they are solved repetitively by the robust techniques used in the previous chapters until the parameters vanish.

The technique explored in section 4.2 discards the use of a neighborhood of the central path to control the iterates. The main reason is that this control is not easy to implement. For example, the stepsize in (2.8), which ensures the iterates to stay in the prescribed neighborhood, is conservative, meaning that one could take a larger stepsize (hence converge more rapidly) and stay in the neighborhood. The presented algorithm has a predictor-corrector nature and has the main goal of being practically efficient (while the theoretical convergence issues are ignored). Algorithms implemented in most available codes follow a similar approach. They use heuristics, but are faster than approaches whose polynomial complexity is well documented.

The big M approach

In all this section, we assume that an initial x 0 ∈ E is known, which satisfies the equality constraint

A(x 0 ) = b. (4.1)
Having a pair (y 0 , s 0 ) satisfying the affine constraint of the dual problem makes no difficulty, since y 0 can be taken arbitrarily and s 0 can be set to c -A * (y 0 ). A way of getting a point x 0 ∈ E satisfying (4.1) is described in section 4.1.1.

We then discuss the so-called big M methods for solving the SDCO problem from a point z 0 = (x 0 , y 0 , s 0 ) satisfying (4.1) and A * (y 0 ) + s 0 = c but not the conditions x 0 ∈ K s and s 0 ∈ K s that would make z 0 strictly feasible and for which the algorithm of chapter 3 could be used. The case when x 0 ∈ K s and s 0 / ∈ K s is considered in section 4.1.2, the case when x 0 / ∈ K s and s 0 ∈ K s is considered in section 4.1.3, and the case when x 0 / ∈ K s and s 0 / ∈ K s is considered in section 4.1.4.

The contents of this section is adapted from [50; § 6].

Getting primal affine feasibility

We assume below that A : E → R m has the following natural extension from E to F, denoted

A F : F → R m .
This one is constructed as follows. First, we take the expression (1.8) of A, with vectors

a i ∈ E, for i ∈ [1 : m]. Then, for x ∈ F, each component [A F (x)] i of A F (x) is defined to be [A F (x)] i = a i , x F . Obviously, the restriction of A F to E is A,
which is what we mean by "extension". Now, if x 0 ∈ E satisfying (4.1) is not given by the user of the solver, this one should compute such a point or claim that (4.1) has no solution x 0 ∈ E. The solver can proceed as follows.

The first step is to check whether b ∈ R(A F ). This can be done using linear algebra techniques. If this is not the case, the primal problem is infeasible and there is no reason to go further. Otherwise, we claim that (4.1) has a solution (in E). This one can be obtained by first computing x ∈ F satisfying the affine constraint

A F (x) = b (4.2a)
and then symmetrizing it by taking

x l 0 := x l and x s,j 0 = 1 2 x s,j + (x s,j ) T , ∀ j ∈ [1 : s]. (4.2b)
It is asked in question ?? to verify that x 0 ∈ E obtained in that manner is indeed a solution to (4.1).

Starting from a strictly feasible primal point

We assume in this section that a strictly feasible primal point x 0 ∈ F s P is known, meaning that

A(x 0 ) = b and x 0 ∈ K s , (4.3)
but that s 0 / ∈ K s . Generally speaking, the big M method to solve the SDCO problem when a strictly feasible point is not available, using the technique of chapter 3, consists in introducing a modified SDCO problem that has the same solutions as the original problem, provided a parameter M in the modified SDCO problem is "large enough" (hence the words "big M "), and for which it is easy to find a strictly feasible point. The technique has therefore an effect that is similar to an exact penalty method.

In case (4.3) holds, the modified SDCO problem is obtained from a modification of the dual problem, for which a strictly feasible point is not known. The dual problem becomes the problem in (y, η, s) ∈ R m × R × E (the modified parts are in red):

       sup b T y -M 1 η A * (y) -ηe + s = c s ∈ K η 0. (4.4)
The new parameter η is used to realize feasibility easily (see below). This nonnegative parameter is then forced to be as small as possible by taking M 1 "sufficiently" large (the rational is that, when η = 0, the original dual problem (D) is recovered). We claim that problem (4.4) can be cast as a standard dual SDCO problem, for which a strictly feasible point is easy to determine, its primal-dual solutions are the same as the original problem, provided M 1 is chosen "sufficiently large" and the dual problem is feasible.

Let us see this.

One can cast problem (4.4) as the standard dual SDCO problem with Ẽ := E × R, K = K × R + , and the dual variable (ỹ, s) ∈ R m+1 × Ẽ that solves It is easy to get a strictly feasible point for this model by taking

   sup bT ỹ Ã * (ỹ) + s = c s ∈ K.
y 0 arbitrary, η 0 = σ 0 > [min K (c -A * (y 0 ))] -, and s 0 = c -A * (y 0 ) + η 0 e,
where min K has been defined by (1.5) and α -:= max(0, -α).

We know that (4.5) is the dual of

     inf c, x Ẽ Ã(x) = b x ∈ K, (4.6)
where

x = (x, ξ 1 ) ∈ E × R. Now, Ã : Ẽ → R m+1 is defined at x by Ã(x) = A(x) -e, x E -ξ 1 .
Therefore, problem (4.6) also reads

           inf c, x E A(x) = b e, x E + ξ 1 = M 1 x ∈ K ξ 1 0. (4.7)
A strictly feasible point (x 0 , ξ 1,0 ) for (4.6)-(4.7) can be the known x 0 and

ξ 1,0 = M 1 -e, x 0 E , provided M 1 > e, x 0 E .
The proposed approach consists in solving (4.7)-(4.4), or equivalently (4.6)-(4.5) in standard form, from a strictly feasible primal-dual pair hoping that at the primal-dual solution (x, ξ 1 , y, η, s, σ) the scalar η vanishes, in which case problem (P ) is actually solved. If η = 0, one increases M 1 and solve problem (4.7)-( 4.4) again. The process is stopped when η = 0 is found or when M 1 is considered as too large. The latter case may occur if (D) is infeasible. Algorithm 4.1.1 (getting a solution from a strictly feasible primal point) Let x ∈ F s P and let M 1 > e, x E . One iteration of the algorithm is as follows.

1. Solve (4.7)-(4.4) from a primal-dual strictly feasible point to get (x, ξ 1 , y, η, s, σ).

2. Successful stopping test. If η = 0, stop and declare that (x, y, s) is a primal-dual solution to (P ).

3. Failure stopping test. If M 1 is too large, stop and declare that it is likely that the dual problem (D) is infeasible.

4.

Increase M 1 .

Starting from a strictly feasible dual point

We assume in this section that a strictly feasible dual pair (y 0 , s 0 ) ∈ F s D is known, meaning that A * (y 0 ) + s 0 = c and s 0 ∈ K s , but that x 0 / ∈ K s , and we determine a modified SDCO problem such that it is easy to maintain (y 0 , s 0 ) strictly feasible and to determine a strictly feasible primal point x 0 .

Consider the following two equivalent modified primal problems in (x, ξ 2 ) ∈ E × R:

       inf c, x E + M 2 ξ 2 A(x) = b x + ξ 2 e ∈ K ξ 2 0 or        inf c, x E + (M 2 -c, e E ) ξ 2 A(x -ξ 2 e) = b x ∈ K ξ 2 0. (4.8)
The problem in the left-hand side in (4.8) shows that one tries to have an as small as possible perturbation of x in K, namely x + ξ e e, by forcing the nonnegative scalar ξ 2 to be as small as possible thanks to the positive number M 2 in the objective, which is taken "sufficiently" large. This is indeed desirable, since if ξ 2 = 0, the original primal problem (P ) is recovered. The equivalent problem in the right-hand side, obtained by the redefinition x + ξ 2 e x, is in a form closer to the standard primal SDCO problem. We claim that this problem can be cast as the standard primal SDCO problem, for which a strictly feasible point is easy to define, its primal-dual solutions are the same as the original problem, provided M 2 is chosen "sufficiently large" and the primal problem is feasible.

Let us see this.

One can cast the problem as the right-hand side of (4.8) in the standard primal SDCO problem in the variable x ∈ Ẽ := E × R:

   inf c, x Ẽ Ã(x) = b x ∈ K. (4.9)
This problem is identical to problem (4.8) provided the vectors c and x ∈ Ẽ, the linear map à : Ẽ → R m , and the cone K are defined by

c = c M 2 -c, e E , x = x ξ 2 , Ã(x) = A(x -ξ 2 e), and
K := K × R + .
A strictly feasible point for (4.9) is given by

x0 := x 0 + ξ 2,0 e ξ 2,0 , where x 0 ∈ E is an arbitrary solution to A(x 0 ) = b, ξ 2,0 > [min K (x 0 )] -.
We know that (4.9) has for dual the following problem in (y, s) ∈ R m × Ẽ: Writing

   sup b T y à * (y) + s = c s ∈ K.
s := s σ 2 ,
we see that problem (4.10) also reads

           sup b T y A * (y) + s = c -A(e) T y + σ 2 = M 2 -c, e E s ∈ K σ 2 0.
Eliminating σ 2 , the dual problem (4.10) becomes the following problem in (y,

s) ∈ R m × E:        sup b T y A * (y) + s = c s ∈ K e, s E M 2 . (4.11)
It can be shown that the variable ξ 2 in (4.8) is a dual variable associated with the constraint e, s E M 2 in (4.11). The given strictly feasible dual pair for (y 0 , s 0 ) for (D) is still a strictly feasible pair for the modified dual problems (4.10) provided

M 2 > e, s 0 E .
Note that this setting ensures that M 2 > 0 (since s 0 ∈ K s ), as desired.

The proposed approach consists in solving (4.8)-(4.11) from a strictly feasible primaldual pair (instead of the original problem (P )) hoping that at the solution (x, ξ 2 , y, s) the scalar ξ 2 vanishes, in which case problem (P ) is actually solved. If ξ 2 = 0, the penalty parameter M 2 is increased and problem (4.8) is solved again. The process is stopped when a solution with ξ 2 = 0 is found or when M 2 is considered as too large, which may occur if (P ) is infeasible. 2. Successful stopping test. If ξ 2 = 0, then stop and declare that (x, y, s) is a primal-dual solution to (P ).

3. Failure stopping test. If M 2 is too large, stop and declare that it is likely that the primal is infeasible.

4. Increase M 2 .

Starting without a strictly feasible point

If no strictly feasible primal and/or dual point is known, then one considers the modified primal problem, obtained by combining the methods of sections 4.1.2 and 4.1.3:

               inf c, x E + M 2 ξ 2 A(x) = b e, x E + ξ 1 = M 1 x + ξ 2 e ∈ K ξ 1 0 ξ 2 0 or                inf c, x E + (M 2 -c, e E )ξ 2 A(x -ξ 2 e) = b e, x E + ξ 1 -n c ξ 2 = M 1 x ∈ K ξ 1 0 ξ 2 0, (4.12)
in which both M 1 and M 2 are taken "sufficiently" large (initial values are given below).

The problem in the left-hand side in (4.12) satisfies the affine constraint A(x) = b (a constraint that can be easily verified by the technique proposed in section 4.1.1) without x ∈ K, but ensures that x + ξ 2 e is in K, where ξ 2 will be small when M 2 is taken large. Its second affine constraint e, x E + ξ 1 = M 1 is directly inspired from the one appearing in (4.7) and will allow us to find easily a strictly feasible point for the associated dual problem. The problem in the left-hand side of (4.12) is not in standard primal form because of the constraint x + ξ 2 e ∈ K, which is the reason why it is rewritten as the one in the right-hand side, which is obtained after the change of variable x + ξ 2 e x. Problem in the right-hand side of (4.12) can be cast as the standard primal SDCO problem in the variable

x ∈ K := K × R 2 + ⊂ Ẽ := E × R 2 :      inf c, x Ẽ Ã(x) = b x ∈ K, (4.13)
where the vectors c and x ∈ Ẽ, the linear map à : Ẽ → R m+1 , and the vector b ∈ R m+1 are defined by

c =   c 0 M 2 -c, e E   , x =   x ξ 1 ξ 2   , Ã(x) = A(x -ξ 2 e) -e, x E -ξ 1 + n c ξ 2 , and b = b -M 1 .
The opposite of the second constraint in (4.12) adopted in the definitions above is motivated by the desire to have -M 1 as the last component of b, which will result in requiring M 1 > 0, which goes in the same sense as the condition M 2 > 0 required above. A strictly feasible point for the modified primal problems (4.12) or (4.13) is

x0 :=   x 0 + ξ 2,0 e ξ 1,0 ξ 2,0   ,
where one takes in order:

x 0 ∈ E is an arbitrary solution to A(x 0 ) = b, ξ 2,0 > [min K (x 0 )] -, M 1 > [ e, x 0 E ] + , ξ 1,0 := M 1 -e, x 0 E , (4.14)
where min K has been defined by (1.5), α + = max(0, α), and α -= max(0, -α).

We know that the dual of (4.13) reads

   sup bT ỹ Ã * (ỹ) + s = c s ∈ K. (4.15)
where à * :

R m+1 → Ẽ is defined at ỹ = (y, η) ∈ R m × R by à * (ỹ) =   A * (y) -ηe -η -e, A * (y) E + n c η   . Therefore, with s = (s, σ 1 , σ 2 ), problem (4.15) reads                    sup b T y -M 1 η A * (y) -ηe + s = c -η + σ 1 = 0 -e, A * (y) E + n c η + σ 2 = M 2 -c, e E s ∈ K σ 1 0 σ 2 0 or                    sup b T y -M 1 η A * (y) -ηe + s = c -η + σ 1 = 0 e, s E + σ 2 = M 2 s ∈ K σ 1 0 σ 2 0. (4.16)
Recall that the parameter M 1 fixed by (4.14) is positive, which is a desired property, since one would like to have η = 0 in the above problem. A strictly feasible point for the modified dual problems (4.15) or (4.16) is ỹ0 = (y 0 , η 0 ) and s0 = (s 0 , σ 1,0 , σ 2,0 ), where one sets successively y 0 arbitrary,

η 0 = σ 1,0 > [min K (c -A * (y 0 ))] -, s 0 := c -A * (y 0 ) + η 0 e, M 2 > [ e, s 0 E ] + , σ 2,0 := M 2 -e, s 0 E , (4.17)
where min K has been defined by (1.5), α -= max(0, -α), and α + = max(0, α). The setting of M 2 above also ensures its positivity, which was a desired property of problems (4.12).

The goal of the following algorithm is to solve a sequence of primal-dual problems (4.13)-(4.15) or (4.12)-(4.16) by increasing M 1 and M 2 at each loop in order to have η = ξ 2 = 0, which guarantees that the original primal-dual SDCO problem has been solved. The form of the problems (4.12) and (4. [START_REF] Helmberg | An interior-point method for semidefinite programming[END_REF]) clearly indicates that one must increase M 1 if η > 0 and one must increase M 2 when ξ 2 > 0.

Algorithm 4.1.3 (getting a solution from a linear feasible point) One loop of the algorithm is as follows.

1. Solve (4.13)-(4.15) or (4.12)-(4.16) from a primal-dual strictly feasible point to get (x, ξ 1 , ξ 2 , y, η, s, σ 1 , σ 2 ).

2. Successful stopping test. If ξ 2 = 0 and η = 0, stop and declare that (x, y, s) is a primal-dual solution to (P ).

3. Failure on dual infeasibility. If M 1 is too large, stop and declare that it is likely that the dual problem is infeasible.

4.

Failure on primal infeasibility. If M 2 is too large, stop and declare that it is likely that the primal problem is infeasible.

5.

Increase M 1 if η > 0.

6. Increase M 2 if ξ 2 > 0.

Implementation

Here are a few recommendations.

Because of the introduction of modified SDCO problems that must be solved by the developed solver, this one has to deal with SDCO problems of various dimensions and data in the same run. It is therefore recommended to have a Matlab function that can find a solution of an SDCO problem when a strictly feasible primal-dual starting point is given (and when the problem has a solution). It is therefore recommended to write this function with the code that has been developed so far in the previous chapter 3. I could have the form

[x,y,s] = sdco_solve (A,b,c,K,x0,y0,s0)
To develop the solver, a good idea is to test it on the easy test cases 1 and 2, in which the initial values of x and y has been discarded. Once the solver works correctly on these test cases, you may want to test it on the more special problems of sections 4.3.1, 4.3.2, and 4.3.3.

The linear operator à does not depend on M 1 and M 2 and can be set outside the loop described in algorithm 4.1.3.

Check whether the linear constraints Ã(x 0 ) = b and à * (ỹ) + s = c are satisfied at the beginning of each cycle, after the setting of Ã, x0 , and b.

A practical infeasible predictor-corrector algorithm

We present in this section an infeasible interior point algorithm, whose first goal is efficiency, letting aside theoretical issues dealing with convergence and complexity. The algorithm uses a predictor-corrector technique of the Mehrotra-type [START_REF] Mehrotra | On the implementation of a primal-dual interior point method[END_REF], which provides superquadratic convergence when it is used to solve monotone linear complementarity problems [START_REF] Wright | A superquadratic infeasible-interior-point method for linear complementarity problems[END_REF]. This is the method implemented in SDPT3 [START_REF] Toh | On the implementation and usage of SDPT3a Matlab software package for semidefinite-quadratic-linear programming, version 4[END_REF].

Second order correction

Direction definition

Here is another approach to get the NT direction. Since the first two equations in the system (1.14) are linear, the difficult part of the derivation of the direction comes from the linearization/symmetrization of its last equation, which reads

x s = µe.

At the current iterate z = (x, y, s) ∈ E × R m × E, satisfying or not this identity, one would like to compute a displacement d

= (d x , d y , d s ) ∈ E × R m × E such that (x + d x ) (s + d s ) = µe.
This also reads

x d s + d x s = µ e -(x s + d x d s ).
Recall the definition of g and g T in (1.24) and (1.3). Multiplying the two sides of the previous identity to the left by g -1 and to the right by g, one gets

g -1 (x d s + d x s) g = µ e -g -1 (x s + d x d s ) g. (4.18)
Given a and b ∈ E, the symmetrization of g -1 a b g ∈ F is obtained by replacing it by

1 2 (g -1 a b g + g T b a g -T ) ∈ E,
where we have set g -T := g -T (recall the rule (1.4a) and that the matrix parts of a vector in E are symmetric). Applying this technique on the previous identity yields

E(d x ) + F(d s ) = µ e -S(x s + d x d s ), (4.19)
where the linear operators E : E → E and F : E → E, and the symmetrization linear operator S : F → E are defined at d x and d s ∈ E, and at z ∈ F by

E(d x ) = 1 2 g -1 d x s g + g T s d x g -T , F(d s ) = 1 2 g -1 x d s g + g T d s x g -T , S(z) = 1 2 g -1 z g + g T z g -T .
Finally, the linearization is obtained by dropping the single nonlinear term d x d s from the right-hand side of (4.19), leading to

E(d x ) + F(d s ) = µ e -S(x s), (4.20)
which is identical to the third equation in the NT system (1.22).

The goal of the so-called second order correction is to define a better direction than the NT direction (1.22) by reconsidering the linearization phase of the last paragraph. Instead of dropping the term d x d s , it is replaced by the one obtained by a prediction phase, say p = (p x , p y , p s ) obtained by setting µ = 0 in (1.22). Hence, (4.20) becomes

E(d x ) + F(d s ) = µe -S(x s + p x p s ). (4.21)

Direction computation

Let us see how to compute the solution to (4.21). Using the notation dx := g -1 d x g -T and ds := g T d s g, the identities

g T s g = g -1 x g -T = σ and g -1 x s g = σ 2 (4.22)
recalled from (1.26), and the fact that σ has its matrix part formed of positive definite diagonal matrices, this identity (4.21) becomes

1 2 dx + ds σ + σ dx + ds = µ e -σ 2 - 1 2 g -1 p x p s g + g T p s p x g -T .
(4.23) This is a Lyapunov equation in the unknown dx + ds , which is particularly easy to solve, since σ has diagonal matrix parts. Taking the (i, j) element of a matrix part in both sides gives

1 2 (Σ ii + Σ jj )( dx + ds ) ij = (µ -Σ 2 ii )δ ij - 1 2 g -1 p x p s g + g T p s p x g -T ij , where δ ij is the Kronecker symbol (δ ij = 1 is i = j, δ ij = 0 if i = j). As a result, dx + ds = µ σ -1 -σ + γ, (4.24)
where the (i, j) element of a matrix part of γ is

γ s ij := -(g s ) -1 p s x p s s g s + (g s ) T p s s p s x (g s ) -T ij /(Σ ii + Σ jj ). (4.25a)
For the vector part of γ, one immediately get from the Lyapunov equation (4.23):

γ l := -(σ l ) -1 • p l x • p l s . (4.25b)
Left and right-multiplying the two sides of the identity (4.24) by g and g T respectively, using w = g g T observed in (1.25), and using (4.22), yield

d x + w d s w = µs -1 -x + g γ g T . (4.26)
We see that we have just added the correction g γ g T to the right-hand side of (1.22) 3 .

To summarize, the NT direction with second order correction is the solution to the system

   A * (d y ) + d s = r d A(d x ) = r p d x + w d s w = r s , (4.27a)
where the residuals r d and r p have been defined by (1.16), and the residual r s is defined by

r s := µs -1 -x + g γ g T . (4.27b)

Stepsize computation

Most implemented algorithms do not use a neighborhood of the central path to control the iterates (this is the case for the algorithms implemented in SDPA [START_REF] Fukuda | Latest developments in the SDPA family for solving large-scale SDPs[END_REF] and SDPT3 [START_REF] Toh | On the implementation and usage of SDPT3a Matlab software package for semidefinite-quadratic-linear programming, version 4[END_REF]). There are at least three reasons for this.

When the iterates are not feasible, the feasible central path cannot be used as a safeguard, so that the neighborhood must take into account infeasibility and are defined with more complexity. It is more easy to determine a stepsize to the boundary of the cone of positive semidefinite matrices than to the boundary of the neighborhood of the central path.

For efficiency reasons, implemented algorithms compute distinct stepsizes for the primal and dual variables. These stepsizes α p and α d are computed from the maximal stepsizes ᾱp and ᾱd maintaining the positive semidefiniteness of matrix parts of x and s ∈ E.

One proceeds as follows. Consider the case of a matrix variable M , to which corresponds the displacement d M (M can be a matrix component of the primal or dual variables x and s ∈ E). The desired stepsize α M from M along d M is determined as large as possible in (0, 1], while ensuring a "sufficient" positive definiteness of M + α M d M . Since M ≻ 0, it makes sense to take

α M = 1, if d M 0. (4.28)
If d M is not positive semidefinite, one looks for the largest α M ∈ (0, 1], with some safeguard, in order to satisfy the equivalent conditions

λ min (M + α M d M ) 0, (4.29a)
λ min (I + α M M -1 d M ) 0, (4.29b)
λ min (I + α M L -1 M P T M d M P M L -T M ) 0. (4.29c)
The equivalence between (4.29a) and (4.29b) comes from the fact that M + α M d M and I + α M M -1 d M have the same eigenvalues, which are therefore in R, although the latter matrix is not symmetric [START_REF] Fujisawa | Numerical evaluation of SDPA (semidefinite programming algorithm)[END_REF]. In (4.29c), L M is a factor of M in the sense that

P T M M P M = L M L T M
for some possible permutation matrix P M (it could be the Cholesky factor with permutation matrix P M ). These conditions can be rewritten as follows

α M λ min (M -1 d M ) -1, (4.30a
)

α M λ min,M (d M ) -1, (4.30b
)

α M λ min (L -1 M P T M d M P M L -T M ) -1, (4.30c)
where we have written λ min,B (A) the minimal λ such that M v = λBv for some nonzero v (the generalized eigenvalues of Av = λBv are real when both A and B are symmetric and B is positive definite, like here). Since, in the present case, d M | 0, hence the minimal eigenvalues above are negative, the previous inequalities, with some safeguard β ∈ (0, 1) becomes respectively

α M = min 1, -β λ min (M -1 d M ) , if d M | 0, (4.31a
)

α M = min 1, -β λ min,M (d M ) , if d M | 0, (4.31b
)

α M = min 1, -β λ min (L -1 M P T M d M P M L -T M ) , if d M | 0. (4.31c)
In the case of the NT direction, the Cholesky factor L M of P T M M P M is computed, so that the use of λ min (L -1

M P T M d M P M L -T M )
is advantageous, in order to have the symmetry of the matrix whose minimum eigenvalue has to be computed, while the computing time for computing the two matrices L -1 M P T M d M P M L -T M and M -1 d M is approximately the same. In Matlab, with eig or eigs, formula (4.31c) seems to be the one that is the most conservative, in the sense that it gives the most negative minimal eigenvalue; for that reason, it is recommended.

The algorithm

The modifications come from the fact that the centering parameter σ is determined by a prediction step (i.e., with µ = 0 in (1.22)), before the computation of the search direction d. Furthermore, the computation of the direction is based on a modification of the residual in the third equatio of (1.22), which is now given by r s in (4.27b). Algorithm 4.2.1 (IPC) One describes one iteration from z = (x, y, s) to z + = (x + , y + , s + ). A stepsize safeguard β ∈ (0, 1) is also updated in the loop to become β + (take β = 0.9 at the first iteration).

1. Stopping criterion. Stop if the feasibility and µ(z) are small enough.

2. Prediction step and its stepsizes. β := 0.9 + 0.09 min(α ′ p , α ′ d ).

6. Update of the variables:

x + = x + α p d x , y + = y + α d d y , s + = s + α d d s .
7. Update of β: β + := 0.9 + 0.09 min(α p , α d ).

The centering parameter σ ∈ (0, 1) set in step 3 is used in step 4 to determine the degree to which the NT direction must aim a central point close to a solution of the problem. If σ is close to 1 a cautious direction is computed, which favors a move to the central path instead of one to the solution. In contrast, a σ close to zero favors an audacious move towards the solution. Assuming that the exponent e = 1, we see that if the prediction step computes a point z ′ providing an important decrease of the duality gap µ(z ′ ), then the centering parameter is taking small in step 3. The exponent e in (4.32) has been introduced to reinforce this behavior of the algorithm. In SDPT3 [START_REF] Todd | SDPT3 -A MATLAB software package for semidefinite programming, version 1.3[END_REF], it has been determined experimentally and set to

e := max[1, 3 min(α ′ p , α ′ d ) 2 ] if µ(z) > 10 -6 , 1
otherwise.

Test cases

Test case 4a

Consider the problem [50; p. 65]:

       sup y 1   0 y 1 0 y 1 y 2 0 0 0 y 1 + 1   0.

Test case 4b

Consider the primal SDO problem written in standard form with the following data:

C = 0 1 1 0 , A 1 = -1 0 0 0 , A 2 = 0 0 0 -1
, and b = -1 0 .

Test case 4c

Consider the primal SDO problem written in standard form with the following data: This chapter describes a few non trivial applications of self-dual conic optimization, which can be used as test-cases for the developed SDCO solver.

C = 0 0 0 0 , A 1 = 1 0 0 0 , A 2 =

Global minimization of a univariate polynomial

Unconstrained minimization

We consider the problem of finding the global minimum of a univariate real polynomial p ∈ R[x] (hence in the variable x ∈ R), which reads

inf x∈R p(x) := α∈N p α x α , (5.1)
where only a finite number of coefficients p α ∈ R are nonzero. The degree d := deg p of the polynomial p is the largest α ∈ N such that p α = 0. Problem (5.1) can be rewritten as an SDCO problem. Its primal form can be deduced from the SOS approach, while its dual form can be obtained by the moment approach. We start with the SOS approach, which is more intuitive for a reader who is not familiar with measure theory.

The SOS approach

The name of this approach comes from the fact that it is grounded on the representation of a nonnegative univariate polynomial as a sum of squares (SOS) of polynomials (proposition 5.1.1).

The starting point consists in reformulating (5.1) as a problem with an infinite number of constraints and a single variable s ∈ R:

sup s∈R s p(x) s, ∀ x ∈ R.
(5.2)

This straightforward reformulation results in an amazing observation: finding the global minimum of a function (not necessarily a polynomial) is a linear optimization problem, but with an infinite number of linear constraints. The question now is to determine whether one can transform this problem in an equivalent one, but with a finite number of constraints. It happens that this is possible when p is a polynomial, essentially because the equality between two polynomials can be expressed by the equality between their coefficients, whose number is finite. Let us introduce the approach progressively. The infinite number of constraints of the previous problem "disappears" if we trivially express them in terms of membership to the cone P of univariate nonnegative polynomials:

sup s∈R s p -s ∈ P.

(5.

3)

It turns out that the membership to P has an LMI representation (this is no longer true for multivariate polynomials, but there are bypasses [START_REF] Parrilo | On a decomposition for multivariable forms via LMI methods[END_REF][START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Lasserre | Moments Positive Polynomials and Their Applications[END_REF][START_REF] Anjos | Introduction to semidefinite, conic and polynomial optimization[END_REF][START_REF] Blekherman | Semidefinite Optimization and Convex Algebraic Geometry[END_REF][START_REF] Lasserre | An Introduction to Polynomial and Semi-Algebraic Optimization[END_REF]). That fact is based on the following two properties: a univariate nonnegative polynomial can be written as a sum of squares (SOS) of polynomials (proposition 5.1.1, which is no longer true for multivariate polynomials) and an SOS polynomial can be represented thanks to a positive semidefinite matrix (proposition 5.1.2, which is still true for multivariate polynomials).

Proposition 5.1.1 (nonnegative univariate polynomial on R) A univariate polynomial p ∈ R[x] is nonnegative on R if and only if it is of even degree, say 2m, and reads p = q 2 + r 2 , with q, r ∈ R[x], deg q = m, and deg r m -1.

Proof. The given conditions are clearly sufficient. Let us show that they are necessary. Being nonnegative on R the polynomial is necessary of even degree, say 2m, and the leading coefficient is positive. There is therefore no loss of generality in supposing that this leading coefficient is 1. Then, the polynomial can be decomposed in m factors of the form (x -a) 2 + b 2 .

It is indeed the form of (x -r)(xr) when r and r are complex conjugate roots a ± ib.

On the other hand, any real root has an even multiplicity (otherwise the polynomial would have positive and negative values around the root) and each double root is of the form above with b = 0. The m factors of the form q 2 + r 2 are then multiplied successively, by using the following formula (q 2 j + r 2 j )(q 2 + r 2 ) = (q j q + r j r) 2 + (q j r -r j q) 2 =: q 2 j+1 + r 2 j+1 .

By induction, we see that deg q j = j and deg r j j -1. It is indeed the case for j = 1 since deg q 1 = 1 and deg r 1 = 0. Now, be r vanishing or not, deg q j+1 = deg q j +1 = j+1. Finally, if r = 0, deg r j+1 = deg r j + 1 j and, if r = 0, deg r j+1 max(deg q j , deg r j + 1) j. ✷ Proposition 5.1.2 (caracterization of SOS polynomials) A polynomial p ∈ R[x] of degree 2m is a sum of r squares of polynomials if and only if there exists a matrix X 0 of order m + 1 and of rank r such that p(x) = v m (x) T Xv m (x), where v m (x) is the vector of monomials (1 x x 2 . . . x m ) T .

Proof.

[⇒] If p ∈ R[x] is of degree 2m and reads r i=1 σ 2 i , where σ i ∈ R[x], the degrees deg σ i m. Therefore, on can find vectors z i ∈ R m+1 such that p(x) = r i=1 (z T i v m (x)) 2 = r i=1 v m (x) T z i z T i v m (x) = v m (x) T Xv m (x),
where

X := r i=1 z i z T i 0 is of rank r. [⇐] Conversely, suppose that p(x) = v m (x) T Xv m (x), with X 0 of rank r. The spectral decomposition of X = r i=1 z i z T i allows us to write p(x) = r i=1 v m (x) T z i z T i v m (x) = r i=1 (z T i v m (x)) 2 ,
showing that p is the sum of the squares of at most r polynomials. ✷

Using the above propositions and setting s = -t, problem (5.3) can be written

   inf (X,t)∈S m+1 ×R t p(x) + t = v m (x) T Xv m (x), ∀ x ∈ R X 0.
(5.4) Knowing p, the equality constraint in (5.4) is an affine constraint on the unknown t ∈ R and the unknown real elements of X, since it is equivalent to the equality between the coefficients of the same monomials in both sides of the identity. To highlight the structure of problem (5.4) and make the determination its dual easier, for α ∈ [0 : 2m], let us introduce the symmetric matrix B α of order m + 1, whose (i, j) ∈ [1 : m + 1] 2 element is defined by

(B α ) ij := 1 if i + j -2 = α 0 otherwise. (5.5) Then, v m (x) T Xv m (x) = i,j∈[1 : m+1] X ij x i+j-2 = α∈[0 : 2m]   i+j-2=α X ij   x α = α∈[0 : 2m] B α , X x α .
Finally, problem (5.4) reads

   inf (X,t)∈S m+1 ×R t B α , X = p α + tδ 0α , ∀ α ∈ [0 : 2m] X 0, (5.6)
where δ 0α = 1 if α = 0 and δ 0α = 0 otherwise. Problem (5.6) is almost in the primal form of an SDCO problem. The difference comes from the free variable t, which is only constrained by the affine constraint in (5.6), not by a nontrivial cone. At least two techniques could be considered to solve this problem, which is not in the standard form of the primal SDCO problem.

(F1) Introduce additional free variables x f ∈ R p in the SDCO model (1.7) (by free variables, we mean variables that are forced to satisfy the affine constraint but not to belong to an additional cone, hence they are not completely free).

(F2) Write t = u -v and impose u 0, v 0. The standard SDCO model (1.7) can then be used with s = 1,

n 1 = m + 1, n l = 2, x =   X u v   , c =   0 S m+1 1 -1   , (5.7)
where 0 S m+1 is the zero matrix of order m + 1, and A is the map linking linearly the coefficients of p(x) + u -v and those of v m (x) T Xv m (x).

We suggest using the second approach, which has the advantage of allowing us to use the so far developed code without modification. Actually, some authors [START_REF] Toh | On the implementation and usage of SDPT3a Matlab software package for semidefinite-quadratic-linear programming, version 4[END_REF] decompose each free variable into the difference of two nonnegative variables, as we did it for t above, since otherwise their numerical treatment seems to introduce difficulties. The final primal formulation of the problem is therefore

   inf (X,u,v)∈S m+1 ×R×R u -v B α , X -(u -v)δ 0α = p α , ∀ α ∈ [0 : 2m] X 0, u 0, v 0, (5.8)
where δ 0α = 1 if α = 0 and δ 0α = 0 otherwise.

As a result, for finding the global minimum value of a polynomial of degree d = 2m (the degree must be even, otherwise the polynomial is not bounded below), the SDCO approach requires determining a positive semidefinite matrix of order m + 1 and a nonnegative vector of size 2, subject to d + 1 affine constraints. For this problem, the linear operator A of the SDCO problem (1.7) has actually a very sparse matrix representation; more specifically, the proportion of nonzero elements is

m 2 + 2m + 1 2m 3 + 5m 2 + 7m + 3 ∼ 1 2m = 1 deg p (for large m),
so that using sparsity techniques allows the approach to compute the global minimal value of polynomials of rather high degree.

The moment approach

Whilst the SOS approach provides the minimal value of a polynomial, the moment approach provides the minimizer of the polynomial. Both approaches are dual to each other in the sense of Lagrange duality, so that an SDCO solver finds a solution to both problems simultaneously. The dual of problem (5.8) is the following problem

       sup (y,S)∈R d+1 ×S m+1 α∈[0 : d] p α y α α∈[0 : d] y α B α + S = 0 S m+1 y 0 = -1 S 0.
(5.9)

Noting ȳα := -y α , the problem also reads The condition ȳ0 = 1 in (5.10) reflects the fact that the sought measure is forced to have a unit mass, i.e., R dµ = 1. Hence, implicitly, problem (5.10) determines the measure µ through its moments ȳα . When the minimizer of the polynomial p is unique, the optimal measure µ is the Dirac measure at the solution so that the minimizer can be recovered in ȳ1 = R x dµ.

   inf (ȳ,S)∈R d+1 ×S m+1 α∈[0 : d] p α ȳα α∈[0 : d] ȳα B α 0 ȳ0 = 1.

Constrained minimization

Assume now that the problem consists in minimizing the polynomial p in (5.11)

Like in section 5.1.1, we introduce the methods by the SOS approach. By the same reasoning as the one leading to (5.2), we see that problem (5.11) can obviously be expressed by sup s∈R s p -s is nonnegative on I.

(5.12)

Like with (5.2), this is a linear optimization problem with an infinite number of constraints, one constraint p(x) s on s for each x in I. From this point of view, we are face to the problem of expressing the nonnegativity of a univariate polynomial on an interval by a finite number of constraints. Remarkably, it turns out that this property can also be expressed by an LMI. One can give a representation of a univariate polynomial that is nonnegative on the interval [0, ∞) without sophisticated algebraic concepts, see proposition below 5.1.3 [START_REF] Pólya | Problems and Theorems in Analysis II[END_REF][START_REF] Powers | Polynomials that are positive on an interval[END_REF]. The case of intervals of the form [a, ∞), (-∞, b], or [a, b], with a and b ∈ R, are then dealt with by a change of variable. A common interest of the obtained representations of univariate polynomials is to be of the form

p = j∈J σ i q j , (5.13)
where |J| 2, the σ j 's are SOS with a degree that can be determined from the one of p, and the q j 's are polynomials of degree 0, 1, or 2 (the constraints q j (x) 0 represent the considered interval). This property is important from a numerical point of view, since the degrees of the σ j 's directly determine the order of the positive semidefinite matrices representing them (proposition 5.1.2).

Nonnegativity on [a, ∞)

We start with proposition 5. Proof.

[⇒] Observe first that if p i ∈ R[x], i ∈ {1, 2}, are two polynomials of the required form σ 0,i + σ 1,i x with σ 0,i and σ 1,i ∈ Σ[x], deg σ 0,i deg p i and deg σ 1,i deg p i -1, then the same property holds for their product, which reads p 1 p 2 = (σ 0,1 σ 0,2 + σ 1,1 σ 1,2 x 2 ) + (σ 0,1 σ 1,2 + σ 0,2 σ 1,1 )x.

Indeed, σ 0,1 σ 0,2 + σ 1,1 σ 1,2 x 2 ∈ Σ[x], with a degree deg p 1 p 2 and σ 0,1 σ 1,2 + σ 0,2 σ 1,1 ∈ Σ[x] with a degree deg p 1 p 2 -1. Now, write p like the product of its factors. The product of two factors corresponding to conjugate roots a ± ib reads (x -a) 2 + b 2 ∈ Σ[x], which is of the required form (with σ 1 = 0). The positive real roots must have an even multiplicity and the product of the corresponding factors is therefore in Σ[x], which is of the required form (with σ 1 = 0). To each nonopsivite real root corresponds a factor of the form a + x, with a 0, which is of the required form (with σ 0 = a and σ 1 = 1).

where each term has a degree d. This problem covers a large family of problems, including algebraic optimization, whose problems have polynomial objective and constraints. It also includes problems with binary constraints since x j ∈ {0, 1} can be expressed by the polynomial constraint x j (x j -1) = 0. Whilst a QCQO problem is NP-hard, its optimal value can be approached from below by a semidefinite optimization problem.

QCQO formulation of an OPF problem

One of the emblematic instances of the so-called Optimal Power Flow (OPF) problem, in alternating current, consists in minimizing the Joule heating losses in an electricity transmission network, while satisfying the electricity demand, by determining appropriately the powers of the generators installed at given buses (or nodes of the network) [START_REF] Carpentier | Contribution à l'étude du dispatching économique[END_REF][START_REF] Capitanescu | State-of-the-art, challenges, and future trends in security-constrained optimal power flow[END_REF][START_REF] Capitanescu | Critical review of recent advances and further developments needed in AC optimal power flow[END_REF]. This optimization problem is NP hard, implying that there is (presently) no algorithm to find the (global) solution in polynomial time (i.e., in a reasonable time when the network size is like the one in most countries or union of countries), while there are several good reasons for being only interested in the global solutions.

Structurally, one can view (or reduce) the problem to a nonconvex quadratic optimization problem with nonconvex quadratic constraints (sometimes abbreviated by the acronym QCQO for Quadratically Constrained Quadratic Programming). Recalling that any polynomial optimization problem can be written that way [START_REF] Shor | Class of global minimum bounds of polynomial functions[END_REF][START_REF] Ferrier | Hilbert's 17th problem and best dual bounds in quadratic minimization[END_REF][START_REF] Madani | Rank-2 matrix solution for semidefinite relaxations of arbitrary polynomial optimization problems[END_REF], one understands the potential difficulty of such a problem.

A QCQO representation of an OPF instance can be obtained by the Matlab function qcqp_opf [START_REF] Josz | The Matlab function qcqp_opf[END_REF], available in Matpower [START_REF] Zimmerman | MATPOWER steady-state operations, planning and analysis tools for power systems research and education[END_REF] (a Matlab package for simulating power systems). This function writes a specified OPF instance as a QCQO problem in complex numbers:

   inf z z H Cz z H A k z = a k , for k ∈ [1 : m] z H B k z b k , for k ∈ [1 : p], (5.24)
where z ∈ C n is a complex vector, C, A k , and B k ∈ C n×n are Hermitian matrices1 , a and b are real vectors, and the exponent • H is used to denote the conjugate transpose. In this representation, n is the number of buses (or nodes of the network).

To solve (5.24) (more precisely a relaxation of it) by an SDCO solver in real numbers, the first task is to rewrite the problem in real numbers (see the discussion in [START_REF] Gilbert | Plea for a semidefinite optimization solver in complex numbers[END_REF] on the interest in having an all-complex approach). For M ∈ C n×n , we write M = ℜ(M ) + iℑ(M ), where ℜ(M ) and ℑ(M ) ∈ R n×n , and i ∈ C the pure imaginary number (i 2 = -1). Similarly, a vector z ∈ C n is decomposed in z = ℜ(z) + iℑ(z), with ℜ(z) and ℑ(z) ∈ R n . It is easy to see that where testcase is a string giving the name of a Matpower test case (see the table below for a short list), A, b, and c stand for the data A, b, and c of the standard primal SDCO problem in (1.7), while x0 and y0 are empty variables (no initial primal-dual point is provided).

M is Hermitian ⇐⇒ ℜ(M ) is symmetric, ℑ (
The rank relaxation of a QCQO can be viewed as the Lasserre [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF] (or moment-sos) relaxation of degree one [START_REF] Lasserre | Convergent LMI relaxations for nonconvex quadratic programs[END_REF]. For the OPF problem, when the degree one relaxation of its QCQO formulation is inexact (i.e., it does not provide the solution to the original QCQO), it is often the case that the degree 2 moment-sos relaxation is exact [START_REF] Josz | Application of the moment-SOS approach to global optimization of the OPF problem[END_REF].

Notes

The proof of proposition 5.1.1 is taken from [START_REF] Pólya | Problems and Theorems in Analysis II[END_REF]§ VI 44]. The proof of proposition 5.1.3 is taken from [39; 2000 ; proposition 2] and is a rewriting of the proof of Pólya and Szegö [38;[START_REF] Robinson | Stability theory for systems of inequalities, part II: differentiable nonlinear systems[END_REF]§ VI 45]. For more advanced technique to minimize polynomials subject to polynomial constraints (the algebraic optimization problem), see [START_REF] Parrilo | On a decomposition for multivariable forms via LMI methods[END_REF][START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Lasserre | Moments Positive Polynomials and Their Applications[END_REF][START_REF] Anjos | Introduction to semidefinite, conic and polynomial optimization[END_REF][START_REF] Blekherman | Semidefinite Optimization and Convex Algebraic Geometry[END_REF][START_REF] Lasserre | An Introduction to Polynomial and Semi-Algebraic Optimization[END_REF].

3

 3 The HKM direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Data representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Calling statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Matlab functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.4 Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Test cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Test case 1a: an easy SDO problem of dimension 3 . . . . . . . . . . . . . . 1.4.2 Test case 1b: a simple LO problem of dimension 2 . . . . . . . . . . . . . . 1.4.3 Test case 1c: two SDO problems and one LO problem in parallel . . Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Moving along the central path . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Algorithmic concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 The closest central point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Neighborhood of the central path . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Two algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 A predictor-corrector algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 A large step algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Test case 2: minimum matrix norm . . . . . . . . . . . . . . . . . . . . . . . . . . . Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . This first session introduces the self-dual conic optimization (SDCO) problem to solve.

1 .

 1 Directions of displacement •, • S n j : (a, b) ∈ S n j × S n j → a, b S n j := tr ab = n j k,l=1a kl b kl ,

1 .

 1 Directions of displacement (v + dx ) (v + ds ) = µe, where dx = w -1/2 d x w -1/2 and ds = w 1/2 d s w 1/2 .

( 4 . 5 )

 45 This problem is identical to problem (4.4) provided the vectors b and ỹ ∈ R m+1 , the adjoint linear map à * : R m+1 → Ẽ, and the vectors s and c ∈ Ẽ are defined by

( 4 . 10 )

 410 Now, à * : R m → Ẽ is defined at x byà * (y) =A * (y) -A(e) T y .

Algorithm 4 .

 4 1.2 (getting a solution from a strictly feasible dual point) Let (y, s) ∈ F s D and let M 2 > e, s E . One iteration of the algorithm is as follows. 1. Solve (4.8)-(4.11) from a primal-dual strictly feasible point to get (x, ξ 2 , y, s).

4 . 5 .

 45 Compute the prediction direction d ′ := (d ′x , d ′ y , d ′ s ) by(1.22) with µ = 0. Compute the primal and dual stepsizes along that direction d ′ by (4.31c), for the current stepsize safeguard β, and denote them by α ′ p and α ′ d (there is one stepsize per matrix component and one for each element of the vector component). Setz ′ := (x + α ′ p d ′ x , y + α ′ d d ′ y , s + α ′ d d ′ s ).3. Centering parameter. For some exponent e discussed below, set the centering parameter to Direction. Compute the search direction d := (d x , d y , d s ) by (4.27) with p = d ′ and µ set to µ := σµ(z). (4.33) Stepsizes. Compute the primal and dual stepsizes α p and α d by (4.31c), with

( 5 . 10 )

 510 This problem has an interpretation in terms of the moments ȳα of a nonnegative measure µ on R, defined by ȳα = R x α dµ.

  (5.1) on an interval I ⊂ R, which may take the form I = R (unconstrained problem), I = [a, ∞), I = (-∞, b], or I = [a, b], for some a and b ∈ R: inf x∈I p(x).

1 . 3 ,

 13 which gives the Pólya and Szegö representation of a univariate polynomial that is nonnegative on R + [38; 1976 ; § VI 45]. Proposition 5.1.3 (nonnegative univariate polynomial on [0, ∞)) A univariate polynomial p ∈ R[x] of degree d is nonnegative on [0, ∞) if and only if there exist σ 0 and σ 1 ∈ Σ[x] such that deg σ 0 d, deg σ 1 d -1 and, for all x ∈ R: p(x) = σ 0 (x) + σ 1 (x)x.(5.14)

Proof. ✷ 5 . 1 . 3

 513 Support toolThe Matlab function [A,b,c,K] = polmin2sdco (p,a,b); transforms in the SDCO format (A, b, c, K) the problem of finding the global minimum of the following 1D polynomial of degree d d+1 i=1 p i x d+1-i , on the inteval [a, b] ⊂ R, along the lines explained in sections 5.1.1 and 5.1.2. In this expression, p i ∈ R is the ith element of the Matlab vector p, which has length d + 1.The function purges p from its first zero elements and makes sure that the minimzation problem is not unbounded, by examining the sign of p 1 , the parity of d, and the possible infinite value of a and b.5.2 Rank relaxation of a QCQO problemA quadratically constrained quadratic optimization (QCQO) problem has the following structure   inf x x T A 0 x x T A k x = b k , ∀ k ∈ E, x T A k x b k , ∀ k ∈ I,where the matrices A k ∈ S n (for k ∈ [0 : m]), the scalars b k ∈ R (for k ∈ [1 : m]), and the index sets E and I form a partition of [1 : m]. Since the matrices are not assumed to be positive semidefinite, the problem is nonconvex. It is actually NP-hard.

  M ) is skew symmetric.With this notation and a Hermitian matrix M , there holdsz H M z = zT M z, ) -ℑ(M ) ℑ(M ) ℜ(M ) ∈ S 2n .Note indeed that M is symmetric. Then, problem (5.24) becomes the QCQO in real numbers[A,b,c,K,x0,y0] = testcase_5 (testcase);

Questions 5 . 1 .

 51 Global minimization of a univariate polynomial. Find the global minimum of a few univariate polynomials p on intervals I ⊂ R, using the moment-SOS technique presented in sections 5.1.1 and 5.1.2, with the help of the Matlab function polmin2sdco described in section 5.1.3:1) p(x) = x 6 -7x 5 + 7x 4 + 35x 3 -56x 2 -28x + 48, I = R [43], 2) p(x) = x 3 + 3x 2 -9x, I = [-6, ∞) and I = [-3, ∞).3) p(x) = -x 3 + 3x 2 + 9x, I = (-∞, 6] and I = (-∞, 3].

Table of Contents

 of 1 Directions of displacement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 The problem to solve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Space structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 The primal and dual SDCO problems . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.3 The central path . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Two directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 The NT direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.

1.2. Two directions

Directions of displacement

Moving along the central path

Finding an appropriate starting point

Starting from an infeasible point

Recall that M ∈ C n×n is Hermitian if M H = M .

[⇐] This is clear since for x ∈ [0, ∞), x 0, while σ 0 and σ 1 are nonnegative polynomials. ✷ we have the following table

Note also that the degrees of σ 0 and σ 1 in corollary 5.1.4 are respectively 2m 0 and 2m 1 . Using proposition 5.1.1, the representation (5.15) of p, and setting s = -t, problem (5.12) can be written

(5.17) Knowing p, the equality constraint in (5.17) is an affine constraint on the unknown t ∈ R and the unknown real elements of X 0 and X 1 , since it is equivalent to the equality between the coefficients of the same monomials in both sides of the identity. Let us introduce the matrices B 0 α ∈ S m 0 +1 and B 1 α ∈ S m 1 +1 by the formula (5.5). Then, one has

Set t := u -v with u and v 0. Then, problem (5.17) becomes

To summarize, for finding the global minimum value of a polynomial of degree d, the SOS approach requires determining two positive semidefinite matrices of order m 0 + 1 and m 1 + 1, as well as a nonnegative vector of size 2, subject to d + 2 affine constraints.

Nonnegativity on (-∞, b]

Proof. Let p be the polynomial defined at x ∈ R by p(x) = p(b -x). Then, p 0 on [0, ∞). By proposition 5.1.3, there exist σ0 and σ1 ∈ Σ[x], with deg σ0 d and ) and use the table that follows. Note that the degrees of σ 0 and σ 1 in corollary 5.1.5 are, like in the previous section, 2m 0 and 2m 1 respectively. Using proposition 5.1.1, the representation (5.19) of p, and setting s = -t, problem (5.12) can be written

(5.20) Knowing p, the equality constraint in (5.20) is an affine constraint on the unknown t ∈ R and the unknown real elements of X 0 and X 1 , since it is equivalent to the equality between the coefficients of the same monomials in both sides of the identity. Let us introduce the matrices B 0 α ∈ S m 0 +1 and B 1 α ∈ S m 1 +1 by the formula (5.5). Then, one has

Set t := u -v with u and v 0. Then, problem (5.20) becomes

is odd X 0 0, X 1 0, u 0, and v 0.

(

To summarize, for finding the global minimum value of a polynomial of degree d, the SOS approach requires determining two positive semidefinite matrices of order m 0 + 1 and m 1 + 1, as well as a nonnegative vector of size 2, subject to d + 2 affine constraints. (ii) there exist σ 0 and σ 1 ∈ Σ[x] such that, for all x ∈ R, one has

where each term has a degree 2⌈d/2⌉, (iii) if d is odd, there exist σ 1 and σ 2 ∈ Σ[x] such that, for all x ∈ R, one has

(5.25)

The dimension of this problem is twice the one of (5.24).

Rank relaxation

The second task is to define an SDCO relaxation of (5.25). Recalling that •, • S 2n denotes the scalar product on the space of symmetric matrices (see section 1.1.1), problem (5.25) reads

Ãk , z zT

The rank constraint of this formulation is very annoying (it takes integer values, hence is discontinuous). The rank relaxation [START_REF] Vandenberghe | Semidefinite programming[END_REF][START_REF] Nesterov | Semidefinite relaxation and nonconvex quadratic optimization[END_REF][START_REF] Nesterov | Squared functional systems and optimization problems[END_REF][START_REF] Nesterov | Semidefinite programming relaxations of nonconvex quadratic optimization[END_REF][START_REF] Lasserre | Convergent LMI relaxations for nonconvex quadratic programs[END_REF][START_REF] Bomze | Solving standard quadratic optimization problems via linear, semidefinite and copositive programming[END_REF] of the problem drops that constraint and therefore reads

We still have to write the relaxed problem in the standard primal SDCO form in (1.7), which looks like the closest to the previous one. This form can be obtained by introducing p slack variables v k ∈ R, for k ∈ [1 : p], and by writing the problem as follows

which is in the primal form (1.7), in the unknown ( X, v).

Test cases

The real SDCO relaxation of the OPF problem, namely the standard primal SDCO form deduced in (5.26), is available by