SuperSCS: fast and accurate large-scale conic optimization

TL;DR

SuperSCS introduces a novel, efficient method for large-scale convex conic problems by integrating SuperMann with Douglas-Rachford splitting on a self-dual embedding, enhancing speed and accuracy.

Contribution

It combines SuperMann and Douglas-Rachford splitting on a self-dual embedding to improve large-scale conic optimization solving capabilities.

Findings

Faster convergence on large-scale problems

High accuracy in solutions

Effective handling of infeasibility and unboundedness

Abstract

We present SuperSCS: a fast and accurate method for solving large-scale convex conic problems. SuperSCS combines the SuperMann algorithmic framework with the Douglas-Rachford splitting which is applied on the homogeneous self-dual embedding of conic optimization problems: a model for conic optimization problems which simultaneously encodes the optimality conditions and infeasibility/unboundedness certificates for the original problem. SuperMann allows the use of fast quasi-Newtonian directions such as a modified restarted Broyden-type direction and Anderson's acceleration.

Figures

Tables

Table 1. TABLE I: Regularized PCA SDP: Solver Statistics

Method	${\mathrm{sgm}}_{10}\left(t^s\right)$	Success
SCS	96.82	74.31%
RB (mem:50)	3.17	100%
RB (mem:100)	3.31	100%
AA (mem:5)	2.13	100%
AA (mem:10)	2.66	100%

Table 2. TABLE II: LASSO: Solver Statistics

Method	${\mathrm{sgm}}_{10}\left(t^s\right)$	Success
SCS	5.97	100%
RB (mem:50)	2.88	100%
RB (mem:100)	2.61	100%
AA (mem:5)	3.38	100%
AA (mem:10)	3.87	100%

Table 3. TABLE III: Sparse ℓ 1 subscript ℓ 1 \ell_{1} -regularized logistic regression: Solver Statistics

Method	${\mathrm{sgm}}_{10}\left(t^s\right)$	Success
SCS	4.85	100%
RB (mem:50)	7.23	100%
RB (mem:100)	7.28	100%
AA (mem:5)	3.00	100%
AA (mem:10)	3.09	100%

Table 4. TABLE IV: Maros-Mészáros QP problems: Solver Statistics

Method	${\mathrm{sgm}}_{10}\left(t^s\right)$	Success
SCS	56.61	83.02%
RB (mem:50)	9.66	90.57%
RB (mem:100)	6.57	90.57%
AA (mem:5)	5.79	90.57%
AA (mem:10)	8.62	91.51%

Equations

$$\operatorname{mat}(x)=\tfrac{1}{\sqrt{2}}\smallmat{\sqrt{2}x_{1}&x_{2}\cdots x_{n}\\ x_{2}\sqrt{2}x_{n+1}\cdots x_{2n-1}\\ \vdots\vdots\ddots\vdots\\ x_{n}x_{2n-1}\cdots\sqrt{2}x_{n(n+1)/2}},$$

$$\begin{array}[]{rl}\displaystyle{\operatorname*{minimize}_{x\in\R^{n}}}&\ \<c,x\\ \operatorname*{subject\ to}&\ b{}-{}Ax{}={}s,\ {}s{}\in{}\mathcal{K},\end{array}$$

$$\begin{array}[]{rl}\displaystyle{\operatorname*{minimize}_{y\in\R^{m}}}&\ \<b,y\\ \operatorname*{subject\ to}{}&\ y\in\mathcal{K}^{*},\ {}A^{{\scriptscriptstyle\top}}y{}+{}c{}={}0\end{array}$$

$$\begin{array}{r}A{x}^{\star }+s^{\star }=b,s^{\star }\in \mathcal{K},y^{\star }\in {\mathcal{K}}^{\ast },\\ A^{\top }{y}^{\star }+c=0,\text{<}{y}^{\star },s^{\star }=0.\end{array}$$

$$Q\text{dfn}\text{smallmat}$$

$$x^{\star }={\chi }^{\star }/{\tau }^{\star },y^{\star }={\psi }^{\star }/{\tau }^{\star },s^{\star }={\varsigma }^{\star }/{\tau }^{\star },$$

