Primal-dual interior-point Methods for Semidefinite Programming from an algebraic point of view, or: Using Noncommutativity for Optimization
Konrad Schrempf




Primal-dual interior-point Methods 

  for Semidefinite Programming from an

  algebraic point of view, or: Using

  Noncommutativity for Optimization
    Konrad Schrempf
  [Contact: math@versibilitas.at (Konrad Schrempf),
    <https://orcid.org/0000-0001-8509-009X>,
    Research Group ASiC, University of Applied Sciences Upper Austria,
    Ringstraße 43a, 4600 Wels;
    Faculty of Mathematics, University of Vienna,
    Oskar-Morgenstern-Platz 1, 1090 Wien; Austria.
    ]
  https://orcid.org/0000-0001-8509-009X
  
    < g r a p h i c s >

  
    December 30, 2023
============================================================================================================================================================================================================================================================================================================================================================================================================




Since more than three decades, interior-point methods proved very
useful for optimization, from linear over semidefinite
to conic (and partly beyond non-con­vex) programming;
despite the fact that already in the semidefinite case
(even when strong duality holds) “hard” problems are known.
We shade a light on a rather surprising restriction in the
non-commutative world (of semidefinite programming),
namely “commutative” paths and propose a new
family of solvers that is able to use the full richness
of “non-commutative” search directions:
(primal) feasible-interior-point methods.
Beside a detailed basic discussion,
we illustrate some variants of “non-commutative” paths
and provide a simple implementation
for further (problem specific) investigations.


Keywords and 2010 Mathematics Subject Classification.
Semidefinite programming, linear programming,
interior-point methods,
primal-dual central path, sum-of-squares;
Primary 90C22, 90C51;
Secondary 90C05



§ INTRODUCTION


Roughly speaking, semidefinite programming (SDP)
is linear programming (LP) where the (non-negative)
vectors are replaced by (positive-semidefinite) matrices.
Indeed, it is the natural generalization of LP
by dropping the restriction to diagonal matrices.
From a purely algebraic point of view, linear programming
is just a commutative version of semidefinite
programming which becomes visible in particular in
(primal-dual) interior-point methods.


