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

TL;DR

This paper introduces a new family of primal feasible-interior-point methods for semidefinite programming that leverage non-commutative search directions, addressing limitations of traditional commutative approaches.

Contribution

It proposes novel non-commutative search directions for interior-point methods, expanding the optimization toolkit for semidefinite programming beyond traditional commutative paths.

Findings

New non-commutative paths demonstrated

Implementation provided for further research

Potential to solve previously hard semidefinite problems

Abstract

Since more than three decades, interior-point methods proved very useful for optimization, from linear over semidefinite to conic (and partly beyond non-convex) 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.

Tables

Table 1. Table 1: The number of iterations (that for finding a feasible starting vector + + that for minimization) for the linear program from Example  2.7 . “PC” stands for predictor-corrector which implies more expensive steps. For SDPA the sparse format with 5 diagonal blocks of size 1 × 1 1 1 1\times 1 is used, for SeDuMi one can specify that the problem is linear. For both solvers their respective standard settings are used. One must be careful with a direct comparison due to the large number of possible tuning parameters. There are almost no parameters in ncminlp .

Solver	Algorithm	Iterations
SeDuMi 1.30 [Stu99]	0	17
	1 v-corr.	5
	2 $xz$-corr (PC)	3
SDPA 7.3.9 [YFK03]	PC, $\tau =0.9$	12
	PC, $\tau =0.99999$	9
ncminlp (Section 7.2)	$\tau =0.99999$	$1+5$
	$\tau =1-{10}^{-10}$	$1+3$

Table 2. Table 2: The number of iterations (that for finding a feasible starting matrix + + that for minimization) for the semidefinite program from [ JCK08 , Section 6] . For SDPA the standard settings did not work for small parameters and therefore λ ∗ = 10 6 superscript 𝜆 superscript 10 6 \lambda^{*}=10^{6} was used in all cases. Iterations in parentheses mean that the result was not very accurate.

		SeDuMi 1.30	SDPA 7.3.9	ncminsdp	Section 7.5
$ϵ$	$\delta$	[Stu99]	[FFK+08]	standard	symmetric
${10}^{-2}$	${10}^{-2}$	9	19	$4+4$	$4+5$
${10}^{-4}$	${10}^{-4}$	16	18	$6+4$	$6+7$
${10}^{-6}$	${10}^{-6}$	23	16	$6+4$	$6+7$
${10}^{-8}$	${10}^{-8}$	(28)	15	$9+\text{failed}$	$9+\left(6\right)$

Table 3. Table 3: The number of iterations (that for finding a feasible starting matrix + + that for minimization) for the semidefinite program from Example  6.1 . “PC” stands for predictor-corrector which implies more expensive steps.

Solver	Algorithm	Iterations
SeDuMi 1.30 [Stu99]	0	19
	1 v-corr.	12
	2 $xz$-corr (PC)	9
SDPA 7.3.9 [FFK+08]	PC	13
ncminsdp (Section 7.5)	standard	$3+5$
	symmetric	$3+5$
	algebraic centered search	$3+5$
	geometric centered search	$3+4$

Equations

$$XZ=\mu I.$$

$$⟨X,Z⟩:=\mathrm{tr}(X^{\top }Z)=\sum _{j=1}^n\sum _{i=1}^n{x}_{ij}{z}_{ij}\text{“=”}(\mathrm{vec}{X}^{\top }{)}^{\top }\mathrm{vec}Z$$

$${\mathcal{M}}_n=\underset{={\mathcal{S}}_n}{\underbrace{{\mathcal{M}}_n^{\xi }\oplus {\mathcal{M}}_n^{\eta }\oplus {\mathcal{M}}_n^{\zeta }}}\oplus {\mathcal{M}}_n^{\nu }$$

$$[0,\infty]\ni\alpha_{\max}^{\pm}=\alpha_{\max}^{\pm}(X,\Delta X)=\sup\bigl{\{}\alpha\in\mathbb{R}\mid X\pm\alpha\Delta X\succ 0\bigr{\}}$$

$$\alpha_{\max}=\bigl{(}\lambda_{\text{max}}(-L^{-1}\Delta XL^{-\!\top})\bigl{)}^{-1}$$

$$X_{\text{cen}}=X_{\text{cen}}(M_{\xi})=\frac{1}{2}\sum_{X_{\xi}\in M_{\xi}}\bigr{(}\alpha_{\max}^{+}X_{\xi}-\alpha_{\max}^{-}X_{\xi}\bigr{)}$$

$$\bigl{\lVert}(X^{\prime})^{-1}\bigr{\rVert}=\min_{X_{\xi}\in\mathcal{M}_{n}^{\xi}}\bigl{\lVert}(X+X_{\xi})^{-1}\bigr{\rVert}$$

$$X=\left[\begin{array}{cc}1& 0.9\\ 0.9& 1\end{array}\right],\Delta X=\left[\begin{array}{cc}-1& .\\ .& -1\end{array}\right],X_{\xi }=\left[\begin{array}{cc}.& 1\\ 1& .\end{array}\right].$$

$$\bigl{\lVert}\bigl{(}X+(1-10^{-k})\alpha_{\max}\Delta X\bigr{)}^{-1}\bigr{\rVert}\approx 10^{k+1}$$

$$X\otimes Z=\left[\begin{array}{cccc}{x}_{1,1}Z& x_{1,2}Z& \dots & x_{1,n}Z\\ x_{2,1}Z& x_{2,2}Z& \dots & x_{2,n}Z\\ ⋮& ⋮& \ddots & ⋮\\ x_{n,1}Z& x_{n,2}Z& \dots & x_{n,n}Z\end{array}\right].$$

$$\left[\begin{array}{cc}(I\otimes \tilde{Z})M_{\xi }& I\otimes Z\end{array}\right]\left[\begin{array}{c}\hat{z}\\ \Delta \tilde{z}\end{array}\right]=(1-\mu )\tilde{z}$$

$$\mathrm{mat}{M}_{\xi }\hat{z}\approx \left[\begin{array}{cc}0& -1\\ -1& 0\end{array}\right]\text{and thus}{X}^{\prime }=\left[\begin{array}{cc}0.9& -0.1\\ -0.1& 0.9\end{array}\right]$$

$$\Vert Ax-b{\Vert }_2=\underset{\bar{x}\in {\mathbb{R}}^n}{\min }\Vert A\bar{x}-b{\Vert }_2.$$

$$\min_{0\leq x\in\mathbb{R}^{n}}\bigl{\{}c^{\!\top}x:Ax=b\bigr{\}}.$$

$$\max_{\begin{subarray}{c}0\leq z\in\mathbb{R}^{n}\\ y\in\mathbb{R}^{m}\end{subarray}}\bigl{\{}b^{\!\top}y:A^{\!\top}y+z=c\bigr{\}}.$$

$$\displaystyle\frac{\mathrm{d}}{\mathrm{d}z_{k}}\biggl{(}b^{\!\top}y+\mu\sum_{j=1}^{n}\log z_{j}\biggr{)}$$

