A Quadratically Convergent Sequential Programming Method for Second-Order Cone Programs Capable of Warm Starts

While a large number of interior-point algorithms for SOCP have been proposed, including some that have been implemented in efficient optimization packages [12, 15, 21], there are only very few approaches for solving SOCPs with an active-set framework. The method proposed by Goldberg and Leyffer [10] is a two-phase algorithm that combines a projected-gradient method with equality-constrained SQP. However, it is limited to instances that have only conic constraints (1c) and no additional linear constraints (1b). Hayashi et al. [13] propose a simplex-type method, where they reformulate the SOCP as a linear semi-infinite program to handle the fact that these instances have infinitely many extreme points. The resulting dual-simplex exchange method shows promising practical behavior. However, in contrast to the method proposed here, the authors conjecture that their method has only an R-linear local convergence rate. Zhadan [30] proposes a similar simplex-type method. Another advantage of the method presented in this paper is that the pivoting algorithm does not need to be designed and implemented from scratch. Instead, it can leverage existing implementations of active-set QP solvers, in particular the efficient handling of linear systems.

The proposed algorithm relies on polyhedral outer approximations based on well-known cutting planes for SOCPs. For instance, the methods for mixed-integer SOCP by Drewes and Ulbrich [7] and Coey et al. [4] use these cutting planes to build LP relaxations of the branch-and-bound subproblems. We note that an LP-based cutting plane algorithm for SOCP could be seen as an active-set method, but it is only linearly convergent. As pointed out in [6], it is crucial to consider the curvature of the conic constraint in the subproblem objective to achieve fast convergence.

The term “SQP method for SOCP” has also been used in the literature to refer to methods for solving nonlinear SOCPs [6, 16, 24, 31]. However, in contrast to the method here, in these approaches, the subproblems themselves are SOCPs (1) and include the linearization of the nonlinear objective and constraints. It will be interesting to explore extensions of the proposed method to nonlinear SOCPs in which feasibility is achieved asymptotically not only for the nonlinear constraints but also for the conic constraints.

For two vectors $x,y\in\mathbb{R}^{n}$, we denote with $x\circ y$ their component-wise product, and the condition $x\perp y$ stands for $x^{T}y=0$. For $x\in\mathbb{R}^{n}$, we define $[x]^{+}$ as the vector with entries $\max\{x_{i},0\}$. We denote by $\|\cdot\|$, $\|\cdot\|_{1}$, $\|\cdot\|_{\infty}$ the Euclidean norm, the $\ell_{1}$-norm, and the $\ell_{\infty}$-norm, respectively. For a cone $\mathcal{K}_{j}$, $e_{ji}\in\mathbb{R}^{n_{j}}$ is the canonical basis vector with $1$ in the element corresponding to $x_{ji}$ for $i\in\{0,\ldots,n_{j}-1\}$, and $\operatorname{int}(\mathcal{K}_{j})$ and $\operatorname{bd}(\mathcal{K}_{j})$ denote the cone’s interior and boundary, respectively.

The definition of the second-order cone in (2) suggests that the conic constraint (1c) can be replaced by the nonlinear constraint

without changing the set of feasible points. Consequently, (1) is equivalent to

Unfortunately, (3) cannot be solved directly with standard gradient-based algorithms for nonlinear optimization, such as SQP methods. The reason is that $r_{j}$ is not differentiable whenever $\bar{x}_{j}=0$. This is particularly problematic when the optimal solution $x^{*}$ of the SOCP lies at the extreme point of a cone, $x^{*}_{j}=0\in\mathcal{K}_{j}$. In that case, the Karush-Kuhn-Tucker (KKT) necessary optimality conditions for the NLP formulation, which are expressed in terms of derivatives, cannot be satisfied. Therefore, any optimization algorithm that seeks KKT points cannot succeed. As a remedy, differentiable approximations of $r_{j}$ have been proposed in the past; see, for example, [28]. However, high accuracy comes at the price of high curvature, which can make finding the numerical solution of the NLP difficult.

