Change to centered difference

Author: ElOceanografo

src/differentiation/compute_hessian_ad.jl

@@ -6,25 +6,25 @@ struct ForwardColorHesCache{THS, THC, TI<:Integer, TD, TGF, TGC, TG}
     buffer::TD
     grad!::TGF
     grad_config::TGC
-    G::TG
-    dG::TG
+    G1::TG
+    G2::TG
 end
 
 function make_hessian_buffers(colorvec, x)
     ncolors = maximum(colorvec)
     D = hcat([float.(i .== colorvec) for i in 1:ncolors]...)
     buffer = similar(D)
-    G = similar(x)
-    dG = similar(x)
-    return (;ncolors, D, buffer, G, dG)
+    G1 = similar(x)
+    G2 = similar(x)
+    return (;ncolors, D, buffer, G1, G2)
 end
 
 function ForwardColorHesCache(f, 
                               x::AbstractVector{<:Number}, 
                               colorvec::AbstractVector{<:Integer}=eachindex(x), 
                               sparsity::Union{AbstractMatrix, Nothing}=nothing,
                               g! = (G, x, grad_config) -> ForwardDiff.gradient!(G, f, x, grad_config))
-    ncolors, D, buffer, G, dG = make_hessian_buffers(colorvec, x)
+    ncolors, D, buffer, G, G2 = make_hessian_buffers(colorvec, x)
     grad_config = ForwardDiff.GradientConfig(f, x)
 
     # If user supplied their own gradient function, make sure it has the right
@@ -42,7 +42,7 @@ function ForwardColorHesCache(f,
     if sparsity === nothing
         sparsity = sparse(ones(length(x), length(x)))
     end
-    return ForwardColorHesCache(sparsity, colorvec, ncolors, D, buffer, g1!, grad_config, G, dG)
+    return ForwardColorHesCache(sparsity, colorvec, ncolors, D, buffer, g1!, grad_config, G, G2)
 end
 
 function numauto_color_hessian!(H::AbstractMatrix{<:Number}, 
@@ -52,11 +52,12 @@ function numauto_color_hessian!(H::AbstractMatrix{<:Number},
                                     safe = true)
     ϵ = cbrt(eps(eltype(x)))
     for j in 1:hes_cache.ncolors
-        hes_cache.grad!(hes_cache.G, x, hes_cache.grad_config)
         x .+= ϵ .* @view hes_cache.D[:, j]
-        hes_cache.grad!(hes_cache.dG, x, hes_cache.grad_config)
-        x .-= ϵ .* @view hes_cache.D[:, j]
-        hes_cache.buffer[:, j] .= (hes_cache.dG .- hes_cache.G) ./ ϵ
+        hes_cache.grad!(hes_cache.G2, x, hes_cache.grad_config)
+        x .-= 2ϵ .* @view hes_cache.D[:, j]
+        hes_cache.grad!(hes_cache.G1, x, hes_cache.grad_config)
+        hes_cache.buffer[:, j] .= (hes_cache.G2 .- hes_cache.G1) ./ 2ϵ
+        x .+= ϵ .* @view hes_cache.D[:, j] #reset to original value
     end
     ii, jj, vv = findnz(hes_cache.sparsity)
     if safe

test/test_sparse_hessian.jl

@@ -13,26 +13,26 @@ colors = matrix_colors(tril(sparsity))
 ncolors = maximum(colors)
 D = hcat([float.(i .== colors) for i in 1:ncolors]...)
 buffer = similar(D)
-G = zero(x)
-dG = zero(x)
+G1 = zero(x)
+G2 = zero(x)
 
 buffers_tup = SparseDiffTools.make_hessian_buffers(colors, x)
 @test buffers_tup.ncolors == ncolors
 @test buffers_tup.D == D
 @test size(buffers_tup.buffer) == size(buffer)
 @test eltype(buffers_tup.buffer) == eltype(buffer)
 @test typeof(buffers_tup.buffer) == typeof(buffer)
-@test size(buffers_tup.G) == size(G)
-@test eltype(buffers_tup.G) == eltype(G)
-@test size(buffers_tup.dG) == size(dG)
-@test eltype(buffers_tup.dG) == eltype(dG)
+@test size(buffers_tup.G1) == size(G1)
+@test eltype(buffers_tup.G1) == eltype(G1)
+@test size(buffers_tup.G2) == size(G2)
+@test eltype(buffers_tup.G2) == eltype(G2)
 
 
 gconfig = ForwardDiff.GradientConfig(fscalar, x)
 g(x) = ForwardDiff.gradient(fscalar, x)           # allocating
 g!(G, x, gconfig) = ForwardDiff.gradient!(G, fscalar, x, gconfig)   # non-allocating
 
-hescache1 = ForwardColorHesCache(sparsity, colors, ncolors, D, buffer, g!, gconfig, G, dG)
+hescache1 = ForwardColorHesCache(sparsity, colors, ncolors, D, buffer, g!, gconfig, G1, G2)
 hescache2 = ForwardColorHesCache(fscalar, x, colors, sparsity, g!)
 hescache3 = ForwardColorHesCache(fscalar, x, colors, sparsity)
 # custom gradient function
@@ -53,7 +53,7 @@ for name in [:sparsity, :colors, :ncolors, :D]
     # hescache5 is the default dense version, so only first axis will match
     @eval @test size(hescache1.$name, 1) == size(hescache5.$name, 1)
 end
-for name in [:buffer, :G, :dG]
+for name in [:buffer, :G1, :G2]
     @eval @test size(hescache1.$name) == size(hescache2.$name)
     @eval @test size(hescache1.$name) == size(hescache3.$name)
     @eval @test size(hescache1.$name) == size(hescache4.$name)
@@ -71,9 +71,9 @@ for (i, hescache) in enumerate([hescache1, hescache2, hescache3, hescache4, hesc
     H = numauto_color_hessian(fscalar, x, colors, sparsity)
     H1 = numauto_color_hessian(fscalar, x, hescache)
     H2 = numauto_color_hessian(fscalar, x)
-    @test all(isapprox.(Hforward, H, rtol=1e-4))
-    @test all(isapprox.(H, H1, rtol=1e-4))
-    @test all(isapprox.(H2, H1, rtol=1e-4))
+    @test all(isapprox.(Hforward, H, rtol=1e-6))
+    @test all(isapprox.(H, H1, rtol=1e-6))
+    @test all(isapprox.(H2, H1, rtol=1e-6))
 
     H1 = similar(H)
     numauto_color_hessian!(H1, fscalar, x, collect(hescache.colors), hescache.sparsity)