# Flexible Differentiable Optimization via Model Transformations

Akshay Sharma<EMAIL_ADDRESS>
Columbia University, NY, USA Mathieu Besançon<EMAIL_ADDRESS>
Zuse Institute Berlin, Germany Joaquim Dias Garcia<EMAIL_ADDRESS>
PSR & PUC-Rio, Rio de Janeiro, Brazil Benoît Legat<EMAIL_ADDRESS>
Massachussets Institute of Technology, MA, USA

###### Abstract

We introduce DiffOpt.jl, a Julia library to differentiate through the solution
of convex optimization problems with respect to arbitrary parameters present
in the objective and/or constraints. The library builds upon MathOptInterface,
thus leveraging the rich ecosystem of solvers and composing well with
modelling languages like JuMP. DiffOpt offers both forward and reverse
differentiation modes, enabling multiple use cases from hyperparameter
optimization to backpropagation and sensitivity analysis, bridging constrained
optimization with end-to-end differentiable programming.

##  1 Introduction

Differentiable Optimization ($\partial\mathcal{O}$) has become a center of
interest in the last years, both as a stand-alone methodology providing
additional information on the sensitivity of general convex optimization
problems and as a principled way to integrate complex convex optimization
components in Machine Learning, following a broader trend of differentiable
programming which extends the set of computer programs for which derivatives
can be computed [Innes et al., 2019]. We present DiffOpt, a library
implementing key methods for $\partial\mathcal{O}$ in Julia [Bezanson et al.,
2017] based on the MathOptInterface (MOI) abstract data structure [Legat et
al., 2021]. Consequently, the extensions to the usual mathematical
optimization machinery are also available in JuMP [Dunning et al., 2017] with
MOI as its backend.

The rest of this paper is structured as follows. Section 2 covers the
background and related research on $\partial\mathcal{O}$. Section 3 presents
the structure of the DiffOpt package and important features. Section 4
highlights some applications illustrating the use of the package.

##  2 Differentiating Convex Optimization Problems

$\partial\mathcal{O}$ regroups a set of methods to compute derivative
information of a function including an optimization problem. If we consider a
generic optimization problem parameterized by
$(\theta_{0},\theta_{1},\dots,\theta_{m})$:

$\displaystyle\mathcal{P}:\min_{x}\,\,\,$ $\displaystyle F(x,\theta_{0})$ s.t.
$\displaystyle G_{i}(x,\theta_{i})\in\mathcal{S}_{i}\,\,\,\forall
i\in\\{1,\dots,m\\}.$

Differentiable optimization in its generic form entails computing the
derivative of the output $x^{*}\in\operatorname*{argmin}_{x}\mathcal{P}$ with
respect to all parameters, i.e., computing the Jacobian matrix:

$\frac{\partial x^{*}}{\partial\theta}\equiv\bigg{(}\frac{\partial
x^{*}_{j}}{\partial\theta_{i}}\bigg{)}_{ij}.$

Instead of the Jacobian matrix of the _solution map_ , $\psi(\theta)=x^{*}$,
we often reason on the the derivative of the solution map as a linear map
denoted as $D\psi(\theta)$, with adjoint $D^{T}\psi(x)$. In general, the
solution map is not smooth, it is even a set-valued map and not a simple
function from parameters to solution. Some lines of work have investigated
non-smooth automatic (sub)differentiation and the behaviour of automatic
differentiation in the non-smooth case [Kakade and Lee, 2018, Bolte and
Pauwels, 2020]. However, in most cases in differentiable programming, a
heuristic quantity is computed when the limit of the directional derivative is
ill-defined111See the ChainRules.jl recommendations
https://juliadiff.org/ChainRulesCore.jl/v1.14/maths/nondiff_points.html,
accessed April 2022. The recent survey [Kotary et al., 2021] on learning with
constrained optimization offers a review of $\partial\mathcal{O}$ techniques.
We briefly go over the lines of work that are closest to the methods
implemented in DiffOpt.jl.

A part of the recent work on $\partial\mathcal{O}$ has focused on specific
optimization problems. In [Blondel et al., 2020], the authors consider a
differentiable optimization method for sorting and ranking problems. In
[Berthet et al., 2020], the differentiation of the solution to a generic
convex optimization problem with respect to a linear objective is considered.
When the feasible set is a polytope, the issue of a Jacobian matrix being zero
almost everywhere arises similarly to the ranking problem from [Blondel et
al., 2020]. The approach followed is that of perturbed optimization,
considering the input cost vector as a random variable centered around a
nominal value. This allows sampling the output solution and Jacobian matrix
from input cost vectors generated from the distribution, yielding unbiased
estimators for both while only requiring access to a linear minimization
oracle i.e. duality information is not required.

For differentiation with respect to arbitrary problem parameters, two main
methods have recently appeared in the literature. The first method, from [Amos
and Kolter, 2019], can be applied to convex quadratic optimization problems
while the second, from [Agrawal et al., 2020], is applied to conic
optimization problems. Both methods are based on rewriting necessary
optimality conditions for the optimization problem as a system of non-linear
equations and then applying the Implicit Function Theorem [Dontchev and
Rockafellar, 2009].

A major step towards integrating differentiable optimization into ML pipelines
was presented in [Amos and Kolter, 2019, Amos, 2019] for convex quadratic
problems (QP). The considered optimization models are of the form:

$\displaystyle\min_{x\in R^{n}}\,\,$ $\displaystyle x^{T}Qx+c^{T}x$
$\displaystyle Gx\leq h:(\lambda)$ $\displaystyle Ax=b:(\mu)$

where $\lambda\in R^{p}$, $\mu\in R^{m}$ denote the dual variables associated
with the inequality, equality constraints respectively. Unlike prior work, the
solution map is differentiated with respect to all problem data $(A\in
R^{m\times n},b\in R^{m},G\in R^{p\times n},h\in R^{p},Q\in R^{n\times n},c\in
R^{n})$. In particular, differentiating the solution with respect to
constraint coefficients opens new applications including learning the
constraints of the convex problem along with solutions as illustrated in
[Amos, 2019] on generic polytopes and Sudoku problems. The solution method
consists in representing the solution process as solving a system of
equations. The solution can then be differentiated with respect to its
parameters using implicit differentiation. In the case of QPs, the KKT
conditions of the system fully describe the optimality conditions for a
primal-dual solution:

$\displaystyle Qx+c+A^{T}\mu+G^{T}\lambda=0$ $\displaystyle(\nabla L)$
$\displaystyle Ax=b$ $\displaystyle(P_{eq})$ $\displaystyle 0\leq
h-Gx\,\,\bot\,\,\lambda\geq 0$ $\displaystyle(C),$

where $(\nabla L)$ represents the gradient of the Lagrangian, $(P_{eq})$
represents the primal feasibility of equality constraints and $(C)$ includes
primal feasibility of inequality constraints, complementarity and dual
feasibility for inequalities. A system of necessary equality constraints can
be derived:

$\displaystyle(\nabla L),\,\,(P_{eq}),$
$\displaystyle\lambda_{i}(h-Gx)_{i}=0\,\,\,\forall i.$

This representation can be implicitly differentiated:

$\displaystyle{{dQ}}x+Q{{dx}}+{{dc}}+{{dA}}^{T}\mu+A^{T}{{d\mu}}+{{dG}}^{T}\lambda+G^{T}{{d\lambda}}=0$
$\displaystyle{{dA}}x+A{{dx}}-{{db}}=0$
$\displaystyle{{d\lambda_{i}}}(h-Gx)_{i}+\lambda_{i}({{dh}}-{{dG}}\ x-G\
{{dx}})_{i}=0\,\,\,\forall i.$