An alternative equivalent reformulation of the conic constraint is given by

In this case, the constraint function is differentiable. But if $x_{j}^{*}=0$, its gradient vanishes, and as a consequence, no constraint qualification applies and the KKT conditions do not hold. Therefore, again, a gradient-based method cannot be employed. By using an outer approximation of the cones that is improved in the course of the algorithm, our proposed variation of the SQP method is able to avoid these kinds of degeneracy.

To facilitate the discussion we define a point-wise partition of the cones.

If the set $\mathcal{E}(x^{*})$ at an optimal solution $x^{*}$ were known in advance, we could compute $x^{*}$ as a solution of (1) by solving the NLP

The constraints involving the linearization of $r_{j}$ are imposed only if $r_{j}$ is differentiable at $x^{*}$, and variables in cones that are extremal-active at $x^{*}$ are explicitly fixed to zero. With this, locally around $x^{*}$, all functions in (4) are differentiable and we could apply standard second-order algorithms to achieve fast local convergence.

In (4), we omitted the cones in $\mathcal{N}(x^{*})$. If $x^{*}$ is feasible for the SOCP and $j\in\mathcal{N}(x^{*})$ we have $\bar{x}^{*}_{j}=0$ and $x^{*}_{j0}>0$, and so $r_{j}(x^{*})<0$. This implies that the nonlinear constraint (4c) for this cone is not active and we can omit it from the problem statement without impacting the optimal solution.

The proposed algorithm is designed to guide the iterates $x^{k}$ into the neighborhood of an optimal solution $x^{*}$. If the optimal solution is not degenerate and the iterates are sufficiently close to $x^{*}$, the steps generated by the algorithm are eventually identical to the steps that the SQP method would take for solving the differentiable optimization problem (4). In this section, we review the mechanisms and convergence results of the basic SQP method [23].

At an iterate $x^{k}$, the basic SQP method, applied to (4), computes a step $d^{k}$ as an optimal solution to the QP subproblem

Here, $H^{k}$ is the Hessian of the Lagrangian function for (4), which in our case is

where $\mu_{j}^{k}\geq 0$ are estimates of the optimal multipliers for the nonlinear constraint (4c), and where ${\rm NewA}^{2}_{xx}r_{j}(x_{j})$ is the $n\times n$ block-diagonal matrix with

in the rows and columns corresponding to $x_{j}$ for $j\in\mathcal{J}$. It is easy to see that ${\rm NewA}^{2}r_{j}(x_{j})$ is positive semi-definite. The estimates $\mu_{j}^{k}$ are updated based on the optimal multipliers $\hat{\mu}_{j}^{k}\geq 0$ corresponding to (5c).

Algorithm 1 formally states the basic SQP method where $\hat{\lambda}^{k}\geq 0$ and $\hat{\eta}^{k}$ denote the multipliers corresponding to (4b) and (4d), respectively. Because we are only interested in the behavior of the algorithm when $x^{k}$ is close to $x^{*}$, we assume here that $\bar{x}_{j}^{k}\neq 0$ for all $j\in\mathcal{D}(x^{*})$ and for all $k$, and hence the gradient and Hessian of $r_{j}$ can be computed. Note that the iterates $\hat{\lambda}^{k}$ and $\hat{\eta}^{k}$ are not explicitly needed in Algorithm 1, but they are necessary to measure the optimality error and define the primal-dual iterate sequence that is analyzed in Theorem 2.2.

A fast rate of convergence can be proven under the following sufficient second-order optimality assumptions [3].

Under these assumptions, the basic SQP algorithm reduces to Newton’s method applied to the optimality conditions of (4) and the following result holds [23].

Theorem 2.2 is a local convergence result. Practical SQP algorithms include mechanisms that make sure that the iterates eventually reach such a neighborhood, even if the starting point is far away. To this end, we employ the exact penalty function

