# GeNIOS: an (almost) second-order operator-splitting solver for large-scale
convex optimization

Theo Diamandis
<EMAIL_ADDRESS>Zachary Frangella
<EMAIL_ADDRESS>Shipu Zhao
<EMAIL_ADDRESS>Bartolomeo Stellato
<EMAIL_ADDRESS>Madeleine Udell
<EMAIL_ADDRESS>
(October 2023)

###### Abstract

We introduce the GEneralized Newton Inexact Operator Splitting solver (GeNIOS)
for large-scale convex optimization. GeNIOS speeds up ADMM by approximately
solving approximate subproblems: it uses a second-order approximation to the
most challenging ADMM subproblem and solves it inexactly with a fast
randomized solver. Despite these approximations, GeNIOS retains the
convergence rate of classic ADMM and can detect primal and dual infeasibility
from the algorithm iterates. At each iteration, the algorithm solves a
positive-definite linear system that arises from a second-order approximation
of the first subproblem and computes an approximate proximal operator. GeNIOS
solves the linear system using an indirect solver with a randomized
preconditioner, making it particularly useful for large-scale problems with
dense data. Our high-performance open-source implementation in Julia allows
users to specify convex optimization problems directly (with or without conic
reformulation) and allows extensive customization. We illustrate GeNIOS’s
performance on a variety of problem types. Notably, GeNIOS is more than ten
times faster than existing solvers on large-scale, dense problems.

##  1 Introduction

Data sets in modern optimization problems are large, motivating the search for
optimization algorithms that scale well with the problem size. The alternating
direction method of multipliers (ADMM) is a particularly powerful algorithm
for tackling these large, data-driven optimization problems. Compared to
interior point methods, ADMM has a modest per-iteration cost and is easy to
parallelize. While the method is slow to produce a high-accuracy solution,
ADMM often finds a low-accuracy solution quickly, which usually suffices for
problems with real-world—and often noisy—data. However, as problem sizes
increase, even ADMM iterations can become unacceptably slow. In particular,
similar to [SBL20], we observe slowdowns for problems with data matrices with
tens of millions of non-zeros or more—a scale easily exceeded by problems with
dense datasets. For example, a dense LASSO problem on, e.g., a gene expression
dataset with $n=1,000$ samples and $p=10,000$ variables is already this large.

In this paper, we introduce the GEneralized Newton Inexact Operator Splitting
solver (GeNIOS, pronounced “genie-ōs”), a new inexact ADMM solver for convex
optimization problems. GeNIOS speeds up ADMM by approximating the ADMM
subproblems at each iteration, and then solving these approximate subproblems
inexactly. Each iteration has two steps: 1) solve a linear system that results
from a second-order approximation to the ADMM subproblem inexactly with a fast
randomized solver, the Nyström preconditioner [FTU21]; 2) compute an
approximate proximal operator. We observe speedups of up to 50$\times$
compared to classic ADMM on large-scale dense optimization problems. GeNIOS
maintains classic ADMM convergence guarantees by controlling subproblem
errors, as outlined in recent theoretical work [Fra+23]. GeNIOS is the first
general-purpose ADMM solver to exploit both forms of inexactness: subproblem
approximation and inexact solves.

###  1.1 The optimization problem

Consider the optimization problem

minimize $\displaystyle f(x)+g(z)$ (1) subject to $\displaystyle Mx-z=c,$

where $x\in\mathbf{R}^{n}$ and $z\in\mathbf{R}^{m}$ are decision variables,
and $f:\mathbf{R}^{n}\to\mathbf{R}$, $g:\mathbf{R}^{m}\to\mathbf{R}$,
$M\in\mathbf{R}^{m\times n}$, and $c\in\mathbf{R}^{m}$ are the problem data.
Assume that the function $f$ is smooth and convex, and that the function $g$
is convex, proper, and lower-semicontinuous. Thus, problem (1) is a convex
optimization problem, which we will refer to as a _convex program_. The
flexibility of this formulation can provide an important speedup for
statistical learning problems relative to conic reformulation. The formulation
can also be specialized to recover many special cases, including quadratic
programs and conic programs.

#### Quadratic programs.

A _quadratic program_ (QP) has the form

minimize $\displaystyle(1/2)x^{T}Px+q^{T}x$ subject to $\displaystyle l\leq
Mx\leq u,$

where $x\in\mathbf{R}^{n}$ is the decision variable and
$P\in\mathbf{S}_{+}^{n}$, $q\in\mathbf{R}^{n}$, $M\in\mathbf{R}^{m\times n}$,
$l\in\mathbf{R}^{m}$, and $u\in\mathbf{R}^{m}$ are the problem data. Linear
equality constraints are encoded by setting $l_{i}=u_{i}$ for some
$i\in\\{1,\dots,m\\}$. A QP is a special case of (1) with

$f(x)=(1/2)x^{T}Px+q^{T}x\qquad\text{and}\qquad
g(z)=I_{[l,u]}(z)=\begin{cases}0&l\leq z\leq u\\\
+\infty&\text{otherwise,}\end{cases}$

and $c=0$. We call $I_{S}(z)$ the _indicator function_ of the set $S$. Linear
programs can be written in this form by setting $P=0$. QPs are ubiquitous in
practice and encompass a wide variety of problems in disparate fields:
portfolio optimization in finance [Mar52, Boy+17]; model predictive control
[GPM89, Raw00]; denoising in signal processing [PE10]; model fitting in
machine learning [SNW12]; and the transport problem in operations research
[Kan48, Dan51], among others.

#### Conic programs.

Let the function $g$ instead be the indicator function of a convex cone to
recover the _conic program_ :

minimize $\displaystyle(1/2)x^{T}Px+q^{T}x$ (2) subject to $\displaystyle
Mx-z=c$ $\displaystyle z\in\cal K,$

where $\mathcal{K}$ is a non-empty, closed, convex cone. Again, we recognize
(2) as a special case of (1) with

$f(x)=(1/2)x^{T}Px+q^{T}x\qquad\text{and}\qquad g(z)=I_{\mathcal{K}}(z).$

Any convex optimization problem can be written in this form [NN92], including
linear programs, quadratic programs, second-order cone programs (SOCPs), and
semidefinite programs (SDPs), and many modeling languages can translate convex
optimization problems to conic form [GB14, DB16, Agr+18, Ude+14]. As a result,
this formulation is used by many of the popular convex optimization solvers
(e.g., SCS [O’D+16, O’D21], COSMO [GCG21], Hypatia [CKV22], and Mosek
[ApS22]). SOCPs appear in robust optimization [BTN98, BTN99, BTEGN09], model
predictive control [BAS10], and many engineering design problems [Lob+98].
Applications of SDPs include convex relaxations of binary optimization
problems [LS91], experiment design [BV04, §7.5], circuit design [VBEG97,
VBEG98], and sum-of-squares programs [BPT12, PL03, Lau09].

#### Machine learning problems.

Consider machine learning problems of the form

minimize
$\displaystyle\sum_{i=1}^{N}\ell(a_{i}^{T}x-b_{i})+(1/2)\lambda_{2}\|x\|_{2}^{2}+\lambda_{1}\|x\|_{1},$

where $\ell:\mathbf{R}\to\mathbf{R}_{+}$ is a convex per-sample loss function,
$\\{(a_{i},b_{i})\\}_{i=1}^{N}\subseteq\mathbf{R}^{n}\times\mathbf{R}$ are
problem data, and $\lambda_{1}$, $\lambda_{2}\geq 0$ are regularization
parameters. In the framework of (1), take

$f(x)=\sum_{i=1}^{N}\ell(a_{i}^{T}x-b_{i})+(1/2)\lambda_{2}\|x\|_{2}^{2}\qquad\text{and}\qquad
g(z)=\lambda_{1}\|z\|_{1},$

and $M=I$, $c=0$. The conic reformulation of this problem typically has many
additional variables. As an example, consider the $\ell_{1}$-regularized
logistic regression problem:

minimize
$\displaystyle\sum_{i=1}^{N}\log\left(1+\exp(b_{i}a_{i}^{T}x)\right)+\lambda_{1}\|z\|_{1}$
subject to $\displaystyle x-z=0.$

The reformulation of logistic regression as a conic program (specifically, an
exponential cone program) has at least $2n+3N$ variables (see appendix B). In
general, this reformulation slows solve time since the per-iteration time
scales superlinearly in the problem size. For example, if the dominant
operation is an $O(n^{2})$ matrix-vector product, increasing the number of
variables by a factor of $5$ can increase solve time by a factor of $25$.
Handling machine learning problems directly via GeNIOS provides a significant
performance boost over general-purpose conic solvers.

###  1.2 Solution methods

A variety of algorithms can be used to solve (1) and its specific subclasses.
We briefly review some of the most popular and relevant methods.

#### Interior point methods and solvers.

Primal-dual interior point methods (IPMs), with roots largely in work by
Karmarkar, Nesterov, Nemirovisky, and Mehrotra [Kar84, NN94, Meh92], have
historically been the method of choice for convex optimization. Commercial
convex optimization solvers, including Gurobi [Gur23], Mosek [ApS22], and
CVXGEN [MB12] use IPMs. Several open-source IPM solvers exist as well,
including Clarabel [CG23], ECOS [DCB13], CVXOPT [ADV+13], and Hypatia [CKV22].
Most similar to the spirit of our work, Hypatia goes beyond the set of
_standard cones_ (usually the positive orthant, second order cone,
semidefinite cone, exponential cone, and power cone) used by other solvers and
implements many _exotic cones_ , which allow the user to more naturally
specify the problem and avoid problem size bloat under reformulation. Still,
this solver uses a conic formulation (2), just with an expanded set of cones
$\mathcal{K}$. IPMs quickly produce very accurate solutions for small and
medium-sized problems. However, they cannot be easily warm-started, and they
do not scale well for very-large problems due to the formation and
factorization of a large matrix at each iteration. A natural alternative to
address the scaling issue, is to consider indirect methods for solving the
linear system, as they only require matrix-vector products. Unfortunately, the
condition number of the associated linear system becomes large as the solver
nears a solution, rendering indirect methods for solving the system
ineffective [Gon12].

#### First-order methods.

The poor scaling of interior point methods and rise of large-scale, data-
driven optimization prompted the resurgence of first-order methods. These
methods can be implemented so that the iteration time is dominated by matrix-
vector product computations, which can be accelerated on modern computing
architectures. Although first-order methods are slow to converge to high
accuracy solutions, they usually produce moderately accurate solutions
quickly. Two of the most common methods in solvers today are ADMM [GM75, GM76,
DR56, LM79, Gab83] (see [Boy+11] for a modern survey) and primal-dual hybrid
gradient (PDHG), also called Chambolle-Pock [CP11]. Both of these methods (and
many others) are special cases of Rockafellar’s proximal point algorithm
[Roc76, EB92] (see [RY22, §3] and [ZLU21] for related derivations).

#### First-order solvers.

In the last decade, many open-source solvers have implemented variants of ADMM
or PDHG. The conic solver SCS [O’D+16, O’D21] applies ADMM directly to the
homogenous self-dual embedding [YTM94, XHY96] of a conic program and was the
first ADMM-based solver that could handle infeasible or unbounded problems.
Later, the QP solver OSQP [Ste+20] and conic solver COSMO [GCG21] used the
results of [Ban+19] to apply ADMM more directly to (2) while still being able
to detect infeasibility and unboundedness by looking at differences of the
iterates. PDHG has been applied to solve large-scale LPs in PDLP [App+21,
App+22] and SDPs in ProxSDP [SGV22]. Notably, PDHG does not require a linear
system solve at each iteration for conic programs of the form (2). Both PDLP
and ProxSDP use the method of [App+21a] to detect infeasibility, but to the
best of our knowledge, it’s unclear to what extent these theoretical results
extend to the case of conic programs with inexact subproblem solves. The ADMM-
based solvers require the solution of a linear system at each iteration, and
although this system can be solved by indirect methods, ill-conditioning of
the problem data—a common phenomenon in real-world data matrices [UT19]—can
slow convergence of these methods. In addition, all of these solvers require
problems to be passed in a conic form resembling (2), which typically leads to
a substantial increase in problem size, especially in machine learning
problems.

#### Beyond conic programs.

All of the solvers discussed solve conic problems of a form similar to (2), or
solve a specialized version of this form more tailored to a specific problem
class (e.g., LPs or QPs). While this form allows solvers to handle essentially
all convex optimization problems provided by a user (usually with the help of
a modeling framework such as JuMP [DHL17] or Convex.jl [Ude+14]), transforming
these problems to conic form can make the problem more difficult to solve by
increasing the size or obscuring the structure of the constraint matrix $M$ or
objective matrix $P$. Another line of work avoids transforming problems into
conic form by building solvers directly on a library of proximal operators,
including Epsilon [WWK15], which is unmaintained, POGS [FB18], which requires
a separable objective, and ProxImaL [Hei+16], which focuses on image
optimization problems. In general, defining $M$ directly as a linear operator
(rather than as a concrete matrix) speeds up computation but can be difficult
to achieve for solvers that require a conic form. SCS does support matrix-free
linear operators but still requires a conic form problem [DB16a].

#### Inexact ADMM.

It is well known that ADMM applied to (1) converges to an optimum at an
$O(1/k)$ rate if both subproblems are solved exactly [HY12, MS13]. Classic
work [EB92] established that ADMM also converges when subproblems are solved
inexactly, provided the errors of subproblems are summable. A more
contemporary line of work [Ouy+15, DY16] considers replacing the $x$\- and
$z$-subproblems of ADMM with approximate subproblems that are easier to solve.
These papers establish convergence of the resulting methods, provided the
approximate subproblems are solved exactly. In recent work [Fra+23], the
authors combined these two approaches—approximate subproblems and inexact
solves of these subproblems—and show that the combination still achieves
ADMM’s $O(1/k)$ convergence rate. This last result motivates GeNIOS, as these
subproblem approximations and inexact solves significantly reduce ADMM
iteration time. (The approximation of the $x$-subproblem also resembles
sequential quadratic programming [Wil63], which is frequently used to
heuristically solve nonconvex problems [BT95].) Recent work [ZFU22] showed
that combining randomized Nyström preconditioning and approximate
$x$-subproblem solves with ADMM yields significant speed-ups over standard
solvers for a variety of unconstrained machine learning problems. Similar
second-order approximations in the subproblems are used in the ALADIN solver
[HFD16]. To the best of the authors’ knowledge, there are no open-source
solvers for general convex programs that implement ADMM with inexact solves of
approximate subproblems.

| GeNIOS | COSMO | OSQP | SCS | ProxSDP
---|---|---|---|---|---
method | ADMM | ADMM | ADMM | ADMM | PDHG
interface | any convex | conic | QP | conic | SOCP, SDP
approximate subproblems | yes | no | no | no | no
linear system solver | indirect | (in)direct | (in)direct | (in)direct | n/a
inexact solves | yes | yes | GPU only | yes | n/a
preconditioning | Nyström | diagonal | diagonal | diagonal | n/a
inexact projection | yes | no | n/a | no | yes
Table 1:  GeNIOS offers a more flexible interface and more inexactness
compared to the most-similar existing solvers. ProxSDP uses an inexact
projection onto the positive semidefinite cone. A conic interface indicates
that the solver solves (2) for the full set of standard cones: the positive
orthant, second-order cone, positive semidefinite cone, exponential cone, and
power cone.

###  1.3 Contributions

This paper showcases, through the GeNIOS solver, how several recent
theoretical ideas can be combined to speed up ADMM while preserving
convergence guarantees. Our contributions can be summarized as follows:

  * •

We develop the GeNIOS.jl solver: the first open-source solver that

    * –

accelerates subproblem solves by inexactly solving approximate ADMM
subproblems using Nyström preconditioning; and

    * –

allows direct specification of ADMM problems via function, gradient and
Hessian oracles for $f$ and function and proximal operator oracles for $g$.

  * •

We show that GeNIOS can dectect infeasibility in conic programs despite
inexact solves using [RGN22], and we show that GeNIOS converges at the
standard $O(1/k)$ rate based on results from the authors’ prior work [Fra+23].

  * •

We show how GeNIOS allows the user to exploit problem structure by specifying
ADMM problems directly (not in conic form) and by using the Julia programming
language’s multiple dispatch to implement efficient linear operators.

  * •

We showcase GeNIOS’s up to 50$\times$ speedup over classic ADMM on a variety
of optimization problems with real-world and simulated data.

The code is available online at

https://github.com/tjdiamandis/GeNIOS.jl

with documentation that includes several examples. Table 1 compares GeNIOS to
the most similar existing solvers.

#### Roadmap.

We overview the GeNIOS algorithm in §2, including the overall method,
convergence guarantees, infeasibility detection, and randomized
preconditioning for the linear system solve at each iteration. In §3, we
discuss the solver interface for general convex optimization problems and the
specialized interfaces for quadratic programming and machine learning
problems. In §4, we numerically demonstrate the performance improvements
gained by leveraging inexactness, randomized preconditioning, and the more
natural problem formulations allowed by GeNIOS. Finally, we point to some
directions of future work in §5.

