## Second-Order Cone Programming Algorithm

### Definition of Second-Order Cone Programming

A second-order cone programming problem has the form

$$\underset{x}{\mathrm{min}}{f}^{T}x$$

subject to the constraints

$$\begin{array}{c}\Vert {A}_{\text{sc}}(i)\cdot x-{b}_{\text{sc}}(i)\Vert \le {d}_{\text{sc}}^{T}(i)\cdot x-\gamma (i)\\ A\cdot x\le b\\ \text{Aeq}\cdot x=\text{beq}\\ \text{lb}\le x\le \text{ub}.\end{array}$$

*f*, *x*, *b*,
*beq*, *lb*, and *ub* are
vectors, and *A* and *Aeq* are matrices. For each
*i*, the matrix
*A*_{sc}(*i*), the vectors
*b*_{sc}(*i*) and
*d*_{sc}(*i*), and the
scalar *γ*(*i*) are in a second-order cone
constraint that you create using `secondordercone`

.

In other words, the problem has a linear objective function and linear constraints, as well as a set of second-order cone constraints of the form $$\Vert {A}_{\text{sc}}(i)\cdot x-{b}_{\text{sc}}(i)\Vert \le {d}_{\text{sc}}^{T}(i)\cdot x-\gamma (i)$$.

`coneprog`

Algorithm

The `coneprog`

solver uses the algorithm described in Andersen, Roos, and
Terlaky [1]. This method is
an interior-point algorithm similar to the Interior-Point linprog Algorithm.

#### Standard Form

The algorithm starts by placing the problem in *standard
form*. The algorithm adds nonnegative slack variables so that the
problem has the form

$$\underset{x}{\mathrm{min}}{f}^{T}x$$

subject to the constraints

$$\begin{array}{l}A\cdot x=b\\ x\in K.\end{array}$$

The solver expands the sizes of the linear coefficient vector
*f* and linear constraint matrix *A* to
account for the slack variables.

The region *K* is the cross product of *Lorentz
cones*
Equation 1 and the
nonnegative orthant. To convert each convex cone

$$\Vert {A}_{\text{sc}}(i)\cdot x-{b}_{\text{sc}}(i)\Vert \le {d}_{\text{sc}}^{T}(i)\cdot x-\gamma (i)$$

to a Lorentz cone Equation 1, create a
column vector of variables *t*_{1},
*t*_{2}, …,
*t*_{n+1}:

$$\begin{array}{c}{t}_{1}={d}^{T}x-\gamma \\ {t}_{2:(n+1)}={A}_{\text{sc}}x-{b}_{\text{sc}}.\end{array}$$

Here, the number of variables *n* for each cone
*i* is the number of rows in
*A*_{sc}(*i*). By its
definition, the variable vector *t* satisfies the
inequality

$$\Vert {t}_{2:(n+1)}\Vert \le {t}_{1}.$$ | (1) |

Equation 1 is the
definition of a Lorentz cone in (*n*+1) variables. The
variables *t* appear in the problem in place of the variables
*x* in the convex region *K*.

Internally, the algorithm also uses a *rotated Lorentz
cone* in the reformulation of cone constraints, but this topic
does not address that case. For details, see Andersen, Roos, and Terlaky [1].

When adding slack variables, the algorithm negates variables, as needed, and adds appropriate constants so that:

Variables with only one bound have a lower bound of zero.

Variables with two bounds have a lower bound of zero and, using a slack variable, have no upper bound.

Variables without bounds are placed in a Lorentz cone with a slack variable as the constrained variable. This slack variable is not part of any other expression, objective or constraint.

#### Dual Problem

The dual cone is

$${K}_{*}=\left\{s\text{\hspace{0.17em}}:\text{\hspace{0.17em}}{s}^{T}x\ge 0\text{\hspace{0.17em}}\forall x\in K\right\}.$$

The dual problem is

$$\underset{y}{\mathrm{max}}{b}^{T}y$$

such that

$${A}^{T}y+s=f$$

for some

$$s\in {K}_{*}.$$

A dual optimal solution is a point (*y*,*s*)
that satisfies the dual constraints and maximizes the dual objective.

#### Homogeneous Self-Dual Formulation

To handle potentially infeasible or unbounded problems, the algorithm adds two
more variables *τ* and *κ* and formulates the
problem as homogeneous (equal to zero) and self-dual.

