When the Jacobian Is Dense but the Program Is Not

23 October 2026 · Download PDF · Download code · Marimo notebook

In machine learning, automatic differentiation is mostly used for products with the Jacobian: the gradient of a loss is a single vector-Jacobian product, and training never forms the Jacobian itself. Many methods in scientific computing need the matrix. Newton’s method for a nonlinear system $𝐹\left(𝑥\right)=0$, at the heart of implicit time stepping for stiff differential equations and of steady-state solvers for discretized partial differential equations, factorizes the Jacobian in every step, and solvers for constrained optimization, such as Ipopt [1] and Uno [2], take the Jacobian of the constraints and the Hessian of the Lagrangian as matrices. These Jacobians are large, with thousands to millions of unknowns, but more often than not sparse, since each equation involves only a few unknowns.

An Illustrated Guide to Automatic Sparse Differentiation [3] explains how standard automatic differentiation tools can be used to compute a sparse Jacobian at a fraction of the cost of a dense one. The idea is fairly simple: first detect which entries can be nonzero, group the columns (or rows) that never share a row (or column), and compute Jacobian-vector (or vector-Jacobian) products in order to generate the sparse Jacobian.¹¹ The grouping is for efficiency: all columns of a group are computed by a single Jacobian-vector product, and all rows of a group by a single vector-Jacobian product. In Backpropagation as a Linear Solve we derived that the Jacobian of any computation can be computed as

where $𝑃$, $𝑀$, $𝐿$ are matrices which more often than not are sparse and with $𝐿$ further being unit lower triangular. Even when $𝐽$ is dense, products with it are cheap: $𝐽𝑣=𝑃\left({𝐿}^{-1}\left({𝑀}^{\top }𝑣\right)\right)$ is evaluated from right to left, with one forward substitution in the sparse, unit lower triangular $𝐿$, at about the cost of running the program.²² Likewise for the transpose: ${𝐽}^{\top }𝑤=𝑀\left({𝐿}^{-\top }\left({𝑃}^{\top }𝑤\right)\right)$ needs one backward substitution in the unit upper triangular ${𝐿}^{\top }$, which is reverse mode. What a dense $𝐽$ makes expensive is forming it, or solving with it.

This post asks what the matrices $𝐿$, $𝑀$, and $𝑃$ add to the understanding of the Jacobian. The answer in short is that they pay off when $𝐿$ is sparse but ${𝐿}^{-1}$ is dense, resulting in a dense Jacobian $𝐽$ (e.g., Example 1).

1. From program to matrices

We can create $𝐿$, $𝑀$, and $𝑃$ automatically by recording the program. A tape stores the partial derivatives of every operation, evaluated as the program runs; the intermediate values are only needed during the run, to compute those partial derivatives and the outputs. The tape of this post is deliberately small. A Tracked number holds a value and the index of its node. Arithmetic on Tracked numbers is overloaded, as for dual numbers: a * b computes the product, gives it the next node index, and stores its partial derivatives, b with respect to a and a with respect to b, as entries of $𝐿$ or ${𝑀}^{\top }$, depending on whether a and b are nodes or inputs. The function trace runs $𝑓$ on Tracked inputs and returns $𝑓\left(𝑥\right)$ together with the sparse matrices $𝐿$, ${𝑀}^{\top }$, and $𝑃$.

1
# A scalar tape. Every operation on a Tracked number is one node Φᵢ
2
# of the forward relation u = Φ(u, x), and it stores its partial
3
# derivatives, the entries L[i, j] = -∂Φᵢ/∂uⱼ and Mᵀ[i, k] = ∂Φᵢ/∂xₖ.
4
using LinearAlgebra, SparseArrays
5
6
# The tape: rows, columns and values of the entries of L - I and Mᵀ,
7
# and the number of nodes
8
const L_entries = (Int[], Int[], Float64[])
9
const Mt_entries = (Int[], Int[], Float64[])
10
const nodes = Ref(0)
11
12
# index > 0 is the node uᵢ, index < 0 is the input x₋ᵢ
13
struct Tracked
14
    value::Float64