in which $\rho>0$ is a penalty parameter. Note that we define $\varphi$ in terms of all conic constraints $\mathcal{J}$, even though $r_{j}$ appears in (4c) only for $j\in\mathcal{D}(x^{*})$. We do this because the proposed algorithm does not know $\mathcal{D}(x^{*})$ in advance and the violation of all cone constraints needs to be taken into account when the original problem (1) is solved. Nevertheless, in this section, we may safely ignore the terms for $j\not\in\mathcal{D}(x^{*})$ because for $j\in\mathcal{E}(x^{*})$ we have $x^{k}_{j}=0$ and hence $[r_{j}(x^{k})]^{+}=0$ for all $k$ due to (5d), and when $j\in\mathcal{N}(x^{*})$, we have $r_{j}(x_{j}^{k})<0$ when $x^{k}$ is close to $x^{*}$ since $r_{j}(x_{j}^{*})<0$.

It can be shown, under suitable assumptions, that the minimizers of $\varphi(\,\cdot\,;\rho)$ over the set defined by the linear constraints (4b),

coincide with the minimizers of (4) when $\rho$ is chosen sufficiently large. Because it is not known upfront how large $\rho$ needs to be, the algorithm uses an estimate, $\rho^{k}$, in iteration $k$, which might be increased during the course of the algorithm.

To ensure that the iterates eventually reach a minimizer of $\varphi(\,\cdot\,;\rho)$, and therefore a solution of (4), we require that the decrease of $\varphi(\,\cdot\,;\rho)$ is at least a fraction of that achieved in the piece-wise linear model of $\varphi(\,\cdot\,;\rho)$ given by

constructed at $x^{k}$. More precisely, the algorithm accepts a trial point $\hat{x}^{k+1}=x^{k}+d$ as a new iterate only if the sufficient decrease condition

holds with some fixed constant $c_{\text{dec}}\in(0,1)$. The trial iterate $\hat{x}^{k+1}=x^{k}+d^{k}$ with $d^{k}$ computed from (5) might not always satisfy this condition. The proposed algorithm generates a sequence of improved steps of which one is eventually accepted.

However, to apply Theorem 2.2, it would be necessary that the algorithm take the original step $d^{k}$ computed from (5); see Step 5 of Algorithm 1. Unfortunately, $\hat{x}^{k+1}=x^{k}+d^{k}$ might not be acceptable even when the iterate $x^{k}$ is arbitrarily close to a non-degenerate solution $x^{*}$ satisfying Assumption 1 (a phenomenon called the Maratos effect [19]). Our remedy is to employ the second-order correction step [9], $s^{k}$, which is obtained as an optimal solution of the QP

For later reference, let $\hat{\lambda}^{S,k}$, $\hat{\mu}^{S,k}$ and $\hat{\eta}^{S,k}$ denote optimal multiplier vectors corresponding to (12b)–(12d), respectively. The algorithm accepts the trial point $\hat{x}^{k+1}=x^{k}+d^{k}+s^{k}$ if it yields sufficient decrease (11) with respect to the original SQP step $d=d^{k}$. Note that (12) is a variation of the second-order correction that is usually used in SQP methods, for which (12c) reads

and avoids the evaluation of ${\rm NewA}r_{j}(x_{j}^{k}+d_{j}^{k})$. In our setting, however, evaluating ${\rm NewA}r_{j}(x_{j}^{k}+d_{j}^{k})$ takes no extra work and (12c) is equivalent to a supporting hyperplane, see Section 3.1. As the following theorem shows (see, e.g., [9] or [5, Section 15.3.2.3]), this procedure computes steps with sufficient decrease (11) and results in quadratic convergence.

The proposed algorithm aims at finding an optimal solution of the SOCP (1), or equivalently, values of the primal variables, $x^{*}\in\mathbb{R}^{n}$, and the dual variables, $\lambda^{*}\in\mathbb{R}^{m}$ and $z_{j}^{*}\in\mathbb{R}^{n_{j}}$ for $j\in\mathcal{J}$, that satisfy the necessary and sufficient optimality conditions [1, Theorem 16]

