Bug(FastDifferentiation): hessian returns wrong values when a variable appears inside a nonlinear op and as a multiplier

[bdrhill]

Description

hessian with AutoFastDifferentiation() returns incorrect values when a variable appears both inside a nonlinear function (e.g. exp, sin, log, ^) and as a separate multiplier outside. The same incorrect result is produced by native FastDifferentiation.hessian and FastDifferentiation.sparse_hessian, so the root cause is in FastDifferentiation.jl. The chained FastDifferentiation.derivative(FastDifferentiation.derivative([f], x), x) path produces the correct result, suggesting the bug is in the multivariate Hessian assembly rather than the underlying derivative rules.

MWE

using DifferentiationInterface
using FastDifferentiation: FastDifferentiation
using ForwardDiff: ForwardDiff

f(v) = exp(v[1] - v[2]) * v[2]
x = [-0.7, 2.3]

# Analytical: ∂²f/∂x[2]² = exp(x[1]-x[2])*(x[2]-2) = exp(-3)*0.3 ≈ 0.014936

hessian(f, AutoForwardDiff(), x)
# 2×2 Matrix{Float64}:
#   0.11451    -0.0647232
#  -0.0647232   0.0149361      ← correct (matches analytical)

hessian(f, AutoFastDifferentiation(), x)
# 2×2 Matrix{Float64}:
#   0.11451    -0.0647232
#  -0.0647232   0.129446       ← H[2,2] is off by ~9x, equal to 2*(x[2]-1)*exp(x[1]-x[2])

H[1,1] and H[1,2]=H[2,1] are correct; only H[2,2] is wrong, off by exp(x[1]-x[2]) * x[2] = f(x) itself.

Range of affected functions

The bug appears whenever a variable y is both inside a nonlinear operator and multiplied outside it; another variable in the nonlinear operator is necessary. The pattern is not specific to exp.

Function	H[2,2] diff (FD vs FastDiff)
exp(x1-x2)*x2	0.115
x2*exp(x1-x2)	0.115
exp(x1+x2)*x2	11.39
sin(x1-x2)*x2	0.325
log(x1+x2)*x2 (x>0)	0.082
(x1-x2)^2*x2	4.6
exp(-x2)*x2 (only one var inside)	0.0 (correct)
exp(x1)*x2 (no shared var)	0.0 (correct)
(x1-x2)*x2 (no nonlinearity)	0.0 (correct)

Native Backend Comparison

The same wrong value is produced by FastDifferentiation.hessian, FastDifferentiation.sparse_hessian, and FastDifferentiation.jacobian(FastDifferentiation.jacobian([f], xv), xv). Manual FastDifferentiation.derivative(FastDifferentiation.derivative([f], xv[2]), xv[2]) gives the correct value, so the bug is in the higher-level Hessian/Jacobian assembly rather than the derivative rules themselves.

using FastDifferentiation
xv = make_variables(:x, 2)
fexpr = exp(xv[1] - xv[2]) * xv[2]

# FastDifferentiation.hessian: WRONG H[2,2]
H_sym = FastDifferentiation.hessian(fexpr, xv)
make_function(H_sym, xv; in_place=false)([-0.7, 2.3])
# 2×2 Matrix{Float64}:
#   0.11451    -0.0647232
#  -0.0647232   0.129446        ← wrong

# Manual chained derivative: CORRECT
d1 = FastDifferentiation.derivative([fexpr], xv[2])
d2 = FastDifferentiation.derivative(d1, xv[2])
make_function(d2, xv; in_place=false)([-0.7, 2.3])
# 1-element Vector{Float64}:
#  0.014936120510359176          ← correct

The symbolic Hessian FD produces for the (2,2) entry is:

(-((exp((x1 - x2)) + (exp((x1 - x2)) * -(x2)))) + ((-(x2) + 1) * -(exp((x1 - x2)))))

This simplifies algebraically to 2*(x[2]-1)*exp(x[1]-x[2]). The correct answer is (x[2]-2)*exp(x[1]-x[2]). The second term in FD's expression, ((-(x[2])+1) * -(exp((x[1] - x[2])))), carries a spurious factor of (1-x[2]); without it, the formula matches the correct manual result.

Expected Behavior

hessian(f, AutoFastDifferentiation(), x) matches hessian(f, AutoForwardDiff(), x) for these functions, to within numerical precision.

Actual Behavior

H[2,2] disagrees, with relative error of order f(x)/H_true[2,2] (~9x in the exp(x1-x2)*x2 example).

Backend

• Backend: AutoFastDifferentiation()
• Works with other backends: Yes — AutoForwardDiff(), AutoHyperHessians(), AutoSymbolics() all match the analytical answer
• Native API gives same result: No — FastDifferentiation.hessian produces the same incorrect value

This may be related to FastDifferentiation issues #65 (factorization algorithm) or #76 (closed; simplification rule for matching terms).

Environment

• Julia 1.12.6
• DifferentiationInterface v0.7.18
• FastDifferentiation v0.4.5
• ForwardDiff v1.3.3

Full environment

Status `/tmp/di-bughunt/Project.toml`
  [47edcb42] ADTypes v1.22.0
  [a0c0ee7d] DifferentiationInterface v0.7.18
  [eb9bf01b] FastDifferentiation v0.4.5
  [f6369f11] ForwardDiff v1.3.3

Julia Version 1.12.6
Commit 15346901f00 (2026-04-09 19:20 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × Intel(R) Core(TM) i5-2520M CPU @ 2.50GHz
  WORD_SIZE: 64

---

🤖 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]

@brianguenter maybe this is of interest to you

[brianguenter]

@gdalle I'll take a look. How hard was it to set up the AI bug-catching? That's something I would like to learn. Got any pointers?

[adrhill]

Hi Brian, there currently is not to much to it: main...bdrhill:DifferentiationInterface.jl:bh/find-bugs

The /find-bugs/SKILL.md is a reusable prompt which keeps a LOG.md as a "memory" of things that have already been tested for in previous runs. I know some people have much more sophisticated and automated setups (e.g., @ChrisRackauckas-Claude), but I just call this skill manually, as I don't want it to accidentally run through my usage limits. I run Claude with --dangerously-skip-permissions on a separate headless computer that I set up as a sandbox and ssh into. The @bdrhill account functions as a sandbox as well.