Regrouping the differential forms of parameters and solution variables results
in the following system:

$\displaystyle\begin{bmatrix}Q&G^{T}&A^{T}\\\ D(\lambda)G&D(h-Gx)&0\\\
A&0&0\end{bmatrix}\begin{bmatrix}dx\\\ d\lambda\\\
d\mu\end{bmatrix}=-\begin{bmatrix}dQ\ x+dc+dG^{T}\lambda+dA^{T}\mu\\\
D(\lambda)dh-D(\lambda)dG\ x\\\ dAx-db\end{bmatrix}.$

where $D(\cdot)$ denotes the diagonal matrix formed from a given vector. The
system above can be used to compute several quantities of interest. At a given
point

$\centering(\theta=(c,Q,A,G,b,h),x,\lambda,\mu),\@add@centering$

solving the system can let us derive the full Jacobian matrix. As pointed out
in [Amos, 2019], the highly sparse structure of the matrix on the right-hand
side can be exploited to compute some partial derivatives first by selecting
only some rows, which are then back-substituted. However, in most applications
where derivative information is needed, we only need to perform Jacobian-
vector products (JVP) or vector-transpose-Jacobian products (VJP), thus
exposing the two principal modes of differentiation present in differentiable
libraries or Automatic Differentiation (AD) tools: forward- and reverse-mode
AD.

Although QPs capture numerous problems of interest, many problems require a
richer set of constraints which can be modelled as conic constraints. Convex
conic optimization problems have standard closed-form primal and dual
expressions,

Primal Problem Dual Problem $\displaystyle\min_{x\in R^{n}}\quad$
$\displaystyle c^{T}x\quad$ $\displaystyle\min_{y\in R^{m}}\quad$
$\displaystyle b^{T}y$ s.t. $\displaystyle Ax+s=b\quad$ s.t. $\displaystyle
A^{T}y+c=0$ $\displaystyle s\in\mathcal{K}$ $\displaystyle
y\in\mathcal{K}^{*}$

where $x$ is the primal variable, $y$ is the dual variable, $s\in R^{m}$ is
the primal slack variable, $\mathcal{K}\subseteq R^{m}$ is a closed convex
cone and $\mathcal{K}^{*}\subseteq R^{m}$ is the corresponding dual cone. Note
that in above form, $A\in R^{m\times n}$, $b\in R^{m}$, $c\in R^{n}$ are
problem data.

The methodology of $\partial\mathcal{O}$ has been developed for conic problems
in [Amos, 2019] and [Agrawal et al., 2020] in parallel and extended in
[Agrawal et al., 2019]. The conic case is more involved due to the conic
constraints, we follow closely the notation of [Agrawal et al., 2020] for a
complete and self-contained description of the method. Instead of deriving the
necessary system of equations from the KKT conditions, the proposed approach
will rely on the Homogeneous Self-Dual Embedding (HSDE) [O’Donoghue et al.,
2016]. In this case the solution map from the problem data to the primal dual
pair, $(x,y,s)=\psi(A,b,c)$, is represented as a composition of three
functions: $\psi=\phi\circ\mathcal{S}\circ\mathcal{Q}$. The definition of each
of these functions and their derivatives are given by:

  1. 1.

$\mathcal{Q}$ maps the problem data, $(A,b,c)$ to the a skew-symmetric matrix
$Q\in\mathbb{R}^{n+m+1}$:

$Q=\begin{bmatrix}0&A^{T}&c\\\ -A&0&b\\\
-c^{T}&-b^{T}&0\end{bmatrix},\text{and its differential is:
}d\mathcal{Q}=\begin{bmatrix}0&dA^{T}&dc\\\ -dA&0&db\\\
-dc^{T}&-db^{T}&0\end{bmatrix}$

  2. 2.

$\mathcal{S}$ maps the $Q$ matrix from the HSDE into its solution $z=(u,v,w)$.
In order to compute the derivative of this function, the author [Agrawal et
al., 2020] use the normalized residual map, defined in [Busseti et al., 2019]:

$\mathcal{N}(z,Q)=((Q-I)\Pi+I)\left(\frac{z}{|z_{n+m+1}|}\right),$

where $\Pi$ is the projection operator onto the Cartesian product
$\mathbb{R}^{n}\times\mathcal{K}^{*}\times\mathbb{R}_{+}$. They key point is
that, given a matrix $Q$, the solution $z$ of $\mathcal{N}(z,Q)=0$ is the
solution of the HSDE. Hence, it is possible to implicitly differentiate the
equation $\mathcal{N}(\mathcal{S}(Q),Q)=0$ and obtain:

$\displaystyle
D\mathcal{S}(Q)=-(D_{z}\mathcal{N}(\mathcal{S}(Q),Q))^{-1}D_{Q}\mathcal{N}(\mathcal{S}(Q),Q),$

where

$\displaystyle D_{Q}\mathcal{N}(z,Q)\left[U\right]=U\Pi(\frac{z}{z_{n+m+1}})$
$\displaystyle
D_{z}\mathcal{N}(z,Q)=\frac{((Q-I)D\Pi(z)+I)}{z_{n+m+1}}-\operatorname*{sign}(z_{n+m+1})((Q-I)\Pi+I)(\frac{z}{z_{n+m+1}^{2}})e_{n+m+1}^{T}$

where $D_{x}f(x,y)[v]$ denotes the directional derivative of $f$ at $x$ in
direction $v$ and $e_{n+m+1}$ the basis vector with 1 at index $n+m+1$. Note
that the second term of $D_{z}\mathcal{N}(z,Q)$ vanishes if $z$ is a solution
of the HSDE.

  3. 3.

$\phi$ maps the solution of the HSDE to the primal-dual pair:

$(x,y,s)=\phi(z)=\phi(u,v,w)=(u,\Pi_{\mathcal{K}^{*}}(v),\Pi_{\mathcal{K}^{*}}(v)-v)/z_{n+m+1}.$

Its derivative is given by:

$D\phi(z)=\begin{bmatrix}I&0&-x\\\ 0&D\Pi_{\mathcal{K}^{*}}(v)&-y\\\
0&D\Pi_{\mathcal{K}^{*}}(v)-I&s\end{bmatrix}$

The derivatives of projections onto classic cones like the positive orthant,
second-order , positive semidefinite, and exponential cones are given in
[Busseti et al., 2019] and implemented for DiffOpt in
MathOptSetDistances.jl[Besançon et al., 2022]. Exponential cone projections
are performed using the technique from [Friberg, 2021].

The methodology of differentiable conic problems has been extended to log-log
convex problems in [Agrawal and Boyd, 2020], using the grammar from
Disciplined Geometric Programming and the rules established for parameterized
disciplined convex problems in [Agrawal et al., 2019]. A differentiable method
has been developed for submodular functions in [Djolonga and Krause, 2018],
opening $\partial\mathcal{O}$ to optimization problems with discrete
structures. In [Gould et al., 2019], deep learning models are studied with
nonlinear optimization problems as nodes instead of closed-form functions,
defining a differentiable optimization method for nonlinear problems without
requiring convexity. The authors leverage implicit differentiation of the
Lagrangian reformulation at the optimal point to estimate derivative
information of the output solution with respect to input parameters of the
node.

The framework proposed in [Paulus et al., 2021] extends $\partial O$ to
optimization problems including integer constraints. Similar to previous work
in the convex setting and unlike prior models tackling combinatorial problems,
their method handles the differentiation of constraints using an estimation of
active constraints at the optimum in the backpropagation phase.

