Epigraph projections for fast general convex programming

Po-Wei Wang
Machine Learning Department, Carnegie Mellon University, Pittsburgh, PA 15213 USA
Matt Wytock
Electrical Engineering Department, Stanford University, Stanford, CA 94305 USA
J. Zico Kolter
School of Computer Science, Carnegie Mellon University, Pittsburgh, PA 15213 USA

Abstract
This paper develops an approach for efficiently
solving general convex optimization problems
specified as disciplined convex programs (DCP),
a common general-purpose modeling framework. Specifically we develop an algorithm
based upon fast epigraph projections, projections
onto the epigraph of a convex function, an approach closely linked to proximal operator methods. We show that by using these operators, we
can solve any disciplined convex program without transforming the problem to a standard cone
form, as is done by current DCP libraries. We
then develop a large library of efficient epigraphprojection operators, mirroring and extending
work on fast proximal algorithms, for many common convex functions. Finally, we evaluate theperformance of the algorithm, and show it often
achieves order of magnitude speedups over existing general-purpose optimization solvers.

1. Introduction
Although convex optimization techniques underly a large
number of machine learning algorithms, there has been a
traditional tension between general purpose optimization
methods and specialized algorithms. General purpose algorithms, exemplified by standard cone form solvers likelinear, second order, and semidefinite cone solvers, withthe addition of modeling languages such as cvx (Grant &Boyd, 2008) or cvxpy (Diamond & Boyd, 2015) (whichconvert problems to these forms), provide a very flexible “rapid prototyping” framework for convex optimizaProceedings of the 33 rd International Conference on Machine
Learning, New York, NY, USA, 2016. JMLR: W&CP volume
48. Copyright 2016 by the author(s).

POWEIW @ CS . CMU . EDU

MWYTOCK @ STANFORD . EDUZKOLTER @ CS . CMU . EDU

tion. However, these methods typically do not scale beyond medium-size problems, and are not well-suited to thelarger problems that make up much current machine learning work. This has lead to the development of a number of specialized solvers for many problems of interest
including support vector machines (Platt, 1999; ShalevShwartz et al., 2011; Chang & Lin, 2011), ranking approaches (Joachims, 2002), sparse inverse covariance estimation (Friedman et al., 2008; Hsieh et al., 2014), etc,
just to name a very small number of examples. Recently,
many special purpose algorithms have started to use proximal operators as a key building block.
This current work looks to begin bridging the gap between
specialized and general purposes solvers. The first main
contribution of this work is a method for converting any
disciplined convex program (DCP) into a form that can bedirectly solved by proximal methods, without converting
to cone form. Doing so requires a new set of operatorsfor epigraph projections, an operator related to but distinct from most existing proximal operators; specifically,
for some convex function f : Rn → R, epigraph projections are optimization problems of the form
epif (v, s) = argmin {kx−vk22 +(t−s)2 | f (x) ≤ t} (1)x,ti.e., they are a projection of some (v ∈ Rn , s ∈ R) onto theepigraph set {(x, t) | f (x) ≤ t}. Thus, the second main
contribution of this work is the development of a wide set of
epigraph projection algorithms (plus some additional proximal operators), which in turn can solve a very broad rangeof DCPs without ever resorting to standard cone transformations. To build these fast epigraph projections operators we develop new optimization approaches, including
an implicit dual Newton method and a O(n) “sum-of-clip”
solver, both of which we detail in later sections. Finally, we
implement these methods in a generic optimization solver,
and show significant speedup over existing general-purpose
approaches, often an order of magnitude or more.

Epigraph projections for fast general convex programming

2. Background
Disciplined convex programming Disciplined convex
programing (DCP) (Grant et al., 2006) is a modeling framework for convex programs that allows very general convex
problems to be specified using a relatively small set of convex atomic function and set of compositional rules that preserve convexity. For instance, a basic rule states that for
h : Rk → R, gi : Rn → R, and one of the following holds
for each i = 1, . . . , k:
• gi convex, h convex nondecreasing in argument i• gi concave, h concave nonincreasing in argument i• gi affine, h convex,
then f = h(g1 (x), . . . , gk (x)), is convex (see e.g. Boyd &Vandenberghe (2004, §3.2.4) for details). The DCP ruleset
is sufficient but not necessary for a problem to be convex:
the “log-sum-exp” function is convex, but cannot be verified as such by the above rules. However, log-sum-exp canbe represented as a separate atomic function that is convex(and monotonic in its arguments). In practice, most convex problems can be written using the DCP ruleset with arelatively small set of atoms.
In addition to this verification, DCP libraries also provide ameans of converting the problems to a standard form cone
problem. Each functional atom also provides an graph implementation of the function, a representation of the function as the solution to a linear cone problem. For example,
the `1 norm has the graph form
kxk1 = min 1T y, subject to −y ≤ x ≤ y.
y

(2)Once an optimization problem has been verified as convex
by the DCP rules, we replace all instances of DCP atoms
with their corresponding epigraph implementation (introducing new variables as needed). The resulting problem,
which must be itself a linear cone problem, is guaranteed
to be equivalent to the original optimization problem, andcan then be solved by standard form solvers.
Proximal methods A recent trend in machine learning
optimization methods has been the increased developmentof algorithms based upon proximal operators. Given a convex function f : Rn → R, the proximal operator is definedas1
proxλf (v) = argmin λf (x) + kx − vk22 .
(3)2x
Crucially, for many functions f , including many nonsmooth functions, we can compute the proximal operator
in closed form (or if not, at least compute it to numericalprecision very efficiently). For instance, the proximal operator for the `1 norm is given by soft thresholding (Donoho,1995), proxλk·k1 (v) = max{|v| − λ, 0} · sign(v), where all
operators are applied elementwise.

