# Laplacian for a minimized function

**URL:** <https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836>\
**Category:** General Usage\
**Tags:** optimization, forwarddiff\
**Created:** [September 14, 2023, 5:10am UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836 "2023-09-14T05:10:36Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![yosuga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yosuga/32/52638_2.png) [@yosuga](https://discourse.julialang.org/u/yosuga)\
**Post date:** [September 14, 2023, 5:10am UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836/1 "2023-09-14T05:10:36Z")

</div>

I’m trying to evaluate the following `laplacian(b)`.

```julia
using Optim
using ForwardDiff

function f(a,b)
    sin(a[1])*(cos(a[2]) + cos(b[1]+b[2]))
end

function g(b)
    h(a) = f(a,b)
    res = optimize(h,zeros(2))
    return Optim.minimum(res)
end

function laplacian(b)
    tr(ForwardDiff.hessian(g,[b[1],b[2]]))
end

```

`f(a,b)` and `g(b)` work but `laplacian(b)` dosen’t.

```julia
MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(g), Int64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(g), Int64}, Float64, 2}, 2})

```

How should I change the code?  
Any help would be greatly appreciated. Thanks!

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [September 14, 2023, 7:37am UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836/2 "2023-09-14T07:37:04Z")

</div>

The issue is that the optimization starts from `zeros(2)`, which is implicitly equivalent to `zeros(Float64, 2)`. If you want ForwardDiff to be able to differentiate through your optimization algorithm, you have to allow arbitrary types everywhere, like this:

```julia
function g(b::Vector{T}) where {T}
    h(a) = f(a, b)
    res = optimize(h, zeros(T, 2))
    return Optim.minimum(res)
end

```

I checked and the Laplacian works now.  
However, can I ask if your actual problem really looks like that MWE? Depending on the size of the input space and the autodiff mode (forward or reverse), you may want to take a look at other packages like [DiffOpt.jl](https://github.com/jump-dev/DiffOpt.jl) or [ImplicitDifferentiation.jl](https://github.com/gdalle/ImplicitDifferentiation.jl).

---

<div class="post-metadata">

**Author:** ![yosuga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yosuga/32/52638_2.png) [@yosuga](https://discourse.julialang.org/u/yosuga)\
**Post date:** [September 14, 2023, 8:44am UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836/3 "2023-09-14T08:44:40Z")

</div>

Thank you for your reply. I checked the Laplacian works.

My actual setup is the following.

```julia
function R(μ,E)
    function dR(θ) 
        T = exp(μ[1]+im*θ[1]) + exp(-μ[1]-im*θ[1]) + exp(μ[2]+im*θ[2]) + exp(-μ[2]-im*θ[2])
        Td = (exp(μ[1]+im*θ[1]) + exp(-μ[1]-im*θ[1]))*(exp(μ[2]+im*θ[2]) + exp(-μ[2]-im*θ[2]))
        Γ = exp(μ[1]+im*θ[1]) - exp(-μ[1]-im*θ[1]) + exp(μ[2]+im*θ[2]) - exp(-μ[2]-im*θ[2])
        return (1/2π)*(1/2π)*log(abs(t*T+td*Td+γ*Γ-E[1]-im*E[2]))
    end
    hcubature(dR,[0,0],[2π,2π],rtol=1e-5)[1]
end

function Φ(E::Vector{T}) where {T}
    Ronkin(μ) = R(μ,E)
    res = optimize(Ronkin,zeros(T,2))
    return Optim.minimum(res)
end

function DOS(E)
    tr(ForwardDiff.hessian(Φ,E))
end

```

here `t=1.0, td=0.5, γ=0.2`.  
I need to evaluate `DOS(E)` for a specific region `-L<E[1],E[2}<L`.  
I’m not familiar with autodiff. What approach is appropriate for this kind of problem?  
As far as I checked, the above `DOS(E)` did not return the correct value, while [the difference approximation](https://en.wikipedia.org/wiki/Five-point_stencil) for `Φ(E::Vector{T})` provided an insufficient but nearly accurate value.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [September 14, 2023, 8:57am UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836/4 "2023-09-14T08:57:39Z")

</div>

If I understand correctly you only have 2 parameters? In that case, forward-mode autodiff (with ForwardDiff.jl) should be your best bet.

I’m a bit surprised that it does not return the “right” value, can you provide a fully runnable example (with imports and numerical values) where the results differ?

---

<div class="post-metadata">

**Author:** ![yosuga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yosuga/32/52638_2.png) [@yosuga](https://discourse.julialang.org/u/yosuga)\
**Post date:** [September 14, 2023, 9:12am UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836/5 "2023-09-14T09:12:36Z")

</div>

Yes, I have 2 parameters, E[1] and E[2], for `DOS(E)`.  
Here is my code.

```julia
using LinearAlgebra
using HCubature
using Optim
using ForwardDiff

function R(μ,E)
    function dR(θ) 
        T = exp(μ[1]+im*θ[1]) + exp(-μ[1]-im*θ[1]) + exp(μ[2]+im*θ[2]) + exp(-μ[2]-im*θ[2])
        Td = (exp(μ[1]+im*θ[1]) + exp(-μ[1]-im*θ[1]))*(exp(μ[2]+im*θ[2]) + exp(-μ[2]-im*θ[2]))
        Γ = exp(μ[1]+im*θ[1]) - exp(-μ[1]-im*θ[1]) + exp(μ[2]+im*θ[2]) - exp(-μ[2]-im*θ[2])
        return (1/2π)*(1/2π)*log(abs(t*T+td*Td+γ*Γ-E[1]-im*E[2]))
    end
    hcubature(dR,[0,0],[2π,2π],rtol=1e-5)[1]
end

function Φ(E::Vector{T}) where {T}
    Ronkin(μ) = R(μ,E)
    res = optimize(Ronkin,zeros(T,2))
    return Optim.minimum(res)
end

function DOS1(E)
    (1/2π)*tr(ForwardDiff.hessian(Φ,E))
end

function DOS2(E)
    δx = δy = 0.1
    Exy = E
    Exp = E.+[δx,0]
    Exm = E.+[-δx,0]
    Eyp = E.+[0,δy]
    Eym = E.+[0,-δy]
    return (1/2π)*((Φ(Exp)-2Φ(Exy)+Φ(Exm))/(δx^2) + (Φ(Eyp)-2Φ(Exy)+Φ(Eym))/(δy^2))
end

```

```julia
t = 1.0
td = 0.5
γ = 0.2
E = [0.0,0.0]

```

```julia
DOS1(E),DOS2(E)
(-1.4135798584282297e-16, 0.14388277024025198)

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [September 14, 2023, 9:17am UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836/6 "2023-09-14T09:17:08Z")

</div>

It looks like either your optimization algorithm or your integration algorithm are nearly constant around your solution, which means the algorithmic derivatives will be close to zero. Gonna take a look

---

<div class="post-metadata">

**Author:** ![yosuga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yosuga/32/52638_2.png) [@yosuga](https://discourse.julialang.org/u/yosuga)\
**Post date:** [September 14, 2023, 5:57pm UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836/7 "2023-09-14T17:57:44Z")

</div>

FYI, it may relate to the nature of `R(μ,E)`. Namely, `R(μ,E)` could take a minimum in a certain region (not at a single point) as a function of `(μ[1],μ[2])` for a fixed `E`, say `E=E0`. See [this topic](https://discourse.julialang.org/t/finding-the-minimizer-region-of-a-function/103819) as an example.  
The important features are

1. When `R(μ,E0)` has a minimizer region, the integrand of `R(μ,E0)` is well-defined in there.
2. When `R(μ,E0)` has a minimizer point, the integrand of `R(μ,E0)` shows a logarithmic singularity at that point.

although, in both cases, `R(μ,E0)` itself is well-defined and takes a finite value.

My observation tells that `DOS1(E)` doesn’t work when a minimizer point exists. For instance,

1. At `E0 = [0.0,0.0]`, which provides a minimizer point,  
`(DOS1(E0),DOS2(E0))=(2.8271597168564594e-16, 0.1418357989357821)`.

2. At `E0 = [-3.0,3.0]`, which gives a minimizer region,  
`(DOS1(E0),DOS2(E0))=(0.0, -0.00018425790497728805)`.

For a given `E0`, you can see whether `R(μ,E0)` has a minimizer point or a minimizer region by calculating [Amoeba](https://en.wikipedia.org/wiki/Amoeba_(mathematics)) of the polynomial, `p(z,w) = t*T+td*Td+γ*Γ-E0[1]-im*E0[2]`, where `z = exp(μ[1]+im*θ[1])` and `w = exp(μ[2]+im*θ[2])`.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [September 14, 2023, 6:52pm UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836/8 "2023-09-14T18:52:39Z")

</div>

But then how do you interpret the result of finite differences in the case of a minimizer point? Does it make any more mathematical sense?

---

<div class="post-metadata">

**Author:** ![yosuga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yosuga/32/52638_2.png) [@yosuga](https://discourse.julialang.org/u/yosuga)\
**Post date:** [September 14, 2023, 7:57pm UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836/9 "2023-09-14T19:57:14Z")

</div>

Roughly speaking, `DOS(E)` behaves as `DOS(E) ~ Σ_i δ(E-Ei)`, here `R(μ,Ei)` gives a minimizer point and `δ(E)` is a 2d delta function. This means that if `E=E0` corresponds to a minimizer region, `DOS(E0)` is zero. On the other hand, in the case of a minimizer point, `DOS(E0)` diverges. Therefore, the result of finite differences can be understood as the suppressed divergence by the discretization of Laplacian.  
So it is expected that `DOS1(E0)` may not match `DOS2(E0)` but it returns non-zero (maybe quite large) value in the case of a minimizer point.

Actually, I’m trying to reproduce Fig.4(a) in [this paper](https://arxiv.org/abs/2212.11743). See Eq.(7),(8),(18),(19),(20) and Sec.IV.

---

<div class="post-metadata">

**Author:** ![Reuel\_Dsouza](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/reuel_dsouza/32/207810_2.png) [@Reuel\_Dsouza](https://discourse.julialang.org/u/Reuel_Dsouza)\
**Post date:** [May 8, 2024, 12:50pm UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836/10 "2024-05-08T12:50:33Z")

</div>

Hi, did you have any luck in reproducing these figures? I am trying this too.

---

<div class="post-metadata">

**Author:** ![yosuga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yosuga/32/52638_2.png) [@yosuga](https://discourse.julialang.org/u/yosuga)\
**Post date:** [May 9, 2024, 7:36am UTC](https://discourse.julialang.org/t/laplacian-for-a-minimized-function/103836/11 "2024-05-09T07:36:15Z")

</div>

Hi, Reuel.  
I confirmed that they can be roughly reproduced by straightforward calculation. However, I haven’t found a numerically efficient way with less error to calculate them yet.