$$\left[\begin{array}{ccc}.& A^{\top }& I\\ A& .& .\\ \mathrm{diag}z& .& \mathrm{diag}x\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]=\left[\begin{array}{c}c-z-A^{\top }y\\ b-Ax\\ (\mu -x_j{z}_j{)}_{j=1}^n\end{array}\right]$$

$$\left[\begin{array}{cc}A& .\\ c^{\top }M& .\\ \mathrm{diag}\tilde{x}& \mathrm{diag}x\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \tilde{x}\end{array}\right]=\left[\begin{array}{c}b-Ax\\ -\gamma \\ 0\end{array}\right]$$

$$\left[\begin{array}{cc}{c}^{\top }M& .\\ (\mathrm{diag}\tilde{x})M& \mathrm{diag}x\end{array}\right]\left[\begin{array}{c}\hat{x}\\ \Delta \tilde{x}\end{array}\right]=\left[\begin{array}{c}-\gamma \\ 0\end{array}\right]$$

$$\left[\begin{array}{cc}A& .\\ \mathrm{diag}\tilde{x}& \mathrm{diag}x\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \tilde{x}\end{array}\right]=\left[\begin{array}{c}b-Ax\\ 0\end{array}\right]$$

$$A=\left[\begin{array}{ccccc}-2& 1& 1& .& .\\ -1& 2& .& 1& .\\ 1& .& .& .& 1\end{array}\right].$$

$$\min_{0\preceq X\in\mathcal{S}_{n}}\Bigl{\{}\langle C,X\rangle:\langle A_{1},X\rangle=b_{1},\langle A_{2},X\rangle=b_{2},\ldots,\langle A_{m},X\rangle=b_{m}\Bigr{\}}.$$

$$\max_{\begin{subarray}{c}0\preceq Z\in\mathcal{S}_{n},\\ \atop y\in\mathbb{R}^{m}\end{subarray}}\Bigl{\{}b^{\!\top}y:y_{1}A_{1}+y_{2}A_{2}+\ldots+y_{m}A_{m}+Z=C\Big{\}}.$$

$$⟨A_i,X⟩$$

$$\sum _{i=1}^m{y}_i{A}_i+Z$$

$$XZ$$

$$⟨C,X⟩-b^{\top }y=⟨X,Z⟩\ge 0.$$

$$\left[\begin{array}{ccc}.& A^{\top }& I_{n^2}\\ A& .& .\\ I_n\otimes Z& .& I_n\otimes X\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]=\left[\begin{array}{c}c-z-A^{\top }y\\ b-Ax\\ \mathrm{vec}(\mu I_n-XZ)\end{array}\right].$$

$$(X+\Delta X)(Z+\Delta Z)$$

$$⟨C,X_{\text{opt}}⟩-b^{\top }{y}_{\text{opt}}=⟨X,Z_{\text{opt}}⟩=0,$$

$$(X+\Delta X)(\tilde{X}+\Delta \tilde{X})$$

Full text

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.

Remark. This is clearly a simplified point of view since semidefinite programming is deeply connected with linear matrix inequalities (LMI’s) [Vin12]. 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 3.9) 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 $\mu>0$ we consider the equation

[TABLE]

Then clearly $Z=\mu X^{-1}$ and since $XX^{-1}=X^{-1}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 [Tod01, Section 6.1], namely the need for symmetrization (of search directions). For a more detailed discussion we recommend [Tod99]. (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.” [HP06, Section 7.1]

Those who are new to semidefinite programming and want to grow their own plant are highly encouraged to start with [HP06, Section 6.1] and take the iteration scheme from [AHO98, Page 6] as a template for an implementation. In parallel —for a detailed introduction, a bigger context and further references— we recommend to have [Nem07], [Tod01], [VB96] and/or [WSV00] at hand. Additionally a classical book on linear programming like [Dan16] can be very useful.

In Section 1 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 2 before we prepare the framework for the investigation of path-following methods for semidefinite programming in Section 3 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 4 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 5. The illustrations are based on a concrete problem discussed in Section 6 and a simple implementation (in Octave) in Section 7 which is intended to serve as a starting point for further investigations. Section 8 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 9.

Remark. 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 2 (page 2) and Figure 3 (page 3).

Remark. 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.

Remark. 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 5 (page 5) a lot of questions arise immediately. Some of them are adressed in Section 9.

1 (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 [Lap18] which are heavily used in Octave [Oct18] and Matlab [Mat18]. For a thorough discussion of applied numerical linear algebra we refer to [Dem97], for linear algebra in semidefinite programming to [WSV00, Chapter 2]. In Section 9 we point out more specialized literature.

Except for the notion of “centering a matrix” —which has a natural geometric interpretation (see for example Figure 3) 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 $\mathbb{R}^{n\times n}\simeq\mathbb{R}^{n^{2}}$ and the algebra of matrices (with the non-commutative multiplication) $\mathcal{M}_{n}(\mathbb{R})$.

The set of the natural numbers is denoted by $\mathbb{N}=\{1,2,\ldots\}$, that of the real numbers by $\mathbb{R}$. 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 $\mathcal{M}_{n}=\mathcal{M}_{n}(\mathbb{R})$ and use their respective lower case letters usually to denote their “vectorized” version in $\mathbb{R}^{n^{2}}$, that is, $x=\operatorname{vec}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=\operatorname{mat}x$. There is one exception to be consistent with the classical notation in linear programming in Section 2, namely $X=\operatorname{\mathrm{diag}}x$ for the diagonal matrix constructed by a vector $x\in\mathbb{R}^{n\times 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 $\operatorname{\mathrm{diag}}x\operatorname{\mathrm{diag}}z=\operatorname{\mathrm{diag}}z\operatorname{\mathrm{diag}}x\in\mathcal{M}_{n}$.

Remark. 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 3) 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 $\tilde{x}\in\mathbb{R}^{n}$ such that $\operatorname{\mathrm{diag}}x\operatorname{\mathrm{diag}}\tilde{x}=I_{n}$.

We are mainly working with symmetric matrices $X=X^{\!\top}\in\mathcal{S}_{n}=\mathcal{S}_{n}(\mathbb{R})\subsetneq\mathcal{M}_{n}\simeq\mathbb{R}^{n\times n}$. A (not necessarily symmetric) matrix $X\in\mathcal{M}_{n}$ is called positive definite (written as $X\succ 0$) if $y^{\!\top}Xy>0$ for all $y\in\mathbb{R}^{n}$ and positive semidefinite (written as $X\succeq 0$) if $y^{\!\top}Xy\geq 0$ for all $y\in\mathbb{R}^{n}$. The trace $\operatorname{\mathrm{tr}}:\mathcal{M}_{n}\to\mathbb{R}$ induces the inner product $\langle.,.\rangle:\mathcal{M}_{n}\times\mathcal{M}_{n}\to\mathbb{R}$,

[TABLE]

for $X=(x_{ij})$ and $Z=(z_{ij})$. By $X\bot Z$ we denote that $X$ and $Z$ are orthogonal, that is, $\langle X,Z\rangle=0$.

A symmetric matrix $X=X^{\!\top}$ has $n$ (not necessarily distinct) real eigenvalues $\lambda_{\min}(X)=\lambda_{1}\leq\lambda_{2}\leq\ldots\leq\lambda_{n}=\lambda_{\max}(X)\in\mathbb{R}$. Our main concern is the decomposition of $\mathcal{M}_{n}$ (as a vector space) in $\ell$ subspaces $\mathcal{M}_{n}^{i}$ of dimension $m_{i}$, that is, $\mathcal{M}_{n}=\mathcal{M}_{n}^{1}\oplus\mathcal{M}_{n}^{2}\oplus\ldots\oplus\mathcal{M}_{n}^{\ell}$ (with $m_{1}+m_{2}+\ldots+m_{\ell}=n^{2}$). As indices we will use mostly $i\in\{\xi,\eta,\zeta,\nu\}$ and define the (vectorsub-)space of non-symmetric matrices as the orthogonal complement of the symmetric matrices, that is, $\mathcal{M}_{n}^{\nu}=(\mathcal{S}_{n})^{\bot}=\operatorname{span}M_{\nu}$ with basis $M_{\nu}\in\mathbb{R}^{n^{2}\times m_{\nu}}$ (of $m_{\nu}=n(n-1)/2$ column vectors). With a decomposition of $\mathcal{S}_{n}$ we have

[TABLE]

with basis $M=[M_{\xi},M_{\eta},M_{\zeta},M_{\nu}]$ and $n^{2}=m_{\xi}+m_{\eta}+m_{\zeta}+m_{\nu}$. Notice that some $m_{i}$ can be zero and recall that $\operatorname{rank}M_{i}=m_{i}$. By abuse of notation we write X\in M_{i}=\bigl{\{}M_{i}^{(1)},\ldots,M_{i}^{(m_{\xi})}\bigr{\}} for $X=\operatorname{mat}M_{i}^{(j)}$ for some $1\leq j\leq m_{\xi}$.

Remark. The decomposition of the vector space of symmetric matrices into “degrees of freedom” $\mathcal{M}_{n}^{\xi}$, “minimization directions” $\mathcal{M}_{n}^{\eta}$ and “constraints” $\mathcal{M}_{n}^{\zeta}$ is used in particular in Section 4, illustrated in Figure 3 (page 3). 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})\in\mathcal{M}_{n}$ we denote by $\lVert X\rVert=\sum x_{ij}^{2}=\langle X,X\rangle$ the Frobenius norm which is just the “classical” vector norm $\lVert X\rVert=\lVert\operatorname{vec}X\rVert_{2}$. Very often we need to “update” $X$, written as $X+\alpha\Delta X$ with “direction” $\Delta X\in\mathcal{M}_{n}$ (“$\Delta X$” is one symbol and we write $\Delta x$ for “$\Delta\operatorname{vec}X$”) and “steplength” $\alpha\in\mathbb{R}$. When $X$ is symmetric and positive definite, that is $X=X^{\!\top}\succ 0$ and $\Delta X\in\mathcal{S}_{n}$ we denote by