Generally speaking, “proximal algorithms” refer to any optimization method that uses a proximal operator in its iteration. Such algorithms are not new, with the original proximal point algorithm proposed in 1976 (Rockafellar, 1976),
but they have seen increased usage in recent years, often
in conjunction with the increasing use of `1 regularization;
see e.g. (Parikh & Boyd, 2013) for recent survey of several
such methods.
Operator splitting techniques are one class of proximal algorithm, which solve a composite optimization problem
minimize f (x) + g(x),x(4)

typically by exploiting fast proximal operators for f and g.
A general review of operator splitting algorithms is givenin (Ryu & Boyd, 2016), and two algorithms of particular
recent interest are Douglas-Rachford splitting (Douglas &Rachford, 1956) and the alternating direction method ofmultipliers (ADMM) (Gabay & Mercier, 1976; Boyd et al.,2011). Ultimately, our algorithm uses ADMM to solve theresulting problems after transforming them to a suitable
form, though other operator splitting methods can be applied as well. Of particular relevance to our problem is the
splitting cone solver (SCS) (O’Donoghue et al., 2013), an
application of ADMM to linear cone programs; together
with existing DCP libraries, this system is very comparable to our own (a solver capable of solving DCP problems
using proximal methods), and comparisons to SCS will be
one of the primary focuses in later sections. Also related
is the TFOCS algorithm (Becker et al., 2011), which uses
first order methods to solve several classes of cone problems, though it does not solve general DCPs.
In the general theme of proximal algorithms, there havebeen a few papers that do consider epigraph projections
in a limited manner. For example, Tofighi et al. (2014)use a epigraph projection onto a total-variation-type norm
described via a linear program for image deconvolution;
Chierchia et al. (2012) describe projections on to the epigraph of the `2 and `∞ norms (in the later case using a naiveO(n log n) algorithm), for a few specialized optimization
problems; and Harizanov et al. (2013) use an epigraph projection onto a very specific function useful for image processing. None of these approaches relate to solving general
optimization problems, nor do they the develop the rangeof epigraph projections that we cover.

3. Converting DCPs to proximal form with
epigraph projections
In this section, we show that any DCP can be solved using a proximal algorithm that employs standard proximal
operators plus epigraph projections. In particular we showthat any DCP composed entirely of atoms from some setof functions G (but which can be composed according toEpigraph projections for fast general convex programming

the DCP ruleset to form much more complex functions that
do not admit efficient proximal operators or epigraph projections) can be solved using only proximal operators andepigraph projections for functions g ∈ G; thus, if we implement proximal and epigraph projection operators for every
element in a DCP atom library, we can directly solve theproblem without ever using any of the DCP graph implementations, but just the relevant operators directly. This isformalized as follows
Theorem 1. Consider the DCP optimization problem

+
×
λ

sum
∥ · ∥22

max{0, ·}θ+1−diag(y)X×∥ · ∥2×θPTθ

minimize f0 (x)xsubject to fi (x) ≤ 0, i = 1, . . . , p
Ax = b(5)

with optimization variable x, affine constraints A ∈ Rm×n
and b ∈ Rm , and convex functions fi : Rn → R built via
the DCP ruleset from convex atoms in a set G. Then thisoptimization problem can be written equivalently as
minimizex̄NXf¯i (x̄)

i=1

for some expanded set of variables x̄ ∈ Rn̄ such that x =x̄I0 for some set I0 ; affine constraints Ā ∈ Rm̄×n̄ andb̄ ∈ Rm̄ ; and convex functions f¯ : Rn̄ → R of the formg(x̄Ii )f¯i (x̄) =or(7)I{g(x̄Ii ) ≤ x̄ti }for some atomic function g ∈ G where Ii selects some subject of the indices and ti 6∈ I selects the epigraph variable.
Again, the key point here is that although atom g ∈ G admits efficient operators, the functions fi made up of compositions of these functions can be much more complex(one simple example for a robust SVM is below). But thetheorem states that the problem can be transformed into an
equivalent one that only requires proximal operators andepigraph projections for functions in the atom set G.
The proof of Therem 1 is straightforward, but requiresa slightly more formal definition of the representation of
DCP functions; we describe this only briefly here, withmore detail about DCPs in general available in Grant et al.
(2006). Each DCP function is represented as a expression
tree, where each non-leaf node in the tree must be a DCP
atom (which itself can be either convex, concave, or affine),and each leaf node must be a variable or a constant. For example, an `∞ robust SVM (see Appendix A.4 for details
about the derivation of this objective term) can be written
as the optimization problem
θ

Figure 1 shows the DCP expression tree for this objective:k · k22 , k · k1 , and max{0, ·} atoms are convex, the sum, +and × atoms are affine (the × atom under DCP must havea constant on one side), the leaf node θ is a variable, andthe λ, 1, −diag(y)X and P T leaf nodes are constants.

(6)

subject to Āx̄ = b̄minimize

Figure 1. A DCP expression tree representing an `∞ robust SVM
optimization objective.

mX
λkθk22 +
max{0, 1 − yi · θT xi + kP T θk1 }.
2
i=1
(8)

Proof of Theorem 1. The proof proceeds by considering
each DCP function fi as an expression tree of depth n,and employs a type of “bottom up” reduction to reduce
the tree to one of level n − 1, plus some additional equality and epigraph constraints. We select some leaf nodeat level n, which we denote l1 and consider its immediate parent g and siblings, l2 , . . . , ln . We introduce a newvariable ˆl and, based upon whether g is affine, convex, orconcave, we add the constraint g(l1 , . . . , ln ) {=, ≤, ≥} ˆl:we add an equality constraint for g affine and place these
constraints into the Āx̄ = b̄ matrix; and we add ≤ constraints for g convex and ≥ constraints for g concave,and introduce the appropriate epigraph indicator functionsf¯i = I{g(l1 , . . . , ln ) {≥, ≤} ˆl} in these two cases. Finally,
we replace the entire subtree with the ˆl variable. We repeat
this process for all leaves at level n, resulting in an equivalent expression tree of depth n − 1 (plus additional equality
constraints and epigraph indicators).
When we reach the root node of the tree for the objectivefunction f0 , we simply add the expression itself as a function f¯i = g (or by way of optimization, we immediately
terminate if the root node is an addition operator and thenodes at depth 2 already have efficient proximal operators
or epigraph projections). For the root nodes of expression trees for fi , i ≥ 1, we add the epigraph constraintf¯i = I(g(l1 , . . . , ln ) ≤ ˆl) plus the constraint ˆl = 0.
Example The theorem is best illustrated by example, so
we consider here how we can use this approach to transform the robust SVM problem (8). This objective composes the hinge loss, a linear function of the parameters θ,
and an `1 norm inside the hinge loss, and there is no prox-

Epigraph projections for fast general convex programming

imal operator that can be applied directly to the entire toplevel objective terms. However, the objective still admits a
DCP formulation, as shown in Figure 1, and applying theoperations from Theorem 1 (using atoms for the squared`2 norm, hinge function max{0, ·}, and `1 norm) results inthe new variables and constraints
ˆl1 ⇒ P T x = ˆl1
ˆl2 ⇒ kˆl1 k1 ≤ ˆl2(9)ˆl3 ⇒ 1 − diag(y)Xθ + 1ˆl2 = ˆl3
and the optimization problem over new variable x̄ =(x, ˆl1 , ˆl2 , ˆl3 )X
λkθk22 +
max{ˆl3 , 0} + I{kˆl1 k1 ≤ ˆl2 }2subject to P T x = ˆl1
1 − diag(y)Xθ + 1ˆl2 = ˆl3 .
minimizex̄(10)This is an optimization problem in the form (6): each termin the objective is either one of the atomic functions, or anindicator of the epigraph of an atomic function, plus additional linear constraints over these variables.
Optimization approaches Finally, as an immediate
corollary of the above theorem comes the fact that we cansolve the modified DCP problem using just a sequenceof proximal and epigraph projection operators of just theatomic functions. The proof here is algorithmic, and comesby simply applying an operator splitting algorithm like
ADMM directly to the modified problem (6). Indeed, thesolver we use in Section 5 uses an ADMM solver, thoughmany operator splitting approaches can solve problems of
the form (6) equally well. We omit the full derivation forbrevity (see e.g. Boyd et al. (2011) for a full derivation of
ADMM in this so-called consensus setting), but briefly, tosolve (6) we can introduce N copies of the x̄ variable, plus
a consensus variable z, and solve the optimization problem
minimize,x̄1 ,...,x̄N ,zNX

