Year: 2021

Abstract:
Automatic differentiation (autodiff) has revolutionized machine learning. It allows expressing complex computations by composing elementary ones in creative ways and removes the burden of computing their derivatives by hand. More recently, differentiation of optimization problem solutions has attracted widespread attention with applications such as optimization layers, and in bi-level problems such as hyper-parameter optimization and meta-learning. However, so far, implicit differentiation remained difficult to use for practitioners, as it often required case-by-case tedious mathematical derivations and implementations. In this paper, we propose an efficient and modular approach for implicit differentiation of optimization problems. In our approach, the user defines directly in Python a function F capturing the optimality conditions of the problem to be differentiated. Once this is done, we leverage autodiff of F and implicit differentiation to automatically differentiate the optimization problem. Our approach thus combines the benefits of implicit differentiation and autodiff. It is efficient as it can be added on top of any state-of-the-art solver and modular as the optimality condition specification is decoupled from the implicit differentiation mechanism. We show that seemingly simple principles allow to recover many existing implicit differentiation methods and create new ones easily. We demonstrate the ease of formulating and solving bi-level optimization problems using our framework. We also showcase an application to the sensitivity analysis of molecular dynamics.

heading: 1 INTRODUCTION
text: Automatic differentiation (autodiff) is now an inherent part of machine learning software. It allows to express complex computations by composing elementary ones in creative ways and removes the tedious burden of computing their derivatives by hand. In parallel, the differentiation of optimization problem solutions has found many applications. A classical example is bi-level optimization, which typically involves computing the derivatives of a nested optimization problem in order to solve an outer one. Examples of applications in machine learning include hyper-parameter optimization (Chapelle et al., 2002; Seeger, 2008; Pedregosa, 2016; Franceschi et al., 2017; Bertrand et al., 2020; 2021), neural networks (Lorraine et al., 2020), and meta-learning (Franceschi et al., 2018; Rajeswaran et al., 2019). Another line of active research involving differentiation of optimization problem solutions are optimization layers (Kim et al., 2017; Amos & Kolter, 2017; Niculae & Blondel, 2017; Djolonga & Krause, 2017; Gould et al., 2019), which can be used to encourage structured outputs, and implicit deep networks (Bai et al., 2019; El Ghaoui et al., 2019), which have a smaller memory footprint than backprop-trained networks.
Since optimization problem solutions typically do not enjoy an explicit formula in terms of their inputs, autodiff cannot be used directly to differentiate these functions. In recent years, two main approaches have been developed to circumvent this problem. The first one consists of unrolling the iterations of an optimization algorithm and using the final iteration as a proxy for the optimization problem solution (Wengert, 1964; Domke, 2012; Deledalle et al., 2014; Franceschi et al., 2018; Ablin et al., 2020). This allows to explicitly construct a computational graph relating the algorithm output to the inputs, on which autodiff can then be used transparently. However, this requires a reimplementation of the algorithm using the autodiff system, and not all algorithms are necessarily autodiff friendly. Moreover, forward-mode autodiff has time complexity that scales linearly with the number of variables and reverse-mode autodiff has memory complexity that scales linearly with the number of algorithm iterations. In contrast, a second approach consists in implicitly relating an optimization problem solution to its inputs using optimality conditions. In a machine learning context, such implicit differentiation has been used for stationarity conditions (Bengio, 2000;
Lorraine et al., 2020), KKT conditions (Chapelle et al., 2002; Gould et al., 2016; Amos & Kolter, 2017; Niculae et al., 2018; Niculae & Martins, 2020) and the proximal gradient fixed point (Niculae & Blondel, 2017; Bertrand et al., 2020; 2021). An advantage of implicit differentiation is that a solver reimplementation is not needed, allowing to build upon decades of state-of-the-art software. Although implicit differentiation has a long history in numerical analysis (Griewank & Walther, 2008; Bell & Burke, 2008; Krantz & Parks, 2012; Bonnans & Shapiro, 2013), so far, it remained difficult to use for practitioners, as it required a case-by-case tedious mathematical derivation and implementation. CasADi (Andersson et al., 2019) allows to differentiate various optimization and root finding problem algorithms provided by the library. However, it does not allow to easily add implicit differentiation on top of existing solvers from optimality conditions expressed by the user, as we do. A recent tutorial explains how to implement implicit differentiation in JAX (Duvenaud et al., 2020). However, the tutorial requires the user to take care of low-level technical details and does not cover a large catalog of optimality condition mappings as we do. Other work (Agrawal et al., 2019a) attempts to address this issue by adding implicit differentiation on top of cvxpy (Diamond & Boyd, 2016). This works by reducing all convex optimization problems to a conic program and using conic programming’s optimality conditions to derive an implicit differentiation formula. While this approach is very generic, solving a convex optimization problem using a conic programming solver—an ADMM-based splitting conic solver (O’Donoghue et al., 2016) in the case of cvxpy—is rarely the state-of-the-art approach for each particular problem instance.
In this work, we adopt a different strategy that makes it easy to add implicit differentiation on top of any existing solver. In our approach, the user defines directly in Python a mapping function F capturing the optimality conditions of the problem solved by the algorithm. Once this is done, we leverage autodiff of F combined with implicit differentiation to automatically differentiate the optimization problem solution. In this way, our approach is generic, yet it can exploit the efficiency of state-of-the-art solvers. It therefore combines the benefits of implicit differentiation and autodiff. To summarize, we make the following contributions.
• We describe our framework and its JAX implementation (provided in the supplementary material). Our framework significantly lowers the barrier to use implicit differentiation thanks to the use of autodiff of the optimality conditions and the seamless integration in JAX. Our framework significantly extends JAX for numerical optimization, with low-level details all abstracted away.
• We instantiate our framework on a large catalog of optimality conditions (Table 1), recovering existing schemes and obtaining new ones, such as the mirror descent fixed point based one.
• On the theoretical side, we provide new bounds on the Jacobian error when the optimization problem is only solved approximately, and empirically validate them.
• We implement four illustrative applications, demonstrating our framework’s ease of use.
Beyond our software implementation, we hope this paper provides a self-contained blueprint for creating an efficient and modular implementation of implicit differentiation.
Notation. We denote the gradient and Hessian of f : Rd → R evaluated at x ∈ Rd by∇f(x) ∈ Rd and ∇2f(x) ∈ Rd×d. We denote the Jacobian of F : Rd → Rp evaluated at x ∈ Rd by ∂F (x) ∈ Rp×d. When f or F have several arguments, we denote the gradient, Hessian and Jacobian in the ith argument by ∇i, ∇2i and ∂i, respectively. The standard probability simplex is denoted by 4d := {x ∈ Rd : ‖x‖1 = 1, x ≥ 0}. For any set C ⊂ Rd, we denote the indicator function IC : Rd → R ∪ {+∞} where IC(x) = 0 if x ∈ C, IC(x) = +∞ otherwise. For a vector or matrix A, we note ‖A‖ the Frobenius (or Euclidean) norm, and ‖A‖op the operator norm.

heading: 2 COMBINING IMPLICIT DIFFERENTIATION AND AUTODIFF

