SuperSCS
1.3.2

SCS and SuperSCS solve the following problem
\begin{eqnarray*} &&\mathrm{Minimize}\ \langle c, x \rangle\\ &&Ax + s = b\\ &&s\in\mathcal{K}, \end{eqnarray*}
where \(x\in\mathbb{R}^n\), \(s\in\mathbb{R}^m\) are the primal variables and \(\mathcal{K}\subseteq\mathbb{R}^m\) is a nonempty, closed, convex cone.
The matrix \(A\in\mathbb{R}^{m\times n}\) and the vector \(b\in\mathbb{R}^m\) are the problem data.
The algorithm makes use of the homogeneous selfdual embedding (HSDE) which is the problem of finding a \(u=(x,y,\tau)\in\mathbb{R}^{n+m+1}\) so that
\begin{eqnarray*} &&u\in\mathcal{C}\\ &&Qu\in\mathcal{C}^*\\ &&\langle u, Qu \rangle = 0, \end{eqnarray*}
where
\begin{eqnarray*} Q = \begin{bmatrix} 0 & A^* & c\\ A & 0 & b\\ c^* & b^* & 0 \end{bmatrix} \end{eqnarray*}
and \(\mathcal{C} = \mathbb{R}^n\times \mathcal{K}^* \times \mathbb{R}_+\) is a cone.
Equivalently, the HSDE can be written as a variational inequality:
\begin{eqnarray*} 0 \in Qu + N_{\mathcal{C}}(u), \end{eqnarray*}
where \(N_{\mathcal{C}}\) is the normal cone of \(\mathcal{C}\).
The DouglasRachford splitting can be applied to the above variational problem.
Operator \(T\) is a (firmly) nonexpansive operator on which we may apply the SuperMann algorithmic scheme.
Since \(Q\) is a skewsymmetric linear operator, it is maximally monotone.
Since \(N_{\mathcal{C}}\) is the subdifferential of the proper convex function \(\delta_{\mathcal{C}}\), it is maximally monotone.
Additionally, because of Cor. 24.4(i) in [BauCom], \(Q+N_{\mathcal{C}}\) is maximally monotone.
Therefore, we may apply the DouglasRachford splitting on this monotone inclusion.
The SCS algorithm [ODon16] is precisely the application of the DouglasRachford splitting (DRS) to \(N_{\mathcal{C}}{}+{}Q\); this leads to the following iterations; see Sec. 7.3 in [RyuBoy16]
\begin{eqnarray*} \tilde{u}^{\nu} {}={}& (I+Q)^{1}(u^{\nu}) \\ \bar{u}^{\nu} {}={}& \Pi_{\mathcal{C}}(2\tilde{u}^{\nu}  u^{\nu}) \\ u^{\nu+1} {}={}& u^{\nu} {}+{} \bar{u}^{\nu} {}{} \tilde{u}^{\nu} \end{eqnarray*}
Here, on the other hand, we consider the opposite DRS splitting, \(Q+N_{\mathcal{C}}\), which leads to the following DRS iterations
\begin{eqnarray*} \tilde{u}^{\nu} {}={}& \Pi_{\mathcal{C}}(u^{\nu}) \\ \bar{u}^{\nu} {}={}& (I+Q)^{1}(2\tilde{u}^{\nu}  u^{\nu}) \\ u^{\nu+1} {}={}& u^{\nu} {}+{} \bar{u}^{\nu} {}{} \tilde{u}^{\nu} \end{eqnarray*}
For any initial guess \(u^0\), the iterates \(u^\nu\) converge to a point \(u^\star\) which satisfies the above monotone inclusion; see Thm. 25.6(i), (iv) in [BauCom11].
The linear system above can be either solved "directly" using a sparse LDL factorization or "indirectly" by means of the conjugate gradient method [ODon16].
The projection on \(\mathcal{C}\) essentially requires that we be able to project on the dual cone, \(\mathcal{K}^*\)
The iterative method can be concisely written as
\begin{eqnarray*} u^{\nu+1} = Tu^\nu, \end{eqnarray*}
where \(T:\mathbb{R}^{N}\to\mathbb{R}^{N}\) is a (firmly) nonexpansive operator [RyuBoy16]. As such it fits the Krasnosel'skiiMann framework, because of Sec. 5.2 in [BauCom11], leading to the relaxed iterations
\begin{eqnarray*} u^{\nu+1} {}={} (1\lambda)u^\nu + \lambda Tu^{\nu}, \end{eqnarray*}
with \(\lambda \in (0,2)\) since \(T\) is firmly nonexpansive and, as a result, it fits the SuperMann framework [ThePat16].
SuperMann [ThePat16] reduces the problem of finding a fixedpoint \(x^\star{}\in{}\operatorname{fix} T\) to that of finding a zero of the residual operator
\begin{eqnarray*} R = IT. \end{eqnarray*}
SuperMann, instead of applying Krasnosel'skiiManntype updates as discussed above, takes extragradienttype updates of the general form
\begin{eqnarray*} w^{\nu} {}={}& u^{\nu} {}+{} \tau_\nu d^\nu, \\ u^{\nu+1} {}={}& u^\nu  \zeta_{\nu} Rw^\nu, \end{eqnarray*}
where \( d^{\nu} \) are fast, e.g., quasiNewtonian, directions and scalar parameters \(\tau_{\nu}\) and \(\zeta_{\nu}\) are appropriately chosen so as to guarantee global convergence.
At each step either trigger fast convergence (K1 steps) or ensure global convergence (K2 steps) as shown below.
Alongside, a sufficient decrease of the norm of the residual, \(\Ru^{\nu}\\), may trigger a "blind update" of the form \( u^{\nu+1} {}={} u^{\nu} {}+{} d^{\nu}. \)
SuperSCS Algorithm
By exploiting the structure of \(T\) and, in particular, the fact that \((I+Q)^{1}\) is a linear operator, we may avoid applying \(T\) at every iteration of the line search loop.
Instead, at every line search iteration we only need to apply \(\Pi_{\mathcal{C}}\) once.
SuperSCS does not need to solve a linear system within the line search loop.
Overall, save the computation of the residuals, at every iteration of SuperSCS we need to solve the linear system twice and invoke \(\Pi_{\mathcal{C}}\) only \(1+l_\nu\) times, where \(l_{\nu}\) is the number of line search iterations.