Refactor Newton postfilter (#287)

* Refactor Newton postfilter

* Fix format

* Fixes

* Fix tests

* Fix test

* Fix format

Author: blegat

docs/src/tutorials/Getting started/sum-of-squares_matrices.jl

@@ -69,13 +69,15 @@ p = vec(y)' * P * vec(y)
 # without exploiting this multipartite structure gives the following 6 monomials.
 
 X = monomials(p)
-@test Certificate.monomials_half_newton_polytope(X, tuple(), apply_post_filter = false) == [x * y[1], x * y[2], y[1] * y[2], x, y[1], y[2]] #src
-Certificate.monomials_half_newton_polytope(X, tuple(), apply_post_filter = false) #!jl
+unipartite = Certificate.NewtonDegreeBounds(tuple())
+@test Certificate.monomials_half_newton_polytope(X, unipartite) == [x * y[1], x * y[2], y[1] * y[2], x, y[1], y[2]] #src
+Certificate.monomials_half_newton_polytope(X, unipartite) #!jl
 
 # Exploiting the multipartite structure gives 4 monomials.
 
-@test Certificate.monomials_half_newton_polytope(X, ([x], y), apply_post_filter = false) == [x * y[1], x * y[2], y[1], y[2]] #src
-Certificate.monomials_half_newton_polytope(X, ([x], y), apply_post_filter = false) #!jl
+multipartite = Certificate.NewtonDegreeBounds(([x], y))
+@test Certificate.monomials_half_newton_polytope(X, multipartite) == [x * y[1], x * y[2], y[1], y[2]] #src
+Certificate.monomials_half_newton_polytope(X, multipartite) #!jl
 
 # In the example above, there were only 3 monomials, where does the difference come from ?
 # Using the monomial basis, the only product of two monomials that is equal to
@@ -84,12 +86,12 @@ Certificate.monomials_half_newton_polytope(X, ([x], y), apply_post_filter = fals
 # hence the whole column and row will be zero as well.
 # Therefore, we can remove this monomial.
 
-@test Certificate.monomials_half_newton_polytope(X, ([x], y)) == [x * y[1], x * y[2], y[1]] #src
-Certificate.monomials_half_newton_polytope(X, ([x], y)) #!jl
+@test Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(multipartite)) == [x * y[1], x * y[2], y[1]] #src
+Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(multipartite)) #!jl
 
 # The same reasoning can be used for monomials `y[1]y[2]` and `x` therefore whether
 # we exploit the multipartite structure or not, we get only 3 monomials thanks
 # to this post filter.
 
-@test Certificate.monomials_half_newton_polytope(X, tuple()) == [x * y[1], x * y[2], y[1]] #src
-Certificate.monomials_half_newton_polytope(X, tuple()) #!jl
+@test Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(unipartite)) == [x * y[1], x * y[2], y[1]] #src
+Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(unipartite)) #!jl

src/Certificate/ideal.jl

@@ -96,19 +96,25 @@ with all variables in the same part.
 struct Newton{
     CT<:SumOfSquares.SOSLikeCone,
     BT<:MB.AbstractPolynomialBasis,
-    NPT<:Tuple,
+    N<:AbstractNewtonPolytopeApproximation,
 } <: SimpleIdealCertificate{CT,BT}
     cone::CT
     basis::Type{BT}
-    variable_groups::NPT
+    newton::N
 end
+
+function Newton(cone, basis, variable_groups::Tuple)
+    return Newton(
+        cone,
+        basis,
+        NewtonFilter(NewtonDegreeBounds(variable_groups)),
+    )
+end
+
 function gram_basis(certificate::Newton{CT,B}, poly) where {CT,B}
     return MB.basis_covering_monomials(
         B,
-        monomials_half_newton_polytope(
-            MP.monomials(poly),
-            certificate.variable_groups,
-        ),
+        monomials_half_newton_polytope(MP.monomials(poly), certificate.newton),
     )
 end
 function gram_basis_type(::Type{<:Newton{CT,BT}}) where {CT,BT}

src/Certificate/newton_polytope.jl

@@ -73,21 +73,26 @@ function is_commutative(vars)
     return length(vars) < 2 || prod(vars[1:2]) == prod(reverse(vars[1:2]))
 end
 
+abstract type AbstractNewtonPolytopeApproximation end
+
+# `maxdegree` and `mindegree` on each variable separately
+# and on each group of variable.
+struct NewtonDegreeBounds{NPT} <: AbstractNewtonPolytopeApproximation
+    variable_groups::NPT
+end
+
 # Multipartite
 # TODO we might do this recursively : do 2 parts, merge them, merge with next
 #      one and so on so that the filter at the end prunes more.
