Fixes

Author: blegat

src/Certificate/newton_polytope.jl

@@ -209,70 +209,6 @@ function half_newton_polytope(basis::SA.ExplicitBasis, newton::NewtonFilter)
     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
-# will be zero hence we can remove this monomial.
-# See [Proposition 3.7, CLR95], [Theorem 2, L09] or [Section 2.4, BKP16].
-
-# [CLR95] Choi, M. D. and Lam, T. Y. and Reznick, B.
-# *Sum of Squares of Real Polynomials*.
-# Proceedings of Symposia in Pure mathematics (1995)
-#
-# [L09] Lofberg, Johan.
-# *Pre-and post-processing sum-of-squares programs in practice*.
-# IEEE transactions on automatic control 54.5 (2009): 1007-1011.
-#
-# [BKP16] Sabine Burgdorf, Igor Klep, and Janez Povh.
-# *Optimization of polynomials in non-commuting variables*.
-# Berlin: Springer, 2016.
-function post_filter(monos, X)
-    num = Dict(mono => 1 for mono in X)
-    function _increase(mono)
-        return num[mono] = get(num, mono, 0) + 1
-    end
-    function _decrease(mono)
-        value = num[mono] - 1
-        if iszero(value)
-            delete!(num, mono)
-        else
-            num[mono] = value
-        end
-        return iszero(value)
-    end
-    for a in monos
-        for b in monos
-            if a != b
-                _increase(a * b)
-            end
-        end
-    end
-    back = Dict{eltype(monos),Int}()
-    keep = ones(Bool, length(monos))
-    function _delete(i)
-        keep[i] = false
-        a = monos[i]
-        for (j, b) in enumerate(monos)
-            if keep[j] && a != b
-                for w in [a * b, b * a]
-                    if _decrease(w) && haskey(back, w)
-                        _delete(back[w])
-                    end
-                end
-            end
-        end
-    end
-    for (i, mono) in enumerate(monos)
-        w = mono^2
-        if haskey(num, w)
-            back[w] = i
-        else
-            _delete(i)
-        end
-    end
-    return monos[findall(keep)]
-end
-
 struct DegreeBounds{M}
     mindegree::Int
     maxdegree::Int
@@ -305,7 +241,6 @@ function min_degree(p::SA.AlgebraElement, v)
     return mapreduce(Base.Fix2(min_degree, v), min, SA.supp(p))
 end
 min_degree(p::MB.Polynomial{MB.Monomial}, v) = MP.degree(p.monomial, v)
-minus_min_degree(p, v) = -min_degree(p, v)
 
 function _monomial(vars, exps)
     if any(Base.Fix2(isless, 0), exps)
@@ -389,13 +324,21 @@ function _interval(a, b)
     end
 end
 
+function _gram_shift_min(d_g, d_s, g::SA.AlgebraElement)
+    if _is_monomial_basis(typeof(g))
+        return d_g + 2minimum(d_s)
+    else
+        return 0
+    end
+end
+
 function deg_range(deg, p, gs, gram_deg, truncation)
     d_max = min(deg(p), truncation)
     sign = deg_sign(deg, p, d_max)
     for g in gs
         d_g, sign_g = deg_sign(deg, g)
         d_s = gram_deg(g)
-        if isempty(d_s) || d_g + 2minimum(d_s) > truncation
+        if isempty(d_s) || _gram_shift_min(d_g, d_s, g) > truncation
             continue
         end
         d = d_g + 2maximum(d_s)
@@ -429,7 +372,7 @@ end
 
 #    deg_range(deg, p, gs, gram_deg, range)
 #
-# Maximum value of `deg(s_0 = p - sum s_i g_i for g in gs) in range` where
+# Maximum value of `deg(s_0 = p - sum s_i g_i for i in eachindex(gs)) in range` where
 # `s_0, s_i` are SOS and `deg(s_i) <= gram_deg(g)`.
 # Note that `range` should be in increasing order.
 function deg_range(deg, p, gs, gram_deg, range::UnitRange)
