Laplacian for a minimized function

General Usage
,
1

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

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.

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!

1 Like
2

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:

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 or ImplicitDifferentiation.jl.

3

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

My actual setup is the following.

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 for Φ(E::Vector{T}) provided an insufficient but nearly accurate value.

4

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?

5

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

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

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

DOS1(E),DOS2(E)
(-1.4135798584282297e-16, 0.14388277024025198)
6

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

7

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 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 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]).

8

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?

9

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. See Eq.(7),(8),(18),(19),(20) and Sec.IV.

10

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

11

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.