Is there any function of TypedPolynomials to get a polynomial vector with fixed degree?

[maihoanganh]

In DynamicPolynomials, I can use function monomials to get a polynomial vector with fixed degree, For example, run monomials([x y], 2) to get [x^2 xy y^2]. Could you show me such function if it is available in TypedPolynomials? Otherwise could you investigate this? In polynomial optimization, such function is necessary to generate sum-of-squares.

[blegat]

You can do

monos = monomials([x, y], 2)
using SumOfSquares
model = Model() 
@variable(model, p, Poly(monos))

See the SumOfSquares documentation for more information.

[maihoanganh]

I tried:

using TypedPolynomials

@polyvar x y
monos = monomials([x, y], 2)

then received the feedback:

UndefVarError: monomials not defined

[blegat]

You need to do using MultivariatePolynomials

[maihoanganh]

After using MultivariatePolynomials, ...

MethodError: no method matching monomials(::Array{Monomial{(x, y),2},1}, ::Int64)

[blegat]

This should work:

julia> using SumOfSquares

julia> using TypedPolynomials

julia> @polyvar x y
(x, y)

julia> monos = monomials((x, y), 2)
3-element Array{Monomial{(x, y),2},1}:
 x²
 xy
 y²

julia> model = Model()
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> @variable(model, p, Poly(monos))
(noname)x² + (noname)xy + (noname)y²

[maihoanganh]

Thank for your answers. But I prefer using only TypedPolynomials to combining with SumOfSquares.

[blegat]

Then what should the coefficients of the polynomial be ?

[maihoanganh]

I use JuMP to generate Gram matrix for sum of square. See PV_TSSOS in my github if you care.

[blegat]

You can simply do

using JuMP
model = Model()
@variable(model, coefs[1:length(monos)])
p = polynomial(coefs, monos)

Note that you can generate a Gram Matrix with SumOfSquares using

using SumOfSquares
model = Model()
@variable(model, q, SOSPoly(monos))

For instance, instead of @constraint(model, p in SOSCone()), if you already know the monomials for the gram matrix (e.g. and want to skip the Newton Polytope step), you can do @variable(model, q in SOSPoly(monos)) and @constraint(model, p == q) instead.

[maihoanganh]

Thank you for your answer. But I had known about it before.

There are several reasons because of that I have not used SumOfSquares.jl (as well as MomentOpt.jl) in my paper https://arxiv.org/abs/1911.11428 with Lasserre and Magron and also PV_TSSOS:

N. H. A. Mai, J. Wang, V. Magron, and J.-B. Lasserre. Exploiting term sparsity for polynomial optimization on non-compact semialgebraic sets, 2020. Forthcoming.

For instant, a bug is related with SemialgebraicSets.jl when I write codes relying on SumOfSquares package for the first paper. The following is an excerpt of my email to Weisser (for an example of polynomial optimization problems in the first paper):

"... We found some wrong output data of POPs with equality constraints using MomentOpt while JuMP and Yalmip scripts gives us true result. MomentOpt script below is an example when solving the POP:

-1 = min x1^6+x2^6+x3^6+3x1^2x2^2x3^2-x1^2(x2^4+x3^4)-x2^2(x3^4+x1^4)-x3^2(x1^4+x2^4) s.t. x1+x2+x3-1 = 0

and the wrong approximation of optimal value is -1.2421e-9. Could you find the reason of this issue and investigate it? ..."

I do not know whether this bug is fixed now. If it has been fixed, I will use SumOfSquares.jl in the forthcoming paper.

[maihoanganh]

Let us consider Id 13 of Table 2 and 3 in my paper https://arxiv.org/abs/1911.11428 with Lasserre and Magron. The exact optimal value is -1. Models by JuMP and Yalmip with order 2 provide a same approximate optimal value -0.9999. But I just run:

using DynamicPolynomials, MomentOpt, JuMP, SemialgebraicSets, MosekTools, CPUTime

start = time()

# Define polnomial variables
@polyvar x1 x2 x3

# Polynomial to optimize 
f = (x1^2+x2^2-2)*(x1^2+x2^2)

# Define semi algebraic support for the measure 
K = @set (x1^2+x2^2-1)*(x1-3)==0

# Define quadratic polynomial
theta=1 + x1^2 + x2^2 +x3^2 

