Exploit symmetry in the Hessian

[devmotion]

Fixes #253. Fixes #836. Fixes #845. Fixes #846.

Relation to #837, #840 and #843

• exploit symmetry in the hessian #837 (KristofferC) is the original. This PR carries its two commits unmodified and answers its second review round.
• Exploit symmetry in the Hessian, on top of #840 #843 is that same work rebuilt on Extract derivatives at the indices of x, and store those indices in the config #840, whose structured-array design is still under discussion.
• This PR is Exploit symmetry in the Hessian, on top of #840 #843 without Extract derivatives at the indices of x, and store those indices in the config #840. Same sweep, same fixes. gradient, jacobian, GradientConfig and JacobianConfig are untouched, and gradient! writes to the wrong entries for structured inputs #838, Chunked jacobian throws for Diagonal/LowerTriangular/UpperTriangular inputs #839 and Reusing a GradientConfig/JacobianConfig with a differently structured input silently computes wrong derivatives #842 are out of scope — reusing a config across structures stays exactly as wrong as on master.

kc/symmetric_hessian is untouched.

Numbers

Suite 9759/9759. hessian! into a preallocated matrix, best of 200, Julia 1.12.7:

n	rosenbrock	ackley	allocations (rosenbrock)
10 (single chunk)	1.11x	1.15x	0 vs 1040 B
30	1.59x	1.58x	0 vs 9456 B
100	1.80x	1.84x	0 vs 111216 B

Over 4 functions × n ∈ {1,2,3,5,8,13} × chunk sizes × 4 result kinds, 372 of 432 results are bitwise identical to master. The 60 that differ are exactly the 60 master returned asymmetric, by at most one ulp.

#840's seeding speedup, confined to the Hessian. seed_hessian_chunk! takes its window from stored positions rather than walking Iterators.drop, so its cost stops growing with n. Seeding the final chunk:

x	before	after
Vector n=40 / n=200	25.5 / 93.8 ns	11.1 / 11.0 ns
UpperTriangular 13² / 25²	99.4 / 337.7 ns	13.8 / 13.9 ns
Diagonal 40 / 200	32.1 / 105.9 ns	18.0 / 18.1 ns

End-to-end that is 1.15x at n = 40 and 1.21x at n = 80 with chunk = 2, and nothing at chunk = 12. #840's own gradient!/jacobian! gains are not included — those paths are untouched here.

Fixes

