
\section{Introduction}\label{s-introduction}
Our story starts with a well studied problem, minimizing
a convex function that is the sum of two convex functions
with different access methods, referred to as a \emph{composite function}
\cite{nesterov2013gradient}.
The first function is smooth, and we can access it only by a few methods, such
as evaluating its value and gradient at a given point.
The other function is not necessarily smooth, but is structured, and 
we can access it only by solving a convex optimization problem that involves it.
In the typical setting, the second function is one for which we can 
efficiently compute its proximal operator, usually analytically.
Here we assume a bit more for the second function, specifically, 
that we can minimize it plus a structured function of modest complexity.

Our goal is to minimize such functions automatically, with a method 
that works well across a large variety of problem instances using its
default parameters.
We leverage new technological developments: systems for 
automatic differentiation that automatically and effficiently compute
gradients given a computation graph description, and
systems for solving convex optimization problems
described in a domain specific language for multiple parameter values.

\subsection{Oracle-structured composite function} \label{s-assumptions}
We seek to minimize $h(x) = f(x)+g(x)$ over $x\in {\mbox{\bf R}}^n$.
We assume that
\begin{itemize}
\item $f:\Omega\to{\mbox{\bf R}}$ is convex and differentiable, where 
$\Omega \subseteq {\mbox{\bf R}}^n$ is convex and open.
We assume that a point $x^0 \in \Omega$ is known.
\item $g: {\mbox{\bf R}}^n \to {\mbox{\bf R}}\cup \{+\infty\}$ is convex, with closed sublevel
sets, and not necessarily differentiable.  Infinite values of $g$ encode
constraints on $x$.
\end{itemize}
Our method will access these two functions in very specific ways.
\begin{itemize}
\item We can evaluate $f(x)$ and $\nabla f(x)$ at any $x$.
For $x\not\in \Omega$, our oracle returns the
value $+\infty$ for $f(x)$.
(We will discuss a few other possible access methods for $f$ in the sequel.)
\item We can minimize $g(x)$ plus another structured function.
Here we use the term structured loosely, to mean in the sense of Nesterov,
or as a practical matter, in a disciplined convex programming (DCP)
description.
This is an extension of the usual assumption in composite function
minimization, where the usual assumption is that the proximal 
operator of $g$ can be evaluated analytically.
\end{itemize}
We refer to $f$ as the oracle part of the objective,
$g$ as the structured part of the objective,
and $h=f+g$ as an oracle-structured composite objective.
We denote the optimal value of the problem as
\[
h^\star = \inf_x \left( f(x) + g(x) \right).
\]

\paragraph{Assumptions.}
We assume that the sublevel sets of $h$ are compact, 
and $h^\star < \infty$ (so at least some sublevel sets are nonempty).
These assumptions imply that $f(x)+g(x)$ has a minimizer, {\it i.e.}, 
a point $x^\star$ with $h(x^\star)=h^\star$.
Our convergence proofs make the typical assumption that $\nabla f$ is Lipschitz 
continuous with constant $L$, but we stress that this is not used in the 
algorithm itself.

\paragraph{Optimality condition.}
The optimality condition is \begin{equation}
\label{e-opt-cond}
\nabla f(x)+ q = 0, \qquad q \in\partial g(x),
\end{equation}
where $\partial g(x)$ denotes the subdifferential of $g$ at $x$.
For $x \in \Omega$ and $q \in \partial g(x)$, 
we can interpret $\nabla f(x)+q$ as a
residual in the optimality condition.
(We will use this in the stopping criterion of our algorithm.)
Our access to $f$ directly gives us $\nabla f(x)$; we will see that our
access method to $g$ indirectly produces a subgradient $q \in \partial g(x)$.

\subsection{Practical considerations}
We are motivated by two technological considerations related to 
our access methods to $f$ and $g$, which we mention briefly here.

\paragraph{Handling $f$.}
To handle $f$ we can rely on automatic differentiation systems
that have been developed in recent years, such as PyTorch \cite{NEURIPS2019_9015},
TensorFlow \cite{tensorflow2015-whitepaper}, and Zygote/Flux/Autograd \cite{maclaurin2015autograd,Flux.jl-2018,innes2019differentiable}.
Automatic differentiation is an old topic \cite{nolan1953analytical, adamson1969slang}
(see \cite{baydin2018automatic} for a recent review), but these
recent systems go way beyond the basic algorithms for automatic
differentiation, in terms of ease of use and run-time efficiency,
across multiple computation platforms.
We describe $f$ (but not its
gradient) using existing libraries and languages;
thereafter, $f(x)$ and $\nabla f(x)$ can be 
evaluated very efficiently, on many computation platforms, 
ranging from single CPU to multiple GPUs.
These systems are widely used throughout machine learning,
mostly for fitting deep neural networks.

\paragraph{Handling $g$.}
To handle $g$ we make use of domain specific languages (DSLs) 
for convex optimization, such as CVX~\cite{grant2014cvx}, 
CVXPY~\cite{cvxpy_paper,agrawal2018rewriting}, 
and Convex.jl~\cite{cvxjl}. These systems take a description of $g$
in a special language based on disciplined convex programming (DCP)~\cite{GBY:06}.
This description of $g$ is then automatically transformed into a 
standard form, such as a cone program, and then solved.
Recently, such systems have been enhanced to include parameters,
which are constants in the problem each time it is solved, 
but can be changed and the problem re-solved efficiently, skipping
the compilation process.
These systems are reasonably 
good at preserving structure in the problem during the compilation,
so the solve times can be quite small when exploitable structure 
is present, which we will see is the case in our method.

We mention that $g$ can contain hidden additional variables.
By this we mean that $g$ has the form $g(x) = \inf_z G(x,z)$,
where $G$ is convex in $(x,z)$.  Roughly speaking, $z$ contains
hidden variables that do not appear in $f$.  Such functions
are immediately handled by structured systems, without 
any additional effort.   In particular, we do not need to work out
an analytical form for $g(x)$ in this case.
In this case just evaluating $g$ requires solving an optimization problem
(over $z$).  Our method will avoid any evaluations of $g$.

\paragraph{When to just use a structured solver.}
Finally, we mention that if $f$ is simple enough to be handled by
a DCP-based system, then simply minimizing $f+g$ using such a system
is the preferred method of solution.
We are interested here in problems where this is not the case.
Typically this means that $f$ is complex in the sense of 
involving substantial data, for example, a sample average of some
function with $10^6$ or more samples.
(We will see this phenomenon in the numerical examples given in 
\S\ref{s-experiments}.)

\paragraph{Contribution.}
The method we propose in the next section is in the family of
variable metric bundle methods, and closely related to a number of 
other methods found in the literature (we review related 
work in \S\ref{s-related-work}).  As an algorithm, our method is not 
particularly novel; we consider our contribution to be its careful 
design to be compatible with our access methods.

\section{Oracle-structured minimization method} \label{s-method}
In this section we propose a generic method for solving the oracle-structured
minimization problem,
which we call \emph{oracle-structured minimization method} (OSMM).
OSMM combines several well known methods from optimization,
including variable metric or quasi-Newton curvature estimates 
to accelerate convergence,
bundle methods that build up a piecewise affine model, and two types of 
damping, based on a trust penalty and a line search.
These are chosen to be compatible with our access methods.

We will denote the iterates with a superscript, so $x^k$ denotes the $k$th
iterate of the algorithm.
We will let $x^{k+1/2}$ denote the tentative iterate at the $(k+1)$st iteration,
before the line search.
We will assume that $x^0\in\Omega$, {\it i.e.}, $f(x^0)<\infty$.  Our algorithm
will guarantee that $x^k \in \Omega$ for all $k$.
It is a descent method, {\it i.e.}, $h(x^{k+1}) < h(x^k)$.  While $h(x^0)=\infty$
is possible, we will see that $h(x^k)< \infty$ for $k\geq 1$.

As with many other optimization algorithms, OSMM is 
based on forming an approximation of the function $f$ in each
iteration.

\subsection{Approximation of the oracle function}\label{s-oracle_approximation}
In iteration $k$, we form a convex approximation of $f$, given by 
$\hat f_k: {\mbox{\bf R}}^n \to {\mbox{\bf R}} \cup \{\infty\}$, based 
on information obtained from previous iterations and possibly prior 
knowledge of $f$.  
The approximation has the specific form
\begin{equation}\label{e-fhat}
\hat f_k(x) = l_k(x) + (1/2)(x-x^k)^T H_k (x-x^k).
\end{equation}
Here $H_k$ is positive semidefinite,
and $l_k:{\mbox{\bf R}}^n \to {\mbox{\bf R}} \cup \{\infty\}$ is a convex
minorant of $f$, {\it i.e.}, $l_k(x) \leq f(x)$ for all $x$.

\paragraph{Assumptions on $H_k$.}
The only assumption we make about $H_k$ is that it is positive
semidefinite and bounded, {\it i.e.}, there exists a $C$ such that
$\|H_k \|_2 \leq C$ for all $k$.
In practice, $H_k$ accelerates convergence 
by serving as an estimate of the curvature of $f$.  
The simplest choice, $H_k=0$, results in algorithm that converges
but does not offer the practical benefit of 
convergence acceleration.

\paragraph{Choice of $H_k$.}
There are many ways of choosing a positive semi-definite $H_k$ 
to approximate the curvature of $f$ at $x^k$.
One obvious choice is the Hessian $H_k = \nabla^2 f(x^k)$,
but this requires that $f$ be twice differentiable, and also 
violates our assumption about how we access $f$.
A lesser violation of the access method might use an
approximation of the Hessian based on evaluations of 
the mapping $z \mapsto \nabla^2f(x^k) z$ ({\it i.e.}, Hessian-vector multiplication),
which can be practical in many cases
\cite{erdogdu2015convergence}.
A simple and effective choice is
$H_k = (a_k/n) I$, where $a_k$ is an approximation of $\mathop{\bf Tr} \nabla^2 f(x^k)$,
obtained for example by the Hutchinson 
method \cite{hutchinson1989stochastic,meyer2021hutch++}.