[TABLE]

the maximal “steplength(s)”. We write $\alpha_{\max}$ for $\alpha_{\max}^{+}$. Furthermore $X=X^{\!\top}\succ 0$ admits the Cholesky factorization $X=LL^{\!\top}$ (with lower triangular $L$). It can be used to compute the maximal steplength [AHO98, Section 2]:

[TABLE]

which is simply $\alpha_{\max}=1/\max\operatorname{\mathrm{diag}}(-X^{-1}\Delta X)$ for diagonal matrices $X$ and $\Delta X$. For $\tau<1$ we can ensure $X+\tau\alpha_{\max}\Delta X\succ 0$. However, from a numerical point of view, we could run into troubles if $\tau=1-\varepsilon$ for $\varepsilon\ll 1$ since $\lVert X\rVert\ll\lVert X^{-1}\rVert$. Now let $\mathcal{M}_{n}^{\xi}=\operatorname{span}M_{n}^{\xi}$ such that $X_{\xi}\bot\Delta X$ for all $X_{\xi}\in M_{\xi}$ and $\beta^{\pm}=\alpha_{\max}^{\pm}(X,X_{\xi})$ assuming $0<\beta^{\pm}<\infty$. If we have the freedom to replace $X$ by some $X^{\prime}=X+X_{\xi}$ (strictly) “between” $X-\beta^{-}X_{\xi}$ and $X+\beta^{+}X_{\xi}$ such that $\lVert(X^{\prime})^{-1}\rVert\ll\lVert X^{-1}\rVert$, we should do that. We will call that “centering” of $X$ (with respect to $X_{\xi}$). We will refer to minimizing \bigl{\lVert}(X+X_{\xi})^{-1}\bigr{\rVert} for $X_{\xi}\in\mathcal{M}_{n}^{\xi}$ as “algebraic” centering and $X+\textstyle{\frac{1}{2}}(\beta^{+}-\beta^{-})X_{\xi}$ as “geometric” centering.

Remark. The “symmetrization” $\frac{1}{2}(X+X^{\!\top})$ of a non-symmetric matrix $X\in\mathcal{M}_{n}$ can be interpreted as “arithmetic” centering (with respect to $\mathcal{M}_{n}^{\nu}$). 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”.

Definition 1.2** (Geometric and Algebraic Centering of a Matrix).**

Let $0\prec X\in\mathcal{S}_{n}$ and $M_{\xi}$ a basis of $\mathcal{M}_{n}^{\xi}\subsetneq\mathcal{S}_{n}$ with $\dim\mathcal{M}_{n}^{\xi}=m_{\xi}\geq 1$. If $\alpha_{\max}^{+}(X,X_{\xi})$ and $\alpha_{\max}^{-}(X,X_{\xi})$ are both finite for all $X_{\xi}\in M_{\xi}$ then

[TABLE]

(respectively $X_{\text{cen}}=X$ for $m_{\xi}=0$) is called a geometric center of $X$ with respect to $M_{\xi}$. Any $X^{\prime}$ with

[TABLE]

is called an algebraic center of $X$ with respect to $\mathcal{M}_{n}^{\xi}$.

Example 1.3**.**

Let $n=2$ and

[TABLE]

Then $\alpha_{\max}=0.1$ and $X+\alpha_{\max}\Delta X$ is singular. However

[TABLE]

while \bigl{\lVert}(X+\alpha_{\max}\Delta X-0.9X_{\xi})^{-1}\bigr{\rVert}\approx 1.5713. Now let $k=10$ and $\alpha=(1-10^{-k})\alpha_{\max}$. Then $\beta^{+}=\alpha_{\max}^{+}(X,X_{\xi})\approx 10^{-11}$ and $\beta^{-}=\alpha_{\max}^{-}(X,X_{\xi})\approx 1.8$. Let $\beta=\frac{1}{2}(\beta^{+}-\beta^{-})$. Now the geometric center is $X_{\text{cen}}\approx X+\alpha\Delta X+\beta X_{\xi}$.