This is clearly a simplified point of view since
semidefinite programming is deeply connected with
linear matrix inequalities (LMI's) 
<cit.>
. However, this would go much too far here
and is not necessary in a first reading.
For those who are familiar with LMI's:
Is there an interpretation for the ideas
(for example algebraic barriers
in Remark <ref>)
presented here?


The seed of this work is one trivial observation:
Suppose that X and Z are invertible (square) matrices,
say over the real numbers,
and let I denote the identity matrix (all of the same size).
For a scalar μ > 0 we consider the equation

XZ = μI.

Then clearly Z = μ X and since X X = X X = I we
can conclude that XZ = ZX, that is, the two matrices X and Z
commute.
Now, for diagonal matrices (in linear programs)
assuming commutativity is definitely not a restriction.
But what are the consequences for solving
general semidefinite programs (by interior-point methods)?

One “indirect” implication is discussed in 
<cit.>
, namely the need for symmetrization (of search
directions). For a more detailed discussion we recommend
<cit.>
. (Notice the ambiguity of the word “symmetrization”; here
it is meant as reformulation such that the solution is
symmetric.)
We will go a complete different way here by
“decoupling” the matrices from the primal and the dual problem.

As humus one can use some curiosity and 
the following remark of Helton and Putinar,
which might be a starting point for further (maybe problem specific)
development and in particular for robust and highly efficient implementations.
Once one realizes, why no-one writes such a solver, it is usually
too late …

“A lament is that all current computational semi-algebraic geometry
projects use a packaged semi-definite solver, none write their own. This
limits efficiencies for sum of squares computation.”
<cit.>

Those who are new to semidefinite programming and
want to grow their own plant are highly
encouraged to start with 
<cit.>
 and take the iteration scheme from
<cit.>
 as a template for an implementation.
In parallel —for a detailed introduction, a bigger context
and further references—
we recommend to have 
<cit.>
, <cit.>
, <cit.>
 and/or 
<cit.>
 at hand.
Additionally a classical book on linear programming
like <cit.>
 can be very useful.

In Section <ref> we summarize the very basics
of (numerical) linear algebra and settle a notation which
simplifies the following discussion.
Then we recall briefly linear programming in Section <ref>
before we prepare the framework for the investigation of
path-following methods for
semidefinite programming in Section <ref>
from an applied point of view (almost forgetting how they are derived).
The main contribution is the family of (primal)
feasible-interior-point methods in Section <ref>
which can follow “non-commutative” paths (meaning that
there is no assumption on the commutation of matrices).
How to find a feasible starting point is somewhat
rudimentary
sketched in Section <ref>.
The illustrations are based on a concrete problem
discussed in Section <ref> and a simple implementation
(in Octave) in Section <ref>
which is intended to serve as a starting point for
further investigations.
Section <ref> is an appetizer for semidefinite
programming especially for students.
That there is quite a lot of topics we can hardly touch
becomes clear from the literature which is not cited
here. Since this work is only the very beginning we
point out some directions for further development and
(try to) mention at least starting points to facilitate
digging into different areas in Section <ref>.


Those who expect highly sophisticated analytical tools and
investigations might be disappointed. The only prerequisites
are literacy in linear algebra, some experience
in numerics (mainly from an applied point of view)
and a tiny fraction of algebra.
For an immediate comparison of the basic idea compare
Figure <ref> (page fig:nc.sdp.103)
and Figure <ref> (page fig:nc.sdp.104).



Although we give some information on the number of iterations
in comparison with some classical solvers to stress the
—maybe surprising— simplicity of the presented concept,
we do not claim anything about usability and/or “speed”
in general. The provided code is intended for
educational purposes, in particular for
analyzing concrete optimization problems
to be able to develop more sophisticated (and optimized)
solvers, maybe in combination with existing ones.



When we use the word “path” here, it is usually meant
in a general sense since we leave the “commutative”
(analytic/continuous) central path and need to find some
other orientation in the high-dimensional vector space
of (square) matrices to be able to get to
the minimum.
Having the polyhedron from a linear program
and the path (from vertex to vertex) —generated
by an instance of the Simplex Method— on it
in mind, that generated
by feasible-interior-point methods
can take “shortcuts” through the interior
of the polyhedron (respectively spectrahedron
in the SDP case).
Since it is not clear a priori how to avoid
“detours” (or oscillations)
like those shown in Figure <ref>
(page fig:nc.sdp.202)
a lot of questions arise immediately.
Some of them are adressed in Section <ref>.




§ (NUMERICAL) LINEAR ALGEBRA


As a warm-up we fix the notation and recall
some basic linear algebra. Although there is nothing
really difficult, applied (numerical)
linear algebra can be very subtle.
But we should not worry too much since we
can rely on classical functions and algorithms
from Lapack <cit.>
 which are heavily used in 
Octave <cit.>
 and Matlab <cit.>
. For a thorough discussion of applied numerical linear algebra
we refer to <cit.>
, for linear algebra in semidefinite programming to
<cit.>
. In Section <ref> we point out
more specialized literature.

Except for the notion of “centering a matrix” —which has a
natural geometric interpretation (see for example Figure <ref>)
and is highly relevant from a numerical point of view—
there is nothing new here.
However, one should recall the difference between
the vector space of (real square) matrices ^n × n≃^n^2
and the algebra of matrices (with the non-commutative
multiplication) _n().

The set of the natural numbers is denoted by = { 1, 2, …},
that of the real numbers by . Zero entries in matrices
are usually replaced by (lower) dots. By I_n we denote the
identity matrix (of size n) respectively I if the size is clear
from the context. We use capital letters (mainly X and Z)
for matrices in _n = _n() and use
their respective lower case letters usually to denote their “vectorized”
version in ^n^2, that is,
x =  X (the first n components in x are the first
column of X, that from n+1 to 2n the second column, etc.),
or the other way around, X =  x.
There is one exception to be consistent with the classical
notation in linear programming in Section <ref>,
namely X =  x for the diagonal matrix constructed by
a vector x ∈^n × 1, that is, x_ii=x_i
for X=(x_ij). The advantage is to be able to view a
linear program (LP) as a “commutative”
semidefinite program (SDP)
since x  z =  z  x ∈_n.


At a first glance it seems a little bit strange to talk about
a matrix-matrix multiplication in linear programming (which
we recall in the next section). But since our focus is with
respect to interior-point methods (mainly discussed in Section <ref>)
our (sequence of) vector(s) x is positive, that is,
all n components of x are positive and therefore in particular
invertible. In other words: There exists a vector x̃∈^n
such that x x̃ = I_n.


We are mainly working with symmetric matrices
X = X∈_n = _n() ⊊_n
≃^n× n.
A (not necessarily symmetric) matrix X ∈_n is
called positive definite (written as X  0)
if y X y > 0 for all y ∈^n
and positive semidefinite (written as X  0)
if y X y ≥ 0 for all y ∈^n.
The trace : _n →
induces the inner product
..: _n ×_n →,

XZ := (XZ)
  = ∑_j=1^n ∑_i=1^n x_ij z_ij
    “=”  (X)Z

for X = (x_ij) and Z = (z_ij).
By X  Z we denote that X and Z are orthogonal,
that is, XZ=0.

A symmetric matrix X=X has
n (not necessarily distinct) real eigenvalues
λ_min(X) = λ_1 ≤λ_2 ≤…≤λ_n
= λ_max(X) ∈.
Our main concern is the decomposition of _n (as a vector space)
in ℓ subspaces _n^i of dimension m_i, that is,
_n = _n^1 ⊕_n^2 ⊕…⊕_n^ℓ
(with m_1+m_2+…+m_ℓ=n^2). As indices we will use mostly
i ∈{ξ, η, ζ, ν} and define
the (vectorsub-)space of non-symmetric matrices as
the orthogonal complement of the symmetric matrices,
that is,
_n^ν = (_n)^
=  M_ν
with basis M_ν∈^n^2 × m_ν
(of m_ν=n(n-1)/2 column vectors).
With a decomposition of _n we have

_n = _n^ξ⊕_n^η⊕_n^ζ_=_n ⊕_n^ν
with basis M = [M_ξ,M_η,M_ζ,M_ν]
and n^2 = m_ξ + m_η + m_ζ + m_ν.
Notice that some m_i can be zero
and recall that M_i = m_i.
By abuse of notation we write
X ∈ M_i = { M_i^(1),…,M_i^(m_ξ)}
for X =  M_i^(j) for some 1 ≤ j ≤ m_ξ.


The decomposition of the vector space of symmetric matrices
into “degrees of freedom” _n^ξ,
“minimization directions” _n^η and
“constraints” _n^ζ is used in particular
in Section <ref>,
illustrated in Figure <ref> (page fig:nc.sdp.104).
It is a reminiscence of (the graph of) a concave
function in the x-y-plane (with a possible parameter
in z-direction).


For a matrix X =(x_ij) ∈_n we denote by
X = ∑ x_ij^2 = XX
the Frobenius norm
which is just the “classical” vector norm
X=  X_2.
Very often we need to “update” X,
written as X + αΔ X with “direction”
Δ X ∈_n (“Δ X” is one symbol
and we write Δ x for “Δ X”)
and “steplength” α∈.
When X is symmetric and positive definite,
that is X=X 0 and Δ X ∈_n
we denote by

[0,∞] ∋α_max^±= α_max^±(X,ΔX)
  = sup{α∈|X ±αΔX 0 }

the maximal “steplength(s)”.
We write α_max for α_max^+.
Furthermore X=X 0 admits the Cholesky factorization
X = L L (with lower triangular L).
It can be used to compute the maximal steplength
<cit.>
:

    α_max = ( (-LΔ X L) )

which is simply α_max = 1/max (-XΔ X)
for diagonal matrices X and Δ X.
For τ < 1 we can ensure
X + τα_maxΔ X  0.
However, from a numerical point of view,
we could run into troubles if τ = 1 - for ≪ 1
since X≪X.
Now let _n^ξ =  M_n^ξ
such that X_ξΔ X
for all X_ξ∈ M_ξ
and β^± = α_max^±(X,X_ξ)
assuming 0 < β^± < ∞.
If we have the freedom to replace X by some X' = X + X_ξ
(strictly) “between”
X - β^- X_ξ and
X + β^+ X_ξ
such that (X')≪X,
we should do that. We will call that
“centering” of X (with respect to X_ξ).
We will refer to
minimizing (X+X_ξ) for X_ξ∈_n^ξ
as “algebraic” centering and
X + 12(β^+ - β^-)X_ξ
as “geometric” centering.


The “symmetrization” 1/2(X+X)
of a non-symmetric matrix X ∈_n
can be interpreted as “arithmetic” centering
(with respect to _n^ν).
Since usually there should not be that much confusion
between “purely” non-symmetric matrices and those
which have a non-trivial symmetric part, we use the
adjective “non-symmetric” in a sloppy way to simplify
the wording.
Otherwise one could alternatively use the adjective
“asymmetric”.






Let 0  X ∈_n and M_ξ a basis of
_n^ξ⊊_n
with _n^ξ = m_ξ≥ 1.
If α_max^+(X, X_ξ) and
α_max^-(X, X_ξ) are both finite
for all X_ξ∈ M_ξ then

X_cen = X_cen(M_ξ)
  = 1/2 ∑_X_ξ∈M_ξ
  ( α_max^+X_ξ- α_max^-X_ξ)

(respectively X_cen = X for m_ξ=0)
is called a geometric center of X with
respect to M_ξ.
Any X' with

(X') = min_X_ξ∈_n^ξ (X+X_ξ)

is called an algebraic center of X with
respect to _n^ξ.





Let n=2 and

X =
[   1 0.9; 0.9   1 ],
  ΔX =
[ -1  .;  . -1 ],
  X_ξ=
[ . 1; 1 . ].

Then α_max = 0.1 and X + α_maxΔ X is singular.
However

(X + (1-10^-k) α_max ΔX)≈10^k+1

while
(X + α_maxΔ X - 0.9 X_ξ)≈ 1.5713.
Now let k=10 and α = (1-10^-k)α_max.
Then β^+ = α_max^+(X,X_ξ) ≈ 10^-11
and β^- = α_max^-(X,X_ξ) ≈ 1.8.
Let β = 1/2(β^+ - β^-).
Now the geometric center is
X_cen≈ X + αΔ X + β X_ξ.



Since our main purpose is to stay away from the “boundary”
(of singular matrices) rather than finding a precise center
we did not yet make that rigorous since we have
—for practical reasons—
an orthogonal basis in mind.
One would expect that (geometric) centering is independent
of a particular basis. But this needs to be shown.


Usually we are interested in a convex combination of
a matrix X=(x_ij) and an approximate center X', that is,
μ X + (1-μ)X' for 0 ≤μ≤ 1. This can
be accomplished easily by solving a linear system
of equations (assuming X≫ 1 and X≈ 1).
By “⊗” we denote the Kronecker tensor product,
for example

X ⊗Z = 
[ x_1,1 Z x_1,2 Z       … x_1,n Z; x_2,1 Z x_2,2 Z       … x_2,n Z;       ⋮       ⋮       ⋱       ⋮; x_n,1 Z x_n,2 Z       … x_n,n Z ].

For Z  0 the matrix I ⊗ Z has full rank.
Let Z̃ = Z.
Given a basis M_ξ of _n^ξ (with m_ξ≥ 1)
and assuming Z̃/Z≫ 1,
the least squares solution of

    [ (I ⊗Z̃) M_ξ       I ⊗ Z ][   ẑ; Δz̃ ]
    = (1-μ) z̃

(using for example Lapack/DGELS)
yields the “approximate” algebraic center
Z' = Z +  M_ξx̂ with 
(Z')≪Z for μ = 0
(and Z' = Z for μ=1).



We did not investigate that in detail but it seems that
the Cholesky factorization is not as stable numerically
as computing the inverse. Therefore an “approximate”
algebraic centering can be used if a (positive definite)
matrix is almost singular
(at least if there are some degrees of freedom).



For Example <ref> we set Z = X + αΔ X
and get

M_ξẑ ≈[  0 -1; -1  0 ]
  and thus  X' = 
[  0.9 -0.1; -0.1  0.9 ]

with (X')≈ 1.6008.



How to solve (over- and) underdetermined linear systems
of equations is discussed in detail in
<cit.>
. Given a linear system A x = b,
a least squares solution x satisfies

A x - b_2 = min_x̅ ∈^n A x̅ - b_2.

For sparse matrices (and iterative methods) we refer
to <cit.>
 and their references.
It would be interesting to investigate the
“double approximate” centering compared to the
linear system (<ref>) from before
with approximated inverse (in the context of semidefinite
programming as it is suggested in Section <ref>).


A set 𝒦⊆_n() is called convex
if μ X + (1-μ) Z holds for all X,Z ∈𝒦 and 0 < μ < 1.
A linear equality A x = b with b ∈ (and A∈^1 × n)
defines a hyperplane in ^n, the corresponding
linear inequality Ax ≤ b a halfspace in ^n.
The intersection of several halfspaces given by A x ≤ b and
b ∈^m× 1 for m > 1 is again convex.
Therefore the intersection of 𝒦 and Ax≤ b is convex.




§ LINEAR PROGRAMMING


We briefly recall the basics:
Let 1≤ m<n be positive integers.
Given a (real) matrix A ∈^m × n (with rows A_1,…,A_m)
and two vectors b ∈^m× 1 and c∈^n× 1,
we want to find a vector x∈^n × 1 and
the minimum of the linear objective function
f_0 = c x such that the entries of x are non-negative
—written as x ≥ 0— and
the linear constraints f_i = A_i x = b_i for
i ∈{ 1,2,…, m } —written as linear system Ax=b—
hold. This is called the primal problem

    min_ 0 ≤ x ∈^n { c x : Ax=b }.

The corresponding dual problem is

    max_0 ≤ z ∈^n 
     y ∈^m{ b y : A y + z = c }.

For simplicity we assume that A = m and the
problems are feasible, that is,
there exists an x (respectively y and z)
such that the minimum in (<ref>) (respectively
the maximum in (<ref>)) is attained.




The “classical” solution is via Simplex Method
in finitely many steps
<cit.>
. 
Primal-dual interior-point methods move inside, that is,
x_j >0 and Ax ≤  b (respectively z_j >0 and A y + z ≤ c),
a convex polyhedron.
(See Figure <ref> for an illustration
of the “primal” path.)
Approaching the “boundary”
at x_j=0 (respectively z_j=0) before the optimum is reached
is punished using barrier functions, for example
c x - μ ( log x_1 + … + log x_n )
for the Lagrange multiplier μ > 0.
From the necessary condition

    / z_k( b y + μ∑_j=1^n log z_j )
         = / z_k xA y_=c-z + μ/z_k
       = μ/z_k - x_k != 0

for an optimum we get the
complementary condition x_j z_j = μ (for j=1,2,…,n).
Now, starting at the triple (x,y,z), we find the new search direction
(Δ x, Δ y, Δ z) by solving the
linearized system of equations

    [ . A I; A . .; z . x ][ Δ x; Δ y; Δ z ]
    =
    [         c - z - A y;             b - A x; (μ - x_j z_j)_j=1^n ]

with 2n+m unknowns. Given a 0 < τ < 1,
the steplength α≤ 1 is chosen
such that x+ταΔ x > 0 and z + ταΔ z > 0.
To make this procedure work practically,
μ must decrease (in each step) “in the right way”.

Since we will discuss the general case (for semidefinite
programming) in Section <ref> and develop
a new approach in Section <ref> we only
comment a little concerning the implementation in
Section <ref>.

Taking the classical Lagrange multiplier approach
resulting in x_j z_j = μ one can ask what would happen
if we take n independent μ_j which could
go to zero with different “speed”. In principle the
condition x_j z_j = μ > 0 just tells us that
both components, x_j and z_j, should be invertible.
Focusing on the primal problem and introducing “dummy”
variables x̃_j for the respective inverses of x_j
we need to solve the linear system

    [   A   .; c M   .;  x̃   x ][ Δ x; Δx̃ ]
    =
    [ b - A x;      -γ;       0 ]

for some γ > 0. (For the general primal-dual
formulation see (<ref>) in the
following section.) Assuming feasibility of x
and a solution Δ x of (<ref>)
there is no restriction for the steplength
as long as x + ταΔ x is positive.
Thus for a feasible (positive) starting vector
x one can use a greedy method. This is
implemented in  in Section <ref>.
Recall that we assume A = m.
By M we denote the orthogonal complement of A in
^n. A search direction can be computed by solving
the linear system

    [    c M      .; (x̃) M      x ][  x̂; Δx̃ ]
    =
    [ -γ;  0 ]

which is much smaller than (<ref>).
The search direction is Δ x = M x̂.
And if (<ref>) is solved only
approximately —say by some iterative method for sparse
matrices—, it does not have any influence
on feasibility.


A variant of the linear system (<ref>) can be used
to find an initial feasible vector x_ini > 0
by starting with x = x̃ = [1,…,1] and solving

[  A  .; x̃  x ]
[  Δx; Δx̃ ]
=
[ b - A x;       0 ]

to get a search direction Δ x.
If α = 1/λ_n(-Δ x  ( x)) > 1
we are done by using x_ini = x + Δ x.
Otherwise we set x: = x + ταΔ x and x̃ := x
(meaning componentwise inverse respectively the inverse as
diagonal matrix) for some 0 < τ < 1 and iterate.


Table <ref> shows the number of iterations
for a small linear program compared with classical
solvers. Although that is rather promising, large scale
tests needs to be done. It is not yet clear how to
prove convergence since we “lost” the
analytical setting. See also Section <ref>.
For further information on Karmarkar's algorithm
<cit.>
 in linear programming we refer to
<cit.>
.




Let b = [2,7,3], c = [-1,-2,0,0,0] and

A =
[ -2  1  1  .  .; -1  2  .  1  .;  1  .  .  .  1 ].

Then x_opt = [3,5,3,0,0].
The following output shows that only one iteration is necessary
to find a (numerically) feasible starting vector and 3 iterations
to find the minimum (up to an error smaller than 10^-9).
Notice in particular the right column with increasing norm of
the “inverse” of x due to the last two vanishing components
in x_opt.











§ SEMIDEFINITE PROGRAMMING


Semidefinite programming is the natural generalization
of linear programming by going over from (the cone of)
non-negative vectors to the cone of positive semidefinite
symmetric matrices.
The primal problem (in standard form) is

    min_0  X ∈_n{CX :
      A_1X = b_1, A_2X = b_2, …, A_mX = b_m }.

The corresponding dual problem is

    max_0  Z ∈_n,
     y ∈^m { b y : y_1 A_1 + y_2 A_2 + … + y_m A_m + Z = C }.

Using diagonal matrices X =  x, C =  c, Z =  z
and (by abuse of notation) A_k =  A_k, a linear program
can directly be formulated as semidefinite program (SDP).
For the connection with linear matrix inequalities (LMI's)
we refer to 
<cit.>
, for matrix inequalities (MI's) to 
<cit.>
.


The inner product CX is also often denoted by
C ∙ X (for symmetric matrices C and X).


For convenience we state at least the definition of the
central path and recall weak duality.
Notice that strong duality, that is,
XZ = 0, implies XZ=0
<cit.>
. Why? Does strong duality also imply X +  Z = n?
For an overview about the differences between linear and
semidefinite programming see in particular
<cit.>
. For details we refer to 
<cit.>
.


The set of solutions to the equations

    A_iX   =  b_i   for all i ∈{1,2,…, m}, 
    ∑_i=1^m y_i A_i + Z    = C   and
    
    XZ    = μ I

for all μ > 0 such that X  0 and Z  0
is called central path.




If X is feasible in (<ref>) and
(y,Z) is feasible in (<ref>), then

CX - by = XZ ≥0.




Starting from the triple (X,y,Z), a new search direction
(Δ X, Δ y, Δ Z)
can be computed by solving a linear system of equations
similar to (<ref>)
with 2n^2 +m unknowns:

    [       .       A   I_n^2;       A       .       .; I_n ⊗ Z       . I_n ⊗ X ][ Δ x; Δ y; Δ z ]
    =
    [  c - z - A y;      b - A x; (μ I_n - XZ) ].

Here we are only interested in the principal idea
of finding a (matrix-valued) search direction.
For detailed discussions and references we refer to
<cit.>
 and <cit.>
, with respect to (many subtle technical details for)
an implementation to
<cit.>
. There is a whole zoo of search directions
discussed in 
<cit.>
.

The last matrix equation of (<ref>) comes from the
linearization of the centering condition XZ = μ I
(for μ > 0) 
<cit.>
:

    (X + Δ X)(Z + Δ Z) 
         = XZ + XΔ Z + Δ X Z + Δ X Δ Z !=μ I.

This condition implies that Z is a scalar multiple of X and
therefore XZ = ZX, that is, X and Z commute.
Although this is a classical assumption and a “natural”
generalization of the complementary condition x_j z_j = μ
from (primal-dual in­te­rior-point methods for) linear programming
we have not found much discussions about its implications
except for the need for symmetrization in 
<cit.>
.
Even if strong duality holds, that is,
the duality gap is zero, 

CX_opt - by_opt = XZ_opt = 0,

it is far from clear what happens in practical
situations (where the optimum is never reached exactly).
We conjecture that the origin of hardness of SDP programs in
the sense of <cit.>
 is the restriction to “commutative” (central) paths.
While in principle this would not be any problem
for linear programming (and diagonal matrices)
one can use a much weaker assumption, namely that of
invertibility (of the components), leading to
a new class of interior-point methods.
For a short discussion in the context of LP see
Section <ref>.

When we go over from the analytic formulation (in particular
of some barrier functions 
<cit.>
) to the discrete (algebraic) —not necessarily linear—
formulation, some information is “lost”.
What we use in fact in interior-point methods
(for convex problems)
is nothing more than the possibility to find “good”
search directions
for minimization (respectively maximization) —often by solving
linear systems of equations—
and the maximal steplength
for a given search direction such that we do not
leave the “interior”.
For semidefinite programming, both steps are rather simple
(at least if we ignore fast convergence issues),
the tricky part is numerics (and the implementation).

Now we assume that X is positive definite,
in particular invertible, that is,
there exists an invertible X̃ such that
XX̃ = X̃X = I
<cit.>
.
Since we are looking for an update Δ X such that
X +Δ X is still invertible, we get the condition

    (X + Δ X)(X̃ + ΔX̃)
         = X X̃ + XΔX̃ + Δ X X̃ + Δ X ΔX̃!= I

which we linearize (assuming “small” Δ X and ΔX̃)
and reformulate (assuming only Δ X X̃ = X̃Δ X) as

    X ΔX̃ + X̃Δ X = I - XX̃ = 0.

Thus —compared to (<ref>)—
we can find also “non-commutative” primal-dual search directions by solving
the following linear system of equations
(with 4n^2+m unknowns):

    [       .       A   I_n^2       .       .;       A       .       .       .       .;       z       .       x       .       .; I_n ⊗X̃       .       . I_n ⊗ X       .;       .       . I_n ⊗Z̃       . I_n ⊗ Z;         ][ Δ x; Δ y; Δ z; Δx̃; Δz̃ ]
    =
    [ c - z - A y;     b - A x;      μ - XZ;   0^n^2 × 1;   0^n^2 × 1 ].

Recall that x= X, etc. Its system matrix has 3n^2+m+1 rows,
that is, the system is underdetermined
(we always assume A = m).
The third block row reads
ZΔ X + XΔ Z = μ - XZ.
If one is familiar with the “classical” issues of
primal-dual interior-point methods like symmetrization,
step-length computation, update of μ, etc. then an
implementation is —at least for moderate sized problems—
straight forward.


An ad hoc implementation applied to the problem
from Section <ref> showed that the primal path
“converges” much faster than the dual path.
We have not investigated that in detail (for different
problem classes) since working with the primal problem
only has the advantage of a much smaller linear system of equations
(to compute a search direction).
It would be interesting to see what happens if one uses alternating
“primal” and “dual” steps,
or how they can be used to accelerate “classical”
interior-point methods.



It is a funny coincidence that the
approach we tried 
to ensure invertibility of an (unknown) matrix X,
namely introducing another matrix X̃ and add
the equations X X̃ = I,
as an alternative to the classical X ≠ 0 
seems to have a drawback when a linearization (like we can use here)
is not appropriate (and Groebner bases are needed) because
the number of unknowns doubles
<cit.>
.
The important observation is the origin of
XZ=μ I from the classical barrier -log X
(respectively log Z)
<cit.>
, 
<cit.>
. So it is rather natural to carry that over directly to
our discrete (iterative) setting as “algebraic”
barrier X X̃ = I (respectively Z Z̃ = I).
Recall that X > 0 does not imply
positive eigenvalues λ_i > 0.






§ FEASIBLE-INTERIOR-POINT METHODS




Now we start (somewhat naively) from scratch
with the basic linear algebra from Section <ref>
in mind and recall the relevant setup for the
semidefinite program (<ref>),

min_0 X ∈_n { CX :
  A_1X = b_1, A_2X = b_2, …, A_mX = b_m }.

Given m+1 square matrices C,A_i ∈_n() = ^n × n
and the vector b ∈^m × n we want to find the
minimum of CX = (C X) such that
X is symmetric and positive semidefinite
and the m constraints A_iX = b_i are fulfilled.
We write x =  X and
A = [ A_1,  A_2, …,  A_m]
and use Example <ref> (from the following section)
for illustration:

A =
[ . 1 . 1 . . . . .; . . 1 . 1 . 1 . .; . . . . . 1 . 1 .; . . . . . . . . 1 ],
  b = [    0; 13/4; 15/4;    1 ]
  and  C =
[ 1 . .; . 0 .; . . 0 ].


The vector space of n × n matrices _n = _n()
decomposes (almost) naturally in
degrees of freedom _n^ξ,
(pure) minimization directions _n^η,
constraint directions _n^ζ and
non-symmetric directions _n^ν,
namely

_n = _n^ξ⊕_n^η⊕_n^ζ_=_n ⊕_n^ν
with m_ν = n(n-1)/2 and
m_ξ + m_η + m + m_ν = n^2 and
bases M^ξ, M^η, M^ζ and M^ν respectively.


However it is not so clear, how to choose M^η.
In some sense m_η =  M^η should be maximal
for (M^η) in row echelon form.
But on the other hand, if there are minimization
directions which “correlate” with other directions
it might be helpful to have more degrees of freedom
for numerical stability (at least in the beginning).
Notice that in the family of examples from
<cit.>
 m_ξ = 0 (see also Table <ref>).



Figure <ref> shows the cone of
symmetric positive definite matrices X = X 0
intersected with the “hyperplane” Ax = b for x =  X.
For a feasible
matrix X_k = X_k 0,
a (symmetric) direction Δ X_k ∈_n^ξ⊕_n^η
and a scalar α_k ∈,
the matrix X_k+1 := X_k + α_k Δ X_k is again feasible,
that is

A(x_k + α_k Δx_k) = b.

If Δ X_k is non-symmetric, that is,
Δ X_k ∈_n^ξ⊕_n^η⊕_n^ν,
a symmetric direction is given by
Δ X_k := 1/2(Δ X_k + (Δ X_k)).
Before we discuss how to find such search directions Δ X_k
(given a feasible positive definite X_k),
we have a look on how to compute the maximal steplength α_k^max
such that X_k + α_k^maxΔ X_k  0.
(Without loss of generality we assume that Δ X_k is such that
α_k^max < ∞.)

To simplify notation we drop the iteration index in the following
and write X = X_k.
Since X is positive definite it admits a Cholesky factorization
X = L L. The maximal steplength
is (<ref>),
α_max = ( (-LΔ X L) )
<cit.>
. So for 0 < τ < 1 we have

X + τα_max ΔX 0.

Now we take a (not necessarily orthonormal) basis M^ξ
(respectively M^η and M^ν)
for _n^ξ (respectively _n^η and _n^ν)
and write M^ξ,η (respectively M^ξ,η,ν)
for [M^ξ,M^η] (respectively [M^ξ,M^η,M^ν])
or simply just M (for M^ξ,η and M^ξ,η,ν).


If the basis M^ξ is not orthogonal one
might be a little bit more careful with respect to
geometric centering, see Definition <ref>
and the following remark.



Steplength scaling can be quite subtle in practice.
For LP's, tpically τ=0.99995 is used
which is too aggressive for “predictor-corrector”
methods 
<cit.>
 therefore τ=0.98 (and adaptive steplength scaling)
is used.
Since here we restrict ourself to feasible paths,
steplength scaling seems to be —except close to
the “boundary” of singular matrices— not a difficult
issue.



From Section <ref> we recall the linearized
condition
(<ref>) for the invertibility
of X + Δ X, namely

X ΔX̃ + X̃ ΔX = 0

for given X and X̃ = X.
Since we start with a feasible X and restrict the search
directions to Δ X ∈ M, this invertibility condition
now reads

    X ΔX̃ + X̃M x̂_=Δ X = 0

for x̂∈^m_ξ+m_η or x̂∈^m_ξ+m_η+m_ν.
Additionally, since we want to minimize,
we need a condition of the form
CΔ X = CMx̂ < 0.
Thus, a minimization search direction can be computed as
a least squares solution of the (underdetermined)
linear system of equations

    [ (I_n ⊗X̃) M     I_n ⊗ X;      ( C) M           . ][  x̂; Δx̃ ]
    =
    [  0; -γ ]

for γ > 0. Now we can use (<ref>)
to compute the steplength α_max and iterate,
at least in principal.
Notice that the search direction computed in (<ref>)
is independent of γ > 0 (modulo scaling).




Although the path (X_i) we would get —at least if
we ignore the numerical issues—
is in the “hyperplane” Ax=b and “positive definite”,
we need to ensure that we indeed approach the minimum.
That this is not clear can be seen in Figure <ref>
(cyan line): There is a first huge step, but then it seems
that the path “converges” to some X' such that CX'≈ 1.005 > 1,
that is, the minimum is not reached.
One has to be very careful here because it might be that
—using multiprecision arithmetics— just a very high number
of iterations is necessary to get arbitrary close to the minimum.
At the end this does not really matter since we want to
reduce the number of iterations as much as possible.



For a case where centering is not needed (or even
must not be used) —and how to adapt the implementation
from Section <ref>—
see (the end of) Section <ref>.
We have not yet investigated that in detail.


There are roughly three different variants of
(primal) feasible-interior-point methods, namely:



 §.§ All Directions (STD)


Use all “feasible” search directions (including the non-symmetric),
symmetrize and center mainly for “staying away from
the boundary” by using a centering parameter 0.65 ≤μ≤ 0.85
and X_new := μ X + (1-μ) X_cen.
This corresponds to  in line 7
in  (Section <ref>)
for algebraic centering.

Some paths for different feasible initial matrices
are shown in Figure <ref>,
some for different centering parameters in
Figure <ref>.
For a typical problem run see (the end of)
Section <ref>.

Instead of the “approximated” algebraic centering,
“exact” geometric centering could be used.
It seems that the optimal centering parameter (at least for
algebraic centering) is around μ≈ 0.75.
But so far it is neither clear if it is (within some range)
problem dependent nor how it can be interpreted.
A dynamical tuning (depending on the relative error)
might be possible.




 §.§ Symmetric Directions (SYM)


Use only symmetric search directions M^ξ,η
and center mainly for “staying close to the center”
by using a centering parameter 0 ≤μ≤ 0.25.
This corresponds to  in line 7
in  (Section <ref>)
for algebraic centering.

Some paths for different feasible initial matrices are
shown in Figure <ref>.
For a problem run (in the context of the relaxation
of a combinatorial problem) without centering
see (the end of) Section <ref>.
For a comparison of the number of iterations
for the example from 
<cit.>
 see Table <ref>
(page tab:nc.jck).

Instead of the “approximated” algebraic centering,
“exact” geometric centering could be used.
The latter could be used to estimate a “good”
centering parameter μ for the former.
On the other hand, the former might be necessary
to get away from the boundary for an “aggressive”
steplength scaling parameter τ < 1 for the latter
(see Remark <ref>).





 §.§ Constructed Directions (ACE or GCE)


If degrees of freedom are available, that is,
m_ξ≥ 1, search directions can be constructed
by using two
different centered matrices X_1 and X_2
with CX_2 < CX_1
and setting Δ X := X_2 - X_1.

For some paths using “approximate” algebraic
centering see Figure <ref>, for
some using “exact” geometric centering
Figure <ref>. An investigation for
m_ξ≥ 2 is open.
For a typical problem run (and some further comments)
see Section <ref>.
Notice however, that the usability in the case of
several minimization directions, that is, m_η≥ 2,
is not yet clear.



















§ FINDING A FEASIBLE INITIAL MATRIX




Even in “classical” interior-point methods (for semidefinite programming)
the role of “good” initial iterates should not be underestimated
<cit.>
. For SDPA a “sufficiently large” multiple of the identity matrix
should be used <cit.>
, the default parameter is λ^*=100.
Since “λ” is heavily used for eigenvalues,
we use “γ” here, for example X_0 = γ I.
Unfortunately we have not found that much beyond short remarks in the literature.
For linear programming one could start with 
<cit.>
.
Since in our case finding a feasible initial iterate is crucial,
we tried several approaches and leave some comments for further investigation.

Actually it turned out that this task seems to be highly non-trivial.
Several of our ad hoc approaches —using variants of the techniques from
the minimization in the previous section— failed in some cases.
So far only one survived and is presented in detail in Section <ref>.
The variant for linear programs in Section <ref> is just
a simplified version. If this also fails in general, a clever
use of semidefinite programming in the context of “norm-minimization”
might work,
maybe in combination with Geršgorin circles
<cit.>
.

Still, starting with a least squares solution to A x =b and using
repeatedly a spectral shift γ I and a correction
A Δ x = b - A(x+(γ I))
seems to work for practical problems.
For small parameters  and δ in the example
<cit.>
 one can increase the shift γ dynamically.

From a different point of view one can ask how the
“projection” of the central ray γ I onto the
constraints A x = b looks like. Not directly but in
terms of the maximal steplength α(γ)
as a function of γ, see Figure <ref>.
If we find any γ such that α(γ) > 1 we
are done. Assuming the existence, starting with a small
γ, it is not difficult to find it with only a few
iterations (increasing γ if the gradient is positive
and decreasing γ if it is negative respectively
continuing with the old iterate and increase γ
less “aggressive”).

The search for an initial iterate is done in two steps:
Firstly, given X and setting X̃ = X,
we get a direction “towards” Ax=b as a least
squares solution to the linear system of equations

    [       A       .; I_n ⊗X̃ I_n ⊗ X ][ Δ x; Δx̃ ]
    =
    [ b-A x;     0 ].

Notice that a scaling of the first block row might
be necessary depending on the norm of X.
Since Δ X is not necessarily symmetric,
we set Δ X := 1/2(Δ X + (Δ X)).
(In the implementation it is always done for numerical reasons.)
Now we can calculate α = α_max(X,Δ X)
and update X := X + min{ 1, τα}Δ X
for some τ < 1.
Secondly, for numerical stability we use
either algebraic centering (<ref>)
if m_ξ≥ 1 or a modified version by computing
a least squares solution Δ X = (M x̂) for
M = [ M_η, M_ν ]
to the linear system of equations

    [ I_n ⊗X̃ M   I_n ⊗ X ][  x̂; Δx̃ ]
    =
    [ (1-μ) x̃ ]

for 0 ≤μ≤ 1.
In both cases we update
X := X + 1/2( Δ X + (Δ X)).
Notice that X̃ = X needs to be updated in between.

Using M= M_η only in the last step did not work
for the example from 
<cit.>
. Some numbers for the initial iterates (for the presented approach)
are in Table <ref>.
For the problem from the next section see also
Table <ref>. Notice that in this case
there is a degree of freedom for “approximated” algebraic
centering.






§ SUM OF SQUARES


Suppose, that we want to find the global minimum of the
polynomial p(x) = 2 + 13/4 x^2 + 15/4 x^3 + x^4
(its graph is the thick line in Figure <ref>).
What is easy (in this case) using high-school
mathematics will turn out to be rather involved using sum-of-squares
and semidefinite programming 
<cit.>
, although the main idea is so simple: Minimize η such that
p_η(x) = η + 13/4 x^2 + 15/4 x^3 + x^4
can be written as

p_η(x) = (q_η,1(x))^2 + (q_η,2(x))^2
    + …+ (q_η,n_η(x))^2
    = ( q̅(x))q̅(x)

for some n_η∈ and polynomials q_η,i∈[x].




Not every non-negative polynomial can be written as a sum of squares
<cit.>
. However, an “approximation” might be possible
<cit.>
. 



Let p_η(x) = η + 13/4 x^2 + 15/4 x^3 + x^4 ∈[x].
It can be written as

p_η(x) =
[   1   x x^2 ]
[        η        0     -ξ/2;        0 13/4 + ξ     15/8;     -ξ/2     15/8        1 ]_=:Q(ξ,η)
[   1;   x; x^2 ].

For η=2, and (approximately) 0.30775 < ξ < 1.39437,
the matrix Q(ξ,2) is positive definite.
For η=1, the matrix Q(0.5,1) is positive semidefinite,
and for η<1 it has at least one negative eigenvalue.
Thus the global minimum of p_2(x) is 1 (at x=-2)
which we can find by solving the semidefinite program

min_0 X ∈_3 
  { CX :
    A_1X = b_1, A_2X = b_2
    A_3X = b_3, A_4X = b_4 }

with

    A_1    = 
    [ . 1 .; 1 . .; . . . ],  
    A_2 =
    [ . . 1; . 1 .; 1 . . ],  
    A_3 =
    [ . . .; . . 1; . 1 . ],  
    A_4 = 
    [ . . .; . . .; . . 1 ],
    
    b    =
    [    0; 13/4; 15/4;    1 ]  and  
    C = 
    [ 1 . .; . . .; . . . ].

We write x =  X and denote by
A∈^m× n^2 (m=4, n=3)
the matrix with rows
( A_i) to be able to write
A (X) = b for the constraints.


The following problem runs indicates that (in principal)
quadratic convergence is possible,
independent of the starting point.
The resulting path for the initial matrix with ξ=0.75 and η=2
is shown in Figure <ref> (blue line). 
Table <ref> shows the number of iterations compared
with some classical solvers and that finding an initial feasible
matrix is an important issue.
For a basic discussion see Section <ref>.
It might be possible to develop problem specific methods
for initial iterates.
Notice in particular the right column in the outputs with
the (Frobenius) norm of X
and the advantage of centered search directions
for the (approximate) solution.













§ IMPLEMENTATION NCMIN


The following implementation in Octave <cit.>
 is meant for educational purposes only without any warranties.
By changing the function  to 
it can be used also in Matlab <cit.>
. (Notice however that the output can be slightly different due
to rounding, etc.)
It consists of a simple standalone linear solver 
with net less than 70 lines (see Section <ref>
for the code and Section <ref> for a typical application)
and a semidefinite solver  (Section <ref>)
building on  (Section <ref>) and
 (Section <ref>)
with net less than 210 lines together.
For a typical problem run of the latter see Section <ref>.


Please keep in mind that the code is not streamlined to be
able to follow each step without jumping into a cascade of
functions, in particular to serve as a thread through the
theory. Resisting the temptation to further “optimize”
it was not that easy. But already for sparse matrices
some algorithms for “simple” tasks like solving a
underdetermined linear system
of equations change significantly and hand­ling the differences
between Octave and Matlab would make the
code unnecessary complicated (to read). Those who are interested in
solving large scale problems might find 
<cit.>
 helpful. For a general background in numerical linear
algebra we refer to 
<cit.>
.
One typical “inefficiency” is that the inverse of X is computed
twice, for example in lines 26 and 34 (Section <ref>)
or lines 80 and 88 (Section <ref>). But to really
attack that seriously one needs to think about approximating the
inverse (and the consequences for the solution) in the case of
sparse matrices (“fill in”, etc.)



If something goes wrong or the result is not like
expected then one should check every single input argument
and try different parameters. Is the matrix C symmetric?
Is the (primal) problem feasible? Is it bounded?
What is the rank (what are the singular values) of A?
Is the tolerance appropriate (especially for a badly conditioned
problem)? Is the steplength scaling τ too “aggressive”?



If everything looks nice, have a look on the
dual problem. Is the duality gap zero?
What are the ranks (with respect to some tolerance) of
the (almost) optimal pair X and Z?
Is the complementary primal-dual (almost) optimal
solution strict complementary
<cit.>
?
Is the result (within some tolerance) independent
of the initial starting value?
Yield the other search strategies the same result?
Is the result “stable” with respect to a change of
the parameters (although the number of iterations changes
significantly)? Is the maximal number of iterations
sufficiently large? Are there some “problem specific”
(plausibility) checks?



language=Octave,basicstyle=
numberblanklines=true, basewidth=0.5em,0.45em
morecomment=[l][]%
numbers=left, numberstyle=, stepnumber=5, firstnumber=1
frame=lines,aboveskip=0.5ex,belowskip=1.5ex



 §.§ NCMIN Examples


These tiny examples are for basic testing and detailled
investigations. The matrix  in line 17 yields
an unbounded problem.
The family of examples from 
<cit.>
 are useful to check the “limits” of a solver.
Although there are some similarities to the
example from Section <ref>, there
are no degrees of freedom and hence
search directions cannot be constructed
(search types 3 and 4 in ).
SeDuMi
<cit.>
 provides an interface to the sparse format of SDPA
<cit.>
.


[firstline=7,label=lst:nc.min.ex,
  title=]/ncminex.m




 §.§ NCMIN Linear Solver


For the theoretical setting see Section <ref>
with the main linear system (<ref>)
corresponding to “” in lines 60–63.
Search type 2 (line 7) is only for illustrating that a naive
approach —not ensuring invertibility of the (diagonal) matrix
as a linearized condition—
can easily trigger wrong results.
There are several parts which need further investigation.

Firstly, finding a feasible initial vector with positive
entries. This is a simplified version of that presented
in Section <ref> and computes iteratively
a direction “towards” Ax=b subject to the linearized
invertibility condition

x̃ Δx + x Δx̃ = 0

until the maximal steplength α = α_max(x,Δ x) > 1
or the residuum b-Ax is smaller than a predefined tolerance
(or the maximal number of iterations is reached).

And secondly, it would be interesting how a modified version
of search type 1 including a correction of the linearization
error (compare with classical predictor-corrector
methods <cit.>
) works. 

For a small linear (scheduling) program with m=24 and n=80
 needs
7 iterations for the initial iterate and
18 iterations for the minimization (with relative error
10^-8).
SeDuMi 1.30 
<cit.>
 needs 8 iterations for the same problem with standard settings
(predictor-corrector method)
and 32 iterations with the “not recommended”
algorithm 0.
Notice however that the latter is a highly developed program,
while  has a lot of potential for improvements:
Reasonable (relative) accuracy is reached quite fast (after a few
iterations) and it is possible to stop at any point
without impact on the feasibility of the (approximate) solution.
Additionally, the linear systems of equations to solve
are smaller than for classical primal-dual interior-point methods …




[firstline=7,label=lst:nc.min.lp,
  title=]/ncminlp.m




 §.§ NCMIN Vector Space Decomposition


The basic theory is explained in Section <ref>.
Here,  corresponds to M_ξ,
 to M_η
and  to M_ν.
The loop from line 21 to 27 takes care of distinguishing
directions between centering and minimization
and does some simple
“normalization” of the vectors in M_ξ.
Notice the limitations of this (naive) approach
with respect to the use of dense matrices
(in total) of size n^2 × n^2.
Although this task is not directly related to “preprocessing”
it is maybe worth to view it from a more general context
<cit.>
.


[firstline=7,label=lst:nc.min.dof,
  title=]/ncmindof.m




 §.§ NCMIN Initial Feasible Matrix


For details see Section <ref>.
For Example <ref> one can subtract  from
a matrix X ∈_3 to read off ξ in entry (2,2)
because M_ξ=[0,0,-1/2,0,1,0,-1/2,0,0]
and η in entry (1,1) because
M_η=[1,0,…,0].
In the following example we have ξ≈ 0.767 and η≈ 1.987:


[firstline=7,label=lst:nc.min.ini,
  title=]/ncminini.m



 §.§ NCMIN Semidefinite Solver


For details see Section <ref>.
For a “typical” path for search type 1 (all directions)
see Figure <ref>,
type 2 (symmetric directions) Figure <ref>,
type 3 (algebraic centered) Figure <ref> and
type 4 (geometric centered) Figure <ref>.
In any case it is a greedy method, that is,
after computing a “feasible” search direction,
it takes (almost) the maximal steplength.
The main ingredient for making that work is
centering, or in other words, “staying away from the boundary (of
singular matrices)”. See Section <ref>,
in particular Example <ref> and <ref>.

[firstline=7,label=lst:nc.min.sdp,
  title=]/ncminsdp.m



§ (NOT) ONLY FOR STUDENTS


Starting with semidefinite programming is difficult,
even with a “linear” background. Maybe because there is
no end in sight and sooner than later one is in the middle
of current research or needs deep results which are hard
to understand even roughly within a few weeks (or months).
Just to get an impression:

“Semidefinite programming (SDP) is probably the most important
new development in optimization in the last two decades. (…)
Notice that semidefinite programming is a far reaching extension
of linear programming (LP) …” 
<cit.>
.

Its importance cannot be overestimated due to the many
practical applications, in particular in engineering
<cit.>
 and the availability of highly developed
solvers. But semidefinite programming is also fascinating
from a purely mathematical point of view because of all
the manifold connections to (at a first glance) different
mathematical areas and the rich “explorable” structure.
The initial hurdle of many definitions and different notations
is admittedly high. Getting confident in (applied)
non-commutative algebra takes some time …

So, where should one start? With small
problems that can be “touched” and from a perspective
one is familiar with, no matter if it is of geometrical,
analytical, algebraic or applied (in the sense of programming)
nature. Almost all here originated from a “tiny” semidefinite
problem with 3 × 3 matrices. And if there were not
Halmos' recommendation to stop
<cit.>
, this paper would have never been finished (and therefore not
exist).

One could start with the simple linear matrix pencil
A = A(x,y) ∈^2 × 2
from <cit.>
 and try to figure out, which “object” it describes
with respect to the real variables/parameters
x and y such that A is positive (semi-)definite:

A = A(x,y) = 
[ 1+x   y;   y 1-x ]
=

[ 1 .; . 1 ]_=:A_0
⊗1 +

[  1  .;  . -1 ]_=:A_x
⊗x +

[ . 1; 1 . ]_=:A_y
⊗y.

Here the zeros are replaced by (lower) dots to emphasize the
structure and the tensor product “⊗” just
means that we “plug in” the variables (into the respective
matrix). Depending on the context/area, matrix pencils are
written in many different ways, for example
A = A(x,y) = A_0 + x A_x + y A_y.
The “constant” coefficient matrix A_0 usually plays
a special role (and is often just the identity matrix).

How can we find a matrix X_0
“inside” this convex object?
And given X_0, how can we find the “center”
algorithmically by using (<ref>)
respectively lines 61ff or 86ff in 
(Section <ref>)?
Which other objects can be represented in such a way
(by 2 × 2 matrices)?
How many variables are needed at most for 3 × 3 matrices?
(Recall that we are talking about symmetric matrices only.)

Questions about the representation of convex sets by
linear matrix pencils (or “linear matrix inequalities”, LMI's)
are very difficult in general. But it is still worth to get
an idea by having a look on the many illustrations in
<cit.>
.

Another possibility is to use the power of semidefinite
programming for concrete problems. One application
is discussed in Section <ref>.
Another is the relaxation of
hard combinatorial problems
<cit.>
. For an illustration we restrict ourselves to 0-1 variables
and a quadratic functional f : ^2 →.
Let C ∈^3 × 3,

A_0 =
[ 1 . .; . . .; . . . ],
  A_1 =
[    . -1/2    .; -1/2    1    .;    .    .    . ]
  and  A_2 = [    .    . -1/2;    .    .    .; -1/2    .    1 ]

and define x̅ = [ 1, x_1, x_2 ]
and X̅ = x̅x̅.
Now the optimum from the (non-convex) combinatorial problem
can be bounded from below by a semidefinite program:

    min_x ∈^2{ f(x) : x_i ∈{ 0,1 }}   
    ≥min_0 X ∈_3{CX.
         : . A_0X = 1, A_1X = 0, A_2X = 0 }.

Without the application of the trace the constraints
reads (typically)

A_1 X = 
[    . -1/2    .; -1/2    1    .;    .    .    . ]
[       1     x_1     x_2;     x_1   x_1^2 x_1 x_2;     x_2 x_1 x_2   x_2^2 ]
=
[         -x_1/2              *              *;                -x_1/2 + x_1^2              *;                             *              0 ].

Some trivial but crucial questions:
What is the rank of X̅? What that
of X  0 in the SDP (<ref>)?
And which linear algebraic concept(s) are useful
to measure the (approximate) “numerical” rank
of some X  0?

The following works only for symmetric
search directions ( in line 7)
and deactivated centering (by uncommenting
line 106 in Section <ref> and
line 45 in Section <ref>).
There might be a simple explanation …










§ EPILOGUE


Or prologue? Is there enough here that makes it worth
to look back? At the end everything grew out naturally 
from the trivial observation that the matrix multiplication is non-commutative
(in general). But what is surprisingly simple from a latter perspective is
not at all obvious in the beginning. And what we know from an old
proverb, namely that the search for a needle in a haystack
does not get easier when one increases the stack, seems to
hold true in particular for a “matrix valued” needle
(pointing towards the minimum). So one challenging question
remains: How can we find “good”
search directions —as solutions to underdetermined
linear systems of equations— along a
feasible interior path?

Here we indicated only that it might be possible
by combining some rather simple techniques. But even if it turns
out that these are not usable in general it could stimulate
new ideas and might explain some “strange” behaviour
<cit.>
 or “hard” problem 
<cit.>
.
In any case one should keep Todd's words
from the abstract of <cit.>
 in mind:
“The most effective computational methods are not always
provable efficient in theory, and vice versa.”
Very interesting to read is also 
<cit.>
 from Forsgren, Gill and Wright about the roots of
non-linear programming:
“Furthermore, a simplex-centric world view had the effect that
even `new' techniques mimicked the motivation of the simplex method
by always staying on a subset of exactly satisfied constraints.”
With other words: When is it necessary to leave a
feasible path?

“Pedagogical and philosophical issues remain about the
best way to motivate interior-point methods —perturbing optimality
conditions? minimizing a barrier function?— and the multiplicity
of viewpoints continues to create new insights and new algorithms.”
<cit.>

There is quite a lot we could not even touch here.
So we just mention a non-complete list
of papers (of quite different flavour)
which might help to find further literature.
For an overview and/or introduction
one could start with
<cit.>
.
Although one should use the latest versions (and the
current documentation) of the implementations, the
underlying publications
<cit.>
 are in some sense timeless.
More on interior-point methods, search directions,
predictor-corrector methods, etc. can be found in
<cit.>
.
Around verification, exact solutions, high precision, etc. there are
<cit.>
 and
<cit.>
.
Around strong duality, regularization, preprocessing,
hard and bad semidefinite programs, etc. one should have a look in 
<cit.>
.

On one hand, semidefinite programming is a special
case of conic programming
<cit.>
.
On the other, it
is very useful for relaxations of hard combinatorial
problems. So there is a lot which goes far beyond
the basics we presented here:
<cit.>
, <cit.>
.
Since semidefinite programming is tight together
with numerics (least squares approximation, iterative
linear solvers, finite precision arithmetics, etc.) one
should mention at least
<cit.>
. But semidefinite programming is also interwoven
with semialgebraic geometry 
<cit.>
.

At a first glance a topic like
“non-commutative sum of squares”
<cit.>
 seems to be far away. But indeed it was
one of the initial questions around
non-commutative convexity 
<cit.>
 which spans the bridge to my
“non-commutative” algebraic research.

There is a lot which is not discussed here (in detail).
One topic is numerics (in general) and sparse matrices
(in particular) in combination with iterative linear
solvers. One starting point for literature is
<cit.>
. Another is that around estimating good (maybe problem
specific) parameters. Since there are not that many in
the presented methods it should not be that difficult
to change them dynamically. One advantage is that there
is no need for a subtle control for μ→ 0^+.
One disadvantage is that there is no continuous path
(and no analysis) anymore, so it might be difficult to
estimate a rate of convergence.

It is clear that the different methods can be combined,
in particular with “classical” methods. Additionally,
information from the dual problem can be used for
verification
<cit.>
. That this should not be underestimated is shown
in particular in 
<cit.>
.

The list of (not just small) improvements is long:
enable the use of a relative error and/or numerical limits 
as stopping criteria, export the initial data or some
intermediate results in SDPA's 
<cit.>
 sparse format, etc.
For the issues with respect to benchmarking
one could start with <cit.>
, for large scale programs with
<cit.>
.

In addition, thorough testing is necessary,
in particular with the problems from SDPLIB
<cit.>
 and the library <cit.>
. However, even if some basic testing could be done
automatically, a deeper analysis of the reasons
where something goes wrong (with the presented approach
or some classical solver) will definitely help to
get further insight in semidefinite programming.
In this context the survey
<cit.>
 could be interesting.




§ ACKNOWLEDGEMENT


I am very grateful for Roland Speicher's invitation to
Saarbrücken and the opportunity to give a talk on this
topic and continue in particular “non-commutative”
algebraic discussions in October 2018.
I thank Gabor Pataki and Michael J. Todd for hints on the
literature and Birgit Janko for very valuable feedback on
several drafts of this work.

This work has been supported by research subsidies granted by
the government of Upper Austria
(research project “Methodenentwicklung für Energieflussoptimierung”).

tocsectionBibliography
alpha