##  2 Method

Our method GeNIOS uses inexact ADMM and techniques from randomized numerical
linear algebra to speed up solve times on large-scale optimization problems.
This method builds on ideas in GeNI-ADMM from the authors’ previous work
[Fra+23], and the implementation is largely inspired by OSQP [Ste+20]. Recall
that GeNIOS solves convex problems in the form (1),

minimize $\displaystyle f(x)+g(z)$ (1) subject to $\displaystyle Mx-z=c,$

where $x\in\mathbf{R}^{n}$ and $z\in\mathbf{R}^{m}$ are decision variables,
and $f:\mathbf{R}^{n}\to\mathbf{R}$, $g:\mathbf{R}^{m}\to\mathbf{R}$,
$M\in\mathbf{R}^{m\times n}$ and $c\in\mathbf{R}^{m}$ are the problem data.
The function $f$ is a smooth and convex, and the function $g$ is a convex,
proper, and lower-semicontinuous. GeNIOS requires the ability to evaluate $f$
and $g$, the gradient $\nabla f$, a Hessian-vector product (HVP)
$v\mapsto\nabla^{2}f(x)v$, and the proximal operator of $g$,

$\mathbf{prox}_{g/\rho}(v)=\mathop{\rm
argmin}_{\tilde{z}}g(\tilde{z})+(\rho/2)\|\tilde{z}-v\|_{2}^{2}.$

The HVP and proximal operator may be approximate, subject to a condition on
the incurred errors. This flexibility allows GeNIOS to easily model a variety
of problems of interest, as gradients and HVPs can be easily
specified—including with automatic differentiation—and proximal operators can
be efficiently computed for many functions (see, e.g., [CP11a, PB+14, WWK15]
and references therein). As a result, GeNIOS not only handles conic programs
but also specializes to problems with bespoke objective functions, such as
robust regression problems in machine learning.

GeNIOS replaces the exact subproblem solutions of classic ADMM with inexact
ones. The standard (scaled) ADMM algorithm applied to (1) consists of the
iterations

$\displaystyle x^{k+1}$ $\displaystyle=\mathop{\rm
argmin}_{x}\left(f(x)+(\rho/2)\|Mx-z^{k}-c+u^{k}\|_{2}^{2}\right)$ (3)
$\displaystyle z^{k+1}$ $\displaystyle=\mathop{\rm
argmin}_{z}\left(g(z)+(\rho/2)\|Mx^{k+1}-z-c+u^{k}\|_{2}^{2}\right)$ (4)
$\displaystyle u^{k+1}$ $\displaystyle=u^{k}+Mx^{k+1}-z^{k+1}-c.$ (5)

GeNIOS replaces the function $f$ in the $x$-subproblem with a second-order
approximation,

$f(x)\approx f(x^{k})+\nabla
f(x^{k})^{T}(x-x^{k})+(\eta^{k}/2)\|x-x^{k}\|^{2}_{\Theta_{k}},$

where $\Theta_{k}$ is a positive semidefinite matrix and $\eta^{k}$ is a
positive parameter. GeNIOS sets $\Theta^{k}=\nabla^{2}f(x^{k})+\sigma I$ with
$\sigma\geq 0$, and $\eta^{k}=1$, but these choices can be modified based on
computational considerations [Fra+23]. For example, the Hessian-vector-product
may be hard to compute exactly, in which case $\Theta^{k}$ can be some
approximation (e.g., the Hessian may only be updated every 10 iterations).

After the substitution of an approximate $f$, which may be viewed as a second-
order Taylor expansion, the $x$-subproblem becomes

$x^{k+1}=\mathop{\rm argmin}_{x}\Big{(}f(x^{k})+\nabla
f(x)^{T}(x-x^{k})+(1/2)\|x-x^{k}\|^{2}_{\nabla^{2}f(x^{k})+\sigma
I}+(\rho/2)\|Mx-z^{k}-c+u^{k}\|_{2}^{2}\Big{)}.$ (6)

Minimizing this unconstrained convex quadratic is equivalent to solving a
linear system. Thus, the approximation of $f$ simplifies and may significantly
speed up the solution to the $x$-subproblem. The $z$-subproblem (4) is
unchanged, except instead of requiring optimality, we only require a $z^{k+1}$
such that the $z$-subproblem objective value at this point is close to the
optimal value. In other words,

$\displaystyle g(z^{k+1})+(\rho/2)\|Mx^{k+1}+z^{k+1}-c+u^{k}\|_{2}^{2}$
$\displaystyle-$ $\displaystyle
g(z^{k+1,\star})+(\rho/2)\|Mx^{k+1}+z^{k+1,\star}-c+u^{k}\|_{2}^{2}<\varepsilon_{z}^{k},$

where

$z^{k+1,\star}=\mathbf{prox}_{g/\rho}(Mx^{k+1}-c+u^{k}).$

The subproblem errors must decay quickly to retain convergence guarantees
[Fra+23]: defining $x^{k,\star}$ as the exact solution to the original
$x$-subproblem, $p_{z}^{k,\star}$ as the optimal value of the $z$ subproblem,
and $p_{z}^{k}$ as the value when the $z$-subproblem objective is evaluated at
the approximate solution $z^{k+1}$, the method provably converges if the
errors

$\varepsilon_{x}^{k}=\|x^{k}-x^{k,\star}\|\qquad\text{and}\qquad\varepsilon_{z}^{k}=p_{z}^{k}-p_{z}^{k,\star}$
(7)

are summable and square root summable respectively, i.e.,

$\sum_{k=1}^{\infty}\varepsilon_{x}^{k}<\infty\qquad\text{and}\qquad\sum_{k=1}^{\infty}\sqrt{\varepsilon_{z}^{k}}<\infty.$
(8)

For detailed discussion of the theoretical assumptions, see [Fra+23, §4].

###  2.1 Solving the linear system

The $x$-subproblem update after approximation (6) is an unconstrained convex
QP. Thus, its solution is the solution to the linear system

$\left(\nabla^{2}f(x^{k})+\rho M^{T}M+\sigma
I\right)x=(\nabla^{2}f(x^{k})+\sigma I)x^{k}-\nabla f(x^{k})+\rho
M^{T}(z^{k}+c-u^{k}).$ (9)

The term $\sigma I$ ensures that this system is positive definite in the case
that $\nabla^{2}f(x^{k})+\rho M^{T}M$ is rank deficient. Generally, this is
not the case (for example, in machine learning problems, $M=I$), and the user
may set $\sigma=0$. For moderate problem sizes, we can use a direct method to
efficiently solve (9). In fact, if $\nabla^{2}f$ is constant (e.g., as in
conic programs (2)), an initial Cholesky factorization can be cached and
reused for all subsequent subproblems. However, for very large problems, even
the initial factorization can become infeasible to compute. We must resort to
the conjugate gradient method (CG) [HS+52], which only requires matrix-vector-
products. GeNIOS uses the CG implementation in Krylov.jl [MOc20].
Unfortunately, the left-hand-side matrix is often ill-conditioned due to the
low-rank nature of real-world data [UT19], and this ill-conditioning causes CG
to converge slowly. GeNIOS uses a Nyström preconditioner to ameliorate this
behavior.

#### Rank deficiency.

The left-hand-side matrix of (9) is usually full rank. For example,
$\nabla^{2}f(x)$ may be the a sum of a positive semidefinite matrix and a
positive diagonal matrix, or $M$ may include an identity block. In this case,
the user may set $\sigma=0$. If the rank is not obvious, the user can estimate
the minimum eigenvalue of $M^{T}M$ via power-iteration or the randomized
Lanczos method [MT20, Alg. 5]. For a conic program (or QP), the user can also
estimate the minimum eigenvalue of $\nabla^{2}f(x)=P$, which is constant
across iterations. These algorithms are relatively cheap to compute and can be
done as part of the problem setup. In general, however, GeNIOS cannot assume
that the left-hand-side matrix is full-rank and must use $\sigma>0$.

#### Inexact solves.

GeNIOS does not solve the system (9) exactly; it only ensures the $x$-update
solution error in (7) satisfies the error condition (8). For conic programs,
i.e., when the $x$-subproblem is not approximated, choosing the CG termination
criterion (the norms of the residuals) to be summable satisfies this
condition. GeNIOS uses the linear system relative tolerance

$\frac{\min\left(\sqrt{\|r_{p}^{k}\|\|r_{d}^{k}\|},\;1.0\right)}{k^{\gamma}}.$

This stopping tolerance combines strategies of existing conic solvers: SCS
[O’D+16] and COSMO [GCG21], which use $1/k^{\gamma}$ with $\gamma=1.5$ by
default; and OSQP’s GPU implementation [SBL20], which uses an error
proportional to the geometric mean of the residuals
$\sqrt{\|r_{p}^{k}\|\|r_{d}^{k}\|}$. Intuitively, GeNIOS uses a looser
tolerance when the residuals are large but then tightens this tolerance as the
residuals decrease to hasten convergence.

###  2.2 Randomized preconditioning

GeNIOS uses techniques from randomized numerical linear algebra to build a
Nyström preconditioner [FTU21] for the $x$-subproblem linear system (9). The
left-hand-side matrix is often ill-conditioned, resulting in slow convergence
of CG. To build a preconditioner, we want to quickly find an approximate
inverse of the dominant eigenspace of this matrix, which we do by creating a
random sketch that is easy to invert. When the matrix is not-too-sparse and
low rank, which is often true for real-world data (see [UT19] and references
therein), this preconditioner provides a significant speedup over standard CG.

#### Preliminaries.

Let $H$ denote the matrix in the $x$-subproblem linear system solve (9):

$H=\nabla^{2}f(x^{k})+\rho M^{T}M+\sigma I.$

In many practical settings, the optimization problem (1) can be written such
that either the Hessian or $M^{T}M$ is a diagonal matrix. GeNIOS uses a
preconditioner for matrices of the form

$A+D\qquad\text{where}\qquad A=\nabla^{2}f(x)\quad\text{or}\quad A=\rho
M^{T}M$

and $D$ is diagonal. Without loss of generality, we can assume that $D=\nu I$,
for some $\nu>0$ because the linear system

$(A+D)w=b$

is equivalent to the linear system

$\left(D^{-1/2}AD^{-1/2}+I\right)\tilde{w}=D^{-1/2}b,$

where $w=D^{-1/2}\tilde{w}$. For the remainder of this section, let $D=\nu I$
and assume the linear system in (9) has the form

$(A+\nu I)x=b.$ (10)

#### The preconditioner.

To precondition the linear system with left-hand-side matrix $A+\nu I$, GeNIOS
first constructs a randomized Nyström approximation [MT20, Alg. 16] to $A$
using test matrix $\Omega\in\mathbf{R}^{n\times r}$:

$\hat{A}=(A\Omega)(\Omega^{T}A\Omega)^{\dagger}(A\Omega)^{T}=U\hat{\Lambda}U^{T},$

where $U\in\mathbf{R}^{n\times r}$ has orthonormal columns, and
$\hat{\Lambda}\in\mathbf{R}^{r\times r}$ is diagonal. GeNIOS uses a standard
normal111 This choice does not preserve sparsity, and other options like
subsampled trigonometric transforms, may provide better performance for sparse
matrices. See [MT20, §9] for discussion and references.  test matrix $\Omega$,
as extensive theoretical and numerical work has found that this choice yields
an excellent low-rank approximation [AM15, MM17, Tro+17, Tro+19]. Given the
low-rank assumption holds, GeNIOS can take $r\ll n$ without losing much
accuracy in the approximation to $A$. Computing the sketch then approximately
scales as the cost of a matrix-vector-product with $A$. The randomized Nyström
preconditioner is

$L^{-1}=(\hat{\Lambda}_{r,r}+\nu I)U\left(\hat{\Lambda}+\nu
I\right)^{-1}U^{T}+I-UU^{T}.$ (11)

This preconditioner (approximately) inverts the dominant eigenspace of $A$,
while leaving the orthogonal complement unaffected. Additionally, $L^{-1}$
never needs to be formed explicitly; GeNIOS stores it in a factored form with
a light storage footprint. In this form, GeNIOS cheaply applies the
preconditioner in $O(nr)$ time and updates the parameter $\nu$ without
recomputing the Nyström approximation.

#### Adaptive sketch size selection.

A good preconditioner decreases the condition number of the linear system (10)
to a small constant. When this is the case, the standard analysis of CG
[TBI97, §38] implies that after only $O\left(\log(1/\varepsilon)\right)$
iterations, the CG iterate is within an $\varepsilon$-ball of the solution. In
practice, a user does not know the spectrum of $A$, so choosing a good sketch
size a-priori is difficult. However, after building the preconditioner, we can
explicitly bounded the condition number $\kappa$ of the left-hand-side matrix
[FTU21, Prop. 5.3]:

$\kappa(L^{-1/2}(A+\nu I)L^{-1/2})\leq
1+\frac{\hat{\Lambda}_{rr}+\|E\|}{\nu}\qquad\text{where }E=A-\hat{A}.$

This bound motivates an adaptive technique. Starting from some small initial
sketch size, GeNIOS can increase the size of the sketch until $\|E\|$ and
$\hat{\Lambda}_{rr}$ are suitably small, or until the sketch size is
unacceptably large. GeNIOS computes $\hat{\Lambda}_{rr}$ as part of the
Nyström sketch, and it can reliably estimate the error $\|E\|$ using a few
iterations of the randomized power method [MT20, Alg. 4]. See [FTU21, §5.4]
for additional discussion. By default GeNIOS uses a constant sketch size, but
the adaptive technique described above is also implemented.

###  2.3 Convergence

The algorithm implemented by GeNIOS is a special case of the GeNI-ADMM
framework in [Fra+23], which converges at the $O(1/k)$ rate of classic ADMM
despite the subproblem approximations and inexact solves, provided a solution
exists [Fra+23, Thm. 5]. Specifically, given that the subproblem errors (7)
satisfy the condition (8), the averaged iterates

$\bar{x}^{k+1}=\frac{1}{k}\sum^{k+1}_{t=2}x^{t}\qquad\text{ and
}\qquad\bar{z}^{k+1}=\frac{1}{k}\sum^{k+1}_{t=2}z^{t}$

produce objective value error and primal residual error that converge at rate
$O(1/k)$, i.e.,

$f(\bar{x}^{k+1})+g(\bar{z}^{k+1})-p^{\star}=O\left(1/k\right),$

and

$\|M\bar{x}^{k+1}+\bar{z}^{k+1}-c\|=O\left(1/k\right).$

The theory suggests that GeNIOS will require approximately the same number of
iterations as standard ADMM (demonstrated in [Fra+23]) and that each iteration
will be faster due to the approximation of the $x$-subproblem and inexact
subproblem solves, so overall we expect GeNIOS to converge faster than
standard ADMM. Our numerical experiments corroborate these expectations.

#### Optimality conditions.

Optimality conditions for the general problem (1) are primal feasibility,

$Mx^{\star}-z^{\star}=c,$

and that the Lagrangian, which can be written as

$\mathcal{L}(x,z,y)=f(x)+g(z)+\rho u^{T}(Mx-z-c),$

has a vanishing gradient when evaluated at the optimal primal and dual
variables:

$\displaystyle 0$ $\displaystyle=\nabla f(x^{\star})+\rho M^{T}u^{\star}$
$\displaystyle 0$ $\displaystyle\in\partial g(z^{\star})-\rho u^{\star}.$

In classic ADMM, $z^{k+1}$ and $u^{k+1}$ always satisfy the second condition
above exactly [Boy+11, §3.3]. However, in GeNIOS’s approximate version, this
condition no longer holds exactly; instead, $\rho u^{k+1}$ is almost a
subgradient of $g$ at $z^{k+1}$. The meaning of this last statement is made
precise in the following lemma.

######  Lemma 1 (Approximate optimality condition [Fra+23, Lem. 4.1]).

At each iteration $k$, there exists $s^{k}$ with

$\|s^{k}\|\leq\sqrt{\frac{2\varepsilon^{k}_{z}}{\rho}},$

such that the approximate $z$-subproblem solution $z^{k+1}$ satisfies

$\rho u^{k+1}+s^{k}\in\partial_{\varepsilon_{z}^{k}}g(z^{k+1}).$

In words, $\rho u^{k+1}+s^{k}$ belongs to the _$\varepsilon^{k}_{z}$
-subdifferential_ at $z^{k+1}$, which means

$g(z)-g(z^{k+1})\geq\langle\rho
u^{k+1}+s^{k},z-z^{k+1}\rangle-\varepsilon_{z}^{k},\quad\forall
z\in\mathbf{R}^{m}.$

