Support for sparse jacobians

Description

I really want to be able to call pt.jacobian and get back a sparse matrix. For very large systems of equations, the jacobian is essentially never going to be dense, so this would lead to huge performance gains.

I talked with @aseyboldt about it informally once, and he thought it was doable. There is a jax package for it, so that's nice. I confess I don't fully understand what these packages are doing under the hood. From the sparsejac function in the linked package:

    This function uses reverse-mode automatic differentiation to compute the
    Jacobian. The `fn` must accept a rank-1 array and return a rank-1 array, and
    the Jacobian should be sparse with nonzero elements identified by `sparsity`.
    Sparsity is exploited in order to make the Jacobian computation efficient.

    This is done by identifying "structurally independent" groups of output
    elements, which is isomorphic to a graph coloring problem. This allows
    project to a lower-dimensional output space, so that reverse-mode
    differentiation can be more efficiently applied.

Seems like we should be able to do that too!

Start with the pt.jacobian entry point and read the linked sparsejac implementation, especially its description of reverse-mode differentiation and sparsity grouping. Done means supporting sparse Jacobian results for large systems while exploiting the stated sparsity structure; the issue does not identify project files or tests.