heading: 2.1 GENERAL PRINCIPLES
text: Overview. Contrary to autodiff through unrolled algorithm iterations, implicit differentiation typically involves a manual, sometimes complicated, mathematical derivation. For instance, numerous works (Chapelle et al., 2002; Gould et al., 2016; Amos & Kolter, 2017; Niculae et al., 2018; Niculae & Martins, 2020) use Karush–Kuhn–Tucker (KKT) conditions in order to relate a constrained optimization problem’s solution to its inputs, and to manually derive a formula for its derivatives. The derivation and implementation in these works are always case-by-case.
In this work, we propose a generic way to easily add implicit differentiation on top of existing solvers. In our approach, the user defines directly in Python a mapping function F capturing the optimality conditions of the problem solved by the algorithm. We provide reusable building blocks to easily express such F . The provided F is then plugged into our Python decorator @custom_root, which we append on top of the solver declaration we wish to differentiate. Under the hood, we combine implicit differentiation and autodiff of F to automatically differentiate the optimization problem solution. A simple illustrative example is given in Figure 1.
Differentiating a root. Let F : Rd×Rn → Rd be a user-provided mapping, capturing the optimality conditions of a problem. An optimal solution, denoted x?(θ), should be a root of F :
F (x?(θ), θ) = 0 . (1) We can see x?(θ) as an implicitly defined function of θ ∈ Rn, i.e., x? : Rn → Rd. More precisely, from the implicit function theorem (Griewank & Walther, 2008; Krantz & Parks, 2012), we know that for (x0, θ0) satisfying F (x0, θ0) = 0 with a continuously differentiable F , if the Jacobian ∂1F evaluated at (x0, θ0) is a square invertible matrix, then there exists a function x?(·) defined on a neighborhood of θ0 such that x?(θ0) = x0. Furthermore, for all θ in this neighborhood, we have that F (x?(θ), θ) = 0 and ∂x?(θ) exists. Using the chain rule, the Jacobian ∂x?(θ) satisfies
∂1F (x ?(θ), θ)∂x?(θ) + ∂2F (x ?(θ), θ) = 0 .
Computing ∂x?(θ) therefore boils down to the resolution of the linear system of equations −∂1F (x?(θ), θ)︸ ︷︷ ︸
A∈Rd×d
∂x?(θ)︸ ︷︷ ︸ J∈Rd×n = ∂2F (x ?(θ), θ)︸ ︷︷ ︸ B∈Rd×n . (2)
When (1) is a one-dimensional root finding problem (d = 1), (2) becomes particularly simple since we then have ∇x?(θ) = B>/A, where A is a scalar value. We will show that existing and new implicit differentiation methods all reduce to this simple principle. We call our approach hybrid, since it combines implicit differentiation (and as such requires solving a linear system) with autodiff of the optimality conditions F . Our approach is efficient as it can be added on top of any state-of-the-art solver and modular as the optimality condition specification is decoupled from the implicit differentiation mechanism. This contrasts with existing works, where the mathematical derivation and implementation are specific to each optimality condition.
Differentiating a fixed point. We will encounter numerous applications where x?(θ) is implicitly defined through a fixed point:
x?(θ) = T (x?(θ), θ) ,
where T : Rd × Rn → Rd. This can be seen as a particular case of (1) by defining the residual F (x, θ) = T (x, θ)− x . (3)
In this case, using the chain rule, we have

heading: A = −∂1F (x?(θ), θ) = I − ∂1T (x?(θ), θ) and B = ∂2F (x?(θ), θ) = ∂2T (x?(θ), θ).
text: Computing JVPs and VJPs. In most practical scenarios, it is not necessary to explicitly form the Jacobian matrix, and instead it is sufficient to left-multiply or right-multiply by ∂1F and ∂2F . These are called vector-Jacobian product (VJP) and Jacobian-vector product (JVP), and are useful for integrating x?(θ) with reverse-mode and forward-mode autodiff, respectively. Oftentimes, F will be explicitly defined. In this case, computing the VJP or JVP can be done via autodiff. In some cases, F may itself be implicitly defined, for instance when F involves the solution of a variational problem. In this case, computing the VJP or JVP will itself involve implicit differentiation.
The right-multiplication (JVP) between J = ∂x?(θ) and a vector v, Jv, can be computed efficiently by solving A(Jv) = Bv. The left-multiplication (VJP) of v> with J , v>J , can be computed by first solving A>u = v. Then, we can obtain v>J by v>J = u>AJ = u>B. Note that when B changes but A and v remain the same, we do not need to solve A>u = v once again. This allows to compute the VJP w.r.t. different variables while solving only one linear system.
To solve these linear systems, we can use the conjugate gradient method (Hestenes et al., 1952) when A is symmetric positive semi-definite and GMRES (Saad & Schultz, 1986) or BiCGSTAB (Vorst & van der Vorst, 1992) otherwise. These algorithms are all matrix-free: they only require matrix-vector products. Thus, all we need from F is its JVPs or VJPs. An alternative to GMRES/BiCGSTAB is to solve the normal equation AA>u = Av using conjugate gradient. We implement this using JAX’s automatic transpose routine jax.linear_transpose (Frostig et al., 2021).
Pre-processing and post-processing mappings. Oftentimes, the goal is not to differentiate θ per se, but the parameters of a function producing θ. One example of such pre-processing is to convert the parameters to be differentiated from one form to another canonical form, such as a quadratic program (Amos & Kolter, 2017) or a conic program (Agrawal et al., 2019a). Another example is when x?(θ) is used as the output of a neural network layer, in which case θ is produced by the previous layer. Likewise, x?(θ) will often not be the final output we want to differentiate. One example of such post-processing is when x?(θ) is the solution of a dual program and we apply the dual-primal mapping to recover the solution of the primal program. Another example is the application of a loss function, in order to reduce x?(θ) to a scalar value. We leave the differentiation of such pre/post-processing mappings to the autodiff system, allowing to compose functions in complex ways.
Implementation details. When a solver function is decorated with @custom_root, we use jax.custom_jvp and jax.custom_vjp to automatically add custom JVP and VJP rules to the function, overriding JAX’s default behavior. As mentioned above, we use linear system solvers based on matrix-vector products and therefore we only need access to F through the JVP or VJP with ∂1F and ∂2F . This is done by using jax.jvp and jax.vjp, respectively. Note that, as in Figure 1, the definition of F will often include a gradient mapping ∇1f(x, θ). Thankfully, JAX supports second-order derivatives transparently. For convenience, our library also provides a @custom_fixed_point decorator, for adding implicit differentiation on top of a solver, given a fixed point iteration T ; see code examples in Appendix A.