$$0\in Qu+N_{\mathcal{C}}(u),$$

$$N_{\mathcal{C}}(u)=$$

$$=$$

$${\tilde{u}}^{\nu }=$$

$${\bar{u}}^{\nu }=$$

$$u^{\nu +1}=$$

$$u^{\nu +1}=T{u}^{\nu },$$

$$u^{\nu +1}=(1-\lambda )u^{\nu }+\lambda T{u}^{\nu },$$

$$R=I-T.$$

$$w^{\nu }=$$

$$u^{\nu +1}=$$

$${\tilde{w}}^{\nu }=$$

$$=$$

$$\nu$$

$${\mathrm{pr}}_{\nu }=\frac{\Vert A{\bar{x}}^{\nu }+{\bar{s}}^{\nu }-b\Vert }{1+\Vert b\Vert }$$

$${\mathrm{dr}}_{\nu }=\frac{\Vert A^{\top }{\bar{y}}^{\nu }+c\Vert }{1+\Vert c\Vert }$$

$${\mathrm{gap}}_{\nu }=\frac{|\text{<}c,{\bar{x}}^{\nu }+\text{<}b,{\bar{y}}^{\nu }|}{1+|\text{<}c,{\bar{x}}^{\nu }|+|\text{<}b,{\bar{y}}^{\nu }|}$$

