Semidefinite programming

A linear programming problem is one in which we wish to maximize or minimize a linear objective function of real variables over a polytope. In semidefinite programming, we instead use real-valued vectors and are allowed to take the dot product of vectors; nonnegativity constraints on real variables in LP (linear programming) are replaced by semidefiniteness constraints on matrix variables in SDP (semidefinite programming). Specifically, a general semidefinite programming problem can be defined as any mathematical programming problem of the form

${\displaystyle {\begin{array}{rl}{\displaystyle \min _{x^{1},\ldots ,x^{n}\in \mathbb {R} ^{n}}}&{\displaystyle \sum _{i,j\in [n]}c_{i,j}(x^{i}\cdot x^{j})}\\{\text{subject to}}&{\displaystyle \sum _{i,j\in [n]}a_{i,j,k}(x^{i}\cdot x^{j})\leq b_{k}}{\text{ for all }}k\\\end{array}}}$

where the ${\displaystyle c_{i,j},a_{i,j,k}}$, and the ${\displaystyle b_{k}}$ are real numbers and ${\displaystyle x^{i}\cdot x^{j}}$ is the dot product of ${\displaystyle x^{i}}$ and ${\displaystyle x^{j}}$.

An ${\displaystyle n\times n}$ matrix ${\displaystyle M}$ is said to be positive semidefinite if it is the Gramian matrix of some vectors (i.e. if there exist vectors ${\displaystyle x^{1},\ldots ,x^{n}}$ such that ${\displaystyle m_{i,j}=x^{i}\cdot x^{j}}$ for all ${\displaystyle i,j}$). If this is the case, we denote this as ${\displaystyle M\succeq 0}$. Note that there are several other equivalent definitions of being positive semidefinite, for example, positive semidefinite matrices are self-adjoint matrices that have only non-negative eigenvalues.

Denote by ${\displaystyle \mathbb {S} ^{n}}$ the space of all ${\displaystyle n\times n}$ real symmetric matrices. The space is equipped with the inner product (where ${\displaystyle {\rm {tr}}}$ denotes the trace) ${\displaystyle \langle A,B\rangle _{\mathbb {S} ^{n}}={\rm {tr}}(A^{T}B)=\sum _{i=1,j=1}^{n}A_{ij}B_{ij}.}$

We can rewrite the mathematical program given in the previous section equivalently as

${\displaystyle {\begin{array}{rl}{\displaystyle \min _{X\in \mathbb {S} ^{n}}}&\langle C,X\rangle _{\mathbb {S} ^{n}}\\{\text{subject to}}&\langle A_{k},X\rangle _{\mathbb {S} ^{n}}\leq b_{k},\quad k=1,\ldots ,m\\&X\succeq 0\end{array}}}$

where entry ${\displaystyle i,j}$ in ${\displaystyle C}$ is given by ${\displaystyle {\frac {c_{i,j}+c_{j,i}}{2}}}$ from the previous section and ${\displaystyle A_{k}}$ is a symmetric ${\displaystyle n\times n}$ matrix having ${\displaystyle i,j}$th entry ${\displaystyle {\frac {a_{i,j,k}+a_{j,i,k}}{2}}}$ from the previous section. Thus, the matrices ${\displaystyle C}$ and ${\displaystyle A_{k}}$ are symmetric and the above inner products are well-defined.

Note that if we add slack variables appropriately, this SDP can be converted to one of the form

${\displaystyle {\begin{array}{rl}{\displaystyle \min _{X\in \mathbb {S} ^{n}}}&\langle C,X\rangle _{\mathbb {S} ^{n}}\\{\text{subject to}}&\langle A_{k},X\rangle _{\mathbb {S} ^{n}}=b_{k},\quad k=1,\ldots ,m\\&X\succeq 0.\end{array}}}$