Remark 1.4**.**

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^{\prime}$, that is, $\mu X+(1-\mu)X^{\prime}$ for $0\leq\mu\leq 1$. This can be accomplished easily by solving a linear system of equations (assuming $\lVert X^{-1}\rVert\gg 1$ and $\lVert X\rVert\approx 1$). By “$\otimes$” we denote the Kronecker tensor product, for example

[TABLE]

For $Z\succ 0$ the matrix $I\otimes Z$ has full rank. Let $\tilde{Z}=Z^{-1}$. Given a basis $M_{\xi}$ of $\mathcal{M}_{n}^{\xi}$ (with $m_{\xi}\geq 1$) and assuming $\lVert\tilde{Z}\rVert/\lVert Z\rVert\gg 1$, the least squares solution of

[TABLE]

(using for example Lapack/DGELS) yields the “approximate” algebraic center $Z^{\prime}=Z+\operatorname{mat}M_{\xi}\hat{x}$ with $\lVert(Z^{\prime})^{-1}\rVert\ll\lVert Z^{-1}\rVert$ for $\mu=0$ (and $Z^{\prime}=Z$ for $\mu=1$).

Remark 1.6**.**

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).

Example 1.7**.**

For Example 1.3 we set $Z=X+\alpha\Delta X$ and get

[TABLE]

with $\lVert(X^{\prime})^{-1}\rVert\approx 1.6008$.

Remark 1.8**.**

How to solve (over- and) underdetermined linear systems of equations is discussed in detail in [Dem97, Chapter 3]. Given a linear system $Ax=b$, a least squares solution $x$ satisfies

[TABLE]

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

A set $\mathcal{K}\subseteq\mathcal{M}_{n}(\mathbb{R})$ is called convex if $\mu X+(1-\mu)Z$ holds for all $X,Z\in\mathcal{K}$ and $0<\mu<1$. A linear equality $Ax=b$ with $b\in\mathbb{R}$ (and $A\in\mathbb{R}^{1\times n}$) defines a hyperplane in $\mathbb{R}^{n}$, the corresponding linear inequality $Ax\leq b$ a halfspace in $\mathbb{R}^{n}$. The intersection of several halfspaces given by $Ax\leq b$ and $b\in\mathbb{R}^{m\times 1}$ for $m>1$ is again convex. Therefore the intersection of $\mathcal{K}$ and $Ax\leq b$ is convex.

2 Linear Programming

We briefly recall the basics: Let $1\leq m<n$ be positive integers. Given a (real) matrix $A\in\mathbb{R}^{m\times n}$ (with rows $A_{1},\ldots,A_{m}$) and two vectors $b\in\mathbb{R}^{m\times 1}$ and $c\in\mathbb{R}^{n\times 1}$, we want to find a vector $x\in\mathbb{R}^{n\times 1}$ and the minimum of the linear objective function $f_{0}=c^{\!\top}x$ such that the entries of $x$ are non-negative —written as $x\geq 0$— and the linear constraints $f_{i}=A_{i}x=b_{i}$ for $i\in\{1,2,\ldots,m\}$ —written as linear system $Ax=b$— hold. This is called the primal problem

[TABLE]

The corresponding dual problem is

[TABLE]

For simplicity we assume that $\operatorname{rank}A=m$ and the problems are feasible, that is, there exists an $x$ (respectively $y$ and $z$) such that the minimum in $\eqref{eqn:nc.lp.primal}$ (respectively the maximum in $\eqref{eqn:nc.lp.dual}$) is attained.

The “classical” solution is via Simplex Method in finitely many steps [Dan16]. Primal-dual interior-point methods move inside, that is, $x_{j}>0$ and $Ax\leq b$ (respectively $z_{j}>0$ and $A^{\!\top}y+z\leq c$), a convex polyhedron. (See Figure 1 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^{\!\top}x-\mu(\log x_{1}+\ldots+\log x_{n})$ for the Lagrange multiplier $\mu>0$. From the necessary condition

[TABLE]

for an optimum we get the complementary condition $x_{j}z_{j}=\mu$ (for $j=1,2,\ldots,n$). Now, starting at the triple $(x,y,z)$, we find the new search direction $(\Delta x,\Delta y,\Delta z)$ by solving the linearized system of equations

[TABLE]

with $2n+m$ unknowns. Given a $0<\tau<1$, the steplength $\alpha\leq 1$ is chosen such that $x+\tau\alpha\Delta x>0$ and $z+\tau\alpha\Delta z>0$. To make this procedure work practically, $\mu$ must decrease (in each step) “in the right way”.

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

Taking the classical Lagrange multiplier approach resulting in $x_{j}z_{j}=\mu$ one can ask what would happen if we take $n$ independent $\mu_{j}$ which could go to zero with different “speed”. In principle the condition $x_{j}z_{j}=\mu>0$ just tells us that both components, $x_{j}$ and $z_{j}$, should be invertible. Focusing on the primal problem and introducing “dummy” variables $\tilde{x}_{j}$ for the respective inverses of $x_{j}$ we need to solve the linear system

[TABLE]

for some $\gamma>0$. (For the general primal-dual formulation see (3.7) in the following section.) Assuming feasibility of $x$ and a solution $\Delta x$ of (2.4) there is no restriction for the steplength as long as $x+\tau\alpha\Delta x$ is positive. Thus for a feasible (positive) starting vector $x$ one can use a greedy method. This is implemented in ncminlp in Section 7.2. Recall that we assume $\operatorname{rank}A=m$. By $M$ we denote the orthogonal complement of $A^{\!\top}$ in $\mathbb{R}^{n}$. A search direction can be computed by solving the linear system

[TABLE]

which is much smaller than (2.3). The search direction is $\Delta x=M\hat{x}$. And if (2.5) 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 (2.5) can be used to find an initial feasible vector $x_{\text{ini}}>0$ by starting with $x=\tilde{x}=[1,\ldots,1]^{\!\top}$ and solving

[TABLE]

to get a search direction $\Delta x$. If \alpha=1/\lambda_{n}\bigl{(}-\operatorname{\mathrm{diag}}\Delta x\,(\operatorname{\mathrm{diag}}x)^{-1}\bigr{)}>1 we are done by using $x_{\text{ini}}=x+\Delta x$. Otherwise we set $x:=x+\tau\alpha\Delta x$ and $\tilde{x}:=x^{-1}$ (meaning componentwise inverse respectively the inverse as diagonal matrix) for some $0<\tau<1$ and iterate.

Remark 2.6**.**

Table 1 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 7.2. For further information on Karmarkar’s algorithm [Kar84] in linear programming we refer to [ARVK89].

Example 2.7**.**

Let $b=[2,7,3]^{\!\top}$, $c=[-1,-2,0,0,0]^{\!\top}$ and

[TABLE]

Then $x_{\text{opt}}=[3,5,3,0,0]^{\!\top}$. 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_{\text{opt}}$.

octave:1> ncminex octave:2> x = ncminlp(A_lp, b_lp, c_lp)

NCMINLP Search STD Version 0.99 December 2018 (C) KS

Linear Program: n=5, m=3 1-tau=1.00e-10, tol=2.00e-09, maxit=30 ini alpha min(x) tr(cx) ||Ax-b|| ||x.^-1|| 1 1.80e+15 1.000e+00 -7.000000000e+00 2.55e-15 1.6e+00 cnt alpha tr(cdx) tr(cx) ||A*x-b|| ||x.^-1|| 0 0.00e+00 -1.000e-01 -7.000000000e+00 2.55e-15 1.6e+00 1 3.17e+01 -3.171e+00 -1.017140874e+01 2.64e-15 5.0e+09 2 2.83e+01 -2.829e+00 -1.300000000e+01 2.95e-15 8.7e+09 3 3.54e+10 -4.243e-10 -1.300000000e+01 2.55e-15 7.1e+19 -1.299999999994142e+01

x = 3.0000e+00 5.0000e+00 3.0000e+00 5.8587e-11 1.4143e-20

octave:3> norm(x-x_lp) ans = 7.1754e-11

3 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

[TABLE]

The corresponding dual problem is

[TABLE]

Using diagonal matrices $X=\operatorname{\mathrm{diag}}x$, $C=\operatorname{\mathrm{diag}}c$, $Z=\operatorname{\mathrm{diag}}z$ and (by abuse of notation) $A_{k}=\operatorname{\mathrm{diag}}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 [Vin12], for matrix inequalities (MI’s) to [CHS06].

Remark. The inner product $\langle C,X\rangle$ is also often denoted by $C\bullet 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, $\langle X,Z\rangle=0$, implies $XZ=0$ [Tod01]. Why? Does strong duality also imply $\operatorname{rank}X+\operatorname{rank}Z=n$? For an overview about the differences between linear and semidefinite programming see in particular [Ali95, Section 3.5]. For details we refer to [WW10].

Definition 3.3** **(Central Path

[TTT98, Section 3.1], [WSV00, Section 10.2]).

The set of solutions to the equations

[TABLE]

for all $\mu>0$ such that $X\succ 0$ and $Z\succ 0$ is called central path.

Proposition 3.4** **(Weak Duality

[Tod01, Proposition 2.1]).

If $X$ is feasible in (3.1) and $(y,Z)$ is feasible in (3.2), then

[TABLE]

Starting from the triple $(X,y,Z)$, a new search direction $(\Delta X,\Delta y,\Delta Z)$ can be computed by solving a linear system of equations similar to (3.5) with $2n^{2}+m$ unknowns:

[TABLE]

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

The last matrix equation of (3.5) comes from the linearization of the centering condition $XZ=\mu I$ (for $\mu>0$) [AHO98]:

[TABLE]

This condition implies that $Z$ is a scalar multiple of $X^{-1}$ 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}=\mu$ from (primal-dual interior-point methods for) linear programming we have not found much discussions about its implications except for the need for symmetrization in [Tod99, Tod01]. Even if strong duality holds, that is, the duality gap is zero,