A thorough discussion of SOCPs is given in the comprehensive review by Alizadeh and Goldfarb [1]. The authors consider the formulation in which the linear constraints (1b) are equality constraints, but the results in [1] can be easily extended to inequalities.

The primal-dual solution $(x^{*},\lambda^{*},z^{*})$ is unique under the following assumption.

The definition of non-degeneracy for SOCP is somewhat involved and we refer the reader to [1, Theorem 21]. Strict complementarity holds if $x_{j}^{*}+z_{j}^{*}\in\operatorname{int}(\mathcal{K}_{j})$ and implies that: (i) $x^{*}_{j}\in\operatorname{int}(K_{j})\Longrightarrow z^{*}_{j}=0$; (ii) $z^{*}_{j}\in\operatorname{int}(K_{j})\Longrightarrow x^{*}_{j}=0$; (iii) $x^{*}_{j}\in\operatorname{bd}(K_{j})\setminus\{0\}\Longleftrightarrow z^{*}_{j}\in\operatorname{bd}(K_{j})\setminus\{0\}$; and (iv) not both $x^{*}_{j}$ and $z_{j}^{*}$ are zero.

In the following, consider a particular cone $\mathcal{K}_{j}$ and let $\mathcal{Y}_{j}$ be a finite subset of $\{y_{j}\in\mathbb{R}^{n_{j}}:\bar{y}_{j}\neq 0,y_{j0}\geq 0\}$. We define the cone

generated by the points in $\mathcal{Y}_{j}$. For each $x_{j}\in\mathcal{K}_{j}$ we have $r_{j}(x_{j})\leq 0$, and using

we obtain for any $y_{j}\in\mathcal{Y}_{j}$ that

Therefore $\mathcal{C}_{j}(\mathcal{Y}_{j})\supseteq\mathcal{K}_{j}$. Also, for $y_{j}\in\mathcal{Y}_{j}$, consider $x_{j}=(1,\bar{y}_{j}^{T}/\|\bar{y}_{j}\|)^{T}.$ Then

and also $r_{j}(x_{j})=\|\bar{x}_{j}\|-x_{j0}={\bar{y}_{j}}/{\|\bar{y}_{j}\|}-1=0.$ Hence $x_{j}\in\mathcal{C}_{j}(\mathcal{Y}_{j})\cap\mathcal{K}_{j}$. Therefore, for any $y_{j}\in\mathcal{Y}_{j}$, the inequality

defines a hyperplane that supports $\mathcal{K}_{j}$ at $(1,{\bar{y}_{j}}/{\|\bar{y}_{j}\|})$. In summary, $\mathcal{C}_{j}(\mathcal{Y}_{j})$ is a polyhedral outer approximation of $\mathcal{K}_{j}$, defined by supporting hyperplanes.

In addition, writing $\mathcal{Y}_{j}=\{y_{j,1},\ldots,y_{j,m}\}$, we also define the cone

For all $x_{j}\in\mathcal{C}_{j}(\mathcal{Y}_{j})$ and $z_{j}=-\sum_{l=1}^{m}\sigma_{j,l}{\rm NewA}r_{j}(y_{j,l})+\eta_{j}e_{j0}\in\mathcal{C}_{j}^{\circ}(\mathcal{Y}_{j})$, we have

because ${\rm NewA}r_{j}(y_{j,l})^{T}x_{j}\leq 0$ and $x_{j0}\geq 0$ from the definition of $\mathcal{C}_{j}(\mathcal{Y}_{j})$. Therefore $\mathcal{C}_{j}^{\circ}(\mathcal{Y}_{j})$ is included in the dual of the cone $\mathcal{C}_{j}(\mathcal{Y}_{j})$.