For convenience, an SDP may be specified in a slightly different, but equivalent form. For example, linear expressions involving nonnegative scalar variables may be added to the program specification. This remains an SDP because each variable can be incorporated into the matrix ${\displaystyle X}$ as a diagonal entry (${\displaystyle X_{ii}}$ for some ${\displaystyle i}$). To ensure that ${\displaystyle X_{ii}\geq 0}$, constraints ${\displaystyle X_{ij}=0}$ can be added for all ${\displaystyle j\neq i}$. As another example, note that for any positive semidefinite matrix ${\displaystyle X}$, there exists a set of vectors ${\displaystyle \{v_{i}\}}$ such that the ${\displaystyle i}$, ${\displaystyle j}$ entry of ${\displaystyle X}$ is ${\displaystyle X_{ij}=(v_{i},v_{j})}$ the scalar product of ${\displaystyle v_{i}}$ and ${\displaystyle v_{j}}$. Therefore, SDPs are often formulated in terms of linear expressions on scalar products of vectors. Given the solution to the SDP in the standard form, the vectors ${\displaystyle \{v_{i}\}}$ can be recovered in ${\displaystyle O(n^{3})}$ time (e.g., by using an incomplete Cholesky decomposition of X).

Analogously to linear programming, given a general SDP of the form

${\displaystyle {\begin{array}{rl}{\displaystyle \min _{X\in \mathbb {S} ^{n}}}&\langle C,X\rangle _{\mathbb {S} ^{n}}\\{\text{subject to}}&\langle A_{i},X\rangle _{\mathbb {S} ^{n}}=b_{i},\quad i=1,\ldots ,m\\&X\succeq 0\end{array}}}$

(the primal problem or P-SDP), we define the dual semidefinite program (D-SDP) as

${\displaystyle {\begin{array}{rl}{\displaystyle \max _{y\in \mathbb {R} ^{m}}}&\langle b,y\rangle _{\mathbb {R} ^{m}}\\{\text{subject to}}&{\displaystyle \sum _{i=1}^{m}}y_{i}A_{i}\preceq C\end{array}}}$

where for any two matrices ${\displaystyle P}$ and ${\displaystyle Q}$, ${\displaystyle P\succeq Q}$ means ${\displaystyle P-Q\succeq 0}$.

The weak duality theorem states that the value of the primal SDP is at least the value of the dual SDP. Therefore, any feasible solution to the dual SDP lower-bounds the primal SDP value, and conversely, any feasible solution to the primal SDP upper-bounds the dual SDP value. This is because

${\displaystyle \langle C,X\rangle -\langle b,y\rangle =\langle C,X\rangle -\sum _{i=1}^{m}y_{i}b_{i}=\langle C,X\rangle -\sum _{i=1}^{m}y_{i}\langle A_{i},X\rangle =\langle C-\sum _{i=1}^{m}y_{i}A_{i},X\rangle \geq 0,}$

where the last inequality is because both matrices are positive semidefinite, and the result of this function is sometimes referred to as duality gap.

Under a condition known as Slater's condition, the value of the primal and dual SDPs are equal. This is known as strong duality. Unlike for linear programs, however, not every SDP satisfies strong duality; in general, the value of the dual SDP may lie strictly below the value of the primal.

(i) Suppose the primal problem (P-SDP) is bounded below and strictly feasible (i.e., there exists ${\displaystyle X_{0}\in \mathbb {S} ^{n},X_{0}\succ 0}$ such that ${\displaystyle \langle A_{i},X_{0}\rangle _{\mathbb {S} ^{n}}=b_{i}}$, ${\displaystyle i=1,\ldots ,m}$). Then there is an optimal solution ${\displaystyle y^{*}}$ to (D-SDP) and

${\displaystyle \langle C,X^{*}\rangle _{\mathbb {S} ^{n}}=\langle b,y^{*}\rangle _{\mathbb {R} ^{m}}.}$

(ii) Suppose the dual problem (D-SDP) is bounded above and strictly feasible (i.e., ${\displaystyle \sum _{i=1}^{m}(y_{0})_{i}A_{i}\prec C}$ for some ${\displaystyle y_{0}\in \mathbb {R} ^{m}}$). Then there is an optimal solution ${\displaystyle X^{*}}$ to (P-SDP) and the equality from (i) holds.

Consider three random variables ${\displaystyle A}$, ${\displaystyle B}$, and ${\displaystyle C}$. By definition, their correlation coefficients ${\displaystyle \rho _{AB},\ \rho _{AC},\rho _{BC}}$ are valid if and only if

