Second-Order Cone Programming Algorithm

Definition of Second-Order Cone Programming

A second-order cone programming problem has the form

$\min_{x} f^{T} x$

subject to the constraints

$\begin{matrix} ‖ A_{sc} (i) \cdot x - b_{sc} (i) ‖ \leq d_{sc}^{T} (i) \cdot x - γ (i) \\ A \cdot x \leq b \\ Aeq \cdot x = beq \\ lb \leq x \leq ub . \end{matrix}$

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 $‖ A_{sc} (i) \cdot x - b_{sc} (i) ‖ \leq d_{sc}^{T} (i) \cdot x - γ (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

$\min_{x} 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

$‖ A_{sc} (i) \cdot x - b_{sc} (i) ‖ \leq d_{sc}^{T} (i) \cdot x - γ (i)$

to a Lorentz cone Equation 1, create a column vector of variables t₁, t₂, …, t_n+1:

$\begin{matrix} t_{1} = d^{T} x - γ \\ t_{2 : (n + 1)} = A_{sc} x - b_{sc} . \end{matrix}$

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

‖ t_{2 : (n + 1)} ‖ \leq 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_{*} = {s : s^{T} x \geq 0 \forall x \in K} .$

The dual problem is

$\max_{y} 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{matrix} A x - b τ = 0 \\ A^{T} y + s - f τ = 0 \\ - f^{T} x + b^{T} y - κ = 0 \end{matrix}

(2)

along with the constraints

(x; τ) \in \bar{K}, (s; κ) \in {\bar{K}}_{*} .

(3)

Here, $\bar{K}$ is the cone K adjoined with the nonnegative real line, which is the space for (x;τ). Similarly ${\bar{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^Ts + τκ = 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^Ty > 0
f^Tx < 0.
If 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{matrix} A x - b τ = γ (A x_{0} - b τ_{0}) \\ A^{T} y + s - c τ = γ (A^{T} y_{0} + s_{0} - f τ_{0}) \\ - f^{T} x + b^{T} y - κ = γ (- f^{T} x_{0} + b^{T} y_{0} - κ_{0}) \\ X S e = γ μ_{0} e \\ τ κ = γ μ_{0} . \end{matrix}

(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₁,x₂,…,x_n], the arrow head matrix X has the definition
$X = mat (x) = [\begin{matrix} x_{1} & x_{2 : n}^{T} \\ x_{2 : n} & x_{1} I \end{matrix}] .$
By its definition, X is symmetric.
The variable e is the vector with a 1 in each cone coordinate corresponding to the x₁ Lorentz cone coordinate.
The variable μ₀ has the definition
$μ_{0} = \frac{x_{0}^{T} s_{0} + τ_{0} κ_{0}}{k + 1},$
where k is the number of nonzero elements in x₀.

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^Ts = 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_i is the variable associated with the Lorentz cone i, S_i = mat(s_i), and e_i 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.

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{matrix} A d_{x} - b d_{τ} = (γ - 1) (A x_{0} - b τ_{0}) \\ A^{T} d_{y} + d_{s} - f d_{τ} = (γ - 1) (A^{T} y_{0} + s_{0} - f τ_{0}) \\ - f^{T} d_{x} + b^{T} d_{y} - d_{κ} = (γ - 1) (- f^{T} x_{0} + b^{T} y_{0} - κ) \\ X_{0} d_{s} + S_{0} d_{x} = - X_{0} S_{0} e + γ μ_{0} e \\ τ_{0} d_{κ} + κ_{0} d_{t} a u = - τ_{0} κ_{0} + γ μ_{0} . \end{matrix}$

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

$[\begin{matrix} x_{1} \\ τ_{1} \\ y_{1} \\ \begin{array}{l} s_{1} \\ κ_{1} \end{array} \end{matrix}] = [\begin{matrix} x_{0} \\ τ_{0} \\ y_{0} \\ \begin{array}{l} s_{0} \\ κ_{0} \end{array} \end{matrix}] + α [\begin{matrix} d_{x} \\ d_{τ} \\ d_{y} \\ \begin{array}{l} d_{s} \\ d_{κ} \end{array} \end{matrix}]$

for some step $α \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
${Infeas}_{Primal}^{k} = \frac{‖ A x_{k} - b τ_{k} ‖}{\max (1, ‖ A x_{0} - b τ_{0} ‖)} .$
Dual infeasibility
${Infeas}_{Dual}^{k} = \frac{‖ A^{T} y_{k} + s_{k} - f τ_{k} ‖}{\max (1, ‖ A^{T} y_{0} + s_{0} - f τ_{0} ‖)} .$
Gap infeasibility
${Infeas}_{Gap}^{k} = \frac{| - f^{T} x_{k} + b^{T} y_{k} - κ_{k} |}{\max (1, | - f^{T} x_{0} + b^{T} y_{0} - κ_{0} |)} .$

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 τ_k approaches a positive constant.

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

${Optimality}^{k} = \frac{| f^{T} x_{k} - b^{T} y_{k} |}{τ_{k} + | b^{T} y_{k} |} = \frac{| f^{T} x_{k} / τ_{k} - b^{T} y_{k} / τ_{k} |}{1 + | b^{T} y_{k} / τ_{k} |} .$

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

$τ_{k} \leq c \max (1, κ_{k}) .$

If b^Ty_k > 0, then the solver declares that the primal problem is infeasible. If f^Tx_k < 0, then the solver declares that the dual problem is infeasible.

The algorithm also stops when

$μ_{k} \leq c μ_{0}$

and

$τ_{k} \leq c \max (1, κ_{k}) .$

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

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

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