[TABLE]

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 [WW10] 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 2.

When we go over from the analytic formulation (in particular of some barrier functions [Nem07, Section 3]) 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 $\tilde{X}$ such that $X\tilde{X}=\tilde{X}X=I$ [Sch18, Remark 4.13]. Since we are looking for an update $\Delta X$ such that $X+\Delta X$ is still invertible, we get the condition

[TABLE]

which we linearize (assuming “small” $\Delta X$ and $\Delta\tilde{X}$) and reformulate (assuming only $\Delta X\tilde{X}=\tilde{X}\Delta X$) as

[TABLE]

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

[TABLE]

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

Remark 3.8**.**

An ad hoc implementation applied to the problem from Section 6 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.

Remark 3.9**.**

It is a funny coincidence that the approach we tried to ensure invertibility of an (unknown) matrix $X$, namely introducing another matrix $\tilde{X}$ and add the equations $X\tilde{X}=I$, as an alternative to the classical $\det X\neq 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 [Sch18, Remark 4.13]. The important observation is the origin of $XZ=\mu I$ from the classical barrier $-\log\det X$ (respectively $\log\det Z$) [Nem07], [TTT98, Section 2]. So it is rather natural to carry that over directly to our discrete (iterative) setting as “algebraic” barrier $X\tilde{X}=I$ (respectively $Z\tilde{Z}=I$). Recall that $\det X>0$ does not imply positive eigenvalues $\lambda_{i}>0$.

4 Feasible-interior-point Methods

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

[TABLE]

Given $m+1$ square matrices $C,A_{i}\in\mathcal{M}_{n}(\mathbb{R})=\mathbb{R}^{n\times n}$ and the vector $b\in\mathbb{R}^{m\times n}$ we want to find the minimum of $\langle C,X\rangle=\operatorname{\mathrm{tr}}(C^{\!\top}\!X)$ such that $X$ is symmetric and positive semidefinite and the $m$ constraints $\langle A_{i},X\rangle=b_{i}$ are fulfilled. We write $x=\operatorname{vec}X$ and $A^{\!\top}=[\operatorname{vec}A_{1},\operatorname{vec}A_{2},\ldots,\operatorname{vec}A_{m}]$ and use Example 6.1 (from the following section) for illustration:

[TABLE]

The vector space of $n\times n$ matrices $\mathcal{M}_{n}=\mathcal{M}_{n}(\mathbb{R})$ decomposes (almost) naturally in degrees of freedom $\mathcal{M}_{n}^{\xi}$, (pure) minimization directions $\mathcal{M}_{n}^{\eta}$, constraint directions $\mathcal{M}_{n}^{\zeta}$ and non-symmetric directions $\mathcal{M}_{n}^{\nu}$, namely

[TABLE]

with $m_{\nu}=n(n-1)/2$ and $m_{\xi}+m_{\eta}+m+m_{\nu}=n^{2}$ and bases $M^{\xi}$, $M^{\eta}$, $M^{\zeta}$ and $M^{\nu}$ respectively.

Remark 4.1**.**

However it is not so clear, how to choose $M^{\eta}$. In some sense $m_{\eta}=\operatorname{rank}M^{\eta}$ should be maximal for $(M^{\eta})^{\!\top}$ 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 [JCK08] $m_{\xi}=0$ (see also Table 2).

Figure 3 shows the cone of symmetric positive definite matrices $X=X^{\!\top}\succeq 0$ intersected with the “hyperplane” $Ax=b$ for $x=\operatorname{vec}X$. For a feasible matrix $X_{k}=X_{k}^{\!\top}\succ 0$, a (symmetric) direction $\Delta X_{k}\in\mathcal{M}_{n}^{\xi}\oplus\mathcal{M}_{n}^{\eta}$ and a scalar $\alpha_{k}\in\mathbb{R}$, the matrix $X_{k+1}:=X_{k}+\alpha_{k}\Delta X_{k}$ is again feasible, that is

[TABLE]