In [Blondel et al., 2021], differentiable optimization is viewed as a problem
of implicit differentiation. By expressing the solution map of the
optimization problem as the root of a system of equation (such as the KKT
conditions) or the solution to a fixed point equation, the system can leverage
automatic differentiation (AD) of the implicit equations to differentiate the
solution map. The approach generalizes the differentiation of quadratic and
conic optimization from [Amos and Kolter, 2019] and [Agrawal et al., 2019] but
requires a user-defined system of equations that are necessary for optimality.

##  3 Package structure

DiffOpt is a Julia package that offers differentiable optimization algorithms
to the JuMP ecosystem. In order to integrate seamlessly with other packages,
DiffOpt is built on top of MathOptInterface.jl (MOI), a foundational unifying
package for constrained optimization, designed to be a backend for modelling
interfaces such as JuMP.jl or Convex.jl. MOI allows the user to describe
structured optimization problems in a unified format based on the constraint
representation of functions of variables belonging to sets. This abstraction
covers a wide variety of problems including linear, quadratic, and conic
constraints as well as more specific sets, such as Special Ordered Sets or
complementarity constraints. Moreover, MOI includes a bridging mechanism
turning the user-provided problem into a problem structure that the chosen
solver accepts through successive transformations (or bridges) of the
function-set pairs.

The design of MOI enables extending its interface though various mechanisms.
First, constraints are defined as pairs func-in-set where func is a function
of the decision variables and set is a set in which the value of the function
should belong when evaluated at feasible points. Solvers that can exploit
special problem structure can therefore allow the user to communicate it by
defining new function or set types.

Second, most of the MOI interface is built on top of _attributes_. This
enables both the user to communicate custom information to the solver such as
starting values or callbacks but also the solver to communicate custom results
such as basis status for simplex solvers or Irreducible Inconsistent
Subsystems (IIS) for infeasible instances.

Third, MOI optimizers support a layered structure to combine several features.
These layers also supports custom constraints and attributes that are
typically defined by an inner layer thanks to the well-posed API of MOI. These
layers are commonly referred to as _meta-solvers_ as they form solvers
parameterized by other solvers. Bridges are defined as MOI layers that
transform constraints into constraints of different types for its inner layer.
When defining a new constraint type, defining bridges to transform it into
classical constraint types allows the user to encode this special structure in
models while still being able to use solvers not supporting this structure.
Below are a few examples of MOI extensions that illustrate this.

SumOfSquares.jl defines the Sum-of-Squares cone as a new set type [Weisser et
al., 2019] and _Gram matrices_ , _Moment matrices_ and _Sum-of-Squares
decompositions_ as new attributes. It then defines a bridge for transforming
Sum-of-Squares constraints into semidefinite constraints. The bridges is
applied for all semidefinite programming solvers except for the solver
Hypatia.jl [Coey et al., 2020] which supports this custom structure.

Dualization.jl [Bodin et al., 2021] offers a dualization meta-solver. The
optimization problem is automatically converted to its dual form and reaches
the internal solver only in dual form, of which the solution is then mapped
back to the user.

Other extensions have been developed, like QuadraticToBinary.jl [Garcia,
2021], another meta-solver that converts quadratically-constrained problems
into Mixed Integer Linear Programs. Also, constraint programming solvers have
defined new sets and functions [Kröger, 2020].

DiffOpt can be decomposed into 5 components and aspects that are covered in
the following five sections. First, it extends MOI by creating new attributes
allowing the user and solver to communicate forward or reverse differentiation
input and output. Second, it implements the quadratic and conic problem
differentiation rules described in Section 2 as MOI models implementing these
attributes. Third, it implements the communication of these attributes through
the MOI caching and bridging layers. Fourth, it implements a meta-solver that
implements the computation of the differentiation attributes and is
parametrized by a solver that should support solving the problem. Fifth, it
allows to integrate optimization layers into AD systems in Julia using the
ChainRulesCore.jl package.

###  3.1 Interface

After the problem is solved, the user can pass parameter perturbations in case
of forward differentiation or sensitivities with respect to solution variable
values in the case of reverse differentiation (or even both). These are passed
to the DiffOpt.Optimizer using the attributes detailed in Table 1. After
sensitivities are loaded, the user might call DiffOpt.forward_differentiate!
or DiffOpt.reverse_differentiate! to compute the derivatives with respect to
the input sensitivities. The resulting derivatives are queried again as
typical solver attributes, detailed in Table 2.

ForwardObjectiveFunction | Forward-mode tangent for the objective function
---|---
ForwardConstraintFunction | Forward-mode tangent for a constraint function
ReverseVariablePrimal | Reverse-mode tangent for a variable value
Table 1: Differential Optimization attributes (DiffOpt attributes for short) for passing perturbations and sensitivities. ForwardVariablePrimal | Forward-mode tangent for a variable value
---|---
ReverseObjectiveFunction | Reverse-mode tangent for the objective function
ReverseConstraintFunction | Reverse-mode tangent for a constraint function
Table 2: Differential Optimization attributes (DiffOpt attributes for short)
for querying resulting derivatives.

Perturbations are passed and queried in form of MOI functions (affine
functions and quadratic functions) associated with constraints or the
objective. The coefficients of variables (or quadratic terms) are the
perturbations related to those variables (or quadratic terms) and the
respective associated constraint or objective. This API allows the
perturbation to flow through model transformations defined by the above
mentioned bridges. To keep the API efficient the returned functions are lazily
computed so that the only the coefficients required by the final user are
actually evaluated.

DiffOpt integrates smoothly with JuMP just like any solver implementing the
MOI interface. Because derivative-related attributes are slightly more complex
than traditional attributes, JuMP was added as a dependency so that we could
overload a small subset of the JuMP API to achieve a better use experience.
Moreover, DiffOpt can communicate with any solver that has an MOI interface,
more than 30 of them are listed in the JuMP Manual [JuMP Developers, 2021].

###  3.2 Differentiation rules as MOI models

DiffOpt implements the differentiation rules for QPs (resp. CPs) described in
Section 2 as a _differentiation optimization MOI model_ (_DiffOpt model_ for
short) QuadraticDiffProblem (resp. ConicDiffProblem). These DiffOpt models
represent QPs (resp. CPs) in the matrix standard form described in Section 2.
They do not support solving the QP (resp. CP) but they support having the
primal and dual solution being set by the user. Then, they support the
differentiation API described in Section 3.1.

The QP and CP standard forms described in Section 2 seems quite restrictive.
For instance, the QP form supports inequality constraints $Gx\leq h$,
i.e.,MOI.ScalarAffineFunction-in-MOI.LessThan, but not $Gx\geq h$,i.e.,
MOI.ScalarAffineFunction-in-MOI.GreaterThan,. This requires the user to
transform their model in order to fit the solver-compatible representation. It
is then quite tedious and error-prone to map both the primal and dual
solutions through these transformations as well as the DiffOpt attributes.
Fortunately, as described in the next section, the transformation of these
attributes through bridges is implemented. Therefore, adding a bridging outer-
layer on top of DiffOpt.QuadraticDiffProblem or DiffOpt.ConicDiffProblem
allows the user to model the QP or CP in the most convenient form while all
these transformations are carried out transparently. Moreover, it also allows
extending DiffOpt to new problem classes through bridges. For instance, as
SumOfSquares.jl defines the transformation from a Sum-of-Squares constraint to
a semidefinite constraint using a bridge, defining how to transform the
DiffOpt attributes through these bridges automatically broadens the class of
programs supported by DiffOpt.ConicDiffProblem to Sum-of-Squares programs.

###  3.3 Model transformations

