In [1]:
using Polyopt
We define two variables $x$ and $z$,
In [2]:
x, z = variables(["x", "z"]);
and the function $f(x,z)=4 x^2 + xz - 4z^2 - \frac{21}{10}x^4 + 4z^4 + \frac{1}{3}x^6$.
In [3]:
f = 4*x^2 + x*z - 4*z^2 - 21//10*x^4 + 4*z^4 + 1//3*x^6;
To minimize $f(x,z)$ we form a moment relaxation of order 3,
In [4]:
prob = momentprob(3, f);
and we solve it,
In [5]:
X, Z, t, y, solsta = solve_mosek(prob);
The optimal lower bound is
In [6]:
t
Out[6]:
If an optimal polynomial solution is found, then $y$ corresponds to the monomial basis vector
In [7]:
[prob.basis y]
Out[7]:
Because of multiple global minima of $f(x,z)$ this is not the case.
We perturb $f(x,z)$ to find a global minimizer,
In [8]:
fp = f + 1e-4*(x+z)
Out[8]:
In [9]:
probp = momentprob(3, fp);
In [10]:
Xp, Zp, tp, yp, solsta = solve_mosek(probp);
We get (approximately) the same optimal lower bound
In [11]:
tp
Out[11]:
In [12]:
[probp.basis yp]
Out[12]:
Now we have found a global minimizer, and can we verify optimality by evaluating $f(x,z)$ at that point
In [13]:
xo = [yp[8], yp[2]]
Out[13]:
In [14]:
Polyopt.evalpoly(f, xo)
Out[14]:
In [15]:
using PyPlot
n = 100
X = linspace(-1.5, 1.5, n)
Z = linspace(-1, 1, n)
Xgrid = repmat(X',n,1)
Zgrid = repmat(Z,1,n)
T = zeros(n,n)
for i in 1:n
for j in 1:n
T[i:i,j:j] = Polyopt.evalpoly(f, [X[j];Z[i]])
end
end
fig = figure(figsize=(15,10))
ax = fig[:gca](projection = "3d")
ax[:plot_surface](Xgrid, Zgrid, T, rstride=2, edgecolors="k", cstride=2, cmap=ColorMap("gray"), linewidth=0.25)
ax[:view_init](20.0, 55.0)
ax[:contour](Xgrid, Zgrid, T, 15, offset=-2.0, colors="black")
ax[:plot]( [xo[1]], [xo[2]], -2.0, "ro")
xlabel("x")
ylabel("z")
Out[15]:
The monomial basis vector defining the problem stored is stored in prob.basis,
In [16]:
prob.basis
Out[16]:
The coefficients for $f$ are stored in prob.obj,
In [17]:
dot(prob.obj, prob.basis)
Out[17]:
The symmetric matrices defining the moment matrix are stored in prob.mom[1] column by column,
In [18]:
size(prob.mom[1])
Out[18]:
In [19]:
reshape(prob.mom[1]*prob.basis, 10, 10)
Out[19]:
and similarly the moment matrix evaluated at yp is
In [20]:
M=reshape(prob.mom[1]*yp, 10, 10)
Out[20]:
In [21]:
eigvals(Symmetric(M,:L))
Out[21]:
The (primal) solution is a semidefinite matrix
In [22]:
Xp[1]
Out[22]:
In [23]:
eigvals(Xp[1])
Out[23]:
If we store Xp[1] in vectorized form (conforming with prob.mom),
In [24]:
xp = vec(Xp[1])
Out[24]:
then we can verify that
In [25]:
prob.mom[1]'*xp - (prob.obj - tp*eye(28,1))
Out[25]:
In [ ]: