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

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.

Start with test/Back/FiniteDiff/test.jl and the HVP path described as the generic JVP-of-gradient fallback; compare how fdtype and fdhtype are used for AutoFiniteDiff. Check backends.md for the documented support and limitations. Done means adding regression coverage for the reported HVP cases and either aligning accuracy with hessian or documenting the chosen limitation and workaround.