heading: 2.2 EXAMPLES
text: We now give various examples of mapping F or fixed point iteration T , recovering existing implicit differentiation methods and creating new ones. Each choice of F or T implies different trade-offs in terms of computational oracles; see Table 1. Source code examples are given in Appendix A.
Stationary point condition. The simplest example is to differentiate through the implicit function x?(θ) = argmin
x∈Rd f(x, θ),
where f : Rd × Rn → R is twice differentiable. In this case, F is simply the gradient mapping F (x, θ) = ∇1f(x, θ). (4)
We then have ∂1F (x, θ) = ∇21f(x, θ) and ∂2F (x, θ) = ∂2∇1f(x, θ), the Hessian of f in its first argument and the Jacobian in the second argument of ∇1f(x, θ). In practice, we use autodiff to compute Jacobian products automatically. Equivalently, we can use the gradient descent fixed point
T (x, θ) = x− η∇1f(x, θ), (5) for all η > 0. Using (3), it is easy to check that we obtain the same linear system since η cancels out.
KKT conditions. We now show that the KKT conditions, manually differentiated in several works (Chapelle et al., 2002; Gould et al., 2016; Amos & Kolter, 2017; Niculae et al., 2018; Niculae & Martins, 2020), fit our framework. As we will see, the key will be to group the optimal primal and dual variables as our x?(θ). Let us consider the general problem
argmin z∈Rp
f(z, θ) subject to G(z, θ) ≤ 0, H(z, θ) = 0,
where z ∈ Rp is the primal variable, f : Rp×Rn → R, G : Rp×Rn → Rr and H : Rp×Rn → Rq . The stationarity, primal feasibility and complementary slackness conditions give
∇1f(z, θ) + [∂1G(z, θ)]>λ+ [∂1H(z, θ)]>ν = 0 H(z, θ) = 0
λ ◦G(z, θ) = 0, (6) where ν ∈ Rq and λ ∈ Rr+ are the dual variables, also known as KKT multipliers. The primal and dual feasibility conditions can be ignored almost everywhere. The system of (potentially nonlinear) equations (6) fits our framework, as we can group the primal and dual solutions as x?(θ) = (z?(θ), ν?(θ), λ?(θ)) to form the root of a function F (x?(θ), θ), where F : Rd×Rn → Rd and d = p + q + r. The primal and dual solutions can be obtained from a generic solver, such as an interior point method. In practice, the above mapping F will be defined directly in Python (see Figure 7 in Appendix A) and F will be differentiated automatically via autodiff.
Proximal gradient fixed point. Unfortunately, not all algorithms return both primal and dual solutions. Moreover, if the objective contains non-smooth terms, proximal gradient descent may be more efficient. We now discuss its fixed point (Niculae & Blondel, 2017; Bertrand et al., 2020; 2021). Let x?(θ) be implicitly defined as
x?(θ) := argmin x∈Rd f(x, θ) + g(x, θ),
where f : Rd × Rn → R is twice-differentiable convex and g : Rd × Rn → R is convex but possibly non-smooth. Let us define the proximity operator associated with g by
proxg(y, θ) := argmin x∈Rd
1 2 ‖x− y‖22 + g(x, θ).
To implicitly differentiate x?(θ), we use the fixed point mapping (Parikh & Boyd, 2014, p.150)
T (x, θ) = proxηg(x− η∇1f(x, θ), θ), (7) for any step size η > 0. The proximity operator is 1-Lipschitz continuous (Moreau, 1965). By Rademacher’s theorem, it is differentiable almost everywhere. Many proximity operators enjoy a closed form and can easily be differentiated, as discussed in Appendix B.
Projected gradient fixed point. As a special case, when g(x, θ) is the indicator function IC(θ)(x), where C(θ) is a convex set depending on θ, we obtain
x?(θ) = argmin x∈C(θ) f(x, θ). (8)
The proximity operator proxg becomes the Euclidean projection onto C(θ)
proxg(y, θ) = projC(y, θ) := argmin x∈C(θ) ‖x− y‖22
and (7) becomes the projected gradient fixed point
T (x, θ) = projC(x− η∇1f(x, θ), θ). (9)
Compared to the KKT conditions, this fixed point is particularly suitable when the projection enjoys a closed form. We discuss how to compute the JVP / VJP for a wealth of convex sets in Appendix B.
Mirror descent fixed point. We again consider the case when x?(θ) is implicitly defined as the solution of (8). We now generalize the projected gradient fixed point beyond Euclidean geometry. Let the Bregman divergence Dϕ : dom(ϕ)× relint(dom(ϕ))→ R+ generated by ϕ be defined by
Dϕ(x, y) := ϕ(x)− ϕ(y)− 〈∇ϕ(y), x− y〉.
We define the Bregman projection of y onto C(θ) ⊆ dom(ϕ) by
projϕC (y, θ) := argmin x∈C(θ) Dϕ(x,∇ϕ∗(y)). (10)
Definition (10) includes the mirror map∇ϕ∗(y) for convenience. It can be seen as a mapping from Rd to dom(ϕ), ensuring that (10) is well-defined. The mirror descent fixed point mapping is then
x̂ = ∇ϕ(x) y = x̂− η∇1f(x, θ)
T (x, θ) = projϕC (y, θ). (11)
Because T involves the composition of several functions, manually deriving its JVP/VJP is error prone. This shows that our approach leveraging autodiff allows to handle more advanced fixed point mappings. A common example of ϕ is ϕ(x) = 〈x, log x− 1〉, where dom(ϕ) = Rd+. In this case, Dϕ is the Kullback-Leibler divergence. An advantage of the Kullback-Leibler projection is that it sometimes easier to compute than the Euclidean projection, as we detail in Appendix B.
Other fixed points. More fixed points are described in Appendix C.

heading: 3 JACOBIAN PRECISION GUARANTEES
text: In practice, either by the limitations of finite precision arithmetic or because we perform a finite number of iterations, we rarely reach the exact solution x?(θ), but instead only reach an approximate solution x̂ and apply the implicit differentiation equation (2) at this approximate solution. This motivates the need for precision guarantees of this approach. We introduce the following formalism. Definition 1. Let F : Rd × Rn → Rd be an optimality criterion mapping. Let A := −∂1F and B := ∂2F . We define the Jacobian estimate at (x, θ) as the solution to the linear equation A(x, θ)J(x, θ) = B(x, θ). It is a function J : Rd × Rn → Rd×n.
It holds by construction that J(x?(θ), θ) = ∂x?(θ). Computing J(x̂, θ) for an approximate solution x̂ of x?(θ) therefore allows to approximate the true Jacobian ∂x?(θ). In practice, an algorithm used to solve (1) depends on θ. Note however that, what we compute is not the Jacobian of x̂(θ), unlike works differentiating through unrolled algorithm iterations, but an estimate of ∂x?(θ). We therefore use the notation x̂, leaving the dependence on θ implicit.
We develop bounds of the form ‖J(x̂, θ)− ∂x?(θ)‖ < C‖x̂− x?(θ)‖, hence showing that the error on the estimated Jacobian is at most of the same order as that of x̂ as an approximation of x?(θ). These bounds are based on the following main theorem, whose proof is included in Appendix D.
Theorem 1 (Jacobian estimate). Let F : Rd × Rn → Rd. Assume that there exist α, β, γ, ε, R > 0 such that A = −∂1F and B = ∂2F satisfy, for all v ∈ Rd, θ ∈ Rn and x such that ‖x−x?(θ)‖ ≤ ε:

heading: A is well-conditioned, Lipschitz: ‖A(x, θ)v‖ ≥ α‖v‖ , ‖A(x, θ)−A(x?(θ), θ)‖op ≤ γ‖x− x?(θ)‖.
text: B is bounded and Lipschitz: ‖B(x?(θ), θ)‖ ≤ R , ‖B(x, θ)−B(x?(θ), θ)‖ ≤ β‖x− x?(θ)‖. Under these conditions, when ‖x̂− x?(θ)‖ ≤ ε, we have
‖J(x̂, θ)− ∂x?(θ)‖ ≤ ( βα−1 + γRα−2 ) ‖x̂− x?(θ)‖ .
This result is inspired by (Higham, 2002, Theorem 7.2), that is concerned with the stability of solutions to inverse problems. Here we consider that A(·, θ) is uniformly well-conditioned, rather than only at x?(θ). This does not affect the first order in ε of this bound, and makes it valid for all x̂. It is also more tailored to applications to equation-specific cases.
Indeed, Theorem 1 can be applied to specific cases. In particular, for the gradient descent fixed point, where T (x, θ) = x− η∇1f(x, θ) and F (x, θ) = T (x, θ)− x, this yields
A(x, θ) = η∇21f(x, θ) and B(x, θ) = −η∂2∇1f(x, θ) . The guarantees on Jacobian precision under regularity conditions rely on f directly; see Corollary 1 in Appendix D. This reveals in particular that Jacobian estimation by implicit differentiation gains a factor of t compared to automatic differentiation, after t iterations of gradient descent in the strongly-convex setting (Ablin et al., 2020, Proposition 3.2). While our guarantees concern the Jacobian of x?(θ), we note that other studies (Grazzi et al., 2020; Ji et al., 2021; Bertrand et al., 2021) give guarantees on hypergradients (i.e., the gradient of an outer objective).
We illustrate these results on ridge regression, where x?(θ) = argminx ‖Φx− y‖22 + ∑ i θix 2 i . This problem has the merit that the solution x?(θ) and its Jacobian ∂x?(θ) are available in closed form. By running gradient descent for t iterations, we obtain an estimate x̂ of x?(θ) and an estimate J(x̂, θ) of ∂x?(θ); cf. Definition 1. By doing so for different numbers of iterations t, we can graph the relation between the error ‖x?(θ) − x̂‖2 and the error ‖∂x?(θ) − J(x̂, θ)‖2, as shown in Figure 2, empirically validating Theorem 1. The results in Figure 2 were obtained using the diabetes dataset from Efron et al. (2004), with other datasets yielding a qualitatively similar behavior. We derive similar guarantees in Corollary 2 in Appendix D for proximal gradient descent.

heading: 4 EXPERIMENTS
text: To conclude this work, we demonstrate the ease of solving bi-level optimization problems with our framework. We also present an application to the sensitivity analysis of molecular dynamics.

