Sparse Jacobians

A dense Jacobian with computeJacobian (XAD/Jacobian.hpp) takes one forward-mode pass per input, or one adjoint sweep per output. When most entries of the Jacobian are known to be zero, far fewer passes are enough: columns that never affect the same row can be seeded together in one forward-mode pass, and their contributions stay separate. This is called column colouring. It needs no special API: it uses the forward mode type and derivative(), as described in the C++ tutorial.

When to Use It

Use it when each output depends on only a few inputs and the sparsity pattern is known, for example banded Jacobians from discretised differential equations or local interactions. The number of passes is the number of colours, which depends on the pattern and not on the number of inputs. For dense Jacobians, use computeJacobian.

Example: a Tridiagonal Jacobian

The function below, with x_{-1} = x_n = 0, has outputs that depend only on their neighbours:

template <typename T>
std::vector<T> f(const std::vector<T>& x)
{
    using std::exp;
    const std::size_t n = x.size();
    std::vector<T> y(n);
    for (std::size_t i = 0; i < n; ++i)
    {
        T left = (i > 0) ? x[i - 1] : T(0.0);
        T right = (i + 1 < n) ? x[i + 1] : T(0.0);
        y[i] = left - 2.0 * x[i] + right + left * right + 0.1 * exp(x[i]);
    }
    return y;
}

Row i depends on columns i-1, i and i+1, which have different residues modulo 3. Giving column j the colour j % 3, three forward passes compute the whole Jacobian, whatever its size:

typedef xad::fwd<double>::active_type AD;
const std::size_t n = x0.size();
std::vector<double> lower(n), diag(n), upper(n);   // the three bands
for (std::size_t colour = 0; colour < 3; ++colour)
{
    // seed all the columns of this colour at once
    std::vector<AD> x(x0.begin(), x0.end());
    for (std::size_t j = 0; j < n; ++j)
        derivative(x[j]) = (j % 3 == colour) ? 1.0 : 0.0;
    std::vector<AD> y = f(x);
    // only one of the columns i-1, i, i+1 has this colour,
    // so derivative(y[i]) is exactly that one entry of row i
    for (std::size_t i = 0; i < n; ++i)
    {
        double d = derivative(y[i]);
        if (i > 0 && (i - 1) % 3 == colour)
            lower[i] = d;
        else if (i % 3 == colour)
            diag[i] = d;
        else if (i + 1 < n && (i + 1) % 3 == colour)
            upper[i] = d;
    }
}

With the vector forward mode type xad::fwd<double, 3>, all colours can be computed in a single pass by seeding derivative(x[j])[j % 3] = 1.0.

Sample

The SparseJacobian sample contains the full code and checks the result against the dense Jacobian from computeJacobian in adjoint mode.