Quasi-Newton methods are a general class of curvature approximations 
that are compatible with our assumptions on access method for $f$.
These methods, which have a very long history,
build up an approximation of $H_k$ using only the current 
and previously evaluated gradients
\cite{davidon1959variable,fletcher1963rapidly,
broyden1965class,dennis1977quasi,byrd1994representations}.
When $H_k$ is low rank, or diagonal plus low rank, the method is 
practical even for large values of $n$.  (Such methods are often called
limited memory, since they do not require the storage of 
an $n \times n$ matrix.)
For OSMM we propose to use the
low-rank quasi-Newton choice given in \cite{fletcher2005new},
described in detail in \S\ref{appendix:s-curvature}.
We can express $H_k$ as
\begin{equation}\label{PSD_low_rank}
H_k = G_k G_k^T,
\end{equation}
where $G_k\in{\mbox{\bf R}}^{n\times r}$, and $r$ is a chosen (maximum) rank
for $H_k$.

As many others have observed, limited memory quasi-Newton methods
deliver most of their benefit for relatively small values of $r$, 
like $r=10$ or $r=20$.   These values allow the methods to be used 
even when $n$ is very large (say, $10^5$), since the storage requirement
(specifically, of $G_k$) grows linearly with $r$, and the computational
cost of evaluating $x^{k+1/2}$ grows quadratically in $r$,
and only linearly in $n$.

\paragraph{Assumptions on $l_k$.}
We make the usual assumption on the minorant $l_k$ that it is tight
at $x^k$, {\it i.e.}, $l_k(x^k)=f(x^k)$.
It follows that $\hat f_k(x^k) = f(x^k)$.
It also follows that $l_k$ is differentiable at $x^k$, and 
$\nabla l_k(x^k) = \nabla f(x^k)$.
To see this, we note that since $l_k$ is a minorant of $f$,
tight at $x^k$, we have 
\[
\partial l_k(x^k) \subseteq \partial f(x^k) = \{ \nabla f(x^k) \}.
\]
The first inclusion can be seen since any affine lower bound on $l_k$,
tight at $x^k$, is also an affine lower bound on $f$, tight at $x^k$,
so its linear part is a subgradient of $f$ at $x^k$.
The right-hand equality holds since $f$ is differentiable, so its 
subdifferential contains only one element, its gradient.
Finally, since $\partial l_k(x^k)$ contains only $\nabla f(x^k)$,
we conclude it is differentiable at $x^k$, with gradient $\nabla f(x^k)$.

In addition to the mathematical assumptions about $l_k$ described above,
we will assume that $l_k$ has a structured description.
This implies that $\hat f_k$ has a structured description.

\paragraph{Minorants.}
The simplest minorant is the first order Taylor approximation
\[
l_k(x) = f(x^k) + \nabla f(x^k)^T(x-x^k).
\]
A more complex minorant is the piecewise affine minorant
\[
l_k(x) = \max_{i=1,\ldots, k} \left( f(x^i) + \nabla f(x^i)^T(x-x^i)\right),
\]
which uses all previously evaluated gradients of $f$.

For OSMM we propose the piecewise affine minorant
\begin{equation}\label{l_k_PWA}
l_k(x) = \max_{i = \max\{0, k-M+1\},\ldots,k} \left(
f(x^i) + \nabla f(x^i)^T(x-x^i) \right),
\end{equation}
the pointwise maximum of the affine minorants from the previous
$M$ gradient evaluations, where $M$ is the memory.
With memory $M=1$, this reduces to the Taylor approximation.

The problem of choosing the memory $M$ is very similar to the problem of
choosing $r$, the rank of the curvature approximation.  The storage 
requirements grows linearly with $M$, and the computational cost
of evaluating $x^{k+1/2}$ grows quadratically with $M$.
As with the choice of $r$, small values such as $M=10$ or $M=20$
seem to work well in practice.

We mention a few additional useful minorants that use additional 
prior information about $f$ or its domain $\Omega$.
First, we can add to $l_k$ constraints that contain $\Omega$.
Suppose we know that $\tilde \Omega \supset \Omega$, where 
$\tilde \Omega$ has a structured description, {\it e.g.}, a box.
We can then use the minorant
\[
l_k(x) = \max_{i = \max\{0, k-M+1\},\ldots,k} \left(
f(x^i) + \nabla f(x^i)^T(x-x^i) \right)
+ I_{\tilde \Omega}(x),
\]
where $I_{\tilde \Omega}$ is the indicator function of $\tilde \Omega$.
In a similar way, if a (constant) lower bound $\ell$ on $f$ is known, 
we can replace any minorant $l_k$ with $\max\{l_k(x),\ell\}$.

If it is known that $f$ is $\rho$-strongly convex, 
we can strengthen the piecewise affine minorant \eqref{l_k_PWA} 
to the piecewise quadratic minorant
\[
l_k(x) = \max_{i=\max\{0, k-M+1\},\ldots,k} \left( f(x^i) 
+ \nabla f(x^i)^T(x-x^i) + (\rho/2)\|x-x^i\|_2^2 \right).
\]
(Each term in the maximum has the same quadratic part $(\rho/2)\|x\|_2^2$,
so $l_k$ can be expressed as a piecewise affine function plus 
$(\rho/2)\|x\|_2^2$.)

\paragraph{Lower bound.}
We observe that
\[
\ell_k = \inf_{x} \left( l_k(x) + g(x) \right)
\]
is a lower bound on the optimal value $h^\star$.
It can be computed using a system for structured optimization.
At iteration $k$ we let $L_k$ denote the best (largest) lower bound
found so far,
\begin{equation}\label{e-Lk}
L_k = \max\{\ell_1, \ldots, \ell_k\}.
\end{equation}


\subsection{Tentative update}\label{s-tentative-update}
At iteration $k$,
our tentative next iterate $x^{k+1/2}$ is obtained by minimizing
our approximation of $f$, plus $g$ and a trust penalty term:
\begin{equation}\label{e-xk+half}
x^{k+1/2} = \mathop{\rm argmin}_{x} \left( 
\hat f_k(x) + g(x) + \frac{\lambda_k}{2} \|x-x^k\|_2^2
\right).
\end{equation}
The last term is a (Levenberg-Marquardt or proximal) trust penalty, 
which penalizes deviation from $x^k$;
the positive parameter $\lambda_k$ scales the trust penalty.
We assume that $x^{k+1/2}$ in \eqref{e-xk+half} can be computed
using a system for structured optimization.
(The minimizer in \eqref{e-xk+half} exists and is unique,
so $x^{k+1/2}$ is well defined.
To see this we observe that the objective is finite for $x=x^k$, and the
function being minimized is strictly convex.)
We note that while $\hat f_k(x^{k+1/2})$ and $g(x^{k+1/2})$ are finite,
$f(x^{k+1/2})=\infty$ (and therefore also $h(x^{k+1/2})=\infty$) is possible.

The two quadratic terms in the objective in \eqref{e-xk+half}
can be combined to express the tentative update as
\begin{equation}\label{e-xk+half2}
x^{k+1/2} = \mathop{\rm argmin}_{x} \left( 
l_k(x) + g(x) + \frac{1}{2}(x-x^k)^T (H_k + \lambda_k I)(x-x^k)
\right),
\end{equation}
which shows that the trust penalty term can be interpreted as a
regularizer for the curvature estimate $H_k$.

In the DCP description of the problem \eqref{e-xk+half2},
using \eqref{PSD_low_rank} we express the last term as
\[
\frac{1}{2}(x-x^k)^T (H_k + \lambda_k I)(x-x^k) = 
\frac{1}{2}\|G_k^T(x - x^k)\|_2^2 + \frac{\lambda_k}{2} \|x - x^k\|_2^2.
\]
This keeps the problem \eqref{e-xk+half2} tractable 
when $r \ll n$ and $n$ is large.
In particular, there is no need to form the $n \times n$ matrix $H_k$.

\paragraph{Tentative update optimality condition.}
For future reference, we note that the optimality 
condition for the minimization in \eqref{e-xk+half2} that 
defines $x^{k+1/2}$ is
\begin{equation}\label{e-xk+half-opt}
0 \in \partial l_k(x^{k+1/2}) + \partial g(x^{k+1/2})
+ (H_k+\lambda_k I) (x^{k+1/2} - x^k).
\end{equation}
When $x^{k+1/2}$ is computed using a structured solver,
we can recover specific subgradients in
the subdifferentials $\partial l_k(x^{k+1/2})$ and $\partial g(x^{k+1/2})$
that satisfy \eqref{e-xk+half-opt}.
(How to do this is explained in \S\ref{appendix-q_k}.
We only need to find one of the two subgradients, since 
the three terms in \eqref{e-xk+half-opt} sum to zero,
and we know the third term.)
We will denote the specific subgradient in $\partial g(x^{k+1/2})$
in \eqref{e-xk+half-opt} as $q^{k+1}$.

When $l_k$ is the piecewise affine minorant \eqref{l_k_PWA},
its subdifferential $\partial l_k(x^{k+1/2})$ has the form
\begin{equation} \label{e-subdiff-lk}
\partial l_k(x^{k+1/2}) = {\mathop {\bf Co}} \{ \nabla f(x^i) \mid l_k(x^{k+1/2}) =
f(x^i)+ \nabla f(x^i)^T(x^{k+1/2}-x^i) \},
\end{equation}
the convex hull of the gradients associated with the 
active terms in maximum defining $l_k$.  
In the simplest case when $l_k$ is differentiable at 
$x^{k+1/2}$, {\it i.e.}, only one term is active, this reduces to
$\{ \nabla l_k(x^{k+1/2})\} = 
\{ \nabla f(x^i)\}$, where $i$ is the (unique) index for which
$l_k(x^{k+1/2}) = f(x^i)+ \nabla f(x^i)^T(x^{k-1/2}-x^i)$.

As explained in \S\ref{appendix-q_k}, we can compute 
a specific subgradient
$\sum_i \gamma_i \nabla f(x^i) \in \partial l_k(x^{k+1/2})$
for which \eqref{e-xk+half-opt} holds,
where $\gamma_i $ are nonnegative and sum to one,
and positive only for $i$ associated with 
active terms in the maximum that defines $l_k(x^{k+1/2})$.
We then have
\begin{equation}\label{e-qk}
q^{k+1} = 
- \sum_i \gamma_i \nabla f(x^i)
-(H_k+\lambda_k I) (x^{k+1/2}-x^k) \in \partial g(x^{k+1/2}).
\end{equation}

