add HesVecGrad

Author: ChrisRackauckas

src/SparseDiffTools.jl

@@ -16,7 +16,9 @@ export  contract_color,
         num_hesvec,num_hesvec!,
         numauto_hesvec,numauto_hesvec!,
         autonum_hesvec,autonum_hesvec!,
-        JacVec,HesVec
+        num_hesvecgrad,num_hesvecgrad!,
+        auto_hesvecgrad,auto_hesvecgrad!,
+        JacVec,HesVec,HesVecGrad
 
 include("coloring/high_level.jl")
 include("coloring/contraction_coloring.jl")

src/differentiation/jaches_products.jl

@@ -108,6 +108,42 @@ function autonum_hesvec(f,x,v)
     partials.(g(Dual{DeivVecTag}.(x, v)), 1)
 end
 
+function num_hesvecgrad!(du,g,x,v,
+                     cache2 = similar(v),
+                     cache3 = similar(v))
+    T = eltype(x)
+    # Should it be min? max? mean?
+    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
+    @. x += ϵ*v
+    g(cache2,x)
+    @. x -= 2ϵ*v
+    g(cache3,x)
+    @. du = (cache2 - cache3)/(2ϵ)
+end
+
+function num_hesvecgrad(g,x,v)
+    T = eltype(x)
+    # Should it be min? max? mean?
+    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
+    x += ϵ*v
+    gxp = g(x)
+    x -= 2ϵ*v
+    gxm = g(x)
+    (gxp - gxm)/(2ϵ)
+end
+
+function auto_hesvecgrad!(du,g,x,v,
+                     cache2 = ForwardDiff.Dual{DeivVecTag}.(x, v),
+                     cache3 = ForwardDiff.Dual{DeivVecTag}.(x, v))
+    cache2 .= Dual{DeivVecTag}.(x, v)
+    g(cache3,cache2)
+    du .= partials.(cache3, 1)
+end
+
+function auto_hesvecgrad(g,x,v)
+    partials.(g(Dual{DeivVecTag}.(x, v)), 1)
+end
+
 ### Operator Forms
 
 mutable struct JacVec{F,T1,T2,uType}
@@ -174,3 +210,34 @@ function LinearAlgebra.mul!(du::AbstractVector,L::HesVec,v::AbstractVector)
         num_hesvec!(du,L.f,L.u,v,L.cache1,L.cache2,L.cache3)
     end
 end
+
+mutable struct HesVecGrad{G,T1,T2,uType}
+    g::G
+    cache1::T1
+    cache2::T2
+    u::uType
+    autodiff::Bool
+end
+
+function HesVecGrad(g,u::AbstractArray;autodiff=true)
+    if autodiff
+        cache1 = ForwardDiff.Dual{DeivVecTag}.(u, u)
+        cache2 = ForwardDiff.Dual{DeivVecTag}.(u, u)
+    else
+        cache1 = similar(u)
+        cache2 = similar(u)
+    end
+    HesVecGrad(g,cache1,cache2,u,autodiff)
+end
+
+Base.size(L::HesVecGrad) = (length(L.cache2),length(L.cache2))
+Base.size(L::HesVecGrad,i::Int) = length(L.cache2)
+Base.:*(L::HesVecGrad,x::AbstractVector) = L.autodiff ? auto_hesvecgrad(L.g,L.u,x) : num_hesvecgrad(L.g,L.u,x)
+
+function LinearAlgebra.mul!(du::AbstractVector,L::HesVecGrad,v::AbstractVector)
+    if L.autodiff
+        auto_hesvecgrad!(du,L.g,L.u,v,L.cache1,L.cache2)
+    else
+        num_hesvecgrad!(du,L.g,L.u,v,L.cache1,L.cache2)
+    end
+end

test/test_jaches_products.jl

@@ -31,6 +31,20 @@ f(u) = sum(u.^2)
 @test autonum_hesvec!(du, f, x, v, similar(v), cache1, cache2) ≈ ForwardDiff.hessian(f,x)*v rtol=1e-2
 @test autonum_hesvec(f, x, v) ≈ ForwardDiff.hessian(f,x)*v rtol=1e-8
 
+function g(x)
+      DiffEqDiffTools.finite_difference_gradient(f,x)
+end
+function g(dx,x)
+      DiffEqDiffTools.finite_difference_gradient!(dx,f,x)
+end
+@test num_hesvecgrad!(du, g, x, v) ≈ ForwardDiff.hessian(f,x)*v rtol=1e-2
+@test num_hesvecgrad!(du, g, x, v, similar(v), similar(v)) ≈ ForwardDiff.hessian(f,x)*v rtol=1e-2
+@test num_hesvecgrad(g, x, v) ≈ ForwardDiff.hessian(f,x)*v rtol=1e-2
+
+@test auto_hesvecgrad!(du, g, x, v) ≈ ForwardDiff.hessian(f,x)*v rtol=1e-2
+@test auto_hesvecgrad!(du, g, x, v, cache1, cache2) ≈ ForwardDiff.hessian(f,x)*v rtol=1e-2
+@test auto_hesvecgrad(g, x, v) ≈ ForwardDiff.hessian(f,x)*v rtol=1e-2
+
 f(du,u) = mul!(du,A,u)
 f(u) = A*u
 L = JacVec(f,x)
@@ -72,3 +86,17 @@ L.u .= v
 ### Integration test with IterativeSolvers
 out = similar(v)
 gmres!(out, L, v)
+
+L = HesVecGrad(g,x,autodiff=false)
+@test L*x ≈ num_hesvec(f, x, x)
+@test L*v ≈ num_hesvec(f, x, v)
+@test mul!(du,L,v) ≈ num_hesvec(f, x, v) rtol=1e-2
+L.u .= v
+@test mul!(du,L,v) ≈ num_hesvec(f, v, v) rtol=1e-2
+
+L = HesVecGrad(g,x)
+@test L*x ≈ numauto_hesvec(f, x, x)
+@test L*v ≈ numauto_hesvec(f, x, v)
+@test mul!(du,L,v) ≈ numauto_hesvec(f, x, v) rtol=1e-8
+L.u .= v
+@test mul!(du,L,v) ≈ numauto_hesvec(f, v, v) rtol=1e-8