# Define small parameter 
eps = 1e-5

# Define the order of relaxation
k=2

gmp = GMPModel()

# Add measure to the model
@measure gmp μ [x1,x2,x3] support=K

# Define the degree d
d = 1 + floor(Int,degree(leadingmonomial(f))/2)

# Define the objective 
@objective gmp Min Mom(theta^k*(f+eps*theta^d),μ)

# Constrain μ to be a probablity measure
@constraint gmp Mom(theta^k,μ) == 1

# We solve the relaxation
relax!(gmp, d+k, with_optimizer(Mosek.Optimizer, QUIET=true))

println("Relaxation order: $(k)")
println("Objective value: $(objective_value(gmp))")

elapsed = time() - start
println("elapsed time = ",elapsed)

and receive:

Relaxation order: 2
Objective value: -2.828137730150422e-9
elapsed time = 0.24699997901916504

We know that MomentOpt is based on SumOfSquares. Although MomentOpt (as well as SumOfSquares) is more effective than the optimization model relying only on JuMP, it is hard to be acceptable to provide false results by MomentOpt (as well as SumOfSquares).

[blegat]

The difference might be due to the fact that SumOfSquares computes Groebner basis and computes the reduced form instead of creating multipliers λ. Let

using DynamicPolynomials, SemialgebraicSets

# Define polnomial variables
@polyvar x1 x2 x3

# Polynomial to optimize 
f = (x1^2 + x2^2 - 2) * (x1^2 + x2^2)

# Define semi algebraic support for the measure
h = (x1^2 + x2^2 - 1) * (x1 - 3)
K = @set h == 0

# Define quadratic polynomial
θ = 1 + x1^2 + x2^2 + x3^2 

# Define small parameter 
ϵ = 1e-5

# Define the order of relaxation
k = 2

# Define the degree d
d = 1 + div(degree(leadingmonomial(f)), 2)

p = f + ϵ * θ^d

The MomentOpt problem is

using MomentOpt, JuMP, MosekTools

gmp = GMPModel()

# Add measure to the model
@measure(gmp, μ, [x1,x2,x3], support=K)

# Define the objective 
@objective(gmp, Min, Mom(θ^k * p, μ))

# Constrain μ to be a probablity measure
@constraint(gmp, Mom(θ^k, μ) == 1)

# We solve the relaxation
relax!(gmp, d + k, with_optimizer(Mosek.Optimizer, QUIET=true))
println(gmp.dstatus)
println(objective_value(gmp))

which gives

INFEASIBLE
6.114033382041494e-9

MomentOpt writes the SumOfSquares problem:

using SumOfSquares
model = Model(with_optimizer(Mosek.Optimizer, QUIET=true))
@variable(model, γ)
@objective(model, Max, γ)
@constraint(model, θ^k * (p - γ) in SOSCone(), domain=K)
optimize!(model)
println(termination_status(model))
println(objective_value(model))

which is equivalent to

using SumOfSquares
model = Model(with_optimizer(Mosek.Optimizer, QUIET=true))
@variable(model, γ)
@objective(model, Max, γ)
q = rem(θ^k * (p - γ), ideal(K))
monos = SumOfSquares.monomials_half_newton_polytope(monomials(q), tuple())
@variable(model, s, SOSPoly(monos))
@constraint(model, rem(q - s, ideal(K)) == 0)
optimize!(model)
println(termination_status(model))
println(objective_value(model))

where rem takes the reduction using the Groebner basis of ideal(K).

The approach used by YALMIP should be close to the following

using SumOfSquares
model = Model(with_optimizer(Mosek.Optimizer, QUIET=true))
@variable(model, γ)
@objective(model, Max, γ)
q = θ^k * (p - γ)
gram_monos = monomials(variables(q), 0:div(maxdegree(q) + 1, 2))
@variable(model, gram, SOSPoly(gram_monos))
d = q - gram
λ_monos = monomials(variables(d), 0:maxdegree(d))
@variable(model, λ, Poly(λ_monos))
@constraint(model, d == λ * h)
optimize!(model)
println(termination_status(model))
println(objective_value(model))

which gives

OPTIMAL
-0.9999199975755827

it is hard to be acceptable to provide false results by MomentOpt (as well as SumOfSquares)