\subsection{Descent direction}\label{s-descent_direction}
If $x^k$ is a fixed point of the tentative update,
{\it i.e.}, $x^{k+1/2} =x^k$, 
then $x^k$ is optimal. To see this,
if $x^{k+1/2} =x^k$, from \eqref{e-xk+half-opt} we have
\begin{equation}\label{e-xk+half-opt-cond}
0\in \partial l_k(x^{k}) +  \partial g(x^{k}) = \nabla f(x^k) +  \partial g(x^{k}),
\end{equation}
so $x^k$ is optimal. 
From the first inclusion in \eqref{e-xk+half-opt-cond}, we can also
conclude that $x^k$ minimizes $l_k(x)+g(x)$, so
$L_k = l_k(x^k) + g(x^k) = f(x^k) + g(x^k)$,
{\it i.e.}, the lower bound in \eqref{e-Lk} is tight when $x^{k+1/2} = x^k$.

Otherwise, the tentative step
\begin{equation}\label{e-vk}
v^k = x^{k+1/2}-x^k
\end{equation} 
is a descent direction for $h$ at $x^k$, {\it i.e.},
for small enough $t>0$ we have $h(x^k + t v^k)<h(x^k)$.
That is, the directional derivative $h'(x^k; v^k)$ is negative.

To see this, we first observe that by \eqref{e-xk+half2},
\begin{equation}\label{desc_dir_step1}
l_k(x^{k+1/2}) + g(x^{k+1/2}) + \frac{1}{2}(v^k)^T(H_k+\lambda_k I)v^k
<
l_k(x^k) + g(x^{k}),
\end{equation}
since $x^{k+1/2}$ minimizes the left-hand side, and the right-hand side
is the same expression, evaluated at $x^k$.
We also have
\[
l_k(x^{k+1/2}) + g(x^{k+1/2}) \geq
l_k(x^{k}) + g(x^{k}) + (l_k+g)'(x^k; v^k).
\]
Combining these two inequalities we get
\[
(l_k+g)'(x^k; v^k) <
- \frac{1}{2}(v^k)^T(H_k+\lambda_k I)v^k.
\]
Finally, we observe that
\[
(l_k+g)'(x^k; v^k) = (f+g)'(x^k; v^k) = h'(x^k; v^k),
\]
since $l_k$ is differentiable at $x^k$, with $\nabla l_k(x^k) = \nabla f(x^k)$.
So we have
\begin{equation}\label{e-descent}
h'(x^k; v^k) < -\frac{1}{2} (v^k)^T (H_k+\lambda_k I)v^k ,
\end{equation}
which shows that $v^k$ is a descent direction for $h$ at $x^k$.

\subsection{Line search}\label{s-line-search}
The next iterate $x^{k+1}$ is found as 
\begin{equation}\label{e-xk+1}
x^{k+1} = x^k + t_k v^k = x^k + t_k (x^{k+1/2}-x^k),
\end{equation}
where $t_k \in (0,1]$ is the step size.
When $t_k=1$, we say the step is un-damped.
We will choose $t_k$ using a variation on 
a traditional Armijo-type line search \cite{armijo1966minimization} 
that avoids additional evaluations of $g$.

For $t \in [0,1]$ we define
\begin{equation}\label{e_phi_k}
\phi_k(t) = f(x^k + t v^k) + t g(x^{k+1/2}) + (1-t)g(x^k)).
\end{equation}
Since the second and third terms are the chord above $g$, 
we have, for $t \in [0,1]$,
\begin{equation}\label{e-chord}
\phi_k(t) \geq h(x^k+t v^k).
\end{equation}
Evidently $\phi_k(0)=h(x^k)$, and $\phi_k$ is differentiable, with
\[
\phi_k'(0) = \nabla f(x^k)^T v^k + g(x^{k+1/2})-g(x^k).
\]
Since $\nabla f(x^k) = \nabla l_k(x^k)$ and $l_k$ is convex, we get
\[
\nabla f(x^k)^T v^k = \nabla l_k(x^k)^T v^k \leq l_k(x^{k+1/2}) - l_k(x^k),
\]
so
we have
\[
\phi_k^\prime(0) \leq l_k(x^{k+1/2}) - l_k(x^k)+ g(x^{k+1/2})-g(x^k).
\]
Combining this with \eqref{desc_dir_step1}, we obtain
\begin{equation}\label{phi_prime_bound}
\phi_k^\prime(0) < - \frac{1}{2}(v^k)^T ( H_k + \lambda_kI) v^k.
\end{equation}
Thus for $t>0$ small,
\begin{equation}\label{e-phik-bnd}
\phi_k(t) = \phi_k(0) + \phi_k(0)^\prime t + o(t^2) 
< h(x^k)  -\frac{t}{2}(v^k)^T ( H_k + \lambda_kI) v^k + o(t^2).
\end{equation}

\paragraph{Step length.}
Let $\alpha, \beta \in (0,1)$.  We take $t_k = \beta^{j}$, 
where $j$ is the smallest nonnegative integer for which
\begin{equation}\label{armijo_step}
\phi_k(t_k) \leq h(x^k)
- \frac{\alpha t_k}{2} (v^k)^T (H_k+\lambda_k I)v^k
\end{equation}
holds.  
(The condition~\eqref{armijo_step} holds for some $j$ by \eqref{e-phik-bnd}.)
A nice feature of this line search is that it does not require any additional
evaluations of the function $g$ (which can be expensive),
since we already know $g(x^k)$ and $g(x^{k+1/2})$.
As has been noted by many authors, the choice of the line 
search parameters $\alpha$ and $\beta$ is not critical.
Traditional default values such as 
\begin{equation}\label{e-line-search-params}
\alpha=0.05, \qquad \beta=0.5
\end{equation}
work well in practice.

\subsection{Adjusting the trust parameter}
We have already observed that $\lambda_{k}$ is a regularizer for $H_{k}$.
A natural choice is to choose the regularizer parameter roughly 
proportional to $\tau_k= \mathop{\bf Tr} H_{k}/n$,
since $\tau_k I$ is the minimum Frobenius norm approximation 
of $H_{k}$ by a multiple of the identity. 
Thus we take
\begin{equation}\label{e-lamk}
\lambda_{k} = \mu_k \left(\tau_k + \tau_{\min} \right),
\end{equation}
where $\mu_k$ gives the trust parameter relative to $\tau_k$,
and $\tau_{\min}$ is a positive lower limit.