f¯i (x̄i ) + I{Āz = b̄}i=1

terms that are indicators of an epigraph set is exactly
proxf¯i (v) = epig (vIi , vti ).

Thus, the solution above requires computing only proximaland epigraph projection operators for the underlying atoms
in the DCP, as stated originally. In practice, we make a few
further refinements of this algorithm (the splitting consensus variable z can be created jointly with the constraints
Ax̄ = b, and some prox operators can directly incorporate
affine terms), but the above method is suffices to show the
key point: that we can solve any DCP problem built fromatoms for which we can implement efficient proximal operators and epigraph projections.

4. Fast epigraph projection solvers
While proximal operators for many functions have beenstudied extensively in the machine learning and optimization literature, there has been very little treatment of epigraph projections. Therefore, because we need these operators to solve general DCP problems using the method
described previously, in this section we develop a numberof methods for efficiently computing the epigraph projection for a wide set of convex atoms.
Consider the optimization problem for epigraph projection
1
k(x, t) − (v, s)k22
2
subject to f (x) ≤ t.
minimizex,twhich results in the update rules (assuming for simplicity
an augmented Lagrangian parameter ρ = 1)x̄k+1
← proxf¯i (uki − z k )i
−1  1 PN
I ĀT(x̄k+1+ uk )
k+1
i
i=1
N
(12)
z
←
Ā 0
b̄
uk+1
← uki + x̄k+1− z k+1 .
iiSince the proximal operator for the indicator of a set is simply the projection onto that set, the proximal operator for f¯i(14)

The Lagrangian of this problem is given byL(x, t, λ) =1k(x, t) − (v, s)k22 + λ(f (x) − t).
2(15)

Minimizing over x and t gives the solutions
x? = proxλf (v), t? = s + λ,(16)which correspond to the dual problem
1
maximize Pλf (v) − λ2 − λs,
λ≥0
2

(11)

subject to x̄i = z, i = 1, . . . , N

(13)

where

(17)

1
Pλf (v) = min λf (x) + kx − vk22
(18)
x
2
denotes the objective obtained by the proximal operator.
The dual problem is concave with a single variable λ, so
it can always be solved via bisection. However, the bisection algorithm takes O(log 1/) to iterations to reach accuracy , which can itself become a bottleneck. A main
contribution of this paper, therefore, is the development of
more efficient epigraph projection algorithms in many settings, often with practical (and theoretical) efficiency thatis nearly that of the equivalent proximal operator itself. In

Epigraph projections for fast general convex programming

developing these methods, we also develop some numberof new proximal algorithms as well, for which we know ofno existing methods in the literature. A summary of theoperator we develop is shown in the Appendix, Table 2.
In the remainder of this section, we detail several of theapproaches that we use to develop fast epigraph projection operators (and sometimes also proximal operators) forthe convex atoms listed in this table. Although many ofthe epigraph projection operators require specialized algorithms, they broadly fall into five categories: those solved
by an exact analytical solution, by a primal-dual Newton
method, by a new optimization approach we refer to as theimplicit dual Newton method, by a new approach we callthe sum-of-max algorithm for finding zeros of piecewise
linear functions, and orthogonal matrix methods. As a final
note, for all functions that apply elementwise to vectors,(i.e., the absolute value f (x) = |x|), we will be constructing epigraph projections for the sum over all the entries
epif (v, s) =1
k(x, t) − (v, s)k2
2
i f (xi )≤t

argmin
P

x,t:

(19)

as this strictly generalizes the case of constructing an epigraph for each element individually.
4.1. Exact methods
While many proximal operators have a closed form solution, a relatively smaller number of epigraph projections
have them. A few notable exceptions are the elementwisesquare f (x) = x2 and the `2 norm. For the square function, for instance, note that the dual epigraph problem canbe solved by the following derivation
1v
proxλ(·)2 (v) = argmin λx + (x − v)2 =2
1 + 2λx
λ=⇒ Pλf (v) =v2(18)
1 + 2λ


d
λv 2
1 2
− λ − λs = 0=⇒(17) dλ1 + 2λ 2
1=⇒ ( λ + s)(1 + 2λ)2 − v 2 = 0
2
(20)
2

which is a cubic equation in λ. Similarly the epigraph projection for the `2 norm is simply the well-known projection
on to the second order cone, which has a standard closed
form solution.
4.2. Primal-dual Newton method
For cases where the proximal operator itself has no closed
form, and the domain of the function is unconstrained, we
can employ a primal-dual Newton method to solve the resulting epigraph form. Intuitively, this simply involves us-

ing Newton’s method to directly solve the Lagrangian minimax problem for the epigraph projection operator, resulting
in the optimality conditions

x − v + λ∇x f (x) = 0,t−s−λr(x, t, λ) = (21)f (x) − twhich can be solved by computing the Newton direction ∆I + λ∇2x f (x) 0 ∇x f (x) ∆ = −r(x, t, λ)

0
1
−1
T
∇x f (x)−10(22)
and taking step sizes according to the standard line searchfor primal-dual methods (Boyd & Vandenberghe, 2004, pg.
612). Note that although we are optimizing over λ ≥ 0,
this constraint does not affect the optimization problem,
because we know the optimal solution must have λ > 0unless f (v) ≤ s, so we can treat the overall maximization
as an unconstrained problem. If f is an elementwise function, or if its Hessian has special structure (as in the caseof log-sum-exp, where the Hessian is diagonal plus rankone), then the Newton step can be computed in O(n) time,
and the overall complexity of the algorithm is O(n) times
the number of Newton iterations (which in theory can vary
between O(log 1/) or O(log log 1/), but which in practice is very small, often around 10 iterations to reach machine precision for these problems). Epigraph projectionsof the exponential, logistic, and log-sum-exp functions alltake this form.
4.3. Implicit dual Newton method
Sometimes the domain of the function is constrained andthe proximal operator has no closed-form solution, but canstill be solved efficiently. In these cases, the primal-dual
Newton method for epigraph projection is not usable, as
we would need to explicitly include the primal constraints.
However, it is still possible to solve the dual problem evenif we do not have a explicit dual formulation. Let
L(x, t, λ) ≡1k(x, t) − (v, s)k2 + λ(f (x) − t)2(23)be the Lagrangian for the epigraph projection problem.
Consider ∇(x,t) L(x) = 0 to be a constraint, whose solution
is the prox operator. Further, the solution (x(λ), t(λ)) =(proxλf (v), s + λ) can be regarded as a function of λ.
Thus, if we apply the implicit function theorem (see e.g.
Dontchev & Rockafellar, 2009) on ∇x L = 0, we have thatdx
= −(I + λ∇2x f (x))−1 ∇x f (x),
dλ

(24)