heading: 4.1 HYPERPARAMETER OPTIMIZATION OF MULTICLASS SVMS
text: In this example, we consider the hyperparameter optimization of multiclass SVMs (Crammer & Singer, 2001) trained in the dual. Here, x?(θ) is the optimal dual solution, a matrix of shape m × k, where m is the number of training examples and k is the number of classes, and θ ∈ R+ is the regularization parameter. The challenge in differentiating x?(θ) is that each row of x?(θ) is constrained to belong to the probability simplex4k. More formally, let Xtr ∈ Rm×p be the training feature matrix and Ytr ∈ {0, 1}m×k be the training labels (in row-wise one-hot encoding). Let W (x, θ) := X>tr(Ytr − x)/θ ∈ Rp×k be the dual-primal mapping. Then, we consider the following bi-level optimization problem
min θ=exp(λ)
1
2 ‖XvalW (x?(θ), θ)− Yval‖2F︸ ︷︷ ︸
outer problem
subject to x?(θ) = argmin x∈C
f(x, θ) := θ
2 ‖W (x, θ)‖2F + 〈x, Ytr〉︸ ︷︷ ︸
inner problem
,
where C = 4k × · · · ×4k is the Cartesian product of m probability simplices. We apply the change of variable θ = exp(λ) in order to guarantee that the hyper-parameter θ is positive. The matrix W (x?(θ), θ) ∈ Rp×k contains the optimal primal solution, the feature weights for each class. The outer loss is computed against validation data Xval and Yval.
While KKT conditions can be used to differentiate x?(θ), a more direct way is to use the projected gradient fixed point (9). The projection onto C can be easily computed by row-wise projections on the simplex. The projection’s Jacobian enjoys a closed form (Appendix B). Another way to differentiate x?(θ) is using the mirror descent fixed point (11). Under the KL geometry, projϕC (y, θ) corresponds to a row-wise softmax. It is therefore easy to compute and differentiate. Figure 3 compares the runtime performance of implicit differentiation vs. unrolling for the latter two fixed points.

heading: 4.2 DATASET DISTILLATION
text: Dataset distillation (Wang et al., 2018; Lorraine et al., 2020) aims to learn a small synthetic training dataset such that a model trained on this learned data set achieves a small loss on the original training set. Formally, let Xtr ∈ Rm×p and ytr ∈ [k]m denote the original training set. The distilled dataset will contain one prototype example for each class and therefore θ ∈ Rk×p. The dataset distillation problem can then naturally be cast as a bi-level problem, where in the inner problem we estimate a logistic regression model x?(θ) ∈ Rp×k trained on the distilled images θ ∈ Rk×p, while in the outer problem we want to minimize the loss achieved by x?(θ) over the training set:
min θ∈Rk×p f(x?(θ), Xtr; ytr)︸ ︷︷ ︸ outer problem subject to x?(θ) ∈ argmin x∈Rp×k f(x, θ; [k]) + ε‖x‖2︸ ︷︷ ︸ inner problem , (12)
where f(W,X; y) := `(y,XW ), ` denotes the multiclass logistic regression loss, and ε = 10−3 is a regularization parameter that we found had a very positive effect on convergence.
In this problem, and unlike in the general hyperparameter optimization setup, both the inner and outer problems are high-dimensional, making it an ideal test-bed for gradient-based bi-level optimization methods. For this experiment, we use the MNIST dataset. The number of parameters in the inner problem is p = 282 = 784. while the number of parameters of the outer loss is k × p = 7840. We solve this problem using gradient descent on both the inner and outer problem, with the gradient of the outer loss computed using implicit differentiation, as described in §2. This is fundamentally different from the approach used in the original paper, where they used differentiation of the unrolled iterates instead. For the same solver, we found that the implicit differentiation approach was 4
times faster than the original one. The obtained distilled images θ are visualized in Figure 4.

heading: 4.3 TASK-DRIVEN DICTIONARY LEARNING
text: Task-driven dictionary learning was proposed to learn sparse codes for input data in such a way that the codes solve an outer learning problem (Mairal et al., 2012; Sprechmann et al., 2014; Zarka et al., 2019). Formally, given a data matrix Xtr ∈ Rm×p and a dictionary of k atoms θ ∈ Rk×p, a sparse code is defined as a matrix x?(θ) ∈ Rm×k that minimizes in x a reconstruction loss f(x, θ) := `(Xtr, xθ) regularized by a sparsity-inducing penalty g(x). Instead of optimizing the dictionary θ to minimize the reconstruction loss, Mairal et al. (2012) proposed to optimize an outer problem that depends on the code. Given a set of labels Ytr ∈ {0, 1}m, we consider a logistic regression problem which results in the bilevel optimization problem:
min θ∈Rk×p,w∈Rk,b∈R σ(x?(θ)w + b; ytr)︸ ︷︷ ︸ outer problem subject to x?(θ) ∈ argmin x∈Rm×k f(x, θ) + g(x)︸ ︷︷ ︸ inner problem . (13)
When ` is the squared Frobenius distance between matrices, and g the elastic net penalty, Mairal et al. (2012, Eq. 21) derive manually, using optimality conditions (notably the support of the codes selected at the optimum), an explicit re-parameterization of x?(θ) as a linear system involving θ. This closed-form allows for a direct computation of the Jacobian of x? w.r.t. θ. Similarly, (Sprechmann et al., 2014) derive first order conditions in the case where ` is a β-divergence, while (Zarka et al., 2019) propose to use unrolling of ISTA iterations. Our approach bypasses all of these manual derivations, giving the user more leisure to focus directly on modeling (loss, regularizer) aspects.
We illustrate this on breast cancer survival prediction from gene expression data. We frame it as a binary classification problem to discriminate patients who survive longer than 5 years (m1 = 200) vs patients who die within 5 years of diagnosis (m0 = 99), from p = 1, 000 gene expression values. As shown in Table 2, solving (13) (Task-driven DictL) reaches a classification performance competitive with state-of-the-art L1 or L2 regularized logistic regression with 100 times fewer variables.

heading: 4.4 SENSITIVITY ANALYSIS OF MOLECULAR DYNAMICS
text: Many physical simulations require solving optimization problems, such as energy minimization in molecular (Schoenholz & Cubuk, 2020) and continuum (Beatson et al., 2020) mechanics, structural optimization (Hoyer et al., 2019) and data assimilation (Frerix et al., 2021). In this experiment, we revisit an example from JAXMD (Schoenholz & Cubuk, 2020), the problem of finding energy minimizing configurations to a system of k packed particles in a 2-dimensional box of size `
x?(θ) = argmin x∈Rk×2 f(x, θ) := ∑ i,j U(xi,j , θ),
where x?(θ) ∈ Rk×2 are the optimal coordinates of the k particles, U(xi,j , θ) is the pairwise potential energy function, with half the particles at diameter 1 and half at diameter θ = 0.6, which we optimize with a domain-specific optimizer (Bitzek et al., 2006). Here we consider sensitivity of particle position with respect to
diameter ∂x?(θ), rather than sensitivity of the total energy from the original experiment. Figure 5 shows results calculated via forward-mode implicit differentiation (JVP). Whereas differentiating the unrolled optimizer happens to work for total energy, here it typically does not even converge (see Appendix Figure 17), due the discontinuous optimization method.

heading: Appendix

heading: A CODE EXAMPLES

heading: A.1 CODE EXAMPLES FOR OPTIMALITY CONDITIONS
text: Our library provides several reusable optimality condition mappings F or fixed points T . We nevertheless demonstrate the ease of writing some of them from scratch.
Proximal gradient fixed point. The proximal gradient fixed point (7) with step size η = 1 is T (x, θ) = proxg(x−∇1f(x, θf ), θg). It can be implemented as follows.
We recall that when the proximity operator is a projection, we recover the projected gradient fixed point as a special case. Therefore, this fixed point can also be used for constrained optimization. We provide numerous proximal and projection operators in the library.
KKT conditions. As a more advanced example, we now describe how to implement the KKT conditions (6). The stationarity, primal feasibility and complementary slackness conditions read
∇1f(z, θf ) + [∂1G(z, θG)]>λ+ [∂1H(z, θH)]>ν = 0 H(z, θH) = 0
λ ◦G(z, θG) = 0. Using jax.vjp to compute vector-Jacobian products, this can be implemented as
Similar mappings F can be written if the optimization problem contains only equality constraints or only inequality constraints.
Mirror descent fixed point. Letting η = 1 and denoting θ = (θf , θproj), the fixed point (11) is
x̂ = ∇ϕ(x) y = x̂−∇1f(x, θf )
T (x, θ) = projϕC (y, θproj).
We can then implement it as follows.
Although not considered in this example, the mapping∇ϕ could also depend on θ if necessary.