Most model transformations rely in an affine relation between two sets
$\mathcal{S}_{1}\subseteq\mathbb{R}^{n}$, $\mathcal{S}_{2}$ of the form

$\mathcal{S}_{1}=\\{\,x\in\mathbb{R}^{n}\mid\exists u\in\mathbb{R}^{m}\text{
s.t. }Ax+Bu+c\in\mathcal{S}_{2}\,\\}.$

In fact, at the time of writing, the only bridge in MOI not based on such
relation is the bridge from a convex quadratic constraint to a rotated second
order cone constraint as it relies on the Cholesky decomposition which is not
a linear map.

As detailed in [Legat, 2020, Section 2.1.2], there is an automated way to
implement the transformation for primal and dual results given the
transformation data $A$, $B$ and $c$. Similarly, we develop in this section
the transformation for the DiffOpt attributes.

Given a constraint $f_{1}(x)\in\mathcal{S}_{1}$, the bridges transforms it
into a constraint $f_{2}(x,u)\in\mathcal{S}_{2}$ by adding variables
$u\in\mathbb{R}^{m}$ where $f_{2}(x,u)=Af_{1}(x)+Bu+c$. We have
$\differential{f_{2}}=A\differential{f_{1}}$. Therefore, we see that, the
ForwardConstraintFunction forward-mode tangent $\Delta f_{1}$ should be mapped
to $A\Delta f_{1}$ and the ReverseConstraintFunction reverse-mode tangent
$\Delta f_{2}$ should be mapped to $A^{*}\Delta f_{2}$, with $A^{*}$ the
adjoint matrix of $A$ satisfying

$\langle Ax,y\rangle_{2}=\langle x,A^{*}y\rangle_{1}$

for all $x$ in the space of $\mathcal{S}_{1}$ and $y$ in the dual space of
$\mathcal{S}_{2}$ where $\langle\cdot,\cdot\rangle_{i}$ is the scalar product
between the space of $\mathcal{S}_{i}$ and its dual space for $i=1,2$.

###  3.4 Meta-solver

DiffOpt is designed as a meta-solver with a structure illustrated in Figure 1.
The main structure made available by the package is DiffOpt.Optimizer that is
parametrized by a mathematical programming inner solver. This inner solver may
be any object supporting storing and solving the problem provided by the user.
As shown in Figure 1, the entry point is a cache that is added to ensure
efficient storage and access of the user model. Also, bridge layers are added
as the user model may use constraint types that are not natively supported by
the inner solver (or the DiffOpt models) and need to be transformed. The
solution found by the inner solver is communicated to the DiffOpt model as
described in Section 3.2. This model is then used to compute forward and/or
reverse differentiation tangents.

Note that the bridges used for the inner solver and for the DiffOpt model are
completely independent. They can be entirely different as the primal and dual
results are automatically transformed through the inner solver bridge layer
from the solver solution into the solution corresponding to the user model and
then through the DiffOpt model bridges into the solution corresponding to the
DiffOpt model standard form.

This design of DiffOpt.Optimizer also enables adding new DiffOpt models in
addition to DiffOpt.QuadraticProgram.Model and DiffOpt.ConicProgram.Model in
order to further broaden the class of supported model.

[ [ [ Inner solver ].Bridges [ QuadraticProgram.Model ConicProgram.Model
].Bridges ].Cache ].DiffOpt.Optimizer

Figure 1: Design of the DiffOpt.Optimizer structure

Because MOI does not differentiate between the two considered classes QP and
CP, the user can pass any of the two problem classes to DiffOpt without ever
having to know which of the two methods will be used. DiffOpt automatically
selects the appropriate problem class based on the type of constraints and
objective function.

In the following code excerpt, we demonstrate the usage of DiffOpt, starting
from a simple JuMP model, then going through a reverse differentiation
procedure.

⬇

u

model = JuMP.Model(() -> diff_optimizer(Clp.Optimizer))

@variable(model, x)

@constraint(model, cons, x >= 3)

@objective(

model,

Min,

2x,

)

optimize!(model) # solve

MOI.set.( # set pertubations / gradient inputs

model,

DiffOpt.ReverseVariablePrimal(),

x,

1.0,

)

DiffOpt.reverse_differentiate!(model) # differentiate

# fetch expression of the gradient of constraint

grad_exp = MOI.get( # -3x+1

model,

DiffOpt.ReverseConstraintFunction(),

cons

)

JuMP.constant(grad_exp) # 1

JuMP.coefficient(grad_exp, x) # -3

###  3.5 ChainRules integration

Automatic Differentiation (AD) has become a cornerstone of machine learning,
evolving from the static transformation of program sources to a fully dynamic
process [Innes et al., 2019]. The Julia AD landscape has been evolving rapidly
[Schäfer et al., 2021] and is now converging towards a shared set of
derivative primitives implemented for elementary functions and exploited by AD
libraries to compute derivatives of full programs. The methods to implement to
provide derivative information are defined in ChainRulesCore.jl [JuliaDiff,
2021] and implemented in ChainRules.jl for Julia Base and standard library
functions.

One drawback is that the ChainRules.jl system reasons on the differentiation
of a function’s output with respect to its inputs. The DiffOpt interface
however is based on MOI and thus on the incremental construction of a model
represented as a single mutable object. One approach is to construct an
implementation of the solution map which takes as argument the model
parameters, builds and optimizes the model object, and returns the optimal
solution. The solution map is a pure function and its derivatives can be
expressed in terms of ChainRules.jl primitives and implemented using
derivative information from DiffOpt.

The implementation allows external users to effortlessly bring any MOI model
built directly or through a modelling interface like JuMP or Convex.jl to a
differentiable pipeline, regardless of the underlying solver used to solve the
main model. However, the user still needs to write down the solution map and
implement the ChainRules.jl interface functions frule and rrule for forward
and reverse differentiation respectively. We will demonstrate how DiffOpt and
ChainRules.jl can be combined in the following section.

##  4 Application examples

We discuss how differentiating an optimization program allows a variety of
applications for different computational tasks. These examples use various
convex solvers to compute the primal and dual solutions, in order to highlight
the ease to swap a solver for another in a single line.

All the following examples can be seen in detail in the Tutorials section of
the DiffOpt manual. The code and data use version 0.4 of the package.

###  4.1 Sensitivity Analysis

Sensitivity analysis [Saltelli et al., 2004] focuses on studying how the
changes in the inputs of a mathematical model affects its output.
Sensitivities of a JuMP model can be computed automatically using the
previously described methods. We illustrate sensitivity analysis for a
classification and a regression task.

####  4.1.1 Classification using SVM

Support vector machines (SVM) classify labelled data points with the
hyperplane minimizing the norm of classification errors on all points (or
achieve the largest margin if the two classes are separable). Assuming
$X\in\mathbb{R}^{n\times d}$ the feature matrix of $n$ data points with $d$
features and $y$ their labels, a soft-margin $\ell_{1}$-SVM with $\ell_{2}$
regularization can be modelled as:

$\displaystyle\min_{\xi,w,b}\,$
$\displaystyle\sum_{i=1}^{n}\xi_{i}+\lambda\|w\|^{2}$ s.t. $\displaystyle
y_{i}(w^{T}X_{i}+b)\geq 1-\xi_{i}\,\,\,\forall i\in 1..n$
$\displaystyle\xi_{i}\geq 0\,\,\,\forall i\in 1..n,$

where $e_{i}$ is the soft margin loss on the $i-th$ data point, $w^{T}x=b$ is
the SVM hyperplane, $\lambda$ is the regularization parameter.

The plots and code transcripts are from Sensitivity Analysis of SVM tutorial.
The model is implemented below and solved using Ipopt [Wächter and Biegler,
2006].

