AutoSparse(AutoFiniteDiff()) inaccurate for variables of different magnitudes #1085

Goal

We compute Jacobians and Hessians of nonlinear models by finite differences, for use in nonlinear model-predictive control. They are sparse, so we would like to use AutoSparse(AutoFiniteDiff()). The models are badly scaled in the sense that the variables differ in magnitude by several orders, for example a pressure of order 1e5 Pa next to a variable of order one or smaller.

Observation

jacobian with AutoFiniteDiff() calls FiniteDiff.finite_difference_jacobian, which takes a step per column relative to that column's variable, and is accurate. jacobian with AutoSparse(AutoFiniteDiff()) computes the compressed Jacobian through pushforward with the seed of each color, which calls FiniteDiff.finite_difference_jvp! (ext/DifferentiationInterfaceFiniteDiffExt/twoarg.jl#L111). There the step is max(relstep * norm(x), absstep) / norm(v) (FiniteDiff src/jvp.jl#L99-L105), i.e. relative to the norm of the whole point. For a small variable next to a large one, this step is orders of magnitude larger than the step of the dense Jacobian, and the sparse Jacobian is wrong in the leading digits. With central differences, where the relative step is larger, the result for the small variable is entirely wrong.

The sparse Hessian is computed through hvp, i.e. finite differences of a finite-difference gradient, with the same norm-relative outer step, and is likewise inaccurate, whereas hessian with AutoFiniteDiff() (which calls FiniteDiff.finite_difference_hessian) is accurate. The forward-over-forward part of this is #1012; the example below shows that the sparse Hessian is also affected by the scaling of the variables.

Example

The Jacobian and the Hessian are diagonal, so both columns are given one color.

using DifferentiationInterface, FiniteDiff, SparseConnectivityTracer, SparseMatrixColorings, LinearAlgebra

f!(y, x) = (y[1] = x[1]^3; y[2] = x[2]^3; y)   # diagonal Jacobian, both columns in one color
Jexact(x) = Diagonal(3 .* x .^ 2)

for x in ([1e5, 1e-2], [1e5, 1.0]), fdtype in (Val(:forward), Val(:central))
    dense = AutoFiniteDiff(; fdtype, fdjtype = fdtype)
    sparse_backend = AutoSparse(dense; sparsity_detector = TracerSparsityDetector(),
        coloring_algorithm = GreedyColoringAlgorithm())
    y = zeros(2)
    Jd = jacobian(f!, y, dense, x)
    Js = jacobian(f!, y, sparse_backend, x)
    println("x = $x, $fdtype")
    println("  AutoFiniteDiff           relative error: ", diag(Jd) ./ diag(Jexact(x)) .- 1)
    println("  AutoSparse(AutoFiniteDiff) relative error: ", diag(Js) ./ diag(Jexact(x)) .- 1)
end

# Hessian of a separable function, scaled so that both terms are of order one
g(x) = (x[1] / 1e5)^3 + x[2]^3
Hexact(x) = Diagonal(6 .* [x[1] / 1e15, x[2]])
x = [1e5, 1e-2]
sb = AutoSparse(AutoFiniteDiff(); sparsity_detector = TracerSparsityDetector(),
    coloring_algorithm = GreedyColoringAlgorithm())
println("Hessian, x = $x")
println("  AutoFiniteDiff           relative error: ", diag(hessian(g, AutoFiniteDiff(), x)) ./ diag(Hexact(x)) .- 1)
println("  AutoSparse(AutoFiniteDiff) relative error: ", diag(hessian(g, sb, x)) ./ diag(Hexact(x)) .- 1)

Output:

x = [100000.0, 0.01], Val{:forward}()
  AutoFiniteDiff           relative error: [1.495313073363036e-8, 1.4900886828872473e-6]
  AutoSparse(AutoFiniteDiff) relative error: sparsevec([1, 2], [7.650372779721692e-9, 0.10906786469265284], 2)
x = [100000.0, 0.01], Val{:central}()
  AutoFiniteDiff           relative error: [1.7613910330283034e-11, 1.2222844958387213e-7]
  AutoSparse(AutoFiniteDiff) relative error: sparsevec([1, 2], [1.100008972798605e-12, 611.142143750178], 2)
x = [100000.0, 1.0], Val{:forward}()
  AutoFiniteDiff           relative error: [1.495313073363036e-8, 1.4901161193847656e-8]
  AutoSparse(AutoFiniteDiff) relative error: sparsevec([1, 2], [7.600377660565982e-9, 0.0010540412872301896], 2)
x = [100000.0, 1.0], Val{:central}()
  AutoFiniteDiff           relative error: [1.7613910330283034e-11, 8.275380380950992e-12]
  AutoSparse(AutoFiniteDiff) relative error: sparsevec([1, 2], [1.4356515976032824e-11, 0.061114214381128296], 2)
Hessian, x = [100000.0, 0.01]
  AutoFiniteDiff           relative error: [1.2417633588057697e-9, 3.9736429924275285e-8]
  AutoSparse(AutoFiniteDiff) relative error: sparsevec([1, 2], [-0.5571909695914141, 0.05264629492636863], 2)

What we tried

• relstep and absstep of AutoFiniteDiff: these set one scalar step for the whole point, so no choice is accurate for both a variable of order 1e5 and a variable of order 1e-2.
• Forward and central differences: both are affected.
• FiniteDiff's own colored Jacobian (finite_difference_jacobian! with colorvec and sparsity) also takes one step per color and is inaccurate for this example as well; reported in Colored Jacobian uses one step per color, inaccurate for variables of different magnitudes FiniteDiff.jl#241.
• As a workaround we wrote our own colored finite differences, in which every column of a color is perturbed by the step that the dense Jacobian takes for that column, and each entry is divided by the step of its column.

Question

Is there an intended way to obtain sparse finite-difference Jacobians and Hessians through DifferentiationInterface whose accuracy matches the dense AutoFiniteDiff() results when the variables differ in magnitude? If not, is the use of a JVP with a norm-relative step for the compressed Jacobian a deliberate design choice?

Versions

Julia 1.12.7, DifferentiationInterface v0.7.21, FiniteDiff v2.33.0, SparseConnectivityTracer v1.2.4, SparseMatrixColorings v0.4.29, ADTypes v1.24.0.

Host: demeter2 · Session: c59a36f8-c4c7-4a84-9858-91adc076743c · Working directory: ~/.julia/dev/MPCComponents

🤖 Generated with Claude Code