jump-dev / SumOfSquares.jl

Sum of Squares Programming for Julia
Other
114 stars 24 forks source link

Empty ideal #107

Open chriscoey opened 5 years ago

chriscoey commented 5 years ago

running this code (replace Hypatia with your favorite SDP solver; sorry I could not figure out how to make the example more minimal) on Julia 1.3 master and JuMP 0.19, I get the following error, which seems to come from DynamicPolynomials due to an empty set of monomials getting passed from _sub_half_newton_polytope:

ERROR: Error while loading expression starting at /home/coey/.julia/dev/HermitianWSOS/scratch.jl:61
caused by [exception 1]
BoundsError: attempt to access 0-element Array{Int64,1} at index [1]
Stacktrace:
 [1] setindex! at ./array.jl:766 [inlined]
 [2] fillZfordeg!(::Array{Array{Int64,1},1}, ::Int64, ::Int64, ::Type{Val{true}}, ::getfield(DynamicPolynomials, Symbol("##12#14")){Array{PolyVar{true},1},getfield(SumOfSquares, Symbol("##13#16")){getfield(SumOfSquares, Symbol("##28#29")),Int64,Array{Int64,1},Array{Int64,1}}}) at /home/coey/.julia/packages/DynamicPolynomials/NSiPE/src/monovec.jl:76
 [3] getZfordegs(::Int64, ::UnitRange{Int64}, ::Type{Val{true}}, ::Function) at /home/coey/.julia/packages/DynamicPolynomials/NSiPE/src/monovec.jl:119
 [4] Type at /home/coey/.julia/packages/DynamicPolynomials/NSiPE/src/monovec.jl:126 [inlined]
 [5] monomials at /home/coey/.julia/packages/DynamicPolynomials/NSiPE/src/monovec.jl:140 [inlined]
 [6] _sub_half_newton_polytope(::MonomialVector{true}, ::Tuple{Int64,Int64}, ::getfield(SumOfSquares, Symbol("##28#29")), ::Array{PolyVar{true},1}) at /home/coey/.julia/packages/SumOfSquares/x4OcB/src/certificate.jl:44
 [7] half_newton_polytope at /home/coey/.julia/packages/SumOfSquares/x4OcB/src/certificate.jl:87 [inlined]
 [8] monomials_half_newton_polytope at /home/coey/.julia/packages/SumOfSquares/x4OcB/src/certificate.jl:93 [inlined]
 [9] SumOfSquares.SOSPolynomialBridge{Float64,MathOptInterface.VectorAffineFunction{Float64},AlgebraicSet{Float64,Polynomial{true,Float64},Buchberger,SemialgebraicSets.SolverUsingMultiplicationMatrices{SemialgebraicSets.GröbnerBasisMultiplicationMatricesAlgorithm,ReorderedSchurMultiplicationMatricesSolver{Float64,Random.MersenneTwister}}},Union{SumOfSquares.GenericVariableBridge{Float64,MathOptInterface.Nonnegatives}, SumOfSquares.GenericVariableBridge{Float64,MathOptInterface.PositiveSemidefiniteConeTriangle}, SumOfSquares.GenericVariableBridge{Float64,SumOfSquares.EmptyCone}, SumOfSquares.GenericVariableBridge{Float64,SumOfSquares.PositiveSemidefinite2x2ConeTriangle}},MathOptInterface.PositiveSemidefiniteConeTriangle,MonomialBasis,Monomial{true},MonomialVector{true}}(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Hypatia.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::MathOptInterface.VectorAffineFunction{Float64}, ::SumOfSquares.SOSPolynomialSet{AlgebraicSet{Float64,Polynomial{true,Float64},Buchberger,SemialgebraicSets.SolverUsingMultiplicationMatrices{SemialgebraicSets.GröbnerBasisMultiplicationMatricesAlgorithm,ReorderedSchurMultiplicationMatricesSolver{Float64,Random.MersenneTwister}}},NonnegPolyInnerCone{MathOptInterface.PositiveSemidefiniteConeTriangle},MonomialBasis,Monomial{true},MonomialVector{true},Tuple{}}) at /home/coey/.julia/packages/SumOfSquares/x4OcB/src/sos_polynomial_bridge.jl:25
 [10] add_constraint(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Hypatia.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::MathOptInterface.VectorAffineFunction{Float64}, ::SumOfSquares.SOSPolynomialSet{AlgebraicSet{Float64,Polynomial{true,Float64},Buchberger,SemialgebraicSets.SolverUsingMultiplicationMatrices{SemialgebraicSets.GröbnerBasisMultiplicationMatricesAlgorithm,ReorderedSchurMultiplicationMatricesSolver{Float64,Random.MersenneTwister}}},NonnegPolyInnerCone{MathOptInterface.PositiveSemidefiniteConeTriangle},MonomialBasis,Monomial{true},MonomialVector{true},Tuple{}}) at /home/coey/.julia/packages/MathOptInterface/C3lip/src/Bridges/bridgeoptimizer.jl:320
 [11] SumOfSquares.SOSPolynomialInSemialgebraicSetBridge{Float64,MathOptInterface.VectorAffineFunction{Float64},AlgebraicSet{Float64,Polynomial{true,Float64},Buchberger,SemialgebraicSets.SolverUsingMultiplicationMatrices{SemialgebraicSets.GröbnerBasisMultiplicationMatricesAlgorithm,ReorderedSchurMultiplicationMatricesSolver{Float64,Random.MersenneTwister}}},NonnegPolyInnerCone{MathOptInterface.PositiveSemidefiniteConeTriangle},Union{SumOfSquares.GenericVariableBridge{Float64,MathOptInterface.Nonnegatives}, SumOfSquares.GenericVariableBridge{Float64,MathOptInterface.PositiveSemidefiniteConeTriangle}, SumOfSquares.GenericVariableBridge{Float64,SumOfSquares.EmptyCone}, SumOfSquares.GenericVariableBridge{Float64,SumOfSquares.PositiveSemidefinite2x2ConeTriangle}},MonomialBasis,Monomial{true},MonomialVector{true},Tuple{}}(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Hypatia.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::MathOptInterface.VectorAffineFunction{Float64}, ::SumOfSquares.SOSPolynomialSet{BasicSemialgebraicSet{Float64,Polynomial{true,Float64},AlgebraicSet{Float64,Polynomial{true,Float64},Buchberger,SemialgebraicSets.SolverUsingMultiplicationMatrices{SemialgebraicSets.GröbnerBasisMultiplicationMatricesAlgorithm,ReorderedSchurMultiplicationMatricesSolver{Float64,Random.MersenneTwister}}}},NonnegPolyInnerCone{MathOptInterface.PositiveSemidefiniteConeTriangle},MonomialBasis,Monomial{true},MonomialVector{true},Tuple{}}) at /home/coey/.julia/packages/SumOfSquares/x4OcB/src/sos_polynomial_in_semialgebraic_set_bridge.jl:54
 [12] add_constraints(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Hypatia.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::Array{MathOptInterface.VectorAffineFunction{Float64},1}, ::Array{SumOfSquares.SOSPolynomialSet{BasicSemialgebraicSet{Float64,Polynomial{true,Float64},AlgebraicSet{Float64,Polynomial{true,Float64},Buchberger,SemialgebraicSets.SolverUsingMultiplicationMatrices{SemialgebraicSets.GröbnerBasisMultiplicationMatricesAlgorithm,ReorderedSchurMultiplicationMatricesSolver{Float64,Random.MersenneTwister}}}},NonnegPolyInnerCone{MathOptInterface.PositiveSemidefiniteConeTriangle},MonomialBasis,Monomial{true},MonomialVector{true},Tuple{}},1}) at /home/coey/.julia/packages/MathOptInterface/C3lip/src/Bridges/bridgeoptimizer.jl:338
 [13] copyconstraints!(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Hypatia.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}, ::Bool, ::MathOptInterface.Utilities.IndexMap, ::Type{MathOptInterface.VectorAffineFunction{Float64}}, ::Type{SumOfSquares.SOSPolynomialSet{BasicSemialgebraicSet{Float64,Polynomial{true,Float64},AlgebraicSet{Float64,Polynomial{true,Float64},Buchberger,SemialgebraicSets.SolverUsingMultiplicationMatrices{SemialgebraicSets.GröbnerBasisMultiplicationMatricesAlgorithm,ReorderedSchurMultiplicationMatricesSolver{Float64,Random.MersenneTwister}}}},NonnegPolyInnerCone{MathOptInterface.PositiveSemidefiniteConeTriangle},MonomialBasis,Monomial{true},MonomialVector{true},Tuple{}}}) at /home/coey/.julia/packages/MathOptInterface/C3lip/src/Utilities/copy.jl:157
 [14] default_copy_to(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Hypatia.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}, ::Bool) at /home/coey/.julia/packages/MathOptInterface/C3lip/src/Utilities/copy.jl:200
 [15] #automatic_copy_to#61 at /home/coey/.julia/packages/MathOptInterface/C3lip/src/Utilities/copy.jl:15 [inlined]
 [16] #automatic_copy_to at ./none:0 [inlined]
 [17] #copy_to#1 at /home/coey/.julia/packages/MathOptInterface/C3lip/src/Bridges/bridgeoptimizer.jl:91 [inlined]
 [18] (::getfield(MathOptInterface, Symbol("#kw##copy_to")))(::NamedTuple{(:copy_names,),Tuple{Bool}}, ::typeof(MathOptInterface.copy_to), ::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Hypatia.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}) at ./none:0
 [19] attach_optimizer(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}}) at /home/coey/.julia/packages/MathOptInterface/C3lip/src/Utilities/cachingoptimizer.jl:130
 [20] optimize!(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{JuMP._MOIModel{Float64}}}) at /home/coey/.julia/packages/MathOptInterface/C3lip/src/Utilities/cachingoptimizer.jl:166
 [21] #optimize!#77(::Bool, ::Bool, ::typeof(optimize!), ::Model, ::Nothing) at /home/coey/.julia/packages/JuMP/jnmGG/src/optimizer_interface.jl:132
 [22] optimize! at /home/coey/.julia/packages/JuMP/jnmGG/src/optimizer_interface.jl:105 [inlined] (repeats 2 times)