If $\Delta X_{k}$ is non-symmetric, that is, $\Delta X_{k}\in\mathcal{M}_{n}^{\xi}\oplus\mathcal{M}_{n}^{\eta}\oplus\mathcal{M}_{n}^{\nu}$, a symmetric direction is given by \Delta X_{k}:=\frac{1}{2}\bigr{(}\Delta X_{k}+(\Delta X_{k})^{\!\top}\bigl{)}. Before we discuss how to find such search directions $\Delta X_{k}$ (given a feasible positive definite $X_{k}$), we have a look on how to compute the maximal steplength $\alpha_{k}^{\text{max}}$ such that $X_{k}+\alpha_{k}^{\text{max}}\Delta X_{k}\succeq 0$. (Without loss of generality we assume that $\Delta X_{k}$ is such that $\alpha_{k}^{\text{max}}<\infty$.)

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=LL^{\!\top}$. The maximal steplength is (1.1), \alpha_{\max}=\bigl{(}\lambda_{\text{max}}(-L^{-1}\Delta XL^{-\!\top})\bigl{)}^{-1} [AHO98, Section 2]. So for $0<\tau<1$ we have

[TABLE]

Now we take a (not necessarily orthonormal) basis $M^{\xi}$ (respectively $M^{\eta}$ and $M^{\nu}$) for $\mathcal{M}_{n}^{\xi}$ (respectively $\mathcal{M}_{n}^{\eta}$ and $\mathcal{M}_{n}^{\nu}$) and write $M^{\xi,\eta}$ (respectively $M^{\xi,\eta,\nu}$) for $[M^{\xi},M^{\eta}]$ (respectively $[M^{\xi},M^{\eta},M^{\nu}]$) or simply just $M$ (for $M^{\xi,\eta}$ and $M^{\xi,\eta,\nu}$).

Remark. If the basis $M^{\xi}$ is not orthogonal one might be a little bit more careful with respect to geometric centering, see Definition 1.2 and the following remark.

Remark. Steplength scaling can be quite subtle in practice. For LP’s, tpically $\tau=0.99995$ is used which is too aggressive for “predictor-corrector” methods [TTT98, Section 4.4] therefore $\tau=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 3 we recall the linearized condition (3.6) for the invertibility of $X+\Delta X$, namely

[TABLE]

for given $X$ and $\tilde{X}=X^{-1}$. Since we start with a feasible $X$ and restrict the search directions to $\Delta X\in\operatorname{span}M$, this invertibility condition now reads

[TABLE]

for $\hat{x}\in\mathbb{R}^{m_{\xi}+m_{\eta}}$ or $\hat{x}\in\mathbb{R}^{m_{\xi}+m_{\eta}+m_{\nu}}$. Additionally, since we want to minimize, we need a condition of the form $\langle C,\Delta X\rangle=\langle C,M\hat{x}\rangle<0$. Thus, a minimization search direction can be computed as a least squares solution of the (underdetermined) linear system of equations

[TABLE]

for $\gamma>0$. Now we can use (1.1) to compute the steplength $\alpha_{\max}$ and iterate, at least in principal. Notice that the search direction computed in (4.3) is independent of $\gamma>0$ (modulo scaling).

Remark. 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 5 (cyan line): There is a first huge step, but then it seems that the path “converges” to some $X^{\prime}$ such that $\langle C,X^{\prime}\rangle\approx 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.

Remark 4.4**.**

For a case where centering is not needed (or even must not be used) —and how to adapt the implementation from Section 7— see (the end of) Section 8. 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\leq\mu\leq 0.85$ and $X_{\text{new}}:=\mu X+(1-\mu)X_{\text{cen}}$. This corresponds to typ=1 in line 7 in ncminsdp (Section 7.5) for algebraic centering.

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

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 $\mu\approx 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^{\xi,\eta}$ and center mainly for “staying close to the center” by using a centering parameter $0\leq\mu\leq 0.25$. This corresponds to typ=2 in line 7 in ncminsdp (Section 7.5) for algebraic centering.

Some paths for different feasible initial matrices are shown in Figure 6. For a problem run (in the context of the relaxation of a combinatorial problem) without centering see (the end of) Section 8. For a comparison of the number of iterations for the example from [JCK08, Section 6] see Table 2 (page 2).

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

Constructed Directions (ACE or GCE)

If degrees of freedom are available, that is, $m_{\xi}\geq 1$, search directions can be constructed by using two different centered matrices $X_{1}$ and $X_{2}$ with $\langle C,X_{2}\rangle<\langle C,X_{1}\rangle$ and setting $\Delta X:=X_{2}-X_{1}$.

For some paths using “approximate” algebraic centering see Figure 7, for some using “exact” geometric centering Figure 8. An investigation for $m_{\xi}\geq 2$ is open. For a typical problem run (and some further comments) see Section 6. Notice however, that the usability in the case of several minimization directions, that is, $m_{\eta}\geq 2$, is not yet clear.

5 Finding a feasible initial matrix

Even in “classical” interior-point methods (for semidefinite programming) the role of “good” initial iterates should not be underestimated [TTT99, Section 4]. For SDPA a “sufficiently large” multiple of the identity matrix should be used [FFK*+*08], the default parameter is $\lambda^{*}=100$. Since “$\lambda$” is heavily used for eigenvalues, we use “$\gamma$” here, for example $X_{0}=\gamma I$. Unfortunately we have not found that much beyond short remarks in the literature. For linear programming one could start with [ARVK89, Section 4]. 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 7.4. The variant for linear programs in Section 7.2 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 [Var04].

Still, starting with a least squares solution to $Ax=b$ and using repeatedly a spectral shift $\gamma I$ and a correction A\Delta x=b-A\bigl{(}x+\operatorname{vec}(\gamma I)\bigr{)} seems to work for practical problems. For small parameters $\varepsilon$ and $\delta$ in the example [JCK08, Section 6] one can increase the shift $\gamma$ dynamically.

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

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

[TABLE]

Notice that a scaling of the first block row might be necessary depending on the norm of $X^{-1}$. Since $\Delta X$ is not necessarily symmetric, we set \Delta X:=\frac{1}{2}\bigl{(}\Delta X+(\Delta X)^{\!\top}\bigr{)}. (In the implementation it is always done for numerical reasons.) Now we can calculate $\alpha=\alpha_{\max}(X,\Delta X)$ and update $X:=X+\min\{1,\tau\alpha\}\Delta X$ for some $\tau<1$. Secondly, for numerical stability we use either algebraic centering (1.5) if $m_{\xi}\geq 1$ or a modified version by computing a least squares solution $\Delta X=\operatorname{mat}(M\hat{x})$ for $M=[M_{\eta},M_{\nu}]$ to the linear system of equations

[TABLE]

for $0\leq\mu\leq 1$. In both cases we update X:=X+\frac{1}{2}\bigl{(}\Delta X+(\Delta X)^{\!\top}\bigr{)}. Notice that $\tilde{X}=X^{-1}$ needs to be updated in between.

Using $M=M_{\eta}$ only in the last step did not work for the example from [JCK08]. Some numbers for the initial iterates (for the presented approach) are in Table 2. For the problem from the next section see also Table 3. Notice that in this case there is a degree of freedom for “approximated” algebraic centering.

6 Sum of Squares

Suppose, that we want to find the global minimum of the polynomial $p(x)=2+\frac{13}{4}x^{2}+\frac{15}{4}x^{3}+x^{4}$ (its graph is the thick line in Figure 10). What is easy (in this case) using high-school mathematics will turn out to be rather involved using sum-of-squares and semidefinite programming [Las01, PS01, HP06], although the main idea is so simple: Minimize $\eta$ such that $p_{\eta}(x)=\eta+\frac{13}{4}x^{2}+\frac{15}{4}x^{3}+x^{4}$ can be written as