The result isn't incorrect, as the termination status is INFEASIBLE, the objective value is computed with the infeasibility certificate returned by Mosek so it cannot be interpreted as a lower bound for the Moment problem. Maybe MomentOpt should print a warning but not printing any warning and requiring the user to always check the status is consistent with what we do in JuMP.

[maihoanganh]

Thank you very much. Your explanation is absolutely correct. I was convinced of using SumOfSquares. Moreover, the following codes show the efficiency of SumOfSquares compared to the model relying only on JuMP:

using DynamicPolynomials, SemialgebraicSets, SumOfSquares, CPUTime, MosekTools

start = time()

# Define polnomial variables
@polyvar x1 x2 x3

# Polynomial to optimize 
f = (x1^2 + x2^2 - 2) * (x1^2 + x2^2)

# Define semi algebraic support for the measure
h = (x1^2 + x2^2 - 1) * (x1 - 3)
K = @set h == 0

# Define quadratic polynomial
θ = 1 + x1^2 + x2^2 + x3^2 

# Define small parameter 
ϵ = 1e-5

# Define the order of relaxation
k = 2

# Define the degree d
d = 1 + div(degree(leadingmonomial(f)), 2)

p = f + ϵ * θ^d

model = Model(with_optimizer(Mosek.Optimizer, QUIET=true))
@variable(model, γ)
@objective(model, Max, γ)
q = θ^k * (p - γ)
gram_monos = monomials(variables(q), 0:div(maxdegree(q) + 1, 2))
@variable(model, gram, SOSPoly(gram_monos))
d = q - gram
λ_monos = monomials(variables(d), 0:maxdegree(d))
@variable(model, λ, Poly(λ_monos))
@constraint(model, d == λ * h)
optimize!(model)
println(termination_status(model))
println(objective_value(model))

elapsed = time() - start
println("elapsed time = ",elapsed)

gives:

OPTIMAL
-0.9999199975677612
elapsed time = 0.25999999046325684

while

using DynamicPolynomials, JuMP, MosekTools, CPUTime

start = time()

# define polnomial variables
@polyvar x1 x2 x3
x=[x1;x2;x3];n=length(x)

# Polynomial to optimize 
f = (x1^2 + x2^2 - 2) * (x1^2 + x2^2)

#inequalities polynomial
g = [] ; m=length(g)

#equalities polynomial
h = [(x1^2 + x2^2 - 1) * (x1 - 3)] ; l=length(h)

# small parameter
eps = 1e-5

# index of relaxation
k = 2

# quadraic polynomial
theta=1+x'*x

# Degree of objective polynomial
if m==0 && l==0
    df = ceil(Int,degree(leadingmonomial(f))/2)
else
    df = floor(Int,degree(leadingmonomial(f))/2) + 1
end

# Degree of inequalities polynomials
dg = []
for i = 1:m
    dg = [dg; ceil(Int,degree(leadingmonomial(g[i]))/2)]
end

# Degree of inequalities polynomials
dh = []
for j = 1:l
    dh = [dh; ceil(Int,degree(leadingmonomial(h[j]))/2)]
end

# Define vetor of monomials
v0= monomials(x, 0)
for j in 1:k+df
    v0= [v0;monomials(x, j)]
end

w0= monomials(x, 0)
for j in 1:2*(k+df)
    w0= [w0;monomials(x, j)]
end

length_v=[];
for i in 1:m
    length_v=[length_v;binomial(k+df-dg[i]+n,n)]
end

length_v_max=length(v0)

length_w=[];
for j in 1:l  
    length_w=[length_w;binomial(2*(k+df-dh[j])+n,n)]
end

# Define sum of square cone
model = Model(with_optimizer(Mosek.Optimizer, QUIET=true))

# Weighted SOS matrix
@variable(model, G0[1:length_v_max, 1:length_v_max],PSD)

# Weighted SOS decomposition
wSOS=v0'*G0*v0

for i in 1:m
     G = @variable(model, [1:length_v[i], 1:length_v[i]],PSD)
    wSOS += g[i]*v0[1:length_v[i]]'*G*v0[1:length_v[i]]
end

for j in 1:l
    q = @variable(model, [1:length_w[j]])
    wSOS += h[j]*w0[1:length_w[j]]'*q