for all x = proxλf . Then, for the dual objectiveD(λ) ≡ min L(x, t, λ) = L(x(λ), t(λ), λ),
x,t

(25)

Epigraph projections for fast general convex programming

Algorithm 1 Implicit dual newton method for epigraphprojection
λ := 1
while not yet converged do
x := proxλf (v);dD(λ)
dλ := f (x) − λ − s;d2 D(λ):= −∇f (x)T (I + λ∇2 f (x))−1 ∇f (x) − 1;
dλ2
dD(λ)
if | dλ | ≤  then break;d2 D(λ)λ := max(0, λ − dD(λ)
dλ / dλ2 );end while
return (proxλf (v), s + λ).

we can derive the first and second derivative by the chain
rule and equation (25):
dD(λ)
d
=
L(x(λ), t(λ), λ) = f (x) − s − λ,
dλ
dλ

T
d2 D(λ)df (x)dx
d dD(λ)
=
=
− 1.
dλ2
dλ dλ
dx
dλ

(26)

Note that the relation holds only on x = proxλf (v). By
using these first and second derivatives, we can apply Newton’s method on the dual problem without computing it
directly, an approach we call the implicit dual Newton
method. If the proximal operator is also solved by a Newton method, the total number of iteration performed can be(#N ewton)2 , but in practice is still usually very small.
Negative
log epigraph The proximal operator ofP− i log(xi ) has the analytic solution
xi =

vi +

p

vi2 + 4λ,2where i = 1, . . . , n.

Algorithm 2 Linear time algorithm for absolute epigraphy := u, a = −s, b = 0;
while y is not empty do
Choose a mid point ym in y, assign λ = |ym |;Denote y(op) as the subvector satisfying [y (op) λ],and #(y(op) ) as the length of the vector;
Partition y to y< , y> , and y= by λ,g := a + ky> k1 − λ · (1 + b + #(y> ));if g < 0 thena := a+ky= k1 +ky> k1 , b := b+#(y= )+#(y> );y := y≤ ;else if g > 0 then y := y≥ ;else break;end while
xi := sign(xi ) max(0, |xi | − λ), ∀i,t = s + λ;return (x, t);

all piecewise linear functions. By inspecting the KKT conditions, we find that the dual solution λ of these problems
all satisfies the equation
F (λ) =nX

max(0, ai λ + yi ) + bλ + c = 0,

(30)

i=1

in which F (λ) is also a piecewise linear function. The
above equation can be solved in O(n) time by enumerating all the knot points and performing a quick-select ormedian-of-median algorithm (Cormen, 2009). The sum-klargest proximal operator can also be solved by a similar
equation with the max functions replaced by the clipping
functions. We present the simplified sum-of-max algorithm
for the absolute value epigraph in Algorithm 2.

(27)
Absolute value epigraph (`1 norm) As a concrete example, the epigraph projection problem of f (x) = kxk1
can be formulated as

By the implicit function theorem, we have
1

dxi
xi
=
, ∀x = proxλf (v).
dλ
1 + λ x12

(28)

i

n
X
d2
1
D
=
−λ
2
dλ2
x
i=1 i

1
1 + λ x12

such that kxk1 ≤ t.

The above problem is equivalent to solvingX
F (λ) =
max(0, |vi | − λ) − λ − s = 0

Thus, the derivatives of D can be computed asXd
D = −λ − s −log xi ,dλi1
min k(x, t) − (v, s)k2 ,
x,t 2

(31)

(32)

i

(29)

!
− 1.

which can be solved by the sum-of-max method. We thenrecover the epigraph projection (x, t) byi

xi = sign(vi ) max(0, |vi | − λ), t = λ + s.

(33)

4.4. Sum of max solvers
Several atomic convex functions are non-smooth, presenting problems for the above approaches. For example, thehinge loss, the absolute value, and the max of elements are

4.5. Orthogonal matrix epigraph solvers
Finally, we consider the case of functions that take matrix
inputs f : Rm×n → R. Although many such functions canEpigraph projections for fast general convex programming(34)

for all orthogonal U, V . The immediate implication of thisfact is that f can depend only on the singular values σ(X)of X, since given the singular value decomposition X =U SV TF (X) = F (U T XV ) = F (S).
(35)
A further consequence is that if f (X) = g(σ(X)) for some
vector (or elementwise) function g, then the epigraph projection of f can be computed from the epigraph projectionof g (a similar property holds for proximal operators)
T

epif (W, s) = ( U0 01 ) epig (σ(W ), s) ( V0 10 )

10 2
10 1

10 3

EpsilonSCS
ECOSRunning time (seconds)f (U XV T ) = f (X)

10 3
Running time (seconds)

be expressed just as equivalent vector/elementwise operations, we focus here on the class of orthogonally invariant
matrix functions, meaning that

10 0
10 -1
10 -2

10 4
10 3
Number of variables

10 2

10 2
10 1

Epsilon
SCS
ECOS

10 0
10 -1
10 -2
10 1

10 2
10 3
Number of variables

Figure 2. Scaling comparison on Lasso (left) for problems with nexamples and 10n variables, and for support vector data description (right) with 10n examples and n variables.

problems that have been reduced to cone form (with operators for cone projections) while Epsilon solves the higherlevel representation with a larger library of proximal andepigraph operators.

(36)

for W = U σ(W )V T . This property follows from the factthat the squared Frobenius norm is also orthogonally invariant. Using this property, we can immediately develop epigraph operators for matrix expressions such as the negative
logdet term (− log applied to singular values), the nuclear
norm (`1 norm applied to singular values), and the operator
norm (maximum singular value).

5. Experiments
In this section we compare our method to existing methods for general convex programming based on conic
solvers. Our approach is implemented in Epsilon (Epigraph
Proximal Solver), a library based on the ideas describedin this paper: Epsilon accepts general convex programs
specified according to the DCP ruleset, transforms them
to a sum of proximal and epigraph operators as in Theorem 1, and solves them by employing the library of operator implementations described previously. Epsilon is open
source and available at http://epopt.io, and all examples here are included in the distribution. We highlightthe benchmark problems briefly in this section, and include
a full description in Appendix A. Epsilon integrates directly with cvxpy (Diamond & Boyd, 2015) and thus we
make the natural comparison between Epsilon and the existing solvers which solve the conic form. In particular, we
compare Epsilon to ECOS (Domahidi et al., 2013), an interior point method, and SCS (O’Donoghue et al., 2013), the“splitting conic solver”. In general, interior point methods
achieve highly accurate solutions but have trouble scaling
to larger problems and so it is unsurprising that Epsilon isable to solve problems to moderate accuracy several ordersof magnitude faster than ECOS. On the other hand, SCS
employs an operator splitting method that is similar in spirit
to Epsilon, both being variants of ADMM. The main difference between Epsilon and SCS is in the intermediate representation to which operator splitting is applied: SCS solves

Standard problems We begin with a few standard machine learning problems: Lasso (Tibshirani, 1996), sparse
inverse covariance estimation (Banerjee et al., 2008) andimage classification on MNIST (LeCun et al., 1998). Table
1 shows that on these problems Epsilon achieves the samelevel of accuracy as existing approaches often several orders of magnitude faster. These are moderately-sized problems by ML standards, i.e. 1500 examples, 5000 features
for Lasso; 2000 examples, 1000 features for MNIST; and a
200 × 200 covariance estimation problem. However, interior point methods (ECOS) take minutes to hours on these
problems, demonstrating the lack of scalability traditionally associated with general convex programming. SCS(which itself is relatively recent work) does improve on
ECOS by an order of magnitude, but Epsilon is still 1-2
orders of magnitude faster. Figure 2 compares the runningtimes of each approach on the Lasso problem as problem
size increases: Epsilon solves larger problems in minutes
which require hours for SCS and ECOS.
Next, we consider several problems involving more complex objectives that go beyond the standard regularized
loss functions for regression and classification found in traditional machine learning models (again, described much
more fully in Appendix A).
Robust SVM We consider the “`∞ ” robust SVM (based
upon those in (Lanckriet et al., 2003; Shivaswamy et al.,2006) but with an `∞ uncertainy ball) described above by
the optimization problem (8).
Robust regression We consider the Chebyshev-like regression problem (Boyd & Vandenberghe, 2004, pg. 323),
which solves a least-squares problem with A in some unPk
certainty set described by Ā + i=1 ui Ai for kuk2 ≤ 1.
This leads to the optimization problem
minimize max kAi xk2 + |āTi x − bi | .
(37)x

i=1,...,kEpigraph projections for fast general convex programming

Lasso
MNIST (2000 examples)Sparse inverse covariance
Robust SVM
Robust regression
SV data description
Sum-k-largest softmax

Time
3.7s
2.53s
2.09s
52.9s
27.3s
40.5s
0.94s

EpsilonObjective
3.21 × 101
1.72 × 103
3.73 × 102
2.46 × 101
3.95 × 10−1
2.23 × 102
2.14 × 101

Time
20.5s
219.6s25.1s174.6s66.5s443.6s325.28s

SCS
Objective
3.21 × 101
1.72 × 103
3.73 × 102
2.48 × 103
3.95 × 10−1
2.23 × 102
2.14 × 101

Time215.7s1753.0s
141.6s67.9s893.3s
13.55sECOS
Objective
3.21 × 101
1.72 × 103
2.44 × 101
3.95 × 10−1
2.23 × 102
2.14 × 101

Table 1. Comparison of running time and objective value for Epsilon, SCS and ECOS, “-” indicates no result after 1 hour.

minimizeΘz[i] + λkΘk22 ,

(38)

i=1

10 2

EpsilonSCS

10 -2

Relative objective

k
X

10 1
Relative objective

Sum-k-largest softmax We consider a form of “worstcase” softmax regression, where we minimize the loss suffered by only the top k elements in the data set. This canbe written as the optimization problem

10 -5
10 -8
10 -11
10 -140

1020
Running time (seconds)

30

10 0

Epsilon
SCS

10 -2
10 -4
10 -6
10 -80200 400 600 800
Running time (seconds)

1000

where z[i] is the i-largest element of vector z, and zi =− log softmax(xi , yi , Θ) is the multiclass softmax loss.

Figure 3. Comparison on time and primal objective suboptimality
for Lasso problem (top) and SVDD problem (bottom).

Support vector data description (SVDD) Given a setof points, x1 , . . . , xm ∈ Rn , SVDD (Tax & Duin, 2004;
Chang et al., 2007) describes those points with an ndimensional Euclidean ball by solving

superior given reasonable running times, determining if we
can produce guaranteed bounds in our setting is an important question for future work. Finally, although it is not thefocus of this paper, we also compare Epsilon to several specialized solvers developed for specific classes of problems,and find that in many cases it is quite competitive, despite
being a fully generic approach; these results are presented
in Appendix B.

minimize
ρ,amX[kxi − ak22 − ρ)]+ + λ[ρ]+

