  We consider the problem of efficiently computing the derivative of the
solution map of a convex cone program, when it exists. We do this by implicitly
differentiating the residual map for its homogeneous self-dual embedding, and
solving the linear systems of equations required using an iterative method.
This allows us to efficiently compute the derivative operator, and its adjoint,
evaluated at a vector. These correspond to computing an approximate new
solution, given a perturbation to the cone program coefficients (i.e.,
perturbation analysis), and to computing the gradient of a function of the
solution with respect to the coefficients. Our method scales to large problems,
with numbers of coefficients in the millions. We present an open-source Python
implementation of our method that solves a cone program and returns the
derivative and its adjoint as abstract linear maps; our implementation can be
easily integrated into software systems for automatic differentiation.