end

@variable(model, lambda)
@constraint(model, coefficients(theta^k*(f-lambda + eps*theta^df) - wSOS) .== 0)
@objective(model, Max, lambda)

optimize!(model)

opt_val = value(lambda)
println("termination status = ", termination_status(model))
println("opt_val = ",opt_val)

elapsed = time() - start
println("elapsed time = ",elapsed)

gives:

termination status = OPTIMAL
opt_val = -0.9999200025317224
elapsed time = 1.186999797821045

[blegat]

Note that you can still use the reduced form with Groebner basis instead of creating multipliers to reduce solve time:

using SumOfSquares
model = Model(with_optimizer(Mosek.Optimizer))
@variable(model, γ)
@objective(model, Max, γ)
q = θ^k * (p - γ)
gram_monos = monomials(variables(q), 0:div(maxdegree(q) + 1, 2))
@variable(model, gram, SOSPoly(gram_monos))
@constraint(model, q == gram, domain = K)
optimize!(model)
println(termination_status(model))
println(objective_value(model))
println(MOI.get(model, MOI.SolveTime()))

gives

Problem
  Name                   :                 
  Objective sense        : max             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 166             
  Cones                  : 0               
  Scalar variables       : 1               
  Matrix variables       : 1               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Problem
  Name                   :                 
  Objective sense        : max             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 166             
  Cones                  : 0               
  Scalar variables       : 1               
  Matrix variables       : 1               
  Integer variables      : 0               

Optimizer  - threads                : 4               
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 166
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 2                 conic                  : 2               
Optimizer  - Semi-definite variables: 1                 scalarized             : 1596            
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 1.39e+04          after factor           : 1.39e+04        
Factor     - dense dim.             : 0                 flops                  : 1.09e+07        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   1.0e+00  1.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  0.00  
1   3.7e-01  3.7e-01  7.0e-02  4.68e-01   -9.500526849e-01  -7.219866307e-01  3.7e-01  0.01  
2   7.8e-02  7.8e-02  8.3e-03  2.80e+00   -1.006874189e+00  -9.945732188e-01  7.8e-02  0.01  
3   1.4e-02  1.4e-02  3.7e-04  1.72e+00   -1.003180151e+00  -1.000938091e+00  1.4e-02  0.02  
4   1.6e-03  1.6e-03  9.1e-06  1.52e+00   -1.000021740e+00  -9.998175813e-01  1.6e-03  0.02  
5   3.2e-04  3.2e-04  6.6e-07  1.31e+00   -9.998846320e-01  -9.998486390e-01  3.2e-04  0.03  
6   6.5e-05  6.5e-05  4.8e-08  1.39e+00   -9.999189547e-01  -9.999129053e-01  6.4e-05  0.03  
7   1.0e-05  1.0e-05  2.4e-09  1.14e+00   -9.999199692e-01  -9.999190684e-01  9.2e-06  0.04  
8   3.9e-07  3.9e-07  1.8e-11  9.89e-01   -9.999199987e-01  -9.999199648e-01  3.5e-07  0.04  
9   6.1e-09  9.7e-10  2.2e-15  1.00e+00   -9.999200000e-01  -9.999199999e-01  8.6e-10  0.04  
Optimizer terminated. Time: 0.04    

OPTIMAL
-0.999920000016737
0.044381141662597656

While with a multiplier you get 536 constraints, 287 scalar variables and a solve time of 0.0731 seconds.

[maihoanganh]

By the way, please show me how to get the Gram matrix of gram and the corresponding moment matrix (the dual solution of Gram matrix) in your latter example.

[blegat]

Assuming the code is

using SumOfSquares, MosekTools
model = Model(with_optimizer(Mosek.Optimizer, QUIET=true))
@variable(model, γ)
@objective(model, Max, γ)
q = θ^k * (p - γ)
gram_monos = monomials(variables(q), 0:div(maxdegree(q) + 1, 2))
@variable(model, gram, SOSPoly(gram_monos))
c = @constraint(model, q == gram, domain = K)
optimize!(model)
println(termination_status(model))
println(objective_value(model))
println(MOI.get(model, MOI.SolveTime()))

how to get the Gram matrix of gram and the corresponding moment matrix

value(gram) and moment_matrix(dual(c), gram_monos).