blegat commented 5 years ago

It seems that the algebraic set is empty, it is defined by

-0.961538zr₁² - 0.961538zi₁² == 3.5 && -4.80769zr₁² - 4.80769zi₁² == -3.5

which is rewritten into

1 == 0

When when taking the polynomial modulo the ideal there is no term left.

blegat commented 5 years ago

I agree that the fix would be to fix the newton polytope computation for empty monomial vector. We should both add a test of newton polytope for empty monomial vector and a test with empty algebraic set.

chriscoey commented 5 years ago

thanks @blegat. if the ideal is empty should there be an error or a warning or something? should the WSOS constraint just be relaxed/dropped? what should the newton polytope function return for an empty monomial vector?

blegat commented 5 years ago

It should return an empty monomial vector. Then it will create a 0x0 PSD matrix: https://github.com/JuliaOpt/SumOfSquares.jl/blob/fdb6adbd87e17e777b7f6abdb124bc8ddaef01fb/src/psd2x2.jl#L3 which is bridged into adding no constraint: https://github.com/JuliaOpt/SumOfSquares.jl/blob/fdb6adbd87e17e777b7f6abdb124bc8ddaef01fb/src/empty_bridge.jl#L1-L8 so it is equivalent to ignoring the constraint. It makes sense since every polynomial is nonnegative over an empty set.

chriscoey commented 5 years ago

OK I can have a try at fixing this. Thanks

lkapelevich commented 4 years ago

I was getting the same error (slightly different error trace now) with a single equality constraint, e.g.:

model = SOSModel(with_optimizer(Mosek.Optimizer))
@variable(model, gamma)
@constraint(model, y - gamma >= 0, domain = @set(y == 1))
@objective(model, Max, gamma)
optimize!(model)

is that expected?