From 0091edd55c4ea8fd0acd93cdfc883e357e963058 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Sun, 16 Aug 2026 10:05:43 +0200 Subject: [PATCH 1/6] exploit symmetry in the hessian instead of relying on the jacobian of gradient for the hessian explicitly seed dual numbers and only calculate the upper triangular part when chunking, gives ~2x speedup as input length becomes big --- src/hessian.jl | 151 ++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 125 insertions(+), 26 deletions(-) diff --git a/src/hessian.jl b/src/hessian.jl index 9c755c9a..59c99e94 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -5,7 +5,7 @@ """ ForwardDiff.hessian(f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, x), check=Val{true}()) -Return `H(f)` (i.e. `J(∇(f))`) evaluated at `x`, assuming `f` is called as `f(x)`. +Return `H(f)` evaluated at `x`, assuming `f` is called as `f(x)`. This method assumes that `isa(f(x), Real)`. @@ -14,8 +14,8 @@ Set `check` to `Val{false}()` to disable tag checking. This can lead to perturba function hessian(f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {F, T,CHK} require_one_based_indexing(x) CHK && checktag(T, f, x) - ∇f = y -> gradient(f, y, cfg.gradient_config, Val{false}()) - return jacobian(∇f, x, cfg.jacobian_config, Val{false}()) + H, _ = symmetric_hessian(f, x, cfg, nothing) + return H end """ @@ -31,29 +31,12 @@ Set `check` to `Val{false}()` to disable tag checking. This can lead to perturba function hessian!(result::AbstractArray, f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {F,T,CHK} require_one_based_indexing(result, x) CHK && checktag(T, f, x) - ∇f = y -> gradient(f, y, cfg.gradient_config, Val{false}()) - jacobian!(result, ∇f, x, cfg.jacobian_config, Val{false}()) + xlen = structural_length(x) + H = result isa AbstractMatrix && size(result) == (xlen, xlen) ? result : reshape(result, xlen, xlen) + symmetric_hessian!(H, f, x, cfg, nothing) return result end - -# We use this struct below instead of an -# equivalent closure in order to avoid -# JuliaLang/julia#15276-related performance -# issues. See #316. -mutable struct InnerGradientForHess{R,C,F} - result::R - cfg::C - f::F -end - -function (g::InnerGradientForHess)(y, z) - inner_result = DiffResult(zero(eltype(y)), y) - gradient!(inner_result, g.f, z, g.cfg.gradient_config, Val{false}()) - g.result = DiffResults.value!(g.result, value(DiffResults.value(inner_result))) - return y -end - """ ForwardDiff.hessian!(result::DiffResult, f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, result, x), check=Val{true}()) @@ -64,8 +47,124 @@ because `isa(result, DiffResult)`, `cfg` is constructed as `HessianConfig(f, res Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ function hessian!(result::DiffResult, f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, result, x), ::Val{CHK}=Val{true}()) where {F,T,CHK} + require_one_based_indexing(x) CHK && checktag(T, f, x) - ∇f! = InnerGradientForHess(result, cfg, f) - jacobian!(DiffResults.hessian(result), ∇f!, DiffResults.gradient(result), x, cfg.jacobian_config, Val{false}()) - return ∇f!.result + xlen = structural_length(x) + hess = DiffResults.hessian(result) + H = hess isa AbstractMatrix && size(hess) == (xlen, xlen) ? hess : reshape(hess, xlen, xlen) + _, ydual = symmetric_hessian!(H, f, x, cfg, DiffResults.gradient(result)) + result = DiffResults.value!(result, value(T, value(T, ydual))) + return result +end + +############################ +# symmetric Hessian kernel # +############################ + +const HESSIAN_ERROR = DimensionMismatch("hessian(f, x) expects that f(x) is a real number. Perhaps you meant jacobian(f, x)?") + +# Seed a chunk in either layer of the nested duals. A `nothing` seed clears that layer. +function seed_hessian_chunk!(duals::AbstractArray{Dual{T,Dual{T,V,N},N}}, x, index, + iseeds::Union{Nothing,NTuple{N,Partials{N,V}}}, + oseeds::Union{Nothing,NTuple{N,Partials{N,Dual{T,V,N}}}}, + chunksize = N) where {T,V,N} + izero = zero(Partials{N,V}) + ozero = zero(Partials{N,Dual{T,V,N}}) + idxs = Iterators.drop(structural_eachindex(duals, x), index - 1) + if isbitstype(V) + for (i, idx) in zip(1:chunksize, idxs) + inner = Dual{T,V,N}(x[idx], iseeds === nothing ? izero : iseeds[i]) + duals[idx] = Dual{T,Dual{T,V,N},N}(inner, oseeds === nothing ? ozero : oseeds[i]) + end + else + for (i, idx) in zip(1:chunksize, idxs) + if isassigned(x, idx) + inner = Dual{T,V,N}(x[idx], iseeds === nothing ? izero : iseeds[i]) + duals[idx] = Dual{T,Dual{T,V,N},N}(inner, oseeds === nothing ? ozero : oseeds[i]) + else + Base._unsetindex!(duals, idx) + end + end + end + return duals +end + +# Copy a block from the nested partials and fill its transpose. On diagonal blocks, read +# only the upper triangle so the result is exactly symmetric. +function extract_hessian_chunk!(::Type{T}, H, ydual, roffset, coffset, rsize, csize) where {T} + for r in 1:rsize + drow = partials(T, ydual, r) + cstart = roffset == coffset ? r : 1 + for c in cstart:csize + h = partials(T, drow, c) + H[roffset + r, coffset + c] = h + H[coffset + c, roffset + r] = h + end + end + return H +end + +# The inner partials of a diagonal block contain the corresponding gradient chunk. +extract_hessian_gradient_chunk!(::Type{T}, ::Nothing, ydual, index, chunksize) where {T} = nothing +extract_hessian_gradient_chunk!(::Type{T}, grad, ydual, index, chunksize) where {T} = + extract_gradient_chunk!(T, grad, value(T, ydual), index, chunksize) + +# Evaluate one pair of chunks at a time using nested duals. Only one triangle of block +# pairs is evaluated; the other is filled by symmetry (see #836). +function symmetric_hessian_expr(result_definition::Expr) + return quote + xlen = structural_length(x) + if xlen < N + throw(ArgumentError(lazy"chunk size cannot be greater than ForwardDiff.structural_length(x) ($(N) > $(structural_length(x)))")) + end + + nblocks = xlen == 0 ? 1 : div(xlen + N - 1, N) + + xdual = cfg.gradient_config.duals + iseeds = cfg.jacobian_config.seeds + oseeds = cfg.gradient_config.seeds + + # Keep all unseeded blocks at zero between evaluations. + seed_hessian_chunk!(xdual, x, 1, nothing, nothing, xlen) + + # The first evaluation determines the output type. + seed_hessian_chunk!(xdual, x, 1, iseeds, oseeds) + ydual1 = f(xdual) + ydual1 isa Real || throw(HESSIAN_ERROR) + $(result_definition) + extract_hessian_chunk!(T, H, ydual1, 0, 0, N, N) + extract_hessian_gradient_chunk!(T, grad, ydual1, 1, N) + seed_hessian_chunk!(xdual, x, 1, nothing, nothing) + + for q in 2:nblocks + qoffset = (q - 1) * N + qsize = min(N, xlen - qoffset) + # Off-diagonal blocks: p seeds columns and q seeds rows. + for p in 1:(q - 1) + poffset = (p - 1) * N + seed_hessian_chunk!(xdual, x, poffset + 1, iseeds, nothing) + seed_hessian_chunk!(xdual, x, qoffset + 1, nothing, oseeds, qsize) + ydual = f(xdual) + extract_hessian_chunk!(T, H, ydual, qoffset, poffset, qsize, N) + seed_hessian_chunk!(xdual, x, poffset + 1, nothing, nothing) + seed_hessian_chunk!(xdual, x, qoffset + 1, nothing, nothing, qsize) + end + # Diagonal blocks seed both layers. + seed_hessian_chunk!(xdual, x, qoffset + 1, iseeds, oseeds, qsize) + ydual = f(xdual) + extract_hessian_chunk!(T, H, ydual, qoffset, qoffset, qsize, qsize) + extract_hessian_gradient_chunk!(T, grad, ydual, qoffset + 1, qsize) + seed_hessian_chunk!(xdual, x, qoffset + 1, nothing, nothing, qsize) + end + + return H, ydual1 + end +end + +@eval function symmetric_hessian(f::F, x, cfg::HessianConfig{T,V,N}, grad) where {F,T,V,N} + $(symmetric_hessian_expr(:(H = similar(x, typeof(value(T, value(T, ydual1))), xlen, xlen)))) +end + +@eval function symmetric_hessian!(H, f::F, x, cfg::HessianConfig{T,V,N}, grad) where {F,T,V,N} + $(symmetric_hessian_expr(:())) end From 0c04e90a91aa52e02c3ee90dc7ee4e2ec52db5d3 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Mon, 17 Aug 2026 11:57:17 +0200 Subject: [PATCH 2/6] Address symmetric Hessian review feedback --- ext/ForwardDiffStaticArraysExt.jl | 32 +++++++++++++-- src/apiutils.jl | 63 +++++++++++++++-------------- src/config.jl | 11 ++--- src/hessian.jl | 67 +++++++++++-------------------- test/AllocationsTest.jl | 10 +++++ test/GradientTest.jl | 5 +++ test/HessianTest.jl | 66 ++++++++++++++++++++++++++++++ test/JacobianTest.jl | 4 ++ test/SeedTest.jl | 18 +++++++++ 9 files changed, 192 insertions(+), 84 deletions(-) diff --git a/ext/ForwardDiffStaticArraysExt.jl b/ext/ForwardDiffStaticArraysExt.jl index bf0ef99a..26abf43b 100644 --- a/ext/ForwardDiffStaticArraysExt.jl +++ b/ext/ForwardDiffStaticArraysExt.jl @@ -7,7 +7,7 @@ using ForwardDiff: Dual, partials, npartials, Partials, GradientConfig, Jacobian gradient, hessian, jacobian, gradient!, hessian!, jacobian!, extract_gradient!, extract_jacobian!, extract_value!, vector_mode_gradient, vector_mode_gradient!, - vector_mode_jacobian, vector_mode_jacobian!, valtype, value + vector_mode_jacobian, vector_mode_jacobian!, HESSIAN_ERROR, valtype, value using DiffResults: DiffResult, ImmutableDiffResult, MutableDiffResult @generated function dualize(::Type{T}, x::StaticArray) where T @@ -107,11 +107,34 @@ end end # Hessian -ForwardDiff.hessian(f::F, x::StaticArray) where {F} = jacobian(Base.Fix1(gradient, f), x) +@inline function extract_hessian(::Type{T}, ydual::Partials, x::StaticArray) where {T} + H = extract_jacobian(T, ydual, x) + return typeof(H)(Symmetric(H, :U)) +end + +@inline function extract_hessian(::Type{T}, ydual::Partials{0}, x::S) where {T,S<:StaticArray} + R = StaticArrays.similar_type(S, valtype(T, eltype(ydual)), Size(length(x), length(x))) + return zero(R) +end + +@inline function ForwardDiff.hessian(f::F, x::StaticArray) where {F} + T = typeof(Tag(f, eltype(x))) + ydual = f(dualize(T, dualize(T, x))) + ydual isa Real || throw(HESSIAN_ERROR) + return extract_hessian(T, partials(T, ydual), x) +end + ForwardDiff.hessian(f::F, x::StaticArray, cfg::HessianConfig) where {F} = hessian(f, x) ForwardDiff.hessian(f::F, x::StaticArray, cfg::HessianConfig, ::Val) where {F} = hessian(f, x) -ForwardDiff.hessian!(result::AbstractArray, f::F, x::StaticArray) where {F} = jacobian!(result, Base.Fix1(gradient, f), x) +@inline function ForwardDiff.hessian!(result::AbstractArray, f::F, x::StaticArray) where {F} + T = typeof(Tag(f, eltype(x))) + ydual = f(dualize(T, dualize(T, x))) + ydual isa Real || throw(HESSIAN_ERROR) + H = result isa AbstractMatrix ? result : reshape(result, length(x), length(x)) + ForwardDiff.extract_hessian_chunk!(T, H, ydual, 0, 0, length(x), length(x)) + return result +end ForwardDiff.hessian!(result::MutableDiffResult, f::F, x::StaticArray) where {F} = hessian!(result, f, x, HessianConfig(f, result, x)) @@ -123,9 +146,10 @@ function ForwardDiff.hessian!(result::ImmutableDiffResult, f::F, x::StaticArray) d1 = dualize(T, x) d2 = dualize(T, d1) fd2 = f(d2) + fd2 isa Real || throw(HESSIAN_ERROR) val = value(T,value(T,fd2)) grad = extract_gradient(T,value(T,fd2), x) - hess = extract_jacobian(T,partials(T,fd2), x) + hess = extract_hessian(T,partials(T,fd2), x) result = DiffResults.hessian!(result, hess) result = DiffResults.gradient!(result, grad) result = DiffResults.value!(result, val) diff --git a/src/apiutils.jl b/src/apiutils.jl index 0615fdb3..1d54d7bb 100644 --- a/src/apiutils.jl +++ b/src/apiutils.jl @@ -88,14 +88,22 @@ end function _seed_zero_partials!(duals::AbstractArray{Dual{T,V,N}}, x, idxs) where {T,V,N} seed = zero(Partials{N,V}) + return _seed!(duals, x, idxs) do value, _ + Dual{T,V,N}(value, seed) + end +end + +# Write a sequence of duals while preserving unassigned entries in arrays whose element type is not +# stored inline. `make_dual` receives the primal value and its one-based position in `idxs`. +@inline function _seed!(make_dual::F, duals::AbstractArray{Dual{T,V,N}}, x, idxs) where {F,T,V,N} if isbitstype(V) - for idx in idxs - duals[idx] = Dual{T,V,N}(x[idx], seed) + for (i, idx) in enumerate(idxs) + duals[idx] = make_dual(x[idx], i) end else - for idx in idxs + for (i, idx) in enumerate(idxs) if isassigned(x, idx) - duals[idx] = Dual{T,V,N}(x[idx], seed) + duals[idx] = make_dual(x[idx], i) else Base._unsetindex!(duals, idx) end @@ -106,38 +114,31 @@ end function seed!(duals::AbstractArray{Dual{T,V,N}}, x, seeds::NTuple{N,Partials{N,V}}) where {T,V,N} - if isbitstype(V) - for (i, idx) in zip(1:N, structural_eachindex(duals, x)) - duals[idx] = Dual{T,V,N}(x[idx], seeds[i]) - end - else - for (i, idx) in zip(1:N, structural_eachindex(duals, x)) - if isassigned(x, idx) - duals[idx] = Dual{T,V,N}(x[idx], seeds[i]) - else - Base._unsetindex!(duals, idx) - end - end + idxs = Iterators.take(structural_eachindex(duals, x), N) + return _seed!(duals, x, idxs) do value, i + Dual{T,V,N}(value, seeds[i]) end - return duals end function seed!(duals::AbstractArray{Dual{T,V,N}}, x, index, seeds::NTuple{N,Partials{N,V}}, chunksize = N) where {T,V,N} offset = index - 1 - idxs = Iterators.drop(structural_eachindex(duals, x), offset) - if isbitstype(V) - for (i, idx) in zip(1:chunksize, idxs) - duals[idx] = Dual{T,V,N}(x[idx], seeds[i]) - end - else - for (i, idx) in zip(1:chunksize, idxs) - if isassigned(x, idx) - duals[idx] = Dual{T,V,N}(x[idx], seeds[i]) - else - Base._unsetindex!(duals, idx) - end - end + idxs = Iterators.take(Iterators.drop(structural_eachindex(duals, x), offset), chunksize) + return _seed!(duals, x, idxs) do value, i + Dual{T,V,N}(value, seeds[i]) + end +end + +# Seed a chunk in either layer of nested duals. A `nothing` seed clears that layer. +function seed_hessian_chunk!(duals::AbstractArray{Dual{T,Dual{T,V,N},N}}, x, index, + iseeds::Union{Nothing,NTuple{N,Partials{N,V}}}, + oseeds::Union{Nothing,NTuple{N,Partials{N,Dual{T,V,N}}}}, + chunksize = N) where {T,V,N} + izero = zero(Partials{N,V}) + ozero = zero(Partials{N,Dual{T,V,N}}) + idxs = Iterators.take(Iterators.drop(structural_eachindex(duals, x), index - 1), chunksize) + return _seed!(duals, x, idxs) do value, i + inner = Dual{T,V,N}(value, iseeds === nothing ? izero : iseeds[i]) + Dual{T,Dual{T,V,N},N}(inner, oseeds === nothing ? ozero : oseeds[i]) end - return duals end diff --git a/src/config.jl b/src/config.jl index 3c6c97e3..58c145f3 100644 --- a/src/config.jl +++ b/src/config.jl @@ -207,10 +207,9 @@ Return a `HessianConfig` instance based on the type of `f` and type/shape of the vector `x`. The returned `HessianConfig` instance contains all the work buffers required by -`ForwardDiff.hessian` and `ForwardDiff.hessian!`. For the latter, the buffers are -configured for the case where the `result` argument is an `AbstractArray`. If -it is a `DiffResult`, the `HessianConfig` should instead be constructed via -`ForwardDiff.HessianConfig(f, result, x, chunk)`. +`ForwardDiff.hessian` and `ForwardDiff.hessian!`, including when the latter stores into a +`DiffResult`. The `ForwardDiff.HessianConfig(f, result, x, chunk)` constructor may also be +used with any of these methods. If `f` is `nothing` instead of the actual target function, then the returned instance can be used with any target function. However, this will reduce ForwardDiff's ability to catch @@ -234,7 +233,9 @@ Return a `HessianConfig` instance based on the type of `f`, types/storage in `re type/shape of the input vector `x`. The returned `HessianConfig` instance contains all the work buffers required by -`ForwardDiff.hessian!` for the case where the `result` argument is an `DiffResult`. +`ForwardDiff.hessian` and `ForwardDiff.hessian!`. It is interchangeable with a config +constructed via `ForwardDiff.HessianConfig(f, x, chunk)`; this constructor retains the +result-aware form for compatibility. If `f` is `nothing` instead of the actual target function, then the returned instance can be used with any target function. However, this will reduce ForwardDiff's ability to catch diff --git a/src/hessian.jl b/src/hessian.jl index 59c99e94..489b9014 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -6,6 +6,8 @@ ForwardDiff.hessian(f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, x), check=Val{true}()) Return `H(f)` evaluated at `x`, assuming `f` is called as `f(x)`. +The returned Hessian is exactly symmetric: its two triangles are filled from the same +derivative values. This method assumes that `isa(f(x), Real)`. @@ -21,8 +23,9 @@ end """ ForwardDiff.hessian!(result::AbstractArray, f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, x), check=Val{true}()) -Compute `H(f)` (i.e. `J(∇(f))`) evaluated at `x` and store the result(s) in `result`, -assuming `f` is called as `f(x)`. +Compute `H(f)` evaluated at `x` and store the result(s) in `result`, assuming `f` is +called as `f(x)`. The stored Hessian is exactly symmetric: its two triangles are filled +from the same derivative values. This method assumes that `isa(f(x), Real)`. @@ -32,7 +35,7 @@ function hessian!(result::AbstractArray, f::F, x::AbstractArray, cfg::HessianCon require_one_based_indexing(result, x) CHK && checktag(T, f, x) xlen = structural_length(x) - H = result isa AbstractMatrix && size(result) == (xlen, xlen) ? result : reshape(result, xlen, xlen) + H = result isa AbstractMatrix ? result : reshape(result, xlen, xlen) symmetric_hessian!(H, f, x, cfg, nothing) return result end @@ -40,9 +43,10 @@ end """ ForwardDiff.hessian!(result::DiffResult, f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, result, x), check=Val{true}()) -Exactly like `ForwardDiff.hessian!(result::AbstractArray, f, x::AbstractArray, cfg::HessianConfig)`, but -because `isa(result, DiffResult)`, `cfg` is constructed as `HessianConfig(f, result, x)` instead of -`HessianConfig(f, x)`. +Exactly like `ForwardDiff.hessian!(result::AbstractArray, f, x::AbstractArray, cfg::HessianConfig)`, +but also stores the value and gradient in `result`. The default `cfg` is constructed as +`HessianConfig(f, result, x)`, though a config constructed as `HessianConfig(f, x)` may also +be used. Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ @@ -51,7 +55,7 @@ function hessian!(result::DiffResult, f::F, x::AbstractArray, cfg::HessianConfig CHK && checktag(T, f, x) xlen = structural_length(x) hess = DiffResults.hessian(result) - H = hess isa AbstractMatrix && size(hess) == (xlen, xlen) ? hess : reshape(hess, xlen, xlen) + H = hess isa AbstractMatrix ? hess : reshape(hess, xlen, xlen) _, ydual = symmetric_hessian!(H, f, x, cfg, DiffResults.gradient(result)) result = DiffResults.value!(result, value(T, value(T, ydual))) return result @@ -63,32 +67,6 @@ end const HESSIAN_ERROR = DimensionMismatch("hessian(f, x) expects that f(x) is a real number. Perhaps you meant jacobian(f, x)?") -# Seed a chunk in either layer of the nested duals. A `nothing` seed clears that layer. -function seed_hessian_chunk!(duals::AbstractArray{Dual{T,Dual{T,V,N},N}}, x, index, - iseeds::Union{Nothing,NTuple{N,Partials{N,V}}}, - oseeds::Union{Nothing,NTuple{N,Partials{N,Dual{T,V,N}}}}, - chunksize = N) where {T,V,N} - izero = zero(Partials{N,V}) - ozero = zero(Partials{N,Dual{T,V,N}}) - idxs = Iterators.drop(structural_eachindex(duals, x), index - 1) - if isbitstype(V) - for (i, idx) in zip(1:chunksize, idxs) - inner = Dual{T,V,N}(x[idx], iseeds === nothing ? izero : iseeds[i]) - duals[idx] = Dual{T,Dual{T,V,N},N}(inner, oseeds === nothing ? ozero : oseeds[i]) - end - else - for (i, idx) in zip(1:chunksize, idxs) - if isassigned(x, idx) - inner = Dual{T,V,N}(x[idx], iseeds === nothing ? izero : iseeds[i]) - duals[idx] = Dual{T,Dual{T,V,N},N}(inner, oseeds === nothing ? ozero : oseeds[i]) - else - Base._unsetindex!(duals, idx) - end - end - end - return duals -end - # Copy a block from the nested partials and fill its transpose. On diagonal blocks, read # only the upper triangle so the result is exactly symmetric. function extract_hessian_chunk!(::Type{T}, H, ydual, roffset, coffset, rsize, csize) where {T} @@ -118,38 +96,39 @@ function symmetric_hessian_expr(result_definition::Expr) throw(ArgumentError(lazy"chunk size cannot be greater than ForwardDiff.structural_length(x) ($(N) > $(structural_length(x)))")) end - nblocks = xlen == 0 ? 1 : div(xlen + N - 1, N) + # `N == 0` only for empty inputs, which still need one evaluation to determine the + # output type and value. + nblocks = xlen == 0 ? 1 : cld(xlen, N) xdual = cfg.gradient_config.duals iseeds = cfg.jacobian_config.seeds oseeds = cfg.gradient_config.seeds - # Keep all unseeded blocks at zero between evaluations. - seed_hessian_chunk!(xdual, x, 1, nothing, nothing, xlen) - - # The first evaluation determines the output type. + # The first evaluation determines the output type. Seeding the first block and clearing + # the untouched tail partitions the fresh buffer, so every element is initialized once. seed_hessian_chunk!(xdual, x, 1, iseeds, oseeds) + seed_hessian_chunk!(xdual, x, N + 1, nothing, nothing, xlen - N) ydual1 = f(xdual) ydual1 isa Real || throw(HESSIAN_ERROR) $(result_definition) extract_hessian_chunk!(T, H, ydual1, 0, 0, N, N) extract_hessian_gradient_chunk!(T, grad, ydual1, 1, N) - seed_hessian_chunk!(xdual, x, 1, nothing, nothing) + nblocks > 1 && seed_hessian_chunk!(xdual, x, 1, nothing, nothing) for q in 2:nblocks qoffset = (q - 1) * N qsize = min(N, xlen - qoffset) - # Off-diagonal blocks: p seeds columns and q seeds rows. + # Off-diagonal blocks: p seeds columns and q seeds rows. The outer seeds for q + # remain unchanged throughout this loop. + seed_hessian_chunk!(xdual, x, qoffset + 1, nothing, oseeds, qsize) for p in 1:(q - 1) poffset = (p - 1) * N seed_hessian_chunk!(xdual, x, poffset + 1, iseeds, nothing) - seed_hessian_chunk!(xdual, x, qoffset + 1, nothing, oseeds, qsize) ydual = f(xdual) extract_hessian_chunk!(T, H, ydual, qoffset, poffset, qsize, N) seed_hessian_chunk!(xdual, x, poffset + 1, nothing, nothing) - seed_hessian_chunk!(xdual, x, qoffset + 1, nothing, nothing, qsize) end - # Diagonal blocks seed both layers. + # The diagonal block adds q's inner seeds while retaining its outer seeds. seed_hessian_chunk!(xdual, x, qoffset + 1, iseeds, oseeds, qsize) ydual = f(xdual) extract_hessian_chunk!(T, H, ydual, qoffset, qoffset, qsize, qsize) @@ -162,7 +141,7 @@ function symmetric_hessian_expr(result_definition::Expr) end @eval function symmetric_hessian(f::F, x, cfg::HessianConfig{T,V,N}, grad) where {F,T,V,N} - $(symmetric_hessian_expr(:(H = similar(x, typeof(value(T, value(T, ydual1))), xlen, xlen)))) + $(symmetric_hessian_expr(:(H = similar(x, valtype(T, valtype(T, typeof(ydual1))), xlen, xlen)))) end @eval function symmetric_hessian!(H, f::F, x, cfg::HessianConfig{T,V,N}, grad) where {F,T,V,N} diff --git a/test/AllocationsTest.jl b/test/AllocationsTest.jl index 94e7cddd..3a59a5ad 100644 --- a/test/AllocationsTest.jl +++ b/test/AllocationsTest.jl @@ -29,6 +29,16 @@ convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,F allocs_szp!(duals, x, 1, 4) @test iszero(allocs_szp!(duals, x, 1, 4)) + hcfg = ForwardDiff.HessianConfig(nothing, x) + hduals = hcfg.gradient_config.duals + iseeds = hcfg.jacobian_config.seeds + oseeds = hcfg.gradient_config.seeds + allocs_hseed!(args...) = @allocated ForwardDiff.seed_hessian_chunk!(args...) + allocs_hseed!(hduals, x, 1, iseeds, oseeds) + @test iszero(allocs_hseed!(hduals, x, 1, iseeds, oseeds)) + allocs_hseed!(hduals, x, 1, nothing, nothing, 4) + @test iszero(allocs_hseed!(hduals, x, 1, nothing, nothing, 4)) + allocs_convert_test_574() = @allocated convert_test_574() allocs_convert_test_574() @test iszero(allocs_convert_test_574()) diff --git a/test/GradientTest.jl b/test/GradientTest.jl index bf121239..c9967812 100644 --- a/test/GradientTest.jl +++ b/test/GradientTest.jl @@ -56,6 +56,7 @@ end cfgx = ForwardDiff.GradientConfig(sin, x) @test_throws ForwardDiff.InvalidTagException ForwardDiff.gradient(f, x, cfgx) @test ForwardDiff.gradient(f, x, cfgx, Val{false}()) == ForwardDiff.gradient(f,x) +@test_throws ArgumentError ForwardDiff.gradient(f, x, ForwardDiff.GradientConfig(f, x, ForwardDiff.Chunk{length(x) + 1}())) ######################## @@ -115,6 +116,10 @@ end ForwardDiff.gradient!(out, prod, sx, scfg) @test out == actual + out = similar(x) + ForwardDiff.gradient!(out, prod, sx, scfg, Val{false}()) + @test out == actual + result = DiffResults.GradientResult(x) result = ForwardDiff.gradient!(result, prod, x) diff --git a/test/HessianTest.jl b/test/HessianTest.jl index 8be72ee5..119fd14d 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -51,11 +51,21 @@ h = [-66.0 -40.0 0.0; @test isapprox(DiffResults.value(out), v) @test isapprox(DiffResults.gradient(out), g) @test isapprox(DiffResults.hessian(out), h) + + # The result-aware and result-independent config constructors are interchangeable. + out = DiffResults.HessianResult(x) + ForwardDiff.hessian!(out, f, x, cfg) + @test isapprox(DiffResults.value(out), v) + @test isapprox(DiffResults.gradient(out), g) + @test isapprox(DiffResults.hessian(out), h) end cfgx = ForwardDiff.HessianConfig(sin, x) @test_throws ForwardDiff.InvalidTagException ForwardDiff.hessian(f, x, cfgx) @test ForwardDiff.hessian(f, x, cfgx, Val{false}()) == ForwardDiff.hessian(f,x) +@test_throws ArgumentError ForwardDiff.hessian(f, x, ForwardDiff.HessianConfig(f, x, ForwardDiff.Chunk{length(x) + 1}())) +@test_throws DimensionMismatch ForwardDiff.hessian(identity, x) +@test_throws DimensionMismatch ForwardDiff.hessian!(similar(x, 3, 3), identity, x) ######################## @@ -108,10 +118,22 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) @test ForwardDiff.hessian(prod, sx, scfg, Val{false}()) == actual @test ForwardDiff.hessian(prod, sx, scfg, Val{false}()) isa StaticArray + symmetry_f(z) = sum(sin(z[i]) / (1 + z[mod1(i + 1, length(z))]^2) for i in eachindex(z)) + symmetric_static = ForwardDiff.hessian(symmetry_f, sx) + @test symmetric_static == transpose(symmetric_static) + @test symmetric_static == ForwardDiff.hessian(symmetry_f, x) + @test all(iszero, ForwardDiff.hessian(Returns(2.0), sx)) + @test_throws DimensionMismatch ForwardDiff.hessian(identity, sx) + out = similar(x, 9, 9) ForwardDiff.hessian!(out, prod, sx) @test out == actual + out = similar(x, 9, 9) + ForwardDiff.hessian!(out, symmetry_f, sx) + @test out == symmetric_static + @test out == transpose(out) + out = similar(x, 9, 9) ForwardDiff.hessian!(out, prod, sx, cfg) @test out == actual @@ -156,6 +178,50 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) @test DiffResults.hessian(sresult3) == DiffResults.hessian(result) end +@testset "LowerTriangular, UpperTriangular and Diagonal" begin + for n in (3, 5), T in (LowerTriangular, UpperTriangular, Diagonal) + x = T(randn(n, n)) + xlen = ForwardDiff.structural_length(x) + weights = reshape(collect(1.0:n^2), n, n) + objective = x -> dot(weights, abs2.(x)) + expected = diagm(2 .* [weights[idx] for idx in ForwardDiff.structural_eachindex(x)]) + + H = ForwardDiff.hessian(objective, x) + @test size(H) == (xlen, xlen) + @test H == expected + + out = fill(NaN, xlen, xlen) + ForwardDiff.hessian!(out, objective, x) + @test out == expected + + flat = fill(NaN, xlen^2) + ForwardDiff.hessian!(flat, objective, x) + @test reshape(flat, xlen, xlen) == expected + end +end + +@testset "BigFloat with an unassigned input entry" begin + x = Vector{BigFloat}(undef, 10) + hole = 5 + for i in eachindex(x) + i == hole || (x[i] = BigFloat(i)) + end + used = [i for i in eachindex(x) if i != hole] + f(x) = sum(abs2(x[i]) for i in used) + expected = zeros(BigFloat, 10, 10) + for i in used + expected[i, i] = 2 + end + + @test !isassigned(x, hole) + for chunksize in (1, 2, 10) + cfg = ForwardDiff.HessianConfig(f, x, ForwardDiff.Chunk{chunksize}()) + H = ForwardDiff.hessian(f, x, cfg) + @test H isa Matrix{BigFloat} + @test H == expected + end +end + @testset "branches in dot" begin # https://github.com/JuliaDiff/ForwardDiff.jl/issues/551 H = [1 2 3; 4 5 6; 7 8 9]; diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index b6d36180..adc63cd7 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -198,6 +198,10 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) ForwardDiff.jacobian!(out, _diff, sx, scfg) @test out == actual + out = similar(x, 6, 9) + ForwardDiff.jacobian!(out, _diff, sx, scfg, Val{false}()) + @test out == actual + result = DiffResults.JacobianResult(similar(x, 6), x) result = ForwardDiff.jacobian!(result, _diff, x) diff --git a/test/SeedTest.jl b/test/SeedTest.jl index 02b821c3..90ef1858 100644 --- a/test/SeedTest.jl +++ b/test/SeedTest.jl @@ -90,4 +90,22 @@ end end end +@testset "seed_hessian_chunk!: $(nameof(typeof(x)))" for (x, sidx) in SEED_CASES + cfg = ForwardDiff.HessianConfig(nothing, x, ForwardDiff.Chunk{3}()) + duals = cfg.gradient_config.duals + iseeds = cfg.jacobian_config.seeds + oseeds = cfg.gradient_config.seeds + nstruct = length(sidx) + + ForwardDiff.seed_hessian_chunk!(duals, x, 1, nothing, nothing, nstruct) + ForwardDiff.seed_hessian_chunk!(duals, x, 4, iseeds, oseeds) + @test [i for (i, idx) in enumerate(sidx) if !iszero(ForwardDiff.partials(ForwardDiff.value(duals[idx])))] == collect(4:6) + @test [i for (i, idx) in enumerate(sidx) if !iszero(ForwardDiff.partials(duals[idx]))] == collect(4:6) + @test all(idx -> ForwardDiff.value(ForwardDiff.value(duals[idx])) == x[idx], eachindex(x)) + + ForwardDiff.seed_hessian_chunk!(duals, x, 4, nothing, nothing) + @test all(idx -> iszero(ForwardDiff.partials(ForwardDiff.value(duals[idx]))), sidx) + @test all(idx -> iszero(ForwardDiff.partials(duals[idx])), sidx) +end + end # module From 0b480dd3c0fcabdac9b02bbe9f88c0d8615981eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Fri, 21 Aug 2026 15:03:58 +0200 Subject: [PATCH 3/6] Hold in the Hessian config what the symmetric sweep reads `HessianConfig` wrapped a `JacobianConfig` and a `GradientConfig` because the Hessian was `jacobian(gradient(f), x)`: the outer sweep seeded the Jacobian config's buffer and the inner `gradient` seeded the gradient config's. The symmetric sweep seeds both layers of the one nested buffer, so it reads four things -- the two seed tuples, the nested buffer and its positions -- and never touches the Jacobian config's buffer again after the constructor derives the nested element type from it. Holding those four directly drops the dead buffer, and `checkstructure(cfg, x)` now resolves through the generic `AbstractConfig` method like every other config's. At `length(x) == 1000` and a chunk size of 12: HessianConfig(f, x) 1.472 MB -> 1.368 MB HessianConfig(f, result, x) 1.576 MB -> 1.368 MB The result-aware constructor allocated two dead buffers rather than one, since it built the `f!(y, x)` `JacobianConfig`. Nothing about the work buffers depends on `result`, so it forwards to the plain constructor and the two now return the same type -- which is what the tests asserting the two configs interchangeable already implied. Co-Authored-By: Claude Opus 5 (1M context) --- src/config.jl | 45 +++++++++++++++++++---------------------- src/hessian.jl | 14 ++++++------- test/AllocationsTest.jl | 34 ++++++++++++++++++------------- test/SeedTest.jl | 5 +---- 4 files changed, 49 insertions(+), 49 deletions(-) diff --git a/src/config.jl b/src/config.jl index fade0939..6f167500 100644 --- a/src/config.jl +++ b/src/config.jl @@ -206,9 +206,11 @@ Base.eltype(::Type{JacobianConfig{T,V,N,D,I}}) where {T,V,N,D,I} = Dual{T,V,N} # HessianConfig # ################# -struct HessianConfig{T,V,N,DG,DJ,IG,IJ} <: AbstractConfig{N} - jacobian_config::JacobianConfig{T,V,N,DJ,IJ} - gradient_config::GradientConfig{T,Dual{T,V,N},N,DG,IG} +struct HessianConfig{T,V,N,D,I} <: AbstractConfig{N} + iseeds::NTuple{N,Partials{N,V}} + oseeds::NTuple{N,Partials{N,Dual{T,V,N}}} + duals::D + indices::I end """ @@ -230,11 +232,13 @@ This constructor does not store/modify `x`. """ function HessianConfig(f::F, x::AbstractArray{V}, - chunk::Chunk = Chunk(x), - tag = Tag(f, V)) where {F,V} - jacobian_config = JacobianConfig(f, x, chunk, tag) - gradient_config = GradientConfig(f, jacobian_config.duals, chunk, tag) - return HessianConfig(jacobian_config, gradient_config) + ::Chunk{N} = Chunk(x), + ::T = Tag(f, V)) where {F,V,N,T} + iseeds = construct_seeds(Partials{N,V}) + oseeds = construct_seeds(Partials{N,Dual{T,V,N}}) + duals = similar(x, Dual{T,Dual{T,V,N},N}) + indices = structural_indices(duals) + return HessianConfig{T,V,N,typeof(duals),typeof(indices)}(iseeds, oseeds, duals, indices) end """ @@ -243,27 +247,20 @@ end Return a `HessianConfig` instance based on the type of `f`, types/storage in `result`, and type/shape of the input vector `x`. -The returned `HessianConfig` instance contains all the work buffers required by -`ForwardDiff.hessian` and `ForwardDiff.hessian!`. It is interchangeable with a config -constructed via `ForwardDiff.HessianConfig(f, x, chunk)`; this constructor retains the -result-aware form for compatibility. +Equivalent to `ForwardDiff.HessianConfig(f, x, chunk)`: the work buffers do not depend on +`result`. The result-aware form is retained for compatibility. If `f` is `nothing` instead of the actual target function, then the returned instance can be used with any target function. However, this will reduce ForwardDiff's ability to catch and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83). -This constructor does not store/modify `x`. +This constructor does not store/modify `result` or `x`. """ -function HessianConfig(f::F, - result::DiffResult, - x::AbstractArray{V}, - chunk::Chunk = Chunk(x), - tag = Tag(f, V)) where {F,V} - jacobian_config = JacobianConfig((f,gradient), DiffResults.gradient(result), x, chunk, tag) - gradient_config = GradientConfig(f, jacobian_config.duals[2], chunk, tag) - return HessianConfig(jacobian_config, gradient_config) -end +HessianConfig(f::F, + ::DiffResult, + x::AbstractArray{V}, + chunk::Chunk = Chunk(x), + tag = Tag(f, V)) where {F,V} = HessianConfig(f, x, chunk, tag) checktag(::HessianConfig{T},f,x) where {T} = checktag(T,f,x) -Base.eltype(::Type{HessianConfig{T,V,N,DG,DJ,IG,IJ}}) where {T,V,N,DG,DJ,IG,IJ} = - Dual{T,Dual{T,V,N},N} +Base.eltype(::Type{HessianConfig{T,V,N,D,I}}) where {T,V,N,D,I} = Dual{T,Dual{T,V,N},N} diff --git a/src/hessian.jl b/src/hessian.jl index 57270812..2b3a927d 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -16,7 +16,7 @@ Set `check` to `Val{false}()` to disable tag checking. This can lead to perturba function hessian(f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {F, T,CHK} require_one_based_indexing(x) CHK && checktag(T, f, x) - checkstructure(cfg.gradient_config, x) + checkstructure(cfg, x) H, _ = symmetric_hessian(f, x, cfg, nothing) return H end @@ -35,7 +35,7 @@ Set `check` to `Val{false}()` to disable tag checking. This can lead to perturba function hessian!(result::AbstractArray, f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {F,T,CHK} require_one_based_indexing(result, x) CHK && checktag(T, f, x) - checkstructure(cfg.gradient_config, x) + checkstructure(cfg, x) hlen = length(x) H = result isa AbstractMatrix ? result : reshape(result, hlen, hlen) symmetric_hessian!(H, f, x, cfg, nothing) @@ -55,7 +55,7 @@ Set `check` to `Val{false}()` to disable tag checking. This can lead to perturba function hessian!(result::DiffResult, f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, result, x), ::Val{CHK}=Val{true}()) where {F,T,CHK} require_one_based_indexing(x) CHK && checktag(T, f, x) - checkstructure(cfg.gradient_config, x) + checkstructure(cfg, x) hlen = length(x) hess = DiffResults.hessian(result) H = hess isa AbstractMatrix ? hess : reshape(hess, hlen, hlen) @@ -107,10 +107,10 @@ function symmetric_hessian_expr(result_definition::Expr) # output type and value. nblocks = xlen == 0 ? 1 : cld(xlen, N) - xdual = cfg.gradient_config.duals - indices = cfg.gradient_config.indices - iseeds = cfg.jacobian_config.seeds - oseeds = cfg.gradient_config.seeds + xdual = cfg.duals + indices = cfg.indices + iseeds = cfg.iseeds + oseeds = cfg.oseeds # The first evaluation determines the output type. Seeding the first block and clearing # the untouched tail partitions the fresh buffer, so every element is initialized once. diff --git a/test/AllocationsTest.jl b/test/AllocationsTest.jl index 65776de1..a4f87d6c 100644 --- a/test/AllocationsTest.jl +++ b/test/AllocationsTest.jl @@ -11,9 +11,7 @@ convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,F @testset "Test seed!/seed_zero_partials! allocations" begin x = rand(1000) cfg = ForwardDiff.GradientConfig(nothing, x) - duals = cfg.duals - seeds = cfg.seeds - indices = cfg.indices + (; duals, seeds, indices) = cfg allocs_seed!(args...) = @allocated ForwardDiff.seed!(args...) allocs_seed!(duals, x, indices, seeds) @@ -29,22 +27,30 @@ convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,F allocs_szp!(duals, x, indices, 1, 4) @test iszero(allocs_szp!(duals, x, indices, 1, 4)) - hcfg = ForwardDiff.HessianConfig(nothing, x) - hduals = hcfg.gradient_config.duals - hindices = hcfg.gradient_config.indices - iseeds = hcfg.jacobian_config.seeds - oseeds = hcfg.gradient_config.seeds - allocs_hseed!(args...) = @allocated ForwardDiff.seed_hessian_chunk!(args...) - allocs_hseed!(hduals, x, hindices, 1, iseeds, oseeds) - @test iszero(allocs_hseed!(hduals, x, hindices, 1, iseeds, oseeds)) - allocs_hseed!(hduals, x, hindices, 1, nothing, nothing, 4) - @test iszero(allocs_hseed!(hduals, x, hindices, 1, nothing, nothing, 4)) - allocs_convert_test_574() = @allocated convert_test_574() allocs_convert_test_574() @test iszero(allocs_convert_test_574()) end +@testset "Test seed_hessian_chunk! allocations" begin + x = rand(1000) + cfg = ForwardDiff.HessianConfig(nothing, x) + (; duals, indices, iseeds, oseeds) = cfg + + allocs_hseed!(args...) = @allocated ForwardDiff.seed_hessian_chunk!(args...) + # all four seed combinations, the mixed ones being what the off-diagonal blocks use + @testset "iseeds=$(i === nothing) oseeds=$(o === nothing)" for (i, o) in + ((iseeds, oseeds), + (iseeds, nothing), + (nothing, oseeds), + (nothing, nothing)) + allocs_hseed!(duals, x, indices, 1, i, o) + @test iszero(allocs_hseed!(duals, x, indices, 1, i, o)) + allocs_hseed!(duals, x, indices, 1, i, o, 4) + @test iszero(allocs_hseed!(duals, x, indices, 1, i, o, 4)) + end +end + @testset "Test jacobian! allocations" begin # jacobian! should not allocate when called with a pre-allocated result Matrix. # Previously, reshape() inside extract_jacobian! allocated a wrapper diff --git a/test/SeedTest.jl b/test/SeedTest.jl index 13076bae..8fc55d39 100644 --- a/test/SeedTest.jl +++ b/test/SeedTest.jl @@ -134,10 +134,7 @@ end @testset "seed_hessian_chunk!: $(nameof(typeof(x)))" for (x, sidx) in SEED_CASES cfg = ForwardDiff.HessianConfig(nothing, x, ForwardDiff.Chunk{3}()) - duals = cfg.gradient_config.duals - indices = cfg.gradient_config.indices - iseeds = cfg.jacobian_config.seeds - oseeds = cfg.gradient_config.seeds + (; duals, indices, iseeds, oseeds) = cfg nstruct = length(sidx) ForwardDiff.seed_hessian_chunk!(duals, x, indices, 1, nothing, nothing, nstruct) From dcdf668fd20e14328b7b42963964119405d2fc2e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Fri, 21 Aug 2026 15:53:44 +0200 Subject: [PATCH 4/6] Restore the result shape check dropped with `extract_jacobian!` The symmetric sweep writes the result entry by entry, which lost the validation that `reshape_jacobian` and the broadcast in `extract_jacobian!` used to provide. A vector result was still checked by `reshape`, a matrix one no longer was: hessian!(fill(NaN, 4, 4), f, rand(3)) # no error, row/col 4 left NaN hessian!(fill(NaN, 4, 4), f, SVector(1., 2., 3.)) `reshape_hessian` mirrors `reshape_jacobian`, down to its `DiffResult` method, so the `DiffResult` path is checked too -- it holds a buffer no entry point ever passes to `require_one_based_indexing`, hence the extra call here. Co-Authored-By: Claude Opus 5 (1M context) --- ext/ForwardDiffStaticArraysExt.jl | 2 +- src/hessian.jl | 19 ++++++++++++------- test/HessianTest.jl | 8 ++++++++ 3 files changed, 21 insertions(+), 8 deletions(-) diff --git a/ext/ForwardDiffStaticArraysExt.jl b/ext/ForwardDiffStaticArraysExt.jl index 9e52cc44..ae970102 100644 --- a/ext/ForwardDiffStaticArraysExt.jl +++ b/ext/ForwardDiffStaticArraysExt.jl @@ -131,7 +131,7 @@ ForwardDiff.hessian(f::F, x::StaticArray, cfg::HessianConfig, ::Val) where {F} = T = typeof(Tag(f, eltype(x))) ydual = f(dualize(T, dualize(T, x))) ydual isa Real || throw(HESSIAN_ERROR) - H = result isa AbstractMatrix ? result : reshape(result, length(x), length(x)) + H = ForwardDiff.reshape_hessian(result, x) ForwardDiff.extract_hessian_chunk!(T, H, ydual, structural_indices(x), 0, 0, length(x), length(x)) return result end diff --git a/src/hessian.jl b/src/hessian.jl index 2b3a927d..e273bd9a 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -36,9 +36,7 @@ function hessian!(result::AbstractArray, f::F, x::AbstractArray, cfg::HessianCon require_one_based_indexing(result, x) CHK && checktag(T, f, x) checkstructure(cfg, x) - hlen = length(x) - H = result isa AbstractMatrix ? result : reshape(result, hlen, hlen) - symmetric_hessian!(H, f, x, cfg, nothing) + symmetric_hessian!(reshape_hessian(result, x), f, x, cfg, nothing) return result end @@ -56,10 +54,8 @@ function hessian!(result::DiffResult, f::F, x::AbstractArray, cfg::HessianConfig require_one_based_indexing(x) CHK && checktag(T, f, x) checkstructure(cfg, x) - hlen = length(x) - hess = DiffResults.hessian(result) - H = hess isa AbstractMatrix ? hess : reshape(hess, hlen, hlen) - _, ydual = symmetric_hessian!(H, f, x, cfg, DiffResults.gradient(result)) + _, ydual = symmetric_hessian!(reshape_hessian(result, x), f, x, cfg, + DiffResults.gradient(result)) result = DiffResults.value!(result, value(T, value(T, ydual))) return result end @@ -70,6 +66,15 @@ end const HESSIAN_ERROR = DimensionMismatch("hessian(f, x) expects that f(x) is a real number. Perhaps you meant jacobian(f, x)?") +function reshape_hessian(result::AbstractMatrix, x) + require_one_based_indexing(result) + size(result) == (length(x), length(x)) || throw(DimensionMismatch( + lazy"cannot store the $(length(x))×$(length(x)) Hessian in a result of size $(size(result))")) + return result +end +reshape_hessian(result::AbstractArray, x) = reshape(result, length(x), length(x)) +reshape_hessian(result::DiffResult, x) = reshape_hessian(DiffResults.hessian(result), x) + # Copy a block from the nested partials and fill its transpose. On diagonal blocks, read # only the upper triangle so the result is exactly symmetric. Both axes are indexed by the linear # indices of `x`, as the columns of a Jacobian are, so `indices` gives the row and column of a block. diff --git a/test/HessianTest.jl b/test/HessianTest.jl index e686d3ac..8a55d807 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -67,6 +67,14 @@ cfgx = ForwardDiff.HessianConfig(sin, x) @test_throws DimensionMismatch ForwardDiff.hessian(identity, x) @test_throws DimensionMismatch ForwardDiff.hessian!(similar(x, 3, 3), identity, x) +@testset "wrongly sized result: $(nameof(typeof(z)))" for z in ([1.0, 2.0, 3.0], + SVector(1.0, 2.0, 3.0)) + msg = "DimensionMismatch: cannot store the 3×3 Hessian in a result of size (4, 4)" + @test_throws msg ForwardDiff.hessian!(fill(NaN, 4, 4), prod, z) + @test_throws msg ForwardDiff.hessian!(DiffResults.DiffResult(0.0, zeros(3), + fill(NaN, 4, 4)), prod, z) +end + ######################## # test vs. Calculus.jl # From c9d6a014bed82c416659505df0cc44bd905680f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Fri, 21 Aug 2026 16:28:15 +0200 Subject: [PATCH 5/6] Read the same triangle at every chunk size Diagonal blocks read `outer <= inner`, off-diagonal blocks read `outer > inner`, and the StaticArrays path reads `outer <= inner` throughout. Reading an entry with i in the outer layer rounds differently from reading it with j there, so the result was not reproducible: at `n = 16`, chunk sizes 1, 2, 3, 5, 7, 11 each differed from chunk 16 and from the `SVector` path by ~1e-16. Swapping which layer block q carries fixes it. Same number of evaluations and seed writes, and q is still seeded once outside the loop; afterwards every chunk size is bitwise identical to the `SVector` path. `log(sum(exp, z))` is the objective in the new test because its mixed partials actually round differently in the two orders -- `sum(z)^3`, `exp(sum(z))`, `prod(z)` and `sum(sin, z) * sum(cos, z)` all give bitwise equal results either way, so none of them would have caught this. It also makes `symmetric_static == hessian(symmetry_f, x)` robust rather than accidental: that only passed because `n = 9` is below `DEFAULT_CHUNK_THRESHOLD`, so the array path ran a single block. Co-Authored-By: Claude Opus 5 (1M context) --- src/hessian.jl | 12 ++++++------ test/HessianTest.jl | 12 ++++++++++++ 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/src/hessian.jl b/src/hessian.jl index e273bd9a..a5394600 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -134,17 +134,17 @@ function symmetric_hessian_expr(result_definition::Expr) for q in 2:nblocks qoffset = (q - 1) * N qsize = min(N, xlen - qoffset) - # Off-diagonal blocks: p seeds columns and q seeds rows. The outer seeds for q - # remain unchanged throughout this loop. - seed_hessian_chunk!(xdual, x, indices, qoffset + 1, nothing, oseeds, qsize) + # Outer-i inner-j and outer-j inner-i round differently, so the outer layer always + # takes the earlier position -- else the result would depend on the chunk size. + seed_hessian_chunk!(xdual, x, indices, qoffset + 1, iseeds, nothing, qsize) for p in 1:(q - 1) poffset = (p - 1) * N - seed_hessian_chunk!(xdual, x, indices, poffset + 1, iseeds, nothing) + seed_hessian_chunk!(xdual, x, indices, poffset + 1, nothing, oseeds) ydual = f(xdual) - extract_hessian_chunk!(T, H, ydual, indices, qoffset, poffset, qsize, N) + extract_hessian_chunk!(T, H, ydual, indices, poffset, qoffset, N, qsize) seed_hessian_chunk!(xdual, x, indices, poffset + 1, nothing, nothing) end - # The diagonal block adds q's inner seeds while retaining its outer seeds. + # The diagonal block adds q's outer seeds while retaining its inner seeds. seed_hessian_chunk!(xdual, x, indices, qoffset + 1, iseeds, oseeds, qsize) ydual = f(xdual) extract_hessian_chunk!(T, H, ydual, indices, qoffset, qoffset, qsize, qsize) diff --git a/test/HessianTest.jl b/test/HessianTest.jl index 8a55d807..ccd231e7 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -186,6 +186,18 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) @test DiffResults.hessian(sresult3) == DiffResults.hessian(result) end +# `log(sum(exp, z))` rounds differently in the two nesting orders, hence the bitwise comparison +@testset "chunk size independence" begin + n = 16 + x = randn(n) + f = z -> log(sum(exp, z)) + expected = ForwardDiff.hessian(f, SVector{n}(x)) + @testset "chunk size = $c" for c in (1, 2, 3, 5, 7, 11, n) + cfg = ForwardDiff.HessianConfig(f, x, ForwardDiff.Chunk{c}()) + @test ForwardDiff.hessian(f, x, cfg) == expected + end +end + # `n = 5` is the only structured case whose default `Chunk(x)` reaches the off-diagonal blocks. @testset "structured inputs: $(nameof(W)) of size $n" for n in (3, 5), (W, sidx) in ( From a82a9e6141183596a76ac46ac39c19ea0cf71288 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Fri, 21 Aug 2026 18:54:05 +0200 Subject: [PATCH 6/6] Build only the zero a layer needs, and cover the paths left untested MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `seed_hessian_chunk!` materialised both zeros even when both seeds were supplied. Free for an isbits value type, but for `BigFloat` at `N = 3` it cost 288 bytes per call, identical for all four seed combinations. Each `=== nothing` is a compile-time constant, so deciding per layer folds away: iseeds oseeds before after given given 288 0 given - 288 192 - given 288 96 - - 288 288 Also documents the result shape, in the wording `jacobian` uses for its own, and adds a sentence on why the clearing calls cannot go through `seed_zero_partials!` and why the `Partials{0}` method of `extract_hessian` is load-bearing -- for a constant `f` the generic method builds a `0 × length(x)` result, not `length(x) × length(x)`. New coverage: the extension's `reshape` branch and its two `HESSIAN_ERROR` throws, and empty inputs on both paths. Not added: the mixed seed forms in `SeedTest`. Which layer a seed lands in is enforced by the types -- `oseeds` only fits the outer `Dual` -- and both mixed forms run in every multi-block sweep, so the bitwise chunk-size test covers them with a failure mode that a unit test would only relocate. This is unlike `seed_zero_partials!`, whose testset exists because over-clearing is invisible through the public API. Co-Authored-By: Claude Opus 5 (1M context) --- ext/ForwardDiffStaticArraysExt.jl | 2 ++ src/apiutils.jl | 7 ++++--- src/hessian.jl | 3 +++ test/AllocationsTest.jl | 8 ++++++++ test/HessianTest.jl | 12 ++++++++++++ 5 files changed, 29 insertions(+), 3 deletions(-) diff --git a/ext/ForwardDiffStaticArraysExt.jl b/ext/ForwardDiffStaticArraysExt.jl index ae970102..1ef9dbc0 100644 --- a/ext/ForwardDiffStaticArraysExt.jl +++ b/ext/ForwardDiffStaticArraysExt.jl @@ -112,6 +112,8 @@ end return typeof(H)(Symmetric(H, :U)) end +# An `f` ignoring its argument returns no partials at all, not `length(x)` zero ones, so the method +# above would build a result with no rows. Reached for an empty `x` too. @inline function extract_hessian(::Type{T}, ydual::Partials{0}, x::S) where {T,S<:StaticArray} R = StaticArrays.similar_type(S, valtype(T, eltype(ydual)), Size(length(x), length(x))) return zero(R) diff --git a/src/apiutils.jl b/src/apiutils.jl index f9a1724b..3e749156 100644 --- a/src/apiutils.jl +++ b/src/apiutils.jl @@ -175,13 +175,14 @@ function seed!(duals::AbstractArray{Dual{T,V,N}}, x, indices, index, end end -# Seed a chunk in either layer of nested duals. A `nothing` seed clears that layer. +# Seed a chunk in either layer of nested duals. A `nothing` seed clears that layer; +# `seed_zero_partials!` cannot, as it would pass the primal where a nested `Dual` is wanted. function seed_hessian_chunk!(duals::AbstractArray{Dual{T,Dual{T,V,N},N}}, x, indices, index, iseeds::Union{Nothing,NTuple{N,Partials{N,V}}}, oseeds::Union{Nothing,NTuple{N,Partials{N,Dual{T,V,N}}}}, chunksize = N) where {T,V,N} - izero = zero(Partials{N,V}) - ozero = zero(Partials{N,Dual{T,V,N}}) + izero = iseeds === nothing ? zero(Partials{N,V}) : nothing + ozero = oseeds === nothing ? zero(Partials{N,Dual{T,V,N}}) : nothing return _seed!(duals, x, structural_chunk(indices, index, chunksize)) do value, i inner = Dual{T,V,N}(value, iseeds === nothing ? izero : iseeds[i]) Dual{T,Dual{T,V,N},N}(inner, oseeds === nothing ? ozero : oseeds[i]) diff --git a/src/hessian.jl b/src/hessian.jl index a5394600..ddd59043 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -6,6 +6,9 @@ ForwardDiff.hessian(f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, x), check=Val{true}()) Return `H(f)` evaluated at `x`, assuming `f` is called as `f(x)`. +Multidimensional arrays are flattened in iteration order: the array +`H(f)` has shape `length(x) × length(x)`, and its elements are +`H(f)[j,k] = ∂²f(x)/∂x[j]∂x[k]`. The returned Hessian is exactly symmetric: its two triangles are filled from the same derivative values. diff --git a/test/AllocationsTest.jl b/test/AllocationsTest.jl index a4f87d6c..3108039d 100644 --- a/test/AllocationsTest.jl +++ b/test/AllocationsTest.jl @@ -49,6 +49,14 @@ end allocs_hseed!(duals, x, indices, 1, i, o, 4) @test iszero(allocs_hseed!(duals, x, indices, 1, i, o, 4)) end + + # a zero of a non-isbits value type does allocate, so supplying both seeds must not build one + @testset "BigFloat" begin + y = BigFloat[1, 2, 3] + cfg = ForwardDiff.HessianConfig(nothing, y, ForwardDiff.Chunk{3}()) + allocs_hseed!(cfg.duals, y, cfg.indices, 1, cfg.iseeds, cfg.oseeds) + @test iszero(allocs_hseed!(cfg.duals, y, cfg.indices, 1, cfg.iseeds, cfg.oseeds)) + end end @testset "Test jacobian! allocations" begin diff --git a/test/HessianTest.jl b/test/HessianTest.jl index ccd231e7..47c1446b 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -132,6 +132,12 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) @test symmetric_static == ForwardDiff.hessian(symmetry_f, x) @test all(iszero, ForwardDiff.hessian(Returns(2.0), sx)) @test_throws DimensionMismatch ForwardDiff.hessian(identity, sx) + @test_throws DimensionMismatch ForwardDiff.hessian!(similar(x, 9, 9), identity, sx) + @test_throws DimensionMismatch ForwardDiff.hessian!(DiffResults.HessianResult(sx), identity, sx) + + flat = fill(NaN, 81) + @test ForwardDiff.hessian!(flat, prod, sx) === flat + @test reshape(flat, 9, 9) == actual out = similar(x, 9, 9) ForwardDiff.hessian!(out, prod, sx) @@ -186,6 +192,12 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) @test DiffResults.hessian(sresult3) == DiffResults.hessian(result) end +@testset "empty input: $(nameof(typeof(z)))" for z in (Float64[], SVector{0,Float64}()) + @test ForwardDiff.hessian(sum, z) == zeros(0, 0) + out = fill(NaN, 0, 0) + @test ForwardDiff.hessian!(out, sum, z) === out +end + # `log(sum(exp, z))` rounds differently in the two nesting orders, hence the bitwise comparison @testset "chunk size independence" begin n = 16