The lemma shows that at each iteration, $\rho u^{k+1}$ is nearly a subgradient
of $g$ at $z^{k+1}$. As $\varepsilon_{z}^{k}$ decays quickly to zero, the
error in this second optimality condition quickly becomes small. Hence for
purposes of termination criteria, GeNIOS only needs to monitor the first
condition: $\nabla f(x^{k})+\rho M^{T}u^{k}$.

#### Termination criteria.

Based on the above discussion, we define the primal and dual residuals of
problem (1) as

$\displaystyle r_{\mathrm{p}}^{k}$ $\displaystyle=Mx^{k}-z^{k}-c,$ (12)
$\displaystyle r_{\mathrm{d}}^{k}$ $\displaystyle=\nabla f(x^{k})+\rho
M^{T}u^{k}.$ (13)

Under the inexactness assumptions discussed above, these residuals converge to
$0$ as $k\to\infty$. GeNIOS’s termination criteria are a mix of an absolute
and relative criterion based on these residuals:

$\displaystyle\|r_{\mathrm{p}}^{k}\|$
$\displaystyle\leq\sqrt{m}\varepsilon_{\mathrm{abs}}+\varepsilon_{\mathrm{rel}}\max\\{\|Mx^{k}\|,\|z^{k}\|,\|c\|\\},$
$\displaystyle\|r_{\mathrm{d}}^{k}\|$
$\displaystyle\leq\sqrt{n}\varepsilon_{\mathrm{abs}}+\varepsilon_{\mathrm{rel}}\|\rho
M^{T}u^{k}\|,$

where $\|\cdot\|$ is any norm. In some cases, the user may wish to use another
convergence criterion, which GeNIOS supports. For example, in machine learning
problems, GeNIOS supports a duality gap criterion, provided the user specifies
the conjugate function of the per-sample loss (see §3.3 for details).

###  2.4 Infeasibility detection

For generic convex optimization problems (1), GeNIOS’s iterates will diverge
if a solution does not exist; however, GeNIOS cannot produce a certificate of
infeasibility. For conic programs (2) on the other hand, GeNIOS can detect and
certify infeasiblity despite inexact solves. Our result follows [Ban+19] and
is a direct consequence of [RGN22, Thm. 2.1], which considered inexactness in
subproblem solves previously. For conic programs (2), the $x$-subproblem finds
an approximate solution, $x^{k}$ satisfying (8), to the linear system

$(P+\sigma I+\rho M^{T}M)x=\sigma x^{k}+\rho M^{T}(z^{k}-u^{k}-c)-q.$

The $z$-subproblem is simply a projection onto the cone $\mathcal{K}$,

$\Pi_{\mathcal{K}}(Mx^{k+1}+u^{k}-c),$

and GeNIOS produces an iterate $z^{k+1}$ satisfying

$\|\tilde{z}^{k+1}-(Mx^{k+1}+u^{k}-c)\|^{2}-\|\Pi_{\mathcal{K}}(Mx^{k+1}+u^{k}-c)-(Mx^{k+1}+u^{k}-c)\|^{2}\leq\varepsilon_{z}^{k},$

where we absorbed the constant $\rho/2$ into the error $\varepsilon_{z}^{k}$.
GeNIOS detects infeasibility by monitoring the sequences $\delta
x^{k}=x^{k+1}-x^{k}$ and $\delta u^{k}=u^{k+1}-u^{k}$.

######  Proposition 2 (Infeasibility certificate).

If the error sequences $\\{\varepsilon^{k}_{x}\\}_{k}$ and
$\\{\sqrt{\varepsilon^{k}_{z}}\\}_{k}$ are summable, then as $k\to\infty$, we
have $\delta x^{k}\to\delta x$ and $\delta u^{k}\to\delta u$. Further,

  1. 1.

If $\delta u\neq 0$, then (2) is primal infeasible. Additionally $\delta u$
provides a certificate of primal infeasibility, in that it satisfies

$M^{T}\delta u=0\quad\text{and}\quad S_{\mathcal{K}}(\delta u)<0,$ (14)

where $S_{\mathcal{K}}$ denotes the support function of $\mathcal{K}$, defined
as

$S_{\mathcal{K}}(\delta u)=\sup_{y\in\mathcal{K}}y^{T}\delta u.$

  2. 2.

If $\delta x\neq 0$, then (2) is dual infeasible. Additionally $\delta x$
provides a certificate of dual infeasibility, in that it satisfies

$P\delta x=0,\quad M\delta x\in\mathcal{K}^{\infty},\quad\text{and}\quad
q^{T}\delta x<0,$ (15)

where $\mathcal{K}^{\infty}$ is the recession cone of $\mathcal{K}$, defined
as

$\mathcal{K}^{\infty}=\\{y\in\mathbf{R}^{n}\mid x+\tau
y\in\mathcal{K}~{}\text{for all}~{}x\in\mathcal{K},\,\tau\geq 0\\}$

  3. 3.

If $\delta x\neq 0$ and $\delta y\neq 0$, then (2) is both primal and dual
infeasible, and $\delta x$ and $\delta u$ provide certificates of
infeasibility as above.

###### Proof.

$\rho$-Strong-convexity of the $z$-subproblem, along with the $z$-subproblem
inexactness condition implies that

$\|z^{k+1}-\Pi_{\mathcal{K}}(Mx^{k+1}+u^{k}-c)\|\leq\sqrt{\frac{2}{\rho}\varepsilon_{z}^{k}}.$

As $\sum_{k}\sqrt{\varepsilon_{z}^{k}}<\infty$ by the construction of the
GeNIOS algorithm, it follows that the errors in the solution to the
$z$-subproblem are summable. The same holds for the $x$-subproblem errors,
again by the construction of the GeNIOS algorithm. The desired result then
follows immediately from Theorem 2.1 in [RGN22], which proves the result in
the case where the $x$ and $z$-subproblem solution errors are summable. ∎

#### Algorithmic certificates.

Following [Ban+19, RGN22, Ste+20], we translate Proposition 2 into simple
algorithmic certificates of infeasibility. GeNIOS declares a problem to be
primal infeasible if the primal infeasibility certificate (14) holds
approximately:

$\|M^{T}\delta u^{k}\|<\varepsilon_{\textrm{inf}}\|\delta u^{k}\|,\quad
S_{\mathcal{K}}(\delta u^{k})<\varepsilon_{\textrm{inf}}\|\delta u^{k}\|,$

where $\varepsilon_{\textrm{inf}}$ is a positive tolerance. Similarly, GeNIOS
declares a problem to be dual infeasible if the dual infeasibility certificate
(15) holds approximately:

$\|P\delta x^{k}\|<\varepsilon_{\textrm{inf}}\|\delta
x^{k}\|,\quad\textrm{dist}_{\mathcal{K}^{\infty}}(M\delta
x^{k})<\varepsilon_{\textrm{inf}}\|\delta x^{k}\|,\quad q^{T}\delta
x^{k}<\varepsilon_{\textrm{inf}}\|\delta x^{k}\|.$

For the case of QPs, the support functions and recession cone for the
hyperrectangle $[l,u]\subseteq\mathbf{R}^{n}$ are well-defined, even though
this set is not necessarily a cone.

###  2.5 Performance improvements

GeNIOS includes performance improvements which are known to speed up
convergence in practice and are implemented in many ADMM solvers.

#### Over-relaxation.

In the $z$\- and $u$-updates, GeNIOS replaces the quantity $Mx^{k+1}$ with

$\alpha Mx^{k+1}+(1-\alpha)(z^{k}+c),$

where $\alpha\in(0,2)$ is a relaxation parameter. Experiments in the
literature [Eck94, EF98] show empirically that $\alpha$ in the range
$[1.5,1.8]$ can improve convergence.

#### Adjusting the penalty parameter.

First-order algorithms, like ADMM, are sensitive to scaling of the problem
data and to the penalty parameter $\rho$. GeNIOS uses preconditioning to
moderate the impact of the problem data scaling and selects $\rho$ using the
simple rule [HYW00, WL01, Boy+11]

$\rho^{k+1}=\begin{cases}\tau\rho^{k}&\|r^{k}_{\mathrm{p}}\|>\mu\|r^{k}_{\mathrm{d}}\|\\\
\rho^{k}/\tau&\|r^{k}_{\mathrm{d}}\|>\mu\|r^{k}_{\mathrm{p}}\|\\\
\rho^{k}&\mathrm{otherwise}.\end{cases}$

Since GeNIOS uses an indirect method, these updates can be applied cheaply;
GeNIOS does not need to refactor a matrix. GeNIOS can also update the
preconditioner by simply changing a scalar parameter (see §2.2). In practice,
we find that the update only makes sense to apply every 25 iterations or so
(this parameter is adjustable by the user). In the scaled version of ADMM, the
(scaled) dual variable $u^{k}$ must also be updated when $\rho$ changes. These
penalty parameter updates must stop after some finite number of iterations for
convergence guarantees to hold.

##  3 Applications

In this section, we detail the three problem classes GeNIOS handles and
provide an example of the interface for each class: general convex programs
(as in (1)) with the GenericSolver, quadratic programs with the QPSolver, and
regularized machine learning problems with the MlSolver.

###  3.1 General convex programs

Recall that GeNIOS solves convex optimization problems of the form

minimize $\displaystyle f(x)+g(z)$ subject to $\displaystyle Mx-z=c,$

where the variables are $x\in\mathbf{R}^{n}$ and $z\in\mathbf{R}^{m}$, and the
problem data are the functions $f$ and $g$, the linear operator $M$, and the
vector $c$. GeNIOS, heavily leverages multiple dispatch in the Julia
programming language [Bez+17] to implement fast primitives for inexact ADMM.
The base GeNIOS implementation uses a fully generic interface, accessible
through GeNIOS’s GenericSolver, which we describe in this section. In the
subsequent sections, we will detail how this interface is specialized for
quadratic programs and machine learning problems, each of which has its own
interface. Multiple dispatch allows GeNIOS to make performance optimizations
for these problem subclasses while still using the infrastructure of the
GenericSolver.

To define an optimization problem, GeNIOS’s GenericSolver needs the following
ingredients:

  * •

The function $f:\mathbf{R}^{n}\to\mathbf{R}$, its gradient $\nabla
f:\mathbf{R}^{n}\to\mathbf{R}^{n}$, and, given a point $x\in\mathbf{R}^{n}$,
its Hessian-vector product (HVP)
$\nabla^{2}f(x):\mathbf{R}^{n}\to\mathbf{R}^{n}$.

  * •

The function $g:\mathbf{R}^{m}\to\mathbf{R}$ and its proximal operator,
$\mathbf{prox}_{g/\rho}:\mathbf{R}^{m}\to\mathbf{R}^{m}$.

  * •

The linear operator $M:\mathbf{R}^{n}\to\mathbf{R}^{m}$ and the vector
$c\in\mathbf{R}^{m}$.

Much of the performance optimization, beyond the algorithm described in §2,
comes from an efficient implementation of the HVP for the $x$-subproblem and
of the proximal operator for the $z$-subproblem. Often, these operators can
take advantage of problem structure. For example, the Hessian of a diagonal-
plus-low-rank matrix can be computed much more efficiently by keeping it in a
factored form (e.g., as in §4.5). Julia has a number of useful packages for
defining efficient linear operators for the HVP and $M$, including
LinearMaps.jl [KHc23], and for defining proximal operators, including
ProximalOperators.jl [Ste+22]. GeNIOS’s documentation includes examples using
these packages. The fact that GeNIOS approximates the subproblems has the
added benefit of often allowing the user to avoid a conic reformulation, which
may significantly increase the size of the problem.

#### Example.

The Lasso regression problem is

minimize $\displaystyle(1/2)\|Ax-b\|^{2}_{2}+\lambda\|x\|_{1},$ (16)

with variable $x\in\mathbf{R}^{n}$, and problem data $A\in\mathbf{R}^{N\times
n}$, $b\in\mathbf{R}^{N}$, and $\lambda\in\mathbf{R}_{+}$. In the form of (1),
this problem becomes

minimize $\displaystyle(1/2)\|Ax-b\|^{2}_{2}+\lambda\|z\|_{1}$ subject to
$\displaystyle x-z=0,$

where we have introduced a new variable $z\in\mathbf{R}^{n}$. To put this
problem into the general form, take

$f(x)=(1/2)\|Ax-b\|^{2}_{2}\qquad\text{and}\qquad g(z)=\lambda\|z\|_{1}.$

The gradient of $f$ is $\nabla f(x)=A^{T}(Ax-b)$, and its HVP is
$\nabla^{2}f(x):z\mapsto A^{T}Az$. (In this problem, the HVP is independent of
the current iterate $x$.) Often in Lasso problems, $A\in\mathbf{R}^{N\times
n}$ with $N\ll n$. Thus, it is more efficient to compute the Hessian by first
multiplying by $A$ then by $A^{T}$ instead of forming the $n\times n$ matrix
$A^{T}A$. The proximal operator of $g$ is the soft-thresholding operator,

$\mathbf{prox}_{g/\rho}(v)_{i}=\left(\mathop{\rm
argmin}_{\tilde{z}}\lambda\|x\|_{1}+(\rho/2)\|\tilde{z}-v\|_{2}^{2}\right)_{i}=\begin{cases}v_{i}-\lambda/\rho&v_{i}>\lambda/\rho\\\
0&\lvert{v_{i}}\rvert\leq\lambda/\rho\\\
v_{i}+\lambda/\rho&v_{i}<-\lambda/\rho.\end{cases}$

Finally, take $M=I$ and $c=0$. Julia represents the identity matrix $I$ as a
special UniformScaling type to efficiently dispatch linear algebra subroutines
at compile time. In other words, these routines recognize that $I$ is the
identity operator and can be computed in $O(1)$ time.

#### Code example.

To use GeNIOS’s GenericSolver, the user first defines the function $f$, its
gradient, the function $g$, and its proximal operator directly in the Julia
language.222 This is not optimized code! Performance in the examples appearing
in §3 has been sacrificed for the sake of clarity. For performant examples,
please see the GeNIOS documentation and the code for the experiments in §4.

⬇

# Assume A, b, N, n, lambda have been defined

p = (; A=A, b=b, lambda=lambda)

f(x, p) = 0.5 * norm(p.A*x .- p.b)^2

function grad_f!(g, x, p)

g .= p.A’*(p.A*x .- p.b)

end

\parg(z, p) = p.lambda*sum(abs, z)

soft_thresh(zi, kappa) = sign(zi) * max(0.0, abs(zi) - kappa)

function prox_g!(v, z, rho, p)

v .= soft_thresh.(z, p.lambda/rho)

end

Julia style guidelines use an exclamation point ! at the end of function names
that modify their arguments. The dot . next to a scalar function is used to
_broadcast_ this function to act elementwise over vector arguments. To access
problem data inside the functions, we define the named tuple p, which is
passed to each function as the last argument and will be later used to
construct the solver.

To define the HVP, the user implements GeNIOS’s HessianOperator type:

⬇

struct HessianLasso{T, S<:AbstractMatrix{T}} <: HessianOperator

A::S

vN::Vector{T}

end

function LinearAlgebra.mul!(y, H::HessianLasso, x)

mul!(H.vN, H.A, x)

mul!(y, H.A’, H.vN)

end

update!(::HessianLasso, ::Solver) = nothing

\parHf = HessianLasso(A, zeros(N))

The update! function updates the internal data of the HessianOperator object
based on the current iterate. For the Lasso problem, the Hessian is
independent of the current iterate $x^{k}$, but we define this function for
completeness. Here, the HVP uses two multiplications: one by $A$, followed by
one by $A^{T}$. It avoids additional allocation of memory by storing a cache
vector vN. (Other functions should do the same; these examples do not for the
sake of simplicity.) Finally, the user defines a GenericSolver by simply
putting all these ingredients together, and then the user solves the problem
with the solve! function.

⬇

solver = GeNIOS.GenericSolver(

f, grad_f!, Hf, # f(x)

g, prox_g!, # g(z)

I, zeros(n); # M, c: Mx - z = c

params=p

)

res = solve!(solver)

Here, $M$ does not need to be a matrix; it only needs to define its operation
(mul!) on a vector. Additional details are provided in the documentation. In
the rest of this section, we discuss how this solver is specialized for
certain subclasses of problems. We leave detailed examination of specific
problems and problem-specific performance improvements to §4.

###  3.2 Quadratic Programs

GeNIOS’s QPSolver solves constrained QPs of the form

minimize $\displaystyle(1/2)x^{T}Px+q^{T}x$ (17) subject to $\displaystyle
Mx=z$ $\displaystyle l\leq z\leq u,$

where $x\in\mathbf{R}^{n}$ and $z\in\mathbf{R}^{m}$ are variables, and
$P\in\mathbf{S}_{+}^{n}$, $q\in\mathbf{R}^{n}$, $M\in\mathbf{R}^{m\times n}$,
$l\in\mathbf{R}^{m}$, and $u\in\mathbf{R}^{m}$ are the problem data. This
problem can be cast in the form (1) by taking $f(x)=(1/2)x^{T}Px+q^{T}x$,
$g(z)=I_{[l,u]}(z)$, and $c=0$:

