# Evaluating the Hessian of the Lagrangian efficiently

**URL:** <https://discourse.julialang.org/t/evaluating-the-hessian-of-the-lagrangian-efficiently/54659>\
**Category:** Performance\
**Tags:** modelingtoolkit\
**Created:** [February 4, 2021, 11:04pm UTC](https://discourse.julialang.org/t/evaluating-the-hessian-of-the-lagrangian-efficiently/54659 "2021-02-04T23:04:10Z")\
**Posts on this page:** 1\
**Showing post:** 2

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [February 5, 2021, 12:44am UTC](https://discourse.julialang.org/t/evaluating-the-hessian-of-the-lagrangian-efficiently/54659/2 "2021-02-05T00:44:21Z")

</div>

Maybe [ModelingToolkit.jl](https://mtk.sciml.ai/stable/) would be helpful? For a system like this that consists of simple equations, it can do symbolic operations rather than AD.

```julia
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

```julia
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

```julia
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

```

---

_[View the full topic](https://discourse.julialang.org/t/evaluating-the-hessian-of-the-lagrangian-efficiently/54659)._