${\displaystyle {\begin{pmatrix}1&\rho _{AB}&\rho _{AC}\\\rho _{AB}&1&\rho _{BC}\\\rho _{AC}&\rho _{BC}&1\end{pmatrix}}\succeq 0,}$

in which case this matrix is called the correlation matrix. Suppose that we know from some prior knowledge (empirical results of an experiment, for example) that ${\displaystyle -0.2\leq \rho _{AB}\leq -0.1}$ and ${\displaystyle 0.4\leq \rho _{BC}\leq 0.5}$. The problem of determining the smallest and largest values that ${\displaystyle \rho _{AC}\ }$ can take is given by:

minimize/maximize ${\displaystyle x_{13}}$

subject to

${\displaystyle -0.2\leq x_{12}\leq -0.1}$

${\displaystyle 0.4\leq x_{23}\leq 0.5}$

${\displaystyle {\begin{pmatrix}1&x_{12}&x_{13}\\x_{12}&1&x_{23}\\x_{13}&x_{23}&1\end{pmatrix}}\succeq 0.}$

We set ${\displaystyle \rho _{AB}=x_{12},\ \rho _{AC}=x_{13},\ \rho _{BC}=x_{23}}$ to obtain the answer. This can be formulated by an SDP. We handle the inequality constraints by augmenting the variable matrix and introducing slack variables, for example

${\displaystyle \mathrm {tr} \left(\left({\begin{array}{cccccc}0&1&0&0&0&0\\0&0&0&0&0&0\\0&0&0&0&0&0\\0&0&0&1&0&0\\0&0&0&0&0&0\\0&0&0&0&0&0\end{array}}\right)\cdot \left({\begin{array}{cccccc}1&x_{12}&x_{13}&0&0&0\\x_{12}&1&x_{23}&0&0&0\\x_{13}&x_{23}&1&0&0&0\\0&0&0&s_{1}&0&0\\0&0&0&0&s_{2}&0\\0&0&0&0&0&s_{3}\end{array}}\right)\right)=x_{12}+s_{1}=-0.1}$

Solving this SDP gives the minimum and maximum values of ${\displaystyle \rho _{AC}=x_{13}\ }$ as ${\displaystyle -0.978}$ and ${\displaystyle 0.872}$ respectively.

Consider the problem

minimize ${\displaystyle {\frac {(c^{T}x)^{2}}{d^{T}x}}}$

subject to ${\displaystyle Ax+b\geq 0}$

where we assume that ${\displaystyle d^{T}x>0}$ whenever ${\displaystyle Ax+b\geq 0}$.

Introducing an auxiliary variable ${\displaystyle t}$ the problem can be reformulated:

minimize ${\displaystyle t}$

subject to ${\displaystyle Ax+b\geq 0,\,{\frac {(c^{T}x)^{2}}{d^{T}x}}\leq t}$

In this formulation, the objective is a linear function of the variables ${\displaystyle x,t}$.

The first restriction can be written as

${\displaystyle {\textbf {diag}}(Ax+b)\geq 0}$

where the matrix ${\displaystyle {\textbf {diag}}(Ax+b)}$ is the square matrix with values in the diagonal equal to the elements of the vector ${\displaystyle Ax+b}$.

The second restriction can be written as

${\displaystyle td^{T}x-(c^{T}x)^{2}\geq 0}$

Defining ${\displaystyle D}$ as follows

${\displaystyle D=\left[{\begin{array}{cc}t&c^{T}x\\c^{T}x&d^{T}x\end{array}}\right]}$

We can use the theory of Schur Complements to see that

${\displaystyle D\succeq 0}$

(Boyd and Vandenberghe, 1996)

The semidefinite program associated with this problem is

minimize ${\displaystyle t}$

subject to ${\displaystyle \left[{\begin{array}{ccc}{\textbf {diag}}(Ax+b)&0&0\\0&t&c^{T}x\\0&c^{T}x&d^{T}x\end{array}}\right]\succeq 0}$