$$\begin{array}{c}Ax-b\tau =0\\ {A}^{T}y+s-f\tau =0\\ -{f}^{T}x+{b}^{T}y-\kappa =0\end{array}$$ | (2) |

along with the constraints

$$(x;\tau )\in \overline{K},\text{\hspace{1em}}(s;\kappa )\in {\overline{K}}_{*}\text{\hspace{0.17em}}.$$ | (3) |

Here, $$\overline{K}$$ is the cone *K* adjoined with the nonnegative
real line, which is the space for (*x*;*τ*).
Similarly $${\overline{K}}_{*}$$ is the cone $${K}_{*}$$ adjoined with the nonnegative real line, which is the space
for (*s*;*κ*). In this formulation, the
following lemma shows that *τ* is the scaling for feasible
solutions, and *κ* is the indicator of an infeasible
problem.

**Lemma** ([1] Lemma
2.1)

Let (*x*, *τ*, *y*,
*s*, *κ*) be a feasible solution of Equation 2 along with
the constraints in Equation 3.

*x*+^{T}s*τκ*= 0.If

*τ*> 0, then (*x*,*y*,*s*)/*τ*is a primal-dual optimal solution of the standard form second-order cone problem.If

*κ*> 0, then at least one of these strict inequalities holds:*b*> 0^{T}y*f*< 0.^{T}xIf the first inequality holds, then the standard form, primal second-order cone problem is infeasible. If the second inequality holds, then the standard form, dual second-order cone problem is infeasible.

In summary, for feasible problems, the variable *τ* scales
the solution between the original standard form problem and the homogeneous
self-dual problem. For infeasible problems, the final iterate
(*x*, *y*, *s*,
*τ*, *κ*) provides a certificate of
infeasibility for the original standard form problem.

#### Start Point

The start point for the iterations is the feasible point:

*x*= 1 for each nonnegative variable, 1 for the first variable in each Lorentz cone, and 0 otherwise.*y*= 0.*s*= (1,0,…,0) for each cone, 1 for each nonnegative variable.*τ*= 1.*κ*= 1.

#### Central Path

The algorithm attempts to follow the *central path*,
which is the parameterized solution to the following equations for
*γ* decreasing from 1 toward 0.

$$\begin{array}{c}Ax-b\tau =\gamma \left(A{x}_{0}-b{\tau}_{0}\right)\\ {A}^{T}y+s-c\tau =\gamma \left({A}^{T}{y}_{0}+{s}_{0}-f{\tau}_{0}\right)\\ -{f}^{T}x+{b}^{T}y-\kappa =\gamma \left(-{f}^{T}{x}_{0}+{b}^{T}{y}_{0}-{\kappa}_{0}\right)\\ XSe=\gamma {\mu}_{0}e\\ \tau \kappa =\gamma {\mu}_{0}.\end{array}$$ | (4) |

Each variable with a 0 subscript indicates the start point of the variable.

The variables

*X*and*S*are*arrow head*matrices formed from the*x*and*s*vectors, respectively. For a vector*x*= [*x*_{1},*x*_{2},…,*x*], the arrow head matrix_{n}*X*has the definition$$X=\text{mat}(x)=\left[\begin{array}{cc}{x}_{1}& {x}_{2:n}^{T}\\ {x}_{2:n}& {x}_{1}I\end{array}\right].$$

By its definition,

*X*is symmetric.The variable

*e*is the vector with a 1 in each cone coordinate corresponding to the*x*_{1}Lorentz cone coordinate.The variable

*μ*_{0}has the definition$${\mu}_{0}=\frac{{x}_{0}^{T}{s}_{0}+{\tau}_{0}{\kappa}_{0}}{k+1},$$

where

*k*is the number of nonzero elements in*x*_{0}.

The central path begins at the start point and ends at an optimal solution to the homogeneous self-dual problem.

Andersen, Roos, and Terlaky [1] show in
Lemma 3.1 that the complementarity condition
*x ^{T}s* = 0, where

*x*and

*s*are in a product of Lorentz cones

*L*, is equivalent to the condition

$${X}_{i}{S}_{i}{e}_{i}={S}_{i}{X}_{i}{e}_{i}=0$$

for every cone *i*. Here
*X _{i}* =
mat(

*x*),

_{i}*x*is the variable associated with the Lorentz cone

_{i}*i*,

*S*= mat(

_{i}*s*), and

_{i}*e*is the unit vector [1,0,0,…,0] of the appropriate dimension. This discussion shows that the central path satisfies the complementarity condition at its end point.

_{i}#### Search Direction

To obtain points near the central path as the parameter *γ*
decreases from 1 toward 0, the algorithm uses Newton's method. The variables to
find are labeled (*x*, *τ*,
*y*, *s*, *κ*). Let
*d _{x}* represent the search
direction for the

*x*variables, and so on. Then the Newton step solves the following linear system, derived from Equation 4.

$$\begin{array}{c}A{d}_{x}-b{d}_{\tau}=(\gamma -1)\left(A{x}_{0}-b{\tau}_{0}\right)\\ {A}^{T}{d}_{y}+{d}_{s}-f{d}_{\tau}=(\gamma -1)\left({A}^{T}{y}_{0}+{s}_{0}-f{\tau}_{0}\right)\\ -{f}^{T}{d}_{x}+{b}^{T}{d}_{y}-{d}_{\kappa}=(\gamma -1)\left(-{f}^{T}{x}_{0}+{b}^{T}{y}_{0}-\kappa \right)\\ {X}_{0}{d}_{s}+{S}_{0}{d}_{x}=-{X}_{0}{S}_{0}e+\gamma {\mu}_{0}e\\ {\tau}_{0}{d}_{\kappa}+{\kappa}_{0}{d}_{t}au=-{\tau}_{0}{\kappa}_{0}+\gamma {\mu}_{0}.\end{array}$$

The algorithm obtains its next point by taking a step in the
*d* direction.

$$\left[\begin{array}{c}{x}_{1}\\ {\tau}_{1}\\ {y}_{1}\\ \begin{array}{l}{s}_{1}\\ {\kappa}_{1}\end{array}\end{array}\right]=\left[\begin{array}{c}{x}_{0}\\ {\tau}_{0}\\ {y}_{0}\\ \begin{array}{l}{s}_{0}\\ {\kappa}_{0}\end{array}\end{array}\right]+\alpha \left[\begin{array}{c}{d}_{x}\\ {d}_{\tau}\\ {d}_{y}\\ \begin{array}{l}{d}_{s}\\ {d}_{\kappa}\end{array}\end{array}\right]$$

for some step $$\alpha \in [0,1]$$.

For both numerical stability and accelerated convergence, the algorithm scales the step according to a suggestion in Nesterov and Todd [8]. Also, the algorithm corrects the step according to a variant of Mehrotra's predictor-corrector [7]. (For further details, see Andersen, Roos, and Terlaky [1].)

#### Step Solver Variations

The preceding discussion relates to the `LinearSolver`

option
with the value `'augmented'`

specified. The solver has other
values that change the step calculation to suit different types of
problems.

`'auto'`

(default) —`coneprog`

chooses the step solver:If the problem is sparse, the step solver is

`'prodchol'`

.Otherwise, the step solver is

`'augmented'`

.

`'normal'`

— The solver uses a variant of the`'augmented'`

step that is suitable when the problem is sparse. See Andersen, Roos, and Terlaky [1].`'schur'`

— The solver uses a modified Schur complement method for handling a sparse problem with a few dense columns. This method is also suitable for large cones. See Andersen [2].`'prodchol'`

— The solver uses the methods described in Goldfarb and Scheinberg ([4] and [5]) for handling a sparse problem with a few dense columns. This method is also suitable for large cones.

#### Iterative Display and Stopping Conditions

At each iteration *k*, the algorithm computes three relative
convergence measures:

Primal infeasibility

$${\text{Infeas}}_{\text{Primal}}^{k}=\frac{\Vert A{x}_{k}-b{\tau}_{k}\Vert}{\mathrm{max}\left(1,\Vert A{x}_{0}-b{\tau}_{0}\Vert \right)}.$$

Dual infeasibility

$${\text{Infeas}}_{\text{Dual}}^{k}=\frac{\Vert {A}^{T}{y}_{k}+{s}_{k}-f{\tau}_{k}\Vert}{\mathrm{max}\left(1,\Vert {A}^{T}{y}_{0}+{s}_{0}-f{\tau}_{0}\Vert \right)}.$$

Gap infeasibility

$${\text{Infeas}}_{\text{Gap}}^{k}=\frac{\left|-{f}^{T}{x}_{k}+{b}^{T}{y}_{k}-{\kappa}_{k}\right|}{\mathrm{max}\left(1,\left|-{f}^{T}{x}_{0}+{b}^{T}{y}_{0}-{\kappa}_{0}\right|\right)}.$$

You can view these three statistics at the command line by specifying iterative display.

options = optimoptions('coneprog','Display','iter');

All three should approach zero when the problem is feasible and the solver
converges. For a feasible problem, the variable
*κ _{k}* approaches zero, and the
variable

*τ*approaches a positive constant.

_{k}One stopping condition is somewhat related to the gap infeasibility. The stopping condition is when the following optimality measure decreases below the optimality tolerance.

$${\text{Optimality}}^{k}=\frac{\left|{f}^{T}{x}_{k}-{b}^{T}{y}_{k}\right|}{{\tau}_{k}+\left|{b}^{T}{y}_{k}\right|}=\frac{\left|{f}^{T}{x}_{k}/{\tau}_{k}-{b}^{T}{y}_{k}/{\tau}_{k}\right|}{1+\left|{b}^{T}{y}_{k}/{\tau}_{k}\right|}.$$

This statistic measures the precision of the objective value.

The solver also stops and declares the problem to be infeasible under the
following conditions. The three relative infeasibility measures are less than
*c* = `ConstraintTolerance`

, and

$${\tau}_{k}\le c\text{\hspace{0.17em}}\mathrm{max}\left(1,{\kappa}_{k}\right).$$

If *b ^{T}y_{k}*
> 0, then the solver declares that the primal problem is
infeasible. If

*f*< 0, then the solver declares that the dual problem is infeasible.

^{T}x_{k}The algorithm also stops when

$${\mu}_{k}\le c{\mu}_{0}$$

and

$${\tau}_{k}\le c\text{\hspace{0.17em}}\mathrm{max}\left(1,{\kappa}_{k}\right).$$

In this case, `coneprog`

reports that the problem is
numerically unstable (exit flag `-10`

).

The remaining stopping condition occurs when at least one infeasibility
measure is greater than `ConstraintTolerance`

and the computed
step size is too small. In this case, `coneprog`

reports that
the search direction became too small and no further progress could be made
(exit flag `-7`

).

## References

[1] Andersen, E. D., C. Roos,
and T. Terlaky. *On implementing a primal-dual interior-point method for
conic quadratic optimization.* Math. Program., Ser. B **95**, pp. 249–277 (2003). https://doi.org/10.1007/s10107-002-0349-3

[2] Andersen, K. D.
*A modified schur-complement method for handling dense columns in
interior-point methods for linear programming.* ACM Transactions on
Mathematical Software (TOMS), 22(3):348–356, 1996.

[3] Ben-Tal, Aharon, and
Arkadi Nemirovski. *Convex Optimization in Engineering: Modeling,
Analysis, Algorithms.* (1998). Available at https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.455.2733&rep=rep1&type=pdf.

[4] Goldfarb, D. and K.
Scheinberg. *A product-form cholesky factorization method for handling
dense columns in interior point methods for linear programming.*
Mathematical Programming, 99(1):1–34, 2004.

[5] Goldfarb, D. and K.
Scheinberg. *Product-form cholesky factorization in interior point methods
for second-order cone programming.* Mathematical Programming,
103(1):153–179, 2005.

[6] Luo, Zhi-Quan, Jos F.
Sturm, and Shuzhong Zhang. *Duality and Self-Duality for Conic Convex
Programming.* (1996). Available at https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.48.6432

[7] Mehrotra, Sanjay. “On the
Implementation of a Primal-Dual Interior Point Method.” *SIAM Journal on
Optimization 2*, no. 4 (November 1992): 575–601. https://doi.org/10.1137/0802028.

[8] Nesterov, Yu. E., and M.
J. Todd. “Self-Scaled Barriers and Interior-Point Methods for Convex Programming.”
*Mathematics of Operations Research 22*, no. 1 (February
1997): 1–42. https://doi.org/10.1287/moor.22.1.1.