Finding the minimizer region of a function

[yosuga]

As an example, let f be a polynomial,

f(z,w) = z + z^-1 + w - w^-1 + 2zw - 5,

where

z = exp(x[1]+iθ[1]), w = exp(x[2]+iθ[2]), 
x,θ ∈ R^2

Then we can define the Ronkin function as follows.

function R(x)
    function dR(θ) 
        #z = exp(x[1]+im*θ[1]),w = exp(x[2]+im*θ[2])
        f = exp(x[1]+im*θ[1])+exp(-x[1]-im*θ[1])+exp(x[2]+im*θ[2])-exp(-x[2]-im*θ[2])+2exp(x[1]+x[2]+im*θ[1]+im*θ[2])-5
        return (1/2π)*(1/2π)*log(abs(f))
    end
    hcubature(dR,[0,0],[2π,2π],rtol=1e-6)[1]
end

Some general features of the Ronkin function are

1. It is convex in entire (x[1],x[2])-space.
2. There would exist an (finite or infinite) area where the Ronkin function takes its minimum.

For the above polynomial, we can see the minimizer region via a contour plot.

The pink region is called Amoeba, where the Ronkin function is strongly convex. And the Ronkin function takes the minimum value at any points in the “vacuole” of Amoeba.

Question:
How can we get the informations of the minimizer region?
Such as sufficient data points in the region, the area of the region, the boundary of the region, etc.

I tried Optim.jl but it returned only a single point in the vacuole (the red point in the figure).

[jd-foster]

An initial thought would be that if you know the coordinates (x_1^*, x_2^*) of the red point, then you would need to solve the implicit curve problem R(x_1, x_2) - R(x_1^*, x_2^*) = 0?

[yosuga]

Mathematically, yes. Then one gets the solutions filling the “vacuole” of Amoeba.
I need a highly accurate and efficient numerical method.

[jd-foster]

I have a tracing method for the boundary, but wouldn’t go so far as to say it’s

@yosuga

a highly accurate and efficient numerical method

[yosuga]

Good!
How did you do that?

[kellertuer]

I fear that any optimizer usually returns a minimiser but not all (i.e. the set). So What you are actually looking for is the indicator function of the set as @jd-foster described it.

Concerning the highly accurate and efficient – can you describe a bit what actual outcome you want to have? A function that answers fast (and accurate) whether a point x1,x2 is in the set?
Such a function would be very easy: Take the minimiser you have from Optim.jl and compare function values.

edit: A wanted to give that idea a short try, but your example above is not reproducible, where is hcubature from?

[jd-foster]

MWE:

using LinearAlgebra
using HCubature
using Optim
using ForwardDiff
using PlotlyJS

function R(x)
    function dR(θ)
        z = exp(x[1]+im*θ[1])
        w = exp(x[2]+im*θ[2])
        f = z + z^-1 + w - w^-1 + 2z*w - 5
        return (1/2π)*(1/2π)*log(abs(f))
    end
    hcubature(dR,[0,0],[2π,2π],rtol=1e-6)[1]
end

## Plots:

xs = ys = -6:0.1:6;
x_y = [ [x,y] for x in xs, y in ys];
zs = R.(x_y);

trace1 = contour(; z=zs, x=xs, y=ys,
        line_color=:white,
        contours=Dict(:start=>1, :end=>13, :size=>0.02))

fig = plot(trace1, Layout( yaxis=attr(scaleanchor="x")))

[kellertuer]

Yes, thanks, and you could combine this with an Optim run to just contour the minimiser region, or use other contour plot libraries – if just the figure is the goal.

Or to be precise one could also use the same set of functions to trigger a contour as in Contour Plots · Plots

[jd-foster]

Consider one-dimensional “slices” across the function (using the above MWE as a starting point), say along x1 = 0:

Q(x2) = t -> R([ t, x2])
S(x1) = t -> R([x1,  t])

xc = 0.0

F = S(xc)
p1 = plot(scatter(; x=xs, y=F.(xs), name="Function at $(xc)"));
p1

G = t -> ForwardDiff.derivative(F, t)
p2 = plot(scatter(; x=xs, y=G.(xs), name="Gradient at $(xc)"));

H = t -> ForwardDiff.derivative(G, t)
p3 = plot(scatter(; x=xs, y=H.(xs), name="2nd deriv. at $(xc)"));

[p1; p2; p3]

By convexity, the gradient is monotonic non-decreasing, where the “vacuole” is the “step” in the centre. If we know an interior point of the middle step, then tracking when the function changes indicates the boundary of those regions. Of course, the numerical behaviour means that tuning tolerances for distinguishing “zero” from “non-zero” becomes dependent on the function.

nonzero_tol = 0.1
E_L = Tuple[]
E_R = Tuple[] # {Float64,Float64}

x_scan = -2:0.05:1;
y_scan = -2:0.05:2;

for y1 in y_scan
    G = t -> ForwardDiff.derivative(S(y1), t)
    D = collect(zip(x_scan, G.(x_scan)))
    left_edge = D[findlast(x -> x[2] < -nonzero_tol, collect(D))]
    right_edge = D[findfirst(x -> x[2] > nonzero_tol, collect(D))]
    push!(E_L, (y1, left_edge...))
    push!(E_R, (y1, right_edge...))
end

traceE_L = scatter(x=getindex.(E_L,2), y=getindex.(E_L,1))
traceE_R = scatter(x=getindex.(E_R,2), y=getindex.(E_R,1))

fig = plot(trace1, Layout( yaxis=attr(scaleanchor="x"), showlegend=false))
add_trace!(fig, traceE_L)
add_trace!(fig, traceE_R)
display(fig)

[yosuga]

Thank you!
Does it work for a very small “vacuole”? In other words, can we detect a narrow middle “step” of a gradient? Such a situation can be obtained by f->f(a=-2.05) = z + z^-1 + w - w^-1 + 2z*w - a. Of course, I know it depends on the tolerance as you mentioned, but my final goal is in there.

[yosuga]

By introducing a parameter to R(x) as f->f(a) = z + z^-1 + w - w^-1 + 2z*w - a, the shape of Amoeba become to depend on a. My goal is to distinguish the existence and the absence of “vacuole”. For example, I want to differentiate a=2.0 (no vacuole) and a=2.05.

One of the simplest way is to see a picture of R(x) directly, but this is hard when we have many a’s. The other way is collecting the points near the minimum obtained from Optim.jl as you said, but this is also hard for a small “vacuole”. So I’m asking here to get a hint to establish a criterion for detecting the existence and the absence of “vacuole” such as the boundary of Amoeba.

[kellertuer]

Oh, that I would approach analytically maybe? In 1D on a plateau you would have the second derivative to be zero, for a strictly convex function, you would not have that. You can do the same in the multivariate case looking at Hessians.

[yosuga]

I understand the theoretical expectations, but I’m not sure whether a numerical calculation works or not. Since the integrand of the Ronkin function has a logarithmic singularity on Amoeba, one may get an incorrect result about optimization. Actually, I have failed to evaluate Laplacian (of a-space) for a different model.

But I tried to calculate Laplacian (of x-space) for the present cases, the results look good. So if Optim.jl always gives a correct minimizer, namely a point in the vacuole, one would distinguish the difference. Thank you.

[kellertuer]

That looks good. Depending on what is more stable, you could also check the determinant of the Hessian.

[yosuga]

Yes, I’ll consider that. Thanks again.