minimize $\displaystyle(1/2)x^{T}Px+q^{T}x+I_{[l,u]}(z)$ subject to
$\displaystyle Mx-z=0.$

#### ADMM subproblems.

In the quadratic program case, the $x$-subproblem requires the (approximate)
solution to the linear system

$(P+\rho M^{T}M+\sigma I)x=\sigma x^{k}-q+\rho M^{T}(z^{k}+c-u^{k}).$

Because $f$ is quadratic, the second-order approximation is exact when
$\sigma=0$. The $z$-subproblem is simply a projection onto a hyperrectangle
defined by $l$ and $u$, which can be computed exactly in $O(m)$ time. Thus,
the main computational bottleneck at each iteration is the solution to the
linear system defining the $x$-subproblem.

#### Solving the linear system.

In many problems, the left-hand-side matrix of this linear system is full-
rank, even with $\sigma=0$, and the problem data appears in the objective
matrix $P$. (Even in a standard form linear program, where $P=0$, the
constraint that $x\geq 0$ will force $M$ to be full column rank. Of course,
there are reformulations such that the matrix becomes rank deficient.) Thus,
by default, GeNIOS simply sketches $P$ and sets the regularization parameter
in the preconditioner to be $\rho$. If the problem data appears in the
constraints, GeNIOS could instead divide the system by $\rho$ and then sketch
$M^{T}M$, but this feature is not implemented as of this writing. In most
cases, problems can be reformulated so that the data appear either in the
constraints or the objective but not both.

#### Performance improvements.

The unique features of GeNIOS and the Julia programming language allow for
several performance optimizations. Often, the matrix $P$ has a more efficient
form for both storage and computation. For example, $P$ may be the sum of a
diagonal and low rank component, i.e., $P=D+FF^{T}$, where $D$ is diagonal and
$F\in\mathbf{R}^{n\times k}$ with $k\ll n$. In this case, the user can create
a custom object which stores $P$ using only $O(nk)$ numbers and computes a
matrix-vector product in $O(nk)$ time.

#### Interface.

QPs are a restricted class of problems compared to (1). This allows for a
simple interface, accessed through QPSolver: the user only needs to specify
$P$, $q$, $M$, $l$, and $u$ in (17). Alternatively, the user can access this
interface through the JuMP modeling language [DHL17]. The linear operators $P$
and $M$ only need matrix-vector products defined, so they can be implemented
to exploit problem structure. An example appears in §4.5.

#### Code example.

Reformulate the lasso problem (16) as a QP with the same solution (cf. (17)):

minimize $\displaystyle(1/2)\begin{bmatrix}x\\\
t\end{bmatrix}^{T}\begin{bmatrix}A^{T}A&\\\ &0\end{bmatrix}\begin{bmatrix}x\\\
t\end{bmatrix}+\begin{bmatrix}-A^{T}b\\\
\lambda\mathbf{1}\end{bmatrix}^{T}\begin{bmatrix}x\\\ t\end{bmatrix}$ subject
to $\displaystyle\begin{bmatrix}0\\\
-\infty\end{bmatrix}\leq\begin{bmatrix}I&I\\\
I&-I\end{bmatrix}\begin{bmatrix}x\\\
t\end{bmatrix}\leq\begin{bmatrix}\infty\\\ 0\end{bmatrix}.$

From this formulation, identify $P$, $q$, $M$, $l$, and $u$. The code to solve
this QP is below.

⬇

# Assume that A, b, m, n, lambda have all been defined

P = blockdiag(sparse(A’*A), spzeros(n, n))

q = vcat(-A’*b, lambda*ones(n))

M = [

sparse(I, n, n) sparse(I, n, n);

sparse(I, n, n) -sparse(I, n, n)

]

l = [zeros(n); -Inf*ones(n)]

u = [Inf*ones(n); zeros(n)]

solver = GeNIOS.QPSolver(P, q, M, l, u)

res = solve!(solver)

###  3.3 Machine Learning

GeNIOS’s MlSolver solves convex machine learning problem of the form

minimize
$\displaystyle\sum_{i=1}^{N}\ell(a_{i}^{T}x-b_{i})+\lambda_{1}\|x\|_{1}+(1/2)\lambda_{2}\|x\|_{2}^{2},$
(18)

where $\ell:\mathbf{R}\to\mathbf{R}_{+}$ is some convex, per-sample loss
function, and $\lambda_{1}$, $\lambda_{2}\in\mathbf{R}_{+}$ are regularization
parameters. The variable $x\in\mathbf{R}^{n}$ is often called the (learned)
model weights, $a_{i}\in\mathbf{R}^{n}$ the feature vectors, and
$b_{i}\in\mathbf{R}$ the responses for $i=1,\dots,N$. Put this problem in the
general form (1) by setting
$f(x)=\sum_{i=1}^{N}\ell(a_{i}^{T}x-b_{i})+(1/2)\lambda_{2}\|x\|_{2}^{2}$ and
$g(z)=\lambda_{1}\|z\|_{1}$, and then adding the constraint that $x=z$.

#### Defining the problem.

The problem (18) is defined by the per-sample loss function $\ell$, the
feature matrix $A\in\mathbf{R}^{N\times n}$ with rows $a_{i}^{T}$ for
$i=1,\dots,N$, the response vector $b\in\mathbf{R}^{N}$, and the
regularization parameters $\lambda_{1}$, $\lambda_{2}\in\mathbf{R}_{+}$. The
gradient and Hessian of $f(x)$ have the form

$\nabla
f(x)=A^{T}\ell^{\prime}(Ax-b)+\lambda_{2}x\quad\text{and}\quad\nabla^{2}f(x)=A^{T}\mathop{{\bf
diag}\left({\ell^{\prime\prime}(Ax-b)}\right)}A+\lambda_{2}I,$

where $\ell^{\prime}$ and $\ell^{\prime\prime}$ are applied elementwise.

#### Solving the linear system.

Since $M=I$, the linear system associated with this problem (9) is always
positive definite; at each iteration, GeNIOS solves the system

$\left(A^{T}\mathop{{\bf
diag}\left({{\ell^{\prime\prime}}^{k}}\right)}A+(\lambda_{2}+\rho)I\right)x=\left(A^{T}\mathop{{\bf
diag}\left({{\ell^{\prime\prime}}^{k}}\right)}A+\lambda_{2}I\right)x^{k}-A^{T}{\ell^{\prime}}^{k}-\lambda_{2}x^{k}-\rho(z^{k}-c+u^{k}),$

where ${\ell^{\prime}}^{k}$, ${\ell^{\prime\prime}}^{k}\in\mathbf{R}^{N}$ are
shorthand for $\ell^{\prime}$ and $\ell^{\prime\prime}$ applied elementwise to
the vector $Ax^{k}-b$. The matrix in this system clearly has minimum
eigenvalue at least $\rho+\lambda_{2}$. Thus, to construct the preconditioner,
GeNIOS sketches $A^{T}\mathop{{\bf
diag}\left({{\ell^{\prime\prime}}^{k}}\right)}A$ and adds a regularization
parameter of $\rho+\lambda_{2}$. For some problems, including lasso
regression, $\ell^{\prime\prime}$ is constant, and GeNIOS only need to sketch
the left-hand-side matrix once. However, for others, including logistic
regression (see §4.2), GeNIOS needs to re-sketch the matrix as the weights are
updated. We found that a preconditioner works well for many subproblem solves
after it is created, so GeNIOS only resketches the matrix occasionally.

#### Custom convergence criterion.

For machine learning problems, GeNIOS uses a bound on the duality gap of (18)
to determine convergence. In particular, the solver constructs a dual feasible
point $\nu$ for the dual function of an equivalent reformulation of (18) and
then uses the termination criterion

$\frac{\ell(x)-g(\nu)}{\min\left(\ell(x),\lvert{g(\nu)}\rvert\right)}\leq\varepsilon_{\mathrm{dual}}.$

We call the quantity on the left the _relative duality gap_ , as it is clearly
a bound for the same quantity evaluated at the optimal dual variable
$\nu^{\star}$. For a full derivation, see appendix A. GeNIOS also allows the
user to specify other custom criteria, such as the relative change in the loss
function or the norm of the gradient.

#### Interface.

To use GeNIOS’s machine learning interface MlSolver, the user defines the per-
sample loss function $\ell:\mathbf{R}\to\mathbf{R}_{+}$ and the nonnegative
regularization parameters $\lambda_{1}$ and $\lambda_{2}$, as well as the data
matrix $A$ and response vector $b$. GeNIOS can use forward-mode automatic
differentiation (using ForwardDiff.jl [RLP16]) to define the first and second
derivatives of $\ell$, or these can also be supplied directly by the user.
GeNIOS defaults to using the relative duality gap convergence criterion for
lasso, elastic net, and logistic regression problems, which have their own
interfaces. The user must provide the conjugate function of the per-sample
loss, $\ell^{*}$, to compute the duality gap for a custom loss $\ell$.

#### Code example.

The lasso problem (16) is easily recognized as a machine learning problem (18)
without transformation. The user may call the solver directly:

⬇

# Assume that A, b, m, n, lambda have all been defined

f(x) = 0.5*x^2

reg_l1 = lambda

reg_l2 = 0.0

fconj(x) = 0.5*x^2

solver = MLSolver(f, reg_l1, reg_l2, A, b; fconj=fconj)

res = solve!(solver; options=SolverOptions(use_dual_gap=true))

For the lasso, elastic net, and logistic regression problems, GeNIOS also
provides specialized interfaces that only require the regularization
parameter(s) to be specified. The Lasso interface is:

⬇

solver = GeNIOS.LassoSolver(lambda, A, b)

res = solve!(solver; options=SolverOptions(use_dual_gap=true))

##  4 Numerical experiments

parameter | default value
---|---
linear system offset $\sigma$ (§2.1) | 1e-6
linear system tolerance exponent $\gamma$ (§2.1) | $1.2$
sketch size (§2.2) | $\min(50,\,n/20)$
resketching frequency (§2.2) | every $20$ iterations
norm for residuals (§2.3) | $\ell_{2}$
stopping tolerances $\varepsilon_{\mathrm{abs}}$ and $\varepsilon_{\mathrm{rel}}$ (§2.3) | 1e-4
infeasibility tolerance $\varepsilon_{\mathrm{inf}}$ (§2.4) | 1e-8
over-relaxation parameter $\alpha$ (§2.5) | $1.6$
penalty update factor $\tau$ (§2.5) | $2$
penalty update threshold $\mu$ (§2.5) | $10$
Table 2: Default parameters for GeNIOS

In this section, we showcase the performance of GeNIOS on a variety of
problems with simulated and real-world data. The first four examples are all
machine learning problems. The final two examples come from finance and signal
processing respectively. Each highlights one or more specific features of
GeNIOS. For some examples, we compare GeNIOS against popular open-source, ADMM
solvers OSQP and COSMO, and against commercial, interior-point solver Mosek.
For other examples, we are primarily interested in comparing different GeNIOS
interfaces or options against each other. In brief, the numerical examples are
outlined below:

  * •

The _elastic net_ problem demonstrates the impact of GeNIOS’s Nyström
preconditioning and inexact subproblem solves on a dense, low-rank machine
learning problem.

  * •

The _logistic regression_ problem shows the impact of GeNIOS’s approximation
of the $x$-subproblem and of avoiding conic reformulation.

  * •

The _Huber regression_ problem compares GeNIOS’s MlSolver to a conic
reformulation.

  * •

The _constrained least squares_ problem compares GeNIOS to OSQP, COSMO, and
Mosek on a dense, low-rank QP.

  * •

The _portfolio optimization_ problem compares GeNIOS’s QPSolver and
GenericSolver to OSQP, COSMO, and Mosek on a sparse, structured QP.

  * •

Finally, the _signal decomposition_ problem demonstrates that GeNIOS’s
GenericSolver can be used to solve nonconvex problems as well. (Of course,
GeNIOS loses the convergence guarantees of the convex case.)

All numerical examples can be found in the benchmarks folder of the solver
repository,

https://github.com/tjdiamandis/GeNIOS.jl.

All experiments were run using GeNIOS v0.2.0 on a MacBook Pro with a M1 Max
processor (8 performance cores) and 64GB of RAM. Unless otherwise stated, we
use default parameters for GeNIOS, listed in table 2. The experiments employ
no parallelization except for multithreaded BLAS. Of course, GeNIOS’s use of
CG to solve the linear system can naturally benefit from additional
parallelization, and many proximal operators can be efficiently parallelized
as well. We leave exploration of GPU acceleration and other forms of
parallelism to future work. Similar examples to those in this section, with
toy data, can be found in the examples section of the documentation.

###  4.1 Elastic net

We first use the elastic net problem, a standard quadratic machine learning
problem, to highlight the advantages of randomized preconditioning and of
inexact subproblem solves. (Because the objective is quadratic, the
$x$-subproblem is not approximated here.) The elastic net problem is

minimize
$\displaystyle(1/2)\|Ax-b\|_{2}^{2}+\lambda_{1}\|x\|_{1}+(\lambda_{2}/2)\|x\|_{2}^{2},$

where $x\in\mathbf{R}^{n}$ is the variable, $A\in\mathbf{R}^{N\times n}$ is
the feature matrix, $b\in\mathbf{R}^{N}$ is the response vector, and
$\lambda_{1}$, $\lambda_{2}$ are regularization parameters. The lasso problem
can be recovered by setting $\lambda_{2}=0$. Clearly, this problem is a
special case of (1) with

$f(x)=(1/2)\|Ax-b\|_{2}^{2}+(\lambda_{2}/2)\|x\|_{2}^{2}\quad\text{and}\quad
g(x)=\lambda_{1}\|x\|_{1},$

and of the machine learning problem formulation (18) with
$\ell(w)=(1/2)w^{2}$. The second derivative is fixed across iterations, so
GeNIOS never needs to update the sketch used to build the preconditioner. (The
preconditioner itself may be updated if the penalty parameter $\rho$ changes.)
Furthermore, since the Hessian of $f$ is constant, the $x$-subproblem is not
approximated—the linear system is simply solved inexactly.

#### Problem data.

We solve the elastic net problem with both the sparse real-sim dataset [CL11]
and the dense YearMSD dataset [BM11], which we augment with random features
[RR07, RR08] as in [ZFU22]. The dataset statistics are summarized in table 3.
We set $\lambda_{2}=\lambda_{1}=0.1\|A^{T}b\|_{\infty}$. Both of these
datasets are approximately low rank, and we estimate the maximum and minimum
eigenvalues of the Gram matrix $A^{T}A$ using the randomized Lanczos method
[MT20, Alg. 5]. In both cases, the minimum eigenvalue is significantly less
than the penalty parameter $\rho=1$. As a result, the condition number of the
$\rho$-regularized linear system is approximately the maximum eigenvalue
$\lambda_{\mathrm{max}}$.

| samples $N$ | features $n$ | nonzeros | density |  $\lambda_{\mathrm{max}}$ (est.) |  $\lambda_{\mathrm{min}}$ (est.)
---|---|---|---|---|---|---
real-sim | 72.3k | 21.0k | 3.709M | 0.24% | 920.8 | 0.0020
YearMSD | 10k | 20k | 200M | 100% | 4450 | 0.0002
Table 3: Dataset statistics.

#### Experiments.

We examine the impact of the preconditioner and the inexact solves for the
elastic net problem with each of these datasets. Figure 1 shows the
convergence on the sparse real-sim dataset, and figure 2 shows the convergence
on the dense YearMSD dataset. For both datasets, inexact subproblem solves
speed up convergence by about 2$\times$. For the dense dataset, the
preconditioner reduces the time to solve the linear system by approximately
50%. For the sparse dataset, the preconditioner provides only a modest
reduction in the time to solve the linear system and so a modest improvement
in the overall solve time: the dataset is so sparse that applying the dense
preconditioner introduces considerable overhead compared to the low cost of a
CG iteration. Detailed timings are in tables 4 and 5. We also show a high
precision solve in figures 3(a) and 3(b), which illustrates the linear
convergence of GeNIOS past the stopping tolerances used in our comparisons.

(a) Relative duality gap

(b) Objective value

(c) Primal residual

(d) Dual residual

Figure 1: Inexact subproblem solves improve GeNIOS’s convergence time on the
elastic net problem with the real-sim dataset, but the preconditioner has
little effect since the dataset is very sparse.

(a) Relative duality gap

(b) Objective value

(c) Primal residual

(d) Dual residual