⬇

#

$\lambda$ = 0.05

model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))

# Add the variables

@variable(model, $\xi$[1:N] >= 0)

@variable(model, w[1:D])

@variable(model, b)

# Add the constraints.

@constraint(model,

con[i in 1:N],

y[i] * (dot(X[i,:], w) + b) >= 1 - $\xi$[i]

)

# Define the objective and solve

@objective(model, Min, $\lambda$ * dot(w, w) + sum($\xi$))

optimize!(model)

Using the forward differentiation mode of DiffOpt, we compute the partial
Jacobian with respect to each $i$-th individual data point and display its
norm:

$\left\|\frac{\partial w}{\partial X_{i}}\right\|+\left|\frac{\partial
b}{\partial X_{i}}\right|.$

The computation is performed as follows:

⬇

∇

for i in 1:N

for j in 1:N

if i == j

# identical perturbations on all x_i

MOI.set(

model,

DiffOpt.ForwardConstraintFunction(),

con[j],

y[j] * sum(w),

)

else

MOI.set(

model,

DiffOpt.ForwardConstraintFunction(),

con[j],

0.0,

)

end

end

DiffOpt.forward_differentiate!(model)

dw = MOI.get.(

model,

DiffOpt.ForwardVariablePrimal(),

w,

)

db = MOI.get(

model,

DiffOpt.ForwardVariablePrimal(),

b,

)

${\scriptstyle\nabla}$[i] = norm(dw) + norm(db)

end

A perturbation of the feature matrix $X$ affects the solution and can induce a
change in the separating hyperplane decisions $(w,b)$. Unlike other
classification models, not all points affect the hyperplane with small enough
perturbations, the optimal solution depends only on a few data points, the
support vectors that name the method. The impact of these perturbations are
displayed in Figure 2.

Figure 2: Learned SVM and sensitivities. Data points colors indicate the class
and marker size denoting the sensitivity of the hyperplane to $x_{i}$.

####  4.1.2 Ridge regression sensitivity

Ridge regression avoids overfitting with an $\ell_{2}$-norm penalty added to a
linear regression model and is particularly advantageous when the number of
features is as large as the number of observations. Assume
$X=\\{(x,y)\\}\subset\mathbb{R}^{d+1}$ to be the set of $n$ data points. Then
a ridge regression fitting problem can be modelled as:

$\displaystyle\min_{w,b}\quad$
$\displaystyle\sum_{i=1}^{N}(y_{i}-w^{T}x_{i}-b)^{2}+\alpha(\|w\|_{2}^{2}+b^{2})$

with $\alpha$ the regularization constant.

The plots and code transcripts are from Sensitivity Analysis of Ridge
Regression tutorial. We implement and solve below a univariate example with
DiffOpt using Ipopt as the underlying QP solver.

⬇

#

model = Model(() -> diff_optimizer(Ipopt.Optimizer))

@variable(model, w)

@variable(model, b)

alpha = 0.8 # regularization constant

@expression(model, e[i=1:N], Y[i] - w * X[i] - b)

@objective(

model,

Min,

1 / N * dot(e, e) + alpha * (w^2 + b^2),

)

optimize!(model)

Similar to the SVM example, a change in a single independent or dependent
variable value $x_{i}$ or $y_{i}$ for a given data point can affect the
learned model. We use DiffOpt in forward mode to quantify these sensitivities
to individual data points:

$\frac{\partial w}{\partial x_{i}}\text{ and }\frac{\partial w}{\partial
y_{i}}.$

A ForwardInObjective attribute can be set for the perturbation of the
objective function. It takes as input the expression proportionally dependent
on the perturbed parameter $\theta$. Given a generic expression
$f(x;\theta)=\theta g(x)$ with parameter $\theta$, the expected input is
$g(x)$. A particular aspect here is that the $x_{i}$ and $y_{i}$ values appear
as linear and quadratic term in the loss function. If the parameter $\theta$
appears both linearly and quadratically, the corresponding objective
perturbation can be derived with a first-order Taylor expansion:

$\displaystyle f(x;\theta)$ $\displaystyle=\theta g(x)+\theta^{2}h(x)$
$\displaystyle f(x;\theta+\delta)$
$\displaystyle=(\theta+\delta)g(x)+\theta^{2}h(x)+2\delta\theta
h(x)+\delta^{2}h(x)$ $\displaystyle f(x;\theta+\delta)$ $\displaystyle\approx
f(x;\theta)+\delta(g(x)+2\theta h(x)).$

When applied to the loss function, the perturbation $\delta^{(x)}_{i}$,
$\delta^{(y)}_{i}$ on $x_{i}$, $y_{i}$ respectively results in first-order
perturbation approximations:

$\displaystyle\delta^{(x)}_{i}(2w^{2}x_{i}+2bw-2wy_{i}),$
$\displaystyle\delta^{(y)}_{i}(2y_{i}-2b-2wx_{i})$

respectively. Using these input perturbations we can extract the output
sensitivity of the slope $w$ using the DiffOpt.ForwardOutVariablePrimal
variable attribute:

⬇

∇

${\scriptstyle\nabla}$x = zero(X)

for i in 1:N

MOI.set(

model,

DiffOpt.ForwardObjectiveFunction(),

2w^2 * X[i] + 2b * w - 2 * w * Y[i]

)

DiffOpt.forward_differentiate!(model)

${\scriptstyle\nabla}$x[i] = MOI.get(

model,

DiffOpt.ForwardVariablePrimal(),

w,

)

MOI.set(

model,

DiffOpt.ForwardObjectiveFunction(),

(2Y[i] - 2b - 2w * X[i]),

)

DiffOpt.forward_differentiate!(model)

${\scriptstyle\nabla}$y[i] = MOI.get(model,
DiffOpt.ForwardOutVariablePrimal(), w)

end

Figure 3 shows how sensitive in amplitude and direction is the slope $w$ to
perturbations of each data point.

(a) $x$ sensitivity of the regression slope $w$.

(b) $y$ sensitivity of the regression slope $w$.

Figure 3: Sensitivity analysis of the data points in ridge regression. Radius
of the markers is proportional to the sensitivity of the slope to x or y
perturbations. Blue markers indicate a negative sensitivity, red markers a
positive one.

###  4.2 Convex Optimization for Neural Network Layers

Most of the common neural network layers are either simple closed-form
operators or a composition of several operators. $\partial\mathcal{O}$ opens
new possibilities by providing derivatives for layers defined as solutions to
optimization problems. The derivatives computed by DiffOpt can be used to
backpropagate through the convex layer.

####  4.2.1 Custom ReLU layer

This example will follow the tutorial Custom ReLU layer of the documentation.

The Rectified Linear Unit or ReLU, a commonly used linear layer in machine
learning networks, is defined as $f(x)=\max\\{x,0\\}$. It can be interpreted
as projecting a point $x\in\mathbb{R}^{n}$ onto the non-negative orthant.

$\displaystyle\min_{y\geq 0}\quad$ $\displaystyle\|x-y\|_{2}^{2}$

where $y$ is the optimization variable and $x$ is the output.

Using the above definition, we model ReLU as a layer in a neural network
created in Flux.jl [Innes, 2018]. The model was trained on MNIST image dataset
[LeCun et al., 2010] with $60,000$ greyscale training bitmaps of size
$28\times 28$. We define the function matrix_relu with a matrix input because
it allows training in batches, the first dimension of the matrix is the layer
size, while the second dimension is the size of the batch.

⬇

f

y::Matrix;

model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))

)

layer_size, batch_size = size(y)

empty!(model)

set_silent(model)

@variable(model, x[1:layer_size, 1:batch_size] >= 0)

@objective(model, Min, x[:]'x[:] -2y[:]'x[:])

optimize!(model)

return value.(x)

end

Using our function as a neural network layer requires a corresponding
derivative implementation to differentiate the model and propagate the
gradients backward. This can be achieved using DiffOpt in conjunction with
ChainRules by defining a method for the reverse-mode primitive function rrule.

⬇

f

model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))

pv = matrix_relu(y, model = model)

function pullback_matrix_relu(dl_dx)

x = model[:x] # load decision variable `x` into scope

dl_dy = zeros(T, size(dl_dx))

dl_dq = zeros(T, size(dl_dx))

# set sensitivities

MOI.set.(model, DiffOpt.ReverseVariablePrimal(), x[:], dl_dx[:])

# compute grad

DiffOpt.reverse_differentiate!(model)

# return gradient wrt objective parameters

obj_exp = MOI.get(model, DiffOpt.ReverseObjectiveFunction())

# coeff of `x` in q'x = -2y'x

dl_dq[:] .= JuMP.coefficient.(obj_exp, x[:])

dq_dy = -2 # dq/dy = -2

dl_dy[:] .= dl_dq[:] * dq_dy

return (ChainRulesCore.NoTangent(), dl_dy)

end

return pv, pullback_matrix_relu

end

Note the ChainRulesCore.NoTangent term which corresponds to the derivative of
the output with respect to the function matrix_relu itself. The next example
will demonstrate how to handle stateful functions. Since this is a pure
function, no derivative is defined. We can now define the neural network
architecture including our custom layer and train it on the MNIST data.

⬇

u

using Flux

# neural network definition

layer_size = 10

m = Flux.Chain(

Flux.Dense(784, layer_size), # 784 being image linear dimension (28 x 28)

matrix_relu,

Flux.Dense(layer_size, 10), # 10 being the number of outcomes (0 to 9)

Flux.softmax,

)

# dataset preprocessing

N = 1000 # batch size

imgs = MLDatasets.MNIST.traintensor(1:N)

labels = MLDatasets.MNIST.trainlabels(1:N)

train_X = float.(reshape(imgs, size(imgs, 1) * size(imgs, 2), N))

train_Y = Flux.onehotbatch(labels, 0:9);

epochs = 50

dataset = repeated((train_X, train_Y), epochs)

# neural network optimization

custom_loss(x, y) = crossentropy(m(x), y)

opt = Flux.ADAM()

Flux.train!(custom_loss, params(m), dataset, opt)

Implementing a custom layer for a known closed-form function is not directly
useful, and solving a quadratic problem is costly in contrast with a simple
ReLU operation. However, it opens the door to more flexible variations of the
layer.

####  4.2.2 Polyhedral projection layer

We generalize the custom layer defined as an optimization problem from the
ReLU example. This use case is available in the tutorials as Polyhedral QP
layer.

Given $m$ vector-scalar pairs $(w_{i},b_{i})\forall i\in 1..m$, we define the
layer taking $y$ as input and projecting it on the polytope defined by the $m$
hyperplanes:

$\displaystyle\min_{x}$ $\displaystyle\|x-y\|^{2}_{2}$ s.t. $\displaystyle
w_{i}^{T}x\geq b_{i}\,\,\,\forall i\in 1..m.$

Instead of a function, we will represent the layer with a functor (or callable
object).

⬇

s

w::NTuple{N, Vector{Float64}}

b::Vector{Float64}

end

Polytope(w::NTuple{N}) where {N} = Polytope{N}(w, randn(N))

We define a "call" operation on the polytope, making it a so-called functor.
Calling the polytope with a matrix y operates an Euclidean projection of each
of the matrix columns onto the polytope.

⬇

f

y::AbstractMatrix;

model = direct_model(DiffOpt.diff_optimizer(Ipopt.Optimizer))

) where {N}

layer_size, batch_size = size(y)

empty!(model)

set_silent(model)

@variable(model, x[1:layer_size, 1:batch_size])

@constraint(model,

greater_than_cons[idx in 1:N, sample in 1:batch_size],

dot(polytope.w[idx], x[:, sample]) ${\scriptstyle\geq}$ polytope.b[idx]

)

@objective(model, Min, dot(x - y, x - y))

optimize!(model)

return JuMP.value.(x)

end

Flux.@functor Polytope

The @functor macro from Flux implements auxiliary functions for collecting the
parameters of our custom layer and operating backpropagation. Similarly to the
ReLU example, ChainRulesCore.rrule is used to implement the reverse-mode
differentiation of the layer.

⬇

f

polytope::Polytope{N},

y::AbstractMatrix) where {N}

model = direct_model(DiffOpt.diff_optimizer(Ipopt.Optimizer))

xv = polytope(y; model = model)

function pullback(dl_dx)

layer_size, batch_size = size(dl_dx)

dl_dx = ChainRulesCore.unthunk(dl_dx)

# `dl_dy` is the derivative of `l` wrt `y`

x = model[:x]

# grad wrt input parameters

dl_dy = zeros(size(dl_dx))

# grad wrt layer parameters

dl_dw = zero.(polytope.w)

dl_db = zero(polytope.b)

# set sensitivities

MOI.set.(model, DiffOpt.ReverseVariablePrimal(), x, dl_dx)

# compute grad

DiffOpt.reverse_differentiate!(model)

# compute gradient wrt objective function parameter y

obj_expr = MOI.get(model, DiffOpt.ReverseObjectiveFunction())

dl_dy .= -2 * JuMP.coefficient.(obj_expr, x)

greater_than_cons = model[:greater_than_cons]

for idx in 1:N, sample in 1:batch_size

cons_expr = MOI.get(model,

DiffOpt.ReverseConstraintFunction(),

greater_than_cons[idx, sample])

dl_db[idx] -= JuMP.constant(cons_expr)/batch_size

dl_dw[idx] .+= JuMP.coefficient.(cons_expr, x[:,sample])/batch_size

end

dself = ChainRulesCore.Tangent{Polytope{N}}(; w = dl_dw, b = dl_db)

return (dself, dl_dy)

end

return xv, pullback

end

Note that the inner pullback returns a ChainRulesCore.Tangent that represents
the tangent of a composite type. This will allow Flux to operate gradient
descent on the parameters of the Polytope struct directly. Similarly to the
previous example, we can now build and train the network (we omit other
details like dataset preprocessing):

⬇

l

m = Flux.Chain(

Flux.Dense(784, layer_size), # 784 being image linear dimension (28 x 28)

Polytope((randn(layer_size), randn(layer_size), randn(layer_size))),

Flux.Dense(layer_size, 10), # 10 being the number of outcomes (0 to 9)

Flux.softmax,

)

Flux.train!(custom_loss, Flux.params(m), dataset, opt)

The capacity to embed a convex problem as a neural network layer enables
formulations of SVMs or regression layers with custom constraints which have
direct applications in meta-learning [Lee et al., 2019].

###  4.3 Hyperparameter optimization

Most machine learning algorithms involve hyperparameters that require tuning
to accelerate the training process and reach good out-of-sample performance,
helping achieve balance between model variance and bias. In the past few
years, many developments were made in gradient-based methods [Maclaurin et
al., 2015] where researchers typically introduced extra hyperparameters for
tuning. Recent interest in automated machine learning (AutoML) has resulted in
a resurgence of research in this field [Feurer and Hutter, 2019].

When the learning process can be formulated as solving an optimization problem
of which the hyperparameters are parameters, the problem of out-of-sample
optimization can be formulated as a bilevel optimization problem [Guyon et
al., 2019]:

