Non-commutative optimization under constraints fails

[araujoms]

The following code gives me a stack overflow error:

using DynamicPolynomials
using SumOfSquares
import SCS

model = Model(SCS.Optimizer)

@ncpolyvar a0 a1 b0 b1
p = a0*b0 + a0*b1 + a1*b0 - a1*b1
@variable(model, λ)
@objective(model, Min, λ)
S = @set a0^2 == 1 && a1^2 == 1 && b0^2 == 1 && b1^2 == 1
con_ref = @constraint(model, λ - p in SOSCone(), domain = S)

optimize!(model)

Here I'm optimizing the CHSH functional, the simplest problem in nonlocality. The documentation is not clear whether this is supposed to work. If I change @ncpolyvar to @polyvar then it works, and the answer is 2*sqrt(2) as expected.

[odow]

A more minimal example is

julia> using DynamicPolynomials

julia> using SumOfSquares

julia> import SCS

julia> model = Model(SCS.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: SCS

julia> @ncpolyvar a
(a,)

julia> @constraint(model, a in SOSCone(), domain = @set a^2 == 1)
(1)a is SOS

julia> optimize!(model)
ERROR: StackOverflowError:
Stacktrace:
 [1] mapexponents(f::Function, m1::Monomial{false}, m2::Monomial{false}) (repeats 79984 times)
   @ MultivariatePolynomials ~/.julia/packages/MultivariatePolynomials/sWAOE/src/monomial.jl:128

[blegat]

Thanks for reporting this, the noncommutative part does not have enough tests unfortunately, especially with constraints

[sebastiendesignolle]

I'm running into the exact same problem, and I fail to understand the logic behind map_exponents to try to fix things. In particular, it's not clear where the function targeted by map_exponents(f, monomial(m1), monomial(m2)) in monomial.jl:147 is supposed to be. For the moment, this calls itself indefinitely, hence the error. Any thoughts?

Also, I've noticed that the minimal example above is now throwing a different error as it's trying to divide polynomials with non-commutative variables. Here is the beginning of this error.

ERROR: Not implemented yet
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] divides(m1::DynamicPolynomials.Monomial{…}, m2::DynamicPolynomials.Monomial{…})
    @ DynamicPolynomials ~/.julia/packages/DynamicPolynomials/UsmGj/src/div.jl:2
  [3] divides(t1::DynamicPolynomials.Monomial{…}, t2::Term{…})
    @ MultivariatePolynomials ~/.julia/packages/MultivariatePolynomials/iTfZE/src/division.jl:25
  [4] divrem(f::Polynomial{…}, g::Vector{…}; kwargs::@Kwargs{})
    @ MultivariatePolynomials ~/.julia/packages/MultivariatePolynomials/iTfZE/src/division.jl:511
  [5] divrem
    @ ~/.julia/packages/MultivariatePolynomials/iTfZE/src/division.jl:485 [inlined]
  [6] #rem#115
    @ ~/.julia/packages/MultivariatePolynomials/iTfZE/src/division.jl:196 [inlined]
  [7] rem
    @ ~/.julia/packages/MultivariatePolynomials/iTfZE/src/division.jl:195 [inlined]
  [8] rem
    @ ~/.julia/packages/SemialgebraicSets/BoKut/src/ideal.jl:60 [inlined]
  [9] bridge_constraint(::Type{…}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, f::MathOptInterface.VectorAffineFunction{…}, s::PolyJuMP.ZeroPolynomialSet{…})
    @ PolyJuMP.Bridges.Constraint ~/.julia/packages/PolyJuMP/aU72s/src/Bridges/Constraint/zero_polynomial_in_algebraic_set.jl:31
[...]

[blegat]

In the commutative case, several operations can fall back to mapexponents. For instance, the product of monomials is mapexponents(+, m1, m2), the division is mapexponents(-, m1, m2), the gcd is mapexponents(min, m1, m2), etc...
For the noncommutative case, you can't really fall back to mapexponents or I agree with you that it's not so clear what mapexponents actually means in that case.

I refactored MultivariatePolynomials quite a lot so that if you just define a new type of monomial, you can easily make it work with SumOfSquares. See for instance https://github.com/blegat/CondensedMatterSOS.jl/

On SumOfSquares master, we refactored things to work with arbitrary basis defining the StarAlgebras interface and I showcased an example in https://youtu.be/CGPHaHxCG2w so in the future, you won't even have to implement the MultivariatePolynomials API

So these are other options than DynamicPolynomials.@ncpolyvar.

In any case, the issue here is that we trying to take the remainder modulo the ideal (which doesn't work in the noncommutative case) instead of just instantiating multipliers which is something I'd like to fix soon. In the meantime, you can create the multipliers for the equalities of the ideal manually with @variable(model, ..., Poly(...))

[sebastiendesignolle]

Thanks a lot, I'll try to use this!