(39)

i=1

with optimization variables ρ ∈ R and a ∈ Rn .
The appeal of general convex programming is that models such as these can be rapidly developed and prototyped;
however, as with simpler problems, existing solvers often have difficultly scaling to larger problem sizes. Table
1 shows that for these problems Epsilon is at least 2-3x
faster than existing approaches and in some cases reduces
running times by 5-10x. For support vector data description, Figure 2 shows that for larger problems Epsilon issignificantly faster than SCS and ECOS, e.g. for a dataset
x1 , . . . , x20000 ∈ R2000 , Epsilon requires roughly 6.5 minutes, SCS 32 minutes and ECOS more than 1 hour.
In Figure 3 we also highlight, for the Lasso and SVDD
problems, the evolution of primal objective versus time.
We see that in both cases, Epsilon consistently has a lower
objective value after the same amount of computation time;
this holds for all the examples in Table 1, though we do notinclude all plots for brevity. One element to note, however,is that unlike interior points methods, ADMM-based methods cannot easily produce guaranteed bounds on the suboptimality of a solution; although we do not pursue this issue
here, and simply show that the resulting solutions are far

6. Conclusion
This paper has focused on the development of fast optimization methods for general convex programs, here specified using a DCP modeling framework. We show that byimplementing just proximal and epigraph projection operators for a set of atomic convex functions, we can solve any
DCP built from those atoms. We have then implemented
a wide set of such operators and we show that the resulting solution methods can substantially outperform existing
approaches.
Acknowledgements This work was supported in part by
the National Science Foundation under award IIS-1320402,
and by support from Google. Matt Wytock was partly
funded by a scholarship from the Siebel Foundation.

References
Banerjee, Onureena, El Ghaoui, Laurent, andd’Aspremont, Alexandre. Model selection throughsparse maximum likelihood estimation for multivariate

Epigraph projections for fast general convex programming

gaussian or binary data. The Journal of Machine
Learning Research, 9:485–516, 2008.
Becker, Stephen R, Candès, Emmanuel J, and Grant,
Michael C. Templates for convex cone problems withapplications to sparse signal recovery. Mathematical
Programming Computation, 3(3):165–218, 2011.
Boyd, Stephen and Vandenberghe, Lieven. Convex optimization. Cambridge university press, 2004.
Boyd, Stephen, Parikh, Neal, Chu, Eric, Peleato, Borja, andEckstein, Jonathan. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends R in Machine Learning, 3(1):1–122, 2011.
Chang, Chien-Chung, Tsai, Hsi-Chen, and Lee, Yuh-Jye.
A minimum enclosing balls labeling method. 2007.
Chang, Chih-Chung and Lin, Chih-Jen. LIBSVM: A library for support vector machines. ACM Transactions
on Intelligent Systems and Technology, 2:27:1–27:27,
2011. Software available at http://www.csie.
ntu.edu.tw/˜cjlin/libsvm.
Chierchia, Giovanni, Pustelnik, Nelly, Pesquet, JeanChristophe, and Pesquet-Popescu, Béatrice. Epigraphical Projection and Proximal Tools for Solving Constrained Convex Optimization Problems: Part I. Technical report, October 2012.
Cormen, Thomas H. Introduction to algorithms. MITpress, 2009.
Diamond, Steven and Boyd, Stephen. CVXPY: A Pythonembedded modeling language for convex optimization.
Journal of Machine Learning Research, 2015. To appear.
Domahidi, Alexander, Chu, Eric, and Boyd, Stephen. Ecos:
An socp solver for embedded systems. In Control Conference (ECC), 2013 European, pp. 3071–3076. IEEE,
2013.
Donoho, David L. De-noising by soft-thresholding. Information Theory, IEEE Transactions on, 41(3):613–627,
1995.
Dontchev, Asen L and Rockafellar, R Tyrrell. Implicit functions and solution mappings. Springer, 2009.

