From 978826543df373ef1be4fd01da1468f669d1b7c4 Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Sun, 22 May 2022 14:58:02 -0700 Subject: [PATCH 01/19] start adding sparse Hessian stuff --- src/SparseDiffTools.jl | 3 ++ src/differentiation/compute_hessian.jl | 58 ++++++++++++++++++++++++++ test/test_sparse_hessian.jl | 36 ++++++++++++++++ 3 files changed, 97 insertions(+) create mode 100644 src/differentiation/compute_hessian.jl create mode 100644 test/test_sparse_hessian.jl diff --git a/src/SparseDiffTools.jl b/src/SparseDiffTools.jl index 07ca008f..1293f2ce 100644 --- a/src/SparseDiffTools.jl +++ b/src/SparseDiffTools.jl @@ -28,6 +28,9 @@ export contract_color, forwarddiff_color_jacobian!, forwarddiff_color_jacobian, ForwardColorJacCache, + forwarddiff_color_hessian!, + forwarddiff_color_hessian, + ForwardColorHesCache, auto_jacvec,auto_jacvec!, num_jacvec,num_jacvec!, num_hesvec,num_hesvec!, diff --git a/src/differentiation/compute_hessian.jl b/src/differentiation/compute_hessian.jl new file mode 100644 index 00000000..754aa54b --- /dev/null +++ b/src/differentiation/compute_hessian.jl @@ -0,0 +1,58 @@ + + + +# ForwardColorJacCache(f,x,_chunksize = nothing; +# dx = nothing, +# colorvec=1:length(x), +# sparsity = nothing) + + +# function HessianConfig(logdensity, imarginal, ijoint, forwarddiff_sparsity=false) +# if forwarddiff_sparsity +# println("Detecting Hessian sparsity via ForwardDiff...") +# H = ForwardDiff.hessian(f, x) +# Hsparsity = sparse(H)[imarginal, imarginal] .!= 0 +# else +# println("Detecting Hessian sparsity via SparsityDetection...") +# Hsparsity = hessian_sparsity(f, x)[imarginal, imarginal] +# end +# Hcolors = matrix_colors(Hsparsity) +# + +struct ForwardColorHesCache{THS, THC, TI<:Integer, TD, TG} + Hsparsity::THS + Hcolors::THC + ncolors::TI + D::TD + Hcomp_buffer::TD + G::TG + δG::TG +end + +function ForwardColorHesCache(f, x, colorvec=1:length(x), sparsity=nothing) + D = hcat([float.(i .== colorvec) for i in 1:maximum(colorvec)]...) + Hcomp_buffer = similar(D) + G = zero(x) + δG = zero(x) + return ForwardColorHesCache(Hsparsity, colorvec, size(Hcolors, 2), D, Hcomp_buffer, G, δG) +end + +function sparse_hessian!(f, g!, θ, hessconfig::ForwardColorHesCache, δ=sqrt(eps(Float64))) + nc = hessconfig.ncolors + for j in one(nc):nc + g!(hessconfig.G, θ) + g!(hessconfig.δG, θ + δ * @view hessconfig.D[:, j]) + hessconfig.Hcomp_buffer[:, j] .= (hessconfig.δG .- hessconfig.G) ./ δ + end + ii, jj, vv = findnz(hessconfig.Hsparsity) + H = sparse(ii, jj, zeros(length(vv))) + for (i, j) in zip(ii, jj) + H[i, j] = hessconfig.Hcomp_buffer[i, hessconfig.Hcolors[j]] + end + return H +end + +function sparse_hessian!(f, θ, hessconfig::ForwardColorHesCache, δ=sqrt(eps(Float64))) + g!(G, θ) = ForwardDiff.gradient!(G, f, θ) + return sparse_hessian!(f, g!, θ, hessconfig, δ) +end diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl new file mode 100644 index 00000000..992bce86 --- /dev/null +++ b/test/test_sparse_hessian.jl @@ -0,0 +1,36 @@ + +## Hessian tests +using SparsityDetection, SparseDiffTools +using LinearAlgebra, SparseArrays + +function fscalar(x) + return -dot(x, x) +end + +x = randn(100) +hs = hessian_sparsity(fscalar, x) +col = matrix_colors(hs) +hescache = ForwardColorHesCache(fscalar, x, col) + +# ForwardColorJacCache(f,x,_chunksize = nothing; +# dx = nothing, +# colorvec=1:length(x), +# sparsity = nothing) + +# forwarddiff_color_hessian!(J::AbstractMatrix{<:Number}, +# f, +# x::AbstractArray{<:Number}; +# dx = nothing, +# colorvec = eachindex(x), +# sparsity = nothing) + +# forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, +# f, +# x::AbstractArray{<:Number}, +# jac_cache::ForwardColorJacCache) + +# jacout = forwarddiff_color_jacobian(g, x, +# dx = similar(x), +# colorvec = 1:length(x), +# sparsity = nothing, +# jac_prototype = nothing) # matrix w/ sparsity pattern \ No newline at end of file From 94adc3fb9e53365ffb78f94cf7dab7a7843254bd Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Mon, 23 May 2022 22:10:29 -0700 Subject: [PATCH 02/19] basic in-place and mutating sparse hessians --- Project.toml | 5 +- src/SparseDiffTools.jl | 1 + src/differentiation/compute_hessian.jl | 58 ----------- src/differentiation/compute_hessian_ad.jl | 112 ++++++++++++++++++++++ 4 files changed, 115 insertions(+), 61 deletions(-) delete mode 100644 src/differentiation/compute_hessian.jl create mode 100644 src/differentiation/compute_hessian_ad.jl diff --git a/Project.toml b/Project.toml index 4c78eb55..660f4184 100644 --- a/Project.toml +++ b/Project.toml @@ -15,26 +15,25 @@ Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Requires = "ae029012-a4dd-5104-9daa-d747884805df" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SparsityDetection = "684fba80-ace3-11e9-3d08-3bc7ed6f96df" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" VertexSafeGraphs = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" [compat] Adapt = "1, 2.0, 3.0" ArrayInterfaceCore = "0.1.1" -ArrayInterfaceStaticArrays = "0.1" Compat = "2.2, 3" DataStructures = "0.17, 0.18" FiniteDiff = "2.8.1" ForwardDiff = "0.10" -Graphs = "1.4" Requires = "0.5, 1.0" StaticArrays = "1" VertexSafeGraphs = "0.2" julia = "1.6" [extras] -ArrayInterfaceBlockBandedMatrices = "5331f1e9-51c7-46b0-a9b0-df4434785e0a" ArrayInterfaceBandedMatrices = "2e50d22c-5be1-4042-81b1-c572ed69783d" +ArrayInterfaceBlockBandedMatrices = "5331f1e9-51c7-46b0-a9b0-df4434785e0a" BandedMatrices = "aae01518-5342-5314-be14-df237901396f" BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" diff --git a/src/SparseDiffTools.jl b/src/SparseDiffTools.jl index 7d992647..70807db2 100644 --- a/src/SparseDiffTools.jl +++ b/src/SparseDiffTools.jl @@ -51,6 +51,7 @@ include("coloring/greedy_star1_coloring.jl") include("coloring/greedy_star2_coloring.jl") include("coloring/matrix2graph.jl") include("differentiation/compute_jacobian_ad.jl") +include("differentiation/compute_hessian_ad.jl") include("differentiation/jaches_products.jl") include("differentiation/vecjac_products.jl") diff --git a/src/differentiation/compute_hessian.jl b/src/differentiation/compute_hessian.jl deleted file mode 100644 index 754aa54b..00000000 --- a/src/differentiation/compute_hessian.jl +++ /dev/null @@ -1,58 +0,0 @@ - - - -# ForwardColorJacCache(f,x,_chunksize = nothing; -# dx = nothing, -# colorvec=1:length(x), -# sparsity = nothing) - - -# function HessianConfig(logdensity, imarginal, ijoint, forwarddiff_sparsity=false) -# if forwarddiff_sparsity -# println("Detecting Hessian sparsity via ForwardDiff...") -# H = ForwardDiff.hessian(f, x) -# Hsparsity = sparse(H)[imarginal, imarginal] .!= 0 -# else -# println("Detecting Hessian sparsity via SparsityDetection...") -# Hsparsity = hessian_sparsity(f, x)[imarginal, imarginal] -# end -# Hcolors = matrix_colors(Hsparsity) -# - -struct ForwardColorHesCache{THS, THC, TI<:Integer, TD, TG} - Hsparsity::THS - Hcolors::THC - ncolors::TI - D::TD - Hcomp_buffer::TD - G::TG - δG::TG -end - -function ForwardColorHesCache(f, x, colorvec=1:length(x), sparsity=nothing) - D = hcat([float.(i .== colorvec) for i in 1:maximum(colorvec)]...) - Hcomp_buffer = similar(D) - G = zero(x) - δG = zero(x) - return ForwardColorHesCache(Hsparsity, colorvec, size(Hcolors, 2), D, Hcomp_buffer, G, δG) -end - -function sparse_hessian!(f, g!, θ, hessconfig::ForwardColorHesCache, δ=sqrt(eps(Float64))) - nc = hessconfig.ncolors - for j in one(nc):nc - g!(hessconfig.G, θ) - g!(hessconfig.δG, θ + δ * @view hessconfig.D[:, j]) - hessconfig.Hcomp_buffer[:, j] .= (hessconfig.δG .- hessconfig.G) ./ δ - end - ii, jj, vv = findnz(hessconfig.Hsparsity) - H = sparse(ii, jj, zeros(length(vv))) - for (i, j) in zip(ii, jj) - H[i, j] = hessconfig.Hcomp_buffer[i, hessconfig.Hcolors[j]] - end - return H -end - -function sparse_hessian!(f, θ, hessconfig::ForwardColorHesCache, δ=sqrt(eps(Float64))) - g!(G, θ) = ForwardDiff.gradient!(G, f, θ) - return sparse_hessian!(f, g!, θ, hessconfig, δ) -end diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl new file mode 100644 index 00000000..3ac00322 --- /dev/null +++ b/src/differentiation/compute_hessian_ad.jl @@ -0,0 +1,112 @@ + + + +# ForwardColorJacCache(f,x,_chunksize = nothing; +# dx = nothing, +# colorvec=1:length(x), +# sparsity = nothing) + + +# function HessianConfig(logdensity, imarginal, ijoint, forwarddiff_sparsity=false) +# if forwarddiff_sparsity +# println("Detecting Hessian sparsity via ForwardDiff...") +# H = ForwardDiff.hessian(f, x) +# Hsparsity = sparse(H)[imarginal, imarginal] .!= 0 +# else +# println("Detecting Hessian sparsity via SparsityDetection...") +# Hsparsity = hessian_sparsity(f, x)[imarginal, imarginal] +# end +# Hcolors = matrix_colors(Hsparsity) +# + +struct ForwardColorHesCache{THS, THC, TI<:Integer, TD, TG, T} + sparsity::THS + colors::THC + ncolors::TI + D::TD + buffer::TD + G::TG + dG::TG + dx::T +end + +function make_hessian_buffers(colorvec, x) + ncolors = maximum(colorvec) + D = hcat([float.(i .== colorvec) for i in 1:ncolors]...) + buffer = similar(D) + G = zero(x) + dG = zero(x) + return (;ncolors, D, buffer, G, dG) +end + +function ForwardColorHesCache(f, x::AbstractArray{T}, dx=sqrt(eps(T)), colorvec=1:length(x), sparsity=nothing) where {T<:Number} + ncolors, D, buffer, G, dG = make_hessian_buffers(colorvec, x) + return ForwardColorHesCache(sparsity, colorvec, ncolors, D, buffer, G, dG, dx) +end + +# function forwarddiff_color_hessian!(H, f, g!, x, hes_cache::ForwardColorHesCache, d=sqrt(eps(Float64))) +# for j in 1:hes_cache.ncolors +# g!(hes_cache.G, x) +# g!(hes_cache.dG, x + d * @view hes_cache.D[:, j]) +# hes_cache.buffer[:, j] .= (hes_cache.dG .- hes_cache.G) ./ d +# end +# ii, jj, vv = findnz(hes_cache.sparsity) +# for (i, j) in zip(ii, jj) +# H[i, j] = hes_cache.buffer[i, hes_cache.colors[j]] +# end +# return H +# end + +# function forwarddiff_color_hessian(H, f, g!, x, dG, colorvec, sparsity) +# cache = ForwardColorHesCache(f, x, dG, colorvec, sparsity) +# forwarddiff_color_hessian!(H, f, g!, ) +# end + +# function forwarddiff_color_hessian(f, x, hes_cache::ForwardColorHesCache, d=sqrt(eps(Float64))) +# g!(G, x) = ForwardDiff.gradient!(G, f, x) +# ii, jj, vv = findnz(hes_cache.sparsity) +# H = sparse(ii, jj, zeros(eltype(x), length(vv))) +# forwarddiff_color_hessian!(H, f, g!, x, hes_cache, d) +# return H +# end +################################# +function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, + f, + x::AbstractArray{<:Number}, + hes_cache::ForwardColorHesCache) + g!(G, x) = ForwardDiff.gradient!(G, f, x) + for j in 1:hes_cache.ncolors + g!(hes_cache.G, x) + g!(hes_cache.dG, x + hes_cache.dx * @view hes_cache.D[:, j]) + hes_cache.buffer[:, j] .= (hes_cache.dG .- hes_cache.G) ./ hes_cache.dx + end + ii, jj, vv = findnz(hes_cache.sparsity) + for (i, j) in zip(ii, jj) + H[i, j] = hes_cache.buffer[i, hes_cache.colors[j]] + end + return H +end + +function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, + f, + x::AbstractArray{<:Number}; + dx=sqrt(eps()), + colorvec=1:length(x), + sparsity=nothing) + hes_cache = ForwardColorHesCache(f, x, dx, colorvec, sparsity) + forwarddiff_color_hessian!(H, f, x, hes_cache) + return H +end + + +function forwarddiff_color_hessian(f::F, x::AbstractArray{<:Number}, hes_cache::ForwardColorHesCache) where F + ii, jj, vv = findnz(hes_cache.sparsity) + H = sparse(ii, jj, zeros(eltype(x), length(vv))) + forwarddiff_color_hessian!(H, f, x, hes_cache) + return H +end + +# forwarddiff_color_hessian(f::F, x::AbstractArray{<:Number}; colorvec, sparsity, jac_prototype, chunksize, dx) where F +# forwarddiff_color_hessian(f::F, x::AbstractArray{<:Number}, hes_cache::ForwardColorHesCache, hes_prototype) where F +# forwarddiff_color_hessian(H::AbstractMatrix{<:Number}, f::F, x::AbstractArray{<:Number}, hes_cache::ForwardColorHesCache) where F +# forwarddiff_color_hessian(H::AbstractArray{<:Number}, f::F, x::AbstractArray{<:Number}; colorvec, sparsity, hes_prototype, chunksize, dx) where F From c7688e159d6b6ec6adc7201118c1c2c8af41513a Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Wed, 15 Jun 2022 16:12:53 -0700 Subject: [PATCH 03/19] Sort out some method signatures, add tests --- src/differentiation/compute_hessian_ad.jl | 82 ++++++----------------- test/runtests.jl | 1 + test/test_sparse_hessian.jl | 68 +++++++++++-------- 3 files changed, 62 insertions(+), 89 deletions(-) diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index 3ac00322..3e2fbd9c 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -1,25 +1,4 @@ - - - -# ForwardColorJacCache(f,x,_chunksize = nothing; -# dx = nothing, -# colorvec=1:length(x), -# sparsity = nothing) - - -# function HessianConfig(logdensity, imarginal, ijoint, forwarddiff_sparsity=false) -# if forwarddiff_sparsity -# println("Detecting Hessian sparsity via ForwardDiff...") -# H = ForwardDiff.hessian(f, x) -# Hsparsity = sparse(H)[imarginal, imarginal] .!= 0 -# else -# println("Detecting Hessian sparsity via SparsityDetection...") -# Hsparsity = hessian_sparsity(f, x)[imarginal, imarginal] -# end -# Hcolors = matrix_colors(Hsparsity) -# - -struct ForwardColorHesCache{THS, THC, TI<:Integer, TD, TG, T} +struct ForwardColorHesCache{THS, THC, TI<:Integer, TD, TG, T<:Number} sparsity::THS colors::THC ncolors::TI @@ -44,37 +23,13 @@ function ForwardColorHesCache(f, x::AbstractArray{T}, dx=sqrt(eps(T)), colorvec= return ForwardColorHesCache(sparsity, colorvec, ncolors, D, buffer, G, dG, dx) end -# function forwarddiff_color_hessian!(H, f, g!, x, hes_cache::ForwardColorHesCache, d=sqrt(eps(Float64))) -# for j in 1:hes_cache.ncolors -# g!(hes_cache.G, x) -# g!(hes_cache.dG, x + d * @view hes_cache.D[:, j]) -# hes_cache.buffer[:, j] .= (hes_cache.dG .- hes_cache.G) ./ d -# end -# ii, jj, vv = findnz(hes_cache.sparsity) -# for (i, j) in zip(ii, jj) -# H[i, j] = hes_cache.buffer[i, hes_cache.colors[j]] -# end -# return H -# end - -# function forwarddiff_color_hessian(H, f, g!, x, dG, colorvec, sparsity) -# cache = ForwardColorHesCache(f, x, dG, colorvec, sparsity) -# forwarddiff_color_hessian!(H, f, g!, ) -# end -# function forwarddiff_color_hessian(f, x, hes_cache::ForwardColorHesCache, d=sqrt(eps(Float64))) -# g!(G, x) = ForwardDiff.gradient!(G, f, x) -# ii, jj, vv = findnz(hes_cache.sparsity) -# H = sparse(ii, jj, zeros(eltype(x), length(vv))) -# forwarddiff_color_hessian!(H, f, g!, x, hes_cache, d) -# return H -# end -################################# -function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, - f, - x::AbstractArray{<:Number}, - hes_cache::ForwardColorHesCache) - g!(G, x) = ForwardDiff.gradient!(G, f, x) +function forwarddiff_color_hessian!(H, + f, + g!, + x, + hes_cache::ForwardColorHesCache, + d=sqrt(eps(Float64))) for j in 1:hes_cache.ncolors g!(hes_cache.G, x) g!(hes_cache.dG, x + hes_cache.dx * @view hes_cache.D[:, j]) @@ -88,25 +43,32 @@ function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, end function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, - f, + f::F, + x::AbstractArray{<:Number}, + hes_cache::ForwardColorHesCache) where F + g!(G, x) = ForwardDiff.gradient!(G, f, x) + + H = forwarddiff_color_hessian!(H, f, g!, x, hes_cache, hes_cache.dx) + return H +end + +function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, + f::F, x::AbstractArray{<:Number}; dx=sqrt(eps()), colorvec=1:length(x), - sparsity=nothing) + sparsity=nothing) where F hes_cache = ForwardColorHesCache(f, x, dx, colorvec, sparsity) forwarddiff_color_hessian!(H, f, x, hes_cache) return H end -function forwarddiff_color_hessian(f::F, x::AbstractArray{<:Number}, hes_cache::ForwardColorHesCache) where F +function forwarddiff_color_hessian(f::F, + x::AbstractArray{<:Number}, + hes_cache::ForwardColorHesCache) where F ii, jj, vv = findnz(hes_cache.sparsity) H = sparse(ii, jj, zeros(eltype(x), length(vv))) forwarddiff_color_hessian!(H, f, x, hes_cache) return H end - -# forwarddiff_color_hessian(f::F, x::AbstractArray{<:Number}; colorvec, sparsity, jac_prototype, chunksize, dx) where F -# forwarddiff_color_hessian(f::F, x::AbstractArray{<:Number}, hes_cache::ForwardColorHesCache, hes_prototype) where F -# forwarddiff_color_hessian(H::AbstractMatrix{<:Number}, f::F, x::AbstractArray{<:Number}, hes_cache::ForwardColorHesCache) where F -# forwarddiff_color_hessian(H::AbstractArray{<:Number}, f::F, x::AbstractArray{<:Number}; colorvec, sparsity, hes_prototype, chunksize, dx) where F diff --git a/test/runtests.jl b/test/runtests.jl index d80e6e84..869c7878 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,7 @@ if GROUP == "All" @time @safetestset "Acyclic coloring" begin include("test_acyclic.jl") end @time @safetestset "Matrix to graph conversion" begin include("test_matrix2graph.jl") end @time @safetestset "AD using colorvec vector" begin include("test_ad.jl") end + @time @safetestset "Hessian colorvecs" begin include("test_sparse_hessian.jl") end @time @safetestset "Integration test" begin include("test_integration.jl") end @time @safetestset "Special matrices" begin include("test_specialmatrices.jl") end @time @safetestset "Jac Vecs and Hes Vecs" begin include("test_jaches_products.jl") end diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index 992bce86..24d4440e 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -1,36 +1,46 @@ - ## Hessian tests using SparsityDetection, SparseDiffTools +using ForwardDiff using LinearAlgebra, SparseArrays function fscalar(x) - return -dot(x, x) + return -dot(x, x) + 2 * x[2] * x[3]^2 end -x = randn(100) -hs = hessian_sparsity(fscalar, x) -col = matrix_colors(hs) -hescache = ForwardColorHesCache(fscalar, x, col) - -# ForwardColorJacCache(f,x,_chunksize = nothing; -# dx = nothing, -# colorvec=1:length(x), -# sparsity = nothing) - -# forwarddiff_color_hessian!(J::AbstractMatrix{<:Number}, -# f, -# x::AbstractArray{<:Number}; -# dx = nothing, -# colorvec = eachindex(x), -# sparsity = nothing) - -# forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, -# f, -# x::AbstractArray{<:Number}, -# jac_cache::ForwardColorJacCache) - -# jacout = forwarddiff_color_jacobian(g, x, -# dx = similar(x), -# colorvec = 1:length(x), -# sparsity = nothing, -# jac_prototype = nothing) # matrix w/ sparsity pattern \ No newline at end of file +x = randn(5) +sparsity = hessian_sparsity(fscalar, x) +colors = matrix_colors(tril(sparsity)) +dx = sqrt(eps()) +ncolors = maximum(colors) +D = hcat([float.(i .== colors) for i in 1:ncolors]...) +buffer = similar(D) +G = zero(x) +dG = 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 typeof(buffers_tup.buffer) == typeof(buffer) +@test buffers_tup.G == G +@test buffers_tup.dG == dG + +hescache1 = ForwardColorHesCache(sparsity, colors, ncolors, D, buffer, G, dG, dx) +hescache2 = ForwardColorHesCache(fscalar, x, dx, colors, sparsity) + + +Hforward = ForwardDiff.hessian(fscalar, x) +g(x) = ForwardDiff.gradient(fscalar, x) # allocating +g!(G, x) = ForwardDiff.gradient!(G, fscalar, x) # non-allocating + +for hescache in [hescache1, hescache2] + H = forwarddiff_color_hessian(fscalar, x, hescache1) + @test all(isapprox.(Hforward, H, rtol=1e-6)) + + H1 = similar(H) + forwarddiff_color_hessian!(H1, fscalar, g!, x, hescache, dx) + @test all(isapprox.(H1, H)) + + forwarddiff_color_hessian!(H1, fscalar, x, hescache) + @test all(isapprox.(H1, H)) +end From 237aeb816422d3739fcf9d4e0f6212763318291f Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Thu, 16 Jun 2022 14:50:55 -0700 Subject: [PATCH 04/19] fix typo --- test/test_sparse_hessian.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index 24d4440e..9c46e76e 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -34,7 +34,7 @@ g(x) = ForwardDiff.gradient(fscalar, x) # allocating g!(G, x) = ForwardDiff.gradient!(G, fscalar, x) # non-allocating for hescache in [hescache1, hescache2] - H = forwarddiff_color_hessian(fscalar, x, hescache1) + H = forwarddiff_color_hessian(fscalar, x, hescache) @test all(isapprox.(Hforward, H, rtol=1e-6)) H1 = similar(H) From f437dacc41e6787625a3aad738e6c367fa8862d6 Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Thu, 16 Jun 2022 14:52:49 -0700 Subject: [PATCH 05/19] Make perturbation of x in-place --- src/differentiation/compute_hessian_ad.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index 3e2fbd9c..ee0df573 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -32,7 +32,9 @@ function forwarddiff_color_hessian!(H, d=sqrt(eps(Float64))) for j in 1:hes_cache.ncolors g!(hes_cache.G, x) - g!(hes_cache.dG, x + hes_cache.dx * @view hes_cache.D[:, j]) + x .+= hes_cache.dx .* @view hes_cache.D[:, j] + g!(hes_cache.dG, x) + x .-= hes_cache.dx .* @view hes_cache.D[:, j] hes_cache.buffer[:, j] .= (hes_cache.dG .- hes_cache.G) ./ hes_cache.dx end ii, jj, vv = findnz(hes_cache.sparsity) From a1b69f6587b23d3bc4e5e5444194c0d4b018776c Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Fri, 17 Jun 2022 07:47:42 -0700 Subject: [PATCH 06/19] remove unnecessary SparsityDetection dep --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index aff4be61..79ae11b3 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,6 @@ Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Requires = "ae029012-a4dd-5104-9daa-d747884805df" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -SparsityDetection = "684fba80-ace3-11e9-3d08-3bc7ed6f96df" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" VertexSafeGraphs = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" From f926d5dc6ecd8e96cb6e0a5358b701b92826f619 Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Fri, 17 Jun 2022 12:05:24 -0700 Subject: [PATCH 07/19] Cache GradientConfig and simplify call signatures - Include a ForwardDiff.GradientConfig in the ForwardColorHesCache to eliminate allocations - Reduce the number of signatures for forwarddiff_color_hessian to the minimum needed to match forwarddiff_color_jacobian - Get rid of default arguments for colorvec and sparsity for now --- src/differentiation/compute_hessian_ad.jl | 72 ++++++++++++----------- 1 file changed, 37 insertions(+), 35 deletions(-) diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index ee0df573..58a94c56 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -1,41 +1,55 @@ -struct ForwardColorHesCache{THS, THC, TI<:Integer, TD, TG, T<:Number} +struct ForwardColorHesCache{THS, THC, TI<:Integer, TD, TGF, TGC, TG} sparsity::THS colors::THC ncolors::TI D::TD buffer::TD + grad!::TGF + grad_config::TGC G::TG dG::TG - dx::T end function make_hessian_buffers(colorvec, x) ncolors = maximum(colorvec) D = hcat([float.(i .== colorvec) for i in 1:ncolors]...) buffer = similar(D) - G = zero(x) - dG = zero(x) + G = similar(x) + dG = similar(x) return (;ncolors, D, buffer, G, dG) end -function ForwardColorHesCache(f, x::AbstractArray{T}, dx=sqrt(eps(T)), colorvec=1:length(x), sparsity=nothing) where {T<:Number} +function ForwardColorHesCache(f, + x::AbstractVector{<:Number}, + g!, + colorvec::AbstractVector{<:Integer}, + sparsity::AbstractMatrix{<:Integer}) ncolors, D, buffer, G, dG = make_hessian_buffers(colorvec, x) - return ForwardColorHesCache(sparsity, colorvec, ncolors, D, buffer, G, dG, dx) + grad_config = ForwardDiff.GradientConfig(f, x) + return ForwardColorHesCache(sparsity, colorvec, ncolors, D, buffer, g!, grad_config, G, dG) end -function forwarddiff_color_hessian!(H, +function ForwardColorHesCache(f, + x::AbstractVector{<:Number}, + colorvec::AbstractVector{<:Integer}, + sparsity::AbstractMatrix{<:Integer}) + g!(G, x, grad_config) = ForwardDiff.gradient!(G, f, x, grad_config) + return ForwardColorHesCache(f, x, g!, colorvec, sparsity) +end + + +function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, f, - g!, - x, - hes_cache::ForwardColorHesCache, - d=sqrt(eps(Float64))) + x::AbstractArray{<:Number}, + hes_cache::ForwardColorHesCache) + ϵ = sqrt(eps(eltype(x))) for j in 1:hes_cache.ncolors - g!(hes_cache.G, x) - x .+= hes_cache.dx .* @view hes_cache.D[:, j] - g!(hes_cache.dG, x) - x .-= hes_cache.dx .* @view hes_cache.D[:, j] - hes_cache.buffer[:, j] .= (hes_cache.dG .- hes_cache.G) ./ hes_cache.dx + 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) ./ ϵ end ii, jj, vv = findnz(hes_cache.sparsity) for (i, j) in zip(ii, jj) @@ -45,30 +59,18 @@ function forwarddiff_color_hessian!(H, end function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, - f::F, - x::AbstractArray{<:Number}, - hes_cache::ForwardColorHesCache) where F - g!(G, x) = ForwardDiff.gradient!(G, f, x) - - H = forwarddiff_color_hessian!(H, f, g!, x, hes_cache, hes_cache.dx) - return H -end - -function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, - f::F, - x::AbstractArray{<:Number}; - dx=sqrt(eps()), - colorvec=1:length(x), - sparsity=nothing) where F - hes_cache = ForwardColorHesCache(f, x, dx, colorvec, sparsity) + f, + x::AbstractArray{<:Number}, + colorvec::AbstractVector{<:Integer}, + sparsity::AbstractMatrix{<:Integer}) + hes_cache = ForwardColorHesCache(f, x, colorvec, sparsity) forwarddiff_color_hessian!(H, f, x, hes_cache) return H end - -function forwarddiff_color_hessian(f::F, +function forwarddiff_color_hessian(f, x::AbstractArray{<:Number}, - hes_cache::ForwardColorHesCache) where F + hes_cache::ForwardColorHesCache) ii, jj, vv = findnz(hes_cache.sparsity) H = sparse(ii, jj, zeros(eltype(x), length(vv))) forwarddiff_color_hessian!(H, f, x, hes_cache) From 09f8caf2da404696a7066845250bde8b7419f486 Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Fri, 17 Jun 2022 12:09:15 -0700 Subject: [PATCH 08/19] Add more tests --- test/test_sparse_hessian.jl | 39 ++++++++++++++++++++++++++----------- 1 file changed, 28 insertions(+), 11 deletions(-) diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index 9c46e76e..6b9a2445 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -7,10 +7,9 @@ function fscalar(x) return -dot(x, x) + 2 * x[2] * x[3]^2 end -x = randn(5) +x = randn(50) sparsity = hessian_sparsity(fscalar, x) colors = matrix_colors(tril(sparsity)) -dx = sqrt(eps()) ncolors = maximum(colors) D = hcat([float.(i .== colors) for i in 1:ncolors]...) buffer = similar(D) @@ -21,26 +20,44 @@ 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 buffers_tup.G == G -@test buffers_tup.dG == dG +@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) -hescache1 = ForwardColorHesCache(sparsity, colors, ncolors, D, buffer, G, dG, dx) -hescache2 = ForwardColorHesCache(fscalar, x, dx, colors, sparsity) - -Hforward = ForwardDiff.hessian(fscalar, x) +gconfig = ForwardDiff.GradientConfig(fscalar, x) g(x) = ForwardDiff.gradient(fscalar, x) # allocating -g!(G, x) = ForwardDiff.gradient!(G, fscalar, x) # non-allocating +g!(G, x, gconfig) = ForwardDiff.gradient!(G, fscalar, x, gconfig) # non-allocating + +hescache1 = ForwardColorHesCache(sparsity, colors, ncolors, D, buffer, g!, gconfig, G, dG) +hescache2 = ForwardColorHesCache(fscalar, x, g!, colors, sparsity) +hescache3 = ForwardColorHesCache(fscalar, x, colors, sparsity) -for hescache in [hescache1, hescache2] +for name in [:sparsity, :colors, :ncolors, :D] + @eval @test hescache1.$name == hescache2.$name + @eval @test hescache1.$name == hescache3.$name +end +for name in [:buffer, :G, :dG] + @eval @test size(hescache1.$name) == size(hescache2.$name) + @eval @test size(hescache1.$name) == size(hescache2.$name) + @eval @test eltype(hescache1.$name) == eltype(hescache3.$name) + @eval @test eltype(hescache1.$name) == eltype(hescache3.$name) +end + +Hforward = ForwardDiff.hessian(fscalar, x) +for hescache in [hescache1, hescache2, hescache3] H = forwarddiff_color_hessian(fscalar, x, hescache) @test all(isapprox.(Hforward, H, rtol=1e-6)) H1 = similar(H) - forwarddiff_color_hessian!(H1, fscalar, g!, x, hescache, dx) + forwarddiff_color_hessian!(H1, fscalar, x, colors, sparsity) @test all(isapprox.(H1, H)) forwarddiff_color_hessian!(H1, fscalar, x, hescache) @test all(isapprox.(H1, H)) end + + From 34c4a2284f230b28b0819f3cb0ab08d639c7813d Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Mon, 20 Jun 2022 12:00:49 -0700 Subject: [PATCH 09/19] add (inneficient) dense default arguments --- src/differentiation/compute_hessian_ad.jl | 25 +++++++++++++++++------ test/test_sparse_hessian.jl | 19 ++++++++++------- 2 files changed, 31 insertions(+), 13 deletions(-) diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index 58a94c56..8773213a 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -10,6 +10,16 @@ struct ForwardColorHesCache{THS, THC, TI<:Integer, TD, TGF, TGC, TG} dG::TG end +"""ALL-ONE!""" +struct DrBronnerArray{T,N} <: AbstractArray{T,N} + size::NTuple{N, Int} +end +DrBronnerArray(T, tup::NTuple{N, Int}) where N = DrBronnerArray{T, N}(tup) +Base.size(x::DrBronnerArray) = x.size +Base.getindex(x::DrBronnerArray{T, N}, i::Int) where {T, N} = one(T) +Base.getindex(x::DrBronnerArray{T, N}, I::Vararg{Int, N}) where {T, N} = one(T) + + function make_hessian_buffers(colorvec, x) ncolors = maximum(colorvec) D = hcat([float.(i .== colorvec) for i in 1:ncolors]...) @@ -22,18 +32,21 @@ end function ForwardColorHesCache(f, x::AbstractVector{<:Number}, g!, - colorvec::AbstractVector{<:Integer}, - sparsity::AbstractMatrix{<:Integer}) + colorvec::AbstractVector{<:Integer}=eachindex(x), + sparsity::Union{AbstractMatrix, Nothing}=nothing) ncolors, D, buffer, G, dG = make_hessian_buffers(colorvec, x) grad_config = ForwardDiff.GradientConfig(f, x) + if sparsity === nothing + sparsity = sparse(ones(length(x), length(x))) + end return ForwardColorHesCache(sparsity, colorvec, ncolors, D, buffer, g!, grad_config, G, dG) end function ForwardColorHesCache(f, x::AbstractVector{<:Number}, - colorvec::AbstractVector{<:Integer}, - sparsity::AbstractMatrix{<:Integer}) + colorvec::AbstractVector{<:Integer}=eachindex(x), + sparsity::Union{AbstractMatrix, Nothing}=nothing) g!(G, x, grad_config) = ForwardDiff.gradient!(G, f, x, grad_config) return ForwardColorHesCache(f, x, g!, colorvec, sparsity) end @@ -61,8 +74,8 @@ end function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, - colorvec::AbstractVector{<:Integer}, - sparsity::AbstractMatrix{<:Integer}) + colorvec::AbstractVector{<:Integer}=eachindex(x), + sparsity::Union{AbstractMatrix, Nothing}=nothing) hes_cache = ForwardColorHesCache(f, x, colorvec, sparsity) forwarddiff_color_hessian!(H, f, x, hes_cache) return H diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index 6b9a2445..f06db653 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -7,7 +7,7 @@ function fscalar(x) return -dot(x, x) + 2 * x[2] * x[3]^2 end -x = randn(50) +x = randn(5) sparsity = hessian_sparsity(fscalar, x) colors = matrix_colors(tril(sparsity)) ncolors = maximum(colors) @@ -35,29 +35,34 @@ g!(G, x, gconfig) = ForwardDiff.gradient!(G, fscalar, x, gconfig) # non-alloca hescache1 = ForwardColorHesCache(sparsity, colors, ncolors, D, buffer, g!, gconfig, G, dG) hescache2 = ForwardColorHesCache(fscalar, x, g!, colors, sparsity) hescache3 = ForwardColorHesCache(fscalar, x, colors, sparsity) +hescache4 = ForwardColorHesCache(fscalar, x) for name in [:sparsity, :colors, :ncolors, :D] @eval @test hescache1.$name == hescache2.$name @eval @test hescache1.$name == hescache3.$name + # hescache4 is the default dense version, so only first axis will match + @eval @test size(hescache1.$name, 1) == size(hescache4.$name, 1) end for name in [:buffer, :G, :dG] @eval @test size(hescache1.$name) == size(hescache2.$name) - @eval @test size(hescache1.$name) == size(hescache2.$name) - @eval @test eltype(hescache1.$name) == eltype(hescache3.$name) + @eval @test size(hescache1.$name) == size(hescache3.$name) + # hescache4 is the default dense version, so only first axis will match + @eval @test size(hescache1.$name, 1) == size(hescache4.$name, 1) + + @eval @test eltype(hescache1.$name) == eltype(hescache2.$name) @eval @test eltype(hescache1.$name) == eltype(hescache3.$name) + @eval @test eltype(hescache1.$name) == eltype(hescache4.$name) end Hforward = ForwardDiff.hessian(fscalar, x) -for hescache in [hescache1, hescache2, hescache3] +for (i, hescache) in enumerate([hescache1, hescache2, hescache3, hescache4]) H = forwarddiff_color_hessian(fscalar, x, hescache) @test all(isapprox.(Hforward, H, rtol=1e-6)) H1 = similar(H) - forwarddiff_color_hessian!(H1, fscalar, x, colors, sparsity) + forwarddiff_color_hessian!(H1, fscalar, x, collect(hescache.colors), hescache.sparsity) @test all(isapprox.(H1, H)) forwarddiff_color_hessian!(H1, fscalar, x, hescache) @test all(isapprox.(H1, H)) end - - From 46a037544a39e380764ba9a578667fd0b62b52df Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Tue, 21 Jun 2022 15:05:20 -0700 Subject: [PATCH 10/19] add colorvec/sparsity method for allocating forwarddiff_color_hessian --- src/differentiation/compute_hessian_ad.jl | 13 +++++++++++-- test/test_sparse_hessian.jl | 9 ++++++++- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index 8773213a..5d5b7bbf 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -84,8 +84,17 @@ end function forwarddiff_color_hessian(f, x::AbstractArray{<:Number}, hes_cache::ForwardColorHesCache) - ii, jj, vv = findnz(hes_cache.sparsity) - H = sparse(ii, jj, zeros(eltype(x), length(vv))) + H = convert.(eltype(x), hes_cache.sparsity) + forwarddiff_color_hessian!(H, f, x, hes_cache) + return H +end + +function forwarddiff_color_hessian(f, + x::AbstractArray{<:Number}, + colorvec::AbstractVector{<:Integer}=eachindex(x), + sparsity::Union{AbstractMatrix, Nothing}=nothing) + hes_cache = ForwardColorHesCache(f, x, colorvec, sparsity) + H = convert.(eltype(x), hes_cache.sparsity) forwarddiff_color_hessian!(H, f, x, hes_cache) return H end diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index f06db653..1e43ac5a 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -56,13 +56,20 @@ end Hforward = ForwardDiff.hessian(fscalar, x) for (i, hescache) in enumerate([hescache1, hescache2, hescache3, hescache4]) - H = forwarddiff_color_hessian(fscalar, x, hescache) + H = forwarddiff_color_hessian(fscalar, x, colors, sparsity) + H1 = forwarddiff_color_hessian(fscalar, x, hescache) + H2 = forwarddiff_color_hessian(fscalar, x) @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) forwarddiff_color_hessian!(H1, fscalar, x, collect(hescache.colors), hescache.sparsity) @test all(isapprox.(H1, H)) + forwarddiff_color_hessian!(H2, fscalar, x) + @test all(isapprox.(H2, H)) + forwarddiff_color_hessian!(H1, fscalar, x, hescache) @test all(isapprox.(H1, H)) end From fd954f1c1cdf15a1b7c78c2ee1b33266f596a580 Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Tue, 21 Jun 2022 15:12:02 -0700 Subject: [PATCH 11/19] remove DrBronnerArray started implementing an "All-One" array, decided not worth the complexity at this point --- src/differentiation/compute_hessian_ad.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index 5d5b7bbf..a5b5cb53 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -10,14 +10,6 @@ struct ForwardColorHesCache{THS, THC, TI<:Integer, TD, TGF, TGC, TG} dG::TG end -"""ALL-ONE!""" -struct DrBronnerArray{T,N} <: AbstractArray{T,N} - size::NTuple{N, Int} -end -DrBronnerArray(T, tup::NTuple{N, Int}) where N = DrBronnerArray{T, N}(tup) -Base.size(x::DrBronnerArray) = x.size -Base.getindex(x::DrBronnerArray{T, N}, i::Int) where {T, N} = one(T) -Base.getindex(x::DrBronnerArray{T, N}, I::Vararg{Int, N}) where {T, N} = one(T) function make_hessian_buffers(colorvec, x) From 71aedbed2d7a8255c57093d2c6b6d234895c7f03 Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Wed, 6 Jul 2022 17:40:05 -0700 Subject: [PATCH 12/19] Add checks on user-define gradient function Make sure it has the right signature, modify to accept GradientConfig if needed. Also some tests for it. --- src/differentiation/compute_hessian_ad.jl | 15 +++++++++++++- test/test_sparse_hessian.jl | 24 +++++++++++++++++------ 2 files changed, 32 insertions(+), 7 deletions(-) diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index a5b5cb53..65ecde3a 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -28,10 +28,23 @@ function ForwardColorHesCache(f, sparsity::Union{AbstractMatrix, Nothing}=nothing) ncolors, D, buffer, G, dG = make_hessian_buffers(colorvec, x) grad_config = ForwardDiff.GradientConfig(f, x) + + # If user supplied their own gradient function, make sure it has the right + # signature (i.e. g!(G, x) or g!(G, x, grad_config::ForwardDiff.GradientConfig)) + if ! hasmethod(g!, (typeof(G), typeof(G), typeof(grad_config))) + if ! hasmethod(g!, (typeof(G), typeof(G))) + throw(ArgumentError("Signature of `g!` must be either `g!(G, x)` or `g!(G, x, grad_config::ForwardDiff.GradientConfig)`")) + end + # define new method that takes a GradientConfig but doesn't use it + g1!(G, x, grad_config) = g!(G, x) + else + g1! = g! + end + if sparsity === nothing sparsity = sparse(ones(length(x), length(x))) end - return ForwardColorHesCache(sparsity, colorvec, ncolors, D, buffer, g!, grad_config, G, dG) + return ForwardColorHesCache(sparsity, colorvec, ncolors, D, buffer, g1!, grad_config, G, dG) end diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index 1e43ac5a..9bfac11d 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -35,27 +35,39 @@ g!(G, x, gconfig) = ForwardDiff.gradient!(G, fscalar, x, gconfig) # non-alloca hescache1 = ForwardColorHesCache(sparsity, colors, ncolors, D, buffer, g!, gconfig, G, dG) hescache2 = ForwardColorHesCache(fscalar, x, g!, colors, sparsity) hescache3 = ForwardColorHesCache(fscalar, x, colors, sparsity) -hescache4 = ForwardColorHesCache(fscalar, x) +# custom gradient function +hescache4 = ForwardColorHesCache(fscalar, x, (G, x) -> ForwardDiff.gradient!(G, fscalar, x), + colors, sparsity) +hescache5 = ForwardColorHesCache(fscalar, x) +# custom gradient has to have 2 or 3 arguments... +@test_throws ArgumentError ForwardColorHesCache(fscalar, x, (a) -> 1.0, colors, sparsity) +@test_throws ArgumentError ForwardColorHesCache(fscalar, x, (a, b, c, d) -> 1.0, colors, sparsity) +# ...and needs to accept (Vector, Vector, ForwardDiff.GradientConfig) +@test_throws ArgumentError ForwardColorHesCache(fscalar, x, (a::Int, b::Int) -> 1.0, colors, sparsity) +@test_throws ArgumentError ForwardColorHesCache(fscalar, x, (a::Int, b::Int, c::Int) -> 1.0, colors, sparsity) for name in [:sparsity, :colors, :ncolors, :D] @eval @test hescache1.$name == hescache2.$name @eval @test hescache1.$name == hescache3.$name - # hescache4 is the default dense version, so only first axis will match - @eval @test size(hescache1.$name, 1) == size(hescache4.$name, 1) + @eval @test hescache1.$name == hescache4.$name + # 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] @eval @test size(hescache1.$name) == size(hescache2.$name) @eval @test size(hescache1.$name) == size(hescache3.$name) - # hescache4 is the default dense version, so only first axis will match - @eval @test size(hescache1.$name, 1) == size(hescache4.$name, 1) + @eval @test size(hescache1.$name) == size(hescache4.$name) + # hescache5 is the default dense version, so only first axis will match + @eval @test size(hescache1.$name, 1) == size(hescache5.$name, 1) @eval @test eltype(hescache1.$name) == eltype(hescache2.$name) @eval @test eltype(hescache1.$name) == eltype(hescache3.$name) @eval @test eltype(hescache1.$name) == eltype(hescache4.$name) + @eval @test eltype(hescache1.$name) == eltype(hescache5.$name) end Hforward = ForwardDiff.hessian(fscalar, x) -for (i, hescache) in enumerate([hescache1, hescache2, hescache3, hescache4]) +for (i, hescache) in enumerate([hescache1, hescache2, hescache3, hescache4, hescache5]) H = forwarddiff_color_hessian(fscalar, x, colors, sparsity) H1 = forwarddiff_color_hessian(fscalar, x, hescache) H2 = forwarddiff_color_hessian(fscalar, x) From c2a1a6d9573d67c572b2c44c9a2226658e85678f Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Wed, 6 Jul 2022 17:41:20 -0700 Subject: [PATCH 13/19] Add some text and Hessian examples to README --- README.md | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/README.md b/README.md index f53021c4..29c25f83 100644 --- a/README.md +++ b/README.md @@ -207,6 +207,35 @@ function, otherwise it will become a dense matrix. If `jac_prototype` and function has a *square* Jacobian matrix. If it is not the case, please specify the shape of output by giving `dx`. +Similar functionality is available for Hessians. Given a scalar function `f(x)`, +a vector value for `x`, and a color vector and sparsity pattern, this can be +accomplished using `forwarddiff_color_hessian` or its in-place form `forwarddiff_color_hessian!`. + +```julia +H = forwarddiff_color_hessian(fscalar, x, colorvec, sparsity) +forwarddiff_color_hessian!(H, fscalar, x, colorvec, sparsity) +``` + +To avoid unnecessary allocations every time the Hessian is computed, +construct a `ForwardColorHesCache` beforehand: + +```julia +hescache = ForwardColorHesCache(f, x, colorvec, sparsity) +fowrwarddif_color_hessian!(H, f, x, hescache) +``` + +By default, these methods use a mix of numerical and automatic differentiation, +namely by taking finite differences of gradients calculated via ForwardDiff.jl. +Alternatively, if you have your own custom gradient function `g!`, you can specify +it as an argument to `ForwardColorHesCache`: + +```julia +hescache = ForwardColorHesCache(fscalar, x, g!, colorvec, sparsity) +``` +Note that any user-defined gradient needs to have the signature `g!(G, x)`, +i.e. updating the gradient `G` in place. + + ### Jacobian-Vector and Hessian-Vector Products Matrix-free implementations of Jacobian-Vector and Hessian-Vector products is From a30fd7d1c6c5e1332ca8a45a3c56dcfec12bd37b Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Wed, 20 Jul 2022 11:51:44 -0700 Subject: [PATCH 14/19] change sqrt to cbrt --- src/differentiation/compute_hessian_ad.jl | 2 +- test/test_sparse_hessian.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index 65ecde3a..1d2a609f 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -61,7 +61,7 @@ function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, hes_cache::ForwardColorHesCache) - ϵ = sqrt(eps(eltype(x))) + ϵ = 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] diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index 9bfac11d..9a8717ca 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -71,9 +71,9 @@ for (i, hescache) in enumerate([hescache1, hescache2, hescache3, hescache4, hesc H = forwarddiff_color_hessian(fscalar, x, colors, sparsity) H1 = forwarddiff_color_hessian(fscalar, x, hescache) H2 = forwarddiff_color_hessian(fscalar, x) - @test all(isapprox.(Hforward, H, rtol=1e-6)) - @test all(isapprox.(H, H1, rtol=1e-6)) - @test all(isapprox.(H2, H1, rtol=1e-6)) + @test all(isapprox.(Hforward, H, rtol=1e-4)) + @test all(isapprox.(H, H1, rtol=1e-4)) + @test all(isapprox.(H2, H1, rtol=1e-4)) H1 = similar(H) forwarddiff_color_hessian!(H1, fscalar, x, collect(hescache.colors), hescache.sparsity) From 13b2dca2489c2f9442e623edda1c712224014ddc Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Wed, 20 Jul 2022 12:07:44 -0700 Subject: [PATCH 15/19] Condense signature for default gradient function Also change order to put optional argument last --- README.md | 2 +- src/differentiation/compute_hessian_ad.jl | 14 ++------------ test/test_sparse_hessian.jl | 14 +++++++------- 3 files changed, 10 insertions(+), 20 deletions(-) diff --git a/README.md b/README.md index 29c25f83..d196d6a3 100644 --- a/README.md +++ b/README.md @@ -230,7 +230,7 @@ Alternatively, if you have your own custom gradient function `g!`, you can speci it as an argument to `ForwardColorHesCache`: ```julia -hescache = ForwardColorHesCache(fscalar, x, g!, colorvec, sparsity) +hescache = ForwardColorHesCache(fscalar, x, colorvec, sparsity, g!) ``` Note that any user-defined gradient needs to have the signature `g!(G, x)`, i.e. updating the gradient `G` in place. diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index 1d2a609f..1ab0d459 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -23,9 +23,9 @@ end function ForwardColorHesCache(f, x::AbstractVector{<:Number}, - g!, colorvec::AbstractVector{<:Integer}=eachindex(x), - sparsity::Union{AbstractMatrix, Nothing}=nothing) + 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) grad_config = ForwardDiff.GradientConfig(f, x) @@ -47,16 +47,6 @@ function ForwardColorHesCache(f, return ForwardColorHesCache(sparsity, colorvec, ncolors, D, buffer, g1!, grad_config, G, dG) 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) - return ForwardColorHesCache(f, x, g!, colorvec, sparsity) -end - - function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index 9a8717ca..da4598f1 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -33,18 +33,18 @@ 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) -hescache2 = ForwardColorHesCache(fscalar, x, g!, colors, sparsity) +hescache2 = ForwardColorHesCache(fscalar, x, colors, sparsity, g!) hescache3 = ForwardColorHesCache(fscalar, x, colors, sparsity) # custom gradient function -hescache4 = ForwardColorHesCache(fscalar, x, (G, x) -> ForwardDiff.gradient!(G, fscalar, x), - colors, sparsity) +hescache4 = ForwardColorHesCache(fscalar, x, colors, sparsity, + (G, x) -> ForwardDiff.gradient!(G, fscalar, x)) hescache5 = ForwardColorHesCache(fscalar, x) # custom gradient has to have 2 or 3 arguments... -@test_throws ArgumentError ForwardColorHesCache(fscalar, x, (a) -> 1.0, colors, sparsity) -@test_throws ArgumentError ForwardColorHesCache(fscalar, x, (a, b, c, d) -> 1.0, colors, sparsity) +@test_throws ArgumentError ForwardColorHesCache(fscalar, x, colors, sparsity, (a) -> 1.0) +@test_throws ArgumentError ForwardColorHesCache(fscalar, x, colors, sparsity, (a, b, c, d) -> 1.0) # ...and needs to accept (Vector, Vector, ForwardDiff.GradientConfig) -@test_throws ArgumentError ForwardColorHesCache(fscalar, x, (a::Int, b::Int) -> 1.0, colors, sparsity) -@test_throws ArgumentError ForwardColorHesCache(fscalar, x, (a::Int, b::Int, c::Int) -> 1.0, colors, sparsity) +@test_throws ArgumentError ForwardColorHesCache(fscalar, x, colors, sparsity, (a::Int, b::Int) -> 1.0,) +@test_throws ArgumentError ForwardColorHesCache(fscalar, x, colors, sparsity, (a::Int, b::Int, c::Int) -> 1.0) for name in [:sparsity, :colors, :ncolors, :D] @eval @test hescache1.$name == hescache2.$name From d4911a7ae25a68cda9f05a791dfe46eb86058526 Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Wed, 20 Jul 2022 12:24:02 -0700 Subject: [PATCH 16/19] "safe" option to make sure hessian is overwritten --- src/differentiation/compute_hessian_ad.jl | 8 +++++--- test/test_sparse_hessian.jl | 10 ++++++++++ 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index 1ab0d459..9724d622 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -10,8 +10,6 @@ struct ForwardColorHesCache{THS, THC, TI<:Integer, TD, TGF, TGC, TG} dG::TG end - - function make_hessian_buffers(colorvec, x) ncolors = maximum(colorvec) D = hcat([float.(i .== colorvec) for i in 1:ncolors]...) @@ -50,7 +48,8 @@ end function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, - hes_cache::ForwardColorHesCache) + hes_cache::ForwardColorHesCache; + safe = true) ϵ = cbrt(eps(eltype(x))) for j in 1:hes_cache.ncolors hes_cache.grad!(hes_cache.G, x, hes_cache.grad_config) @@ -60,6 +59,9 @@ function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, hes_cache.buffer[:, j] .= (hes_cache.dG .- hes_cache.G) ./ ϵ end ii, jj, vv = findnz(hes_cache.sparsity) + if safe + fill!(H, false) + end for (i, j) in zip(ii, jj) H[i, j] = hes_cache.buffer[i, hes_cache.colors[j]] end diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index da4598f1..c6268795 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -84,4 +84,14 @@ for (i, hescache) in enumerate([hescache1, hescache2, hescache3, hescache4, hesc forwarddiff_color_hessian!(H1, fscalar, x, hescache) @test all(isapprox.(H1, H)) + + forwarddiff_color_hessian!(H1, fscalar, x, hescache, safe=false) + @test all(isapprox.(H1, H)) + + # confirm unsafe is faster + t_safe = minimum(@elapsed(forwarddiff_color_hessian!(H1, fscalar, x, hescache, safe=true)) + for _ in 1:100) + t_unsafe = minimum(@elapsed(forwarddiff_color_hessian!(H1, fscalar, x, hescache, safe=false)) + for _ in 1:100) + @test t_unsafe <= t_safe end From 0a473647a398516e4f279bbf961581418c197165 Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Thu, 28 Jul 2022 16:40:17 -0700 Subject: [PATCH 17/19] Rename forwarddiff_color_hessian to numauto_color_hessian --- README.md | 12 ++++++------ src/SparseDiffTools.jl | 4 ++-- src/differentiation/compute_hessian_ad.jl | 14 +++++++------- test/test_sparse_hessian.jl | 18 +++++++++--------- 4 files changed, 24 insertions(+), 24 deletions(-) diff --git a/README.md b/README.md index d196d6a3..3511ac6a 100644 --- a/README.md +++ b/README.md @@ -207,13 +207,13 @@ function, otherwise it will become a dense matrix. If `jac_prototype` and function has a *square* Jacobian matrix. If it is not the case, please specify the shape of output by giving `dx`. -Similar functionality is available for Hessians. Given a scalar function `f(x)`, -a vector value for `x`, and a color vector and sparsity pattern, this can be -accomplished using `forwarddiff_color_hessian` or its in-place form `forwarddiff_color_hessian!`. +Similar functionality is available for Hessians, using finite differences of forward derivatives. Given a scalar function `f(x)`, a vector value for `x`, +and a color vector and sparsity pattern, this can be accomplished using +`numauto_color_hessian` or its in-place form `numauto_color_hessian!`. ```julia -H = forwarddiff_color_hessian(fscalar, x, colorvec, sparsity) -forwarddiff_color_hessian!(H, fscalar, x, colorvec, sparsity) +H = numauto_color_hessian(fscalar, x, colorvec, sparsity) +numauto_color_hessian!(H, fscalar, x, colorvec, sparsity) ``` To avoid unnecessary allocations every time the Hessian is computed, @@ -221,7 +221,7 @@ construct a `ForwardColorHesCache` beforehand: ```julia hescache = ForwardColorHesCache(f, x, colorvec, sparsity) -fowrwarddif_color_hessian!(H, f, x, hescache) +numauto_color_hessian!(H, f, x, hescache) ``` By default, these methods use a mix of numerical and automatic differentiation, diff --git a/src/SparseDiffTools.jl b/src/SparseDiffTools.jl index 70807db2..b9ae599f 100644 --- a/src/SparseDiffTools.jl +++ b/src/SparseDiffTools.jl @@ -28,8 +28,8 @@ export contract_color, forwarddiff_color_jacobian!, forwarddiff_color_jacobian, ForwardColorJacCache, - forwarddiff_color_hessian!, - forwarddiff_color_hessian, + numauto_color_hessian!, + numauto_color_hessian, ForwardColorHesCache, auto_jacvec,auto_jacvec!, num_jacvec,num_jacvec!, diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index 9724d622..686ae490 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -45,7 +45,7 @@ function ForwardColorHesCache(f, return ForwardColorHesCache(sparsity, colorvec, ncolors, D, buffer, g1!, grad_config, G, dG) end -function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, +function numauto_color_hessian!(H::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, hes_cache::ForwardColorHesCache; @@ -68,30 +68,30 @@ function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, return H end -function forwarddiff_color_hessian!(H::AbstractMatrix{<:Number}, +function numauto_color_hessian!(H::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, colorvec::AbstractVector{<:Integer}=eachindex(x), sparsity::Union{AbstractMatrix, Nothing}=nothing) hes_cache = ForwardColorHesCache(f, x, colorvec, sparsity) - forwarddiff_color_hessian!(H, f, x, hes_cache) + numauto_color_hessian!(H, f, x, hes_cache) return H end -function forwarddiff_color_hessian(f, +function numauto_color_hessian(f, x::AbstractArray{<:Number}, hes_cache::ForwardColorHesCache) H = convert.(eltype(x), hes_cache.sparsity) - forwarddiff_color_hessian!(H, f, x, hes_cache) + numauto_color_hessian!(H, f, x, hes_cache) return H end -function forwarddiff_color_hessian(f, +function numauto_color_hessian(f, x::AbstractArray{<:Number}, colorvec::AbstractVector{<:Integer}=eachindex(x), sparsity::Union{AbstractMatrix, Nothing}=nothing) hes_cache = ForwardColorHesCache(f, x, colorvec, sparsity) H = convert.(eltype(x), hes_cache.sparsity) - forwarddiff_color_hessian!(H, f, x, hes_cache) + numauto_color_hessian!(H, f, x, hes_cache) return H end diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index c6268795..3f2d6086 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -68,30 +68,30 @@ end Hforward = ForwardDiff.hessian(fscalar, x) for (i, hescache) in enumerate([hescache1, hescache2, hescache3, hescache4, hescache5]) - H = forwarddiff_color_hessian(fscalar, x, colors, sparsity) - H1 = forwarddiff_color_hessian(fscalar, x, hescache) - H2 = forwarddiff_color_hessian(fscalar, x) + 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)) H1 = similar(H) - forwarddiff_color_hessian!(H1, fscalar, x, collect(hescache.colors), hescache.sparsity) + numauto_color_hessian!(H1, fscalar, x, collect(hescache.colors), hescache.sparsity) @test all(isapprox.(H1, H)) - forwarddiff_color_hessian!(H2, fscalar, x) + numauto_color_hessian!(H2, fscalar, x) @test all(isapprox.(H2, H)) - forwarddiff_color_hessian!(H1, fscalar, x, hescache) + numauto_color_hessian!(H1, fscalar, x, hescache) @test all(isapprox.(H1, H)) - forwarddiff_color_hessian!(H1, fscalar, x, hescache, safe=false) + numauto_color_hessian!(H1, fscalar, x, hescache, safe=false) @test all(isapprox.(H1, H)) # confirm unsafe is faster - t_safe = minimum(@elapsed(forwarddiff_color_hessian!(H1, fscalar, x, hescache, safe=true)) + t_safe = minimum(@elapsed(numauto_color_hessian!(H1, fscalar, x, hescache, safe=true)) for _ in 1:100) - t_unsafe = minimum(@elapsed(forwarddiff_color_hessian!(H1, fscalar, x, hescache, safe=false)) + t_unsafe = minimum(@elapsed(numauto_color_hessian!(H1, fscalar, x, hescache, safe=false)) for _ in 1:100) @test t_unsafe <= t_safe end From 40e877ccceab3dfc50f3b6bd95c894752ab44443 Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Fri, 29 Jul 2022 14:55:54 -0700 Subject: [PATCH 18/19] Change to centered difference --- src/differentiation/compute_hessian_ad.jl | 23 ++++++++++++----------- test/test_sparse_hessian.jl | 22 +++++++++++----------- 2 files changed, 23 insertions(+), 22 deletions(-) diff --git a/src/differentiation/compute_hessian_ad.jl b/src/differentiation/compute_hessian_ad.jl index 686ae490..f9a07183 100644 --- a/src/differentiation/compute_hessian_ad.jl +++ b/src/differentiation/compute_hessian_ad.jl @@ -6,17 +6,17 @@ 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, @@ -24,7 +24,7 @@ function ForwardColorHesCache(f, 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 diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index 3f2d6086..c8a7c3a8 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -13,8 +13,8 @@ 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 @@ -22,17 +22,17 @@ buffers_tup = SparseDiffTools.make_hessian_buffers(colors, x) @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) From e1c8eafe9224003f5748e2d1078852277b96bcda Mon Sep 17 00:00:00 2001 From: Sam Urmy Date: Fri, 29 Jul 2022 15:36:10 -0700 Subject: [PATCH 19/19] Remove timing test that randomly fails once in a while --- test/test_sparse_hessian.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/test/test_sparse_hessian.jl b/test/test_sparse_hessian.jl index c8a7c3a8..d0308da1 100644 --- a/test/test_sparse_hessian.jl +++ b/test/test_sparse_hessian.jl @@ -88,10 +88,13 @@ for (i, hescache) in enumerate([hescache1, hescache2, hescache3, hescache4, hesc numauto_color_hessian!(H1, fscalar, x, hescache, safe=false) @test all(isapprox.(H1, H)) - # confirm unsafe is faster - t_safe = minimum(@elapsed(numauto_color_hessian!(H1, fscalar, x, hescache, safe=true)) - for _ in 1:100) - t_unsafe = minimum(@elapsed(numauto_color_hessian!(H1, fscalar, x, hescache, safe=false)) - for _ in 1:100) - @test t_unsafe <= t_safe + # the following tests usually pass, but once in a while don't (it's not a big difference + # in timing on these small matrices, and sometimes its less than the timing variability). + # Commenting out for now to avoid rare stochastic test failures. + # # confirm unsafe is faster + # t_safe = minimum(@elapsed(numauto_color_hessian!(H1, fscalar, x, hescache, safe=true)) + # for _ in 1:100) + # t_unsafe = minimum(@elapsed(numauto_color_hessian!(H1, fscalar, x, hescache, safe=false)) + # for _ in 1:100) + # @test t_unsafe <= t_safe end