# Polynomial optimization: not getting the correct solution with extractatoms

**URL:** <https://discourse.julialang.org/t/polynomial-optimization-not-getting-the-correct-solution-with-extractatoms/42635>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, optimization, polynomials\
**Created:** [July 6, 2020, 7:47pm UTC](https://discourse.julialang.org/t/polynomial-optimization-not-getting-the-correct-solution-with-extractatoms/42635 "2020-07-06T19:47:08Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![aliceee](https://avatars.discourse-cdn.com/v4/letter/a/8baadc/32.png) [@aliceee](https://discourse.julialang.org/u/aliceee)\
**Post date:** [July 6, 2020, 7:47pm UTC](https://discourse.julialang.org/t/polynomial-optimization-not-getting-the-correct-solution-with-extractatoms/42635/1 "2020-07-06T19:47:08Z")

</div>

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

```julia
# 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

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

```

When I tried to evaluate the optimal solution

```julia
ν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.

```julia
subs(f, x1=>0, x2=>0)
600

```

I also check

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

```

and get

```julia
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!

---

<div class="post-metadata">

**Author:** ![blegat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blegat/32/217090_2.png) [@blegat](https://discourse.julialang.org/u/blegat)\
**Post date:** [July 7, 2020, 8:55am UTC](https://discourse.julialang.org/t/polynomial-optimization-not-getting-the-correct-solution-with-extractatoms/42635/2 "2020-07-07T08:55:07Z")

</div>

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

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

```

with

```julia
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`.

---

<div class="post-metadata">

**Author:** ![aliceee](https://avatars.discourse-cdn.com/v4/letter/a/8baadc/32.png) [@aliceee](https://discourse.julialang.org/u/aliceee)\
**Post date:** [July 13, 2020, 5:40pm UTC](https://discourse.julialang.org/t/polynomial-optimization-not-getting-the-correct-solution-with-extractatoms/42635/3 "2020-07-13T17:40:14Z")

</div>

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