[TABLE]

for some $n_{\eta}\in\mathbb{N}$ and polynomials $q_{\eta,i}\in\mathbb{R}[x]$.

Remark. Not every non-negative polynomial can be written as a sum of squares [Ble10]. However, an “approximation” might be possible [Las07, LTY17, Par03].

Example 6.1**.**

Let $p_{\eta}(x)=\eta+\frac{13}{4}x^{2}+\frac{15}{4}x^{3}+x^{4}\in\mathbb{R}[x]$. It can be written as

[TABLE]

For $\eta=2$, and (approximately) $0.30775<\xi<1.39437$, the matrix $Q(\xi,2)$ is positive definite. For $\eta=1$, the matrix $Q(0.5,1)$ is positive semidefinite, and for $\eta<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

[TABLE]

with

[TABLE]

We write $x=\operatorname{vec}X$ and denote by $A\in\mathbb{R}^{m\times n^{2}}$ ($m=4$, $n=3$) the matrix with rows $(\operatorname{vec}A_{i})^{\!\top}$ to be able to write $A\operatorname{vec}(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 $\xi=0.75$ and $\eta=2$ is shown in Figure 8 (blue line). Table 3 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 5. 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^{-1}$ and the advantage of centered search directions for the (approximate) solution.

octave:1> ncminex octave:2> X = ncminsdp(A_sdp, b_sdp, C_sdp, Q(0.75,2));

NCMINSDP Search GCE Version 0.99 December 2018 (C) KS

Semidefinite Program: n=3, m=4, m_dof=1, m_min=1 1-tau=5.00e-05, tol=2.00e-08, mu=0.000, maxit=20 cnt alpha tr(CdX) tr(CX) ||A*vec(X)-b|| ||X^-1|| 0 0.00e+00 -1.000e-01 2.000000000e+00 0.00e+00 2.5e+01 1 9.88e+00 -9.876e-01 1.012402018e+00 2.40e-16 1.7e+03 2 1.26e-02 -1.240e-02 1.000002493e+00 5.06e-16 8.4e+06 3 2.01e-04 -2.493e-06 1.000000000e+00 5.06e-16 1.7e+11 4 5.00e-05 -1.247e-10 1.000000000e+00 9.21e-16 3.9e+15 1.000000000000002e+00 octave:3> norm(X-X_sdp) ans = 2.0251e-15

octave:4> X = ncminsdp(A_sdp, b_sdp, C_sdp, Q(0.75,100));

NCMINSDP Search GCE Version 0.99 December 2018 (C) KS

Semidefinite Program: n=3, m=4, m_dof=1, m_min=1 1-tau=5.00e-05, tol=2.00e-08, mu=0.000, maxit=20 cnt alpha tr(CdX) tr(CX) ||A*vec(X)-b|| ||X^-1|| 0 0.00e+00 -1.000e-01 1.000000000e+02 0.00e+00 1.0e+01 1 9.88e+02 -9.878e+01 1.219011619e+00 3.60e-14 9.9e+01 2 2.04e-03 -2.010e-01 1.017994134e+00 3.64e-14 1.2e+03 3 8.95e-02 -1.799e-02 1.000001178e+00 3.64e-14 1.8e+07 4 6.55e-05 -1.178e-06 1.000000000e+00 3.64e-14 3.6e+11 5 5.00e-05 -5.892e-11 1.000000000e+00 3.70e-14 5.9e+15 1.000000000000041e+00 octave:5> norm(X-X_sdp) ans = 4.8726e-14

octave:6> X = ncminsdp(A_sdp, b_sdp, C_sdp)

NCMINSDP Search STD Version 0.99 December 2018 (C) KS

Semidefinite Program: n=3, m=4, m_dof=1, m_min=1 1-tau=1.00e-06, tol=2.00e-08, mu=0.765, maxit=20 ini alpha min(eig(X)) tr(CX) ||Avec(X)-b|| ||X^-1|| 1 5.76e-01 4.247e-01 1.000000000e+00 1.86e+00 2.0e+03 2 9.13e-01 4.673e-02 1.205800795e+00 1.63e-01 4.4e+03 3 1.98e+00 4.037e-02 1.986977186e+00 5.15e-16 3.3e+01 cnt alpha tr(CdX) tr(CX) ||A*vec(X)-b|| ||X^-1|| 0 0.00e+00 -1.000e-01 1.986977186e+00 5.15e-16 2.5e+01 1 2.83e+00 -2.826e-01 1.704376627e+00 4.84e-16 9.4e+01 2 6.72e+00 -6.720e-01 1.032384907e+00 6.38e-16 1.4e+03 3 3.24e-01 -3.238e-02 1.000009826e+00 9.99e-16 4.3e+06 4 9.83e-05 -9.825e-06 1.000000000e+00 4.58e-16 1.1e+11 5 3.88e-09 -3.883e-10 1.000000000e+00 9.99e-16 1.1e+15 1.000000000000039e+00 X = 1.0000e+00 7.0203e-18 -2.5000e-01 7.0203e-18 3.7500e+00 1.8750e+00 -2.5000e-01 1.8750e+00 1.0000e+00

octave:7> norm(X-X_sdp,’fro’) ans = 7.1906-08

7 Implementation NCMIN

The following implementation in Octave [Oct18] is meant for educational purposes only without any warranties. By changing the function mldivide to lsqminnorm it can be used also in Matlab [Mat18]. (Notice however that the output can be slightly different due to rounding, etc.) It consists of a simple standalone linear solver ncminlp with net less than 70 lines (see Section 7.2 for the code and Section 2 for a typical application) and a semidefinite solver ncminsdp (Section 7.5) building on ncmindof (Section 7.3) and ncminini (Section 7.4) with net less than 210 lines together. For a typical problem run of the latter see Section 6.

Remark. 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 handling 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 [MH15] helpful. For a general background in numerical linear algebra we refer to [Dem97]. One typical “inefficiency” is that the inverse of $X$ is computed twice, for example in lines 26 and 34 (Section 7.4) or lines 80 and 88 (Section 7.5). 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.)

Remark. 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 $\tau$ too “aggressive”?

Remark. 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 [WW10]? 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?

7.1 NCMIN Examples

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

7.2 NCMIN Linear Solver

For the theoretical setting see Section 2 with the main linear system (2.5) corresponding to “A_sys*x_sys=b_sys” 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 5 and computes iteratively a direction “towards” $Ax=b$ subject to the linearized invertibility condition

[TABLE]

until the maximal steplength $\alpha=\alpha_{\max}(x,\Delta 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 [Meh92]) works.

For a small linear (scheduling) program with $m=24$ and $n=80$ ncminlp needs 7 iterations for the initial iterate and 18 iterations for the minimization (with relative error $10^{-8}$). SeDuMi 1.30 [Stu99] 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 ncminlp 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 …

7.3 NCMIN Vector Space Decomposition