We update $\mu_k$ by decreasing it when the line search is undamped,
{\it i.e.}, $t_k=1$, and increasing it when the line search is damped,
{\it i.e.}, $t_k<1$.
We do this with
\begin{equation}\label{e-muk}
\mu_{k+1} = \left\{
\begin{array}{ll} 
\max\{ \gamma_\text{dec} \mu_k, \mu_{\min} \} & t_k=1 \\
\min\{ \gamma_\text{inc}\mu_k, \mu_{\max} \} & t_k<1,
\end{array} \right.
\end{equation}
where $\mu_\text{min}$ and $\mu_\text{max}$ are positive lower
and upper limits for $\mu_k$,
$\gamma_\text{dec} \in (0,1)$ is the factor by which we decrease
$\mu_k$,
and $\gamma_\text{inc} \in (1,\infty)$ is the factor by which we increase
$\mu_k$.
The values
\begin{equation}\label{e-lambda-update-params}
\tau_{\min} = 10^{-3}, \quad
\gamma_\text{dec} = 0.8, \quad 
\gamma_\text{inc} = 1.1, \quad
\mu_{\min} = 10^{-4}, \quad
\mu_{\max} = 10^5
\end{equation}
give good results for a wide range of problems.
We can take $\mu_0=1$.

We mention one initialization that is useful when
$f$ is twice differentiable and we have the ability to evaluate
$z \mapsto \nabla^2f(x^k) z$ ({\it i.e.}, Hessian-vector multiplication).
In this case we can replace $\tau_0$ with an estimate of 
$\mathop{\bf Tr} \nabla^2 f(x^0)/n$ obtained using the Hutchinson method \cite{hutchinson1989stochastic}.


\subsection{Stopping criteria}\label{s-stopping-crit}
We use two stopping criteria, one based on a gap 
between upper and lower bounds on the optimal value,
and the other based on an optimality condition residual.
The gap condition is simple:
\begin{equation}\label{e-gap-stop}
h(x^k) - L_k \leq \epsilon_\mathrm{abs}^\mathrm{gap} 
+ \epsilon_\mathrm{rel}^\mathrm{gap}|h(x^k)|,
\end{equation}
where $\epsilon_\mathrm{abs}^\mathrm{gap}$ and 
$\epsilon_\mathrm{rel}^\mathrm{gap}$ are positive
absolute and relative gap tolerances, respectively.
Evaluating $L_k$ can be almost as expensive as evaluating $x^{k+1/2}$,
but it is used only in the stopping criterion.  To reduce
this overhead, we evaluate $L_k$ only every ten iterations.
Reasonable values for the gap tolerances are
$\epsilon_\mathrm{abs}^\mathrm{gap} = 10^{-4}$ and
$\epsilon_\mathrm{rel}^\mathrm{gap} = 10^{-3}$.

The residual based stopping criterion is tested whenever we
take an undamped step, {\it i.e.}, $t_k=1$.  In this case $x^{k+1}=x^{k+1/2}$,
and we obtain $q^{k+1} \in \partial g(x^{k+1})$ in \eqref{e-qk},
so $\nabla f(x^{k+1}) + q^{k+1}$ is a residual for the optimality condition
\eqref{e-opt-cond}.
The stopping criterion is
\begin{equation}\label{stopping_criterion_3}
\frac{1}{\sqrt{n}}\|\nabla f(x^{k+1}) + q^{k+1}\|_2 \leq
\epsilon_\mathrm{abs}^\mathrm{res} 
+ \epsilon_\mathrm{rel}^\mathrm{res} 
\left(\frac{1}{\sqrt{n}}\|\nabla f(x^{k+1})\|_2 + \frac{1}{\sqrt{n}}\|q^{k+1}\|_2 \right),
\end{equation}
where 
$\epsilon_\mathrm{abs}^\mathrm{res}$ and 
$\epsilon_\mathrm{rel}^\mathrm{res}$ are
relative and absolute residual tolerances.
(Dividing the norm expressions above by $\sqrt n$ gives the
root mean square or RMS values of the argument.)
Reasonable values for these parameters are
$\epsilon_\mathrm{abs}^\mathrm{res}=10^{-4}$ and
$\epsilon_\mathrm{rel}^\mathrm{res}=10^{-3}$.

\clearpage
\subsection{Algorithm summary}
We summarize OSMM in algorithm \ref{alg:1}.
\begin{algdesc}\label{alg:1}
\emph{Oracle-structured minimization method.}
\begin{tabbing}
{\bf given} an initial point $x^0\in\Omega$.\\
{\bf for} $k=0, 1, \ldots, k_{\max}$\\
\qquad \= 1.\ \emph{Form a surrogate objective.}
Form $l_k$ and $H_k$.\\
\> 2.\ \emph{Tentative step.} Compute $x^{k+1/2}$ by \eqref{e-xk+half}.\\
\> 3.\ \emph{Line search and update.} Set line search 
step size $t_k$ by \eqref{armijo_step}
and $x^{k+1}$ by \eqref{e-xk+1}.\\
\> 4.\ \emph{Compute lower bound.}  If $k$ is a multiple of $10$,
evaluate $L_k$.\\
\> 5.\ \emph{Check stopping criterion.} Quit if \eqref{e-gap-stop} or 
\eqref{stopping_criterion_3} holds.\\
\> 6.\ \emph{Update trust penalty parameter.} Update $\lambda_{k+1}$ by \eqref{e-lamk}
and \eqref{e-muk}.
\end{tabbing}
\end{algdesc}

The algorithm parameters in OSMM are the memory of the minorant $M$, the rank of 
the curvature estimate $r$, the line search parameters given in
\eqref{e-line-search-params}, the $\lambda_k$ update parameters
given on \eqref{e-lambda-update-params}, and the relative and absolute
gap and residual tolerances, given in \S\ref{s-stopping-crit}.
The practical performance of OSMM is not particularly sensitive to 
the choice of these parameters; our implementation uses as default values
the ones described above, with memory $M=20$ and rank $r=20$.


\subsection{Implementation}
We have implemented OSMM in an open-source Python package, available at
\begin{center}
\url{https://github.com/cvxgrp/osmm}.
\end{center}
The user supplies a PyTorch description of the oracle function $f$, 
a CVXPY description of the structured function
$g$,  and an initial point $x^0\in\Omega$. 
The package invokes PyTorch to evaluate $f$ and its gradient $\nabla f$.  
The convex model $\hat f_k$ is then formed and handed off to CVXPY to 
efficiently compute the next tentative iterate.

\subsection{Related work}\label{s-related-work} There is a lot of prior work
related to the method proposed in this paper.  The work on variable metric
bundle
methods~\cite{lemarechal1978nonsmooth,fukushima1984descent,kiwiel1990proximity,schramm1992version,lemarechal1995new,mifflin1996quasi,lemarechal1997variable,mifflin1998quasi,lukvsan1998bundle,kiwiel2000efficiency,teo2010bundle,yu2010quasi,noll2013bundle,Bagirov2014,van2016strongly,de2016doubly}
and (inexact) proximal Newton-type
methods~\cite{levitin1966constrained,schmidt2009optimizing,sra2012optimization,becker2012quasi,lee2014proximal,scheinberg2016practical,li2017inexact,ghanbari2018proximal,yue2019family,becker2019quasi,lee2019inexact,asi2019modeling,mordukhovich2020globally}
are probably the most closely related, though our method differs
in key ways arising from our assumed access methods.  Proximal Newton methods
are similar in philosophy to our approach, as both types of methods really shine
when the structured part of the objective can be minimized efficiently.
However, on a purely technical level, proximal Newton methods generally form
the (usual) first-order model of the oracle part of the objective, without 
a bundle term, 
which is of course different than the model we use,
and can affect the convergence speed.

A few variable metric bundle methods have been proposed over the years, but
these methods also differ from ours in subtle (but important) ways.  In
general, these methods are not geared towards solving composite minimization
problems.  
Additionally, our line search method is different, 
as ours is carefully designed to satisfy the access conditions.

Finally, there has been a surge of interest in scalable quasi-Newton
methods~\cite{erdogdu2015convergence,pmlr-v70-ye17a,gower2018tracking,huang2020span} recently.
The focus has been on developing low-rank, regularized, and sub-sampled
Newton-type methods.  Of these, the family of low-rank approximations to the
Hessian are most related to the curvature estimate used in our method.

\section{Convergence}
In this section we show the convergence of OSMM,
and give a numerical example to illustrate its typical performance,
as the two most important algorithm parameters, memory $M$ and rank $r$, vary.
(More examples will be given in \S\ref{s-experiments}.)

\subsection{Convergence}\label{s-convergence}
We demonstrate two types of convergence.  First, we show the iterates $x^k$ converge asymptotically in a certain sense, {\it i.e.}, the sequence of tentative steps $v^k \to 0$ as $k \to \infty$.  (From \eqref{e-xk+half-opt}, this implies that when the iterates $x^k$ do converge, the convergence point must be optimal.)  Second, we show the objective values generated by OSMM converge to the optimal value.  
The results additionally require the (standard) assumption that $\nabla f$ be Lipschitz continuous.

\paragraph{Convergence of iterates.}
To see that $v^k \to 0$, we observe from \eqref{e-chord}~and \eqref{armijo_step} that
\begin{equation}\label{e-quad-form-bound}
h(x^{k+1}) - h(x^k) \leq - \frac{\alpha t_k}{2} (v^k)^T (H_k + \lambda_k I) v^k \leq 0.
\end{equation}
Since $h(x^k)$ is decreasing and bounded below by $h^\star$, it converges,
which implies that 
the left-hand side of \eqref{e-quad-form-bound} converges to zero as
$k \to \infty$.  This in turn implies that
\[
t_k (v^k)^T (H_k + \lambda_k I) v^k \to 0 \quad \mathrm{as} \quad k \to \infty.
\]
By construction $\lambda_k$ is lower bounded by $\tau_{\min} \mu_{\min}$,
and $t_k$ is also bounded away from $0$, 
as will be shown later in \eqref{lower-bound-tk},
so we get that $v^k\to 0$, as claimed.


\paragraph{Convergence of objective values.}
We require the basic fact that there exists a subsequence of iterates $(k_\ell)$ accepting unit step length, {\it i.e.}, that $t_{k_\ell} = 1$ for all $\ell$.  We give a proof of this fact in \S\ref{s-tk-one-io}, assuming that $\nabla f$ is Lipschitz continuous.

Now, by convexity,
\[
h(x^{k_\ell+1}) - h^\star \leq (\nabla f(x^{k_\ell+1}) + q^{k_\ell+1})^T (x^{k_\ell+1} - x^\star),
\]
for some $q^{k_\ell} \in \partial g(x^{k_\ell})$.  Adding and subtracting any $l_{k_\ell}^{k_\ell+1} \in \partial l_{k}(x^{k+1})$, we get
\begin{equation}\label{e-func-val-descent0}
h(x^{k_\ell+1}) - h^\star \leq \big( ( \nabla f(x^{k_\ell+1}) - l_{k_\ell}^{k_\ell+1} ) + ( l_{k_\ell}^{k_\ell+1} + q^{k_\ell+1} ) \big)^T (x^{k_\ell+1} - x^\star).
\end{equation}

From our assumptions on $h$, and because OSMM is a descent method, we see that $\| x^{k_\ell+1} - x^\star \|_2$ is bounded.  Moreover, because $H_k$ and $\lambda_k$ are uniformly bounded, we get from \eqref{e-xk+half-opt} that
\[
\| l_{k_\ell}^{k_\ell+1} + q^{k_\ell+1} \|_2 \leq \| H_{k_\ell} 
+ \lambda_{k_\ell} I \|_2 \| v^{k_\ell} \|_2 \lesssim \| v^{k_\ell} \|_2,
\]
where we write $a \lesssim b$ to mean that $a$ and $b$ satisfy $a \leq C b$, for some $C > 0$.  Finally, it can be shown (see \S\ref{s-minorant-accuracy} for details) that the limited memory piecewise affine minorant satisfies
\begin{equation}\label{e-minorant-accuracy}
\| \nabla f(x^{k_\ell+1}) - l_{k_\ell}^{k_\ell+1} \|_2 \lesssim \max_{j=\max\{0,\ldots,k_\ell-M+1\},\ldots,k_\ell} \| v^j \|_2.
\end{equation}
Putting these together with \eqref{e-func-val-descent0}, we obtain that
\[
h(x^{k_\ell+1}) - h^\star \lesssim \max_{j=\max\{0,\ldots,k_\ell-M+1\},\ldots,k} \| v^j \|_2.
\]
Earlier we showed that $v^k \to 0$, so $h(x^{k_\ell+1}) \to h^\star$ as $\ell \to \infty$.  Because the sequence of objective values $h(x^k)$ is convergent, we get $h(x^k) \to h^\star$, as claimed.

\subsection{A numerical example}\label{s-numerical-cvg-example}
Next we investigate the convergence of OSMM numerically.  
Here and throughout this paper, 
we use a Tesla V100-SXM2-32GB-LS GPU with 32 gigabytes of memory 
to evaluate $f$ and $\nabla f$ via PyTorch, 
and an Intel Xeon E5-2698 v4 2.20 GHz CPU 
to compute the tentative iterates via CVXPY and for any baselines.

\paragraph{Problem formulation.}  
We consider an instance of the Kelly gambling problem~\cite[\S 4]{BoydVand04}
\cite{kelly2011new,busseti2016risk},
\begin{equation}\label{e-kelly}
\begin{array}{ll}
\mbox{maximize} & \sum_{i=1}^N \pi_i \log ( r_i^T x ) \\
\mbox{subject to} & x \geq 0, \quad \mathbf 1^T x = 1,
\end{array}
\end{equation}
where $x \in {\mbox{\bf R}}^n$ is the variable, and $r_i \in {\mbox{\bf R}}^n_+$,
$i=1,\ldots,N$, $\pi \in {\mbox{\bf R}}^N_+$ are the problem data, with $\sum_{i=1}^N \pi_i =1$.
Here $x_j$ is the fraction of our wealth we put on bet $j$,
$(r_i)_j$ gives the return of bet $j$ with outcome $i$, and 
$\pi_i$ is the probability of outcome $i$.
To put the Kelly gambling problem into our oracle-structured form, 
we take $f$ to be the objective in~\eqref{e-kelly}, 
and $g$ to be the indicator function of the constraints,
in this case the probability simplex.

\paragraph{Problem instance.}
Our problem instance has $n=1{,}000$ bets and $N = 1{,}000{,}000$ 
possible outcomes.
Note that simply storing the returns $r_i$,
which are dense in this instance, requires roughly four gigabytes, 
causing most existing solvers to struggle, as we will soon see.  
More detail about the Kelly gambling problem, and how the data were chosen,
can be found in~\S\ref{subsec:experiment-kelly}.
This problem instance requires 4,800 seconds to solve using 
MOSEK~\cite{mosek}, a high performance commercial solver, with accuracy set to
its lower value.  

\paragraph{Results.}  We solve~\eqref{e-kelly} by OSMM with rank $r \in
\{0,20,50\}$ and memory $M \in \{1,20,50\}$, a total of nine
choices of algorithm parameters.
The choices $r=0$ and $M=1$ correspond to using
no estimate of curvature, and no memory, respectively.
OSMM is run for 100 iterations for each combination of $r$ and $M$; $x^{k+1/2}$ and $L_k$ are computed using ECOS.
The results are shown
in figure~\ref{fig:convergence_example} and
table~\ref{tab:example_time_details}.  
The optimal objective value is -0.079.
Figure~\ref{fig:convergence_example}
shows the convergence of OSMM, in terms of iterations (the left column)
and elapsed time (the right column).
The top row shows the true suboptimality $h(x^k) - h^\star$, which of course
we do not know when the algorithm is running, since we do not know $h^\star$.
The middle row shows the gap $h(x^k)- L_k$ 
(which we do know as the algorithm runs).   
The bottom row shows the RMS value of 
the optimality condition residual, {\it i.e.}, the left-hand side
of~\eqref{stopping_criterion_3}.

We see that rank values $r=20$ and $r=50$, and memory values 
$M=1$ or $M=20$ yield the fastest convergence.  
We have observed similar results in many other problems.
The default choices in OSMM are $r=20$, $M=20$.

\begin{figure*}
\centering
\includegraphics{figures/large_scale_kelly_subopt_left.pdf}
\includegraphics{figures/large_scale_kelly_subopt_right.pdf}
\includegraphics{figures/large_scale_kelly_gap_left.pdf}
\includegraphics{figures/large_scale_kelly_gap_right.pdf}
\includegraphics{figures/large_scale_kelly_rmsres_left.pdf}
\includegraphics{figures/large_scale_kelly_rmsres_right.pdf}
\caption{Suboptimality (top row), and RMS residual (bottom row), versus
iterations (left column) and run-time in seconds (right column) for
OSMM on the Kelly gambling problem for the nine combinations of 
rank $r$ and memory $M$.}
\label{fig:convergence_example}
\end{figure*}

More detail is given in table~\ref{tab:example_time_details}.
The total solve time and iterations are based on achieving a high accuracy
of $10^{-6}$, which is far more accurate than would be needed in practice.
For our nine choices of algorithm parameters, the total OSMM time 
ranges from around 0.8 to 16 seconds, substantially faster than using MOSEK
to solve the problem.

\begin{table}
  \centering
  \begin{tabular}[]{ll lll ll}
    \toprule
    \multicolumn{2}{l}{Rank and memory} &
    \multicolumn{2}{l}{Compute times (sec.)}&
\multicolumn{3}{l}{Total solve time and iterations} \\
    \cmidrule(lr){1-2}
    \cmidrule(lr){3-4}
	\cmidrule(lr){5-7}
    $r$ & $M$ & $x^{k+1/2}$ & $L_k$ & Time (sec.) & Iterations & $f$ evals./iter.\\
    \midrule
    \midrule
    0 & 1    & 0.029 & 0.012 & 5.1 & 60 & 5.6\\
    0 & 20 & 0.062 & 0.077 & 8.6 & 52 & 4.9\\
    0 & 50 & 0.14 & 0.39 & 16  & 44 & 4.8\\
    \cmidrule(lr){1-2}
    \cmidrule(lr){3-5}
	\cmidrule(lr){6-7}
	20 & 1    & 0.058  &  0.013 & 0.81 & 11 & 1.9\\
	20 & 20 & 0.11  & 0.070  & 1.7 & 11 & 1.9\\
	20 & 50 & 0.20  & 0.34  & 3.8 & 11 & 1.9\\
    \cmidrule(lr){1-2}
    \cmidrule(lr){3-5}
	\cmidrule(lr){6-7}
	50 & 1    & 0.061  & 0.014& 0.85 & 11& 1.9\\
	50 & 20 & 0.11  & 0.093 & 1.7  & 11  & 1.9\\
	50 & 50 & 0.22  & 0.40 & 3.9 & 11  & 1.9\\
    \bottomrule
  \end{tabular}
  \caption{Times to evaluate tentative update $x^{k+1/2}$ and lower bound $L_k$,
total solve time and iterations to reach $10^{-6}$ accuracy, and average number of evaluations of $f$
per iteration, for OSMM with the nine combinations of rank $r$ and memory $M$.}
\label{tab:example_time_details}
\end{table}


\clearpage
\section{Generic applications}\label{s-applications}
In this section we describe some generic applications that reduce to our 
specific oracle-structured composite function minimization problem.
We focus on the function $f$, which should be differentiable, but not so
simple that it can be handled directly in a structured optimization
system.  We will see that this generally occurs when there is a lot of 
data required to specify $f$.

\subsection{Stochastic programming} 
\paragraph{Sample average approximation for stochastic programming.}
We start with the problem of minimizing $\mathop{\bf E{}} F(x,\omega) + g(x)$,
where $F$ is convex in $x$ for each value of the random variable $\omega$.
We will approximate the first objective term using a sample average.
We generate $N$ independent samples $\omega_1, \ldots, \omega_N$,
and take
\begin{equation}\label{e-emp-mean}
f(x) = \frac{1}{N} \sum_{i=1}^N F(x,\omega_i).
\end{equation}
As a variation we can use importance sampling to get a lower variance estimate 
of $\mathop{\bf E{}} F(x,\omega)$. To do this we generate the samples from a 
proposal distribution with density $q$, and form the sample average estimate
\begin{equation}\label{e-imp-sampling}
f(x) = \frac{1}{N} \sum_{i=1}^N \frac{p(\omega_i)}{q(\omega_i)}F(x,\omega_i),
\end{equation}
where $p$ is the density of $\omega$.

In both cases we can take $N$ to be quite large, 
since we will only need to evaluate $f$
and its gradient, and current systems for this are very efficient.
For example, the gradients $\nabla F(x,\omega_i)$ can be computed in
parallel.  As a practical matter, this happens automatically, with no
or very little directive from the user who specifies $f$.

\paragraph{Validation and stopping criterion tolerance.}
The functions $f$ given in \eqref{e-emp-mean} and \eqref{e-imp-sampling}
are only approximations of the true objective $\mathop{\bf E{}} F(x,\omega)$,
though we hope they are good approximations when we take $N$ large,
as we can with OSMM.
To understand how accurate the sample average \eqref{e-emp-mean} is,
we generate another set of independent samples 
$\tilde \omega_1,\ldots, \tilde \omega_N$
and define the validation function as
\[
f^\mathrm{val}(x) = \frac{1}{N} \sum_{i=1}^N 
F(x, \tilde \omega_i)
\]
(and similarly if we use importance sampling).
The magnitude of the difference $|f^\mathrm{val}(x)-f(x)|$ gives us a 
rough idea of the accuracy in approximating $\mathop{\bf E{}} F(x,\omega)$.
(Better estimates of the accuracy can be obtained by repeating
this multiple times, but we are interested in only a crude estimate.)

Solving the oracle-structured problem to an accuracy substantially better than
the Monte Carlo sampling accuracy does not make sense in practice.
This justifies replacing the absolute gap tolerance
$\epsilon_\mathrm{abs}^\mathrm{gap}$ with the maximum of 
a fixed absolute tolerance and 
the sampling error $|f^\mathrm{val}(x^k) - f(x^k)|$.
(We can evaluate the sampling error whenever we evaluate the gap,
{\it i.e.}, every 10 iterations.)
Roughly speaking, we stop when we know we have solved the problem
to an accuracy that is better than our approximation.

\subsection{Utility maximization}
An important special case of stochastic optimization is
utility maximization, where we seek to maximize
\begin{equation}\label{e-util-max}
\mathop{\bf E{}} U(H(x,\omega)) - g(x),
\end{equation}
where $U:{\mbox{\bf R}} \to {\mbox{\bf R}}$ is a concave increasing utility function,
and $H$ is concave in $x$ for each $\omega$.
The first term, the expected utility, is concave.
This is equivalent to the stochastic programming problem of 
minimizing
\[
\mathop{\bf E{}} (-U(H(x,\omega)) + g(x),
\]
which is stochastic programming with $F= -U\circ H$, which is convex
in $x$.  We can replace the expectation with a sample average
using \eqref{e-emp-mean}, or an importance sampling sample average 
using \eqref{e-imp-sampling}.
Utility maximization is a common method for handling the variance or 
uncertainty in a stochastic objective; it introduces risk aversion.

\subsection{Conditional and entropic value-at-risk programming}\label{s-EVaR}
Another method to introduce risk aversion into a stochastic optimization
problem is to mimimize value-at-risk (VaR), or a specific quantile of
$F(x,\omega)$, where
$F$ is convex in $x$ for each value of the random variable $\omega$.
The value-at-risk is defined as
\[
\mathop{\bf VaR}(F(x,\omega);\eta) = \inf \{ \gamma \mid \mathop{\bf Prob}(F(x,\omega) \leq \gamma) \geq \eta \},
\]
where $\eta \in (0,1)$ is a given quantile.
With an additional structured convex function $g$ in the objective, we obtain
the VaR problem
\begin{align}\label{prob-var}
\begin{array}{ll}
\mbox{minimize} & \mathop{\bf VaR}(F(x,\omega); \eta)) + g(x),
\end{array}
\end{align}
which, roughly speaking, is the problem of minimizing the $\eta$-quantile of
the random variable $F(x,\omega)+g(x)$.
Aside from a few special cases, this problem is not convex.  
We proceed by replacing the nonconvex VaR term with a convex upper bound on VaR,
such as conditional value-at-risk (CVaR)
\cite{CVAR_uryasev_rockafellar,rockafellar2002conditional} 
or entropic value-at-risk (EVaR) \cite{ahmadi2011information}.
Beyond resulting in tractable convex problems, these upper bounds possess a number 
of nice properties, such as being coherent risk measures, 
which VaR is not; 
see \cite{rockafellar2002conditional,ahmadi2011information} for a discussion.

CVaR is given by
\begin{equation}\label{e-cvar}
\mathop{\bf CVaR}(F(x,\omega);\eta) = \inf_{\alpha \in \mathbf{R}} \Bigg\{ \frac{\mathop{\bf E{}} (F(x,\omega) - \alpha)_+}{1-\eta} + \alpha \Bigg\},
\end{equation}
where $(z)_+ = \max\{z,0\}$, and EVaR is given by
\begin{equation}\label{e-evar}
\mathop{\bf EVaR}(F(x,\omega);\eta) = \inf_{\alpha > 0} \Bigg\{ \alpha \log \Bigg( \frac{\mathop{\bf E{}} \exp(F(x,\omega)/\alpha)}{1-\eta} \Bigg) \Bigg\}.
\end{equation}
Both of these are convex functions of $x$.
We have, for any $x$, 
\[
\mathop{\bf VaR}(F(x,\omega);\eta) \leq \mathop{\bf CVaR}(F(x,\omega);\eta) \leq \mathop{\bf EVaR}(F(x,\omega); \eta)
\]
(see, {\it e.g.}, \cite{ahmadi2011information}).
With CVaR, we obtain the convex problem
\[
\begin{array}{ll}
\mbox{minimize} & \mathop{\bf CVaR}(F(x,\omega); \eta)) + g(x),
\end{array}
\]
and similarly for EVaR.

We now show how the CVaR and EVaR problems can be approximated as oracle-structured
problems.
We start with CVaR.
We generate independent samples $\omega_i$, $i=1, \ldots, N$,
and replace the expectation with the empirical mean,
\[
f(x,\alpha) = \frac{1}{N} \sum_{i=1}^N
\frac{(F(x,\omega_i) - \alpha)_+}{1-\eta} + \alpha,
\]
with variables $x$ and $\alpha$.   This function is jointly convex in $x$ and $\alpha$;
minimizing over $\alpha$ gives $\mathop{\bf CVaR}(F(x,\omega);\eta)$ for the empirical 
distribution.
We adjoin $\alpha$ to $x$ to obtain a problem in oracle-structured form, {\it i.e.},
minimizing $f(x,\alpha)+g(x)$, over $x$ and $\alpha$.
(That is, we take $(x,\alpha)$ as what we call $x$ in our general form.)
While $f$ is not differentiable in $(x,\alpha)$, we have observed that 
our method still works very well.

For EVaR, we obtain
\[
f(x,\alpha) = \alpha \log \Bigg( \frac{1}{N} \sum_{i=1}^N
\frac{\exp(F(x,\omega_i)/\alpha)}{1-\eta} \Bigg),
\]
which is jointly convex in $x$ and $\alpha$.
(To see convexity, we observe that $f$ is the perspective function
of the log-sum-exp function; see \cite[\S 3.2.6]{BoydVand04}.)
As with CVaR,
minimizing over $\alpha$ yields $\mathop{\bf EVaR}(F(x,\omega);\eta)$ for the empirical
distribution.
Unlike our approximation with CVaR, this function is differentiable.

\subsection{Generic exponential family density fitting} \label{s-exp-fam-fit}
We consider fitting an exponential family of densities, given by
\begin{equation}\label{e-exp-fam}
p_\theta(z) = e^{- \left(\phi(z)^T \theta + A(\theta) \right)},
\end{equation}
to samples $z_1,\ldots, z_m\in\mathcal{S}$.  Here, $\mathcal S$ is the support of the density, $\phi : \mathcal{S} \to {\mbox{\bf R}}^n$ is the sufficient statistic, and $\theta \in \Theta \subseteq {\mbox{\bf R}}^n$ is the canonical parameter (to be fitted).  The density normalizes via the log-partition (or cumulant generating) function
\[
A (\theta) = \log \int_\mathcal{S} e^{-\phi(z)^T \theta} \; dz.
\] 
We assume $\Theta$ is convex, and additionally that it only contains parameters for which the log-partition function is finite.

The negative log-likelihood, given samples $z_1,\ldots,z_m$, is
\[
\sum_{i=1}^m \log p_\theta(z_i) = c^T \theta + m A(\theta),
\]
where $c = \sum_{i=1}^m \phi(z_i)$.  So maximum likelihood estimation of $\theta$ corresponds to solving the density fitting problem
\begin{align}\label{exp_density_prob_0}
\begin{array}{ll}
\mbox{minimize} & \frac{1}{m}c^T \theta + A(\theta) + \lambda r(\theta)\\
\mbox{subject to} 
& \theta\in\Theta,
\end{array}
\end{align}
with variable $\theta$.  We can include a (potentially nonsmooth) convex regularization term $\lambda r(\theta)$ in the objective, where $\lambda \geq 0$ is the regularization strength, and $r : \Theta \to {\mbox{\bf R}}$ is the regularizer.  Since the log-partition function is convex~\cite{wainwright2008graphical}, the density fitting problem~\eqref{exp_density_prob_0} is also convex.

The log-partition function is generally intractable except for a few special cases, so we replace the integral in $A(\theta)$ with a finite sum using importance sampling, {\it i.e.},
\[
\int_\mathcal{S} e^{-\phi(z)^T \theta} \; dz \approx
\frac{1}{N} \sum_{i=1}^N \frac{1(\omega_i\in\mathcal{S})}{q(\omega_i)}e^{-\phi(\omega_i)^T \theta},
\]
where $\omega_i, \; i=1,\ldots,N$, are independent draws from the proposal distribution $q$.  
The number of samples $N$ can be very large, especially when the number of dimensions $n$ is moderate.  
When $\mathcal S$ is bounded and its dimension is small, we can simply use a Riemann sum, so that the samples $\omega_i$ are lattice points in $\mathcal S$ and we have $q(\omega_i) = 1/|\mathcal{S}|$.
The problem~\eqref{exp_density_prob_0} is clearly in oracle-structured form, once we take $A(\theta)$ (with its Monte Carlo approximation) to be $f$, and the rest of the objective plus the indicator of $\Theta$ to be $g$.

A number of interesting regularizers are possible.  The squared $\ell_2$ norm, {\it i.e.}, $r(\theta) = \| \theta \|_2^2$, is of course a natural choice.  
When $\mathcal{S}$ is bounded, a different option is to use the squared $L_2$ norm of the gradient of the log-density
\[
r(\theta) = \int_\mathcal{S} \|\nabla \log p_\theta(z)\|_2^2 \; dz.
\]
This regularizer enforces a certain kind of smoothness: the regularized density $p_\theta$ tends to the uniform distribution on $\mathcal S$, as the regularization strength $\lambda$ grows.  Finally, observe that we can write
\[
\int_\mathcal{S} \|\nabla \log p_\theta(z)\|_2^2 \; dz =
\int_\mathcal{S} \|D\phi(z)^T \theta\|_2^2 \; dz
= \theta^T \left(\int_\mathcal{S} D\phi(z) D\phi(z)^T \; dz \right) \theta,
\]
where $D\phi$ is the Jacobian of sufficient statistic $\phi$.  We can replace the integral with a finite sum again, and obtain the regularizer
\begin{equation*}
r(\theta) \approx \theta^T \left(\frac{1}{N} \sum_{i=1}^N 
\frac{1(\omega_i\in\mathcal{S})}{q(\omega_i)} D\phi(\omega_i) D\phi(\omega_i)^T \right) \theta.
\end{equation*}


\section{Numerical examples}\label{s-experiments} In this section we
demonstrate the performance of OSMM through several numerical examples, 
all taken from the
generic applications described in the previous section.
We start by showing results for two different portfolio selection problems,
the Kelly gambling example (shown earlier in \S\ref{s-numerical-cvg-example}),
and minimizing CVaR.  We then show
a density estimation example.  After that we present results for a
supply chain optimization problem with entropic value-at-risk minimization.  
All of these examples are structured stochastic optimization problems, and we
use simple sample averages to approximate expectations.
OSMM is designed to handle the case when $f$ is complex, 
which in the case of sample averages means $N$ is large.
We will see that when $N$ is small, OSMM is actually slower than just 
solving the problem directly using a structured solver; when $N$ is large,
it is much faster (and in many cases, directly using a structured solver
fails).

We report the time needed for OSMM to reach high accuracy, {\it i.e.}, $h(x^k) - h^\star \leq 10^{-6}$.  
(This makes a fairer comparison with MOSEK, SCS~\cite{ocpb:16,scs}, and ECOS~\cite{Domahidi2013ecos}.)
We also indicate when practical accuracy is reached, using our default 
parameters, and the sampling accuracy.
We use the default parameters in OSMM,
and use ECOS as the solver in CVXPY to compute the tentative 
iterate $x^{k+1/2}$.  
We do not perform any parameter tuning for our method.

\subsection{Kelly gambling}\label{subsec:experiment-kelly}

\paragraph{Problem formulation.} In the Kelly gambling problem there are $n$
bets we can wager on, and $N$ possible outcomes, with probabilities
$\pi_i$, $i=1, \ldots, N$.
The bet returns are given by $r_i\in {\mbox{\bf R}}_+^n$, where $(r_i)_j$ is
the return, {\it i.e.}, the amount you win for each dollar you put on bet $j$
when outcome $i$ occurs.
We are to choose $x\in {\mbox{\bf R}}^n$, with $\mathbf 1^T x= 1$, where $x_j$ is
the fraction of our wealth we place on bet $j$.
We seek to maximize the average log return, which maximizes long-term 
wealth growth if we repeatedly bet.
This leads to the (convex) optimization problem
\begin{equation*}
\begin{array}{ll}
\mbox{maximize} & \sum_{i=1}^N \pi_i \log (r_i^T x) \\
\mbox{subject to} & x \geq 0, \quad \mathbf 1^T x = 1,
\end{array}
\end{equation*}
with variable $x$.

\paragraph{Problem instance.}    
We consider $n=200$ bets. 
The probabilities of the outcomes $\pi_i$ are independent draws 
from a uniform distribution on $[0,1]$, normalized to sum to one.
The returns of the bets in each outcome are independently drawn 
from a log normal distribution, 
{\it i.e.}, $\exp(z)$, $z \sim \mathcal{N}(0, 1)$, 
and then scaled so the expected return
of bet $j$ is $\bar r_j$, {\it i.e.}, $\sum_{i=1}^N \pi_i (r_i)_j = \bar r_j$,
where $\bar r_j$ is drawn from a uniform distribution
on $[0.9, 1.1]$.
For our problem instance, the solution has 57
nonzero entries ranging from 0.001 to 0.04, 
and the optimal mean log return is 0.057.

\paragraph{Results.}  
The run-times for OSMM, MOSEK, SCS, 
and ECOS are shown in table~\ref{tab:runtimes:kelly}.  
When the number of Monte Carlo samples is small, {\it e.g.}, $N = 1{,}000$, 
MOSEK performs the fastest and takes less than a second to attain high accuracy.
OSMM also takes less than one second.
ECOS takes about two seconds,
and SCS takes roughly six seconds.  
However, when $N$ is larger ({\it i.e.}, $N = 10{,}000$, $100{,}000$, 
or $1{,}000{,}000$), OSMM is the fastest method, always taking less than a second to attain the required $10^{-6}$ optimality gap.
When $N=10{,}000$, MOSEK is still competitive, but SCS and ECOS 
are two to four orders of magnitude slower than OSMM.
When $N=100{,}000$, MOSEK and ECOS are two to three 
orders of magnitude slower, and when $N=1{,}000{,}000$, 
MOSEK and ECOS are three to five orders of magnitude slower.
SCS fails for both the two larger values of $N$.
These findings suggest that OSMM is useful when the number of samples is large, 
as it exhibits good scaling with $N$.

OSMM spends 0.0013 and 0.0024 seconds to evaluate $f$ 
and its gradient $\nabla f$, respectively, when $N=1{,}000{,}000$.  
Computing the tentative iterate $x^{k+1/2}$ and the lower bound 
$L_k$ from~\eqref{e-Lk} takes 0.033 and 0.014 seconds, respectively.  
The line search also turns out to be quite efficient here, 
as $f$ is evaluated twice during the line search on average.

\begin{table}
  \centering
  \begin{tabular}[]{l llll}
    \toprule
	$N$ & OSMM & MOSEK & SCS & ECOS \\
    \midrule
    \midrule
    1{,}000   & 0.76 & 0.58&6.4& 2.2   \\
    10{,}000 & 0.64 & 4.5 & 2{,}100 & 50 \\
	100{,}000 & 0.62 & 62 & --- & 910 \\
	1{,}000{,}000 & 0.64 & 890 & --- & 20{,}000 \\
	\bottomrule
	\end{tabular}
	\caption{Solve times in seconds on the Kelly gambling problem.  
		A dash (``---'') means the solver failed, either for numerical reasons,
		or because it did not reach the required $10^{-6}$ suboptimality
		in 24 hours.}\label{tab:runtimes:kelly}
\end{table}

Figure~\ref{fig:kelly_convergence} shows the convergence of OSMM 
with $N = 1{,}000{,}000$.  
In the top panel, we can see that practical accuracy is reached after
14 iterations, and high accuracy is reached after 16 iterations,
as shown by the dotted black and green lines, respectively.

\begin{figure*}
\centering
\includegraphics{figures/kelly_cvg_top.pdf}
\includegraphics{figures/kelly_cvg_bottom.pdf}
\caption{Suboptimality, gap (top row), 
	and RMS residual (bottom row) on the Kelly gambling problem with $N = 1{,}000{,}000$.
	The dotted green and black lines show when high accuracy and 
	practical accuracy are reached, respectively.}
\label{fig:kelly_convergence}
\end{figure*}


\clearpage
\subsection{CVaR portfolio optimization with derivatives}
\label{subsec:experiment-cvar}
\paragraph{Problem formulation.}
We consider making investments in $m$ stocks, and call and put
derivatives on them, with various strike prices.
Our investment will be for one period (of, say, one month).
We let $\omega \in {\mbox{\bf R}}_{++}^m$ denote the (fractional) change in prices
of the $m$ underlying stocks, which we model as random with a log-normal
distribution, {\it i.e.}, $\log \omega \sim \mathcal N(\mu,\Sigma)$,
where the $\log$ is elementwise.
For simplicity we will assume there is one call and one put option available
for each stock.
Let $p_\text{c} \in{\mbox{\bf R}}_{++}^m$ and $s_\text{c}\in{\mbox{\bf R}}_{++}^m$ 
be the call option prices (premia) and strike prices, normalized by the current 
stock price, so, for example, $(s_\text{c})_i = 1.15$ means the strike price
is $1.15$ times the current stock price.
Let $p_\text{p}\in{\mbox{\bf R}}^m_{++}$ and $s_\text{p}\in{\mbox{\bf R}}^m_{++}$ 
denote the corresponding quantities for the $m$ put options.

The amount we receive per dollar of investment in the underlying stocks
is $\omega$, the ratio of the current to final stock price.
For every dollar invested in the call options we receive
$(\omega-s_\text{c})_+ / p_\text{c}$, where the division is 
elementwise, and $(a)_+ = \max\{a,0\}$.
For the put options, the total we receive per dollar of investment is
$(s_\text{p}-\omega)_+ / p_\text{p}$.

We make investments in $n=3m$ different
assets, the $m$ underlying stocks and $m$ associated call and put options.
We let $x \in {\mbox{\bf R}}^n$ denote the fractions of our wealth that we invest in
the assets, so $\mathbf 1^T x =1$.  We consider long and short positions, with
$x_i <0$ denoting a short position.
We consider a simple set of portfolio constraints,
$x \geq x_\text{min}$ ({\it e.g.}, $x_\text{min}=-0.1$ 
limits the maximum short position for 
any asset to not exceed 10\% of the total portfolio value),
and $\|x\|_1 \leq L$, where $L$ is a leverage limit.
(Since $\mathbf 1^T x =1$, this means that the total long position cannot exceed
a multiple $(L+1)/2$ of total wealth, and the total short position
cannot exceed a multiple $(L-1)/2$ of the total wealth.)
We partition $x$ as $x=(x_\text{u},x_\text{c},x_\text{p})$, with
each subvector in ${\mbox{\bf R}}^m$.
Our portfolio has total return
\[
x_\text{u}^T \omega + x_\text{c}^T ((s_\text{c}-\omega)_+ / p_\text{c} ) 
+ x_\text{p}^T ((s_\text{p}-\omega)_+ / p_\text{p}) = r(\omega)^T x,
\]
where $r(\omega) \in {\mbox{\bf R}}^n$ is the total return of the $n=3m$ assets,
{\it i.e.}, stocks and options.

Our problem is to choose the portfolio so as to minimize the conditional
value at risk (described in \eqref{e-cvar}) of the negative total return, {\it i.e.},
\[
\begin{array}{ll}
\mbox{minimize} &  \mathop{\bf CVaR}(-r(\omega)^T x;\eta)\\
\mbox{subject to} & x \geq x_\text{min}, \; \mathbf 1^T x = 1, \; \|x\|_1 \leq L,
\end{array}
\]
where $\eta \in (0,1)$ sets the risk aversion.

We use a sample average approximation of the expectation in CVaR to 
obtain the problem
\[
\begin{array}{ll}
\mbox{minimize} &  
\frac{1}{N} \sum_{i=1}^N
\frac{(-r(\omega_i)^Tx - \alpha)_+}{1-\eta} + \alpha \\
\mbox{subject to} & x \geq x_\text{min}, \; \mathbf 1^T x = 1, \; \|x\|_1 \leq L
\end{array}
\]
with variables $x = (x_\text{u},x_\text{c},x_\text{p}) \in {\mbox{\bf R}}^n$
and $\alpha \in {\mbox{\bf R}}$.
The vectors $\omega_i \in {\mbox{\bf R}}_{++}^n$, $i=1,\ldots,N$, are independent 
samples of the price change $\omega$.

\paragraph{Problem instance.} 
We take the number of stocks as $m=100$, so $n=300$.
We take the minimum position limit $x_\text{min} = -0.1$
and leverage limit $L=1.6$.  We use risk aversion parameter $\eta = 0.8$,
so we are attempting to minimize the 20th percentile of the portfolio loss.
The price change covariance $\Sigma$ is generated according to 
$\Sigma = \sigma^2(I + 0.2 FF^T)$, where $\sigma = 1/\sqrt{2}$, 
and the entries of $F\in{\mbox{\bf R}}^{m \times 5}$ are 
independent draws from a standard normal distribution.
We set the mean price change according to $\mu_i = 0.03 \sqrt{\Sigma_{ii}} - 0.5\Sigma_{ii}$, $i=1,\ldots,m$. 
For each call option, the strike price is set as the 80th percentile
of $\omega$, and for each put option it is the 20th percentile.
The option prices are determined by the Black-Scholes formula
with zero discount.
(These data are approximately consistent with an investment period of 
one month for U.S.~equities.)

When we solve this problem instance,
the optimal portfolio return has mean 1.3\%,
standard deviation 5\%, and loss probability 39\%;
annualized, these correspond to a return of 16\%, 
standard deviation 18\%, and loss probability 20\%.
The optimal portfolio contains as assets the underlying stocks as well as call and put options.

\paragraph{Results.}
Table~\ref{tab:runtimes:cvar} shows the run-times for
$N$ ranging from $10{,}000$ to $1{,}000{,}000$.  
We see again that for small values of $N$, it is more efficient to
solve the problem directly using a structured solver, whereas for 
large values, OSMM is far more efficient.
When $N=1{,}000{,}000$, OSMM (using PyTorch) takes 0.0021 seconds to evaluate $f$ 
and 0.0053 seconds to evaluate $\nabla f$; 
OSMM takes 0.050 seconds to compute the tentative iterate $x^{k+1/2}$, 
and 0.021 seconds to evaluate the lower bound $L_k$ (using CVXPY).
Figure~\ref{fig:cvar_convergence} shows the convergence of OSMM 
with $N=1{,}000{,}000$.  
Practical accuracy are high accuracy are reached at the same time
after 100 iterations.

\begin{table}
  \centering
  \begin{tabular}[]{l llll}
    \toprule
	$N$ & OSMM & MOSEK & SCS & ECOS \\
    \midrule
    \midrule
	10{,}000 & 11 & 6.0 & 250 & 18 \\
	100{,}000 & 6.5 & 64 & 2{,}900 & 310\\
	1{,}000{,}000 & 6.3 & 1{,}900 & 30{,}000 & 5{,}200 \\ 
	\bottomrule
	\end{tabular}
	\caption{Solve times in seconds on the conditional value-at-risk problem.  A dash (``---'') means the solver failed, either for numerical reasons,
	or because it did not reach the required $10^{-6}$ suboptimality
	in 24 hours.}\label{tab:runtimes:cvar}
\end{table}

\begin{figure*}
\centering
\includegraphics{figures/cvar_cvg_top.pdf}
\includegraphics{figures/cvar_cvg_bottom.pdf}
\caption{Suboptimality, gap, sampling accuracy (top row), and RMS residual (bottom row) 
	on the conditional value-at-risk problem with $N=1{,}000{,}000$.
	The dotted green and black lines coincide, indicating that
	high accuracy and practical accuracy are reached at the same time.}
\label{fig:cvar_convergence}
\end{figure*}


\clearpage
\subsection{Exponential series density estimation}\label{subsec:experiment:ESE}
\paragraph{Problem formulation.} We consider an instance of the 
generic exponential family density fitting problem 
described in \S\ref{s-exp-fam-fit}.
We consider data $z_1, \ldots, z_m \in {\mbox{\bf R}}^2$, 
and wish to fit a density $p$ supported on $\mathcal S = [-1,1]^2$.
We let the sufficient statistics $\phi_i, \; i=1,\ldots,n$, be 
the Legendre polynomials up to degree four, so $n=14$.
(This is known as an exponential series 
density estimator~\cite{marsh2007goodness,wu2010exponential,gao2015penalized}.)
We solve the density estimation problem~\eqref{exp_density_prob_0}.

\paragraph{Problem instance.}   We take $m$ data points
sampled from a mixture of three Gaussian densities,
restricted to $\mathcal S$, with means $(1/3,1/3)$, $(1/3,-1/3)$, and
$(-1/3,-1/3)$, weights $0.4$, $0.3$, and $0.3$, 
and common covariance $(1/36) I$.
(So the data do not come from the family of density we use to fit.)
We form a Riemann sum using points in $\mathcal S$ lying on a grid with side lengths $\sqrt N$.
Recall that $N$, the number of samples, here refers to our approximate evaluation
of the integral that arises in the log-partition function, and not 
the number of data samples, which is fixed at $m=2{,}000$.


\paragraph{Results.}  
Table~\ref{tab:runtimes:density} shows the run-times for the various methods.
We see the usual pattern where directly solving the problem can be efficient for small $N$, but OSMM is much faster for large $N$.
When $N=1{,}000{,}000$,
it takes OSMM 0.001 and 0.0014 seconds to evaluate $f$ and $\nabla f$, 
respectively,
and 0.015 and 0.0097 seconds to compute $x^{k+1/2}$ and $L_k$, respectively.

Figure~\ref{fig:exp_density_Legendre_polynomials} shows the convergence 
of OSMM for $N=1{,}000{,}000$.
Practical accuracy is reached after 39 iterations with suboptimality $10^{-5}$, 
and after 42 iterations OSMM reaches high accuracy.
In this instance, our lower bound $L_k = -\infty$, so neither it nor the gap are plotted in the figure.

\begin{table}
  \centering
  \begin{tabular}[]{l llll}
    \toprule
	$N$ & OSMM & MOSEK & SCS & ECOS \\
    \midrule
    \midrule
    10{,}000 & 0.72 & 1.2 & 1{,}100 & 0.92 \\
	100{,}000 & 1.2 & 11 & --- & 14  \\ 
	1{,}000{,}000 &  0.84 &  120  & --- & 190 \\
	\bottomrule
	\end{tabular}
	\caption{Solve times in seconds on the density estimation problem.  A dash (``---'') means the solver failed, either for numerical reasons,
	or because it did not reach the required $10^{-6}$ suboptimality
	in 24 hours.}\label{tab:runtimes:density}
\end{table}

\begin{figure*}
\centering
\includegraphics{figures/exp_cvg_top.pdf}
\includegraphics{figures/exp_cvg_bottom.pdf}
\caption{Suboptimality (top row) and RMS residual (bottom row) 
	on the exponential family density fitting problem with $N=1{,}000{,}000$.
	The dotted green and black lines show when high accuracy and 
	practical accuracy are reached, respectively.}
\label{fig:exp_density_Legendre_polynomials}
\end{figure*}


\clearpage
\subsection{Vector news vendor with entropic value-at-risk}
We consider a variant of the classic news vendor problem that 
involves entropic value-at-risk programming~\cite{ahmed2007coherent}.

\paragraph{Problem formulation.} 
We choose quantities $q \in {\mbox{\bf R}}_+^n$ of $n$
products to produce, at total cost $\phi(q)$, where $\phi : {\mbox{\bf R}}^n_+ \to
{\mbox{\bf R}}_+$.  We have constraints on the quantities we can produce, 
$q \leq q_\mathrm{max}$, and on the total production cost, $\phi(q) \leq 
\phi_\mathrm{max}$.
We sell the amount $\min(q,d)$, where $d\in {\mbox{\bf R}}_+^n$ is the demand,
and the min is taken elementwise, at market prices $p \in {\mbox{\bf R}}_+^n$,
so the total revenue is $p^T \min(q,d)$, and the profit is
\[
P = p^T \min\{q,d\} - \phi(q).
\]
We assume that $\omega = (d,p) \in {\mbox{\bf R}}_+^{2n}$ is a random variable
with known distribution, so $P(q,\omega)$ 
is a random variable that depends on $q$.

We choose the quantities $q$ to minimize the EVaR of the negative
profit,
\[
\begin{array}{ll}
\mbox{minimize} & \mathop{\bf EVaR} (-P(q,\omega); \eta) \\ 
\mbox{subject to} &  \phi(q) \leq \phi_{\max}, \; q \geq 0, \; q \leq q_{\max},
\; \end{array}
\]
where $\eta$ is a specified quantile.

As described in \S\ref{s-EVaR}, 
we approximate this with Monte Carlo samples $(d_1,p_1), \ldots, 
(d_N,p_N)$ to obtain the problem
\[
\begin{array}{ll}
\mbox{minimize} & \alpha \log\left(\frac{1}{(1-\eta)N} \sum_{i=1}^N \exp
\left(\frac{-p_i^T\min(q, d_i) + \phi(q)}{\alpha} \right) \right)\\
\mbox{subject to} &  \phi(q) \leq \phi_{\max}, \; q \geq 0, \; q \leq q_{\max},
\; \alpha \geq 0, \end{array}
\]
with variables $q$ and $\alpha\in{\mbox{\bf R}}_+$.  We
denote the realizations of the demand $d$ and prices $p$ on the $i$th Monte
Carlo simulation by $d_i, p_i \in {\mbox{\bf R}}_+^n, \; i=1,\ldots,N$, respectively.

\paragraph{Problem instance.}  
We take $n=500$ and risk aversion parameter $\eta=0.9$.
We assume the demand and market prices follow 
a joint log-normal distribution, {\it i.e.}, 
$(d,p) = \exp z$, $z \sim \mathcal N(\mu, \Sigma)$, 
and the exponential is elementwise.  
We draw the entries of $\mu\in {\mbox{\bf R}}^{2n}$ independently from 
a uniform distribution on $[-0.2, 0]$,
and set $\Sigma = 0.1FF^T$,
where the entries of $F\in {\mbox{\bf R}}^{2n \times 5}$ are independently 
drawn from a standard normal distribution. 
We use a production cost which is separable and piecewise affine
with one kink point for each entry of $q$,
\[
\phi(q) = a^T q + 0.5 a^T (q - b)_+,
\]
where the elements in $a$ and $b$ are both drawn uniformly at random from $[0.2, 0.9]$ and $[0.01, 0.03]$, respectively.
The maximum production quantities $q_{\max}$ is set as 
$5b$, and the maximum cost is $\phi_\textrm{max} = 1$.
With these parameter values, the optimal
profit has mean 3.2 and standard deviation 0.98.

\paragraph{Results.}  
The run-times for the various methods are in 
table~\ref{tab:runtimes:vecnewsven}.
In this instance, OSMM is the fastest for all values of $N$ ranging from 
$1{,}000$ to $1{,}000{,}000$, and the other solvers fail for nearly all values of $N$.
When $N=1{,}000{,}000$, OSMM takes 0.16 and 0.27 seconds to evaluate 
$f$ and $\nabla f$, respectively; it takes 0.076 seconds to compute 
the tentative iterate, and 0.036 seconds to compute the lower bound. 

Figure~\ref{fig:vecnewsven_convergence} shows 
the convergence of OSMM with $N=1{,}000{,}000$.  
We can see that practical accuracy is reached after 50 
iterations with suboptimality on the order of $10^{-5}$, 
while it takes 58 iterations to reach high accuracy.

\begin{table}
	\centering
	\begin{tabular}[]{l llll}
		\toprule
		$N$ & OSMM & MOSEK & SCS & ECOS \\
		\midrule
		\midrule
		1{,}000 & 6.7 & 120 &  --- & ---\\
		10{,}000 & 11 & 7{,}400 & --- & --- \\
		100{,}000 & 10 & --- & --- & ---  \\ 
		1{,}000{,}000 & 53  &  ---  & --- & --- \\
		\bottomrule
	\end{tabular}
	\caption{Solve times in seconds on the vector news vendor problem.  
		A dash (``---'') means the solver failed, either for numerical reasons,
		or because it did not reach the required $10^{-6}$ suboptimality
		in 24 hours.}\label{tab:runtimes:vecnewsven}
\end{table}

\begin{figure*}
\centering
\includegraphics{figures/vecnewsven_cvg_top.pdf}
\includegraphics{figures/vecnewsven_cvg_bottom.pdf}
\caption{Suboptimality, gap, sampling accuracy (top row), and RMS residual (bottom row) 
	on the vector news vendor problem with $N=1{,}000{,}000$.
	The dotted green and black lines show when high accuracy and 
	practical accuracy are reached, respectively.}
\label{fig:vecnewsven_convergence}
\end{figure*}

\clearpage
\section*{Acknowledgment}
The work was supported in part by the National Key R\&D Program of China 
with grant No. 2018YFB1800800, 
by the Key Area R\&D Program of Guangdong Province 
with grant No. 2018B030338001, 
by Shenzhen Outstanding Talents Training Fund, 
and by Guangdong Research Project No. 2017ZT07X152.