$${\mathrm{ic}}_{\nu }=\{\begin{array}{ll}\Vert b\Vert \Vert A^{\top }{\bar{y}}^{\nu }\Vert /\text{<}b,{\bar{y}}^{\nu },& \text{ if }\text{<}b,{\bar{y}}^{\nu }<0\\ +\infty ,& \text{ else}\end{array}$$

$${\mathrm{uc}}_{\nu }=\{\begin{array}{ll}\Vert c\Vert \Vert A{\bar{x}}^{\nu }+{\bar{s}}^{\nu }\Vert /\text{<}c,{\bar{x}}^{\nu },& \text{ if }\text{<}c,{\bar{x}}^{\nu }<0\\ +\infty ,& \text{ else}\end{array}$$

$$d^{\nu }=-B_{\nu }^{-1}R{u}^{\nu }=-H_{\nu }R{u}^{\nu },$$

$$B_{\nu +1}=B_{\nu }+\frac{1}{\Vert z^{\nu }{\Vert }^2}({\tilde{\xi }}^{\nu }-B_{\nu }{z}^{\nu })z^{\nu \top }$$

$${\theta }_{\nu }=\{\begin{array}{ll}1& \text{if }|{\gamma }_{\nu }|\ge \bar{\vartheta }\\ \frac{1-\mathrm{sgn}({\gamma }_{\nu })\bar{\vartheta }}{1-{\gamma }_{\nu }}& \text{otherwise}\end{array}$$

$$H_{\nu +1}=H_{\nu }+\frac{1}{\text{<}{H}_{\nu }{\tilde{\xi }}^{\nu },z^{\nu }}(z^{\nu }-H_{\nu }{\tilde{\xi }}^{\nu })(z^{\nu \top }{H}_{\nu }),$$

$$d^{\nu }=-R{u}^{\nu }-(Z_{\nu }-{\Xi }_{\nu })t^{\nu },$$

$$\underset{t^{\nu }}{\mathrm{minimize}}\Vert {\Xi }_{\nu }{t}^{\nu }-R{u}^{\nu }{\Vert }^2,$$

$$r_{p,s}=\frac{t_{p,s}}{{\min }_{s^{\prime }\in S}{t}_{p,s^{\prime }}}.$$

$${\rho }_s(\tau )=\text{nicefrac}1|P|\cdot |\{p\in P:r_{p,s}\le \tau \}|,$$

$$t_p^s=\{\begin{array}{ll}{t}_{p,s}& \text{ if }{t}_{p,s}<\infty \\ 100\max \{t_{p,s}\mid p\in P,t_{p,s}<\infty \}& \text{ otherwise}\end{array}$$

$$\operatorname{sgm}_{\sigma}{}={}\exp\Big{[}\textstyle\sum_{p{}\in{}P}\ln\left(\max\{1,\sigma{}+{}t^{s}_{p}\}\right)\Big{]}{}-{}\sigma.$$

Full text

I Introduction

Conic optimization problems are of central importance in convex optimization as several solvers and parsers such as CVX [11], CVXPy [6], YALMIP [14] and MOSEK [17] transform given problems into a conic representation. Indeed, all convex optimization problems can be cast in the standard form of a conic optimization problem.

Various interior point methods have been proposed for solving conic optimization problems [23, 25, 8]. Interior point methods can achieve high accuracy, yet, do not scale well with the problem size. On the other hand, first-order methods have low per-iteration cost and minimum memory requirements, therefore, are better suited for large-scale problems [18]. Recent research has turned to first-order methods for large-scale conic problems such as SDPs [28]. However, their convergence rate is at most Q-linear with a Q-factor often close to one, especially for ill-conditioned problems.

The KKT conditions of a conic optimization problem together with conditions for the detection of infeasibility or unboundedness can be combined in a convex feasibility problem known as the homogeneous self-dual embedding (HSDE) [27]. The HSDE has been used in both interior point [23] and first-order methods [18].

In this paper we present a numerical optimization method for solving the HSDE. The HSDE is first cast as a variational inequality which can be equivalently seen as a monotone inclusion. We observe that the splitting cone solver (SCS) presented in [18] can be interpreted as the application of the Douglas-Rachford splitting (DRS) to that monotone inclusion. We then apply the reverse splitting to that monotone inclusion — which is a firmly nonexpansive operator — and employ the SuperMann scheme [24] which allows the use of quasi-Newtonian directions such as restarted Broyden directions and Anderson’s acceleration [26, 9]. We call the resulting method SuperSCS. This way, SuperSCS can achieve a fast convergence rate while retaining a low per-iteration cost. In fact, SuperSCS uses exactly the same oracle as SCS.

II Mathematical preliminaries

We denote by $\R$, $\R_{{\scriptscriptstyle+}}$, $\R^{n}$ and $\R^{m\times n}$ the sets of real numbers, non-negative reals, $n$-dimensional real vectors and $m$-by-$n$ real matrices respectively. We denote the transpose of a matrix $A$ by $A^{{\scriptscriptstyle\top}}$. For two vectors $x,y\in\R^{n}$, we denote by $\<x,y\>=x^{{\scriptscriptstyle\top}}y$ their standard inner product. Let $E$ be a vector space in $\R^{n}$. We define the orthogonal complement of $E$ in $\R^{n}$ to be the vector space $E^{{\scriptscriptstyle\perp}}=\{y\in\R^{n}{}\mid{}\<y,x\>=0,\forall x\in E\}$.

A set $\mathcal{K}\subseteq\R^{n}$ is called a convex cone if it is convex and $\lambda x{}\in{}\mathcal{K}$ for every $x{}\in{}\mathcal{K}$ and $\lambda>0$. The binary relation $x\succcurlyeq_{\mathcal{K}}y$ is interpreted as $x-y\in\mathcal{K}$. The dual cone of $\mathcal{K}$ is defined as $\mathcal{K}^{*}=\{x^{*}{}\mid{}\<x^{*},x\>\geq 0,\forall x\in\mathcal{K}\}.$

A few examples of convex cones of interest are: (i) the zero cone $\mathcal{K}_{n}^{\mathrm{f}}=\{0\}^{n}$, (ii) the cone of symmetric positive semidefinite matrices $\mathcal{K}^{\mathrm{s}}_{n}{}={}\{x\in\R^{n(n+1)/2}{}\mid{}\operatorname{mat}(x):\text{pos. definite}\},$ where $\operatorname{mat}:\R^{n(n+1)/2}\to\R^{n\times n}$ is defined by

[TABLE]

(iii) the second-order cone $\mathcal{K}^{\mathrm{q}}_{n}{}={}\{z=(x,t){}:{}x\in\R^{n-1},t{}\in{}\R{}\mid{}\|x\|_{2}\leq t\},$ (iv) the positive orthant $\mathcal{K}^{\mathrm{l}}_{n}{}={}\{x\in\R^{n}{}\mid{}x\geq 0\}$ and (v) the three-dimensional exponential cone, $\mathcal{K}^{\mathrm{e}}{}={}\operatorname{cl}\{(x_{1},x_{2},x_{3}){}\mid{}x_{1}{}\geq{}x_{2}e^{x_{3}/x_{2}},{}x_{2}>0\}.$

The normal cone of a nonempty closed convex set $C$ is the set-valued mapping $N_{C}(x)=\{g{}\mid{}\<g,y-x\>\leq 0,\forall y\in C\}$ for $x\in C$ and $N_{C}(x)=\emptyset$ for $x\notin C$. The Euclidean projection of $x$ on $C$ is denoted by $\Pi_{C}$.

III Conic programs

A cone program is an optimization problem of the form

[TABLE]

where $A\in\R^{m\times n}$ is a possibly sparse matrix and $\mathcal{K}$ is a nonempty closed convex cone.

The vast majority of convex optimization problems of practical interest can be represented in the above form [18, 3]. Indeed, cone programs can be thought of as a universal representation for all convex problems of practical interest and many convex optimization solvers first transform the given problem into this form.

The dual of ($\mathcal{P}$) is given by [3, Sec. 1.4.3]

[TABLE]

Let $p^{\star}$ be the optimal value of ($\mathcal{P}$) and $d^{\star}$ be the optimal value of ($\mathcal{D}$). Strong duality holds ($p^{\star}=-d^{\star}$) if the primal or the dual problem are strictly feasible [3, Thm. 1.4.2].

Whenever strong duality holds, the following KKT conditions are necessary and sufficient for optimality of $(x^{\star},s^{\star},y^{\star})\in\R^{n}\times\R^{m}\times\R^{m}$:

[TABLE]

Infeasibility and unboundedness conditions are provided by the so-called theorems of the alternative. The weak theorems of the alternative state that 1) Either primal feasibility holds, or there is a $y$ with $A^{{\scriptscriptstyle\top}}y=0$, $y\succeq_{\mathcal{K}^{*}}0$ and $\<b,y\><0$, and, similarly, 2) Either dual feasibility holds or there is a $x$ so that $Ax\succeq_{\mathcal{K}}0$ and $\<c,x\>\leq 0$ [3, Sec. 1.4.7].

III-A Homogeneous self-dual embedding

In this section we present a key result which is due to Ye et al. [27]: the HSDE, which is a feasibility problem which simultaneously describes the optimality, (in)feasibility and (un)boundedness conditions of a conic optimization problem. Solving the HSDE yields a solution of the original conic optimization problem, when one exists, or a certificate of infeasibility or unboundedness. We start by considering the following feasibility problem in $(\chi,\varsigma,\psi,\tau,\kappa)$

[TABLE]

Note that for $\tau^{\star}=1$ and $\kappa^{\star}=0$, the above equations reduce to the primal-dual optimality conditions. As shown in [27], the solutions of (2) satisfy $\kappa^{\star}\tau^{\star}{}={}0$, i.e., at least one of $\kappa^{\star}$ and $\tau^{\star}$ must be zero. In particular, if $\kappa=0$ and $\tau>0$, then the triplet $(x^{\star},y^{\star},s^{\star})$ with

[TABLE]

is a primal-dual solution of ($\mathcal{P}$) and ($\mathcal{D}$). If instead $\tau^{\star}=0$ and $\kappa>0$, then the problem is either primal- or dual-infeasible. If $\tau{}={}\kappa{}={}0$, no conclusion can be drawn.

We define $u{}={}(\chi,\psi,\tau)$ and $v=(0,\varsigma,\kappa)$. The self-dual embedding reduces to the problem of determining $u$ and $v$ such that $Qu=v$ with $(u,v)\in\mathcal{C}\times\mathcal{C}^{*}$ where $\mathcal{C}{}\dfn{}\R^{n}{}\times{}\mathcal{K}^{*}{}\times{}\R_{+}.$ This is equivalent to the variational inequality

[TABLE]

Indeed, for all $u\in\mathcal{C}$,

[TABLE]

where the second equality follows by considering $v=\nicefrac{{1}}{{2}}u$ and $v=\nicefrac{{3}}{{2}}u$, which both belong to the cone $\mathcal{C}$ (see also [12, Ex. 5.2.6(a)]). From that and the fact that $Qu\in\{u\}^{{\scriptscriptstyle\perp}}$ since $Q$ is skew-symmetric, the equivalence of the HSDE in (2) and the variational inequality in (2c) follows. Equation (2c) is a monotone inclusion which can be solved using operator theory machinery as we discuss in the following section.

IV SuperSCS

IV-A SCS and DRS

Since $Q$ is a skew-symmetric linear operator, it is maximally monotone. Being the normal cone of a convex set, $N_{\mathcal{C}}$ is maximally monotone as well. Additionally, because of [2, Cor. 24.4(i)], $Q+N_{\mathcal{C}}$ is maximally monotone. Therefore, we may apply the Douglas-Rachford splitting on the monotone inclusion (2c). The SCS algorithm [18] is precisely the application of the DRS to $N_{\mathcal{C}}{}+{}Q$ leading to the iterations discussed in [21, Sec. 7.3]. This observation furnishes a short and elegant interpretation of SCS. Here, on the other hand, we consider the reverse splitting, $Q+N_{\mathcal{C}}$, which leads to the following DRS iterations

[TABLE]

For any initial guess $u^{0}$, the iterates $u^{\nu}$ converge to a point $u^{\star}$ which satisfies the monotone inclusion (2c) [2, Thm. 25.6(i), (iv)]. The linear system in (2da) can be either solved “directly” using a sparse LDL factorization or “indirectly” by means of the conjugate gradient method [18]. The projection on $\mathcal{C}$ in (2db) essentially requires that we be able to project on $\mathcal{K}^{*}$.

The iterative method (2d) can be concisely written as

[TABLE]

where $T:\R^{N}\to\R^{N}$ is given by $Tu=u+\Pi_{\mathcal{C}}(2(I+Q)^{-1}u-u)-(I+Q)^{-1}u$ and is firmly nonexpansive [2, Chap. 26]. As such it fits the Krasnosel’skii-Mann framework [2, Sec. 5.2] leading to the relaxed iterations

[TABLE]

with $\lambda\in(0,2)$ and, as a result, it fits the SuperMann framework [24].

IV-B SuperSCS: SuperMann meets SCS

SuperMann considers the problem of finding a fixed-point $x^{\star}{}\in{}\operatorname{fix}T$ from the viewpoint of finding a zero of the residual operator

[TABLE]

SuperMann, instead of applying Krasnosel’skii-Mann-type updates of the form (2f), takes extragradient-type updates of the general form

[TABLE]

where $d^{\nu}$ are fast, e.g., quasi-Newtonian, directions and scalar parameters $\alpha_{\nu}$ and $\zeta_{\nu}$ are appropriately chosen so as to guarantee global convergence.

At each step we perform backtracking line search on $\alpha_{\nu}$ until we either trigger fast convergence (K1 steps) or ensure global convergence (K2 steps) as shown in Algorithm 1. The K2 step, cf. (2hb), can be interpreted as a projection of the current iterate on a hyperplane generated by $w^{\nu}$, separating the set of fixed points from $u^{\nu}$, and thus guarantees that every iterate comes closer to fixed point set. Alongside, a sufficient decrease of the norm of the residual, $\|Ru^{\nu}\|$, may trigger a “blind update” (K0 steps) of the form $u^{\nu+1}{}={}u^{\nu}{}+{}d^{\nu},$ where no line search iterations need to be executed.

By exploiting the structure of $T$ and, in particular, linearity of $(I+Q)^{-1}$, we may avoid evaluating linear system solves at every backtracking step. Instead, we only need to evaluate $\Pi_{\mathcal{C}}$ once in every backtracking iteration. In particular, for $w^{\nu}=u^{\nu}+\alpha d^{\nu}$, we have

[TABLE]

where $\tilde{u}^{\nu}$ has already been computed, since it is needed in the evaluation of $Ru^{\nu}$, while $\tilde{d}^{\nu}$ solves $(I+Q)\tilde{d}^{\nu}=d^{\nu}$. The computation of $\tilde{d}^{\nu}$, which is the most costly operation, is performed only once, before the backtracking procedure takes place. The computation of the fixed-point residual of $w$ is also easily computed by $Rw^{\nu}{}={}\tilde{w}^{\nu}-\operatorname{\Pi}_{\mathcal{C}}(2\tilde{w}^{\nu}-w^{\nu})$.

Overall, save the computation of the residuals, at every iteration of SuperSCS we need to solve the linear system (2da) twice and invoke $\Pi_{\mathcal{C}}$ exactly $1+l_{\nu}$ times, where $l_{\nu}$ is the number of backtracks.

IV-C Termination

[TABLE]

The relative infeasibility certificate is defined as (note that $\bar{y}^{\nu}\in\mathcal{K}^{*}$ as a result of the projection step)

[TABLE]

Likewise, the relative unboundedness certificate is defined as

[TABLE]

Provided that $\bar{u}^{\nu}$ is not a feasible $\epsilon$-optimal point, it is a certificate of unboundedness if $\mathrm{uc}_{\nu}<\epsilon$ and it is a certificate of infeasibility if $\mathrm{ic}_{\nu}<\epsilon$.

IV-D Quasi-Newtonian directions

Quasi-Newtonian directions $d^{\nu}$ can be computed according to the general rule

[TABLE]

where invertible linear operators $H_{\nu}$ are updated according to certain low-rank updates so as to satisfy certain secant conditions starting from an initial operator $H_{0}$.

IV-D1 Restarted Broyden directions

Here we make use of Powell’s trick to update linear operators $B_{\nu}$ in such a way so as to enforce nonsignularity using the recursive formula [20, 24]:

[TABLE]

where $z^{\nu}{}={}w^{\nu}{}-{}u^{\nu}$, $\xi^{\nu}{}={}Rw^{\nu}{}-{}Ru^{\nu}$ and for a fixed parameter $\bar{\vartheta}\in(0,1)$ and $\gamma_{\nu}{}={}\<H_{\nu}\xi^{\nu},z^{\nu}\>/\|z^{\nu}\|^{2}$ we have $\tilde{\xi}^{\nu}{}={}(1-\theta_{\nu})B_{\nu}z^{\nu}+\theta_{\nu}\xi^{\nu}$ with

[TABLE]

and the convention $\operatorname{sgn}(0)=1$. Using the Sherman-Morisson formula, operators $H_{\nu}$ are updated as follows:

[TABLE]

thus lifting the need to compute and store matrices $B_{\nu}$.

The Broyden method requires that we store matrices of dimension $(m{}+n{}+1){}\times{}(m{}+n{}+1)$. Here we employ a limited-memory restarted Broyden (RB) method which affords us a computationally favorable implementation using buffers of length $\mathrm{mem}$, that is $Z_{\nu}{}={}[z^{\nu}~{}z^{\nu-1}~{}\cdots~{}z^{\nu-\mathrm{mem}+1}],$ and $\tilde{Z}_{\nu}{}={}[\tilde{z}^{\nu}~{}\tilde{z}^{\nu-1}~{}\cdots~{}\tilde{z}^{\nu-\mathrm{mem}+1}],$ where $\tilde{z}^{i}$ are the auxiliary variables $\tilde{z}^{i}{}\dfn{}\nicefrac{{z^{i}{}-{}H_{i}\tilde{\xi}^{i}}}{{\<s^{i},H_{i}\tilde{\xi}^{i}\>}}$.

We have observed that SuperSCS performs better when deactivating K0 steps or using a small value for $c_{0}$ (e.g., $c_{0}=0.1$) when using restarted Broyden directions.

IV-D2 Anderson’s acceleration

Anderson’s acceleration (AA) imposes a multi-secant condition [26, 9]. In particular, at every iteration $\nu$ we update a buffer of $\mathrm{mem}$ past values of $z$ and $\xi$, that is we construct a buffer $Z_{\nu}$ as above and a buffer $\Xi_{\nu}{}={}\begin{bmatrix}\xi^{\nu}&\xi^{\nu-1}&\cdots&\xi^{\nu-{\mathrm{mem}}+1}\end{bmatrix}$. Directions are computed as

[TABLE]

where $t^{\nu}$ is a least-squares solution of the linear system $\Xi_{\nu}t^{\nu}=Ru^{\nu},$ that is $t^{\nu}$ solves

[TABLE]

and can be solved using the singular value decomposition of $\Xi_{\nu}$, or a QR factorization which may be updated at every iteration [26]. In practice Anderson’s acceleration works well for short memory lengths, typically between 3 and 10, and, more often than not, outperforms the above restarted Broyden directions.

IV-E Convergence

The convergence properties of SuperSCS are inherited by those of the general SuperMann scheme [24]. In particular, under a weak boundedness assumption for the quasi-Newtonian directions $d^{\nu}$, $Ru^{\nu}$ converges to zero and $u^{\nu}$ converges to a $u^{\star}$ which satisfies the HSDE (2c). If, additionally, $R$ is metrically subregular at $u^{\star}$ — a weak assumption — then, the convergence is R-linear. Under additional assumptions, SuperSCS with the full-memory counterpart of the above restarted Broyden scheme can be shown to converge superlinearly. The restarted Broyden directions of Algorithm 2 and Anderson’s acceleration, though not proven superlinear directions, exhibit steep linear convergence as shown in the next section (see Fig. 2) and have low memory requirements.

V Benchmarks and results

V-A Benchmarking methodology

In order to compare different solvers in a statistically meaningful way, we use the Dolan-Moré (DM) plot [7] and the shifted geometric mean [16]. The DM plot allows us to compare solvers in terms of their relative performance (e.g., computation time, flops, etc) and robustness, i.e., their ability to successful solve a given problem up to a certain tolerance.

Let $P$ be a finite set of test problems and $S$ a finite set of solvers we want to compare to one another. Let $t_{p,s}$ denote the computation time that solver $s$ needs to solve problem $p$. We define the ratio between $t_{p,s}$ and the lowest observed cost to solve this problem using a solver from $S$ as

[TABLE]

If $s$ cannot solve $p$ at all, we define $t_{p,s}=+\infty$ and $r_{p,s}=+\infty$. The cumulative distribution of the performance ratio is the DM performance profile plot. In particular, define

[TABLE]

for $\tau{}\geq{}1$. The DM performance profile plot is the plot of $\rho_{s}(\tau)$ versus $\tau$ on a logarithmic $x$-axis. For every $s\in S$, the value $\rho_{s}(1)$ is the probability of solver $s$ to solve a given problem faster than all other solvers, while $\lim_{\tau\to\infty}\rho_{s}(\tau)$ is the probability that solver $s$ solves a given problem at all.

As demonstrated in [10], DM plots aim at comparing multiple solvers to the best one (cf. (2p)), than to one another. Therefore, alongside we shall report the shifted geometric mean of computation times for each solver following [16]. For solver $s{}\in{}S$, define the vector $t^{s}{}\in{}\R^{|P|}$ with

[TABLE]

The shifted geometric mean of $t^{s}$ with shifting parameter $\sigma{}\geq{}0$ is defined as

[TABLE]

Hereafter we use $\sigma=10\,\mathrm{s}$.

In what follows we compare SuperSCS with the quasi-Newtonian direction methods presented in Section IV-D against SCS [19, 18]. All tolerances are fixed to $10^{-4}$. In order to allow for a fair comparison among algorithms with different per-iteration cost, we do not impose a maximum number of iterations; instead, we consider that an algorithm has failed to produce a solution — or a certificate of unboundedness/infeasibility — if it has not terminated after a certain (large) maximum time. All benchmarks were executed on a system with a quad-core i5-6200U CPU at $2.30\,\mathrm{GHz}$ and $12\,\mathrm{GB}$ RAM running Ubuntu 14.04.

V-B Semidefinite programming problems

Let us consider the problem of sparse principal component analysis with an $\ell_{1}$-regularizer which has the form [5].

[TABLE]

A total of 288 randomly generated problems was used

for benchmarking with $d\in\{50,120,$ $140,180\}$ and $\lambda\in$ $\{0.1,2,5\}$, where $d$ is the dimension of $Z$. The DP plot in Fig. 1 shows that SuperSCS is consistently faster and more robust compared to SCS. In Table I we see that SuperSCS is faster than SCS by more than an order of magnitude; in particular, SuperSCS with AA directions and memory $5$ was found to perform best.

Additionally, for this benchmark we tried the interior-point solvers of SDPT3 [25] and Sedumi [23]. Both exhibited similar performance; in particular, for problems of dimension $Z\in\R^{140\times 140}$, SDPT3 and Sedumi required $6500\,\mathrm{s}$ to $7500\,\mathrm{s}$ and for problems of dimension $Z\in\R^{180\times 180}$ they required $19000\,\mathrm{s}$ to $21500\,\mathrm{s}$. Note that SuperSCS (with RB and memory 50) solves all problems in no more than $11.7\,\mathrm{s}$. Additionally, interior-point methods have an immense memory footprint of several GB, in this example whereas SuperSCS with RB and memory 100 needs as little as $211.1\,\mathrm{MB}$ and with AA and memory 5 consumes just $46.2\,\mathrm{MB}$.

V-C LASSO problems

Regularized least-squares problems with the $\|{}\cdot{}\|_{1}$-regularizer, also known as LASSO problems, are optimization problems of the form

[TABLE]

where $A\in\R^{m\times n}$ is a (sparse) matrix, and $\mu>0$ is the regularization weight. LASSO problems find applications in statistics and compressed sensing [22]. LASSO problems are cast as second-order cone programs [13].

LASSO problems aim at finding a sparse vector $x$ which minimizes $\|Ax-b\|^{2}$. The sparseness of the minimizer $x^{\star}$ can be controlled by $\mu\in\{0.01,0.1,1\}$. We tested 1152 randomly generated LASSO problems with $n\in\{631,1000,1585,2512\}$, $m=\lceil n/5\rceil$, and matrices $A$ with condition numbers $\kappa_{A}\in\{10,215,4600,10^{5}\}$. The DM plot in Fig. 3 and the statistics presented in Table II demonstrate that SuperSCS, both with RB and AA directions, outperforms SCS.

V-D Sparse $\ell_{1}$-regularized logistic regression

Logistic regression is a regression model where dependent variables are binary [4]. The $\ell_{1}$ regularized variant of logistic regression aims at performing simultaneous regression and feature selection and amounts to solving an optimization problem of the following form [1]:

[TABLE]

Similar to LASSO (Sec. V-C), parameter $\lambda>0$ controls the sparseness of the solution $w^{\star}$. Such problems can be cast as conic programs with the exponential cone, which is not self-dual. A total of 288 randomly generated problems was used for benchmarking

with $p\in\{80,100\}$, $q\in\{50,100,120\}$ and $\lambda\in\{10,20,50\}$. In Fig. 4 we observe that SuperSCS with RB directions is slower compared to SCS, however SuperSCS with AA is noticeably faster.

V-E Maros-Mészáros QP problems

Here we present performance of SCS and SuperSCS (with RB and AA directions) on this collection of problems on the Maros-Mészáros collection of problems [15]. As shown in Fig. 5 SuperSCS, both with RB directions and Anderson’s acceleration, is faster and more robust compared to SCS. The two quasi-Newtonian directions appear to be on a par.

VI Conclusions

In this work we introduced SuperSCS: a first-order method for large-scale conic optimization problems which combines the low iteration cost of SCS and the fast convergence of SuperMann. We have compared SuperSCS with SCS on a broad collection of conic optimization problems of practical interest. Using Dolan-Moré plots and runtime statistics, we demonstrated that SuperSCS with Anderson’s acceleration is faster and more robust than SCS.

The C implementation of SuperSCS builds up on SCS and is a free open-source software. SuperSCS can be interfaced from MATLAB, Python, can be invoked via CVX, CVXPy and YALMIP and is also available as a Docker image (see https://kul-forbes.github.io/scs/).