15
    index::Int
16
end
17
18
# Append the node with the given value and (argument => partial) pairs
19
function record!(value, args...)
20
    i = nodes[] += 1
21
    for (a, ∂a) in args
22
        rows, cols, vals = a.index > 0 ? L_entries : Mt_entries
23
        push!(rows, i)
24
        push!(cols, abs(a.index))
25
        push!(vals, a.index > 0 ? -∂a : ∂a)
26
    end
27
    return Tracked(value, i)
28
end
29
30
Base.:+(a::Tracked, b::Tracked) = record!(a.value + b.value, a => 1.0, b => 1.0)
31
Base.:-(a::Tracked, b::Tracked) = record!(a.value - b.value, a => 1.0, b => -1.0)
32
Base.:*(a::Tracked, b::Tracked) = record!(a.value * b.value, a => b.value, b => a.value)
33
Base.:/(a::Tracked, b::Tracked) =
34
    record!(a.value / b.value, a => 1 / b.value, b => -a.value / b.value^2)
35
Base.:+(a::Tracked, c::Real) = record!(a.value + c, a => 1.0)
36
Base.:+(c::Real, a::Tracked) = a + c
37
Base.:-(a::Tracked, c::Real) = a + (-c)
38
Base.:-(c::Real, a::Tracked) = record!(c - a.value, a => -1.0)
39
Base.:-(a::Tracked) = record!(-a.value, a => -1.0)
40
Base.:*(a::Tracked, c::Real) = record!(a.value * c, a => Float64(c))
41
Base.:*(c::Real, a::Tracked) = a * c
42
Base.:/(a::Tracked, c::Real) = a * inv(c)
43
Base.:/(c::Real, a::Tracked) = record!(c / a.value, a => -c / a.value^2)
44
Base.exp(a::Tracked) = (v = exp(a.value); record!(v, a => v))
45
Base.sin(a::Tracked) = record!(sin(a.value), a => cos(a.value))
46
Base.:^(a::Tracked, p::Integer) = record!(a.value^p, a => p * a.value^(p - 1))
47
48
# Run f on tracked inputs and return f(x) with the sparse matrices
49
# L = I - ∂Φ/∂u, Mᵀ = ∂Φ/∂x and the selection P of the outputs, so
50
# that the Jacobian of f at x is P L⁻¹ Mᵀ.
51
function trace(f, x)
52
    foreach(empty!, (L_entries..., Mt_entries...))
53
    nodes[] = 0
54
    y = f([Tracked(xk, -k) for (k, xk) in enumerate(x)])
55
    d, n, m = nodes[], length(x), length(y)
56
    L = I + sparse(L_entries..., d, d) # sums the entries of x * x
57
    Mt = sparse(Mt_entries..., d, n)
58
    P = sparse(1:m, [yi.index for yi in y], ones(m), m, d)
59
    return [yi.value for yi in y], L, Mt, P
60
end

With the matrices at hand, the Jacobian itself is a triangular solve with a sparse right-hand side, $𝐽=𝑃\left({𝐿}^{-1}{𝑀}^{\top }\right)$. Column $𝑘$ of ${𝐿}^{-1}{𝑀}^{\top }$ is nonzero only at the intermediate variables that depend on ${𝑥}_{𝑘}$, and a sparse triangular solve finds those by a search through the graph of $𝐿$ before it does any arithmetic. SparseSparse.jl does this when both sides of \ are sparse, whereas the standard library returns a dense matrix.

1
using SparseSparse # Required for \ with a sparse right-hand side
2
include("tape.jl")
3
4
# u''(xᵢ) by central differences, with u = 0 at both ends
5
function cdiff2(x, i, h)
6
    left = i > 1 ? x[i-1] : 0.0
7
    right = i < length(x) ? x[i+1] : 0.0
8
    return (left - 2x[i] + right) / h^2
9
end
10
11
# Bratu, u'' + λ exp(u) = 0: each output reads three inputs
12
function bratu(x; λ=1.0)
13
    h = 1 / (length(x) + 1)
