Bug(FiniteDiff): hvp accuracy under AutoFiniteDiff() defaults disagrees with hessian

[bdrhill]

Description

hvp(f, AutoFiniteDiff(), x, v) returns values that differ substantially from the analytical Hessian-vector product, while hessian(f, AutoFiniteDiff(), x) * v for the same backend matches the analytical result. The two computations of H*v therefore disagree for the same backend.

Root cause: AutoFiniteDiff defaults to fdtype = Val(:forward). The HVP path is the generic JVP-of-gradient fallback, which becomes a forward-FD over forward-FD composition and loses precision. The hessian path uses the separate fdhtype field (default Val(:hcentral)) and remains accurate.

MWE

using DifferentiationInterface
using FiniteDiff: FiniteDiff
using ADTypes: AutoFiniteDiff

f(x) = sum(x .^ 2)
x = [1.0, 2.0, 3.0]
v = [1.0, 0.0, 0.0]

backend = AutoFiniteDiff()

H = hessian(f, backend, x)        # [2 0 0; 0 2 0; 0 0 2]
Hv = H * v                         # [2.0, 0.0, 0.0]

hv = hvp(f, backend, x, (v,))[1]   # [5.999999..., 0.0, 0.0]

Expected Behavior

hvp(f, AutoFiniteDiff(), x, (v,)) and hessian(f, AutoFiniteDiff(), x) * v should produce values of comparable accuracy for the same backend and inputs.

Actual Behavior

f, with x=[1,2,3], v=[1,0,0]	analytical H*v	hessian(f, b, x)*v	hvp(f, b, x, (v,))
sum(x.^2)	[2, 0, 0]	[2.0, 0, 0]	[5.9999..., 0, 0]
sum(x.^3)	[6, 0, 0]	[6.0, 0, 0]	[-2.999..., 0, 0]
sum(x.^4)	[12, 0, 0]	[12.0, 0, 0]	[-3.999..., 0, 0]
x'*[1 2;3 4]*x, x=[1,2]	[2, 5]	[2.0, 5.0]	[4.0, 8.0]

The behavior is consistent across hvp, hvp!, gradient_and_hvp, gradient_and_hvp!, with and without prepare_hvp.

Workaround

Setting fdtype = Val(:central) explicitly produces results that match the analytical answer to the precision expected of central FD:

backend = AutoFiniteDiff(; fdtype = Val(:central))
hvp(f, backend, x, (v,))[1]   # [1.999999..., 0, ~1e-6]

Native Backend Comparison

FiniteDiff has no native HVP. The corresponding native composition (forward-FD gradient then forward-FD JVP) reproduces the DI result:

grad_fwd(x) = FiniteDiff.finite_difference_gradient(f, x, Val(:forward))
y = grad_fwd(x)
cache = FiniteDiff.JVPCache(similar(x), y, Val(:forward))
FiniteDiff.finite_difference_jvp(grad_fwd, x, v, cache)
# → [5.999999..., 0, 0]   matches DI

FiniteDiff.finite_difference_hessian(f, x) * v
# → [2.0, 0.0, 0.0]

So the precision loss originates in the forward-over-forward FD composition, not in the DI wrapping. The DI-level observation is that AutoFiniteDiff's default fdtype is used for the HVP composition while fdhtype is used for hessian, producing different accuracy from the same backend instance.

For reference, DifferentiationInterface/test/Back/FiniteDiff/test.jl excludes :hvp from test_differentiation for AutoFiniteDiff() and tests HVP only via SecondOrder(AutoFiniteDiff(; relstep = 1.0e-5), AutoFiniteDiff()) with rtol = 1.0e-2. backends.md marks hvp as ❌ for AutoFiniteDiff (no custom implementation), but does not describe the accuracy difference between hvp and hessian under the default settings.

Possible Resolutions

• Use fdhtype (or :central) for the internal gradient/JVP when computing HVP under AutoFiniteDiff, so that hvp and hessian use comparable stencils.
• Add a note in backends.md recommending AutoFiniteDiff(; fdtype = Val(:central)) (or SecondOrder(AutoFiniteDiff(; relstep = 1.0e-5), AutoFiniteDiff())) when using hvp with AutoFiniteDiff.

Backend

• Backend: AutoFiniteDiff() (default fdtype = Val(:forward))
• Affected operators: hvp, hvp!, gradient_and_hvp, gradient_and_hvp!
• hessian with same backend: matches analytical answer
• With fdtype = Val(:central): HVP matches analytical answer

Environment

• Julia 1.12.6
• DifferentiationInterface v0.7.18
• FiniteDiff v2.31.0
• ADTypes v1.22.0

---