Figure 2: Both inexact subproblem solves and the preconditioner improve GeNIOS’s convergence time on the elastic net problem with the dense YearMSD dataset. | GeNIOS | GeNIOS (no pc) | ADMM | ADMM (no pc)
---|---|---|---|---
setup time (total) | 0.752s | 0.011s | 0.405s | 0.011s
preconditioner time | 0.741s | 0.000s | 0.394s | 0.000s
solve time | 13.848s | 15.026s | 26.396s | 29.156s
number of iterations | 190 | 190 | 190 | 190
total linear system time | 8.736s | 10.008s | 21.350s | 24.116s
avg. linear system time | 45.978ms | 52.674ms | 112.367ms | 126.925ms
total prox time | 0.006s | 0.005s | 0.005s | 0.005s
avg. prox time | 0.029ms | 0.028ms | 0.028ms | 0.029ms
total time | 14.600s | 15.037s | 26.801s | 29.167s
Table 4: Timings for GeNIOS with and without preconditioning and inexact
solves for the elastic net problem with the sparse real-sim dataset.

GeNIOS GeNIOS (no pc) ADMM ADMM (no pc) setup time (total) 2.195s 0.053s
1.835s 0.052s preconditioner time 2.142s 0.000s 1.786s 0.000s solve time
18.634s 30.491s 38.221s 62.756s number of iterations 50 50 50 50 total linear
system time 12.194s 23.872s 31.940s 56.278s avg. linear system time 243.879ms
477.442ms 638.796ms 1125.564ms total prox time 0.001s 0.001s 0.001s 0.001s
avg. prox time 0.027ms 0.027ms 0.029ms 0.029ms total time 20.830s 30.544s
40.055s 62.808s

Table 5: Timings for GeNIOS with and without preconditioning and inexact
solves for the elastic net problem with the dense YearMSD dataset.

(a) Elastic net, real-sim

(b) Elastic net, YearMSD

(c) Logistic regression, real-sim

Figure 3: GeNIOS converges linearly on the elastic net problem (a, b) and the
logistic regression problem (c) up to high accuracy. Here, eps denotes machine
epsilon.

###  4.2 Logistic Regression

This example highlights the benefits of using an approximate $x$-subproblem,
in addition to an inexact solve. For the logistic regression problem, with
problem data $\tilde{a}_{i}\in\mathbf{R}^{n}$ and $\tilde{b}_{i}\in\mathbf{R}$
for $i=1,\dots,N$, where $\tilde{b}_{i}\in\\{\pm 1\\}$, define

$\mathbb{P}\left({\tilde{b}_{i}\mid\tilde{a}_{i}}\right)=\frac{1}{1+\exp\left(\tilde{b}_{i}(\tilde{a}_{i}^{T}x)\right)}=\frac{1}{1+\exp(a_{i}^{T}x)},$

where $a_{i}=\tilde{b}_{i}\tilde{a}_{i}$. Use the negative of the log
likelihood as the loss, giving the optimization problem

minimize
$\displaystyle\sum_{i=1}^{m}\log\left(1+\exp(a_{i}^{T}x)\right)+\lambda_{1}\|x\|_{1}.$

Recognize that, in the form of (18), the per-sample loss function is
$\ell(w)=\log(1+\exp(w))$.

#### Experiments.

We use the real-sim dataset as in the elastic net experiment (see table 3),
which corresponds to a binary classification problem and take
$\lambda_{1}=(0.1/2)\|A^{T}\mathbf{1}\|_{\infty}$. We solve this problem with
three distinct methods: GeNIOS’s MlSolver, with and without both
preconditioning and inexact solves; the MlSolver with an exact $x$-update,
which we solve with Optim.jl’s [MR18] implementation of L-BFGS [Noc80, LN89]
with default parameters333 The stopping criterion is when the infinity norm of
the gradient is under 1e-8. ; and GeNIOS’s QPSolver, which we modified to
handle exponential cone constraints. In the conic form problem, we solve the
$x$-subproblem inexactly and we use the fast projection onto the exponential
cone from Friberg [Fri23] for the $z$-subproblem. In our formulation, the
conic form problem has $2n+5N$ variables and $2n+9N$ constraints instead of
the $2n$ variables and $n$ constraints in the MlSolver formulation (see
appendix B for the equivalent conic problem). We show convergence in figure 4
and breakdown timing in table 6. The conic form solve is not plotted because
its residuals refer to different quantities, making direct comparison
difficult. In this example, using the preconditioner does not add more of a
speedup than its overhead since the dataset is quite sparse. Standard ADMM is
slower than GeNIOS, but it is faster than using an exact solve for the
approximate $x$-subproblem. The conic form problem has the slowest solve time,
and its per-iteration time is longer than GeNIOS with inexact solves.

#### Discussion.

The results in table 6 tell an interesting story. First, as in the previous
example, inexact solves of the approximate $x$-subproblem have little impact
on the number of iterations GeNIOS takes to converge. However, there is a
clear trade off between approximating the $x$-subproblem and solving the exact
problem: the exact problem takes longer to solve but significantly cuts down
the number of ADMM iterations the algorithm takes to converge. The conic
form’s $x$-update is also faster than standard ADMM with an exact solve, as
the linear system is solved inexactly and the constraint matrix is very sparse
and structured. However, the algorithm pays for this speed in the $z$-update,
which requires computing $2N$ projections onto the exponential cone. Overall,
these results highlight the benefits of solving logistic regression in a more
natural form, such as (18), instead of the conic form, and of approximating
the $x$-subproblem to speed up solve time, even if the number of iterations
required to converge increases.

(a) Relative duality gap convergence

(b) Objective value convergence

(c) Primal residual convergence

(d) Dual residual convergence

Figure 4:  Both subproblem approximation and inexact solves improve GeNIOS’s
convergence time on the logistic regression problem with the real-sim dataset.
Similarly to the elastic net problem, the preconditioner has little effect
since this dataset is very sparse.

GeNIOS G. (no pc) G. (exact) G. (no pc, exact) ADMM LBFGS Conic setup time
(total) 0.405s 0.012s 0.394s 0.011s 0.009s 0.003s preconditioner time 0.383s
0.000s 0.385s 0.000s 0.000s 0.000s solve time 11.692s 9.565s 40.020s 43.109s
23.898s 59.698s number of iterations 136 140 137 137 45 127 total $x$-update
time 5.392s 5.593s 33.521s 39.017s 22.883s 20.491s avg. $x$-update time
39.650ms 39.951ms 244.677ms 284.794ms 508.514ms 161.343ms total prox time
0.005s 0.005s 0.005s 0.005s 0.002s 34.938s avg. prox time 0.038ms 0.037ms
0.039ms 0.040ms 0.034ms 275.099ms total time 12.097s 9.577s 40.414s 43.119s
23.907s 59.701s

Table 6: Timings for GeNIOS with and without preconditioning, subproblem
approximations, and inexact solves for the logistic regression problem with
the real-sim dataset. We also compare with the conic form problem.

###  4.3 Huber fitting

In this example, we showcase GeNIOS’s ability to handle custom loss functions,
which are often not supported by standard machine learning solvers. Handling
these custom loss functions directly, instead of solving the equivalent conic
program, speeds up the solution of machine learning problems, which we
demonstrate for the case of Huber fitting [Hub64]. The Huber fitting problem
replaces the standard squared-error loss function (cf. the elastic net problem
in §4.1) with a loss function that is less sensitive to outliers, defined as

$\ell^{\mathrm{hub}}(w)=\begin{cases}w^{2}&\lvert w\rvert\leq 1\\\
2|w|-1&\lvert w\rvert>1.\end{cases}$

This function is easily verified to be convex and smooth. The
$\ell_{1}$-regularized Huber fitting problem is then

minimize
$\displaystyle\sum_{i=1}^{N}\ell^{\mathrm{hub}}(a_{i}^{T}x-b_{i})+\lambda_{1}\|x\|_{1},$

with problem data $a_{i}\in\mathbf{R}^{n}$ and $b_{i}\in\mathbf{R}$ for
$i=1,\dots,N$, variable $x\in\mathbf{R}^{n}$, and regularization parameter
$\lambda_{1}\geq 0$. Huber fitting, a form of _robust regression_ , often has
superior performance on real-world data and obviates the need for outlier
detection in data pre-processing. However, the Huber loss and other robust
loss functions are usually not supported by standard solvers. GeNIOS’s
MlSolver, on the other hand, provides full support for custom convex loss
functions, as described in §3.3.

#### Equivalent quadratic program.

The $\ell_{1}$-regularized Huber fitting problem can also be written as a
quadratic program (see [MM00] for details) by introducing new variables
$q\in\mathbf{R}^{n}$, $r\in\mathbf{R}^{N}$, $s\in\mathbf{R}^{N}$, and
$t\in\mathbf{R}^{N}$:

minimize $\displaystyle r^{T}r+2\mathbf{1}^{T}(s+t)+\lambda_{1}q$ subject to
$\displaystyle Ax-r-s+t=b$ $\displaystyle-q\leq x\leq q$ $\displaystyle 0\leq
s,t$

where $A$ is a $N\times n$ matrix with rows $a_{i}^{T}$ and $b$ is a vector
with elements $b_{i}$. This problem now has $2n+3N$ variables, and the
constraint matrix in (17) is of size $2n+3N\times 2n+3N$. Putting other robust
regression problems into standard forms can also lead to substantial increases
in problem size.

#### Problem data.

For this problem, we generate random data with $n$ features and $N=n/2$
samples for varying values of $n$. We sample each feature $A_{i,j}$
independently from the standard normal distribution $\mathcal{N}(0,1)$ and
then normalize the columns to have zero mean and unit $\ell_{2}$ norm. The
true value of weights $x^{\star}$ has 10% nonzero entries, each sampled from
$\mathcal{N}(0,1)$. The response vector $b$ is computed as $Ax^{\star}+0.1v$,
where $v\sim\mathcal{N}(0,1)$. Finally, we make 5% of $b$ outliers by adding
$u$ to these entries, where each $u$ is drawn uniformly at random from
$\\{-10,10\\}$. We take $\lambda_{1}=0.1\|A^{T}b\|_{\infty}$. We vary $n$ from
250 to 16,000.

#### Comparing solver interfaces.

First, we examine the difference between using the GeNIOS MLSolver interface
and using the QPSolver interface. Since the residuals for the MLSolver and
QPSolver interfaces are different, we overwrite the convergence criteria
calculation to examine the subdifferential of the original loss function:

$A^{T}{\ell^{\mathrm{hub}}}^{\prime}(Ax-b)+\lambda\partial\|x\|_{1},$

where the derivative of the loss is applied elementwise and

$\partial\|x\|_{1}=\left\\{v\mid
v_{i}=\mathbf{sign}(x_{i})~{}\text{if}~{}\lvert{x_{i}}\rvert>0~{}\text{and
otherwise}~{}v_{i}\in[-1,1]\right\\}.$

We say that the solver has converged to an optimal solution when the distance
from the vector $0$ to the subdifferential set is less than 1e-4, i.e., when
there is some vector in this set with $\ell_{2}$ norm less than 1e-4. We do
not use randomized preconditioning for this problem, as the random data matrix
$A$ is not approximately low rank. We solve the problem for $n$ varied from
$250$ to $16,000$. Figure 5 shows timing for the linear system solve in the
$x$-subproblem and for the overall solve for varying $n$, and it shows the
residuals’ convergence as a function of cumulative time spent on solving the
linear system for $n=16,000$. We provide detailed solve time breakdowns in
table 7 for $n=16,000$.

#### Discussion.

The MlSolver speeds up solve time by over an order-of-magnitude compared to
the QPSolver, even for small problem sizes. The dominant operation per
iteration—the linear system solve—is over an order of magnitude faster (table
7). MlSolver is also faster than QPSolver even after subtracting the time to
solve the linear system. Both methods require about the same number of
iterations to reach the stopping criterion. Multiplication by $A$ and $A^{T}$
are two most expensive operations in the algorithm. MlSolver represents $A$ as
a dense matrix and uses optimized BLAS operations to compute matrix-vector
products. In contrast, QPSolver stores the data matrix $A$ as part of the
large, sparse constraint matrix $M$ (see (17)) and so cannot apply $A$ and
$A^{T}$ efficiently.

(a) Solve time comparisons for $n=250$ to $n=16,000$.

(b) Primal and Dual residual convergence for $n=16,000$.

Figure 5: GeNIOS’s MlSolver avoids conic reformulation and significantly
outperforms the QPSolver (with and without exact subproblem solves) for the
Huber fitting problem.

MLSolver QPSolver QPSolver (exact) setup time (total) 0.017s 0.000s 0.000s
preconditioner time 0.000s 0.000s 0.000s solve time 14.717s 483.966s 697.787s
number of iterations 74 96 92 total linear system time 9.727s 176.815s
402.297s avg. linear system time 131.441ms 1841.822ms 4372.798ms total prox
time 0.002s 0.006s 0.006s avg. prox time 0.025ms 0.067ms 0.066ms total time
14.734s 483.966s 697.788s

Table 7: Timing comparisons of GeNIOS’s MlSolver and QPSolver for the Huber
fitting problem for $n=16,000$.

###  4.4 Bounded Least Squares

We compare GeNIOS one-to-one with other solvers on the bounded least squares
problem. This problem solves a standard least squares problem with upper and
lower bounds on the weights:

minimize $\displaystyle(1/2)\|Ax-b\|^{2}_{2}$ subject to $\displaystyle 0\leq
x\leq 1.$

The problem is clearly a QP (17) with $P=A^{T}A$, and $q=A^{T}b$. Unlike in
the unconstrained machine learning problems, the $x$-subproblem iterates are
not necessarily feasible; however, the $z$-subproblem iterates are.

#### Problem data.

We use the first $N$ samples in the YearMSD dataset and generate $n=N/2$
random features. We precompute $P$ and $q$. The matrix in the linear system
solve, $P+\rho I$, is dense. Since the matrix $P$ comes from real world data,
we expect it to also be low-rank. The fact that $P$ is low-rank and dense
suggest that using a randomized preconditioner will lead to significant
speedups, similar to those we saw when solving the elastic net problem in
§4.1.

#### Ablation study.

First, we examine the impact of the preconditioner and inexact solves for this
problem. (Since this problem is a QP, GeNIOS does not approximate the
$x$-subproblem.) We solve the bounded least squares problem with $N=20,000$
data samples and $n=N/2$ features with and without both the randomized
preconditioner and inexact $x$-subproblem solves. We set the absolute and
relative termination tolerances to be 1e-5. Figure 6 shows the convergence of
the primal residual, the dual residual, and the objective value. The inexact
solves reduce solve time by over 50%, and the preconditioner introduces a very
modest overhead (under 2%) for an approximately 40% linear system solve time
reduction. Table 8 details these solve times. Again, the addition of inexact
solves and the preconditioner do not affect the number of iterations it takes
for the solution to converge. As with many QPs, the linear system solve time
dominates.

(a) Objective value

(b) Primal residual

(c) Dual residual

Figure 6:  Both inexact subproblem solves and the preconditioner improve GeNIOS’s convergence time on the bounded least squares problem with the dense YearMSD dataset (cf. the elastic net problem in figure 2).  | GeNIOS | GeNIOS (no pc) | ADMM | ADMM (no pc)
---|---|---|---|---
setup time (total) | 0.859s | 0.012s | 0.601s | 0.011s
preconditioner time | 0.848s | 0.000s | 0.589s | 0.000s
solve time | 27.651s | 53.465s | 48.549s | 94.418s
number of iterations | 68 | 68 | 67 | 67
total linear system time | 25.507s | 50.800s | 46.489s | 91.781s
avg. linear system time | 375.098ms | 747.057ms | 693.869ms | 1369.870ms
total prox time | 0.002s | 0.002s | 0.002s | 0.002s
avg. prox time | 0.023ms | 0.023ms | 0.025ms | 0.024ms
total time | 28.510s | 53.476s | 49.150s | 94.429s
Table 8: Timing comparisons of GeNIOS’s MlSolver and QPSolver for the
constrained least squares problem with $n=16,000$ features (augmented) and
$N=32,000$ samples from the YearMSD dataset.

#### Comparison with other solvers.

We solve the same problem, with varying value of $n$ from $250$ to $16,000$,
with our solver, the popular QP solver OSQP [Ste+20], the pure Julia conic
solver COSMO [GCG21], and the commercial solver Mosek [ApS22]. Both OSQP and
COSMO are ADMM-based solvers, while Mosek uses an interior point method. For
COSMO, we use both the QDLDL direct solver (the default) and the CG indirect
solver, which solves the same reduced system as GeNIOS’s QPSolver. We set all
ADMM-based solvers to have absolute and relative termination tolerances of
1e-4 and use the infinity norm of the residuals, since this is the default in
OSQP and COSMO. We use Mosek’s default parameters. Both GeNIOS and COSMO’s
indirect solver use inexact solves in this example. Figure 7(a) shows the
results. GeNIOS with and without preconditioning begins to perform much better
than other methods as the problem size becomes large. Preconditoning provides
a nontrivial solve time reduction—about 50%—as the problem sizes grows. While
OSQP and COSMO with a direct linear system solve are faster than GeNIOS for
smaller problem sizes, they have worse scaling as the problem size grows due
to the matrix factorization.