14
    return [cdiff2(x, i, h) + λ * exp(x[i]) for i in eachindex(x)]
15
end
16
17
# Nonlocal Ohmic heating, u'' + λ exp(-u) / (∫₀¹ exp(-u) dx)² = 0, of the
18
# next section: every output reads the sum s, so every input
19
function heating(x; λ=2.0)
20
    h = 1 / (length(x) + 1)
21
    s = sum(exp(-xj) for xj in x) # h * s is the integral
22
    scaling = λ / (h * s)^2
23
    return [cdiff2(x, i, h) + scaling * exp(-x[i]) for i in eachindex(x)]
24
end
25
26
x = 0.1 .* rand(1000)
27
_, L, Mt, P = trace(bratu, x) # nnz(L) = 13000, nnz(M) = 3998
28
J = P * (L \ Mt) # tridiagonal: nnz(J) = 2998
29
_, L, Mt, P = trace(heating, x) # nnz(L) = 27997, nnz(M) = 4998
30
J = P * (L \ Mt) # dense: nnz(J) = 1000000

2. The augmented system

While both the Bratu and the heating residual (Example 3) in Code 2 give sparse $𝑃$, $𝐿$, and $𝑀$, only the Bratu Jacobian is sparse. The cumulative sum of Example 1 already showed how this can happen: a sparse $𝐿$ can have a dense inverse, in which case sparse differentiation cannot help. For products, a dense $𝐽$ is no problem, since $𝐽𝑣=𝑃\left({𝐿}^{-1}\left({𝑀}^{\top }𝑣\right)\right)$ costs one forward substitution. Newton’s method for $𝐹\left(𝑥\right)=0$, however, needs a solve, $𝐽\mathrm{\Delta }=𝐹\left(𝑥\right)$, in every step, and with a dense $𝐽$ that means forming $𝐽$ and factorizing it, at a cost of order ${𝑛}^3$. The cumulative sum also hinted at the way around this: given $𝑃$, $𝐿$, and $𝑀$, the solve can be written as one sparse augmented system,

whose first block row gives $𝑧={𝐿}^{-1}{𝑀}^{\top }\mathrm{\Delta }$ and second $𝑃𝑧=𝐽\mathrm{\Delta }=𝐹\left(𝑥\right)$. Newton’s method can therefore solve with the sparse $𝐴$, whose unknowns are $𝑧$ and the step $\mathrm{\Delta }$, instead of with the dense $𝐽$.⁵⁵ The same factorization of $𝐴$ also solves with ${𝐽}^{\top }$, as in adjoint methods: ${𝐽}^{\top }𝑤=𝑔$ is $\left[\begin{array}{cc}{𝐿}^{\top }& {𝑃}^{\top }\\ -𝑀& 0\end{array}\right]\left[\begin{array}{c}𝜆\\ 𝑤\end{array}\right]=\left[\begin{array}{c}0\\ 𝑔\end{array}\right]$, the system with ${𝐴}^{\top }$.

A = [L -Mt; P spzeros(m, n)]
Δ = (A \ [zeros(d); F])[d+1:end]

3. Lifting

This could seem like the end for using standard AD software for getting the sparse augmented Jacobian, however that is not true. A simple technique called Lifting [5] instead rewrites the program, so that the intermediate variables become unknowns. That is, the lifted program returns both $𝑢$ and $𝑓\left(𝑥\right)$, and the forward relation and the outputs become equations. Said differently: The lifted program has a Jacobian that is exactly the augmented matrix $𝐴$ of (11) ⁹⁹ $𝐺\left(𝑢,𝑥\right)=\left(0,𝑓\left(𝑥\right)\right)$ holds only at $𝑢=𝑢\left(𝑥\right)$, the values the program computes; for other $𝑢$, the first block measures how far $𝑢$ is from satisfying the forward relation.

Now the Jacobian of $𝐺$ is sparse, even when ${𝐿}^{-1}$ is dense. That means that tools such that DifferentiationInterface.jl or adsex can be used to extract the sparse Jacobian. By hand, lifting every operation would mean rewriting the whole program as equations, so in practice only the few intermediate variables that couple many others are lifted, such as the integral in the heating problem.