• Exact symmetry (Hessian should be symmetric #253): both triangles are filled from the same value. Over that sweep master is exactly symmetric in 228/288 results, this branch in 288/288.
• Chunk-size reproducibility: diagonal blocks read one triangle, off-diagonal blocks the other. At n = 16 with log(sum(exp, z)), chunks 1, 2, 3, 5, 7, 11 each differed from chunk 16. Every chunk size is now bitwise equal to the StaticArrays path.
• hessian! stopped checking the result shape when the sweep replaced extract_jacobian!: hessian!(fill(NaN, 4, 4), f, rand(3)) silently left row/column 4 untouched. reshape_hessian restores it.
• hessian!(::DiffResult, …) stopped checking the gradient buffer, which master validated through reshape_jacobian's row count. Worst case was silent and wrong rather than incomplete: with a 4×3 buffer for an UpperTriangular(3, 3) input the six derivatives went to linear positions 1, 4, 5, 7, 8, 9 of a four-row array. structural_eachindex(grad, x) — the check the seeding utilities already apply to a work buffer — restores it, and still accepts a flat gradient buffer for a matrix x as master did.
• The StaticArrays path went through jacobian(gradient(f), x) and was not exactly symmetric either. hessian(prod, ::SMatrix{3,3}): 195 ns → 61 ns.
• The two nested dual layers shared a tag (hessian: the two nested dual layers share a tag, so results depend on the code path and chunk size #845). HessianConfig handed one Tag(f, V) to both, so value/partials inside f could not tell them apart and the ordering machinery had nothing to compare. f(z) = sum(abs2, z) + ForwardDiff.value(z[1]) * z[2] at [1.0, 2.0, 3.0] gave master three different Hessians and two different gradients depending on the path and the chunk size — and the chunk-size reproducibility above does not hold for such an f under one tag. The outer tag is now derived from the inner one as outer_tag(T, Dual{T,V,N}), registered after it so that T ≺ TO. A result carrying only the inner layer then reads as a gradient with a vanishing Hessian, with no special casing.
• extract_hessian dispatched on partials(T, ydual) (hessian!(::ImmutableDiffResult, f, ::StaticArray) errors when f does not depend on its argument #846), which is a Partials only when the result carries T. For a result carrying only an enclosing tag it is a bare Dual and matched no method, so ForwardDiff.derivative(a -> ForwardDiff.hessian(z -> a * 2.0, SVector(1.0, 2.0, 3.0))[1,1], 1.0) errored. Dispatch is now on the result itself, which also subsumes the Partials{0} shortcut for a constant f and an empty x. Zeros written for a derivative that vanishes take their element type from ydual, as extract_gradient! does, so the array and StaticArrays paths treat a plain result buffer alike.
• Blocks whose derivatives the result cannot carry are skipped. For an f whose result has no outer perturbation the off-diagonal blocks are known to be zero, so evaluations drop from nblocks(nblocks+1)/2 to nblocks, and to one when no gradient is requested.

Breaking

Structured inputs return length(x) × length(x), both axes the linear indices of x, zeros in the rows and columns of the structural zeros. hessian(f, UpperTriangular(rand(3,3))) is 9×9 with rows/columns 2, 3, 6 zero.

Master returns 9×6 — linear rows, structural columns — which cannot be indexed without reading structural_eachindex; chunk mode throws DimensionMismatch; and hessian!(DiffResults.HessianResult(x), f, x) fails to broadcast. The new convention is the smaller change (master already had linear rows) and makes HessianResult(x) the right buffer. The cost is Diagonal, now quartic in the diagonal's size: hessian! for Diagonal(40) spends its time on a 1600×1600 result rather than a 40×40 one.

HessianConfig holds iseeds, oseeds, duals instead of a JacobianConfig and a GradientConfig, and carries the outer tag, so HessianConfig{T,V,N,DG,DJ} becomes HessianConfig{T,TO,V,N,D} and eltype becomes Dual{TO,Dual{T,V,N},N} — any code that writes the type out, not just code that reads the fields, has to change. T stays first and stays the inner tag, so cfg::HessianConfig{T} signatures keep working. The sweep never reads the Jacobian config's buffer. HessianConfig(f, x): 1.422 → 1.312 MB at length(x) == 1000; the result-aware constructor 1.531 → 1.312 MB. Separate commit, easy to drop.

Both of these need a 2.0.

Results change for an f that reaches into the layers of its argument, which is the point of the #845 fix rather than a cost of it. ForwardDiff.value(z[1]) now yields the inner-layer value instead of one the outer layer re-absorbs, so the example above gives [4.0, 5.0, 6.0] — what master's array path returned — from every path and every chunk size, where master's StaticArrays path gave [2.0, 5.0, 6.0]. An f that does not reach into layers is bitwise unchanged: over DiffTests.VECTOR_TO_NUMBER_FUNCS × chunk sizes 1, 2, 3, 5, all 60 results are identical to the one-tag computation.

Implementation notes

structural_linearindices gives structural_eachindex's positions as linear indices, computed once per sweep rather than cached on the config — OneTo/StepRange for dense and Diagonal (allocation-free), one small vector per call for the triangles.

Seeding a structured buffer now walks it by linear index, which is what lets a position double as a Hessian row. Base._unsetindex! exists for Array alone — for a linear index its AbstractArray fallback recurses forever, where a CartesianIndex merely had no method — so the unassigned-entry path now raises an ArgumentError naming the entry. That also removes a pre-existing StackOverflowError for a Diagonal of a non-bits element type.

Beyond the suite: 3 wrappers × n ∈ {1,2,3,5,8} × every chunk size × 4 result kinds against a closed-form reference — 1269 checks, all passing.

🤖 Generated with Claude Code

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 92.16%. Comparing base (b742809) to head (b42dca2).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #844      +/-   ##
==========================================
+ Coverage   90.68%   92.16%   +1.48%     
==========================================
  Files          11       11              
  Lines        1052     1200     +148     
==========================================
+ Hits          954     1106     +152     
+ Misses         98       94       -4     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.