Semidefinite programs are important tools for developing approximation algorithms for NP-hard maximization problems. The first approximation algorithm based on an SDP is due to Michel Goemans and David P. Williamson (JACM, 1995). They studied the max cut problem: Given a graph G = (V, E), output a partition of the vertices V so as to maximize the number of edges crossing from one side to the other. This problem can be expressed as an integer quadratic program:

Maximize ${\displaystyle \sum _{(i,j)\in E}{\frac {1-v_{i}v_{j}}{2}},}$ such that each ${\displaystyle v_{i}\in \{1,-1\}}$.

Unless P = NP, we cannot solve this maximization problem efficiently. However, Goemans and Williamson observed a general three-step procedure for attacking this sort of problem:

1. Relax the integer quadratic program into an SDP.
2. Solve the SDP (to within an arbitrarily small additive error ${\displaystyle \epsilon }$).
3. Round the SDP solution to obtain an approximate solution to the original integer quadratic program.

For max cut, the most natural relaxation is

${\displaystyle \max \sum _{(i,j)\in E}{\frac {1-\langle v_{i},v_{j}\rangle }{2}},}$ such that ${\displaystyle \lVert v_{i}\rVert ^{2}=1}$, where the maximization is over vectors ${\displaystyle \{v_{i}\}}$ instead of integer scalars.

This is an SDP because the objective function and constraints are all linear functions of vector inner products. Solving the SDP gives a set of unit vectors in ${\displaystyle \mathbf {R^{n}} }$; since the vectors are not required to be collinear, the value of this relaxed program can only be higher than the value of the original quadratic integer program. Finally, a rounding procedure is needed to obtain a partition. Goemans and Williamson simply choose a uniformly random hyperplane through the origin and divide the vertices according to which side of the hyperplane the corresponding vectors lie. Straightforward analysis shows that this procedure achieves an expected approximation ratio (performance guarantee) of 0.87856 - ε. (The expected value of the cut is the sum over edges of the probability that the edge is cut, which is proportional to the angle ${\displaystyle \cos ^{-1}\langle v_{i},v_{j}\rangle }$ between the vectors at the endpoints of the edge over ${\displaystyle \pi }$. Comparing this probability to ${\displaystyle (1-\langle v_{i},v_{j}\rangle )/{2}}$, in expectation the ratio is always at least 0.87856.) Assuming the unique games conjecture, it can be shown that this approximation ratio is essentially optimal.

Since the original paper of Goemans and Williamson, SDPs have been applied to develop numerous approximation algorithms. Recently, Prasad Raghavendra has developed a general framework for constraint satisfaction problems based on the unique games conjecture.[1]

Most codes are based on interior point methods (CSDP, MOSEK, SeDuMi, SDPT3, DSDP, SDPA). Robust and efficient for general linear SDP problems. Restricted by the fact that the algorithms are second-order methods and need to store and factorize a large (and often dense) matrix.

First-order methods for conic optimization avoid computing, storing and factorizing a large Hessian matrix and scale to much larger problems than interior point methods, at some cost in accuracy. A first-order method is implemented in the Splitting Cone Solver (SCS).[3] Another first-order method is the alternating direction method of multipliers (ADMM).[4] This method requires in every step projection on the cone of semidefinite matrices.

The code ConicBundle formulates the SDP problem as a nonsmooth optimization problem and solves it by the Spectral Bundle method of nonsmooth optimization. This approach is very efficient for a special class of linear SDP problems.

Algorithms based on Augmented Lagrangian method (PENSDP) are similar in behavior to the interior point methods and can be specialized to some very large scale problems. Other algorithms use low-rank information and reformulation of the SDP as a nonlinear programming problem (SDPLR).[5]

Algorithms that solve SDPs approximately have been proposed as well. The main goal of such methods is to achieve lower complexity in applications where approximate solutions are sufficient and complexity must be minimal. A prominent method that has been used for data detection in multiple-input multiple-output (MIMO) wireless systems is Triangular Approximate SEmidefinite Relaxation (TASER),[6] which operates on the Cholesky decomposition factors of the semidefinite matrix instead of the semidefinite matrix. This method calculates approximate solutions for a max-cut-like problem that are often comparable to solutions from exact solvers but in only 10-20 algorithm iterations.