heading: A.2 CODE EXAMPLES FOR EXPERIMENTS
text: We now sketch how to implement our experiments using our framework. In the following, jnp is short for jax.numpy. In all experiments, we only show how to compute gradients with the outer objective. We can then use these gradients with gradient-based solvers to solve the outer objective.

heading: Multiclass SVM experiment.
text: Task-driven dictionary learning experiment.

heading: Dataset distillation experiment.

heading: Molecular dynamics experiment.

heading: B JACOBIAN PRODUCTS
text: Our library provides numerous reusable building blocks. We describe in this section how to compute their Jacobian products. As a general guideline, whenever a projection enjoys a closed form, we leave the Jacobian product to the autodiff system.

heading: B.1 JACOBIAN PRODUCTS OF PROJECTIONS
text: We describe in this section how to compute the Jacobian products of the projections (in the Euclidean and KL senses) onto various convex sets. When the convex set does not depend on any variable, we simply denote it C instead of C(θ).
Non-negative orthant. When C is the non-negative orthant, C = Rd+, we obtain projC(y) = max(y, 0), where the maximum is evaluated element-wise. This is also known as the ReLu function. The projection in the KL sense reduces to the exponential function, projϕC (y) = exp(y).
Box constraints. When C(θ) is the box constraints C(θ) = [θ1, θ2]d with θ ∈ R2, we obtain
projC(y, θ) = clip(y, θ1, θ2) := max(min(y, θ2), θ1).
This is trivially extended to support different boxes for each coordinate, in which case θ ∈ Rd×2.
Probability simplex. When C is the standard probability simplex, C = 4d, there is no analytical solution for projC(y). Nevertheless, the projection can be computed exactly in O(d) expected time or O(d log d) worst-case time (Brucker, 1984; Michelot, 1986; Duchi et al., 2008; Condat, 2016). The Jacobian is given by diag(s)− ss>/‖s‖1, where s ∈ {0, 1}d is a vector indicating the support of projC(y) (Martins & Astudillo, 2016). The projection in the KL sense, on the other hand, enjoys a closed form: it reduces to the usual softmax projϕC (y) = exp(y)/ ∑d j=1 exp(yj).
Box sections. Consider now the Euclidean projection z?(θ) = projC(y, θ) onto the set C(θ) = {z ∈ Rd : αi ≤ zi ≤ βi, i ∈ [d];w>z = c}, where θ = (α, β, w, c). This projection is a singlyconstrained bounded quadratic program. It is easy to check (see, e.g., (Niculae & Martins, 2020)) that an optimal solution satisfies for all i ∈ [d]
z?i (θ) = [L(x ?(θ), θ)]i := clip(wix?(θ) + yi, αi, βi)
where L : R× Rn → Rd is the dual-primal mapping and x?(θ) ∈ R is the optimal dual variable of the linear constraint, which should be the root of
F (x?(θ), θ) = L(x?(θ), θ)>w − c. The root can be found, e.g., by bisection. The gradient∇x?(θ) is given by∇x?(θ) = B>/A and the Jacobian ∂z?(θ) is obtained by application of the chain rule on L.
Norm balls. When C(θ) = {x ∈ Rd : ‖x‖ ≤ θ}, where ‖ · ‖ is a norm and θ ∈ R+, projC(y, θ) becomes the projection onto a norm ball. The projection onto the `1-ball reduces to a projection onto the simplex, see, e.g., (Duchi et al., 2008). The projections onto the `2 and `∞ balls enjoy a closed-form, see, e.g., (Parikh & Boyd, 2014, §6.5). Since they rely on simple composition of functions, all three projections can therefore be automatically differentiated.
Affine sets. When C(θ) = {x ∈ Rd : Ax = b}, where A ∈ Rp×d, b ∈ Rp and θ = (A, b), we get
projC(y, θ) = y −A†(Ay − b) = y −A>(AA>)−1(Ay − b)
where A† is the Moore-Penrose pseudoinverse of A. The second equality holds if p < d and A is full rank. A practical implementation can pre-compute a factorization of the Gram matrix AA>. Alternatively, we can also use the KKT conditions.
Hyperplanes and half spaces. When C(θ) = {x ∈ Rd : a>x = b}, where a ∈ Rd and b ∈ R and θ = (a, b), we get
projC(y, θ) = y − a>y − b ‖a‖22 a.
When C(θ) = {x ∈ Rd : a>x ≤ b}, we simply replace a>y−b in the numerator by max(a>y−b, 0).
Transportation and Birkhoff polytopes. When C(θ) = {X ∈ Rp×d : X1d = θ1, X>1p = θ2, X ≥ 0}, the so-called transportation polytope, where θ1 ∈ 4p and θ2 ∈ 4d are marginals, we can compute approximately the projections, both in the Euclidean and KL senses, by switching to the dual or semi-dual (Blondel et al., 2018). Since both are unconstrained optimization problems, we can compute their Jacobian product by implicit differentiation using the gradient descent fixed point. An advantage of the KL geometry here is that we can use Sinkhorn (Cuturi, 2013), which is a GPU-friendly algorithm. The Birkhoff polytope, the set of doubly stochastic matrices, is obtained by fixing θ1 = θ2 = 1d/d.
Order simplex. When C(θ) = {x ∈ Rd : θ1 ≥ x1 ≥ x2 ≥ · · · ≥ xd ≥ θ2}, a so-called order simplex (Grotzinger & Witzgall, 1984; Blondel, 2019), the projection operations, both in the Euclidean and KL sense, reduce to isotonic optimization (Lim & Wright, 2016) and can be solved exactly in O(d log d) time using the Pool Adjacent Violators algorithm (Best et al., 2000). The Jacobian of the projections and efficient product with it are derived in (Djolonga & Krause, 2017; Blondel et al., 2020).
Polyhedra. More generally, we can consider polyhedra, i.e., sets of the form C(θ) = {x ∈ Rd : Ax = b, Cx ≤ d}, where A ∈ Rp×d, b ∈ Rp, C ∈ Rm×d, and d ∈ Rm. There are several ways to differentiate this projection. The first is to use the KKT conditions as detailed in §2.2. A second way is consider the dual of the projection instead, which is the maximization of a quadratic function subject to non-negative constraints (Parikh & Boyd, 2014, §6.2). That is, we can reduce the projection on a polyhedron to a problem of the form (8) with non-negative constraints, which we can in turn implicitly differentiate easily using the projected gradient fixed point, combined with the projection on the non-negative orthant. Finally, we apply the dual-primal mapping , which enjoys a closed form and is therefore amenable to autodiff, to obtain the primal projection.

heading: B.2 JACOBIAN PRODUCTS OF PROXIMITY OPERATORS
text: We provide several proximity operators, including for the lasso (soft thresholding), elastic net and group lasso (block soft thresholding). All satisfy closed form expressions and can be differentiated automatically via autodiff. For more advanced proximity operators, which do not enjoy a closed form, recent works have derived their Jacobians. The Jacobians of fused lasso and OSCAR were derived in (Niculae & Blondel, 2017). For general total variation, the Jacobians were derived in (Vaiter et al., 2013; Cherkaoui et al., 2020).

heading: C MORE EXAMPLES OF OPTIMALITY CRITERIA AND FIXED POINTS
text: To demonstrate the generality of our approach, we describe in this section more optimality mapping F or fixed point iteration T .
Newton fixed point. Let x be a root of G(·, θ), i.e., G(x, θ) = 0. The fixed point iteration of Newton’s method for root-finding is
T (x, θ) = x− η[∂1G(x, θ)]−1G(x, θ).
By the chain and product rules, we have
∂1T (x, θ) = I − η(...)G(x, θ)− η[∂1G(x, θ)]−1∂1G(x, θ) = (1− η)I.
Using (3), we get A = −∂1F (x, θ) = ηI . Similarly,
B = ∂2T (x, θ) = ∂2F (x, θ) = −η[∂1G(x, θ)]−1∂2G(x, θ).
Newton’s method for optimization is obtained by choosing G(x, θ) = ∇1f(x, θ), which gives
T (x, θ) = x− η[∇21f(x, θ)]−1∇1f(x, θ). (14)
It is easy to check that we recover the same linear system as for the gradient descent fixed point (5). A practical implementation can pre-compute an LU decomposition of ∂1G(x, θ), or a Cholesky decomposition if ∂1G(x, θ) is positive semi-definite.
Proximal block coordinate descent fixed point. We now consider the case when x?(θ) is implicitly defined as the solution
x?(θ) := argmin x∈Rd f(x, θ) + m∑ i=1 gi(xi, θ),
where g1, . . . , gm are possibly non-smooth functions operating on subvectors (blocks) x1, . . . , xm of x. In this case, we can use for i ∈ [m] the fixed point
xi = [T (x, θ)]i = proxηigi(xi − ηi[∇1f(x, θ)]i, θ), (15)
where η1, . . . , ηm are block-wise step sizes. Clearly, when the step sizes are shared, i.e., η1 = · · · = ηm = η, this fixed point is equivalent to the proximal gradient fixed point (7) with g(x, θ) =∑n i=1 gi(xi, θ).
Quadratic programming. We now show how to use the KKT conditions discussed in §2.2 to differentiate quadratic programs, recovering Optnet (Amos & Kolter, 2017) as a special case. To give some intuition, let us start with a simple equality-constrained quadratic program (QP)
argmin z∈Rp
f(z, θ) = 1
2 z>Qz + c>z subject to H(z, θ) = Ez − d = 0,
where Q ∈ Rp×p, E ∈ Rq×p, d ∈ Rq. We gather the differentiable parameters as θ = (Q,E, c, d). The stationarity and primal feasibility conditions give
∇1f(z, θ) + [∂1H(z, θ)]>ν = Qz + c+ E>ν = 0 H(z, θ) = Ez − d = 0.
In matrix notation, this can be rewritten as[ Q E>
E 0 ] [ z ν ] = [ −c d ] . (16)
We can write the solution of the linear system (16) as the root x = (z, ν) of a function F (x, θ). More generally, the QP can also include inequality constraints
argmin z∈Rp
f(z, θ) = 1
2 z>Qz+ c>z subject to H(z, θ) = Ez− d = 0, G(z, θ) = Mz− h ≤ 0.
where M ∈ Rr×p and h ∈ Rr. We gather the differentiable parameters as θ = (Q,E,M, c, d, h). The stationarity, primal feasibility and complementary slackness conditions give
∇1f(z, θ) + [∂1H(z, θ)]>ν + [∂1G(z, θ)]>λ = Qz + c+ E>ν +M>λ = 0 H(z, θ) = Ez − d = 0
λ ◦G(z, θ) = diag(λ)(Mz − h) = 0 In matrix notation, this can be written as Q E> M>E 0 0
diag(λ)M 0 0 [zν λ ] = [ −c d λ ◦ h ]
While x = (z, ν, λ) is no longer the solution of a linear system, it is the root of a function F (x, θ) and therefore fits our framework. With our framework, no derivation is needed. We simply define f , H and G directly in Python.
Conic programming. We now show that the differentiation of conic linear programs (Agrawal et al., 2019b; Amos, 2019), at the heart of differentiating through cvxpy layers (Agrawal et al., 2019a), easily fits our framework. Consider the problem
z?(λ), s?(λ) = argmin z∈Rp,s∈Rm c>z subject to Ez + s = d, s ∈ K, (17)
where λ = (c, E, d), E ∈ Rm×p, d ∈ Rm, c ∈ Rp and K ⊆ Rm is a cone; z and s are the primal and slack variables, respectively. Every convex optimization problem can be reduced to the form (17). Let us form the skew-symmetric matrix
θ(λ) =  0 E> c−E 0 d −c> −d> 0  ∈ RN×N , whereN = p+m+1. Following (Agrawal et al., 2019b;a; Amos, 2019), we can use the homogeneous self-dual embedding to reduce the process of solving (17) to finding a root of the residual map
F (x, θ) = θΠx+ Π∗x = ((θ − I)Π + I)x, (18) where Π = projRp×K∗×R+ and K
∗ ⊆ Rm is the dual cone. The splitting conic solver (O’Donoghue et al., 2016), which is based on ADMM, outputs a solution F (x?(θ), θ) = 0 which is decomposed as x?(θ) = (u?(θ), v?(θ), w?(θ)). We can then recover the optimal solution of (17) using
z?(λ) = u?(θ(λ)) and s?(λ) = projK∗(v ?(θ(λ)))− v?(θ(λ)).
The key oracle whose JVP/VJP we need is therefore Π, which is studied in (Ali et al., 2017). The projection onto a few cones is available in our library and can be used to express F .
Frank-Wolfe. We now consider
x?(θ) = argmin x∈C(θ)⊂Rd f(x, θ), (19)
where C(θ) is a convex polytope, i.e., it is the convex hull of vertices v1(θ), . . . , vm(θ). The FrankWolfe algorithm requires a linear minimization oracle (LMO)
s 7→ argmin x∈C(θ) 〈s, x〉
and is a popular algorithm when this LMO is easier to compute than the projection onto C(θ). However, since this LMO is piecewise constant, its Jacobian is null almost everywhere. Inspired by SparseMAP (Niculae et al., 2018), which corresponds to the case when f is a quadratic, we rewrite (19) as
p?(θ) = argmin p∈4m g(p, θ) := f(V (θ)p, θ),
where V (θ) is a d × m matrix gathering the vertices v1(θ), . . . , vm(θ). We then have x?(θ) = V (θ)p?(θ). Since we have reduced (19) to minimization over the simplex, we can use the projected gradient fixed point to obtain
T (p?(θ), θ) = proj4m(p ?(θ)−∇1g(p∗(θ), θ)).
We can therefore compute the derivatives of p?(θ) by implicit differentiation and the derivatives of x?(θ) by product rule. Frank-Wolfe implementations typically maintain the convex weights of the vertices, which we use to get an approximation of p?(θ). Moreover, it is well-known that after t iterations, at most t vertices are visited. We can leverage this sparsity to solve a smaller linear system. Moreover, in practice, we only need to compute VJPs of x?(θ).

heading: D JACOBIAN PRECISION PROOFS
text: Proof of Theorem 1. To simplify notations, we note A? := A(x?, θ) and Â := A(x̂, θ), and similarly for B and J . We have by definition of the Jacobian estimate function A?J? = B? and ÂĴ = B̂. Therefore we have
J(x̂, θ)− ∂x?(θ) = Â−1B̂ −A−1? B? = Â−1B̂ − Â−1B? + Â−1B? −A−1? B? = Â−1(B̂ −B?) + (Â−1 −A−1? )B? .
For any invertible matrices M1,M2, it holds that M−11 −M −1 2 = M −1 1 (M2 −M1)M −1 2 , so
‖M−12 −M −1 2 ‖op ≤ ‖M −1 1 ‖op‖M2 −M1‖op‖M −1 2 ‖op .
Therefore,
‖Â−1 −A−1? ‖op ≤ 1
α2 ‖Â−A?‖op ≤
γ
α2 ‖x̂− x?(θ)‖ .
As a consequence, the second term in J(x̂, θ)− ∂x?(θ) can be upper bounded and we obtain
‖J(x̂, θ)− ∂x?(θ)‖ ≤ ‖Â−1(B̂ −B?)‖+ ‖(Â−1 −A−1? )B?‖
≤ ‖Â−1‖op‖B̂ −B?‖+ γ
α2 ‖x̂− x?(θ)‖ ‖B?‖ ,
which yields the desired result.
Corollary 1 (Jacobian precision for gradient descent fixed point). Let f be such that f(·, θ) is twice differentiable and α-strongly convex and ∇21f(·, θ) is γ-Lipschitz (in the operator norm) and ∂2∇1f(x, θ) is β-Lipschitz and bounded in norm by R. The estimated Jacobian evaluated at x̂ is then given by J(x̂, θ) = −(∇21f(x̂, θ))−1∂2∇1f(x̂, θ) . For all θ ∈ Rn, and any x̂ estimating x?(θ), we have the following bound for the approximation error of the estimated Jacobian
‖J(x̂, θ)− ∂x?(θ)‖ ≤ ( β
α + γR α2
) ‖x̂− x?(θ)‖ .
Proof of Corollary 1. This follows from Theorem 1, applied to this specificA(x, θ) andB(x, θ).
For proximal gradient descent, where T (x, θ) = proxηg(x− η∇1f(x, θ), θ), this yields
A(x, θ) = I − ∂1T (x, θ) = I − ∂1proxηg(x− η∇1f(x, θ), θ)(I − η∇21f(x, θ)) B(x, θ) = ∂2proxηg(x− η∇1f(x, θ), θ)− η∂1proxηg(x− η∇1f(x, θ), θ)∂2∇1f(x, θ) .
We now focus in the case of proximal gradient descent on an objective f(x, θ) + g(x), where g is smooth and does not depend on θ. This is the case in our experiments in §4.3. Recent work also exploits local smoothness of solutions to derive similar bounds (Bertrand et al., 2021, Theorem 13) Corollary 2 (Jacobian precision for proximal gradient descent fixed point). Let f be such that f(·, θ) is twice differentiable and α-strongly convex and ∇21f(·, θ) is γ-Lipschitz (in the operator norm) and ∂2∇1f(x, θ) is β-Lipschitz and bounded in norm by R. Let g : Rd → R be a twicedifferentiable µ-strongly convex (with special case µ = 0 being only convex), for which the function Γη(x, θ) = ∇2g(proxηg(x−η∇1f(x, θ)) is κη-Lipschitz in it first argument. The estimated Jacobian evaluated at x̂ is then given by
J(x̂, θ) = −(∇21f(x̂, θ) + Γη(x̂, θ))−1∂2∇1f(x̂, θ) .
For all θ ∈ Rn, and any x̂ estimating x?(θ), we have the following bound for the approximation error of the estimated Jacobian
‖J(x̂, θ)− ∂x?(θ)‖ ≤ ( β + κη α+ µ + γR (α+ µ)2 ) ‖x̂− x?(θ)‖ .
Proof of Corollary 2. First, let us note that proxηg(y, θ) does not depend on θ, since g itself does not depend on θ, and is therefore equal to classical proximity operator of ηg which, with a slight overload of notations, we denote as proxηg(y) (with a single argument). In other words,
proxηg(y, θ) = proxηg(y) , ∂1proxηg(y, θ) = ∂proxηg(y) , ∂2proxηg(y, θ) = 0 .
Regarding the first claim (expression of the estimated Jacobian evaluated at x̂), we first have that proxηg(y) is the solution to (x
′ − y) + η∇g(x′) = 0 in x′ - by first-order condition for a smooth convex function. We therefore have that
proxηg(y) = (I + η∇g)−1(y) ∂proxηg(y) = (Id + η∇2g(proxηg(y)))−1 ,
the first I and inverse being functional identity and inverse, and the second Id and inverse being in the matrix sense, by inverse rule for Jacobians ∂h(z) = [∂h−1(h(z))]−1 (applied to the prox).
As a consequence, we have, for Γη(x, θ) = ∇2g(proxηg(x− η∇1f(x, θ)) that
A(x, θ) = Id − (Id + ηΓη(x, θ))−1(Id − η∇21f(x, θ)) = (Id + ηΓη(x, θ))
−1[Id + ηΓη(x, θ)− (Id − η∇21f(x, θ))] = η(Id + ηΓη(x, θ))
−1(∇21f(x, θ) + Γη(x, θ)) B(x, θ) = −η(Id + ηΓη(x, θ))−1∂2∇1f(x, θ) .
As a consequence, for all x ∈ Rd, we have that
J(x, θ) = −(∇21f(x, θ) + Γη(x, θ))−1∂2∇1f(x, θ) .
In the following, we modify slightly the notation of both A and B, writing
Ã(x, θ) = ∇21f(x, θ) + Γη(x, θ) B̃(x, θ) = −∂2∇1f(x, θ) .
With the current hypotheses, following along the proof of Theorem 1, we have that Ã is (α + µ) well-conditioned, and (γ + κη)-Lipschitz in its first argument, and B̃ is β-Lipschitz in its first argument and bounded in norm by R. The same reasoning yields
‖J(x̂, θ)− ∂x?(θ)‖ ≤ ( β + κη α+ µ + γR (α+ µ)2 ) ‖x̂− x?(θ)‖ .

heading: E EXPERIMENTAL SETUP AND ADDITIONAL RESULTS
text: Our experiments use JAX (Bradbury et al., 2018), which is Apache2-licensed and scikit-learn (Pedregosa et al., 2011), which is BSD-licensed.

heading: E.1 HYPERPARAMETER OPTIMIZATION OF MULTICLASS SVMS
text: Experimental setup. Synthetic datasets were generated using scikit-learn’s sklearn.datasets.make_classification (Pedregosa et al., 2011), following a model adapted from (Guyon, 2003). All datasets consist of m = 700 training samples belonging to k = 5 distinct classes. To simulate problems of different sizes, the number of features is varied as p ∈ {100, 250, 500, 750, 1000, 2000, 3000, 4000, 5000, 7500, 10000}, with 10% of features being informative and the rest random noise. In all cases, an additional mval = 200 validation samples were generated from the same model to define the outer problem.
For the inner problem, we employed three different solvers: (i) mirror descent, (ii) (accelerated) proximal gradient descent and (iii) block coordinate descent. Hyperparameters for all solvers were individually tuned manually to ensure convergence across the range of problem sizes. For mirror descent, a stepsize of 1.0 was used for the first 100 steps, following a inverse square root decay afterwards up to a total of 2500 steps. For proximal gradient descent, a stepsize of 5 · 10−4 was used for 2500 steps. The block coordinate descent solver was run for 500 iterations. All solvers used the same initialization, namely, xinit = 1k1m×k, which satisfies the dual constraints.
For the outer problem, gradient descent was used with a stepsize of 5 · 10−3 for the first 100 steps, following a inverse square root decay afterwards up to a total of 150 steps.
Conjugate gradient was used to solve the linear systems in implicit differentiation for at most 2500 iterations.
All results reported pertaining CPU runtimes were obtained using an internal compute cluster. GPU results were obtained using a single NVIDIA P100 GPU with 16GB of memory per dataset. For each dataset size, we report the average runtime of an individual iteration in the outer problem, alongside a 90% confidence interval estimated from the corresponding 150 runtime values.
Additional results Figure 13 compares the runtime of implicit differentiation and unrolling on GPU. These results highlight a fundamental limitation of the unrolling approach in memory-limited systems such as accelerators, as the inner solver suffered from out-of-memory errors for most problem sizes (p ≥ 2000 for mirror descent, p ≥ 750 for proximal gradient and block coordinate descent). While it might be possible to ameliorate this limitation by reducing the maximum number of iterations in the inner solver, doing so might lead to additional challenges (Wu et al., 2018) and require careful tuning.
Figure 14 depicts the validation loss (value of the outer problem objective function) at convergence. It shows that all approaches were able to solve the outer problem, with solutions produced by different approaches being qualitatively indistinguishable from each other across the range of problem sizes considered.
Figure 15 shows the Jacobian error achieved as a function of the solution error, when varying the number of features.

heading: E.2 TASK-DRIVEN DICTIONARY LEARNING
text: We downloaded from http://acgt.cs.tau.ac.il/multi_omic_benchmark/ download.html a set of breast cancer gene expression data together with survival information generated by the TCGA Research Network (https://www.cancer.gov/tcga) and processed as explained by (Rappoport & Shamir, 2018). The gene expression matrix contains the expression value for p=20,531 genes in m=1,212 samples, from which we keep only the primary tumors (m=1,093). From the survival information, we select the patients who survived at least five years after diagnosis (m1 = 200), and the patients who died before five years (m0 = 99), resulting in a cohort of m = 299 patients with gene expression and binary label. Note that non-selected patients are those who are marked as alive but were not followed for 5 years.
To evaluate different binary classification methods on this cohort, we repeated 10 times a random split of the full cohort into a training (60%), validation (20%) and test (20%) sets. For each split and each method, 1) the method is trained with different parameters on the training set, 2) the parameter that maximizes the classification AUC on the validation set is selected, 3) the method is then re-trained on the union of the training and validation sets with the selected parameter, and 4) we measure the AUC of that model on the test set. We then report, for each method, the mean test AUC over the 10 repeats, together with a 95% confidence interval defined a mean ± 1.96 × standard error of the mean. We used Scikit Learn’s implementation of logistic regression regularized by `1 (lasso) and `2 (ridge) penalty from sklearn.linear_model.LogisticRegression, and varied the C regularization parameter over a grid of 10 values: {10−5, 10−3, . . . , 104}. For the unsupervised dictionary learning experiment method, we estimated a dictionary from the gene expression data in the training and validation sets, using sklearn.decomposition.DictionaryLearning(n_components=10, alpha=2.0), which produces sparse codes in k = 10 dimensions with roughly 50% nonzero coefficients by minimizing the squared Frobenius reconstruction distance with lasso regularization on the code. We then use sklearn.linear_model.LogisticRegression to train a logistic regression on the codes, varying the ridge regularization parameter C over a grid of 10 values {10−1, 100, . . . , 108}. Finally, we implemented the task-driven dictionary learning model (13) with our toolbox, following the pseudo-code in Figure 10. Like for the unsupervised dictionary learning experiment, we set the dimension of the codes to k = 10, and a fixed elastic net regularization on the inner optimization problem to ensure that the codes have roughly 50% sparsity. For the outer optimization problem, we solve an `2 regularized ridge regression problem, varying again the ridge regularization parameter C over a grid of 10 values {10−1, 100, . . . , 108}. Because the outer problem is non-convex, we minimize it using the Adam optimizer (Kingma & Ba, 2014) with default parameters.