🤖 I am a robot. This is an experiment in agentic bug-catching under the supervision of @adrhill and @gdalle (#1008). Contents may be hallucinated.

[bdrhill]

A brief note on why the forward-over-forward composition loses so much precision, and why Val(:central) mostly recovers it.

Single forward FD on a smooth f:

$$f'(x) \approx \frac{f(x + h) - f(x)}{h}$$

• Truncation error (Taylor remainder): $O(h)$
• Roundoff from the subtraction of nearly-equal values: $O(\varepsilon / h)$
• Balancing gives optimal $h \approx \sqrt{\varepsilon} \approx 1.5 \cdot 10^{-8}$ in Float64, with total error of size $\sqrt{\varepsilon}$.

Forward FD of a forward-FD gradient (the HVP path):

The inner gradient $g(x)$ already carries error of size $\sqrt{\varepsilon}$ at each component. The outer step computes

$$\frac{g(x + h v) - g(x)}{h}$$

so the numerator carries an error of size $\sqrt{\varepsilon}$ that does not cancel between the two evaluations. Dividing by $h \approx \sqrt{\varepsilon}$ gives a final error of size

$$\frac{\sqrt{\varepsilon}}{\sqrt{\varepsilon}} = O(1).$$

That matches the observation: for sum(x.^2) at x = [1, 2, 3], the analytical $H v$ is $[2, 0, 0]$ and the returned value is $[5.999\ldots, 0, 0]$, an absolute error of size $\sim 4$.

Central-over-central has the same structure but with better-conditioned ingredients:

• Single central FD has truncation $O(h^2)$ and roundoff $O(\varepsilon/h)$, optimal $h \approx \varepsilon^{1/3}$, error $\approx \varepsilon^{2/3} \approx 4 \cdot 10^{-11}$.
• Central FD of a central FD: numerator error $\sim \varepsilon^{2/3}$, divided by $h \sim \varepsilon^{1/3}$, gives final error $\sim \varepsilon^{1/3} \approx 6 \cdot 10^{-6}$.

That also matches the observation: AutoFiniteDiff(; fdtype = Val(:central)) returns [1.99999930..., 0, ~4e-6] for the same call, with the residual at the predicted $\varepsilon^{1/3}$ scale.

So the workaround mentioned above is not just empirically better — it shifts the asymptotic error from $O(1)$ to $O(\varepsilon^{1/3})$, which is small enough for the result to look "approximately correct" on the test functions in the report. References for the derivations: Trefethen & Bau, Numerical Linear Algebra, Lecture 21; Nocedal & Wright, Numerical Optimization, §8.1.

This is a standard result in finite-difference numerical analysis, but it is not currently surfaced in the DI docs or in FiniteDiff's own README, which is what makes the default AutoFiniteDiff() behavior surprising in practice.

---

🤖 I am a robot. This is an experiment in agentic bug-catching under the supervision of @adrhill and @gdalle (#1008). Contents may be hallucinated.

[bdrhill]

For comparison, AutoFiniteDifferences (FiniteDifferences.jl) does not show the same hvp/hessian split. Both operators route through the single fdm passed to the constructor, so they always use stencils of equal accuracy.

Tested with the same four functions (sum(x.^2), sum(x.^3), sum(x.^4), x' * [1 2; 3 4] * x) at x = [1, 2, 3] (or [1, 2]) and v = [1, 0, 0] (or [1, 0]):

fdm	hvp matches analytical?	hvp vs hessian * v
central_fdm(5, 1)	yes	identical
central_fdm(3, 1)	yes	identical
central_fdm(2, 1)	no	identical
forward_fdm(5, 1)	yes	identical
forward_fdm(3, 1)	yes	identical
forward_fdm(2, 1)	no	identical
backward_fdm(2, 1)	no	identical

The 2-point stencils inherit the same O(1) (forward/backward) and O(ε^{1/3}) (central) accuracy penalty discussed in the previous comment, but hvp and hessian * v agree in every case, because both paths use the user-supplied fdm. AutoFiniteDifferences also has no default — the constructor requires an explicit fdm — so the low-accuracy cases only appear when the caller has already opted into a low-order stencil.

So the inconsistency in this issue appears to be specific to AutoFiniteDiff's separation of fdtype and fdhtype. Routing hvp through fdhtype (or sharing a single field, as AutoFiniteDifferences does) would bring the two backends into alignment.

---

🤖 I am a robot. This is an experiment in agentic bug-catching under the supervision of @adrhill and @gdalle (#1008). Contents may be hallucinated.

[gdalle]

It's funny, I had actually encountered this issue in DI testing but never thought to fix it. Routing hvp through fdhtype inside DI is a good fix. Wanna give it a go @bdrhill? Open a PR