Friedman, Jerome, Hastie, Trevor, and Tibshirani, Robert.
Sparse inverse covariance estimation with the graphicallasso. Biostatistics, 9(3):432–441, 2008.
Friedman, Jerome, Hastie, Trevor, and Tibshirani, Rob.
Regularization paths for generalized linear models via
coordinate descent. Journal of statistical software, 33(1):1, 2010.
Gabay, Daniel and Mercier, Bertrand. A dual algorithm forthe solution of nonlinear variational problems via finite
element approximation. Computers & Mathematics with
Applications, 2(1):17–40, 1976.
Grant, Michael and Boyd, Stephen. Graph implementations for nonsmooth convex programs. In Blondel, V.,
Boyd, S., and Kimura, H. (eds.), Recent Advances inLearning and Control, Lecture Notes in Control and Information Sciences, pp. 95–110. Springer-Verlag Limited, 2008.
Grant, Michael, Boyd, Stephen, and Ye, Yinyu. Disciplined convex programming. In Liberti, Leo and Maculan, Nelson (eds.), Global Optimization, volume 84 of
Nonconvex Optimization and Its Applications, pp. 155–
210. Springer US, 2006. ISBN 978-0-387-28260-2. doi:
10.1007/0-387-30528-9 7.
Harizanov, Stanislav, Pesquet, Jean-Christophe, and Steidl,
Gabriele. Scale Space and Variational Methods in Computer Vision: 4th International Conference, SSVM 2013,
Schloss Seggau, Leibnitz, Austria, June 2-6, 2013. Proceedings, chapter Epigraphical Projection for Solving
Least Squares Anscombe Transformed Constrained Optimization Problems, pp. 125–136. Springer Berlin Heidelberg, 2013.
Hsieh, Cho-Jui, Sustik, Mátyás A., Dhillon, Inderjit S., andRavikumar, Pradeep D. Sparse inverse covariance matrix estimation using quadratic approximation. CoRR,
abs/1306.3212, 2013.
Hsieh, Cho-Jui, Sustik, Mátyás A, Dhillon, Inderjit S, andRavikumar, Pradeep. Quic: quadratic approximation forsparse inverse covariance estimation. The Journal of Machine Learning Research, 15(1):2911–2947, 2014.

Douglas, Jim and Rachford, Henry H. On the numericalsolution of heat conduction problems in two and threespace variables. Transactions of the American mathematical Society, pp. 421–439, 1956.

Joachims, Thorsten. Optimizing search engines using
clickthrough data. In Proceedings of the eighth ACM
SIGKDD international conference on Knowledge discovery and data mining, pp. 133–142. ACM, 2002.

Fan, Rong-En, Chang, Kai-Wei, Hsieh, Cho-Jui, Wang,
Xiang-Rui, and Lin, Chih-Jen. Liblinear: A library
for large linear classification. The Journal of Machine
Learning Research, 9:1871–1874, 2008.

Lanckriet, Gert RG, Ghaoui, Laurent El, Bhattacharyya,
Chiranjib, and Jordan, Michael I. A robust minimax approach to classification. The Journal of Machine Learning Research, 3:555–582, 2003.

Epigraph projections for fast general convex programming

LeCun, Yann, Cortes, Corinna, and Burges, Christopher JC. The mnist database of handwritten digits, 1998.
O’Donoghue, Brendan, Chu, Eric, Parikh, Neal, and Boyd,
Stephen. Operator splitting for conic optimization
via homogeneous self-dual embedding. arXiv preprint
arXiv:1312.3039, 2013.
Optimization, Gurobi et al. Gurobi optimizer reference
manual. URL: http://www. gurobi. com, 2012.
Parikh, Neal and Boyd, Stephen. Proximal algorithms.
Foundations and Trends in optimization, 1(3):123–231,
2013.
Platt, John C. Advances in kernel methods. chapter Fast
Training of Support Vector Machines Using Sequential
Minimal Optimization, pp. 185–208. MIT Press, Cambridge, MA, USA, 1999. ISBN 0-262-19416-3.
Rahimi, Ali and Recht, Benjamin. Random features forlarge-scale kernel machines. In Advances in neural information processing systems, pp. 1177–1184, 2007.
Rockafellar, R Tyrrell. Monotone operators and the proximal point algorithm. SIAM journal on control and optimization, 14(5):877–898, 1976.
Ryu, Ernest K and Boyd, Stephen. Primer on monotone
operator methods. To appear, Appl. Comput. Math., 15(1), 2016.
Shalev-Shwartz, Shai, Singer, Yoram, Srebro, Nathan, andCotter, Andrew. Pegasos: Primal estimated sub-gradientsolver for svm. Mathematical programming, 127(1):3–
30, 2011.
Shivaswamy, Pannagadatta K, Bhattacharyya, Chiranjib,
and Smola, Alexander J. Second order cone programming approaches for handling missing and uncertain
data. The Journal of Machine Learning Research, 7:1283–1314, 2006.
Tax, David MJ and Duin, Robert PW. Support vector data
description. Machine learning, 54(1):45–66, 2004.
Tibshirani, Robert. Regression shrinkage and selection via
the lasso. Journal of the Royal Statistical Society. Series
B (Methodological), pp. 267–288, 1996.
Tofighi, Mohammad, Bozkurt, Alican, Kose, Kivanc, andCetin, A Enis. Deconvolution using projections onto theepigraph set of a convex cost function. In Signal Processing and Communications Applications Conference(SIU), 2014 22nd, pp. 1638–1641. IEEE, 2014.

Epigraph projections for fast general convex programming

A. Experiment details

trivially be put in proximal form withIn this section we provide the details of the experiments
corresponding to the results presented in Section 5. In eachcase we describe the underlying problem, the initial optimization problem, and how this translates into proximal
form for the Epsilon solver (that is, a sum of objective terms
where each term has an efficient proximal operator or epigraph projection, plus a set of linear equality constraints),
using the transformations from Section 3. All these examples are included in the Epsilon distribution, available at
http://epopt.io.

f1 (Θ) = − log |Θ| + tr SΘ,f2 (Θ) = λkΘk1 .

(43)

