Fix Sparse.Putinar with Newton polytope (#278)

* Fix Sparse.Putinar with Newton polytope

* Fixes

Author: blegat

src/Certificate/Sparsity/monomial.jl

@@ -152,18 +152,39 @@ _monos(basis::MB.MonomialBasis) = basis.monomials
 function _gram_monos(vars, certificate::SumOfSquares.Certificate.MaxDegree{CT, MB.MonomialBasis}) where CT
     return _monos(SumOfSquares.Certificate.maxdegree_gram_basis(MB.MonomialBasis, vars, certificate.maxdegree))
 end
+# poly = s0 + sum si gi
+# where `s1` are the multipliers with basis `multiplier_generator_monos`
+# we want to get the gram basis of `s0`
+function _ideal_monos(poly_monos, multiplier_gram_monos)
+    monos_set = Set(poly_monos)
+    for (gram_monos, gen_monos) in multiplier_gram_monos
+        for a in gram_monos
+            for b in gram_monos
+                for m in gen_monos
+                    push!(monos_set, a * b * m)
+                end
+            end
+        end
+    end
+    return MP.monovec(collect(monos_set))
+end
+# The ideal certificate should only ask for `MP.monomial`
+struct DummyPolynomial{M}
+    monomials::M
+end
+MP.monomials(p::DummyPolynomial) = p.monomials
+MP.variables(p::DummyPolynomial) = MP.variables(p.monomials)
 function sparsity(poly::MP.AbstractPolynomial, domain::SemialgebraicSets.BasicSemialgebraicSet, sp::Monomial, certificate::SumOfSquares.Certificate.AbstractPreorderCertificate)
-    gram_monos = _gram_monos(
-        reduce((v, q) -> unique!(sort!([v..., MP.variables(q)...], rev=true)),
-                  domain.p, init = MP.variables(poly)),
-        SumOfSquares.Certificate.ideal_certificate(certificate)
-    )
     processed = SumOfSquares.Certificate.preprocessed_domain(certificate, domain, poly)
     multiplier_generator_monos = [
         (_monos(SumOfSquares.Certificate.multiplier_basis(certificate, index, processed)),
          MP.monomials(SumOfSquares.Certificate.generator(certificate, index, processed)))
         for index in SumOfSquares.Certificate.preorder_indices(certificate, processed)
     ]
+    gram_monos = _monos(SumOfSquares.Certificate.gram_basis(
+        SumOfSquares.Certificate.ideal_certificate(certificate),
+        DummyPolynomial(_ideal_monos(MP.monomials(poly), multiplier_generator_monos)),
+    ))
     cliques, multiplier_cliques = sparsity(MP.monomials(poly), sp, gram_monos, multiplier_generator_monos)
     return MB.MonomialBasis.(cliques), [MB.MonomialBasis.(clique) for clique in multiplier_cliques]
 end

src/constraint.jl

@@ -81,15 +81,15 @@ function default_ideal_certificate(domain::FullSpace, basis, cone, maxdegree, ne
 end
 
 function default_ideal_certificate(
-    domain::AbstractAlgebraicSet, sparsity::Certificate.Sparsity.NoPattern, basis::AbstractPolynomialBasis, cone, args...)
+    ::AbstractAlgebraicSet, ::Certificate.Sparsity.NoPattern, basis::AbstractPolynomialBasis, cone, args...)
     return Certificate.FixedBasis(cone, basis)
 end
 function default_ideal_certificate(
-    domain::AbstractAlgebraicSet, sparsity::Certificate.Sparsity.NoPattern, args...)
+    domain::AbstractAlgebraicSet, ::Certificate.Sparsity.NoPattern, args...)
     return default_ideal_certificate(domain, args...)
 end
 function default_ideal_certificate(
-    domain::AbstractAlgebraicSet, sparsity::Certificate.Sparsity.Pattern, args...)
+    ::AbstractAlgebraicSet, sparsity::Certificate.Sparsity.Pattern, args...)
     return Certificate.Sparsity.Ideal(sparsity, args...)
 end
 

test/sparsity.jl

@@ -18,11 +18,6 @@ end
 
 set_monos(bases::Vector{<:MB.MonomialBasis}) = Set([basis.monomials for basis in bases])
 