heading: E.3 DATASET DISTILLATION
text: Experimental setup. For the inner problem, we used gradient descent with backtracking linesearch, while for the outer problem we used gradient descent with momentum and a fixed step-size. The momentum parameter was set to 0.9 while the step-size was set to 1.
Figure 4 was produced after 4000 iterations of the outer loop on CPU (Intel(R) Xeon(R) Platinum P-8136 CPU @ 2.00GHz), which took 1h55. Unrolled differentiation took instead 8h:05 (4 times more) to run the same number of iterations. As can be seen in Figure 16, the output is the same in both approaches.

heading: E.4 MOLECULAR DYNAMICS
text: Our experimental setup is adapted from the JAX-MD example notebook available at https://github.com/google/jax-md/blob/master/notebooks/meta_ optimization.ipynb.
We emphasize that calculating the gradient of the total energy objective, f(x, θ) = ∑ ij U(xi,j , θ), with respect to the diameter θ of the smaller particles,∇1f(x, θ), does not require implicit differentiation or unrolling. This is because∇1f(x, θ) = 0 at x = x?(θ):
∇θf(x?(θ), θ) = ∂x?(θ)>∇1f(x?(θ), θ) +∇2f(x?(θ), θ) = ∇2f(x?(θ), θ).
This is known as Danskin’s theorem or envelope theorem. Thus instead, we consider sensitivities of position ∂x?(θ) directly, which does require implicit differentiation or unrolling.
Our results comparing implicit and unrolled differentiation for calculating the sensitivity of position are shown in Figure 17. We use BiCGSTAB (Vorst & van der Vorst, 1992) to perform the tangent linear solve. Like in the original JAX-MD experiment, we use k = 128 particles in m = 2 dimensions.