(a) Bounded least squares

(b) Portfolio optimization

Figure 7: GeNIOS outperforms the COSMO, OSQP, and Mosek solvers for the dense
bounded least squares problem (left). For the sparse portfolio optimization
problem, GeNIOS outperforms direct solvers at large problem sizes and enjoys
better scaling with problem data size.

###  4.5 Portfolio optimization

The portfolio optimization problem highlights how GeNIOS’s QPSolver allows the
user to speed up standard QP’s by leveraging Julia’s multiple dispatch, and
how GeNIOS’s GenericSolver provides powerful tools for advanced users to
exploit structure. The portfolio optimization problem finds the fraction of
wealth to be invested across a universe of $n$ assets to maximize risk-
adjusted return. If the portfolio is constrained to be long-only, this problem
can be written as

minimize $\displaystyle-\mu^{T}x+(\gamma/2)x^{T}\Sigma x$ subject to
$\displaystyle\mathbf{1}^{T}x=1$ $\displaystyle x\geq 0,$

where the variable $x\in\mathbf{R}^{n}$ represents the allocation,
$\mu\in\mathbf{R}^{n}$ and $\Sigma\in\mathbf{S}_{+}^{n}$ are the (estimated)
return mean and covariances respectively, and $\gamma\in\mathbf{R}_{++}$ is a
risk-aversion parameter. Often, the covariance matrix $\Sigma$ has a diagonal-
plus-low-rank structure, i.e.,

$\Sigma=FF^{T}+D,$

where $D\in\mathbf{R}^{n\times n}$ is a diagonal matrix indicating asset-
specific risk and $F\in\mathbf{R}^{n\times k}$ is a factor matrix with $k\ll
n$. This problem is clearly a QP as written. To fit into the form of (17), set

$P=\gamma\Sigma,\qquad q=-\mu,\qquad M=\begin{bmatrix}\mathbf{1}^{T}\\\
I\end{bmatrix}\ \qquad l=\begin{bmatrix}1\\\ 0\end{bmatrix},\qquad
u=\begin{bmatrix}1\\\ \infty\end{bmatrix}.$

We solve the portfolio optimization problem using both GeNIOS’s QPSolver and
its GenericSolver, which permits additional performance improvements. We
compare solve times against OSQP and COSMO with absolute and relative
tolerances set to 1e-4, and against Mosek with its default parameters. Again,
we use the infinity norm of the residuals in GeNIOS for the sake of
comparison.

#### Equivalent QP.

In portfolio optimization problem, the right-hand-side matrix of the
$x$-subproblem linear system is a dense matrix; the diagonal-plus-low-rank
structure of $\Sigma$ is lost. To take advantage of structure, the portfolio
optimization problem can be reformulated as the following equivalent quadratic
program:

minimize $\displaystyle(\gamma/2)x^{T}Dx+(\gamma/2)y^{T}y-\mu^{T}x$ subject to
$\displaystyle y=F^{T}x$ $\displaystyle\mathbf{1}^{T}x=1$ $\displaystyle x\geq
0.$

In this equivalent formulation, the matrix $\Sigma$ is no longer formed. We
solve this QP instead of the original QP unless otherwise stated.

#### QP interface.

While we can solve the equivalent problem introduced above, GeNIOS allows us
to avoid this reformulation. Instead, we create types for $P$ and the
constraint matrix $M$ that implement fast multiplication by taking advantage
of structure. We can multiply by $P$ in $O(nk)$ time, and we can multiply by
$M$ or $M^{T}$ in $O(n)$ time. Julia’s multiple dispatch allows the user to
easily define new objects that implement fast operations. We provide a
tutorial in the ‘Markowitz Portfolio Optimization, Three Ways’ example in the
documentation.

#### Generic interface.

The constraint set of the portfolio optimization problem is the
$n$-dimensional simplex, for which there exists a fast projection. We can take
advantage of this fast projection using GeNIOS’s GenericSolver. Instead of
using the QP formulation (17), we formulate the portfolio optimization problem
as

minimize $\displaystyle-\mu^{T}x+(\gamma/2)x^{T}\Sigma x+I_{S}(z)$ subject to
$\displaystyle x-z=0,$

where $I_{S}(z)$ is the indicator function of the simplex
$S=\\{z\mid\mathbf{1}^{T}z=1\text{ and }z\geq 0\\}$. Recognize this
formulation as (1) with $f(x)=(\gamma/2)x^{T}\Sigma x-\mu^{T}x$,
$g(z)=I_{S}(z)$, $M=I$, and $c=0$. The proximal operator of $g$ evaluated at
$v$ is a projection of $v$ onto the set $S$, defined as the solution to the
optimization problem

minimize $\displaystyle(1/2)\|\tilde{z}-v\|_{2}^{2}$ subject to
$\displaystyle\mathbf{1}^{T}\tilde{z}=1$ $\displaystyle\tilde{z}\geq 0.$

This problem can be efficiently solved by considering the optimality condition
for its dual (see appendix C):

$\mathbf{1}^{T}(v-\nu\mathbf{1})_{+}=1.$

The primal solution $\tilde{z}$ can be reconstructed as

$\tilde{z}_{i}=(v_{i}-\nu)_{+}.$

Thus, the proximal operator of $g$ reduces to a single-variable root finding
problem, which can be efficiently solved via bisection search. In this
experiment, we solve this problem to an accuracy of 1e-8. While $z$-subproblem
inexactness is also permitted in GeNIOS (see §2.3), we do not fully explore
the implications in this work.

#### Problem data.

For this problem, we use synthetic data generated as in [Ste+20, App. A]. The
asset specific risk, i.e., the entries of the diagonal matrix $D$, are sampled
independently and uniformly from the interval $[0,\sqrt{k}]$. The factor
loading matrix $F$ is sparse, with half of the entries randomly selected to be
independently sampled standard normals and the other half to be zero. We set
$n=100k$. Finally, the return vector $\mu$ also has independent standard
normal entries. We set the risk parameter $\gamma$ to be $1$. We vary the
number of assets $n$ from $250$ to $256,000$.

#### Numerical results.

For each value of $n$, we solve the portfolio optimization problem with
GeNIOS, COSMO, OSQP, and Mosek. We use GeNIOS’s QPSolver to solve the original
QP, the equivalent QP, and the original problem with custom $P$ and $M$
operators. We use GeNIOS’s GenericSolver to solve the equivalent problem where
we deal with the constraint set directly. We use both COSMO’s indirect and
direct linear system solvers. The results are shown in figure 7(b). Clearly
solving the original, full QP, is a bad idea; the matrix $\Sigma$ becomes
dense and much of the structure is lost, resulting in a slow solve time. When
solving the equivalent QP, GeNIOS outperforms COSMO’s indirect solver (which
also uses inexact solves), but both of these solvers are slower than COSMO’s
direct solver and OSQP. In this problem, the constraint matrix in the
equivalent formulation is sparse (approximately $0.5\%$ nonzeros as $k$ gets
large) and highly structured. As a result, sparse factorization methods
perform relatively well, even for large problem sizes (over $80$M nonzeros).
However, the direct method solve times increase at a faster rate, so for a
large enough problem size, indirect methods will perform better than direct
methods. Solving the original QP using custom operators with GeNIOS becomes
competitive with the direct methods much faster as the problem size gets
large, and it also scales at a slower rate. Finally, GeNIOS’s GenericSolver,
which exploits additional structure, begins to quickly outperform the direct
methods for even moderately-sized problems, while scaling at a rate comparable
to the other indirect methods.

#### Discussion.

The solve time speedup of GeNIOS’s QPSolver with custom operators over the
equivalent formulation is likely, at least in part, a result of how the data
is stored. Using custom operators permits GeNIOS to store and use $F$ as a
dense matrix instead of storing and using $F$ as part of the sparse constraint
matrix. (This phenomenon is similar to what we observed in the Huber fitting
example in §4.3.) The ability to create and use these custom operators in the
original QP highlights the benefits of leveraging Julia’s multiple dispatch to
define optimization problems. Because GeNIOS is written in pure Julia, it does
not need the problem data to be a matrix; it only needs linear operators.
Often, it is much easier for the user to identify a fast way of applying $P$,
$M$, and $M^{T}$ in the original formulation (17) than to identify a good
problem reformulation. Here, for example, we easily defined a fast multiply
with the diagonal-plus-low-rank matrix $\Sigma$ and the original constraint
matrix $M=\begin{bmatrix}\mathbf{1}&I\end{bmatrix}^{T}$. Furthermore, the
Julia ecosystem has many packages that assist with constructing these fast
operators, including LinearMaps.jl [KHc23]. The speedup from GeNIOS’s
GenericSolver further emphasizes the advantage of allowing the user to
naturally specify problem structure to exploit. In this example, we recognized
that projection onto the simplex is fast, and GeNIOS allows us to handle this
part of the problem directly.

###  4.6 Signal decomposition

The signal decomposition problem demonstrates how GeNIOS’s GenericSolver
interface is flexible enough to handle nonconvex problems. Of course, none of
the convergence guarantees apply in this setting; however, ADMM applied to
nonconvex problems has been shown to sometimes work well in practice (see
[DTB18] and references therein). In the signal decomposition problem [MB+23],
a time series $y_{1},\dots,y_{T}\in\mathbf{R}$ is modeled as a sum of
components $x^{1},\dots,x^{K}\in\mathbf{R}^{T}$. Each of these $K$ components
is associated with a loss function
$\phi_{k}:\mathbf{R}^{T}\to\mathbf{R}\cup\\{\infty\\}$, which denotes the
implausibility of the current guess of $x^{k}$. The function $\phi_{k}$ can
take on the value $\infty$ to encode constraints; a signal $x$ is feasible if
$\phi_{k}(x)<\infty$. Given an observed signal $y$ and component classes
$1,\dots,K$ associated with loss functions $\phi_{1},\dots,\phi_{K}$, the
signal decomposition problem is

minimize
$\displaystyle(1/T)\|y-x^{2}-\cdots-x^{K}\|_{2}^{2}+\phi_{2}(x^{2})+\cdots+\phi_{K}(x^{K}),$

with variables $x^{2},\dots,x^{K}$. The first component $x^{1}$ is assumed to
be mean-square small.

In GeNIOS’s framework, this problem can be phrased as

minimize
$\displaystyle(1/T)\|y-x^{2}-\cdots-x^{K}\|_{2}^{2}+\phi_{2}(z^{2})+\cdots+\phi_{K}(z^{K})$
subject to $\displaystyle x-z=0,$

where $x=(x_{2},\dots,x_{K})$ and $z=(z_{2},\dots,z_{K})$. The function $f(x)$
is a quadratic, and the function $g(z)$ is separable across the $K$
components, so the proximal operators for $\phi_{2},\dots,\phi_{K}$ can be
evaluated in parallel.

(a) Observed signal

(b) Reconstructed signal

(c) Individual components of the reconstructed signal.

Figure 8: GeNIOS accurately decomposes the observed signal into its
components.

#### Problem data.

We generate synthetic data similar to the example in [MB+23, §2.9] with
$T=500$, and $K=3$ components. The first component is Gaussian noise. The
second component is a sine wave. The last component is a square wave. We
observe the sum of these components (see figure 8(a)). We choose three
component classes for the signal decomposition problem: mean-square small,
mean-square second-order smooth, and a signal constrained to only take on
values $\pm\theta$:

$\displaystyle\phi_{1}(x)$ $\displaystyle=\frac{1}{T}\|x\|_{2}^{2},$
$\displaystyle\phi_{2}(x)$
$\displaystyle=\frac{1}{T-2}\sum_{t=2}^{T-1}(x_{t+1}-2x_{t}+x_{t-1})^{2},$
$\displaystyle\phi_{3}(x)$ $\displaystyle=I_{\\{\pm\theta\\}^{T}}(x).$

These classes select for a small signal, a smoothly changing signal, and a
binary signal with known amplitude respectively. Here, we assume the amplitude
of the binary signal is known, but in a real example, we would solve the
problem for varying values of the amplitude hyperparameter.

#### Numerical Results.

The proximal operator for $\phi_{2}$ requires the solution of an unconstrained
convex quadratic program. Thus, this operator can be evaluated with a simple
linear system solve, and a banded matrix can be used to take advantage of
structure. The proximal operator for $\phi_{3}$ requires the projection onto
the nonconvex set $\\{\pm\theta\\}^{T}$. This projection can be computed very
quickly, but the nonconvexity of $\phi_{3}$ makes this optimization problem
nonconvex. As a result, GeNIOS is not guaranteed to converge, but it
empirically still works well for this problem. Figure 8 shows that GeNIOS
recovers the true components from the observed signal. This example highlights
the flexibility afforded by GeNIOS’s GenericSolver, which can be used to solve
problems with nonconvex constraint sets.

##  5 Conclusion

This paper introduced the new inexact ADMM solver GeNIOS, implemented in the
Julia language. This solver approximates the ADMM subproblems and solves these
approximate subproblems inexactly. These approximations yield several
benefits. First, the $x$-subproblem becomes a linear system solve, which
GeNIOS solves efficiently even for large problem sizes with the (Nyström
preconditioned) conjugate gradient method. Second, GeNIOS avoids conic
reformulations by handling the objective function directly, reducing the
problem size and improving memory locality. Despite the approximations and
inexact solves, GeNIOS retains the convergence rate of classic ADMM and can
detect infeasibility or unboundedness. Moreover, GeNIOS offers a flexible
interface, allowing the user to specify the objective function with zeroth,
first, and second order oracles. It can also work from just the zeroth order
oracle, computing higher order derivatives with automatic differentiation.

Through examples, we demonstrate that GeNIOS’s algorithmic improvements yield
substantial speedups over classic ADMM and over existing solvers for large
problem sizes. These improvements allow the user to exploit problem structure
using GeNIOS’s GenericSolver interface. GeNIOS also includes a specialized
QPSolver interface for quadratic programs and a MlSolver interface for machine
learning problems. Finally, we show that GeNIOS’s flexible interface allows
the user to specify and solve nonconvex problems as well.

There is more room to speed up GeNIOS. Parallelization of the dominant
operation—matrix-vector multiplies—using GPUs presents an obvious, likely
significant speedup. Future work on the algorithm itself should investigate
speedups from approximations or inexact solutions to the $z$-subproblem.
Additionally, the Nyström preconditioner in §2.2 performs best when the linear
system matrix has the form of low-rank-plus-identity (10). Developing
optimization problem modeling tools to automatically compile problems into a
form with this structure presents another interesting avenue for future work.
Finally, we only use standard normal test matrices to construct the
preconditioner used in our experiments. Other test matrices, such as
subsampled trigonometric transforms, may perform better for sparse
optimization problems.

## Acknowledgements

The authors thank Chris Rackauckas, Gaurav Arya, Axel Feldmann, and Guillermo
Angeris for helpful discussions. T. Diamandis is supported by the Department
of Defense (DoD) through the National Defense Science & Engineering Graduate
(NDSEG) Fellowship Program. B. Stellato is supported by the NSF CAREER Award
ECCS-2239771. Z. Frangella, S. Zhao, and M. Udell gratefully acknowledge
support from the National Science Foundation (NSF) Award IIS-2233762, the
Office of Naval Research (ONR) Award N000142212825 and N000142312203, and the
Alfred P. Sloan Foundation.