-function half_newton_polytope(
-    X::AbstractVector,
-    parts::Tuple;
-    apply_post_filter = true,
-)
+function half_newton_polytope(X::AbstractVector, newton::NewtonDegreeBounds)
     if !is_commutative(MP.variables(X))
         throw(
             ArgumentError(
                 "Multipartite Newton polytope not supported with noncommutative variables.",
             ),
         )
     end
+    parts = newton.variable_groups
     if !all(
         i -> all(j -> i == j || isempty(parts[i] ∩ parts[j]), 1:length(parts)),
         1:length(parts),
@@ -108,25 +113,23 @@ function half_newton_polytope(
     end
     if length(all_parts) == 1
         # all variables on same part, fallback to shortcut
-        return half_newton_polytope(
-            X,
-            tuple();
-            apply_post_filter = apply_post_filter,
-        )
+        return half_newton_polytope(X, NewtonDegreeBounds(tuple()))
     end
     monovecs = map(vars -> sub_half_newton_polytope(X, vars), all_parts)
     # Cartesian product of the newton polytopes of the different parts
     product = [prod(monos) for monos in Iterators.product(monovecs...)]
     mindeg, maxdeg = cfld(MP.extdegree(X), 2)
     # We know that the degree inequalities are satisfied variable-wise and
     # part-wise but for all variables together so we filter with that
-    monos =
-        MP.monovec(filter(mono -> mindeg <= MP.degree(mono) <= maxdeg, product))
-    if apply_post_filter
-        return post_filter(monos, X)
-    else
-        return monos
-    end
+    gram_monos = filter(mono -> mindeg <= MP.degree(mono) <= maxdeg, product)
+    return MP.monovec(gram_monos)
+end
+
+# Filters out points ouside the Newton polytope from the
+# outer approximation given by `outer_approximation`.
+struct NewtonFilter{N<:AbstractNewtonPolytopeApproximation} <:
+       AbstractNewtonPolytopeApproximation
+    outer_approximation::N
 end
 
 function __chip(cur, i, vars, exps, n, op)
@@ -162,16 +165,12 @@ end
 
 # Shortcut for more efficient `extdeg` and `exp` function in case all the
 # variables are in the same part
-function half_newton_polytope(
-    X::AbstractVector,
-    ::Tuple{};
-    apply_post_filter = true,
-)
+function half_newton_polytope(X::AbstractVector, ::NewtonDegreeBounds{Tuple{}})
     vars = MP.variables(X)
     if is_commutative(vars)
         # Commutative variables
         exp(mono, i) = MP.exponents(mono)[i]
-        monos = _sub_half_newton_polytope(X, MP.extdegree(X), exp, vars)
+        return _sub_half_newton_polytope(X, MP.extdegree(X), exp, vars)
     else
         # Non-commutative variables
         # We use Newton chip method of [Section 2.3, BKP16].
@@ -201,15 +200,15 @@ function half_newton_polytope(
                 end
             end
         end
-        monos = MP.monovec(_monos)
-    end
-    if apply_post_filter
-        return post_filter(monos, X)
-    else
-        return monos
+        return MP.monovec(_monos)
     end
 end
 
+function half_newton_polytope(monos::AbstractVector, newton::NewtonFilter)
+    gram_monos = half_newton_polytope(monos, newton.outer_approximation)
+    return post_filter(gram_monos, monos)
+end
+
 # If `mono` is such that there is no other way to have `mono^2` by multiplying
 # two different monomials of `monos` and `mono` is not in `X` then, the corresponding
 # diagonal entry of the Gram matrix will be zero hence the whole column and row
@@ -271,15 +270,10 @@ function post_filter(monos, X)
 end
 
 function monomials_half_newton_polytope(
-    X::AbstractVector,
-    parts;
-    apply_post_filter = true,
+    monos::AbstractVector,
+    newton::AbstractNewtonPolytopeApproximation,
 )
-    return half_newton_polytope(
-        MP.monovec(X),
-        parts;
-        apply_post_filter = apply_post_filter,
-    )
+    return half_newton_polytope(MP.monovec(monos), newton)
 end
 
 struct DegreeBounds{M}
@@ -289,6 +283,27 @@ struct DegreeBounds{M}
     variablewise_maxdegree::M
 end
 
+# TODO add to MP
+function _map_powers(f, mono)
+    exps = map(f, MP.powers(mono))
+    return _monomial(MP.variables(mono), exps)
+end
+function _map_exponents(f, mono)
+    exps = map(f, MP.exponents(mono))
+    return _monomial(MP.variables(mono), exps)
+end
+
+_min_half(d::Integer) = cld(d, 2)
+_max_half(d::Integer) = fld(d, 2)
+function _half(d::DegreeBounds)
+    return DegreeBounds(
+        _min_half(d.mindegree),
+        _max_half(d.maxdegree),
+        _map_exponents(_min_half, d.variablewise_mindegree),
+        _map_exponents(_max_half, d.variablewise_maxdegree),
+    )
+end
+
 function min_degree(p, v)
     return mapreduce(Base.Fix2(MP.degree, v), min, MP.monomials(p))
 end
@@ -306,35 +321,30 @@ function max_degree(p, v)
 end
 
 function min_shift(d, shift)
-    return max(0, cld(d - shift, 2))
-end
-function max_shift(d, shift)
-    return fld(d - shift, 2)
+    return max(0, d - shift)
 end
 
 function minus_shift(
     deg,
-    m::MP.AbstractMonomial,
+    mono::MP.AbstractMonomial,
     p::MP.AbstractPolynomialLike,
     shift,
 )
-    exps = map(MP.powers(m)) do (v, d)
+    return _map_powers(mono) do (v, d)
         return shift(d, deg(p, v))
     end
-    return _monomial(MP.variables(m), exps)
 end
 
 function minus_shift(d::DegreeBounds, p::MP.AbstractPolynomialLike)
     var_mindegree =
         minus_shift(min_degree, d.variablewise_mindegree, p, min_shift)
-    var_maxdegree =
-        minus_shift(max_degree, d.variablewise_maxdegree, p, max_shift)
+    var_maxdegree = minus_shift(max_degree, d.variablewise_maxdegree, p, -)
     if isnothing(var_maxdegree)
         return
     end
     return DegreeBounds(
         min_shift(d.mindegree, MP.mindegree(p)),
-        max_shift(d.maxdegree, MP.maxdegree(p)),
+        d.maxdegree - MP.maxdegree(p),
         var_mindegree,
         var_maxdegree,
     )
@@ -436,8 +446,8 @@ function putinar_degree_bounds(
 )
     mindegree = 0
     # TODO homogeneous case
-    mindeg(g) = min_shift(mindegree, MP.mindegree(g))
-    maxdeg(g) = max_shift(maxdegree, MP.maxdegree(g))
+    mindeg(g) = _min_half(min_shift(mindegree, MP.mindegree(g)))
+    maxdeg(g) = _max_half(maxdegree - MP.maxdegree(g))
     degrange(g) = mindeg(g):maxdeg(g)
     minus_degrange(g) = -maxdeg(g):-mindeg(g)
     # The multiplier will have degree `0:2fld(maxdegree - MP.maxdegree(g), 2)`
@@ -477,3 +487,52 @@ function putinar_degree_bounds(
         _monomial(vars, vars_maxdeg),
     )
 end
+
+function multiplier_basis(g::MP.AbstractPolynomialLike, bounds::DegreeBounds)
+    shifted = minus_shift(bounds, g)
+    if isnothing(shifted)
+        halved = nothing
+    else
+        halved = _half(shifted)
+    end
+    basis = MB.MonomialBasis
+    if isnothing(halved)
+        # TODO add `MB.empty_basis` to API
+        return MB.maxdegree_basis(
+            basis,
+            MP.variables(bounds.variablewise_mindegree),
+            -1,
+        )
+    else
+        return maxdegree_gram_basis(basis, halved)
+    end
+end
+
+function half_newton_polytope(
+    p::MP.AbstractPolynomialLike,
+    gs::AbstractVector{<:MP.AbstractPolynomialLike},
+    vars,
+    maxdegree,
+    ::NewtonDegreeBounds,
+)
+    # TODO take `variable_groups` into account
+    bounds = putinar_degree_bounds(p, gs, vars, maxdegree)
+    return [multiplier_basis(g, bounds) for g in gs]
+end
+
+function half_newton_polytope(
+    p::MP.AbstractPolynomialLike,
+    gs::AbstractVector{<:MP.AbstractPolynomialLike},
+    vars,
+    maxdegree,
+    newton::NewtonFilter{<:NewtonDegreeBounds},
+)
+    # TODO
+    return half_newton_polytope(
+        p,
+        gs,
+        vars,
+        maxdegree,
+        newton.outer_approximation,
+    )
+end