where in the construction of f1 , we have exploited the factthat any proximal operator can easily be combined witha linear function (see e.g. Parikh & Boyd (2013) for details). In our experiments, we construct the sample covariance from 100 samples drawn from N (0, Θ−10 ) whereΘ0 ∈ S200 has uniform random entries and is 1% nonzero.
A.3. MNIST (2000 images)

A.1. Lasso
The Lasso problem isminimize (1/2)kXθ − yk22 + λkθk1 ,θ

(40)with input features X ∈ Rm×n , response variables y ∈Rm , and model parameters θ ∈ Rn . The regularization parameter λ ≥ 0 controls the tradeoff between data fit andthe `1 regularization term which encourages sparsity in themodel parameters. The Lasso is especially useful in thehigh-dimensional case where m < n as the sparsity induced by `1 regularization effectively controls the numberof free parameters in the model, see Tibshirani (1996) for
details. The Lasso problem as written is already in proximal form withf1 (θ) = (1/2)kXθ − yk22 ,f2 (θ) = λkθk1 .

(41)
1500×5000

In our experiments, we generate X ∈ R
froma standard Normal distribution and set y = Xθ0 + with θ0 having 1% nonzero standard Normal entries and ∼ N (0, 0.052 ). We set the regularization parameter
λ = 0.5kX T yk∞ .
A.2. Sparse inverse covariance

The MNIST dataset (LeCun et al., 1998) consists of handwritten digits, constructed with the goal of building a classifier for automatically recognizing each digit ({0, . . . , 9})in each image. As a linear classifier applied directly to theraw pixels performs poorly, we generate random Fourierfeatures (Rahimi & Recht, 2007) and train a classifier using sparse softmax regression
minimize `(X, y; Θ) + λkΘk1
Θwhere X ∈ Rm×n are the image features, y ∈ {0, . . . , 9}m
are the image labels, λ ≥ 0 is the regularization parameter and the softmax loss, parameterized by Θ ∈ Rn×10 , isgiven by
!
m
10
X
X
T
T
`(X, y; Θ) =log
exp(xi θk ) − xi θyi . (45)
i=1

k=1

Unlike the least squares loss employed in the Lasso, thesoftmax loss cannot easily be composed with an arbitrary
linear function and this requires the introduction of an additional auxiliary variable, Z ∈ Rm×10 . With this additionalvariable, the proximal form for this problem is given by
!
m
10
X
X
f1 (Z) =log
exp(Zik )i=1

Sparse inverse covariance estimation models a multivariate
Gaussian distribution over n variables by solving the optimization problem(44)

f2 (Θ) = λkΘk1 −

k=1
m
X

(46)
xTi θyii=1

with equality constraints
minimize − log |Θ| + tr SΘ + λkΘk1
Θ(42)where S ∈ Sn is the sample covariance, Θ ∈ Sn is theestimated inverse covariance and the `1 norm k · k1 isapplied elementwise. As with the Lasso, the `1 penalty
promotes sparse structure and is especially useful in thehigh-dimensional case with more variables than samples(m < n), see Friedman et al. (2008) for details. With respect to our framework, given a proximal operator for the− log | · | term, sparse inverse covariance estimation can

Z = XΘ

(47)

In our experiments, we train the classifier on 2000 images
using 1000 random Fourier features and λ = 0.1.
A.4. Robust SVM
In our experiments we use an `∞ variant of the support vector machine, similar to the formulations in (Lanckriet et al.,
2003; Shivaswamy et al., 2006), but with an `∞ uncertainty

Epigraph projections for fast general convex programming
Function

Matrix X ∈ Rn×n

Vector x ∈ Rn

Elementwise x, y ∈ R

Type

Atom

Definition

Absolute

f (x) = |x| (≡ kxk1 )Squaref (x) = x2 (≡ kxk22 )Hinge

f (x) = max{x, 0}

Deadzonef (x) = max{|x| − , 0},  ≥ 0

Quantilef (x) = max{αx, (α − 1)x}, 0 ≤ α ≤ 1

Logisticf (x) = log(1 + exp(x))

Inverse positivef (x) = 1/x, x ≥ 0

Negative log

f (x) = − log(x), x ≥ 0

Exponential

f (x) = exp(x)Negative entropyf (x) = x · log(x), x ≥ 0

KL Divergencef (x, y) = x · log(x/y), x, y ≥ 0

Quadratic over linear

f (x, y) = x2 /y, y ≥ 0`2 -normf (x) = kxk2

MaximumLog-sum-exp

f (x) = maxi xi (kxk∞ ≡ maxi |xi |)Pf (x) = ki=1 x[i] (x[i] ≥ x[i+1] )Pnf (x) = log
i=1 exp(xi )Negative log det
Nuclear norm
Spectral normf (X) = − log det(X), X ∈ Sn
f (X) = kσ(X)k1 , X ∈ Rm×n
f (X) = kσ(X)k∞ , X ∈ Rm×n

Sum-k-largest

Proximal / epigraph operator
MethodComplexity
soft thresholdingO(n)sum-of-max
O(n)linear equation
O(n)cubic equation
O(n)soft thresholdingO(n)sum-of-max
O(n)soft thresholdingO(n)sum-of-max
O(n)soft thresholdingO(n)sum-of-max
O(n)NewtonO(n) · (# Newton)
Primal-dual Newton
O(n) · (# Newton)
cubic equationO(n)Implicit dual Newton
O(n) · (# Newton)
quadratic equationO(n)Implicit dual Newton
O(n) · (# Newton)
NewtonO(n) · (# Newton)
Primal-dual Newton
O(n) · (# Newton)
O(n) · (# Newton)
Projected Newton
Implicit dual Newton
O(n) · (# Newton)2
O(n) · (# Newton)
Projected Newton
Implicit dual Newton
O(n) · (# Newton)2
O(n) · (# Newton)
Projected Newton
Implicit dual Newton
O(n) · (# Newton)2
group soft thresholding
O(n)analytic projection
O(n)sum-of-max
O(n)sum-of-max
O(n)sum-of-clipO(n)bisectionO(n) · (# Bisection)
NewtonO(n) · (# Newton)
Primal-dual Newton
O(n) · (# Newton)− log on λ(X)k · k1 on σ(X)k · k∞ on σ(X)O(n3 )O(n3 )O(n3 )

Table 2. Complete list of proximal and epigraph projection operators implemented for this work. Most proximal operators (except
sum-k-largest) have appeared in some form in previous literature, but epigraph projections are typically novel to this work.

ball instead of an `2 uncertainty ball. In this setting, given
input data (xi , yi ) we wish to train an SVM by minimizing
the standard regularized hinge loss,mX
λminimizekθk22 +
max{0, 1 − yi · x̄Ti θ}
θ
2
i=1

(48)

but where x̄i lies in some uncertainty set centered at xi ,x̄i = xi + P u where kuk∞ ≤ 1. This can be expressed as
the optimization problem
minimizeθmX

λkθk22 +
sup max{0, 1−yi ·θT (x̄i +P ui )}2ku
i k∞ ≤1i=1(49)

