This came up in a conversation with @seantalts.

When solving an algebraic equation, we find y^* such that f(y^*, \theta) = 0 where \theta are our model parameters. We compute the jacobian matrix of y with respect to \theta using the implicit function theorem:

J = - [J^y(y^*, \theta)]^{-1} J^\theta(y^*, \theta)

where J^y is the Jacobian matrix of f with respect to y, etc.

Currently, I compute J^y and J^\theta separately using `Jacobian()`

, but it could be more efficient to compute \tilde J := [J^y, J^\theta] in one go. The cost should be roughly halved, since it’s the same amount of sweeps in reverse-mode autodiff.

Any objections to this idea? If not, I’ll submit an issue, fix it, and do the PR.