Polynomial optimization: not getting the correct solution with extractatoms

Hi, I’m trying to get the optimal solution of a polynomial optimization problem by calculating the dual.

# Create symbolic variables (not JuMP decision variables)
@polyvar x1 x2
# Create a Sum of Squares JuMP model with the Mosek solver
model = SOSModel(with_optimizer(Mosek.Optimizer))
# Create a JuMP decision variable for the lower bound
@variable(model, γ)
# @variable(model, x)
# f(x) is the Goldstein-Price function
f1 = x1+x2+1
f2 = 19-14*x1+3*x1^2-14*x2+6*x1*x2+3*x2^2
f3 = 2*x1-3*x2
f4 = 18-32*x1+12*x1^2+48*x2-36*x1*x2+27*x2^2
f = (1+f1^2*f2)*(30+f3^2*f4)
@constraint(model, constr, f >= γ)
@objective(model, Max, γ)

The problem is successfully solved to optimal. The optimal objective value can be evaluated by


When I tried to evaluate the optimal solution

ν3 = moment_matrix(constr)

and I got

Atomic measure on the variables x1, x2 with 1 atoms:
at [0.0, 0.0] with weight 0.9999999999872671

However, it is not the optimal solution since its corresponding objective value is 600.

subs(f, x1=>0, x2=>0)

I also check

MultivariateMoments.computesupport!(ν3, 1e-3)

and get

Algebraic Set defined by 12 equalities
 -x2 = 0
 -x1 = 0
 -x2^2 = 0
 -x1*x2 = 0
 -x1^2 = 0
 -x2^3 = 0
 -x1*x2^2 = 0
 -x1^2*x2 = 0
 -x1^3 = 0
 -x1^2*x2^2 + 0.49988336573392234*x1*x2^3 + 1.4998072823025776*x2^4 = 0
 -x1^3*x2 + 1.7498451395778885*x1*x2^3 + 0.7496347168248607*x2^4 = 0
 -x1^4 + 1.6244349952561488*x1*x2^3 + 2.6243893730401324*x2^4 = 0

while clearly, the true optimal solution (0,-1) does not satisfy these constraints.
Can anyone help me? Thanks in advance!

It seems you have chosen a too loose tolerance 1e-3, with 1e-4 you get that the moment matrix is not atomic.
There is no check done in extractatoms to verify that the atoms found have moments close to the moments you gave.
You can check it yourself either by verifying that evaluating the atoms does not give you the same value of the objective of by comparing

ν3 = moment_matrix(constr)
μ3 = measure(ν3)


atoms = extractatoms(ν3, 1e-3)
μ = measure(atoms, μ3.x)

Here, you clearly see that it does not match as all moments are zero for μ except for the moment of the monomial 1 and that’s not the case for μ3.

Got it. Thanks for clarification. In this case, when the moment matrix is not atomic ( the result when tolerance is set to 1e-4), how should I evaluate the optimal solutions for x1 and x2?