$\displaystyle\min_{\theta}\,\,\,$ $\displaystyle f(X_{test},\hat{w},\theta)$
s.t. $\displaystyle\hat{w}\in\operatorname*{argmin}_{w}f(X_{train},w,\theta),$

with $w$ the learned weights of the prediction model, $X_{train/test}$ the
training and testing data, $f$ the loss function and $\theta$ the
hyperparameter. Depending on the expression of
$\nabla_{\theta}f(\cdot,\cdot,\cdot)$, computing (even local) optima of the
bilevel optimization problem can be challenging. We showcase how DiffOpt can
be used to meta-optimize the weights and the hyperparameters following the
DiffOpt tutorial Auto-tuning Hyperparameters.

####  4.3.1 Optimization problem

Let $X=\\{(x,y)|x\in\mathbb{R}^{d},y\in\mathbb{R}\\}$ be the set of $N$ data
points of dimension $D$. The regularized linear model can be modelled as an
optimization problem of the form:

$\min_{w}\quad\frac{1}{2nD}\sum_{i=1}^{n}(y_{i}-w^{T}x_{i})^{2}+\frac{\alpha}{2D}\|w\|_{2}^{2}$
(1)

where $w\in\mathbb{R}^{D}$ are the learned weights and $\alpha$, the
regularization parameter, is the only hyperparameter. Since the problem is
strongly convex, it admits a unique minimum $w^{*}$. Its implementation in
JuMP is given below.

⬇

i

function fit_ridge(model, X, y, $\alpha$)

JuMP.empty!(model)

set_silent(model)

N, D = size(X)

@variable(model, w[1:D])

@expression(model, err_term, X * w - y)

@objective(

model,

Min,

dot(err_term, err_term) / (2 * N * D) + $\alpha$ * dot(w, w) / (2 * D),

)

optimize!(model)

return w

end

####  4.3.2 Model differentiation

We want to find the optimal regularization parameter for the loss on a test
set that was not used to find the optimal weights. We will apply gradient
descent of the unregularized test loss with respect to $\alpha$, for which we
need to compute its gradient $\dfrac{\partial l}{\partial\alpha}$. This can be
achieved using the chain rule:

$\frac{\partial
l}{\partial\alpha}(w,\alpha)=\nabla_{w}l(w,\alpha)\frac{\partial
w}{\partial\alpha}(w,\alpha)$

where $\dfrac{\partial w}{\partial\alpha}$ is the derivative of the optimal
solution of Problem (1) w.r.t. the parameter and can be found using DiffOpt as
follows:

⬇

f

D = length(w)

dw_d$\alpha$ = zeros(D)

MOI.set(

model,

DiffOpt.ForwardObjectiveFunction(),

dot(w, w) / (2 * D),

)

DiffOpt.forward_differentiate!(model)

for i in 1:D

dw_d$\alpha$[i] = MOI.get(

model,

DiffOpt.ForwardVariablePrimal(),

w[i],

)

end

return dw_d$\alpha$

end

function d_testloss_d$\alpha$(model, X_test, y_test, w, ŵ)

N, D = size(X_test)

dw_d$\alpha$ = compute_dw_d$\alpha$(model, w)

err_term = X_test * ŵ - y_test

return sum(eachindex(err_term)) do i

dot(X_test[i,:], dw_d$\alpha$) * err_term[i]

end / (N * D)

end

####  4.3.3 Hyperparameter gradient descent

The value of $\alpha$ is updated using a fixed-step gradient descent scheme
implemented below:

⬇

f

$\alpha$_s = Float64[]

${\scriptstyle\partial}$$\alpha$_s = Float64[]

test_loss = Float64[]

$\alpha$ = $\alpha$0

N, D = size(X_test)

model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))

for iter in 1:max_iters

w = fit_ridge(model, X_train, y_train, $\alpha$)

ŵ = value.(w)

err_term = X_test * ŵ - y_test

${\scriptstyle\partial}$$\alpha$ = d_testloss_d$\alpha$(model, X_test, y_test,
w, ŵ)

push!($\alpha$_s, $\alpha$)

push!(${\scriptstyle\partial}$$\alpha$_s, ${\scriptstyle\partial}$$\alpha$)

push!(test_loss, norm(err_term)^2 / (2 * N * D))

$\alpha$ -= fixed_step * ${\scriptstyle\partial}$$\alpha$

if abs(${\scriptstyle\partial}$$\alpha$) ${\scriptstyle\leq}$ grad_tol

break

end

end

return $\alpha$_s, ${\scriptstyle\partial}$$\alpha$_s, test_loss

end

####  4.3.4 Numerical results

The mean squared error of the regression on the training and test sets is
displayed in Figure 5, normalized for display. For values of $\alpha$ below
$\approx 0.23$, the model is under-regularized i.e. overfitted to the training
data, increasing $\alpha$ improves the test error. After that point,
increasing $\alpha$ excessively shrinks the regression coefficients,
increasing the error.

Figure 4: Mean squared error on the training and test set against $\alpha$

Figure 5: Gradient descent on the test set for $\alpha$ using DiffOpt

Figure 5 displays the trajectory of the gradient descent procedure described
in Section 4.3.3, starting from $\alpha_{0}=0.1$ and performing about 120
iterations.

Cross-validation is the typical procedure applied to tune the hyperparameters
of learning models. Leveraging $\partial\mathcal{O}$ lets us replace black-box
optimization procedures that are typical for hyperparameters with any first-
order method that leverages the gradient obtained from DiffOpt.

##  5 Conclusion

DiffOpt allows users to differentiate through optimization problems
implementing MathOptInterface, formulated with high-level modelling systems
like JuMP orConvex.jl. The package brings a flexible differentiable
optimization framework beyond disciplined convex optimization, separating the
underlying implementation from the MOI idiomatic interface. This separation of
concerns opens the possibility for the integration of novel differentiable
optimization techniques for classes of problems not already covered.

Future work will also consider performance improvements. In particular, using
DiffOpt to tune hyperparameters or construct a convex optimization layer
requires fast computation of derivatives and reoptimization of a perturbed
primal problem. The two steps could be integrated to cache more information
and reduce the computational burden per iteration.

## Acknowledgments

The work of A. Sharma on DiffOpt.jl was funded by the Google Summer of Code
program through NumFocus. M. Besançon was partially supported through the
Research Campus Modal funded by the German Federal Ministry of Education and
Research (fund numbers 05M14ZAM,05M20ZBM). J. Dias Garcia was supported in
part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -
Brasil (CAPES) - Finance Code 001. The authors would like to thank all
contributors and users of the package for their feedback and improvements, and
in particular Invenia Technical Computing for support, feedback on the API and
documentation, and the integration with ChainRules. We also thank Guillaume
Dalle for feedback and discussion on early versions of this manuscript.

