Evaluating the Hessian of the Lagrangian efficiently

General Usage Performance
2

Maybe ModelingToolkit.jl would be helpful? For a system like this that consists of simple equations, it can do symbolic operations rather than AD.

using ModelingToolkit

@parameters t y[1:3] λ[1:3]
@variables x[1:3](t)
D = Differential(t)

eqs = [
    D(x[1]) ~ x[1] * x[2] + x[3] * y[1],  # constraint 1
    D(x[2]) ~ x[1] * y[2] * y[3],         # constraint 2
    D(x[3]) ~ x[1] * x[2] + x[3] + y[2],  # constraint 3
]

ℒ = sum(λ.*ModelingToolkit.rhss(eqs))

ℋ = ModelingToolkit.hessian(ℒ, [x;y])

Produces

6×6 Matrix{Num}:
       0  λ₁ + λ₃   0   0     y₃*λ₂     y₂*λ₂
 λ₁ + λ₃        0   0   0         0         0
       0        0   0  λ₁         0         0
       0        0  λ₁   0         0         0
   y₃*λ₂        0   0   0         0  λ₂*x₁(t)
   y₂*λ₂        0   0   0  λ₂*x₁(t)         0

Which can be compiled to a function

hess, _ = ModelingToolkit.build_function(ℋ, [x;y;λ]; expression=Val{false})
hess(rand(9))

6×6 Matrix{Float64}:
 0.0        1.56516  0.0       0.0       0.0490337   0.0661009
 1.56516    0.0      0.0       0.0       0.0         0.0
 0.0        0.0      0.0       0.927298  0.0         0.0
 0.0        0.0      0.927298  0.0       0.0         0.0
 0.0490337  0.0      0.0       0.0       0.0         0.00383799
 0.0661009  0.0      0.0       0.0       0.00383799  0.0
2 Likes