## References

  * [ADV+13] Martin S Andersen, Joachim Dahl and Lieven Vandenberghe  “CVXOPT: A Python package for convex optimization”  In _Available at cvxopt. org_ 54, 2013
  * [Agr+18] Akshay Agrawal, Robin Verschueren, Steven Diamond and Stephen Boyd  “A rewriting system for convex optimization problems”  In _Journal of Control and Decision_ 5.1, 2018, pp. 42–60
  * [AM15] Ahmed Alaoui and Michael W Mahoney  “Fast Randomized Kernel Ridge Regression with Statistical Guarantees”  In _Advances in Neural Information Processing Systems_ , 2015
  * [App+21] David Applegate, Mateo Díaz, Oliver Hinder, Haihao Lu, Miles Lubin, Brendan O’Donoghue and Warren Schudy  “Practical large-scale linear programming using primal-dual hybrid gradient”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 20243–20257
  * [App+21a] David Applegate, Mateo Díaz, Haihao Lu and Miles Lubin  “Infeasibility detection with primal-dual hybrid gradient for large-scale linear programming”  In _arXiv preprint arXiv:2102.04592_ , 2021
  * [App+22] David Applegate, Oliver Hinder, Haihao Lu and Miles Lubin  “Faster first-order primal-dual methods for linear programming using restarts and sharpness”  In _Mathematical Programming_ Springer, 2022, pp. 1–52
  * [ApS22] MOSEK ApS  “The MOSEK optimization toolbox for MATLAB manual. Version 10.0.”, 2022  URL: http://docs.mosek.com/9.0/toolbox/index.html
  * [Ban+19] Goran Banjac, Paul Goulart, Bartolomeo Stellato and Stephen Boyd  “Infeasibility detection in the alternating direction method of multipliers for convex optimization”  In _Journal of Optimization Theory and Applications_ 183 Springer, 2019, pp. 490–519
  * [BAS10] Lars Blackmore, Behçet Açikmeşe and Daniel P Scharf  “Minimum-landing-error powered-descent guidance for Mars landing using convex optimization”  In _Journal of guidance, control, and dynamics_ 33.4, 2010, pp. 1161–1171
  * [Bez+17] Jeff Bezanson, Alan Edelman, Stefan Karpinski and Viral B Shah  “Julia: A fresh approach to numerical computing”  In _SIAM review_ 59.1 SIAM, 2017, pp. 65–98
  * [BM11] T. Bertin-Mahieux  “YearPredictionMSD”, UCI Machine Learning Repository, 2011  URL: {DOI}:https://doi.org/10.24432/C50K61
  * [Boy+11] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato and Jonathan Eckstein  “Distributed optimization and statistical learning via the alternating direction method of multipliers”  In _Foundations and Trends® in Machine learning_ 3.1 Now Publishers, Inc., 2011, pp. 1–122
  * [Boy+17] Stephen Boyd, Enzo Busseti, Steve Diamond, Ronald N. Kahn, Kwangmoo Koh, Peter Nystrup and Jan Speth  “Multi-Period Trading via Convex Optimization”  In _Foundations and Trends® in Optimization_ 3.1, 2017, pp. 1–76  DOI: 10.1561/2400000023
  * [BPT12] Grigoriy Blekherman, Pablo A Parrilo and Rekha R Thomas  “Semidefinite optimization and convex algebraic geometry”  SIAM, 2012
  * [BT95] Paul T Boggs and Jon W Tolle  “Sequential quadratic programming”  In _Acta numerica_ 4 Cambridge University Press, 1995, pp. 1–51
  * [BTEGN09] Aharon Ben-Tal, Laurent El Ghaoui and Arkadi Nemirovski  “Robust optimization”  Princeton university press, 2009
  * [BTN98] Aharon Ben-Tal and Arkadi Nemirovski  “Robust convex optimization”  In _Mathematics of operations research_ 23.4 INFORMS, 1998, pp. 769–805
  * [BTN99] Aharon Ben-Tal and Arkadi Nemirovski  “Robust solutions of uncertain linear programs”  In _Operations research letters_ 25.1 Elsevier, 1999, pp. 1–13
  * [BV04] Stephen Boyd and Lieven Vandenberghe  “Convex Optimization”  Cambridge, United Kingdom: Cambridge University Press, 2004
  * [CG23] Yuwen Chen and Paul Goulart  “An Efficient IPM Implementation for A Class of Nonsymmetric Cones”  In _arXiv preprint arXiv:2305.12275_ , 2023
  * [CKV22] Chris Coey, Lea Kapelevich and Juan Pablo Vielma  “Solving natural conic formulations with Hypatia.jl”  In _INFORMS Journal on Computing_ 34.5 INFORMS, 2022, pp. 2686–2699  DOI: https://doi.org/10.1287/ijoc.2022.1202
  * [CL11] Chih-Chung Chang and Chih-Jen Lin  “LIBSVM: a library for support vector machines”  In _ACM transactions on intelligent systems and technology (TIST)_ 2.3 Acm New York, NY, USA, 2011, pp. 1–27
  * [CP11] Antonin Chambolle and Thomas Pock  “A first-order primal-dual algorithm for convex problems with applications to imaging”  In _Journal of mathematical imaging and vision_ 40 Springer, 2011, pp. 120–145
  * [CP11a] Patrick L Combettes and Jean-Christophe Pesquet  “Proximal splitting methods in signal processing”  In _Fixed-point algorithms for inverse problems in science and engineering_ Springer, 2011, pp. 185–212
  * [Dan51] George B Dantzig  “Application of the simplex method to a transportation problem”  In _Activity analysis and production and allocation_ Wiley, 1951
  * [DB16] Steven Diamond and Stephen Boyd  “CVXPY: A Python-embedded modeling language for convex optimization”  In _Journal of Machine Learning Research_ 17.83, 2016, pp. 1–5
  * [DB16a] Steven Diamond and Stephen Boyd  “Matrix-free convex optimization modeling”  In _Optimization and Its Applications in Control and Data Sciences: In Honor of Boris T. Polyak’s 80th Birthday_ Springer, 2016, pp. 221–264
  * [DCB13] Alexander Domahidi, Eric Chu and Stephen Boyd  “ECOS: An SOCP solver for embedded systems”  In _2013 European control conference (ECC)_ , 2013, pp. 3071–3076  IEEE
  * [DHL17] Iain Dunning, Joey Huchette and Miles Lubin  “JuMP: A modeling language for mathematical optimization”  In _SIAM review_ 59.2 SIAM, 2017, pp. 295–320
  * [DR56] Jim Douglas and Henry H Rachford  “On the numerical solution of heat conduction problems in two and three space variables”  In _Transactions of the American mathematical Society_ 82.2, 1956, pp. 421–439
  * [DTB18] Steven Diamond, Reza Takapoui and Stephen Boyd  “A general system for heuristic minimization of convex functions over non-convex sets”  In _Optimization Methods and Software_ 33.1 Taylor & Francis, 2018, pp. 165–193
  * [DY16] Wei Deng and Wotao Yin  “On the global and linear convergence of the generalized alternating direction method of multipliers”  In _Journal of Scientific Computing_ 66.3, 2016, pp. 889–916
  * [EB92] Jonathan Eckstein and Dimitri P Bertsekas  “On the Douglas—Rachford splitting method and the proximal point algorithm for maximal monotone operators”  In _Mathematical Programming_ 55.1 Springer, 1992, pp. 293–318
  * [Eck94] Jonathan Eckstein  “Parallel alternating direction multiplier decomposition of convex programs”  In _Journal of Optimization Theory and Applications_ 80.1 Springer, 1994, pp. 39–62
  * [EF98] Jonathan Eckstein and Michael C Ferris  “Operator-splitting methods for monotone affine variational inequalities, with a parallel application to optimal control”  In _INFORMS Journal on Computing_ 10.2 INFORMS, 1998, pp. 218–235
  * [FB18] Christopher Fougner and Stephen Boyd  “Parameter selection and preconditioning for a graph form solver”  In _Emerging Applications of Control and Systems Theory: A Festschrift in Honor of Mathukumalli Vidyasagar_ Springer, 2018, pp. 41–61
  * [Fra+23] Zachary Frangella, Shipu Zhao, Theo Diamandis, Bartolomeo Stellato and Madeleine Udell  “On the (linear) convergence of Generalized Newton Inexact ADMM”  In _arXiv preprint arXiv:2302.03863_ , 2023
  * [Fri23] Henrik A Friberg  “Projection onto the exponential cone: a univariate root-finding problem”  In _Optimization Methods and Software_ Taylor & Francis, 2023, pp. 1–17
  * [FTU21] Zachary Frangella, Joel A Tropp and Madeleine Udell  “Randomized Nyström Preconditioning”  In _arXiv preprint arXiv:2110.02820_ , 2021
  * [Gab83] D Gabay  “Applications of the method of multipliers to variational inequalities, M. Fortin and R. Glowinski (eds.), Augmented Lagrangian methods: Applications to the solution of boundary value problems”  North-Holland, Amsterdam, 1983
  * [GB14] Michael Grant and Stephen Boyd  “CVX: Matlab software for disciplined convex programming, version 2.1”, 2014
  * [GCG21] Michael Garstka, Mark Cannon and Paul Goulart  “COSMO: A Conic Operator Splitting Method for Convex Conic Problems”  In _Journal of Optimization Theory and Applications_ 190.3 Springer, 2021, pp. 779–810  DOI: 10.1007/s10957-021-01896-x
  * [GM75] Roland Glowinski and Americo Marroco  “Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualité d’une classe de problèmes de Dirichlet non linéaires”  In _Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique_ 9.R2 EDP Sciences, 1975, pp. 41–76
  * [GM76] Daniel Gabay and Bertrand Mercier  “A dual algorithm for the solution of nonlinear variational problems via finite element approximation”  In _Computers & mathematics with applications_ 2.1 Elsevier, 1976, pp. 17–40
  * [Gon12] Jacek Gondzio  “Interior point methods 25 years later”  In _European Journal of Operational Research_ 218.3 Elsevier, 2012, pp. 587–601
  * [GPM89] Carlos E Garcia, David M Prett and Manfred Morari  “Model predictive control: Theory and practice—A survey”  In _Automatica_ 25.3 Elsevier, 1989, pp. 335–348
  * [Gur23] Gurobi Optimization, LLC  “Gurobi Optimizer Reference Manual”, 2023  URL: https://www.gurobi.com
  * [Hei+16] Felix Heide, Steven Diamond, Matthias Nießner, Jonathan Ragan-Kelley, Wolfgang Heidrich and Gordon Wetzstein  “Proximal: Efficient image optimization using proximal algorithms”  In _ACM Transactions on Graphics (TOG)_ 35.4 ACM New York, NY, USA, 2016, pp. 1–15
  * [HFD16] Boris Houska, Janick Frasch and Moritz Diehl  “An augmented Lagrangian based algorithm for distributed nonconvex optimization”  In _SIAM Journal on Optimization_ 26.2 SIAM, 2016, pp. 1101–1127
  * [HS+52] Magnus R Hestenes and Eduard Stiefel  “Methods of conjugate gradients for solving linear systems”  In _Journal of research of the National Bureau of Standards_ 49.6, 1952, pp. 409–436
  * [Hub64] Peter J Huber  “Robust Estimation of a Location Parameter”  In _The Annals of Mathematical Statistics_ JSTOR, 1964, pp. 73–101
  * [HY12] Bingsheng He and Xiaoming Yuan  “On the O(1/n) convergence rate of the Douglas–Rachford alternating direction method”  In _SIAM Journal on Numerical Analysis_ 50.2 SIAM, 2012, pp. 700–709
  * [HYW00] Bing-Sheng He, Hai Yang and SL Wang  “Alternating direction method with self-adaptive penalty parameters for monotone variational inequalities”  In _Journal of Optimization Theory and applications_ 106 Springer, 2000, pp. 337–356
  * [Kan48] Leonid V Kantorovich  “On a problem of Monge”  In _CR (Doklady) Acad. Sci. URSS (NS)_ 3, 1948, pp. 225–226
  * [Kar84] Narendra Karmarkar  “A new polynomial-time algorithm for linear programming”  In _Proceedings of the sixteenth annual ACM symposium on Theory of computing_ , 1984, pp. 302–311
  * [KHc23] Daniel Karrasch, Jutho Haegeman and contributors  “LinearMaps.jl: A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently.”, https://github.com/JuliaLinearAlgebra/LinearMaps.jl, 2023
  * [Kim+07] Seung-Jean Kim, Kwangmoo Koh, Michael Lustig, Stephen Boyd and Dimitry Gorinevsky  “An interior-point method for large-scale $\ell\\_1$-regularized least squares”  In _IEEE journal of selected topics in signal processing_ 1.4 IEEE, 2007, pp. 606–617
  * [KKB07] Kwangmoo Koh, Seung-Jean Kim and Stephen Boyd  “An interior-point method for large-scale l1-regularized logistic regression”  In _Journal of Machine Learning Research_ 8.Jul, 2007, pp. 1519–1555
  * [Lau09] Monique Laurent  “Sums of squares, moment matrices and optimization over polynomials”  In _Emerging applications of algebraic geometry_ Springer, 2009, pp. 157–270
  * [LM79] Pierre-Louis Lions and Bertrand Mercier  “Splitting algorithms for the sum of two nonlinear operators”  In _SIAM Journal on Numerical Analysis_ 16.6 SIAM, 1979, pp. 964–979
  * [LN89] Dong C Liu and Jorge Nocedal  “On the limited memory BFGS method for large scale optimization”  In _Mathematical programming_ 45.1-3 Springer, 1989, pp. 503–528
  * [Lob+98] Miguel Sousa Lobo, Lieven Vandenberghe, Stephen Boyd and Hervé Lebret  “Applications of second-order cone programming”  In _Linear algebra and its applications_ 284.1-3 Elsevier, 1998, pp. 193–228
  * [LS91] László Lovász and Alexander Schrijver  “Cones of matrices and set-functions and 0–1 optimization”  In _SIAM journal on optimization_ 1.2 SIAM, 1991, pp. 166–190
  * [Mar52] HM Markowitz  “Portfolio Selection, the journal of finance. 7 (1)”  In _N_ 1, 1952, pp. 71–91
  * [MB12] Jacob Mattingley and Stephen Boyd  “CVXGEN: A code generator for embedded convex optimization”  In _Optimization and Engineering_ 13 Springer, 2012, pp. 1–27
  * [MB+23] Bennet E Meyers and Stephen P Boyd  “Signal decomposition using masked proximal operators”  In _Foundations and Trends® in Signal Processing_ 17.1 Now Publishers, Inc., 2023, pp. 1–78
  * [Meh92] Sanjay Mehrotra  “On the implementation of a primal-dual interior point method”  In _SIAM Journal on optimization_ 2.4 SIAM, 1992, pp. 575–601
  * [MM00] Olvi L Mangasarian and David R. Musicant  “Robust linear and support vector regression”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 22.9 IEEE, 2000, pp. 950–955
  * [MM17] Cameron Musco and Christopher Musco  “Recursive sampling for the Nystrom method”  In _Advances in Neural Information Processing Systems_ 30, 2017
  * [MOc20] A. Montoison, D. Orban and contributors  “Krylov.jl: A Julia Basket of Hand-Picked Krylov Methods”, https://github.com/JuliaSmoothOptimizers/Krylov.jl, 2020  DOI: 10.5281/zenodo.822073
  * [MR18] P Mogensen and A Riseth  “Optim: A mathematical optimization package for Julia”  In _Journal of Open Source Software_ 3.24 Open Journals, 2018
  * [MS13] Renato DC Monteiro and Benar F Svaiter  “Iteration-complexity of block-decomposition algorithms and the alternating direction method of multipliers”  In _SIAM Journal on Optimization_ 23.1 SIAM, 2013, pp. 475–507
  * [MT20] Per-Gunnar Martinsson and Joel A Tropp  “Randomized numerical linear algebra: Foundations and algorithms”  In _Acta Numerica_ 29 Cambridge University Press, 2020, pp. 403–572
  * [NN92] Yu Nesterov and Arkadi Nemirovsky  “Conic formulation of a convex programming problem and duality”  In _Optimization Methods and Software_ 1.2 Taylor & Francis, 1992, pp. 95–115
  * [NN94] Yurii Nesterov and Arkadii Nemirovskii  “Interior-point polynomial algorithms in convex programming”  SIAM, 1994
  * [Noc80] Jorge Nocedal  “Updating quasi-Newton matrices with limited storage”  In _Mathematics of computation_ 35.151, 1980, pp. 773–782
  * [O’D+16] Brendan O’Donoghue, Eric Chu, Neal Parikh and Stephen Boyd  “Conic Optimization via Operator Splitting and Homogeneous Self-Dual Embedding”  In _Journal of Optimization Theory and Applications_ 169.3, 2016, pp. 1042–1068  URL: http://stanford.edu/~boyd/papers/scs.html
  * [O’D21] Brendan O’Donoghue  “Operator Splitting for a Homogeneous Embedding of the Linear Complementarity Problem”  In _SIAM Journal on Optimization_ 31, 2021, pp. 1999–2023
  * [Ouy+15] Yuyuan Ouyang, Yunmei Chen, Guanghui Lan and Eduardo Pasiliao  “An accelerated linearized alternating direction method of multipliers”  In _SIAM Journal on Imaging Sciences_ 8.1, 2015, pp. 644–681
  * [PB+14] Neal Parikh and Stephen Boyd  “Proximal algorithms”  In _Foundations and trends® in Optimization_ 1.3 Now Publishers, Inc., 2014, pp. 127–239
  * [PE10] Daniel P Palomar and Yonina C Eldar  “Convex optimization in signal processing and communications”  Cambridge university press, 2010
  * [PL03] Pablo A Parrilo and Sanjay Lall  “Semidefinite programming relaxations and algebraic optimization in control”  In _European Journal of Control_ 9.2-3 Elsevier, 2003, pp. 307–321
  * [Raw00] James B Rawlings  “Tutorial overview of model predictive control”  In _IEEE control systems magazine_ 20.3 IEEE, 2000, pp. 38–52
  * [RGN22] Nikitas Rontsis, Paul Goulart and Yuji Nakatsukasa  “Efficient semidefinite programming with approximate admm”  In _Journal of Optimization Theory and Applications_ Springer, 2022, pp. 1–29
  * [RLP16] J. Revels, M. Lubin and T. Papamarkou  “Forward-Mode Automatic Differentiation in Julia”  In _arXiv:1607.07892 [cs.MS]_ , 2016  URL: https://arxiv.org/abs/1607.07892
  * [Roc76] R Tyrrell Rockafellar  “Monotone operators and the proximal point algorithm”  In _SIAM journal on control and optimization_ 14.5 SIAM, 1976, pp. 877–898
  * [RR07] Ali Rahimi and Benjamin Recht  “Random features for large-scale kernel machines”  In _Advances in neural information processing systems_ 20, 2007
  * [RR08] Ali Rahimi and Benjamin Recht  “Uniform approximation of functions with random bases”  In _2008 46th annual allerton conference on communication, control, and computing_ , 2008, pp. 555–561  IEEE
  * [RY22] Ernest K Ryu and Wotao Yin  “Large-Scale Convex Optimization: Algorithms & Analyses via Monotone Operators”  Cambridge University Press, 2022
  * [SBL20] Michel Schubiger, Goran Banjac and John Lygeros  “GPU acceleration of ADMM for large-scale quadratic programming”  In _Journal of Parallel and Distributed Computing_ 144 Elsevier, 2020, pp. 55–67
  * [SGV22] Mario Souto, Joaquim D Garcia and Álvaro Veiga  “Exploiting low-rank structure in semidefinite programming by approximate operator splitting”  In _Optimization_ 71.1 Taylor & Francis, 2022, pp. 117–144
  * [SNW12] Suvrit Sra, Sebastian Nowozin and Stephen J Wright  “Optimization for machine learning”  Mit Press, 2012
  * [Ste+20] B. Stellato, G. Banjac, P. Goulart, A. Bemporad and S. Boyd  “OSQP: an operator splitting solver for quadratic programs”  In _Mathematical Programming Computation_ 12.4, 2020, pp. 637–672  DOI: 10.1007/s12532-020-00179-2
  * [Ste+22] Lorenzo Stella, Niccolò Antonello, Mattias Fält, Dmitry Volodin, Domagoj Herceg, Elliot Saba, Fredrik Bagge Carlson, Tony Kelman, sp, Ellis Brown, Julia TagBot and Pantelis Sopasakis  “JuliaFirstOrder/ProximalOperators.jl: v0.15.3”  Zenodo, 2022  DOI: 10.5281/zenodo.6583304
  * [TBI97] Lloyd N Trefethen and David Bau III  “Numerical Linear Algebra”  SIAM, 1997
  * [Tro+17] Joel A Tropp, Alp Yurtsever, Madeleine Udell and Volkan Cevher  “Fixed-rank approximation of a positive-semidefinite matrix from streaming data”  In _Advances in Neural Information Processing Systems_ 30, 2017
  * [Tro+19] Joel A Tropp, Alp Yurtsever, Madeleine Udell and Volkan Cevher  “Streaming low-rank matrix approximation with an application to scientific simulation”  In _SIAM Journal on Scientific Computing_ 41.4 SIAM, 2019, pp. A2430–A2463
  * [Ude+14] Madeleine Udell, Karanveer Mohan, David Zeng, Jenny Hong, Steven Diamond and Stephen Boyd  “Convex optimization in Julia”  In _2014 First Workshop for High Performance Technical Computing in Dynamic Languages_ , 2014, pp. 18–28  IEEE
  * [UT19] Madeleine Udell and Alex Townsend  “Why are big data matrices approximately low rank?”  In _SIAM Journal on Mathematics of Data Science_ 1.1 SIAM, 2019, pp. 144–160
  * [VBEG97] Lieven Vandenberghe, Stephen Boyd and Abbas El Gamal  “Optimal wire and transistor sizing for circuits with non-tree topology”  In _1997 Proceedings of IEEE International Conference on Computer Aided Design (ICCAD)_ , 1997, pp. 252–259  IEEE
  * [VBEG98] Lieven Vandenberghe, Stephen Boyd and Abbas El Gamal  “Optimizing dominant time constant in RC circuits”  In _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ 17.2 IEEE, 1998, pp. 110–125
  * [Wil63] Robert B Wilson  “A simplicial algorithm for concave programming”  In _Ph. D. Dissertation, Graduate School of Bussiness Administration_ Harvard University, 1963
  * [WL01] SL Wang and LZ Liao  “Decomposition method with a variable parameter for a class of monotone variational inequality problems”  In _Journal of optimization theory and applications_ 109 Springer, 2001, pp. 415–429
  * [WWK15] Matt Wytock, Po-Wei Wang and J Zico Kolter  “Convex programming with fast proximal and linear operators”  In _arXiv preprint arXiv:1511.04815_ , 2015
  * [XHY96] Xiaojie Xu, Pi-Fang Hung and Yinyu Ye  “A simplified homogeneous and self-dual linear programming algorithm and its implementation”  In _Annals of Operations Research_ 62.1 Springer, 1996, pp. 151–171
  * [YTM94] Yinyu Ye, Michael J Todd and Shinji Mizuno  “An $O(\sqrt{n}L)$-iteration homogeneous and self-dual linear programming algorithm”  In _Mathematics of operations research_ 19.1 INFORMS, 1994, pp. 53–67
  * [ZFU22] Shipu Zhao, Zachary Frangella and Madeleine Udell  “NysADMM: faster composite convex optimization via low-rank approximation”  In _Proceedings of the 39th International Conference on Machine Learning_ 162, Proceedings of Machine Learning Research  PMLR, 2022  URL: https://proceedings.mlr.press/v162/zhao22a.html
  * [ZLU21] Shipu Zhao, Laurent Lessard and Madeleine Udell  “An automatic system to detect equivalence between iterative algorithms”  In _arXiv preprint arXiv:2105.04684_ , 2021