@@ -451,6 +394,40 @@ _maxdegree(a::SA.AlgebraElement) = maximum(_maxdegree, SA.supp(a))
 _mindegree(p::MB.Polynomial{MB.Monomial}) = MP.mindegree(p.monomial)
 _maxdegree(p::MB.Polynomial{MB.Monomial}) = MP.maxdegree(p.monomial)
 
+_is_monomial_basis(::Type{<:SA.AlgebraElement{<:MB.Algebra{BT,B}}}) = true
+_is_monomial_basis(::Type{<:SA.AlgebraElement}) = false
+
+# Minimum degree of a gram basis for a gram matrix `s`
+# such that the minimum degree of `s * g` is at least `mindegree`.
+function _multiplier_mindegree(
+    mindegree,
+    g::SA.AlgebraElement,
+)
+    if _is_monomial_basis(typeof(g))
+        return _min_half(min_shift(mindegree, _mindegree(g)))
+    else
+        # For instance, with `Chebyshev` a square already produces
+        # a monomial of degree 0 so let's just give up here
+        # The Chebyshev Newton polytope of the product of polynomials
+        # with Chebyshev Newton polytope `N1` and `N2` is a subset
+        # of `(N1 + N2 - R_+^n) ∩ R_+^n` so we compute
+        # `σ(y, cheby(N1) + cheby(N2))` as `σ(max.(y, 0), N1 + N2)`
+        return 0
+    end
+end
+
+# Maximum degree of a gram basis for a gram matrix `s`
+# such that the maximum degree of `s * g` is at most `maxdegree`.
+# Here, we assume that the degree of `*(a::MB.Polynomial, b::MB.Polynomial)`
+# is bounded by the sum of the degree of `a` and `b` for any basis.
+function _multiplier_maxdegree(maxdegree, g::SA.AlgebraElement)
+    return _max_half(maxdegree - _maxdegree(g))
+end
+
+function _multiplier_deg_range(range, g::SA.AlgebraElement)
+    return _multiplier_mindegree(minimum(range), g):_multiplier_maxdegree(maximum(range), g)
+end
+
 function putinar_degree_bounds(
     p::SA.AlgebraElement,
     gs::AbstractVector{<:SA.AlgebraElement},
@@ -459,13 +436,11 @@ function putinar_degree_bounds(
 )
     mindegree = 0
     # TODO homogeneous case
-    mindeg(g) = _min_half(min_shift(mindegree, _mindegree(g)))
-    maxdeg(g) = _max_half(maxdegree - _maxdegree(g))
-    degrange(g) = mindeg(g):maxdeg(g)
-    minus_degrange(g) = -maxdeg(g):-mindeg(g)
+    degrange(g) = _multiplier_deg_range(mindegree:maxdegree, g)
+    minus_degrange(g) = (-).(degrange(g))
     # The multiplier will have degree `0:2fld(maxdegree - _maxdegree(g), 2)`
     mindegree =
-        -deg_range(p -> -_mindegree(p), p, gs, minus_degrange, -maxdegree:0)
+        -deg_range((-) ∘ _mindegree, p, gs, minus_degrange, -maxdegree:0)
     if isnothing(mindegree)
         return
     end
@@ -474,13 +449,17 @@ function putinar_degree_bounds(
         return
     end
     vars_mindeg = map(vars) do v
-        return -deg_range(
-            Base.Fix2(minus_min_degree, v),
-            p,
-            gs,
-            minus_degrange,
-            -maxdegree:0,
-        )
+        return if _is_monomial_basis(eltype(gs))
+            -deg_range(
+                (-) ∘ Base.Fix2(min_degree, v),
+                p,
+                gs,
+                minus_degrange,
+                -maxdegree:0,
+            )
+        else
+            0
+        end
     end
     if any(isnothing, vars_mindeg)
         return
@@ -501,29 +480,24 @@ function putinar_degree_bounds(
     )
 end
 
-function multiplier_basis(g::SA.AlgebraElement, bounds::DegreeBounds)
+function multiplier_basis(g::SA.AlgebraElement{<:MB.Algebra{BT,B}}, bounds::DegreeBounds) where {BT,B}
     shifted = minus_shift(bounds, g)
     if isnothing(shifted)
         halved = nothing
     else
         halved = _half(shifted)
     end
-    basis = MB.FullBasis{MB.Monomial,MP.monomial_type(typeof(g))}()
+    basis = MB.FullBasis{B,MP.monomial_type(typeof(g))}()
     if isnothing(halved)
-        # TODO add `MB.empty_basis` to API
-        return MB.maxdegree_basis(
-            basis,
-            MP.variables(bounds.variablewise_mindegree),
-            -1,
-        )
+        return MB.empty_basis(MB.explicit_basis_type(typeof(basis)))
     else
         return maxdegree_gram_basis(basis, halved)
     end
 end
 
 function half_newton_polytope(
-    p::MP.AbstractPolynomialLike,
-    gs::AbstractVector{<:MP.AbstractPolynomialLike},
+    p::SA.AlgebraElement,
+    gs::AbstractVector{<:SA.AlgebraElement},
     vars,
     maxdegree,
     ::NewtonDegreeBounds,
@@ -554,10 +528,10 @@ function half_newton_polytope(
     gs::AbstractVector{<:SA.AlgebraElement},
     vars,
     maxdegree,
-    ::NewtonFilter{<:NewtonDegreeBounds},
+    filter::NewtonFilter{<:NewtonDegreeBounds},
 )
     bounds = putinar_degree_bounds(p, gs, vars, maxdegree)
-    bases = [multiplier_basis(g, bounds) for g in gs]
+    bases = half_newton_polytope(p, gs, vars, maxdegree, filter.outer_approximation)
     push!(
         bases,
         maxdegree_gram_basis(
@@ -604,6 +578,7 @@ function Base.:+(c::SignCount, a::SignChange{Missing})
     @assert c.unknown >= -a.Δ
     return SignCount(c.unknown + a.Δ, c.positive, c.negative)
 end
+
 function Base.:+(c::SignCount, a::SignChange{<:Number})
     if a.sign > 0
         @assert c.positive >= -a.Δ
@@ -642,6 +617,29 @@ function increase(cache, counter, generator_sign, monos, mult)
     end
 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
+# will be zero hence we can remove this monomial.
+# See [Proposition 3.7, CLR95], [Theorem 2, L09] or [Section 2.4, BKP16].
+
+# This is generalized here to the case with constraints as detailed in [L23].
+
+# [CLR95] Choi, M. D. and Lam, T. Y. and Reznick, B.
+# *Sum of Squares of Real Polynomials*.
+# Proceedings of Symposia in Pure mathematics (1995)
+#
+# [L09] Löfberg, Johan.
+# *Pre-and post-processing sum-of-squares programs in practice*.
+# IEEE transactions on automatic control 54.5 (2009): 1007-1011.
+#
+# [BKP16] Sabine Burgdorf, Igor Klep, and Janez Povh.
+# *Optimization of polynomials in non-commuting variables*.
+# Berlin: Springer, 2016.
+#
+# [L23] Legat, Benoît
+# *Exploiting the Structure of a Polynomial Optimization Problem*
+# SIAM Conference on Applications of Dynamical Systems, 2023
 function post_filter(poly::SA.AlgebraElement, generators, multipliers_gram_monos)
     counter = MB.algebra_element(
         zero(MP.polynomial_type(MP.monomial_type(typeof(SA.basis(poly))), SignCount)),
@@ -737,8 +735,3 @@ end
 function _sub(basis::SubBasis{B}, I) where {B}
     return SubBasis{B}(basis.monomials[I])
 end
-
-function half_newton_polytope(a::SA.AlgebraElement, args...)
-    @assert SA.basis(a) isa MB.SubBasis{MB.Monomial}
-    half_newton_polytope(SA.coeffs(a, MB.FullBasis{MB.Monomial,MP.monomial_type(typeof(a))}()), args...)
-end