Use StarAlgebras.QuadraticForm (#391)

* Use StarAlgebras.QuadraticForm

* Update branch

* Fix format

* Fixes

* Fix format

* Fix

* Back to `StarAlgebras#main`

* Fix

* Fixes

* Fix format

* Remove useless lines

* Missing star

* Fix format

Author: blegat

docs/src/tutorials/Symmetry/cyclic.jl

@@ -90,18 +90,6 @@ model = Model(solver)
 @variable(model, t)
 @objective(model, Max, t)
 pattern = Symmetry.Pattern(G, action)
-basis = MB.explicit_basis(MB.algebra_element(poly - t))
-using SymbolicWedderburn
-summands = SymbolicWedderburn.symmetry_adapted_basis(
-    Float64,
-    pattern.group,
-    pattern.action,
-    basis,
-    semisimple = true,
-)
-
-gram_basis = SumOfSquares.Certificate.Symmetry._gram_basis(pattern, basis, Float64)
-
 con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern)
 optimize!(model)
 @test termination_status(model) == MOI.OPTIMAL #src

src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl

@@ -63,6 +63,8 @@ function MOI.Bridges.Constraint.bridge_constraint(
     λ_constraints = UMCT[]
     preprocessed =
         Certificate.preprocessed_domain(set.certificate, set.domain, p)
+    implicit_basis = MB.implicit_basis(set.basis)
+    cache = zero(MOI.ScalarAffineFunction{T}, MB.algebra(implicit_basis))
     for index in Certificate.preorder_indices(set.certificate, preprocessed)
         λ, λ_variable, λ_constraint, λ_basis = lagrangian_multiplier(
             model,
@@ -79,14 +81,15 @@ function MOI.Bridges.Constraint.bridge_constraint(
         # `Float64` when used with JuMP and the coefficient type is often `Int` if
         # `set.domain.V` is `FullSpace` or `FixedPolynomialSet`.
         g = Certificate.generator(set.certificate, index, preprocessed)
+        MA.operate_to!(cache, +, λ)
         # TODO replace with `MA.sub_mul` when it works.
         p = MA.operate!(
             SA.UnsafeAddMul(*),
             p,
-            λ,
+            cache,
             MB.algebra_element(
                 MB.sparse_coefficients(-one(T) * similar(g, T)),
-                MB.FullBasis{MB.Monomial,MP.monomial_type(g)}(),
+                implicit_basis,
             ),
         )
     end

src/Bridges/Variable/kernel.jl

@@ -12,15 +12,15 @@ function MOI.Bridges.Variable.bridge_constrained_variable(
 ) where {T,M}
     variables = Vector{MOI.VariableIndex}[]
     constraints = MOI.ConstraintIndex{MOI.VectorOfVariables}[]
-    acc = zero(
-        MOI.ScalarAffineFunction{T},
-        MB.algebra(MB.implicit_basis(set.basis)),
-    )
+    algebra = MB.algebra(MB.implicit_basis(set.basis))
+    acc = zero(MOI.ScalarAffineFunction{T}, algebra)
+    cache = zero(MOI.ScalarAffineFunction{T}, algebra)
     for (gram_basis, weight) in zip(set.gram_bases, set.weights)
         gram, vars, con = SOS.add_gram_matrix(model, M, gram_basis, T)
         push!(variables, vars)
         push!(constraints, con)
-        MA.operate!(SA.UnsafeAddMul(*), acc, gram, weight)
+        MA.operate_to!(cache, +, gram)
+        MA.operate!(SA.UnsafeAddMul(*), acc, cache, weight)
     end
     MA.operate!(SA.canonical, SA.coeffs(acc))
     return KernelBridge{T,M}(

src/Certificate/ideal.jl

@@ -20,6 +20,7 @@ function _combine_with_gram(
     weights,
 ) where {B,M}
     p = zero(_NonZero, MB.algebra(MB.FullBasis{B,M}()))
+    cache = zero(_NonZero, MB.algebra(MB.FullBasis{B,M}()))
     for mono in basis
         MA.operate!(
             SA.UnsafeAddMul(*),
@@ -29,12 +30,12 @@ function _combine_with_gram(
         )
     end
     for (gram, weight) in zip(gram_bases, weights)
-        MA.operate!(
-            SA.UnsafeAddMul(*),
-            p,
+        MA.operate_to!(
+            cache,
+            +,
             GramMatrix{_NonZero}((_, _) -> _NonZero(), gram),
-            weight,
         )
+        MA.operate!(SA.UnsafeAddMul(*), p, cache, weight)
     end
     MA.operate!(SA.canonical, SA.coeffs(p))
     return MB.SubBasis{B}(keys(SA.coeffs(p)))

src/Certificate/newton_polytope.jl

@@ -567,6 +567,9 @@ Base.copy(s::SignChange) = s
 Base.iszero(::SignChange) = false
 MA.scaling_convert(::Type, s::SignChange) = s
 Base.:*(s::SignChange, α::Real) = SignChange(s.sign * α, s.Δ)
+Base.:*(α::Real, s::SignChange) = SignChange(α * s.sign, s.Δ)
+#Base.convert(::Type{SignChange{T}}, s::SignChange) where {T} = SignChange{T}(s.sign, s.Δ)
+#Base.:+(a::SignChange, b::SignChange) = convert(SignCount, a) + b
 
 struct SignCount
     unknown::Int
@@ -586,6 +589,9 @@ function _sign(c::SignCount)
         return missing
     end
 end
+function Base.show(io::IO, c::SignCount)
+    return print(io, "[$(c.unknown)|$(c.positive)|$(c.negative)]")
+end
 
 function Base.:*(α, a::SignCount)
     if α > 0
@@ -630,26 +636,16 @@ end
 
 Base.convert(::Type{SignCount}, Δ::SignChange) = SignCount() + Δ
 
-function increase(cache, counter, generator_sign, monos, mult)
-    for a in monos
-        for b in monos
-            MA.operate_to!(
-                cache,
-                *,
-                MB.algebra_element(mult),
-                MB.algebra_element(a),
-                MB.algebra_element(b),
-            )
-            MA.operate!(
-                SA.UnsafeAddMul(*),
-                counter,
-                _term_constant_monomial(
-                    SignChange((a != b) ? missing : generator_sign, 1),
-                    mult,
-                ),
-                cache,
-            )
-        end
+struct SignGram{T,B}
+    sign::T
+    basis::B
+end
+SA.basis(g::SignGram) = g.basis
+function Base.getindex(g::SignGram, i, j)
+    if i == j
+        return SignChange(g.sign, 1)
+    else
+        return SignChange(missing, 1)
     end
 end
 
@@ -715,38 +711,33 @@ function post_filter(
         _DictCoefficients(Dict{MP.monomial_type(typeof(poly)),SignCount}()),
         MB.implicit_basis(SA.basis(poly)),
     )
-    algebra = MB.algebra(MB.implicit_basis(SA.basis(poly)))
-    cache = zero(Float64, algebra)
-    cache3 = zero(SignCount, algebra)
-    cache4 = zero(SignCount, algebra)
+    cache = zero(SignCount, MB.algebra(MB.implicit_basis(SA.basis(poly))))
+    cache2 = zero(SignCount, MB.algebra(MB.implicit_basis(SA.basis(poly))))
     for (mono, v) in SA.nonzero_pairs(SA.coeffs(poly))
         MA.operate!(
-            SA.UnsafeAddMul(*),
+            SA.UnsafeAdd(),
             counter,
             _term(SignChange(_sign(v), 1), SA.basis(poly)[mono]),
         )
     end
     for (mult, gram_monos) in zip(generators, multipliers_gram_monos)
-        for (mono, v) in SA.nonzero_pairs(SA.coeffs(mult))
-            increase(
-                cache,
-                counter,
-                -_sign(v),
-                gram_monos,
-                SA.basis(mult)[mono],
-            )
-        end
+        MA.operate_to!(
+            cache,
+            +,
+            SA.QuadraticForm{SA.star}(SignGram(-1, gram_monos)),
+        )
+        MA.operate!(SA.UnsafeAddMul(*), counter, mult, cache)
     end
     function decrease(sign, a, b, generator)
         MA.operate_to!(
-            cache3,
+            cache,
             *,
             _term(SignChange(sign, -1), a),
             MB.algebra_element(b),
         )
-        MA.operate_to!(cache4, *, cache3, generator)
-        MA.operate!(SA.UnsafeAddMul(*), counter, cache4)
-        for mono in SA.supp(cache4)
+        MA.operate_to!(cache2, *, cache, generator)
+        MA.operate!(SA.UnsafeAdd(), counter, cache2)
+        for mono in SA.supp(cache2)
             count = SA.coeffs(counter)[SA.basis(counter)[mono]]
             count_sign = _sign(count)
             # This means the `counter` has a sign and it didn't have a sign before
@@ -781,16 +772,16 @@ function post_filter(
     for i in eachindex(generators)
         for (j, mono) in enumerate(multipliers_gram_monos[i])
             MA.operate_to!(
-                cache3,
+                cache,
                 *,
                 # Dummy coef to help convert to `SignCount` which is the `eltype` of `cache`
                 _term(SignChange(1, 1), mono),
                 MB.algebra_element(mono),
             )
             # The `eltype` of `cache` is `SignCount`
             # so there is no risk of term cancellation
-            MA.operate_to!(cache4, *, cache3, generators[i])
-            for w in SA.supp(cache4)
+            MA.operate_to!(cache2, *, cache, generators[i])
+            for w in SA.supp(cache2)
                 if ismissing(_sign(SA.coeffs(counter)[SA.basis(counter)[w]]))
                     push!(get!(back, w, Tuple{Int,Int}[]), (i, j))
                 else

src/gram_matrix.jl

@@ -90,6 +90,7 @@ function MP.similar_type(
     return GramMatrix{S,B,US,MS}
 end
 
+SA.basis(g::GramMatrix) = g.basis
 MB.implicit_basis(g::GramMatrix) = MB.implicit_basis(g.basis)
 
 # When taking the promotion of a GramMatrix of JuMP.Variable with a Polynomial JuMP.Variable, it should be a Polynomial of AffExpr
@@ -137,76 +138,13 @@ function GramMatrix(Q::AbstractMatrix, monos::AbstractVector)
     return GramMatrix(Q, MB.SubBasis{MB.Monomial}(sorted_monos), σ)
 end
 
-# `_gram_operate!` is a temporary workaround waiting for a cleaner solution
-# in https://github.com/JuliaAlgebra/StarAlgebras.jl/pull/62
-
-function _gram_operate!(
-    op,
-    p,
-    α,
-    row::MB.Polynomial,
-    col::MB.Polynomial,
-    args::Vararg{Any,N},
-) where {N}
-    return MA.operate!(
-        op,
-        p,
-        MB.term_element(true, row), # TODO `*` shouldn't be defined for group elements so this is a hack
-        MB.term_element(α, col),
-        args...,
-    )
-end
-
-function _gram_operate!(
-    op,
-    p,
-    α,
-    row::SA.AlgebraElement,
-    col::SA.AlgebraElement,
-    args::Vararg{Any,N},
-) where {N}
-    return MA.operate!(
-        op,
-        p,
-        true * row, # TODO `*` shouldn't be defined for group elements so this is a hack
-        α * col,
-        args...,
-    )
-end
-
-function _gram_operate!(
-    op,
-    p,
-    α,
-    row::MB.SemisimpleElement,
-    col::MB.SemisimpleElement,
-    args::Vararg{Any,N},
-) where {N}
-    for (r, c) in zip(row.elements, col.elements)
-        _gram_operate!(op, p, α, r, c, args...)
-    end
-end
-
 function MA.operate!(
     op::SA.UnsafeAddMul{typeof(*)},
     p::SA.AlgebraElement,
     g::GramMatrix,
     args::Vararg{Any,N},
 ) where {N}
-    for row in eachindex(g.basis)
-        row_star = SA.star(g.basis[row])
-        for col in eachindex(g.basis)
-            _gram_operate!(
-                op,
-                p,
-                g.Q[row, col],
-                row_star,
-                g.basis[col],
-                args...,
-            )
-        end
-    end
-    return p
+    return MA.operate!(op, p, SA.QuadraticForm{SA.star}(g), args...)
 end
 
 """
@@ -320,6 +258,51 @@ end
 #    convert(PT, MP.polynomial(p))
 #end
 
+function MA.operate_to!(a::SA.AlgebraElement, ::typeof(+), g::GramMatrix)
+    return MA.operate_to!(a, +, SA.QuadraticForm{SA.star}(g))
+end
+
+function MA.operate!(op::SA.UnsafeAdd, a::SA.AlgebraElement, g::GramMatrix)
+    return MA.operate!(op, a, SA.QuadraticForm{SA.star}(g))
+end
+
+function MA.operate!(
+    op::SA.UnsafeAdd,
+    a::SA.AlgebraElement,
+    g::BlockDiagonalGramMatrix,
+)
+    for block in g.blocks
+        MA.operate!(op, a, block)
+    end
+    return a
+end
+
+function MA.operate_to!(
+    a::SA.AlgebraElement,
+    ::typeof(+),
+    g::BlockDiagonalGramMatrix,
+)
+    MA.operate!(zero, a)
+    MA.operate!(SA.UnsafeAdd(), a, g)
+    MA.operate!(SA.canonical, a)
+    return a
+end
+
+function MB.algebra_element(
+    p::Union{GramMatrix{T,B,U},BlockDiagonalGramMatrix{T,B,U}},
+) where {T,B,U}
+    return MB.algebra_element(p, U)
+end
+
+function MB.algebra_element(
+    g::Union{GramMatrix,BlockDiagonalGramMatrix},
+    ::Type{T},
+) where {T}
+    a = zero(T, MB.algebra(MB.implicit_basis(g)))
+    MA.operate_to!(a, +, g)
+    return a
+end
+
 function MP.polynomial(
     p::Union{GramMatrix{T,B,U},BlockDiagonalGramMatrix{T,B,U}},
 ) where {T,B,U}
@@ -330,10 +313,5 @@ function MP.polynomial(
     g::Union{GramMatrix,BlockDiagonalGramMatrix},
     ::Type{T},
 ) where {T}
-    p = zero(T, MB.algebra(MB.implicit_basis(g)))
-    MA.operate!(SA.UnsafeAddMul(*), p, g)
-    MA.operate!(SA.canonical, SA.coeffs(p))
-    return MP.polynomial(
-        SA.coeffs(p, MB.FullBasis{MB.Monomial,MP.monomial_type(g)}()),
-    )
+    return MP.polynomial(MB.algebra_element(g, T))
 end