Review summaries:
This paper provides a unified tool for combining the implicit differentiation technique and the  automatic differentiation method widely used in existing deep learning packages such as PyTorch and TensorFlow. The proposed implementation is easy to use for numeric optimization such as bilevel optimization, meta-learning and hyperparameter optimization because it covers many existing schemes such as fixed point, KKT point, projected method. In the experiments, the authors illustrate how their tool can be useful to simply the implementation. 

The paper proposes a modular and efficient framework along with its JAX implementation for the implicit differentiation of optimization problems. The user defines the function F capturing the optimality conditions of the problem to be differentiated; then the framework combines implicit differentiation and autodiff of F to automatically differentiate the optimization problem. The proposed framework is labeled as efficient, since it doesn’t have to unroll the computational graph like in autodiff, and modular since it doesn’t require case-by-case mathematical derivation like in implicit differentiation. The authors show that existing implicit differentiation methods can be instantiated in their framework. They provide and empirically validate new bounds on the Jacobian error when the optimization problem is only solved approximately. The authors implemented four illustrative applications of their framework ( Hyperparameter Optimization Of Multiclass SVM; Dataset Distillation; Task-Driven Dictionary Learning; Sensitivity Analysis Of Molecular Dynamics). Code and implementation in JAX are provided along with the paper. 

