# 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, γ)
optimize!(model)
``````

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

``````println(objective_value(model))
value(γ)
3.0000000209530935
``````

When I tried to evaluate the optimal solution

``````ν3 = moment_matrix(constr)
extractatoms(ν3,1e-3)
``````

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)
600
``````

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)
``````

with

``````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`?