Now define $R=[-{\rm NewA}r_{j}(y_{j,1}),\ldots,-{\rm NewA}r_{j}(y_{j,m}),e_{j0}]$ and let $z_{j}\in\mathbb{R}^{n_{j}}$ be in the dual of $\mathcal{C}_{j}(\mathcal{Y}_{j})$. Since this implies that $x_{j}^{T}z_{j}\geq 0$ for all $x\in\mathcal{C}_{j}(\mathcal{Y}_{j})=\{\mathbb{R}^{n_{j}}:R^{T}x_{j}\geq 0\}$, Farkas’ lemma yields that $z_{j}=R\cdot(\sigma^{T},\eta)^{T}$ for some $\sigma_{j}\in\mathbb{R}^{m}_{+}$ and $\eta_{j}\geq 0$, i.e., $z_{j}\in\mathcal{C}_{j}^{\circ}(\mathcal{Y}_{j})$.

Overall we proved that $\mathcal{C}_{j}^{\circ}(\mathcal{Y}_{j})$ defined in (17) is the dual of $\mathcal{C}_{j}(\mathcal{Y}_{j})$, and since $\mathcal{C}_{j}(\mathcal{Y}_{j})\supseteq\mathcal{K}_{j}$, this implies $\mathcal{C}^{\circ}_{j}(\mathcal{Y}_{j})\subseteq\mathcal{K}_{j}$.

In each iteration, at an iterate $x^{k}$, the proposed algorithm computes a step $d^{k}$ as an optimal solution of the subproblem

Here, $H^{k}$ is a positive semi-definite matrix that captures the curvature of the nonlinear constraint (3c), and for each cone, $\mathcal{Y}^{k}_{j}$ is the set of hyperplane-generating points that have been accumulated up to this iteration. From (14), we see that (18d) can be replaced by linear constraints. Consequently, (18) is a QP and can be solved as such.

Algorithm 2 describes a preliminary version of the proposed SQP method based on this subproblem. Observe that the linearization (18c) can be rewritten as

and is equivalent to the hyperplane constraint generated at $x_{j}^{k}$. Consequently, if $x^{k}_{j}\not\in\mathcal{K}_{j}$, then $r_{j}(x^{k}_{j})>0$ and (18c) acts as a cutting plane that excludes $x^{k}_{j}$. Using the update rule

in Step 6 makes sure that $x^{k}_{j}$ is excluded in all future iterations.

In our algorithm, we initialize $\mathcal{Y}_{j}^{0}$ so that

In this way, $x_{j}=0$ is an extreme point of $\mathcal{C}_{j}(\mathcal{Y}_{j}^{0})$, as it is for $\mathcal{K}_{j}$, and the challenging aspect of the cone is already captured in the first subproblem. By choosing the coordinate vectors $e_{ji}$ we have ${\rm NewA}r_{j}(e_{ji})^{T}x_{j}=x_{ji}-x_{j0},$ and the hyperplane constraint (16) becomes a very sparse linear constraint.

When $H^{k}=0$ in each iteration, this procedure becomes the standard cutting plane algorithm for the SOCP (1). It is well-known that the cutting plane algorithm is convergent in the sense that every limit point of the iterates is an optimal solution of the SOCP (1), but the convergence is typically slow. In the following sections, we describe how Algorithm 2 is augmented to achieve fast local convergence. The full method is stated formally in Algorithm 3.

We now describe a strategy that enables our algorithm to identify those cones that are extreme-active at a non-degenerate solution $x^{*}$ within a finite number of iterations, i.e., $\mathcal{E}(x^{k})=\mathcal{E}(x^{*})$ for all large $k$. This will make it possible to apply a second-order method and achieve quadratic local convergence.

Consider the optimality conditions for the QP subproblem (18):

Here, $\hat{\lambda}^{k}$, $\hat{\mu}^{k}_{j}$, and $\hat{\nu}_{j}^{k}$ are the multipliers corresponding to the constraints in (18); for completeness, we define $\hat{\mu}_{j}^{k}=0$ for $j\in\mathcal{J}\setminus\mathcal{D}(x^{k})$. In (21a), ${\rm NewA}_{x}r_{j}(x^{k})$ is the vector in $\mathbb{R}^{n}$ that contains ${\rm NewA}r_{j}(x_{j}^{k})$ in the elements corresponding to $x_{j}$ and is zero otherwise. Similarly, $\hat{\nu}^{k}\in\mathbb{R}^{n}$ is equal to $\hat{\nu}^{k}_{j}$ in the elements corresponding to $x_{j}$ for all $j\in\mathcal{J}$ and zero otherwise.

Let us define

It is easy to verify that, for $j\in\mathcal{D}(x^{k})$, ${\rm NewA}r_{j}(x_{j}^{k})x_{j}^{k}=r_{j}(x_{j}^{k})$ and hence $r_{j}(x_{j}^{k})^{T}(x_{j}^{k}+d^{k})\leq 0$ from (21c). As a consequence we obtain $x^{k}_{j}+d_{j}^{k}\in\mathcal{C}_{j}(\hat{\mathcal{Y}}_{j}^{k})$ for all $j\in\mathcal{J}$. Furthermore, $\hat{\nu}_{j}^{k}\in\mathcal{C}^{\circ}_{j}(\mathcal{Y}_{j}^{k})$ implies that

for suitable values of $\sigma^{k}_{j,l}\geq 0$ and $\eta^{k}_{j}\geq 0$. Then $\hat{z}_{j}^{k}:=-\hat{\mu}_{j}^{k}\nabla r_{j}(x^{k})+\hat{\nu}_{j}^{k}\in\mathcal{C}^{\circ}_{j}(\hat{\mathcal{Y}}_{j}^{k})$ and

from (21a). In conclusion, if $(d,\hat{\lambda}^{k},\hat{\mu}^{k},\hat{\nu}^{k})$ is a primal-dual solution of the QP subproblem (18), then $(d,\hat{\lambda}^{k},\hat{z}^{k})$ satisfies the conditions

which more closely resembles the SOCP optimality conditions (13). Our algorithm maintains primal-dual iterates $(x^{k+1},\hat{\lambda}^{k},\hat{z}^{k})$ that are updated based on (24).

Suppose that strict-complementarity holds at a primal-dual solution $(x^{*},\lambda^{*},z^{*})$ of the SOCP (1) and that $(x^{k+1},\hat{\lambda}^{k},\hat{z}^{k})\to(x^{*},\lambda^{*},z^{*})$. If $j\not\in\mathcal{E}(x^{*})$ then $x_{j}^{*}\in\mathcal{K}_{j}$ implies $x^{*}_{j0}>0$. As $x^{k}_{j}$ converges to $x^{*}_{j}$, we have $x^{k}_{j0}>0$ and therefore $j\not\in\mathcal{E}(x^{k})$ for sufficiently large $k$. This yields $\mathcal{E}(x^{k})\subseteq\mathcal{E}(x^{*})$. We now derive a modification of Algorithm 2 that ensures that $\mathcal{E}(x^{*})\subseteq\mathcal{E}(x^{k})$ for all sufficiently large $k$ under Assumption 2.

Consider any $j\in\mathcal{E}(x^{*})$. We would like to have

for all large $k$, since then complementarity in (24c) implies that $x_{j}^{k+1}=x_{j}^{k}+d^{k}_{j}=0$ and hence $j\in\mathcal{E}(x^{k+1})$ for all large $k$. We will later show that Assumption 2 implies that $\hat{z}_{j}^{k}\to z^{*}_{j}$ and that there exists a neighborhood $N_{\epsilon}(z^{*}_{j})=\{z_{j}\in\mathbb{R}^{n_{j}}:\|z_{j}-z^{*}_{j}\|\leq\epsilon\}$ of $z^{*}_{j}$ so that $z_{j}\in\operatorname{int}(\mathcal{C}^{\circ}_{j}(\hat{\mathcal{Y}}_{j}^{0}\cup\{-y_{j}\}))$ if $z_{j},y_{j}\in N_{\epsilon}(z^{*}_{j})$; see Remark 4.18. This suggests that some vector close to $-z_{j}^{*}$ should eventually be included in $\hat{\mathcal{Y}}_{j}^{k}$ because then (25) holds when $\hat{z}_{j}^{k}$ is close enough to $z^{*}_{j}$. For this purpose, the algorithm computes