##  Appendix A Duality gap bounds for machine learning problems

In this section, we derive bounds for the duality gap for the machine learning
examples we consider in §3.3. The derivation is largely inspired by [Kim+07,
KKB07]. We consider problems of the form

minimize
$\displaystyle\ell(x)=\sum_{i=1}^{m}f(a_{i}^{T}x-b_{i})+\lambda_{1}\|x\|_{1}+(1/2)\lambda_{2}\|x\|_{2}^{2},$
(19)

where $f(\cdot)$ is some loss function, which we assume to be convex and
differentiable and $\lambda_{1},\;\lambda_{2}\geq 0$ are coefficients of the
regularization terms which are assumed to not both be $0$. (We do not need
differentiability of $f$, but is convenient for our purposes.) Tthe optimality
conditions of (19) are that $0$ is contained in the subdifferential, i.e.,

$0\in\sum_{i=1}^{m}f^{\prime}(a_{i}^{T}x-b_{i})\cdot
a_{i}+\lambda_{1}\partial\|x\|_{1}+\lambda_{2}x.$

By rearranging, we have the condition that

$\left(\sum_{i=1}^{m}f^{\prime}(a_{i}^{T}x-b_{i})\cdot
a_{i}\right)_{i}+\lambda_{2}x_{i}\in\begin{cases}\\{+\lambda_{1}\\}&x_{i}<0\\\
\\{-\lambda_{1}\\}&x_{i}>0\\\ [-\lambda_{1},\lambda_{1}]&x_{i}=0\end{cases}$
(20)

for $i=1,\dots,n$. These optimality conditions indicate that $x=0$ is a
solution to (19) if and only if
$\|A^{T}f^{\prime}(-b)+\lambda_{2}x\|_{\infty}\leq\lambda_{1}$, where the
function $f^{\prime}$ is applied elementwise.

#### Lagrangian and primal-dual optimality.

We introduce new variable $y_{i}$ and reformulate the primal problem (19) as

minimize $\displaystyle
l(x)=\sum_{i=1}^{m}f(y_{i})+\lambda_{1}\|x\|_{1}+(1/2)\lambda_{2}\|x\|_{2}^{2}$
(21) subject to $\displaystyle y=Ax-b.$

The Lagrangian is then

$L(x,y,\nu)=\sum_{i=1}^{m}f(y_{i})+\lambda_{1}\|x\|_{1}+(1/2)\lambda_{2}\|x\|_{2}^{2}+\nu^{T}(Ax-
b-y).$ (22)

A primal dual point $(x,y,\nu)$ is optimal when it is feasible, i.e., ,
$y=Ax-b$ and the gradient of the Lagrangian vanishes:

$\nu_{i}=f^{\prime}(a_{i}^{T}x-b_{i})\qquad\text{and}\qquad(A^{T}\nu+\lambda_{2}x)_{i}\in{\begin{cases}\\{+\lambda_{1}\\}&x_{i}<0\\\
\\{-\lambda_{1}\\}&x_{i}>0\\\ [-\lambda_{1},\lambda_{1}]&x_{i}=0.\end{cases}}$
(23)

In deriving the dual function, we will have to consider the case when
$\lambda_{2}=0$, in which case we are solving an $\ell_{1}$-regularized
problem, separately from the case that $\lambda_{2}>0$, in which case the
regularization term is smooth.

#### Optimality gap.

First, we will consider the case when $\lambda_{2}>0$. Partially minimizing
the Lagrangian over $x$ and $y$, we get the dual function

$g(\nu)=-\sum_{i=1}^{m}f^{*}(\nu_{i})-b^{T}\nu-(1/2\lambda_{2})\sum_{i=1}^{n}\left(\left(\lvert{(A^{T}\nu)_{i}}\rvert-\lambda_{1}\right)_{+}\right)^{2},$

where $f^{*}$ is the Fenchel conjugate of $f$ [BV04, §3.3]. If, instead, we
have that $\lambda_{2}=0$, then

$g(v)=\begin{cases}-\sum_{i=1}^{m}f^{*}(\nu_{i})-b^{T}\nu,&\|A^{T}\nu\|_{\infty}\leq\lambda_{1}\\\
-\infty,&\text{otherwise}.\end{cases}$

The dual problem is simply

maximize $\displaystyle g(\nu),$ (24)

where the norm ball constraint is encoded in the objective. Importantly, for
any dual feasible point $\nu$, we have that

$\ell(x)\geq\ell(x^{\star})=g(\nu^{\star})\geq g(\nu),$

where $x^{\star}$ and $\nu^{\star}$ are optimal solutions to (19) and (24)
respectively. It immediately follows that

$\frac{\ell(x)-g(\nu)}{\min\left(\ell(x),\lvert{g(\nu)}\rvert\right)}\geq\frac{\ell(x)-\ell(x^{\star})}{\min\left(\ell(x),\lvert{g(\nu)}\rvert\right)}.$
(25)

We will call the left hand side the _relative duality gap_ , which clearly
bounds the relative suboptimality of a primal-dual feasible point $(x,\nu)$.

#### Dual feasible points.

Now, we must devise a way to construct dual feasible points $\nu$. Inspired by
the optimality conditions (23), we construct a dual feasible point by taking
$\nu_{i}=f^{\prime}(a_{i}^{T}x-b_{i})$ and then, when $\lambda_{2}=0$,
projecting onto the norm ball given by the first optimality condition.
Specifically, in this case, we take

$\nu_{i}=\frac{\lambda_{1}}{\left\|\sum_{i=1}^{m}f^{\prime}(a_{i}^{T}x-b_{i})\cdot
a_{i}+\lambda_{2}x\right\|_{\infty}}\cdot f^{\prime}(a_{i}^{T}x-b_{i}).$

When $x\neq x^{\star}$, this projection ensures that
$\|A^{T}\nu\|_{\infty}\leq\lambda_{1}$, and, therefore, the function $g(\nu)$
is finite-valued if the conjugate function is finite valued. When
$x=x^{\star}$, this dual variable will be optimal, as it satisfies the
optimality conditions by construction; the original optimality condition (20)
ensures that $\nu$ will be unaffected by the projection and, therefore,
satisfy both conditions in (23).

This construction, along with the bound (25) suggests a natural stopping
criterion. For any primal feasible point, we construct a dual feasible point
$\nu$. Using the dual feasible point, we evaluate the duality gap. We then
terminate the solver if the relative gap is less than some tolerance
$\varepsilon$ as

$\frac{l(x)-g(\nu)}{\max\left(\lvert{g(\nu)\rvert},l(x)\right)}\leq\varepsilon.$

If this condition is met, we are guaranteed that the true relative error is
also less than $\varepsilon$ from (25).

##  Appendix B Logistic regression conic form

Here, we consider a modified version of (17),

minimize $\displaystyle(1/2)x^{T}Px+q^{T}x$ (26) subject to $\displaystyle
Mx\in C$

where $C=C_{1}\times\cdots\times C_{d}$ and each $C_{i}$ is either a box or an
exponential cone, defined as

$K_{\mathrm{exp}}=\left\\{(x,y,z)\mid y>0\text{ and }y\exp(x/y)\leq
z\right\\}.$

(We work with this form because it requires only a slight modification of our
QPSolver interface (17).) The logistic regression problem is

minimize
$\displaystyle\sum_{i=1}^{N}\log\left(1+\exp(a_{i}^{T}x)\right)+\lambda_{1}\|x\|_{1},$

with variable $x\in^{n}$. We can reformulate this into an exponential cone
program by introducing new variables $p\in\mathbf{R}^{n}$,
$t\in\mathbf{R}^{N}$, $r\in\mathbf{R}^{N}$, $u\in\mathbf{R}^{N}$,
$s\in\mathbf{R}^{N}$, and $v\in\mathbf{R}^{N}$ via transformation into
epigraph form (see [BV04, §4.2.4]). We also notice that

$\log\left(1+\exp(a_{i}^{T}x)\right)\leq
t_{i}\iff\exp(-t_{i})+\exp(a_{i}^{T}x)\leq 1.$

We can transform convex exponentials with inequalities into exponential cone
constraints by using the fact

$\exp(a)\leq b\iff(a,1,b)\in K_{\mathrm{exp}}.$

Putting these transformations together, we have the equivalent problem

minimize $\displaystyle\sum_{i=1}^{N}t_{i}+\lambda_{1}p_{i}$ subject to
$\displaystyle-q\leq x\leq q$ $\displaystyle u+v\leq 1$
$\displaystyle(-t_{i},r_{i},u_{i})\in K_{\mathrm{exp}},\qquad i=1,\dots N$
$\displaystyle(a_{i}^{T}x-t_{i},s_{i},v_{i})\in K_{\mathrm{exp}},\qquad
i=1,\dots N$ $\displaystyle r=1$ $\displaystyle s=1.$

After completing the transformation and putting it into the conic form (26),
the matrix $P=0$, $q\in\mathbf{R}^{2n+5N}$, and $M\in\mathbf{R}^{2n+9N\times
2n+5N}$. Although the side dimensions of $M$ are large, this matrix is very
sparse and highly structured.

##  Appendix C Projection onto the simplex

In this section, we derive an efficient method to project onto the simplex,
i.e., to solve the optimization problem

minimize $\displaystyle(1/2)\|\tilde{z}-v\|_{2}^{2}$ subject to
$\displaystyle\mathbf{1}^{T}\tilde{z}=1$ $\displaystyle\tilde{z}\geq 0.$

Our approach is similar to the approach used to solve the water-filling
problem in communications. The Lagrangian for this problem is

$\mathcal{L}({\tilde{z}},\nu,\lambda)=(1/2)\|{\tilde{z}}-v\|_{2}^{2}+\nu(\mathbf{1}^{T}{\tilde{z}}-1)-\lambda^{T}{\tilde{z}}.$

Denote the optimal primal and dual variables by ${\tilde{z}}^{\star}$,
$\nu^{\star}$, and $\lambda^{\star}$. The optimality conditions are that
${\tilde{z}}^{\star}$ is primal feasible, $\lambda^{\star}\geq 0$,
complementary slackness,

${\tilde{z}}^{\star}_{i}>0\implies\lambda^{\star}_{i}=0\qquad\text{and}\qquad\lambda_{i}^{\star}>0\implies{\tilde{z}}_{i}^{\star}=0,$

and that the gradient of the Lagrangian vanishes,

${\tilde{z}}^{\star}-v+\nu^{\star}\mathbf{1}-\lambda^{\star}=0.$

From these conditions, we can conclude that

${\tilde{z}}^{\star}_{i}>0\implies\tilde{z}_{i}^{\star}=v_{i}-\nu^{\star}.$

Otherwise, we have that

$v_{i}-\nu^{\star}+\lambda^{\star}_{i}=0\implies v_{i}-\nu^{\star}\leq 0.$

Thus, we can conclude that

$\tilde{z}^{\star}=(v-\nu^{\star}\mathbf{1})_{+}.$

Plugging this expression for $\tilde{z}^{\star}$ into the constraint that the
sum of the entries is equal to $1$, we have that $\nu^{\star}$ must satisfy

$\sum_{i=1}^{n}(v_{i}-\nu^{\star})_{+}=1.$

We can solve this single variable equation with bisection search for
$\nu^{\star}$, and then solve for $\tilde{z}^{\star}$. Easy upper and lower
bounds for $\nu^{\star}$ are $\max_{i}\lvert{v_{i}}\rvert$ and
$-\max_{i}\lvert{v_{i}}\rvert-1$ respectively.