The basic theory is explained in Section 1. Here, M_dof corresponds to $M_{\xi}$, M_min to $M_{\eta}$ and M_nsy to $M_{\nu}$. 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_{\xi}$. Notice the limitations of this (naive) approach with respect to the use of dense matrices (in total) of size $n^{2}\times n^{2}$. Although this task is not directly related to “preprocessing” it is maybe worth to view it from a more general context [CSW13, MPRW09, Ram97, TW12].

7.4 NCMIN Initial Feasible Matrix

For details see Section 5. For Example 6.1 one can subtract Q(0,0) from a matrix $X\in\mathcal{S}_{3}$ to read off $\xi$ in entry $(2,2)$ because $M_{\xi}^{\!\top}=[0,0,-\frac{1}{2},0,1,0,-\frac{1}{2},0,0]^{\!\top}$ and $\eta$ in entry $(1,1)$ because $M_{\eta}^{\!\top}=[1,0,\ldots,0]^{\!\top}$. In the following example we have $\xi\approx 0.767$ and $\eta\approx 1.987$:

octave:1> ncminex octave:2> X_ini = ncminini(A_sdp, b_sdp, C_sdp, 5e-16); octave:3> X_ini - Q(0,0) ans = 1.9870e+00 -1.1829e-16 -3.8329e-01 -1.1829e-16 7.6659e-01 2.2204e-16 -3.8329e-01 2.2204e-16 -1.1102e-16

7.5 NCMIN Semidefinite Solver

For details see Section 4. For a “typical” path for search type 1 (all directions) see Figure 4, type 2 (symmetric directions) Figure 6, type 3 (algebraic centered) Figure 7 and type 4 (geometric centered) Figure 8. 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 1, in particular Example 1.3 and 1.7.

8 (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) …” [Vin12].

Its importance cannot be overestimated due to the many practical applications, in particular in engineering [Nem07, Section 4.3] 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\times 3$ matrices. And if there were not Halmos’ recommendation to stop [Hal70, Section 19], this paper would have never been finished (and therefore not exist).

One could start with the simple linear matrix pencil $A=A(x,y)\in\mathbb{R}^{2\times 2}$ from [Net11, Example 1.1.1] 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:

[TABLE]

Here the zeros are replaced by (lower) dots to emphasize the structure and the tensor product “$\otimes$” 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}+xA_{x}+yA_{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 (1.1) respectively lines 61ff or 86ff in ncminsdp (Section 7.5)? Which other objects can be represented in such a way (by $2\times 2$ matrices)? How many variables are needed at most for $3\times 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 [Net11].

Another possibility is to use the power of semidefinite programming for concrete problems. One application is discussed in Section 6. Another is the relaxation of hard combinatorial problems [Nem07, Section 4.3.2]. For an illustration we restrict ourselves to 0-1 variables and a quadratic functional $f:\mathbb{R}^{2}\to\mathbb{R}$. Let $C\in\mathbb{R}^{3\times 3}$,

[TABLE]

and define $\bar{x}=[1,x_{1},x_{2}]^{\!\top}$ and $\bar{X}=\bar{x}^{\!\top}\bar{x}$. Now the optimum from the (non-convex) combinatorial problem can be bounded from below by a semidefinite program:

[TABLE]

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

[TABLE]

Some trivial but crucial questions: What is the rank of $\bar{X}$? What that of $X\succ 0$ in the SDP (8)? And which linear algebraic concept(s) are useful to measure the (approximate) “numerical” rank of some $X\succ 0$?

The following works only for symmetric search directions (typ=2 in line 7) and deactivated centering (by uncommenting line 106 in Section 7.5 and line 45 in Section 7.4). There might be a simple explanation …

octave:1> ncminex; X = ncminsdp(A_cmb, b_cmb, C_cmb, X_ini)

NCMINSDP Search SYM Version 0.99 December 2018 (C) KS

Semidefinite Program: n=3, m=3, m_dof=1, m_min=2 1-tau=5.00e-05, tol=2.00e-08, mu=0.235, maxit=20 ini alpha min(eig(X)) tr(CX) ||Avec(X)-b|| ||X^-1|| 1 1.10e+00 8.932e-02 1.666666667e-01 5.55e-16 1.2e+01 cnt alpha tr(CdX) tr(CX) ||A*vec(X)-b|| ||X^-1|| 0 0.00e+00 -1.000e-01 1.666666667e-01 5.55e-16 1.2e+01 1 4.95e+00 -4.950e-01 -3.283080272e-01 5.90e-16 6.0e+04 2 1.70e+00 -1.702e-01 -4.985070477e-01 6.05e-16 1.8e+05 3 1.26e-04 -1.258e-05 -4.985196251e-01 6.05e-16 3.5e+09 4 7.03e+10 -6.267e-10 -4.985196257e-01 6.05e-16 7.1e+13 -4.985196257109498e-01 X = 1.0000e+00 9.8043e-04 9.9900e-01 9.8043e-04 9.8043e-04 1.9686e-03 9.9900e-01 1.9686e-03 9.9900e-01

9 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 [WNM12] or “hard” problem [WW10]. In any case one should keep Todd’s words from the abstract of [Tod01] in mind: “The most effective computational methods are not always provable efficient in theory, and vice versa.” Very interesting to read is also [FGW02, Section 1.1] 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.” [FGW02, Section 1.1]

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 [Nem07, Tod01, VB96, Wri05]. Although one should use the latest versions (and the current documentation) of the implementations, the underlying publications [TTT99, Stu99, YFK03] are in some sense timeless. More on interior-point methods, search directions, predictor-corrector methods, etc. can be found in [AHO98, HRVW96, NT08, PW00, TTT98, Toh02]. Around verification, exact solutions, high precision, etc. there are [HNSED16, JCK08, WNM12] and [Rum10, Section 14]. Around strong duality, regularization, preprocessing, hard and bad semidefinite programs, etc. one should have a look in [CSW13, MPRW09, Pat17, RTW97, TW12, WW10].

On one hand, semidefinite programming is a special case of conic programming [BF09, BV04, FGW02, FV99, Jar92, JS95, NT98, Nes12]. 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: [Ali95, BEGFB94, HL03, Mar03, Par03], [Nem07, Section 4]. Since semidefinite programming is tight together with numerics (least squares approximation, iterative linear solvers, finite precision arithmetics, etc.) one should mention at least [CMTH16, MH15, WNM12]. But semidefinite programming is also interwoven with semialgebraic geometry [HN10, Net11, HM12].

At a first glance a topic like “non-commutative sum of squares” [Hel02, CKP11] seems to be far away. But indeed it was one of the initial questions around non-commutative convexity [HMV06] 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 [MH15]. 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 $\mu\to 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 [Rum10, Section 14]. That this should not be underestimated is shown in particular in [WNM12].

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 [FFK*+*08] sparse format, etc. For the issues with respect to benchmarking one could start with [Mit03], for large scale programs with [FNYF07].

In addition, thorough testing is necessary, in particular with the problems from SDPLIB [Bor99] and the library [dKS09]. 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 [dK10] 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”).