-function Certificate.Sparsity.sparsity(monos::AbstractVector{<:MP.AbstractMonomial}, domain::SemialgebraicSets.BasicSemialgebraicSet, sp::Sparsity.Monomial, maxdegree, degs)
-    half_monos = Certificate.maxdegree_gram_basis(MB.MonomialBasis, variables, div(maxdegree, 2))
-    P = Set(monos)
-end
-
 """
     wml19()
 
@@ -82,8 +77,11 @@ function wml19()
         ]))
         @test set_monos(Certificate.Sparsity.sparsity(f, SignSymmetry(), certificate)) == expected
     end
-    @testset "Example 5.4" begin
-        preorder_certificate = Certificate.Putinar(Certificate.MaxDegree(SOSCone(), MB.MonomialBasis, 4), SOSCone(), MB.MonomialBasis, 4)
+    @testset "Example 5.4 $(typeof(ideal_certificate))" for ideal_certificate in [
+        Certificate.MaxDegree(SOSCone(), MB.MonomialBasis, 4),
+        Certificate.Newton(SOSCone(), MB.MonomialBasis, tuple()),
+    ]
+        preorder_certificate = Certificate.Putinar(ideal_certificate, SOSCone(), MB.MonomialBasis, 4)
         @polyvar x[1:2]
         f = x[1]^4 + x[2]^4 + x[1] * x[2]
         K = @set 1 - 2x[1]^2 - x[2]^2 >= 0
@@ -168,9 +166,8 @@ function l09()
         ]))
     end
 end
-function square_domain()
-    d = 6
-    preorder_certificate = Certificate.Putinar(Certificate.MaxDegree(SOSCone(), MB.MonomialBasis, 6), SOSCone(), MB.MonomialBasis, 6)
+function square_domain(ideal_certificate, d)
+    preorder_certificate = Certificate.Putinar(ideal_certificate, SOSCone(), MB.MonomialBasis, d)
     @polyvar x y
     f = x^2*y^4 + x^4*y^2 - 3*x^2*y*2 + 1
     K = @set(1 - x^2 >= 0 && 1 - y^2 >= 0)
@@ -238,22 +235,37 @@ function drop_monomials()
             end
             @test set_monos(Certificate.Sparsity.sparsity(f, Sparsity.Monomial(ChordalCompletion(), k, use_all_monomials), certificate)) == expected
         end
-        preorder_certificate = Certificate.Putinar(Certificate.MaxDegree(SOSCone(), MB.MonomialBasis, 4), SOSCone(), MB.MonomialBasis, 3)
-        f = polynomial(x^3)
-        K = @set x >= 0
-        @testset "$k $use_all_monomials" for k in 0:3, use_all_monomials in [false, true]
-            basis, preorder_bases = Certificate.Sparsity.sparsity(f, K, Sparsity.Monomial(ChordalCompletion(), k, use_all_monomials), preorder_certificate)
-            if k == 1 && !use_all_monomials
-                @test set_monos(basis) == Set(monovec.([[x^2, x]]))
-            elseif (k == 2 && !use_all_monomials) || (k == 1 && use_all_monomials)
-                @test set_monos(basis) == Set(monovec.([[x^2, 1], [x^2, x]]))
-            else
-                @test set_monos(basis) == Set(monovec.([[x^2, x, 1]]))
-            end
-            if k == 1 && !use_all_monomials
-                @test set_monos(preorder_bases[1]) == Set(monovec.([[x]]))
-            else
-                @test set_monos(preorder_bases[1]) == Set(monovec.([[x, 1]]))
+        @testset "$(typeof(ideal_certificate))" for ideal_certificate in [
+            Certificate.MaxDegree(SOSCone(), MB.MonomialBasis, 4),
+            Certificate.Newton(SOSCone(), MB.MonomialBasis, tuple()),
+        ]
+            preorder_certificate = Certificate.Putinar(ideal_certificate, SOSCone(), MB.MonomialBasis, 3)
+            f = polynomial(x^3)
+            K = @set x >= 0
+            @testset "$k $use_all_monomials" for k in 0:3, use_all_monomials in [false, true]
+                basis, preorder_bases = Certificate.Sparsity.sparsity(f, K, Sparsity.Monomial(ChordalCompletion(), k, use_all_monomials), preorder_certificate)
+                if ideal_certificate isa Certificate.Newton
+                    if use_all_monomials
+                        @test set_monos(basis) == Set(monovec.([[x]]))
+                        @test set_monos(preorder_bases[1]) == Set(monovec.([[x, 1]]))
+                    else
+                        @test isempty(basis)
+                        @test set_monos(preorder_bases[1]) == Set(monovec.([[x]]))
+                    end
+                else
+                    if k == 1 && !use_all_monomials
+                        @test set_monos(basis) == Set(monovec.([[x^2, x]]))
+                    elseif (k == 2 && !use_all_monomials) || (k == 1 && use_all_monomials)
+                        @test set_monos(basis) == Set(monovec.([[x^2, 1], [x^2, x]]))
+                    else
+                        @test set_monos(basis) == Set(monovec.([[x^2, x, 1]]))
+                    end
+                    if (k == 1 && !use_all_monomials)
+                        @test set_monos(preorder_bases[1]) == Set(monovec.([[x]]))
+                    else
+                        @test set_monos(preorder_bases[1]) == Set(monovec.([[x, 1]]))
+                    end
+                end
             end
         end
     end
@@ -262,7 +274,8 @@ end
     xor_complement_test()
     wml19()
     l09()
-    square_domain()
+    square_domain(Certificate.MaxDegree(SOSCone(), MB.MonomialBasis, 6), 6)
+    square_domain(Certificate.Newton(SOSCone(), MB.MonomialBasis, tuple()), 6)
     sum_square(8)
     @test Certificate.Sparsity.appropriate_type(32) == Int64
     sum_square(32)