From 0091edd55c4ea8fd0acd93cdfc883e357e963058 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Sun, 16 Aug 2026 10:05:43 +0200 Subject: [PATCH 01/13] 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 02/13] 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 e6c6a45594affbfc151f8f2bddfc73df382332af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Mon, 31 Aug 2026 11:28:19 +0200 Subject: [PATCH 03/13] 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 three things -- the two seed tuples and the nested buffer -- and never touches the Jacobian config's buffer again after the constructor derives the nested element type from it. Holding those three directly drops the dead buffer. At `length(x) == 1000` and a chunk size of 12: HessianConfig(f, x) 1.422 MB -> 1.312 MB HessianConfig(f, result, x) 1.531 MB -> 1.312 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 | 43 ++++++++++++++++++----------------------- src/hessian.jl | 6 +++--- test/AllocationsTest.jl | 6 +++--- test/SeedTest.jl | 5 +---- 4 files changed, 26 insertions(+), 34 deletions(-) diff --git a/src/config.jl b/src/config.jl index 58c145f3..5469b8c9 100644 --- a/src/config.jl +++ b/src/config.jl @@ -195,9 +195,10 @@ Base.eltype(::Type{JacobianConfig{T,V,N,D}}) where {T,V,N,D} = Dual{T,V,N} # HessianConfig # ################# -struct HessianConfig{T,V,N,DG,DJ} <: AbstractConfig{N} - jacobian_config::JacobianConfig{T,V,N,DJ} - gradient_config::GradientConfig{T,Dual{T,V,N},N,DG} +struct HessianConfig{T,V,N,D} <: AbstractConfig{N} + iseeds::NTuple{N,Partials{N,V}} + oseeds::NTuple{N,Partials{N,Dual{T,V,N}}} + duals::D end """ @@ -219,11 +220,12 @@ 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}) + return HessianConfig{T,V,N,typeof(duals)}(iseeds, oseeds, duals) end """ @@ -232,27 +234,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}}) where {T,V,N,DG,DJ} = - Dual{T,Dual{T,V,N},N} +Base.eltype(::Type{HessianConfig{T,V,N,D}}) where {T,V,N,D} = Dual{T,Dual{T,V,N},N} diff --git a/src/hessian.jl b/src/hessian.jl index 489b9014..3dd81511 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -100,9 +100,9 @@ function symmetric_hessian_expr(result_definition::Expr) # 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 + xdual = cfg.duals + 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 3a59a5ad..146b7fc3 100644 --- a/test/AllocationsTest.jl +++ b/test/AllocationsTest.jl @@ -30,9 +30,9 @@ convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,F @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 + hduals = hcfg.duals + iseeds = hcfg.iseeds + oseeds = hcfg.oseeds 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)) diff --git a/test/SeedTest.jl b/test/SeedTest.jl index 90ef1858..d7829e6c 100644 --- a/test/SeedTest.jl +++ b/test/SeedTest.jl @@ -91,10 +91,7 @@ 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 + (; duals, iseeds, oseeds) = ForwardDiff.HessianConfig(nothing, x, ForwardDiff.Chunk{3}()) nstruct = length(sidx) ForwardDiff.seed_hessian_chunk!(duals, x, 1, nothing, nothing, nstruct) From 5d1d4753dd0f61c8c5d84a007e89b166df7eb5b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Mon, 31 Aug 2026 12:00:01 +0200 Subject: [PATCH 04/13] Read the same triangle at every chunk size Diagonal blocks read `outer <= inner`, off-diagonal blocks the other triangle, 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` and `f = log(sum(exp, z))`, chunk sizes 1, 2, 3, 5, 7 and 11 each differed from chunk 16 and from the `SVector` path by ~3e-18. 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 the existing `symmetric_static == hessian(symmetry_f, x)` assertion 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 | 13 +++++++------ test/HessianTest.jl | 15 +++++++++++++++ 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/src/hessian.jl b/src/hessian.jl index 3dd81511..b02ce6c8 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -118,17 +118,18 @@ 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, 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. + # q's inner seeds remain unchanged throughout this loop. + seed_hessian_chunk!(xdual, x, qoffset + 1, iseeds, nothing, 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, poffset + 1, nothing, oseeds) ydual = f(xdual) - extract_hessian_chunk!(T, H, ydual, qoffset, poffset, qsize, N) + extract_hessian_chunk!(T, H, ydual, poffset, qoffset, N, qsize) seed_hessian_chunk!(xdual, x, 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, qoffset + 1, iseeds, oseeds, qsize) ydual = f(xdual) extract_hessian_chunk!(T, H, ydual, qoffset, qoffset, qsize, qsize) diff --git a/test/HessianTest.jl b/test/HessianTest.jl index 119fd14d..b3881f48 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -178,6 +178,21 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) @test DiffResults.hessian(sresult3) == DiffResults.hessian(result) end +@testset "the result does not depend on the chunk size" begin + n = 16 + v = randn(n) + # the mixed partials of `log(sum(exp, ·))` round differently in the two orders; `sum(z)^3`, + # `exp(sum(z))`, `prod(z)` and `sum(sin, z) * sum(cos, z)` do not, and would pass regardless + f = z -> log(sum(exp, z)) + expected = Matrix(ForwardDiff.hessian(f, SVector{n}(v))) + + @testset "chunk size = $c" for c in (1, 2, 3, 5, 7, 11, n) + H = ForwardDiff.hessian(f, v, ForwardDiff.HessianConfig(f, v, ForwardDiff.Chunk{c}())) + @test H == expected + @test H == transpose(H) + end +end + @testset "LowerTriangular, UpperTriangular and Diagonal" begin for n in (3, 5), T in (LowerTriangular, UpperTriangular, Diagonal) x = T(randn(n, n)) From 705e7b141a95ec83297ee2a39cdb150666dfe482 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Mon, 31 Aug 2026 12:20:00 +0200 Subject: [PATCH 05/13] 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 New coverage: the extension's `reshape` branch and its two `HESSIAN_ERROR` throws, and empty inputs on both paths. Also notes why the `Partials{0}` method of `extract_hessian` is load-bearing -- for a constant `f` the generic method would build a `0 × length(x)` result, not `length(x) × length(x)`. 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 ++++--- test/AllocationsTest.jl | 13 +++++++++++-- test/HessianTest.jl | 24 ++++++++++++++++++++++++ 4 files changed, 41 insertions(+), 5 deletions(-) diff --git a/ext/ForwardDiffStaticArraysExt.jl b/ext/ForwardDiffStaticArraysExt.jl index 26abf43b..e3daa34c 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 1d54d7bb..c75bdc9e 100644 --- a/src/apiutils.jl +++ b/src/apiutils.jl @@ -129,13 +129,14 @@ function seed!(duals::AbstractArray{Dual{T,V,N}}, x, 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, 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 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]) diff --git a/test/AllocationsTest.jl b/test/AllocationsTest.jl index 146b7fc3..d773df20 100644 --- a/test/AllocationsTest.jl +++ b/test/AllocationsTest.jl @@ -34,11 +34,20 @@ convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,F iseeds = hcfg.iseeds oseeds = hcfg.oseeds 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)) + for i in (iseeds, nothing), o in (oseeds, nothing) + allocs_hseed!(hduals, x, 1, i, o) + @test iszero(allocs_hseed!(hduals, x, 1, i, o)) + end allocs_hseed!(hduals, x, 1, nothing, nothing, 4) @test iszero(allocs_hseed!(hduals, x, 1, nothing, nothing, 4)) + # a zero is free for an isbits value type, so only a `BigFloat` catches one being built for a + # layer that was given seeds + bx = BigFloat.(x) + bcfg = ForwardDiff.HessianConfig(nothing, bx, ForwardDiff.Chunk{3}()) + allocs_hseed!(bcfg.duals, bx, 1, bcfg.iseeds, bcfg.oseeds) + @test iszero(allocs_hseed!(bcfg.duals, bx, 1, bcfg.iseeds, bcfg.oseeds)) + allocs_convert_test_574() = @allocated convert_test_574() allocs_convert_test_574() @test iszero(allocs_convert_test_574()) diff --git a/test/HessianTest.jl b/test/HessianTest.jl index b3881f48..bb068d9a 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -178,6 +178,30 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) @test DiffResults.hessian(sresult3) == DiffResults.hessian(result) end +@testset "empty input" begin + f = z -> 1.0 + for x in (Float64[], SVector{0,Float64}()) + H = ForwardDiff.hessian(f, x) + @test size(H) == (0, 0) + @test ForwardDiff.hessian!(fill(NaN, 0, 0), f, x) == H + end +end + +@testset "a StaticArray result that is not a matrix" begin + sx = SVector(1.0, 2.0, 3.0) + flat = fill(NaN, 9) + @test ForwardDiff.hessian!(flat, prod, sx) === flat + @test reshape(flat, 3, 3) == ForwardDiff.hessian(prod, sx) +end + +@testset "an array-valued f is not a Hessian" begin + sx = SVector(1.0, 2.0, 3.0) + @test_throws DimensionMismatch ForwardDiff.hessian(identity, sx) + @test_throws DimensionMismatch ForwardDiff.hessian!(fill(NaN, 3, 3), identity, sx) + @test_throws DimensionMismatch ForwardDiff.hessian!(DiffResults.HessianResult(sx), identity, sx) + @test_throws DimensionMismatch ForwardDiff.hessian(identity, [1.0, 2.0, 3.0]) +end + @testset "the result does not depend on the chunk size" begin n = 16 v = randn(n) From 907ccd1d8dfd46c8120edaad95b455b370eb0f68 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Mon, 31 Aug 2026 13:53:01 +0200 Subject: [PATCH 06/13] Index the Hessian by the linear indices of `x` For a structured input the result mixed two coordinate systems: rows were linear indices of `x`, columns were structural positions, so `hessian(f, UpperTriangular(3x3))` was 9x6 and a caller could not say what `H[4, 5]` referred to without reading `structural_eachindex`. Chunk mode did not get that far, throwing a `DimensionMismatch` from `reshape_jacobian`, and `hessian!(DiffResults.HessianResult(x), f, x)` failed to broadcast. Both axes are now the linear indices of `x`, with hard zeros in the rows and columns of the structural zeros. That is the smaller change relative to master, which already used linear indices for the rows and only got the columns wrong, and it makes `DiffResults.HessianResult(x)` the right buffer, since it allocates `length(x)^2`: size(hessian(f, UpperTriangular(randn(3, 3)))) # 9x9, rows/cols 2, 3, 6 zero size(hessian(f, Diagonal(randn(n)))) # n^2 x n^2 The cost is `Diagonal`, whose Hessian is now quartic in the size of the diagonal; differentiating with respect to the diagonal vector is the better choice there. `structural_linearindices` gives the positions of `structural_eachindex` as linear indices, computed once per sweep. `Base.OneTo` and a `StepRange` for the dense and `Diagonal` cases, so those stay allocation-free; the triangles build one position vector per call, against a sweep that evaluates `f` B(B+1)/2 times. Its two-argument form takes the config's buffer and the input, and performs the size check `structural_eachindex` performs for seeding, which the sweep no longer goes through. `extract_hessian_gradient_chunk!` no longer delegates to `extract_gradient_chunk!`, which takes its positions from the result rather than from `x` (#838): `HessianResult` hands back a dense `size(x)` gradient buffer even for a structured `x`, so delegating scattered the derivatives into the first `structural_length(x)` linear positions. Seeding a structured buffer now walks it by linear index rather than by `CartesianIndex`, which is also what lets a position double as a Hessian row. `Base._unsetindex!` is implemented for `Array` alone -- for a linear index its `AbstractArray` fallback recurses forever, where a `CartesianIndex` merely had no method -- so the unassigned-entry path goes through a wrapper that raises an `ArgumentError` naming the entry instead. That also removes a pre-existing `StackOverflowError` for a `Diagonal` of a non-bits element type, whose positions `structural_eachindex` already gave as linear indices. Not addressed here, and unchanged from master: `gradient!`/`jacobian!` extraction (#838, #839) and detecting a config reused across structures (#842). Co-Authored-By: Claude Opus 5 (1M context) --- ext/ForwardDiffStaticArraysExt.jl | 4 +- src/apiutils.jl | 51 ++++++++++++++++++++-- src/hessian.jl | 71 ++++++++++++++++++++----------- test/AllocationsTest.jl | 14 +++--- test/HessianTest.jl | 52 +++++++++++++++------- test/SeedTest.jl | 32 +++++++++----- 6 files changed, 162 insertions(+), 62 deletions(-) diff --git a/ext/ForwardDiffStaticArraysExt.jl b/ext/ForwardDiffStaticArraysExt.jl index e3daa34c..ac7f8215 100644 --- a/ext/ForwardDiffStaticArraysExt.jl +++ b/ext/ForwardDiffStaticArraysExt.jl @@ -5,7 +5,7 @@ using ForwardDiff.LinearAlgebra using ForwardDiff.DiffResults using ForwardDiff: Dual, partials, npartials, Partials, GradientConfig, JacobianConfig, HessianConfig, Tag, Chunk, gradient, hessian, jacobian, gradient!, hessian!, jacobian!, - extract_gradient!, extract_jacobian!, extract_value!, + extract_gradient!, extract_jacobian!, extract_value!, structural_linearindices, vector_mode_gradient, vector_mode_gradient!, vector_mode_jacobian, vector_mode_jacobian!, HESSIAN_ERROR, valtype, value using DiffResults: DiffResult, ImmutableDiffResult, MutableDiffResult @@ -134,7 +134,7 @@ ForwardDiff.hessian(f::F, x::StaticArray, cfg::HessianConfig, ::Val) where {F} = 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)) + ForwardDiff.extract_hessian_chunk!(T, H, ydual, structural_linearindices(x), 0, 0, length(x), length(x)) return result end diff --git a/src/apiutils.jl b/src/apiutils.jl index c75bdc9e..2039d8f3 100644 --- a/src/apiutils.jl +++ b/src/apiutils.jl @@ -70,6 +70,43 @@ function structural_eachindex(x::Diagonal, y::AbstractArray) return diagind(x) end +# The positions of `structural_eachindex`, in the same order, as linear indices of `x`. The two +# argument form is only ever given a config's work buffer and the input it is used with. +function check_structural_size(duals, x) + if size(duals) != size(x) + throw(DimensionMismatch(lazy"the config was built for an array of size $(size(duals)) and cannot be used with an array of size $(size(x))")) + end + return nothing +end + +structural_linearindices(x::AbstractArray) = structural_linearindices(x, x) +function structural_linearindices(duals::AbstractArray, x::AbstractArray) + require_one_based_indexing(duals, x) + check_structural_size(duals, x) + return Base.OneTo(length(duals)) +end +function structural_linearindices(duals::UpperTriangular, x::AbstractArray) + require_one_based_indexing(duals, x) + check_structural_size(duals, x) + n = size(duals, 1) + return [i + n * (j - 1) for j in 1:n for i in 1:j] +end +function structural_linearindices(duals::LowerTriangular, x::AbstractArray) + require_one_based_indexing(duals, x) + check_structural_size(duals, x) + n = size(duals, 1) + return [i + n * (j - 1) for j in 1:n for i in j:n] +end +function structural_linearindices(duals::Diagonal, x::AbstractArray) + require_one_based_indexing(duals, x) + check_structural_size(duals, x) + n = size(duals, 1) + return range(1; step = n + 1, length = n) +end + +# The `count` positions starting at structural position `index`. +structural_chunk(indices, index, count) = view(indices, index:(index + count - 1)) + # Copies the values of `x` into `duals` with zero partials. Used both to remove seeds `duals` is # currently carrying and to initialize a freshly allocated work buffer, whose elements must all be # written before the target function reads them. @@ -93,6 +130,14 @@ function _seed_zero_partials!(duals::AbstractArray{Dual{T,V,N}}, x, idxs) where end end +# `Base._unsetindex!` is implemented for `Array` alone: for a linear index its `AbstractArray` +# fallback recurses forever, and it has no `CartesianIndex` method at all. +_unsetindex!(duals::Array, idx) = Base._unsetindex!(duals, idx) +_unsetindex!(duals::AbstractArray, idx) = throw(ArgumentError(LazyString( + "cannot differentiate at an input with an unassigned entry at index ", idx, + ": that would leave an entry of the ", nameof(typeof(duals)), + " work buffer unassigned, which is only possible for an Array"))) + # 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} @@ -105,7 +150,7 @@ end if isassigned(x, idx) duals[idx] = make_dual(x[idx], i) else - Base._unsetindex!(duals, idx) + _unsetindex!(duals, idx) end end end @@ -131,13 +176,13 @@ end # 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, index, +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 = iseeds === nothing ? zero(Partials{N,V}) : nothing ozero = oseeds === nothing ? zero(Partials{N,Dual{T,V,N}}) : nothing - idxs = Iterators.take(Iterators.drop(structural_eachindex(duals, x), index - 1), chunksize) + idxs = structural_chunk(indices, index, 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]) diff --git a/src/hessian.jl b/src/hessian.jl index b02ce6c8..06d363a9 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)`. +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. @@ -34,8 +36,8 @@ 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) - xlen = structural_length(x) - H = result isa AbstractMatrix ? result : reshape(result, xlen, xlen) + hlen = length(x) + H = result isa AbstractMatrix ? result : reshape(result, hlen, hlen) symmetric_hessian!(H, f, x, cfg, nothing) return result end @@ -53,9 +55,9 @@ 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) - xlen = structural_length(x) + hlen = length(x) hess = DiffResults.hessian(result) - H = hess isa AbstractMatrix ? hess : reshape(hess, xlen, xlen) + H = hess isa AbstractMatrix ? hess : reshape(hess, hlen, hlen) _, ydual = symmetric_hessian!(H, f, x, cfg, DiffResults.gradient(result)) result = DiffResults.value!(result, value(T, value(T, ydual))) return result @@ -68,24 +70,35 @@ end const HESSIAN_ERROR = DimensionMismatch("hessian(f, x) expects that f(x) is a real number. Perhaps you meant jacobian(f, 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. -function extract_hessian_chunk!(::Type{T}, H, ydual, roffset, coffset, rsize, csize) where {T} +# only the upper triangle so the result is exactly symmetric. `indices` maps a block +# position to its row and column, both being linear indices of `x`. +function extract_hessian_chunk!(::Type{T}, H, ydual, indices, roffset, coffset, rsize, csize) where {T} + rows = structural_chunk(indices, roffset + 1, rsize) + cols = structural_chunk(indices, coffset + 1, csize) for r in 1:rsize drow = partials(T, ydual, r) + i = rows[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 + j = cols[c] + H[i, j] = h + H[j, i] = 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) +# TODO: delegate to `extract_gradient_chunk!` once it takes its positions from `x` (#838). +extract_hessian_gradient_chunk!(::Type{T}, ::Nothing, ydual, indices, index, chunksize) where {T} = nothing +function extract_hessian_gradient_chunk!(::Type{T}, grad, ydual, indices, index, chunksize) where {T} + dual = value(T, ydual) + for (i, idx) in enumerate(structural_chunk(indices, index, chunksize)) + grad[idx] = partials(T, dual, i) + end + return grad +end # 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). @@ -103,17 +116,25 @@ function symmetric_hessian_expr(result_definition::Expr) xdual = cfg.duals iseeds = cfg.iseeds oseeds = cfg.oseeds + indices = structural_linearindices(xdual, x) # 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) + seed_hessian_chunk!(xdual, x, indices, 1, iseeds, oseeds) + seed_hessian_chunk!(xdual, x, indices, 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) - nblocks > 1 && seed_hessian_chunk!(xdual, x, 1, nothing, nothing) + # the structural zeros of `x` are not variables, so no block writes their rows and columns + if xlen != length(x) + fill!(H, zero(eltype(H))) + if grad !== nothing + fill!(grad, zero(eltype(grad))) + end + end + extract_hessian_chunk!(T, H, ydual1, indices, 0, 0, N, N) + extract_hessian_gradient_chunk!(T, grad, ydual1, indices, 1, N) + nblocks > 1 && seed_hessian_chunk!(xdual, x, indices, 1, nothing, nothing) for q in 2:nblocks qoffset = (q - 1) * N @@ -121,20 +142,20 @@ function symmetric_hessian_expr(result_definition::Expr) # 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. # q's inner seeds remain unchanged throughout this loop. - seed_hessian_chunk!(xdual, x, qoffset + 1, iseeds, nothing, qsize) + 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, poffset + 1, nothing, oseeds) + seed_hessian_chunk!(xdual, x, indices, poffset + 1, nothing, oseeds) ydual = f(xdual) - extract_hessian_chunk!(T, H, ydual, poffset, qoffset, N, qsize) - seed_hessian_chunk!(xdual, x, poffset + 1, nothing, nothing) + 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 outer seeds while retaining its inner seeds. - seed_hessian_chunk!(xdual, x, qoffset + 1, iseeds, oseeds, qsize) + seed_hessian_chunk!(xdual, x, indices, 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) + extract_hessian_chunk!(T, H, ydual, indices, qoffset, qoffset, qsize, qsize) + extract_hessian_gradient_chunk!(T, grad, ydual, indices, qoffset + 1, qsize) + seed_hessian_chunk!(xdual, x, indices, qoffset + 1, nothing, nothing, qsize) end return H, ydual1 @@ -142,7 +163,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, valtype(T, valtype(T, typeof(ydual1))), xlen, xlen)))) + $(symmetric_hessian_expr(:(H = similar(x, valtype(T, valtype(T, typeof(ydual1))), length(x), length(x))))) 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 d773df20..00dea4d7 100644 --- a/test/AllocationsTest.jl +++ b/test/AllocationsTest.jl @@ -33,20 +33,22 @@ convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,F hduals = hcfg.duals iseeds = hcfg.iseeds oseeds = hcfg.oseeds + hindices = ForwardDiff.structural_linearindices(hduals, x) allocs_hseed!(args...) = @allocated ForwardDiff.seed_hessian_chunk!(args...) for i in (iseeds, nothing), o in (oseeds, nothing) - allocs_hseed!(hduals, x, 1, i, o) - @test iszero(allocs_hseed!(hduals, x, 1, i, o)) + allocs_hseed!(hduals, x, hindices, 1, i, o) + @test iszero(allocs_hseed!(hduals, x, hindices, 1, i, o)) end - allocs_hseed!(hduals, x, 1, nothing, nothing, 4) - @test iszero(allocs_hseed!(hduals, x, 1, nothing, nothing, 4)) + allocs_hseed!(hduals, x, hindices, 1, nothing, nothing, 4) + @test iszero(allocs_hseed!(hduals, x, hindices, 1, nothing, nothing, 4)) # a zero is free for an isbits value type, so only a `BigFloat` catches one being built for a # layer that was given seeds bx = BigFloat.(x) bcfg = ForwardDiff.HessianConfig(nothing, bx, ForwardDiff.Chunk{3}()) - allocs_hseed!(bcfg.duals, bx, 1, bcfg.iseeds, bcfg.oseeds) - @test iszero(allocs_hseed!(bcfg.duals, bx, 1, bcfg.iseeds, bcfg.oseeds)) + bindices = ForwardDiff.structural_linearindices(bcfg.duals, bx) + allocs_hseed!(bcfg.duals, bx, bindices, 1, bcfg.iseeds, bcfg.oseeds) + @test iszero(allocs_hseed!(bcfg.duals, bx, bindices, 1, bcfg.iseeds, bcfg.oseeds)) allocs_convert_test_574() = @allocated convert_test_574() allocs_convert_test_574() diff --git a/test/HessianTest.jl b/test/HessianTest.jl index bb068d9a..df74c068 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -217,25 +217,47 @@ end end 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) +@testset "$(nameof(W)), n = $n" for n in (3, 5), (W, sidx) in ( + (LowerTriangular, [i + n * (j - 1) for j in 1:n for i in j:n]), + (UpperTriangular, [i + n * (j - 1) for j in 1:n for i in 1:j]), + (Diagonal, 1:(n + 1):n^2), +) + x = W(randn(n, n)) + # d²f/dx[a]dx[b] is `1 + (a == b)` on the structural entries and zero everywhere else + f = z -> (sum(abs2, z) + sum(z)^2) / 2 + L = length(x) + + expected = zeros(L, L) + expected[sidx, sidx] .= 1 + for k in sidx + expected[k, k] += 1 + end + grad = zeros(n, n) + grad[sidx] .= x[sidx] .+ sum(x) + + # `length(sidx) - 1` makes the final chunk a partial one + @testset "chunk size = $c" for c in unique((1, 2, length(sidx) - 1, length(sidx))) + cfg = ForwardDiff.HessianConfig(f, x, ForwardDiff.Chunk{c}()) + + H = ForwardDiff.hessian(f, x, cfg) + @test H isa Matrix + @test size(H) == (L, L) @test H == expected - out = fill(NaN, xlen, xlen) - ForwardDiff.hessian!(out, objective, x) + out = fill(NaN, L, L) + @test ForwardDiff.hessian!(out, f, x, cfg) === out @test out == expected - flat = fill(NaN, xlen^2) - ForwardDiff.hessian!(flat, objective, x) - @test reshape(flat, xlen, xlen) == expected + flat = fill(NaN, L^2) + @test ForwardDiff.hessian!(flat, f, x, cfg) === flat + @test reshape(flat, L, L) == expected + + # `DiffResults.HessianResult` allocates a dense gradient buffer even for a structured `x` + result = DiffResults.HessianResult(x) + result = ForwardDiff.hessian!(result, f, x, cfg) + @test DiffResults.value(result) ≈ f(x) + @test DiffResults.gradient(result) == grad + @test DiffResults.hessian(result) == expected end end diff --git a/test/SeedTest.jl b/test/SeedTest.jl index d7829e6c..6a6454ca 100644 --- a/test/SeedTest.jl +++ b/test/SeedTest.jl @@ -16,13 +16,19 @@ include("utils.jl") # The expected structural index sets are written out by hand rather than obtained from # `structural_eachindex`, so a bug in that iterator cannot hide inside the assertions depending on # it; one test ties the two together. Order is significant: `index` and `count` are positions along -# the sequence, not array indices. The sets are heterogeneous by design — `Vector` and `Diagonal` -# enumerate linear indices (the latter via `diagind`), `UpperTriangular` enumerates `CartesianIndex` -# in column-major order. +# the sequence, not array indices. The `sidx` sets are heterogeneous by design — `Vector` and +# `Diagonal` enumerate linear indices (the latter via `diagind`), `UpperTriangular` enumerates +# `CartesianIndex` in column-major order — while `lidx` is the same positions as linear indices. const SEED_CASES = ( - (rand(10), collect(1:10)), - (UpperTriangular(rand(5, 5)), [CartesianIndex(i, j) for j in 1:5 for i in 1:j]), - (Diagonal(rand(6, 6)), collect(1:7:36)), + (rand(10), + 1:10, + 1:10), + (UpperTriangular(rand(5, 5)), + [CartesianIndex(i, j) for j in 1:5 for i in 1:j], + [i + 5 * (j - 1) for j in 1:5 for i in 1:j]), + (Diagonal(rand(6, 6)), + 1:7:36, + 1:7:36), ) # Positions within `sidx` whose partials are zero. @@ -42,7 +48,7 @@ function fill_marker!(duals, x, sidx, marker) return duals end -@testset "seed_zero_partials!: $(nameof(typeof(x)))" for (x, sidx) in SEED_CASES +@testset "seed_zero_partials!: $(nameof(typeof(x)))" for (x, sidx, _) in SEED_CASES cfg = ForwardDiff.GradientConfig(nothing, x, ForwardDiff.Chunk{3}()) duals, seeds = cfg.duals, cfg.seeds N = ForwardDiff.npartials(eltype(duals)) @@ -90,17 +96,21 @@ end end end -@testset "seed_hessian_chunk!: $(nameof(typeof(x)))" for (x, sidx) in SEED_CASES +@testset "seed_hessian_chunk!: $(nameof(typeof(x)))" for (x, sidx, lidx) in SEED_CASES (; duals, iseeds, oseeds) = ForwardDiff.HessianConfig(nothing, x, ForwardDiff.Chunk{3}()) nstruct = length(sidx) + indices = ForwardDiff.structural_linearindices(duals, x) - ForwardDiff.seed_hessian_chunk!(duals, x, 1, nothing, nothing, nstruct) - ForwardDiff.seed_hessian_chunk!(duals, x, 4, iseeds, oseeds) + # the windows below are positions along `indices`, so pin it to the implementation once + @test indices == lidx + + ForwardDiff.seed_hessian_chunk!(duals, x, indices, 1, nothing, nothing, nstruct) + ForwardDiff.seed_hessian_chunk!(duals, x, indices, 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) + ForwardDiff.seed_hessian_chunk!(duals, x, indices, 4, nothing, nothing) @test all(idx -> iszero(ForwardDiff.partials(ForwardDiff.value(duals[idx]))), sidx) @test all(idx -> iszero(ForwardDiff.partials(duals[idx])), sidx) end From 90d949659c46c12114efb9ad7dc4460ca4f482e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Mon, 31 Aug 2026 14:21:20 +0200 Subject: [PATCH 07/13] 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. The non-matrix method checks the length itself rather than leaving it to `reshape`, whose message names neither the Hessian nor the input. Co-Authored-By: Claude Opus 5 (1M context) --- ext/ForwardDiffStaticArraysExt.jl | 2 +- src/apiutils.jl | 4 ++-- src/hessian.jl | 31 +++++++++++++++++++++++-------- test/HessianTest.jl | 10 ++++++++++ 4 files changed, 36 insertions(+), 11 deletions(-) diff --git a/ext/ForwardDiffStaticArraysExt.jl b/ext/ForwardDiffStaticArraysExt.jl index ac7f8215..6eed7f73 100644 --- a/ext/ForwardDiffStaticArraysExt.jl +++ b/ext/ForwardDiffStaticArraysExt.jl @@ -133,7 +133,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_linearindices(x), 0, 0, length(x), length(x)) return result end diff --git a/src/apiutils.jl b/src/apiutils.jl index 2039d8f3..19ce48dd 100644 --- a/src/apiutils.jl +++ b/src/apiutils.jl @@ -70,8 +70,6 @@ function structural_eachindex(x::Diagonal, y::AbstractArray) return diagind(x) end -# The positions of `structural_eachindex`, in the same order, as linear indices of `x`. The two -# argument form is only ever given a config's work buffer and the input it is used with. function check_structural_size(duals, x) if size(duals) != size(x) throw(DimensionMismatch(lazy"the config was built for an array of size $(size(duals)) and cannot be used with an array of size $(size(x))")) @@ -79,6 +77,8 @@ function check_structural_size(duals, x) return nothing end +# The positions of `structural_eachindex`, in the same order, as linear indices of `x`. The two +# argument form is only ever given a config's work buffer and the input it is used with. structural_linearindices(x::AbstractArray) = structural_linearindices(x, x) function structural_linearindices(duals::AbstractArray, x::AbstractArray) require_one_based_indexing(duals, x) diff --git a/src/hessian.jl b/src/hessian.jl index 06d363a9..bb38d110 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -36,9 +36,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) - 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 @@ -55,10 +53,8 @@ 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) - 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 @@ -69,6 +65,23 @@ end const HESSIAN_ERROR = DimensionMismatch("hessian(f, x) expects that f(x) is a real number. Perhaps you meant jacobian(f, x)?") +# Mirrors `reshape_jacobian`. The sweep writes the result entry by entry, so nothing else checks it. +function reshape_hessian(result::AbstractMatrix, x) + require_one_based_indexing(result) + if size(result) != (length(x), length(x)) + throw(DimensionMismatch(lazy"cannot store the $(length(x))×$(length(x)) Hessian in a result of size $(size(result))")) + end + return result +end +function reshape_hessian(result::AbstractArray, x) + require_one_based_indexing(result) + if length(result) != length(x)^2 + throw(DimensionMismatch(lazy"cannot store the $(length(x))×$(length(x)) Hessian in a result of length $(length(result))")) + end + return reshape(result, length(x), length(x)) +end +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. `indices` maps a block # position to its row and column, both being linear indices of `x`. @@ -134,7 +147,9 @@ function symmetric_hessian_expr(result_definition::Expr) end extract_hessian_chunk!(T, H, ydual1, indices, 0, 0, N, N) extract_hessian_gradient_chunk!(T, grad, ydual1, indices, 1, N) - nblocks > 1 && seed_hessian_chunk!(xdual, x, indices, 1, nothing, nothing) + if nblocks > 1 + seed_hessian_chunk!(xdual, x, indices, 1, nothing, nothing) + end for q in 2:nblocks qoffset = (q - 1) * N diff --git a/test/HessianTest.jl b/test/HessianTest.jl index df74c068..cd6ebd5a 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -194,6 +194,16 @@ end @test reshape(flat, 3, 3) == ForwardDiff.hessian(prod, sx) end +@testset "a result of the wrong shape" begin + @testset "$(nameof(typeof(x)))" for x in (randn(3), SVector(1.0, 2.0, 3.0)) + @test_throws "cannot store the 3×3 Hessian in a result of size (4, 4)" ForwardDiff.hessian!(fill(NaN, 4, 4), sum, x) + @test_throws "cannot store the 3×3 Hessian in a result of size (2, 2)" ForwardDiff.hessian!(fill(NaN, 2, 2), sum, x) + @test_throws "cannot store the 3×3 Hessian in a result of length 8" ForwardDiff.hessian!(fill(NaN, 8), sum, x) + end + result = DiffResults.DiffResult(0.0, randn(3), fill(NaN, 4, 4)) + @test_throws "cannot store the 3×3 Hessian in a result of size (4, 4)" ForwardDiff.hessian!(result, sum, randn(3)) +end + @testset "an array-valued f is not a Hessian" begin sx = SVector(1.0, 2.0, 3.0) @test_throws DimensionMismatch ForwardDiff.hessian(identity, sx) From 8070750fb7e868577054ed9d8099d8cbe71e58bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Mon, 31 Aug 2026 17:56:09 +0200 Subject: [PATCH 08/13] Cover the two error paths Codecov flagged, and pin the messages `check_structural_size` and the `_unsetindex!` fallback were the only lines of the sweep left untested: the first needs a `HessianConfig` reused across sizes, the second a work buffer that is not an `Array`, so the existing `BigFloat` test with a `Vector` takes the `Base._unsetindex!` branch instead. Every `@test_throws` the sweep added asserted only the exception type. `HESSIAN_ERROR` is new here, but master threw `DimensionMismatch` at the same calls with the inner gradient's message, so those assertions passed on master too. They now match the message, with the type prefix so the type stays pinned. The chunk size guard named `ForwardDiff.structural_length`, which is internal. --- src/hessian.jl | 2 +- test/GradientTest.jl | 8 ++++++++ test/HessianTest.jl | 27 +++++++++++++++------------ 3 files changed, 24 insertions(+), 13 deletions(-) diff --git a/src/hessian.jl b/src/hessian.jl index bb38d110..73cc5c9f 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -119,7 +119,7 @@ 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)))")) + throw(ArgumentError(lazy"chunk size cannot be greater than the number of differentiated entries of x ($(N) > $(xlen))")) end # `N == 0` only for empty inputs, which still need one evaluation to determine the diff --git a/test/GradientTest.jl b/test/GradientTest.jl index c9967812..50f96656 100644 --- a/test/GradientTest.jl +++ b/test/GradientTest.jl @@ -280,6 +280,14 @@ end end end +# a structured work buffer cannot preserve an unassigned entry: it can only be left unassigned in +# an `Array` +@testset "$(nameof(T)) with an unassigned entry" for T in (LowerTriangular, UpperTriangular) + x = T(Matrix{BigFloat}(undef, 3, 3)) + @test_throws "ArgumentError: cannot differentiate at an input with an unassigned entry at index CartesianIndex(1, 1): that would leave an entry of the $(nameof(T)) work buffer unassigned" ForwardDiff.gradient(sum, x) +end +@test_throws "ArgumentError: cannot differentiate at an input with an unassigned entry at index 1: that would leave an entry of the Diagonal work buffer unassigned" ForwardDiff.gradient(sum, Diagonal(Vector{BigFloat}(undef, 3))) + # issue #769 @testset "functions with `Dual` output" begin x = [Dual{OuterTestTag}(Dual{TestTag}(1.3, 2.1), Dual{TestTag}(0.3, -2.4))] diff --git a/test/HessianTest.jl b/test/HessianTest.jl index cd6ebd5a..a7bef929 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -23,6 +23,8 @@ h = [-66.0 -40.0 0.0; -40.0 130.0 -80.0; 0.0 -80.0 200.0] +hessian_error = "DimensionMismatch: hessian(f, x) expects that f(x) is a real number. Perhaps you meant jacobian(f, x)?" + @testset "running hardcoded test with chunk size = $c and tag = $(repr(tag))" for c in HESSIAN_CHUNK_SIZES, tag in (nothing, Tag((f,ForwardDiff.gradient), eltype(x))) cfg = ForwardDiff.HessianConfig(f, x, ForwardDiff.Chunk{c}(), tag) resultcfg = ForwardDiff.HessianConfig(f, DiffResults.HessianResult(x), x, ForwardDiff.Chunk{c}(), tag) @@ -63,9 +65,10 @@ 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) +@test_throws "ArgumentError: chunk size cannot be greater than the number of differentiated entries of x (4 > 3)" ForwardDiff.hessian(f, x, ForwardDiff.HessianConfig(f, x, ForwardDiff.Chunk{length(x) + 1}())) +@test_throws hessian_error ForwardDiff.hessian(identity, x) +@test_throws hessian_error ForwardDiff.hessian!(similar(x, 3, 3), identity, x) +@test_throws "DimensionMismatch: the config was built for an array of size (3,) and cannot be used with an array of size (4,)" ForwardDiff.hessian(f, rand(4), ForwardDiff.HessianConfig(f, x)) ######################## @@ -123,7 +126,7 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) @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) + @test_throws hessian_error ForwardDiff.hessian(identity, sx) out = similar(x, 9, 9) ForwardDiff.hessian!(out, prod, sx) @@ -196,20 +199,20 @@ end @testset "a result of the wrong shape" begin @testset "$(nameof(typeof(x)))" for x in (randn(3), SVector(1.0, 2.0, 3.0)) - @test_throws "cannot store the 3×3 Hessian in a result of size (4, 4)" ForwardDiff.hessian!(fill(NaN, 4, 4), sum, x) - @test_throws "cannot store the 3×3 Hessian in a result of size (2, 2)" ForwardDiff.hessian!(fill(NaN, 2, 2), sum, x) - @test_throws "cannot store the 3×3 Hessian in a result of length 8" ForwardDiff.hessian!(fill(NaN, 8), sum, x) + @test_throws "DimensionMismatch: cannot store the 3×3 Hessian in a result of size (4, 4)" ForwardDiff.hessian!(fill(NaN, 4, 4), sum, x) + @test_throws "DimensionMismatch: cannot store the 3×3 Hessian in a result of size (2, 2)" ForwardDiff.hessian!(fill(NaN, 2, 2), sum, x) + @test_throws "DimensionMismatch: cannot store the 3×3 Hessian in a result of length 8" ForwardDiff.hessian!(fill(NaN, 8), sum, x) end result = DiffResults.DiffResult(0.0, randn(3), fill(NaN, 4, 4)) - @test_throws "cannot store the 3×3 Hessian in a result of size (4, 4)" ForwardDiff.hessian!(result, sum, randn(3)) + @test_throws "DimensionMismatch: cannot store the 3×3 Hessian in a result of size (4, 4)" ForwardDiff.hessian!(result, sum, randn(3)) end @testset "an array-valued f is not a Hessian" begin sx = SVector(1.0, 2.0, 3.0) - @test_throws DimensionMismatch ForwardDiff.hessian(identity, sx) - @test_throws DimensionMismatch ForwardDiff.hessian!(fill(NaN, 3, 3), identity, sx) - @test_throws DimensionMismatch ForwardDiff.hessian!(DiffResults.HessianResult(sx), identity, sx) - @test_throws DimensionMismatch ForwardDiff.hessian(identity, [1.0, 2.0, 3.0]) + @test_throws hessian_error ForwardDiff.hessian(identity, sx) + @test_throws hessian_error ForwardDiff.hessian!(fill(NaN, 3, 3), identity, sx) + @test_throws hessian_error ForwardDiff.hessian!(DiffResults.HessianResult(sx), identity, sx) + @test_throws hessian_error ForwardDiff.hessian(identity, [1.0, 2.0, 3.0]) end @testset "the result does not depend on the chunk size" begin From 2f11c0b57e815142d754aa9046cbac403f713823 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Tue, 1 Sep 2026 00:50:10 +0200 Subject: [PATCH 09/13] Check the Hessian gradient buffer, not only the Hessian `reshape_hessian` restored the check the sweep lost for the Hessian buffer, but the gradient buffer of a `DiffResult` was still written unchecked. Master validated it indirectly: the `f!(y, x)` `jacobian!` called `require_one_based_indexing` on it, and `reshape_jacobian` compared the Hessian's row count against `length(ydual)`, which followed the buffer. hessian!(DiffResult(0.0, fill(NaN, 4), fill(NaN, 3, 3)), f, rand(3)) # was: no error, gradient == [g1, g2, g3, NaN] hessian!(DiffResult(0.0, fill(NaN, 2), fill(NaN, 3, 3)), f, rand(3)) # was: BoundsError naming neither the gradient nor the input hessian!(DiffResult(0.0, fill(NaN, 4, 3), fill(NaN, 9, 9)), f, UpperTriangular(randn(3, 3))) # was: no error, the six derivatives written to linear positions 1,4,5,7,8,9 of a 4x3 buffer The last one is why this matters: silently wrong values in a buffer the caller believes was filled, and reachable only because indexing the Hessian by the linear indices of `x` made the `DiffResult` path work for structured inputs at all. `structural_eachindex(grad, x)` is the check the seeding utilities already apply to a work buffer, so no new helper is needed and `require_one_based_indexing` comes with it. It compares linear indices where both arrays are `IndexLinear`, which keeps a flat gradient buffer for a matrix `x` working as it did on master -- `length` is the requirement, not `size`, since the sweep writes by linear index of `x`. The new `@test_throws` assert the type alone: the message is Base's, not ours. Co-Authored-By: Claude Opus 5 (1M context) --- src/hessian.jl | 5 ++++- test/HessianTest.jl | 13 +++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/src/hessian.jl b/src/hessian.jl index 73cc5c9f..50c3846f 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -80,7 +80,10 @@ function reshape_hessian(result::AbstractArray, x) end return reshape(result, length(x), length(x)) end -reshape_hessian(result::DiffResult, x) = reshape_hessian(DiffResults.hessian(result), x) +function reshape_hessian(result::DiffResult, x) + structural_eachindex(DiffResults.gradient(result), x) + return reshape_hessian(DiffResults.hessian(result), x) +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. `indices` maps a block diff --git a/test/HessianTest.jl b/test/HessianTest.jl index a7bef929..10bd9594 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -205,6 +205,19 @@ end end result = DiffResults.DiffResult(0.0, randn(3), fill(NaN, 4, 4)) @test_throws "DimensionMismatch: cannot store the 3×3 Hessian in a result of size (4, 4)" ForwardDiff.hessian!(result, sum, randn(3)) + + # the gradient buffer is written by linear index of `x`, so it is checked against `x` as well + @test_throws DimensionMismatch ForwardDiff.hessian!(DiffResults.DiffResult(0.0, fill(NaN, 4), fill(NaN, 3, 3)), sum, randn(3)) + @test_throws DimensionMismatch ForwardDiff.hessian!(DiffResults.DiffResult(0.0, fill(NaN, 2), fill(NaN, 3, 3)), sum, randn(3)) + @test_throws DimensionMismatch ForwardDiff.hessian!(DiffResults.DiffResult(0.0, fill(NaN, 4, 3), fill(NaN, 9, 9)), sum, UpperTriangular(randn(3, 3))) +end + +@testset "a gradient buffer that is flat rather than shaped like x" begin + X = randn(3, 3) + f = z -> sum(abs2, z) + result = ForwardDiff.hessian!(DiffResults.DiffResult(0.0, fill(NaN, 9), fill(NaN, 9, 9)), f, X) + @test DiffResults.gradient(result) == vec(2 .* X) + @test DiffResults.hessian(result) == ForwardDiff.hessian(f, X) end @testset "an array-valued f is not a Hessian" begin From 2ef622c92a05d4563a63b26389a9d62616bad046 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Tue, 1 Sep 2026 00:52:43 +0200 Subject: [PATCH 10/13] Build the triangular position vector at its final size `structural_linearindices` returned a comprehension over a flattened generator for the two triangular wrappers, so the vector grew by `push!`: 512 bytes per sweep at n = 6, where the vector itself is 168. Writing into a `Vector{Int}` of the right length is one allocation of that length, and carrying the linear index along the columns -- advancing by `n - j` for the upper triangle and by `j` for the lower at the end of column j -- drops the multiplication as well. n = 6 512 -> 224 bytes n = 20 6832 -> 1840 bytes Dense inputs and `Diagonal` are unaffected, both still allocation-free. Co-Authored-By: Claude Opus 5 (1M context) --- src/apiutils.jl | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/src/apiutils.jl b/src/apiutils.jl index 19ce48dd..b640272a 100644 --- a/src/apiutils.jl +++ b/src/apiutils.jl @@ -89,13 +89,29 @@ function structural_linearindices(duals::UpperTriangular, x::AbstractArray) require_one_based_indexing(duals, x) check_structural_size(duals, x) n = size(duals, 1) - return [i + n * (j - 1) for j in 1:n for i in 1:j] + indices = Vector{Int}(undef, structural_length(duals)) + k = idx = 0 + for j in 1:n + for _ in 1:j + indices[k += 1] = (idx += 1) + end + idx += n - j + end + return indices end function structural_linearindices(duals::LowerTriangular, x::AbstractArray) require_one_based_indexing(duals, x) check_structural_size(duals, x) n = size(duals, 1) - return [i + n * (j - 1) for j in 1:n for i in j:n] + indices = Vector{Int}(undef, structural_length(duals)) + k = idx = 0 + for j in 1:n + for _ in j:n + indices[k += 1] = (idx += 1) + end + idx += j + end + return indices end function structural_linearindices(duals::Diagonal, x::AbstractArray) require_one_based_indexing(duals, x) From 40c2e0954b3bba6be30b8893e523024552c7f529 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Tue, 1 Sep 2026 00:55:14 +0200 Subject: [PATCH 11/13] Give the structured Hessian testset a reference that can fail `(sum(abs2, z) + sum(z)^2) / 2` has second derivative `1 + (a == b)` on the structural entries, so `expected` was uniform off the diagonal and any bijection of the structural positions maps the diagonal to the diagonal: the Hessian assertion was invariant under every permutation of them, and could not catch a wrong row or column in the off-diagonal blocks -- the one thing the structured testset exists to check. Only the gradient assertion pinned the order, and the gradient comes from the diagonal blocks alone. `dot(w, z)^2 / 2` has second derivative `w[a] * w[b]`, and `w * transpose(w)` for distinct `w` is reproduced by no permutation, since it would need `P * w == w`. With `w` the small integers on the structural positions and zero elsewhere the entries stay exactly representable, so the assertions remain `==` rather than `isapprox`, and `dot(w, x)` reads the same linear indices the sweep writes. Co-Authored-By: Claude Opus 5 (1M context) --- test/HessianTest.jl | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/test/HessianTest.jl b/test/HessianTest.jl index 10bd9594..596a0d8a 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -249,17 +249,14 @@ end (Diagonal, 1:(n + 1):n^2), ) x = W(randn(n, n)) - # d²f/dx[a]dx[b] is `1 + (a == b)` on the structural entries and zero everywhere else - f = z -> (sum(abs2, z) + sum(z)^2) / 2 L = length(x) + # d²f/dx[a]dx[b] is `w[a] * w[b]`, which no permutation of the structural positions reproduces + w = zeros(L) + w[sidx] = 1:length(sidx) + f = z -> dot(w, z)^2 / 2 - expected = zeros(L, L) - expected[sidx, sidx] .= 1 - for k in sidx - expected[k, k] += 1 - end - grad = zeros(n, n) - grad[sidx] .= x[sidx] .+ sum(x) + expected = w * transpose(w) + grad = reshape(w .* dot(w, x), n, n) # `length(sidx) - 1` makes the final chunk a partial one @testset "chunk size = $c" for c in unique((1, 2, length(sidx) - 1, length(sidx))) From 30f6c69f5ad505dbf8fd7a9db06992412a06dba6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Tue, 1 Sep 2026 00:55:41 +0200 Subject: [PATCH 12/13] Guard the allocation-free `hessian!`, and rename the seeding testset The PR's headline is that `hessian!` into a preallocated matrix allocates nothing, and nothing in the suite held it: `hessian_allocs()` from #720 covers the StaticArray `ImmutableDiffResult` path only. The sweep is allocation-free for a dense input because `structural_linearindices` returns a `Base.OneTo` and `structural_chunk` a `UnitRange` rather than a view; either of those turning back into an array would show up here and nowhere else. Both chunk shapes are covered, since a partial final block takes different branches: n = 40 with chunk 6 gives seven blocks with a four-wide last one, n = 10 with chunk 10 a single block. The seeding testset has covered `seed_hessian_chunk!` since the sweep landed, so its name no longer described it. Co-Authored-By: Claude Opus 5 (1M context) --- test/AllocationsTest.jl | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/test/AllocationsTest.jl b/test/AllocationsTest.jl index 00dea4d7..3a51885c 100644 --- a/test/AllocationsTest.jl +++ b/test/AllocationsTest.jl @@ -7,7 +7,7 @@ include(joinpath(dirname(@__FILE__), "utils.jl")) convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,Float64,8},4},2}, 1.3) -@testset "Test seed!/seed_zero_partials! allocations" begin +@testset "Test seeding allocations" begin x = rand(1000) cfg = ForwardDiff.GradientConfig(nothing, x) duals = cfg.duals @@ -71,6 +71,21 @@ end @test iszero(allocs_jacobian!()) end +@testset "Test hessian! allocations" begin + # the sweep is allocation-free only as long as the positions of a dense input stay a range: + # a `structural_linearindices` or `structural_chunk` returning an array would show up here + function allocs_hessian!(n, c) + f(z) = sum(abs2, z) + sum(z)^3 + x = randn(n) + result = zeros(n, n) + cfg = ForwardDiff.HessianConfig(f, x, ForwardDiff.Chunk{c}()) + ForwardDiff.hessian!(result, f, x, cfg) # warmup + return @allocated ForwardDiff.hessian!(result, f, x, cfg) + end + @test iszero(allocs_hessian!(40, 6)) # seven blocks, the last one partial + @test iszero(allocs_hessian!(10, 10)) # a single block +end + @testset "allocation-free nested StaticArray jacobian" begin # test that nested jacobians of StaticArrays do not allocate. # This is a regression test for issue #798, where the inner jacobian was allocating From b42dca24e0e83b42d3bbe2a702212d5817377c3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Tue, 1 Sep 2026 11:42:47 +0200 Subject: [PATCH 13/13] Give the Hessian's two dual layers distinct tags MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `HessianConfig` handed one `Tag(f, V)` to both layers of `Dual{T,Dual{T,V,N},N}`, leaving them type-indistinguishable. `value`/`partials` inside `f` could not tell them apart, the ordering machinery had nothing to compare, and the result depended on the code path and the chunk size (#845). The outer tag is now derived from the inner one as `outer_tag(T, Dual{T,V,N})`, and registered after it so that `T ≺ TO`. A result that carries only the inner layer then reads as a gradient with a vanishing Hessian, with no special casing, and the array path keeps the gradient master produced for such an `f`. Blocks whose derivatives the result cannot carry are no longer evaluated, and `extract_hessian` dispatches on the result rather than on its partials, which no longer matches a method when the result carries only an enclosing tag (#846). Zeros written for a derivative that vanishes take their element type from `ydual`, as `extract_gradient!` does, so both paths treat a plain buffer alike. Fixes #845. Fixes #846. Co-Authored-By: Claude Opus 5 (1M context) --- ext/ForwardDiffStaticArraysExt.jl | 43 ++++++----- src/apiutils.jl | 6 +- src/config.jl | 14 ++-- src/hessian.jl | 85 ++++++++++++++-------- test/HessianTest.jl | 116 ++++++++++++++++++++++++++++++ 5 files changed, 209 insertions(+), 55 deletions(-) diff --git a/ext/ForwardDiffStaticArraysExt.jl b/ext/ForwardDiffStaticArraysExt.jl index 6eed7f73..bafd266d 100644 --- a/ext/ForwardDiffStaticArraysExt.jl +++ b/ext/ForwardDiffStaticArraysExt.jl @@ -6,7 +6,7 @@ using ForwardDiff.DiffResults using ForwardDiff: Dual, partials, npartials, Partials, GradientConfig, JacobianConfig, HessianConfig, Tag, Chunk, gradient, hessian, jacobian, gradient!, hessian!, jacobian!, extract_gradient!, extract_jacobian!, extract_value!, structural_linearindices, - vector_mode_gradient, vector_mode_gradient!, + vector_mode_gradient, vector_mode_gradient!, outer_tag, vector_mode_jacobian, vector_mode_jacobian!, HESSIAN_ERROR, valtype, value using DiffResults: DiffResult, ImmutableDiffResult, MutableDiffResult @@ -107,34 +107,41 @@ end end # Hessian -@inline function extract_hessian(::Type{T}, ydual::Partials, x::StaticArray) where {T} - H = extract_jacobian(T, ydual, x) +@inline function extract_hessian(::Type{T}, ::Type{TO}, ydual::Dual{TO,<:Dual{T}}, x::StaticArray) where {T,TO} + H = extract_jacobian(T, partials(TO, ydual), x) 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))) +# A result that never picked up both perturbations has no second derivatives, and offers neither +# the `length(x)` rows the method above reads nor, for an `f` ignoring its argument, any at all. +@inline function extract_hessian(::Type{T}, ::Type{TO}, ydual, x::S) where {T,TO,S<:StaticArray} + R = StaticArrays.similar_type(S, valtype(T, valtype(TO, typeof(ydual))), + Size(length(x), length(x))) return zero(R) end -@inline function ForwardDiff.hessian(f::F, x::StaticArray) where {F} +# The layers need distinct tags; see `ForwardDiff.outer_tag`. +@inline function hessian_tags(f::F, x::StaticArray) where {F} T = typeof(Tag(f, eltype(x))) - ydual = f(dualize(T, dualize(T, x))) + return T, outer_tag(T, Dual{T,eltype(x),length(x)}) +end + +@inline function ForwardDiff.hessian(f::F, x::StaticArray) where {F} + T, TO = hessian_tags(f, x) + ydual = f(dualize(TO, dualize(T, x))) ydual isa Real || throw(HESSIAN_ERROR) - return extract_hessian(T, partials(T, ydual), x) + return extract_hessian(T, TO, 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) @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))) + T, TO = hessian_tags(f, x) + ydual = f(dualize(TO, dualize(T, x))) ydual isa Real || throw(HESSIAN_ERROR) H = ForwardDiff.reshape_hessian(result, x) - ForwardDiff.extract_hessian_chunk!(T, H, ydual, structural_linearindices(x), 0, 0, length(x), length(x)) + ForwardDiff.extract_hessian_chunk!(T, TO, H, ydual, structural_linearindices(x), 0, 0, length(x), length(x)) return result end @@ -144,14 +151,14 @@ ForwardDiff.hessian!(result::ImmutableDiffResult, f::F, x::StaticArray, cfg::Hes ForwardDiff.hessian!(result::ImmutableDiffResult, f::F, x::StaticArray, cfg::HessianConfig, ::Val) where {F} = hessian!(result, f, x) function ForwardDiff.hessian!(result::ImmutableDiffResult, f::F, x::StaticArray) where {F} - T = typeof(Tag(f, eltype(x))) + T, TO = hessian_tags(f, x) d1 = dualize(T, x) - d2 = dualize(T, d1) + d2 = dualize(TO, 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_hessian(T,partials(T,fd2), x) + val = value(T,value(TO,fd2)) + grad = extract_gradient(T,value(TO,fd2), x) + hess = extract_hessian(T,TO,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 b640272a..c7869c95 100644 --- a/src/apiutils.jl +++ b/src/apiutils.jl @@ -192,15 +192,15 @@ end # 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, +function seed_hessian_chunk!(duals::AbstractArray{Dual{TO,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} + chunksize = N) where {TO,T,V,N} izero = iseeds === nothing ? zero(Partials{N,V}) : nothing ozero = oseeds === nothing ? zero(Partials{N,Dual{T,V,N}}) : nothing idxs = structural_chunk(indices, index, 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]) + Dual{TO,Dual{T,V,N},N}(inner, oseeds === nothing ? ozero : oseeds[i]) end end diff --git a/src/config.jl b/src/config.jl index 5469b8c9..96dd73ee 100644 --- a/src/config.jl +++ b/src/config.jl @@ -195,12 +195,17 @@ Base.eltype(::Type{JacobianConfig{T,V,N,D}}) where {T,V,N,D} = Dual{T,V,N} # HessianConfig # ################# -struct HessianConfig{T,V,N,D} <: AbstractConfig{N} +struct HessianConfig{T,TO,V,N,D} <: AbstractConfig{N} iseeds::NTuple{N,Partials{N,V}} oseeds::NTuple{N,Partials{N,Dual{T,V,N}}} duals::D end +# The layers need distinct tags, or `value`/`partials` inside `f` cannot tell them apart (#845). +# `tagcount` fixes the ordering here rather than at the first comparison, as `Tag` does. +outer_tag(::Type{T}, ::Type{D}) where {T,D} = (tagcount(Tag{T,D}); Tag{T,D}) +outer_tag(::Type{Nothing}, ::Type) = Nothing + """ ForwardDiff.HessianConfig(f, x::AbstractArray, chunk::Chunk = Chunk(x)) @@ -224,8 +229,9 @@ function HessianConfig(f::F, ::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}) - return HessianConfig{T,V,N,typeof(duals)}(iseeds, oseeds, duals) + TO = outer_tag(T, Dual{T,V,N}) + duals = similar(x, Dual{TO,Dual{T,V,N},N}) + return HessianConfig{T,TO,V,N,typeof(duals)}(iseeds, oseeds, duals) end """ @@ -250,4 +256,4 @@ HessianConfig(f::F, 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,D}}) where {T,V,N,D} = Dual{T,Dual{T,V,N},N} +Base.eltype(::Type{HessianConfig{T,TO,V,N,D}}) where {T,TO,V,N,D} = Dual{TO,Dual{T,V,N},N} diff --git a/src/hessian.jl b/src/hessian.jl index 50c3846f..8bc403da 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -50,12 +50,12 @@ be used. 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} +function hessian!(result::DiffResult, f::F, x::AbstractArray, cfg::HessianConfig{T,TO} = HessianConfig(f, result, x), ::Val{CHK}=Val{true}()) where {F,T,TO,CHK} require_one_based_indexing(x) CHK && checktag(T, f, x) _, ydual = symmetric_hessian!(reshape_hessian(result, x), f, x, cfg, DiffResults.gradient(result)) - result = DiffResults.value!(result, value(T, value(T, ydual))) + result = DiffResults.value!(result, value(T, value(TO, ydual))) return result end @@ -88,11 +88,11 @@ 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. `indices` maps a block # position to its row and column, both being linear indices of `x`. -function extract_hessian_chunk!(::Type{T}, H, ydual, indices, roffset, coffset, rsize, csize) where {T} +function extract_hessian_chunk!(::Type{T}, ::Type{TO}, H, ydual::Dual{TO,<:Dual{T}}, indices, roffset, coffset, rsize, csize) where {T,TO} rows = structural_chunk(indices, roffset + 1, rsize) cols = structural_chunk(indices, coffset + 1, csize) for r in 1:rsize - drow = partials(T, ydual, r) + drow = partials(TO, ydual, r) i = rows[r] cstart = roffset == coffset ? r : 1 for c in cstart:csize @@ -105,11 +105,23 @@ function extract_hessian_chunk!(::Type{T}, H, ydual, indices, roffset, coffset, return H end +# Without both perturbations the block is zero. +function extract_hessian_chunk!(::Type{T}, ::Type{TO}, H, ydual, indices, roffset, coffset, rsize, csize) where {T,TO} + rows = structural_chunk(indices, roffset + 1, rsize) + cols = structural_chunk(indices, coffset + 1, csize) + h = zero(valtype(T, valtype(TO, typeof(ydual)))) + for j in cols, i in rows + H[i, j] = h + H[j, i] = h + end + return H +end + # The inner partials of a diagonal block contain the corresponding gradient chunk. # TODO: delegate to `extract_gradient_chunk!` once it takes its positions from `x` (#838). -extract_hessian_gradient_chunk!(::Type{T}, ::Nothing, ydual, indices, index, chunksize) where {T} = nothing -function extract_hessian_gradient_chunk!(::Type{T}, grad, ydual, indices, index, chunksize) where {T} - dual = value(T, ydual) +extract_hessian_gradient_chunk!(::Type{T}, ::Type{TO}, ::Nothing, ydual, indices, index, chunksize) where {T,TO} = nothing +function extract_hessian_gradient_chunk!(::Type{T}, ::Type{TO}, grad, ydual, indices, index, chunksize) where {T,TO} + dual = value(TO, ydual) for (i, idx) in enumerate(structural_chunk(indices, index, chunksize)) grad[idx] = partials(T, dual, i) end @@ -140,16 +152,27 @@ function symmetric_hessian_expr(result_definition::Expr) seed_hessian_chunk!(xdual, x, indices, N + 1, nothing, nothing, xlen - N) ydual1 = f(xdual) ydual1 isa Real || throw(HESSIAN_ERROR) + Vout = valtype(T, valtype(TO, typeof(ydual1))) $(result_definition) - # the structural zeros of `x` are not variables, so no block writes their rows and columns - if xlen != length(x) - fill!(H, zero(eltype(H))) - if grad !== nothing - fill!(grad, zero(eltype(grad))) - end + # A second derivative needs both perturbations, a first derivative only the inner one: + # what the result does not carry vanishes identically. + zero_hessian = !(ydual1 isa Dual{TO,<:Dual{T}}) + zero_gradient = zero_hessian && !(ydual1 isa Dual{T}) + + # Zero what no block writes: a derivative that vanishes, and the rows and columns of + # the structural zeros of `x`, which are not variables. + if zero_hessian || xlen != length(x) + fill!(H, zero(Vout)) + end + if grad !== nothing && (zero_gradient || xlen != length(x)) + fill!(grad, zero(Vout)) + end + # off-diagonal blocks find second derivatives, diagonal ones also the gradient + if zero_hessian && (zero_gradient || grad === nothing) + return H, ydual1 end - extract_hessian_chunk!(T, H, ydual1, indices, 0, 0, N, N) - extract_hessian_gradient_chunk!(T, grad, ydual1, indices, 1, N) + extract_hessian_chunk!(T, TO, H, ydual1, indices, 0, 0, N, N) + extract_hessian_gradient_chunk!(T, TO, grad, ydual1, indices, 1, N) if nblocks > 1 seed_hessian_chunk!(xdual, x, indices, 1, nothing, nothing) end @@ -157,22 +180,24 @@ function symmetric_hessian_expr(result_definition::Expr) for q in 2:nblocks qoffset = (q - 1) * N qsize = min(N, xlen - qoffset) - # 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. - # q's inner seeds remain unchanged throughout this loop. - 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, nothing, oseeds) - ydual = f(xdual) - extract_hessian_chunk!(T, H, ydual, indices, poffset, qoffset, N, qsize) - seed_hessian_chunk!(xdual, x, indices, poffset + 1, nothing, nothing) + if !zero_hessian + # 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. + # q's inner seeds remain unchanged throughout this loop. + 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, nothing, oseeds) + ydual = f(xdual) + extract_hessian_chunk!(T, TO, H, ydual, indices, poffset, qoffset, N, qsize) + seed_hessian_chunk!(xdual, x, indices, poffset + 1, nothing, nothing) + end end # 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) - extract_hessian_gradient_chunk!(T, grad, ydual, indices, qoffset + 1, qsize) + extract_hessian_chunk!(T, TO, H, ydual, indices, qoffset, qoffset, qsize, qsize) + extract_hessian_gradient_chunk!(T, TO, grad, ydual, indices, qoffset + 1, qsize) seed_hessian_chunk!(xdual, x, indices, qoffset + 1, nothing, nothing, qsize) end @@ -180,10 +205,10 @@ function symmetric_hessian_expr(result_definition::Expr) 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, valtype(T, valtype(T, typeof(ydual1))), length(x), length(x))))) +@eval function symmetric_hessian(f::F, x, cfg::HessianConfig{T,TO,V,N}, grad) where {F,T,TO,V,N} + $(symmetric_hessian_expr(:(H = similar(x, Vout, length(x), length(x))))) end -@eval function symmetric_hessian!(H, f::F, x, cfg::HessianConfig{T,V,N}, grad) where {F,T,V,N} +@eval function symmetric_hessian!(H, f::F, x, cfg::HessianConfig{T,TO,V,N}, grad) where {F,T,TO,V,N} $(symmetric_hessian_expr(:())) end diff --git a/test/HessianTest.jl b/test/HessianTest.jl index 596a0d8a..91db825e 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -243,6 +243,122 @@ end end end +# https://github.com/JuliaDiff/ForwardDiff.jl/issues/845 +@testset "f that inspects the layers of its argument" begin + # `ForwardDiff.value` drops the outer perturbation of an intermediate, not of the result + f = z -> sum(abs2, z) + ForwardDiff.value(z[1]) * z[2] + # the mixed derivative of `value(z[1]) * z[2]` is 1 in one order and 0 in the other, and + # only one triangle of block pairs is evaluated + expected_hessian = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 2.0] + + @testset "$(nameof(typeof(x))), chunk size = $c" for x in ( + [1.0, 2.0, 3.0], SVector(1.0, 2.0, 3.0), MVector(1.0, 2.0, 3.0), + ), c in HESSIAN_CHUNK_SIZES + cfg = ForwardDiff.HessianConfig(f, x, ForwardDiff.Chunk{c}()) + + @test ForwardDiff.hessian(f, x) == expected_hessian + @test ForwardDiff.hessian(f, x, cfg) == expected_hessian + + out = fill(NaN, 3, 3) + @test ForwardDiff.hessian!(out, f, x, cfg) === out + @test out == expected_hessian + + result = ForwardDiff.hessian!(DiffResults.HessianResult(x), f, x, cfg) + @test DiffResults.value(result) == 16.0 + @test DiffResults.gradient(result) == [4.0, 5.0, 6.0] + @test DiffResults.hessian(result) == expected_hessian + end +end + +# https://github.com/JuliaDiff/ForwardDiff.jl/issues/845 +# https://github.com/JuliaDiff/ForwardDiff.jl/issues/846 +@testset "a result that does not carry both perturbations" begin + xs = ([1.0, 2.0, 3.0], SVector(1.0, 2.0, 3.0), MVector(1.0, 2.0, 3.0)) + + @testset "the inner perturbation only: $(nameof(typeof(x)))" for x in xs + # the result itself lost its outer layer, so no second derivative survives -- but + # the gradient does, in the inner one + f = z -> ForwardDiff.value(sum(abs2, z)) + @test all(iszero, ForwardDiff.hessian(f, x)) + @test all(iszero, ForwardDiff.hessian!(fill(NaN, 3, 3), f, x)) + result = ForwardDiff.hessian!(DiffResults.HessianResult(x), f, x) + @test DiffResults.value(result) == 14.0 + @test DiffResults.gradient(result) == [2.0, 4.0, 6.0] + @test all(iszero, DiffResults.hessian(result)) + end + + @testset "neither perturbation: $(nameof(typeof(x)))" for x in xs + f = Returns(2.0) + @test all(iszero, ForwardDiff.hessian(f, x)) + @test all(iszero, ForwardDiff.hessian!(fill(NaN, 3, 3), f, x)) + result = ForwardDiff.hessian!(DiffResults.HessianResult(x), f, x) + @test DiffResults.value(result) == 2.0 + @test all(iszero, DiffResults.gradient(result)) + @test all(iszero, DiffResults.hessian(result)) + end + + @testset "an enclosing tag only: $(nameof(typeof(x)))" for x in xs + # `f` does not depend on `z`, so its result carries the `derivative` tag alone + @test ForwardDiff.derivative(a -> ForwardDiff.hessian(z -> a * 2.0, x)[1, 1], 1.0) == 0.0 + + ForwardDiff.derivative(1.0) do a + # a buffer that can hold the enclosing tag is written in full, one that cannot errors + H = fill(a * 111.0, 3, 3) + ForwardDiff.hessian!(H, z -> a * 2.0, x) + @test all(iszero, H) + @test_throws MethodError ForwardDiff.hessian!(fill(111.0, 3, 3), z -> a * 2.0, x) + return zero(a) + end + end +end + +@testset "no block is evaluated for a derivative the result cannot carry" begin + x = [1.0, 2.0, 3.0] + chunk = ForwardDiff.Chunk{1}() + evaluations = Ref(0) + + function second_order(z) + evaluations[] += 1 + return sum(abs2, z) + end + function first_order(z) + evaluations[] += 1 + return ForwardDiff.value(sum(abs2, z)) + end + function constant(z) + evaluations[] += 1 + return 2.0 + end + + # one evaluation per diagonal block and one per pair of distinct blocks + evaluations[] = 0 + ForwardDiff.hessian(second_order, x, ForwardDiff.HessianConfig(second_order, x, chunk)) + @test evaluations[] == 6 + + # without the outer perturbation only the diagonal blocks contribute, and only a gradient + evaluations[] = 0 + ForwardDiff.hessian(first_order, x, ForwardDiff.HessianConfig(first_order, x, chunk)) + @test evaluations[] == 1 + + evaluations[] = 0 + ForwardDiff.hessian!(DiffResults.HessianResult(x), first_order, x, + ForwardDiff.HessianConfig(first_order, x, chunk)) + @test evaluations[] == 3 + + # a constant still needs the evaluation that determines the output type + evaluations[] = 0 + ForwardDiff.hessian!(DiffResults.HessianResult(x), constant, x, + ForwardDiff.HessianConfig(constant, x, chunk)) + @test evaluations[] == 1 +end + +@testset "nested differentiation" begin + f = z -> sum(w -> w^3, z) + @testset "$(nameof(typeof(x)))" for x in ([1.0, 2.0, 3.0], SVector(1.0, 2.0, 3.0)) + @test ForwardDiff.derivative(a -> ForwardDiff.hessian(f, a .* x)[1, 1], 1.0) == 6.0 + end +end + @testset "$(nameof(W)), n = $n" for n in (3, 5), (W, sidx) in ( (LowerTriangular, [i + n * (j - 1) for j in 1:n for i in j:n]), (UpperTriangular, [i + n * (j - 1) for j in 1:n for i in 1:j]),