Sampling basis

docs/src/tutorials/Getting started/sampling.jl

@@ -0,0 +1,23 @@
+# # Sampling basis
+
+#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Getting started/sampling.ipynb)
+#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Getting started/sampling.ipynb)
+# **Contributed by**: Benoît Legat
+
+using Test #src
+using DynamicPolynomials
+using SumOfSquares
+import Hypatia
+
+# In this tutorial, we show how to use a different polynomial basis
+# for enforcing the equality between the polynomial and its Sum-of-Squares decomposition.
+
+@polyvar x
+p = x^4 - 4x^3 - 2x^2 + 12x + 9
+
+model = Model(Hypatia.Optimizer)
+set_silent(model)
+@constraint(model, p in SOSCone(), zero_basis = BoxSampling([-1], [1]))
+optimize!(model)
+solution_summary(model)
+backend(model)

src/Certificate/ideal.jl

@@ -40,12 +40,22 @@ function _combine_with_gram(
     return MB.SubBasis{B}(keys(SA.coeffs(p)))
 end
 
-_reduce_with_domain(basis::MB.SubBasis, ::FullSpace) = basis
+function _reduce_with_domain(basis::MB.SubBasis, zero_basis, ::FullSpace)
+    return MB.explicit_basis_covering(zero_basis, basis)
+end
 
-function _reduce_with_domain(basis::MB.SubBasis{B}, domain) where {B}
-    if B !== MB.Monomial
-        error("Only Monomial basis support with an equalities in domain")
-    end
+function _reduce_with_domain(basis, zero_basis, domain)
+    return __reduce_with_domain(basis, zero_basis, domain)
+end
+
+function __reduce_with_domain(_, _, _)
+    return error("Only Monomial basis support with an equalities in domain")
+end
+function __reduce_with_domain(
+    basis::MB.SubBasis{MB.Monomial},
+    ::MB.FullBasis{MB.Monomial},
+    domain,
+)
     I = ideal(domain)
     # set of standard monomials that are hit
     standard = Set{eltype(basis.monomials)}()
@@ -60,14 +70,15 @@ function _reduce_with_domain(basis::MB.SubBasis{B}, domain) where {B}
 end
 
 function zero_basis(
-    ::AbstractIdealCertificate,
+    cert::AbstractIdealCertificate,
     basis,
     domain,
     gram_bases,
     weights,
 )
     return _reduce_with_domain(
         _combine_with_gram(basis, gram_bases, weights),
+        _zero_basis(cert),
         domain,
     )
 end
@@ -270,15 +281,27 @@ function _quotient_basis_type(
     }
 end
 
+_zero_basis(c::SimpleIdealCertificate) = c.zero_basis
+
+function _zero_basis_type(::Type{<:SimpleIdealCertificate{C,G,Z}}) where {C,G,Z}
+    return Z
+end
+
+_zero_basis(c::Remainder) = _zero_basis(c.gram_certificate)
+
+function _zero_basis_type(::Type{Remainder{C}}) where {C}
+    return _zero_basis_type(C)
+end
+
 function MA.promote_operation(
     ::typeof(zero_basis),
-    ::Type{<:Union{SimpleIdealCertificate,Remainder}},
-    ::Type{B},
+    C::Type{<:Union{SimpleIdealCertificate,Remainder}},
+    ::Type,
     ::Type{SemialgebraicSets.FullSpace},
     ::Type,
     ::Type,
-) where {B}
-    return B
+)
+    return MB.explicit_basis_type(_zero_basis_type(C))
 end
 
 function MA.promote_operation(

src/constraints.jl

@@ -438,6 +438,11 @@ function _default_basis(p::MP.AbstractPolynomialLike, basis)
     )
 end
 
+_default_zero_basis(basis, ::Nothing) = MB.implicit_basis(basis)
+function _default_zero_basis(basis, nodes::MB.AbstractNodes)
+    return MB.ImplicitLagrangeBasis(MP.variables(basis), nodes)
+end
+
 function JuMP.build_constraint(
     _error::Function,
     p,
@@ -446,11 +451,14 @@ function JuMP.build_constraint(
     zero_basis = nothing,
     kws...,
 )
-    __coefs, basis, gram_basis = _default_basis(p, basis)
-    if isnothing(zero_basis)
-        zero_basis = MB.implicit_basis(basis)
-    end
-    set = JuMP.moi_set(cone, basis, gram_basis, zero_basis; kws...)
+    __coefs, basis, gram_basis = _default_gram_basis(p, basis)
+    set = JuMP.moi_set(
+        cone,
+        basis,
+        gram_basis,
+        _default_zero_basis(basis, zero_basis);
+        kws...,
+    )
     _coefs = PolyJuMP.non_constant(__coefs)
     # If a polynomial with real coefficients is used with the Hermitian SOS
     # cone, we want to promote the coefficients to complex