## References

  * [Agrawal et al., 2019] Agrawal, A., Amos, B., Barratt, S., Boyd, S., Diamond, S., and Kolter, Z. (2019).  Differentiable Convex Optimization Layers.  arXiv:1910.12430 [cs, math, stat].  arXiv: 1910.12430.
  * [Agrawal et al., 2020] Agrawal, A., Barratt, S., Boyd, S., Busseti, E., and Moursi, W. M. (2020).  Differentiating Through a Cone Program.  arXiv:1904.09043 [math].  arXiv: 1904.09043.
  * [Agrawal and Boyd, 2020] Agrawal, A. and Boyd, S. (2020).  Differentiating through Log-Log Convex Programs.  arXiv:2004.12553 [math].  arXiv: 2004.12553.
  * [Amos, 2019] Amos, B. (2019).  Differentiable optimization-based modeling for machine learning.  PhD Thesis, PhD thesis. Carnegie Mellon University.
  * [Amos and Kolter, 2019] Amos, B. and Kolter, J. Z. (2019).  OptNet: Differentiable Optimization as a Layer in Neural Networks.  arXiv:1703.00443 [cs, math, stat].  arXiv: 1703.00443.
  * [Berthet et al., 2020] Berthet, Q., Blondel, M., Teboul, O., Cuturi, M., Vert, J.-P., and Bach, F. (2020).  Learning with Differentiable Perturbed Optimizers.  arXiv:2002.08676 [cs, math, stat].  arXiv: 2002.08676.
  * [Besançon et al., 2022] Besançon, M., Diamandis, T., Sharma, A., and Garcia, J. D. (2022).  MathOptSetDistances.jl: v0.2.2.  https://doi.org/10.5281/zenodo.6505452.
  * [Bezanson et al., 2017] Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B. (2017).  Julia: A fresh approach to numerical computing.  SIAM review, 59(1):65--98.
  * [Blondel et al., 2021] Blondel, M., Berthet, Q., Cuturi, M., Frostig, R., Hoyer, S., Llinares-López, F., Pedregosa, F., and Vert, J.-P. (2021).  Efficient and modular implicit differentiation.  arXiv preprint arXiv:2105.15183.
  * [Blondel et al., 2020] Blondel, M., Teboul, O., Berthet, Q., and Djolonga, J. (2020).  Fast differentiable sorting and ranking.  In International Conference on Machine Learning, pages 950--959. PMLR.
  * [Bodin et al., 2021] Bodin, G., Garcia, J. D., Legat, B., Besançon, M., Lubin, M., Dowson, O., and Cornejo, M. (2021).  Dualization.jl: v0.3.4.  https://doi.org/10.5281/zenodo.4718987.
  * [Bolte and Pauwels, 2020] Bolte, J. and Pauwels, E. (2020).  A mathematical model for automatic differentiation in machine learning.  Advances in Neural Information Processing Systems, 33:10809--10819.
  * [Busseti et al., 2019] Busseti, E., Moursi, W. M., and Boyd, S. (2019).  Solution refinement at regular points of conic problems.  Computational Optimization and Applications, 74(3):627--643.  Publisher: Springer.
  * [Coey et al., 2020] Coey, C., Kapelevich, L., and Vielma, J. P. (2020).  Solving natural conic formulations with hypatia. jl.  arXiv preprint arXiv:2005.01136.
  * [Djolonga and Krause, 2018] Djolonga, J. and Krause, A. (2018).  Differentiable learning of submodular models.  Advances in Neural Information Processing Systems 30, 2:1014--1024.
  * [Dontchev and Rockafellar, 2009] Dontchev, A. L. and Rockafellar, R. T. (2009).  Implicit functions and solution mappings, volume 543.  Springer.
  * [Dunning et al., 2017] Dunning, I., Huchette, J., and Lubin, M. (2017).  Jump: A modeling language for mathematical optimization.  SIAM Review, 59(2):295--320.
  * [Feurer and Hutter, 2019] Feurer, M. and Hutter, F. (2019).  Hyperparameter optimization.  Automated Machine Learning, pages 3--33.
  * [Friberg, 2021] Friberg, H. A. (2021).  Projection onto the exponential cone: a univariate root-finding problem.
  * [Garcia, 2021] Garcia, J. D. (2021).  QuadraticToBinary.jl: v0.2.4.  https://doi.org/10.5281/zenodo.4718981.
  * [Gould et al., 2019] Gould, S., Hartley, R., and Campbell, D. (2019).  Deep declarative networks: A new hope.  arXiv preprint arXiv:1909.04866.
  * [Guyon et al., 2019] Guyon, I., Sun-Hosoya, L., Boullé, M., Escalante, H. J., Escalera, S., Liu, Z., Jajetic, D., Ray, B., Saeed, M., Sebag, M., et al. (2019).  Analysis of the automl challenge series.  Automated Machine Learning, page 177.
  * [Innes, 2018] Innes, M. (2018).  Flux: Elegant machine learning with julia.  Journal of Open Source Software.
  * [Innes et al., 2019] Innes, M., Edelman, A., Fischer, K., Rackauckas, C., Saba, E., Shah, V. B., and Tebbutt, W. (2019).  A differentiable programming system to bridge machine learning and scientific computing.  arXiv preprint arXiv:1907.07587.
  * [JuliaDiff, 2021] JuliaDiff (2021).  Chainrules.jl, recipes for ad.  https://www.juliadiff.org/ChainRulesCore.jl/stable/.
  * [JuMP Developers, 2021] JuMP Developers (2021).  JuMP Manual.  https://jump.dev/JuMP.jl/v0.21.9/installation/#Supported-solvers.
  * [Kakade and Lee, 2018] Kakade, S. M. and Lee, J. D. (2018).  Provably correct automatic sub-differentiation for qualified programs.  Advances in neural information processing systems, 31.
  * [Kotary et al., 2021] Kotary, J., Fioretto, F., Van Hentenryck, P., and Wilder, B. (2021).  End-to-end constrained optimization learning: A survey.  arXiv preprint arXiv:2103.16378.
  * [Kröger, 2020] Kröger, O. (2020).  ConstraintSolver.jl.  https://github.com/Wikunia/ConstraintSolver.jl.
  * [LeCun et al., 2010] LeCun, Y., Cortes, C., and Burges, C. (2010).  Mnist handwritten digit database.  ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist, 2\.
  * [Lee et al., 2019] Lee, K., Maji, S., Ravichandran, A., and Soatto, S. (2019).  Meta-learning with differentiable convex optimization.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10657--10665.
  * [Legat, 2020] Legat, B. (2020).  Set programming : theory and computation.  PhD thesis, UCLouvain.
  * [Legat et al., 2021] Legat, B., Dowson, O., Garcia, J., and Lubin, M. (2021).  Mathoptinterface: a data structure for mathematical optimization problems.  INFORMS Journal on Computing (in press).
  * [Maclaurin et al., 2015] Maclaurin, D., Duvenaud, D., and Adams, R. (2015).  Gradient-based hyperparameter optimization through reversible learning.  In International conference on machine learning, pages 2113--2122. PMLR.
  * [O’Donoghue et al., 2016] O’Donoghue, B., Chu, E., Parikh, N., and Boyd, S. (2016).  Conic optimization via operator splitting and homogeneous self-dual embedding.  Journal of Optimization Theory and Applications, 169(3):1042--1068.
  * [Paulus et al., 2021] Paulus, A., Rolínek, M., Musil, V., Amos, B., and Martius, G. (2021).  Comboptnet: Fit the right np-hard problem by learning integer programming constraints.  arXiv preprint arXiv:2105.02343.
  * [Saltelli et al., 2004] Saltelli, A., Tarantola, S., Campolongo, F., and Ratto, M. (2004).  Sensitivity analysis in practice: a guide to assessing scientific models, volume 1.  Wiley Online Library.
  * [Schäfer et al., 2021] Schäfer, F., Tarek, M., White, L., and Rackauckas, C. (2021).  Abstractdifferentiation.jl: Backend-agnostic differentiable programming in julia.  arXiv preprint arXiv:2109.12449.
  * [Wächter and Biegler, 2006] Wächter, A. and Biegler, L. T. (2006).  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  Mathematical programming, 106(1):25--57.
  * [Weisser et al., 2019] Weisser, T., Legat, B., Coey, C., Kapelevich, L., and Vielma, J. P. (2019).  Polynomial and moment optimization in julia and jump.  In JuliaCon.