diff --git a/ext/ForwardDiffStaticArraysExt.jl b/ext/ForwardDiffStaticArraysExt.jl index 010af29b..1ef9dbc0 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!, structural_indices, 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,36 @@ 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 + +# 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) +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 = ForwardDiff.reshape_hessian(result, x) + ForwardDiff.extract_hessian_chunk!(T, H, ydual, structural_indices(x), 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 +148,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 daf9a1f6..3e749156 100644 --- a/src/apiutils.jl +++ b/src/apiutils.jl @@ -174,3 +174,17 @@ function seed!(duals::AbstractArray{Dual{T,V,N}}, x, indices, index, Dual{T,V,N}(value, seeds[i]) end 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, + 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 + return _seed!(duals, x, structural_chunk(indices, index, chunksize)) do value, i + inner = Dual{T,V,N}(value, iseeds === nothing ? izero : iseeds[i]) + Dual{T,Dual{T,V,N},N}(inner, oseeds === nothing ? ozero : oseeds[i]) + end +end diff --git a/src/config.jl b/src/config.jl index 3db333f6..6f167500 100644 --- a/src/config.jl +++ b/src/config.jl @@ -206,9 +206,11 @@ Base.eltype(::Type{JacobianConfig{T,V,N,D,I}}) where {T,V,N,D,I} = Dual{T,V,N} # HessianConfig # ################# -struct HessianConfig{T,V,N,DG,DJ,IG,IJ} <: AbstractConfig{N} - jacobian_config::JacobianConfig{T,V,N,DJ,IJ} - gradient_config::GradientConfig{T,Dual{T,V,N},N,DG,IG} +struct HessianConfig{T,V,N,D,I} <: AbstractConfig{N} + iseeds::NTuple{N,Partials{N,V}} + oseeds::NTuple{N,Partials{N,Dual{T,V,N}}} + duals::D + indices::I end """ @@ -218,10 +220,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 @@ -231,11 +232,13 @@ This constructor does not store/modify `x`. """ function HessianConfig(f::F, x::AbstractArray{V}, - chunk::Chunk = Chunk(x), - tag = Tag(f, V)) where {F,V} - jacobian_config = JacobianConfig(f, x, chunk, tag) - gradient_config = GradientConfig(f, jacobian_config.duals, chunk, tag) - return HessianConfig(jacobian_config, gradient_config) + ::Chunk{N} = Chunk(x), + ::T = Tag(f, V)) where {F,V,N,T} + iseeds = construct_seeds(Partials{N,V}) + oseeds = construct_seeds(Partials{N,Dual{T,V,N}}) + duals = similar(x, Dual{T,Dual{T,V,N},N}) + indices = structural_indices(duals) + return HessianConfig{T,V,N,typeof(duals),typeof(indices)}(iseeds, oseeds, duals, indices) end """ @@ -244,25 +247,20 @@ end Return a `HessianConfig` instance based on the type of `f`, types/storage in `result`, and type/shape of the input vector `x`. -The returned `HessianConfig` instance contains all the work buffers required by -`ForwardDiff.hessian!` for the case where the `result` argument is an `DiffResult`. +Equivalent to `ForwardDiff.HessianConfig(f, x, chunk)`: the work buffers do not depend on +`result`. The result-aware form is retained for compatibility. If `f` is `nothing` instead of the actual target function, then the returned instance can be used with any target function. However, this will reduce ForwardDiff's ability to catch and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83). -This constructor does not store/modify `x`. +This constructor does not store/modify `result` or `x`. """ -function HessianConfig(f::F, - result::DiffResult, - x::AbstractArray{V}, - chunk::Chunk = Chunk(x), - tag = Tag(f, V)) where {F,V} - jacobian_config = JacobianConfig((f,gradient), DiffResults.gradient(result), x, chunk, tag) - gradient_config = GradientConfig(f, jacobian_config.duals[2], chunk, tag) - return HessianConfig(jacobian_config, gradient_config) -end +HessianConfig(f::F, + ::DiffResult, + x::AbstractArray{V}, + chunk::Chunk = Chunk(x), + tag = Tag(f, V)) where {F,V} = HessianConfig(f, x, chunk, tag) checktag(::HessianConfig{T},f,x) where {T} = checktag(T,f,x) -Base.eltype(::Type{HessianConfig{T,V,N,DG,DJ,IG,IJ}}) where {T,V,N,DG,DJ,IG,IJ} = - Dual{T,Dual{T,V,N},N} +Base.eltype(::Type{HessianConfig{T,V,N,D,I}}) where {T,V,N,D,I} = Dual{T,Dual{T,V,N},N} diff --git a/src/hessian.jl b/src/hessian.jl index 9c755c9a..ddd59043 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -5,7 +5,12 @@ """ 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)`. +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. This method assumes that `isa(f(x), Real)`. @@ -14,15 +19,17 @@ 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}()) + checkstructure(cfg, x) + H, _ = symmetric_hessian(f, x, cfg, nothing) + return H 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)`. @@ -31,41 +38,131 @@ 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}()) + checkstructure(cfg, x) + symmetric_hessian!(reshape_hessian(result, x), 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}()) -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. """ 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 + checkstructure(cfg, x) + _, 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 + +############################ +# symmetric Hessian kernel # +############################ + +const HESSIAN_ERROR = DimensionMismatch("hessian(f, x) expects that f(x) is a real number. Perhaps you meant jacobian(f, x)?") + +function reshape_hessian(result::AbstractMatrix, x) + require_one_based_indexing(result) + size(result) == (length(x), length(x)) || throw(DimensionMismatch( + lazy"cannot store the $(length(x))×$(length(x)) Hessian in a result of size $(size(result))")) + return result +end +reshape_hessian(result::AbstractArray, x) = reshape(result, length(x), length(x)) +reshape_hessian(result::DiffResult, x) = reshape_hessian(DiffResults.hessian(result), x) + +# Copy a block from the nested partials and fill its transpose. On diagonal blocks, read +# only the upper triangle so the result is exactly symmetric. Both axes are indexed by the linear +# indices of `x`, as the columns of a Jacobian are, so `indices` gives the row and column of a block. +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] + for c in (roffset == coffset ? r : 1):csize + h = partials(T, drow, c) + 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, indices, index, chunksize) where {T} = nothing +extract_hessian_gradient_chunk!(::Type{T}, grad, ydual, indices, index, chunksize) where {T} = + extract_gradient_chunk!(T, grad, value(T, ydual), indices, 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 + + # `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.duals + indices = cfg.indices + iseeds = cfg.iseeds + oseeds = cfg.oseeds + + # The first evaluation determines the output type. Seeding the first block and clearing + # the untouched tail partitions the fresh buffer, so every element is initialized once. + 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) + # unwrapped once so that the helpers' `valtype` reaches the value type + zero_unseeded_columns!(T, H, value(T, ydual1), x) + grad === nothing || zero_unseeded!(T, grad, value(T, ydual1), x) + 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 + 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. + 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) + 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) + seed_hessian_chunk!(xdual, x, indices, 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, 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} + $(symmetric_hessian_expr(:())) end diff --git a/test/AllocationsTest.jl b/test/AllocationsTest.jl index 5e21ef62..3108039d 100644 --- a/test/AllocationsTest.jl +++ b/test/AllocationsTest.jl @@ -11,9 +11,7 @@ convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,F @testset "Test seed!/seed_zero_partials! allocations" begin x = rand(1000) cfg = ForwardDiff.GradientConfig(nothing, x) - duals = cfg.duals - seeds = cfg.seeds - indices = cfg.indices + (; duals, seeds, indices) = cfg allocs_seed!(args...) = @allocated ForwardDiff.seed!(args...) allocs_seed!(duals, x, indices, seeds) @@ -34,6 +32,33 @@ convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,F @test iszero(allocs_convert_test_574()) end +@testset "Test seed_hessian_chunk! allocations" begin + x = rand(1000) + cfg = ForwardDiff.HessianConfig(nothing, x) + (; duals, indices, iseeds, oseeds) = cfg + + allocs_hseed!(args...) = @allocated ForwardDiff.seed_hessian_chunk!(args...) + # all four seed combinations, the mixed ones being what the off-diagonal blocks use + @testset "iseeds=$(i === nothing) oseeds=$(o === nothing)" for (i, o) in + ((iseeds, oseeds), + (iseeds, nothing), + (nothing, oseeds), + (nothing, nothing)) + allocs_hseed!(duals, x, indices, 1, i, o) + @test iszero(allocs_hseed!(duals, x, indices, 1, i, o)) + allocs_hseed!(duals, x, indices, 1, i, o, 4) + @test iszero(allocs_hseed!(duals, x, indices, 1, i, o, 4)) + end + + # a zero of a non-isbits value type does allocate, so supplying both seeds must not build one + @testset "BigFloat" begin + y = BigFloat[1, 2, 3] + cfg = ForwardDiff.HessianConfig(nothing, y, ForwardDiff.Chunk{3}()) + allocs_hseed!(cfg.duals, y, cfg.indices, 1, cfg.iseeds, cfg.oseeds) + @test iszero(allocs_hseed!(cfg.duals, y, cfg.indices, 1, cfg.iseeds, cfg.oseeds)) + end +end + @testset "Test jacobian! allocations" begin # jacobian! should not allocate when called with a pre-allocated result Matrix. # Previously, reshape() inside extract_jacobian! allocated a wrapper diff --git a/test/GradientTest.jl b/test/GradientTest.jl index 9d7d75c2..230a6798 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 9c72faf3..47c1446b 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -51,11 +51,29 @@ 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) + +@testset "wrongly sized result: $(nameof(typeof(z)))" for z in ([1.0, 2.0, 3.0], + SVector(1.0, 2.0, 3.0)) + msg = "DimensionMismatch: cannot store the 3×3 Hessian in a result of size (4, 4)" + @test_throws msg ForwardDiff.hessian!(fill(NaN, 4, 4), prod, z) + @test_throws msg ForwardDiff.hessian!(DiffResults.DiffResult(0.0, zeros(3), + fill(NaN, 4, 4)), prod, z) +end ######################## @@ -108,10 +126,28 @@ 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) + @test_throws DimensionMismatch ForwardDiff.hessian!(similar(x, 9, 9), identity, sx) + @test_throws DimensionMismatch ForwardDiff.hessian!(DiffResults.HessianResult(sx), identity, sx) + + flat = fill(NaN, 81) + @test ForwardDiff.hessian!(flat, prod, sx) === flat + @test reshape(flat, 9, 9) == actual + out = similar(x, 9, 9) ForwardDiff.hessian!(out, prod, sx) @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,14 +192,33 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) @test DiffResults.hessian(sresult3) == DiffResults.hessian(result) end -# issues #838 and #839, which `hessian` inherits through `jacobian(gradient(f), x)` -@testset "structured inputs: $(nameof(W))" for (W, sidx) in ( +@testset "empty input: $(nameof(typeof(z)))" for z in (Float64[], SVector{0,Float64}()) + @test ForwardDiff.hessian(sum, z) == zeros(0, 0) + out = fill(NaN, 0, 0) + @test ForwardDiff.hessian!(out, sum, z) === out +end + +# `log(sum(exp, z))` rounds differently in the two nesting orders, hence the bitwise comparison +@testset "chunk size independence" begin + n = 16 + x = randn(n) + f = z -> log(sum(exp, z)) + expected = ForwardDiff.hessian(f, SVector{n}(x)) + @testset "chunk size = $c" for c in (1, 2, 3, 5, 7, 11, n) + cfg = ForwardDiff.HessianConfig(f, x, ForwardDiff.Chunk{c}()) + @test ForwardDiff.hessian(f, x, cfg) == expected + end +end + +# `n = 5` is the only structured case whose default `Chunk(x)` reaches the off-diagonal blocks. +@testset "structured inputs: $(nameof(W)) of size $n" for n in (3, 5), + (W, sidx) in ( # both axes are indexed by the linear indices of `x`, hard zeros off the structure - (LowerTriangular, [i + 3 * (j - 1) for j in 1:3 for i in j:3]), - (UpperTriangular, [i + 3 * (j - 1) for j in 1:3 for i in 1:j]), - (Diagonal, 1:4:9), + (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(3, 3)) + x = W(randn(n, n)) # d²f/dx[a]dx[b] is `1 + (a == b)` for structural `a`, `b`, and zero everywhere else f = z -> (sum(abs2, z) + sum(z)^2) / 2 L = length(x) @@ -174,14 +229,16 @@ end expected[k, k] += 1 end val = f(x) - grad = zeros(3, 3) + grad = zeros(n, n) grad[sidx] .= x[sidx] .+ sum(x) - # one chunk size below the full length, so that the final chunk is 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}()) + # `nothing` is the default chunk; `length(sidx) - 1` makes the final chunk a partial one + @testset "chunk size = $c" for c in (nothing, unique((1, 2, length(sidx) - 1, length(sidx)))...) + chunk = c === nothing ? ForwardDiff.Chunk(x) : ForwardDiff.Chunk{c}() + cfg = ForwardDiff.HessianConfig(f, x, chunk) H = ForwardDiff.hessian(f, x, cfg) + @test H isa Matrix @test size(H) == (L, L) @test H == expected @@ -189,10 +246,14 @@ end @test ForwardDiff.hessian!(out, f, x, cfg) === out @test out == 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, - ForwardDiff.HessianConfig(f, result, x, ForwardDiff.Chunk{c}())) + ForwardDiff.HessianConfig(f, result, x, chunk)) @test DiffResults.value(result) ≈ val @test DiffResults.gradient(result) == grad @test DiffResults.hessian(result) == expected @@ -222,6 +283,28 @@ end 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 7cee2b21..51ada2df 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -202,6 +202,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 7c4f1792..8fc55d39 100644 --- a/test/SeedTest.jl +++ b/test/SeedTest.jl @@ -132,4 +132,20 @@ 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, indices, iseeds, oseeds) = cfg + nstruct = length(sidx) + + 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, 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 + end # module