Example 6: Estimating the heating parameter

Example 6 (Estimating the heating parameter)   The same idea regarding lifting can be used in mathematical programming software like JuMP.jl. In this example try to estimate $𝜆$ based on a set of noisy temperature measurements, using the UnoSolver, that is solving the following nonlinear program

$$\begin{array}{cc}\underset{𝑥,𝜆\ge 0}{min}& \sum _{𝑖=1}^{𝑛}{\left({𝑥}_{𝑖}-{𝑦}_{𝑖}\right)}^2\\ \text{subject to}& {𝑥}_{𝑖-1}-2{𝑥}_{𝑖}+{𝑥}_{𝑖+1}+\frac{𝜆exp\left(-{𝑥}_{𝑖}\right)}{{𝑠}^2}=0,𝑖=1,\dots ,𝑛.\end{array}$$

(17)

In JuMP.jl, the sum $𝑠={\sum }_{𝑗}exp\left(-{𝑥}_{𝑗}\right)$ can be written in two ways. Directly defining s as an expression

s = @expression(model, sum(exp(-x[j]) for j in 1:n))

and lifted, it is a variable with a constraint of its own,

@variable(model, s, start = n) # sum(exp(-x[j])) = n at the start x = 0
@constraint(model, s == sum(exp(-x[j]) for j in 1:n))

The expression is not a fixed number but a shorthand, which JuMP.jl substitutes into every constraint that uses s. As a variable, $𝑠$ is an unknown of its own, which Uno updates along with $𝑥$, and its constraint makes it equal the sum only at the solution. With the expression, every constraint reads every ${𝑥}_{𝑗}$, so the Jacobian and the Hessian that JuMP.jl hands to Uno become dense. Lifted, only the constraint of $𝑠$ reads all of $𝑥$, and since JuMP.jl differentiates each constraint on its own, that one full row is cheap and both matrices are sparse.

Table 4: Solving the estimation with $𝑠$ lifted and as an expression, with Uno through JuMP.jl; both estimate $𝜆\approx 2.00$. The memory is that allocated in Julia.

$𝑛$	Lifted		Expression
	Time	Memory	Time	Memory
100	5 ms	2.6 MiB	3.0 s	809 MiB
200	28 ms	5.2 MiB	10.6 s	4.2 GiB
1000	36 ms	25 MiB	—	—
2000	140 ms	50 MiB	—	—
4000	192 ms	102 MiB	—	—

Figure 1: The estimation for $𝑛=50$ with the true $𝜆=2$. The data have noise with standard deviation 0.03 added.

The script is heating_lifted.jl.

4. Summary

Every Jacobian of a program is the product $𝐽=𝑃{𝐿}^{-1}{𝑀}^{\top }$. Even when the three matrices are sparse, the inverse ${𝐿}^{-1}$ can be dense, resulting in a dense Jacobian. Sparse differentiation pays off when $𝐽$ is sparse, but not when it is dense. In this case actually knowing the structure of the Jacobian can be useful. Either directly knowing $𝑃$, $𝐿$, and $𝑀$ gives us a sparse representation of the Jacobian, or lifting the program to expose intermediate variables gives a sparse Jacobian which can be connected to the original dense Jacobian through a Schur complement.¹¹¹¹ Eliminating the intermediate variables from the lifted Jacobian leaves the Schur complement of the upper-left block $𝐿$ in $\left[\begin{array}{cc}𝐿& -{𝑀}^{\top }\\ 𝑃& 0\end{array}\right]$, which is $0-𝑃{𝐿}^{-1}\left(-{𝑀}^{\top }\right)=𝑃{𝐿}^{-1}{𝑀}^{\top }=𝐽$.

@misc{paltorp2026sparsejacobians,
  author = {Paltorp, Mikkel},
  title = {{When the Jacobian Is Dense but the Program Is Not}},
  year = {2026},
  month = oct,
  url = {https://mipals.github.io/pdf/2026/10-23-sparse-jacobians/sparse-jacobians.pdf},
  note = {Blog post}
}