which, using the relation that supkuk∞ ≤1 cT u = kck1
where, is equivalent to the optimization problemmX
λkθk22 +
max{0, 1 − yi · θT xi + kP T θk1 }.
θ
2
i=1
(50)
As discussed in Section 3, this is transformed to proximal
form

minimize

λ
kθk22
2
X
f2 (z3 ) =max{z3 , 0}f1 (θ) =f3 (z1 , z2 ) = I{kz1 k1 ≤ z2 }(51)

Epigraph projections for fast general convex programmingand equality constraints
z1 = P T x
z3 = 1 − diag(y)Xθ + 1z2 .

(52)minimize max kAk xk2 + |āTk x − bk | .

2500×750

In our experiments, we generated X ∈ R
randomuniform [0,1], θ ∈ R750 also random uniform, and set y =sign(xTi θ + N (0, 0.1)). To create well-separated points,
we further added xi ← xi + 0.7yi · θ, and chose P =diag(N(0, 750)).
Support vector data description Given a set of unlabeled points, x1 , . . . , xm ∈ Rn , support vector data description (Tax & Duin, 2004; Chang et al., 2007) describes
those points with an n-dimensional Euclidean ball by solving
minimize
ρ,amX[kxi − ak22 − ρ)]+ + λ[ρ]+

(53)

i=1

with optimization variables ρ ∈ R and a ∈ Rn . The firstterm penalizes points outside a ball centered at a with ra√dius ρ while the second term regularizes the radius withλ ≥ 0 controlling the tradeoff. This problem is transformed
to proximal form with operators
f1 (t, ρ) =m
X[ti ]+ + λ[ρ]+i=1

f3 (a, s) =m
X

For robust regression (Boyd & Vandenberghe, 2004, pg.
323), we wish to find a solution x that minimize this worst
possible error,(54)
I(kxi − ak22 ≤ si )

x

k

(58)

With additional variable t, u, v, p, q ∈ Rk , this problem istransformed to proximal form with operators
f1 (t) = max (ti ),i=1,...,kf2 (p, u) =kXI(kpi k2 ≤ ui ),i=1

f3 (q, v) =

k
X

(59)

I(|qi | ≤ vi ),i=1

with equality constraints
t = u + v,
pi = Ai x,(60)

qi = āTi x − bi , ∀i = 1, . . . , k.
In the experiment, we generate the Ā, Ai , and b from
uniform distribution, then normalize Ā and Ai to a unit
ball. We choose p = 5000 and A, Ai ∈ R10×200 , fori = 1, . . . , p.
sum-k-largest softmax The softmax loss is a multiclassloss function defined asi=1

exp(xT Θyi ),softmax(x, y, Θ) = P
Tk exp(x Θk )and equality constraintt = s − ρ.

(55)

In our experiments, we generate 5000 random points uniformly over the 200-dimensional unit hypersphere and thenchoose 100 outliers at random and add noise,  ∼ N (0, I);
we fit the model with λ = 1.
Robust Regression Consider a noisy matrix bounded bya unit ball over some known perturbation directions Ai ,where Θ ∈ Rn×c is the weight for each class. The softmax
loss is commonly used in, for example, the regularized logistic regression, which can be formulated by the softmax
loss plus a regularization term. Here, we consider to minimize the worst k loss incurred from the regression. I.e., weonly minimize
minimizeΘA = {Ā + c1 A1 + . . . + cp Ap | kpk2 ≤ 1}.

(56)

The maximum error incurred from performing linear regression over the uncertain set can be written as
sup k(Ā + c1 A1 + . . . + cp Ap )x − bk∞ksup ck (Ak x) + (āTk x − bk )kck2 ≤1= max kAk xk2 + |āTk x − bk | .
k

(57)

k
X

z[i] + λkΘk22 ,

(62)

i=1

where z[i] is the i-largest element of vector z, and zi =− log softmax(xi , yi , Θ) is the multiclass softmax loss.
With additional variable ui ∈ Rc , ∀i = 1, . . . , m, thisproblem is transformed to proximal form with operators

kck2 ≤1

= max

(61)

f1 (Θ) = −mX

xTi Θyi + λkΘk22 ,

i=1

f2 (z) = sum-k-largest(z),f3 (Θ, z) = I(log-sum-exp(ui ) ≤ zi ),(63)

Epigraph projections for fast general convex programming

SolverliblinearglmnetGurobi
QUIC

Problemhinge l1
hinge l1 sparse
hinge l2
hinge l2 sparselasso
lasso sparselogreg l1
logreg l1 sparse
mv lasso
lp
qp
covselTime
Epsilon Solver
3.71s
0.49s
14.26s
4.26s
3.58s
0.16s
1.82s
0.83s
3.69s
0.84s
13.58s
0.67s
3.70s
2.31s
6.69s
1.96s
7.14s
7.40s
0.33s
6.02s
1.39s
4.12s
0.93s
6.24s

mature software packages readily available in common
mathematical programming environments (e.g. Matlab, R,
Python). Furthermore, even when specialized solvers areavailable, translating problems to the interface provided bya particular package requires effort to understand and conform to the idiosyncrasies of each implementation. In contrast, general convex programming offers a uniform syntax
and interface allowing problems to be easily formulated,
extended and solved.

Table 3. Comparison of running times between Epsilon and specialized solvers.

with equality constraints
ui = xTi Θ, ∀i = 1, . . . , m.

(64)

In the experiment we choose X ∈ R400×10 , k = 5, andthe number of classes to be 120. We generate the datafrom normalized uniform distribution, and assign classes
uniformly.

B. Comparison with specialized solvers
In this section, we compare Epsilon to an assortment ofspecialized solvers which are available for common problems. Before doing so, we emphasize that the generalconvex programming approach offers many advantages tospecialized algorithms in terms of reuse and extensibility.
In addition, most convex problems do not have dedicated,

In terms of running times, Table 3 compares Epsilon to four
dedicated software packages implementing specialized algorithms: liblinear (Fan et al., 2008), glmnet (Friedman
et al., 2010), Gurobi (Optimization et al., 2012) and QUIC
(Hsieh et al., 2013). The default stopping criteria is used foreach solver, corresponding to moderate accuracy for Epsilon and high accuracy for the specialized solvers. For themost part, Epsilon is competitive, although on a few problems liblinear and glmnet are significantly faster. This isdue to the ability of these specialized algorithms to exploit
sparsity in the solution which arises due to `1 regularization (Lasso problems) or a small number of support vectors
in the dual SVM formulation (hinge problems). At present,
Epsilon does not take advantage of such structure and thus
may have a disadvantage on sparse problems.
On the other hand, Table 3 shows that Epsilon is significantly faster than Gurobi and QUIC in solving linear/quadratic programs and sparse inverse covariance estimation, respectively. However, it is important to highlightthat the specialized algorithms solve these problems to highaccuracy (e.g. tolerances of 10−8 or smaller) while Epsilontargets only moderate accuracy (e.g. 10−3 ). For moderate accuracy, the operator splitting approach can be highly
competitive, allowing Epsilon to be significantly faster onsome problems.