This paper introduces a Jax package for implicitly differentiating various numerical solvers. Concretely, the authors develop a systemic methodology for producing gradients for a variety of optimization problems. Then, the authors prove that the Jacobian solution to the approximate numerical solution produces close enough gradients. Finally, the authors show the power of their framework on four test tasks.

A good paper considers a critical problem, but contributions are somewhat weak.  This paper aims at a unified approach for implicit differentiation of optimization problems. The experiments results seem good, but the proposed methods lack novelty.

This paper presents a module for implicit differentiation in JAX, which allows efficient differentiation through optimization and root-finding problems. The paper presents the general principle of combining the implicit function theorem with autodiff of optimality conditions, and discusses how this can work for a implicit functions such as a range of constrained optimization methods. Finally there are example applications to three machine learning tasks and one molecular dynamics simulation task where this framework improves on currently used methods.

The paper promises extension of Google’s JAX library by automatic and modular differentitation of implicitly defined functions. Given a solver $\theta\mapsto x^*(\theta)$ one asks for its Jacobian $\partial x^*(\theta)$. In cases when $x^*$ does not have an explicit form that can be automatically differentiated, the Jacobian can be computed from the optimality conditions $F$ having the general form $F(x^*(\theta),\theta)=0$. Jacobian $\partial x^*(\theta)$ is then given by a solution of a linear system involving $\partial_1 F$ and $\partial_2 F$. The authors describe a unified implementation in JAX and show the resulting simple interface for end-users. Various use cases are presented on examples and code-snippets figures, including stationary point conditions, KKT conditions, proximal gradient, projective gradient and mirror descent fixed point and more. On the theoretical side, the authors discuss the case when the optimization problem solution is only approximate and derive Jacobian error bounds under certain assumptions on optimality conditions $F$.

The authors propose a unified modular framework for implicit differentiation, which removes the need to manually derive the jacobians for each individual case by leveraging automatic differentiation of the optimality conditions. The proposed framework is shown to encompass various existing implicit differentiation schemes (such as stationary point- or KKT conditions) while also creating new ones (e.g. mirror descent fixed point). To handle the practical case in which the optimization problem can only be solved approximately, a theoretically justified bound on the error of the jacobian estimate is provided, showing that the error on the estimated jacobian is at most of the order of the deviation of the approximate solution from the optimal solution. Additionally, a JAX implementation of the proposed framework is described, with the promise of an open-source JAX library for implicit differentiation upon acceptance. Finally, four applications of the framework are demonstrated, highlighting its flexibility and ease of use.