Hessians

Hessians are a crucial ingredient in NewtonSolvers and SimpleSolvers.NewtonOptimizerStates.

using SimpleSolvers
using LinearAlgebra: norm

x = rand(3)
obj = MultivariateObjective(x -> norm(x - vcat(0., 0., 1.))  ^ 2, x)
hes = HessianAutodiff(obj, x)

HessianAutodiff{Float64, Main.var"#1#2", Matrix{Float64}, ForwardDiff.HessianConfig{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3}, 3}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3}}}}(Main.var"#1#2"(), [NaN NaN NaN; NaN NaN NaN; NaN NaN NaN], ForwardDiff.HessianConfig{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3}, 3}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3}}}(ForwardDiff.JacobianConfig{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3}}}((Partials(1.0, 0.0, 0.0), Partials(0.0, 1.0, 0.0), Partials(0.0, 0.0, 1.0)), ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3}[Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(6.95171462023273e-310,6.95174885566977e-310,6.9517146202335e-310,6.95174885566977e-310), Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(6.9517146202343e-310,6.95174885566977e-310,6.9517444547575e-310,6.95174885567056e-310), Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(6.9517444547575e-310,6.95174885567056e-310,6.95171462024064e-310,6.95174885566977e-310)]), ForwardDiff.GradientConfig{ForwardDiff.Tag{Main.var"#1#2", Float64}, ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3}, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3}, 3}}}((Partials(Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(1.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(0.0,0.0,0.0,0.0)), Partials(Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(1.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(0.0,0.0,0.0,0.0)), Partials(Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(1.0,0.0,0.0,0.0))), ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3}, 3}[Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(0.0,5.0e-324,5.0e-324,1.5e-323),Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(1.5e-323,2.5e-323,2.5e-323,1.5e-323),Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(4.0e-323,4.0e-323,5.0e-323,5.0e-323),Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(4.0e-323,6.4e-323,6.4e-323,7.4e-323)), Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(8.0e-323,8.4e-323,9.0e-323,9.4e-323),Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(1.0e-322,1.04e-322,1.1e-322,1.14e-322),Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(1.2e-322,1.24e-322,1.3e-322,1.33e-322),Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(1.4e-322,1.43e-322,1.4e-322,1.53e-322)), Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(1.4e-322,1.63e-322,1.7e-322,1.73e-322),Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(1.73e-322,1.83e-322,1.73e-322,1.93e-322),Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(1.93e-322,2.03e-322,1.93e-322,2.1e-322),Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}}(2.1e-322,2.2e-322,2.1e-322,2.3e-322))])))

An instance of HessianAutodiff stores a Hessian matrix:

hes.H

3×3 Matrix{Float64}:
 NaN  NaN  NaN
 NaN  NaN  NaN
 NaN  NaN  NaN

The instance of HessianAutodiff can be called:

hes(x)

3×3 Matrix{Float64}:
 2.0          0.0           0.0
 0.0          2.0          -1.11022e-16
 5.55112e-17  1.11022e-16   2.0

Or equivalently with:

update!(hes, x)

This updates hes.H:

hes.H

3×3 Matrix{Float64}:
 2.0          0.0           0.0
 0.0          2.0          -1.11022e-16
 5.55112e-17  1.11022e-16   2.0

BFGS Hessian

using SimpleSolvers: initialize!
hes = HessianBFGS(obj, x)
initialize!(hes, x)

HessianBFGS{Float64, Vector{Float64}, Matrix{Float64}, MultivariateObjective{Float64, Vector{Float64}, Main.var"#1#2", GradientAutodiff{Float64, Main.var"#1#2", ForwardDiff.GradientConfig{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Main.var"#1#2", Float64}, Float64, 3}}}}, Float64, Vector{Float64}}}(MultivariateObjective (for vector-valued quantities only the first component is printed):

    f(x)              = NaN 
    g(x)₁             = 3.82e-01 
    x_f₁              = NaN 
    x_g₁              = 1.91e-01 
    number of f calls = 0 
    number of g calls = 1 
, [NaN, NaN, NaN], [0.19090669902576285, 0.5256623915420473, 0.3905882754313441], [NaN, NaN, NaN], [NaN, NaN, NaN], [0.3818133980515257, 1.0513247830840946, -1.2188234491373118], [NaN, NaN, NaN], [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0])

For computational reasons we save the inverse of the Hessian, it can be accessed by calling inv:

inv(hes)

3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

Similarly to HessianAutodiff we can call SimpleSolvers.update!:

using SimpleSolvers: update!

update!(hes, x)