Sparse submodule (#228)

* sparse -> Sparse

* Sparse module

* Sparse -> Sparsity

* Renamings

* Sparse -> Sparsity

* Renamings

* Update tests

* Add deprecation

* Fixes

* Fix for MOI v0.10.6

* Doc fixes

* Doc fixes

* Fix examples

docs/src/constraints.md

@@ -292,9 +292,9 @@ SumOfSquares.CopositiveInner
 
 Types of sparsity
 ```@docs
-SumOfSquares.Certificate.VariableSparsity
-SumOfSquares.Certificate.MonomialSparsity
-SumOfSquares.Certificate.SignSymmetry
+SumOfSquares.Certificate.Sparsity.Variable
+SumOfSquares.Certificate.Sparsity.Monomial
+SumOfSquares.Certificate.Sparsity.SignSymmetry
 ```
 
 Ideal certificates:
@@ -303,13 +303,14 @@ SumOfSquares.Certificate.MaxDegree
 SumOfSquares.Certificate.FixedBasis
 SumOfSquares.Certificate.Newton
 SumOfSquares.Certificate.Remainder
-SumOfSquares.Certificate.SparseIdeal
+SumOfSquares.Certificate.Sparsity.Ideal
+SumOfSquares.Certificate.Symmetry.Ideal
 ```
 
 Preorder certificates:
 ```@docs
 SumOfSquares.Certificate.Putinar
-SumOfSquares.Certificate.SparsePreorder
+SumOfSquares.Certificate.Sparsity.Preorder
 ```
 
 Attributes
@@ -341,12 +342,12 @@ SumOfSquares.PolyJuMP.bridges
 
 Chordal extension:
 ```@docs
-SumOfSquares.Certificate.ChordalExtensionGraph.neighbors
-SumOfSquares.Certificate.ChordalExtensionGraph.fill_in
-SumOfSquares.Certificate.ChordalExtensionGraph.is_clique
-SumOfSquares.Certificate.ChordalExtensionGraph.LabelledGraph
-SumOfSquares.Certificate.ChordalExtensionGraph.add_node!
-SumOfSquares.Certificate.ChordalExtensionGraph.add_edge!
-SumOfSquares.Certificate.ChordalExtensionGraph.add_clique!
-SumOfSquares.Certificate.ChordalExtensionGraph.completion
+SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.neighbors
+SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.fill_in
+SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.is_clique
+SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.LabelledGraph
+SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.add_node!
+SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.add_edge!
+SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.add_clique!
+SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.completion
 ```

docs/src/tutorials/Sparsity/term_sparsity.jl

@@ -34,7 +34,7 @@ function sos_min(sparsity)
     return value(t), moment_matrix(con_ref)
 end
 
-bound, ν = sos_min(NoSparsity())
+bound, ν = sos_min(Sparsity.NoPattern())
 @test bound ≈ 0 atol=1e-6 #src
 bound
 
@@ -50,7 +50,7 @@ extractatoms(ν, 1e-6)
 
 # Using the monomial/term sparsity method of [WML20a] based on cluster completion, we find the same bound.
 
-bound, ν = sos_min(MonomialSparsity())
+bound, ν = sos_min(Sparsity.Monomial())
 @test bound ≈ 0 atol=1e-6 #src
 bound
 
@@ -62,7 +62,7 @@ bound
 
 # Using the monomial/term sparsity method of [WML20b] based on chordal completion, the lower bound is smaller than 0.
 
-bound, ν = sos_min(MonomialSparsity(ChordalCompletion()))
+bound, ν = sos_min(Sparsity.Monomial(ChordalCompletion()))
 @test bound ≈ -0.00355 rtol=1e-3 #src
 bound
 

examples/chordal_sparsity.jl

@@ -5,8 +5,8 @@ using Test
 include("sparse_polynomials.jl")
 
 function test(test_function, optimizer_constructor)
-    sparse_model = sos_lower_bound(test_function, optimizer_constructor, VariableSparsity())
-    dense_model  = sos_lower_bound(test_function, optimizer_constructor, NoSparsity())
+    sparse_model = sos_lower_bound(test_function, optimizer_constructor, Sparsity.Variable())
+    dense_model  = sos_lower_bound(test_function, optimizer_constructor, Sparsity.NoPattern())
     return abs(objective_value(sparse_model) - objective_value(dense_model))
 end
 

examples/chordal_sparsity_with_domain.jl

@@ -2,7 +2,7 @@ using SumOfSquares
 using DynamicPolynomials
 using Test
 
-function sos_lower_bound(factory, sparsity::Sparsity)
+function sos_lower_bound(factory, sparsity::Sparsity.Pattern)
     @polyvar x y z
     K = @set 1 - x^2 - y^2 >= 0 && 1 - y^2 - z^2 >= 0
     p = 1 - x*y - y*z
@@ -16,10 +16,10 @@ end
 
 function sparse_vs_dense(optimizer_constructor)
     # use chordal sparsity
-    sparse_value = sos_lower_bound(optimizer_constructor, VariableSparsity())
+    sparse_value = sos_lower_bound(optimizer_constructor, Sparsity.Variable())
 
     # no chordal sparsity
-    dense_value = sos_lower_bound(optimizer_constructor, NoSparsity())
+    dense_value = sos_lower_bound(optimizer_constructor, Sparsity.NoPattern())
 
     @test abs(sparse_value - dense_value) < 1e-6
 end

examples/sparse_polynomials.jl

@@ -44,7 +44,7 @@ function generalized_rosenbrock(n::Int)
     return p
 end
 
-function sos_lower_bound(p, factory, sparsity::Sparsity)
+function sos_lower_bound(p, factory, sparsity::Sparsity.Pattern)
     model = Model(factory)
     @variable(model, t)
     @objective(model, Max, t)

src/Certificate/Certificate.jl

@@ -244,11 +244,9 @@ SumOfSquares.matrix_cone_type(::Type{Remainder{GCT}}) where {GCT} = SumOfSquares
 zero_basis(certificate::Remainder) = zero_basis(certificate.gram_certificate)
 zero_basis_type(::Type{Remainder{GCT}}) where {GCT} = zero_basis_type(GCT)
 
-include("sparse/ChordalExtensionGraph.jl")
-using .ChordalExtensionGraph: ChordalCompletion, ClusterCompletion
-export ChordalCompletion, ClusterCompletion
-
-include("sparse/sparse_putinar.jl")
+include("Sparsity/Sparsity.jl")
+using .Sparsity: SignSymmetry, ChordalCompletion, ClusterCompletion
+export Sparsity, SignSymmetry, ChordalCompletion, ClusterCompletion
 
 include("Symmetry/Symmetry.jl")
 export Symmetry

src/Certificate/Sparsity/Sparsity.jl

@@ -0,0 +1,27 @@
+module Sparsity
+
+import MultivariatePolynomials
+const MP = MultivariatePolynomials
+import MultivariateBases
+const MB = MultivariateBases
+using SemialgebraicSets
+
+include("ChordalExtensionGraph.jl")
+using .ChordalExtensionGraph: ChordalCompletion, ClusterCompletion
+export ChordalCompletion, ClusterCompletion
+
+import SumOfSquares
+
+export SignSymmetry
+abstract type Pattern end
+struct NoPattern <: Pattern end
+
+include("xor_space.jl")
+include("sign.jl")
+include("variable.jl")
+include("monomial.jl")
+
+include("preorder.jl")
+include("ideal.jl")
+
+end

src/Certificate/Sparsity/ideal.jl

@@ -0,0 +1,53 @@
+"""
+    struct Sparsity.Ideal{S <: Sparsity.Pattern, C <: AbstractIdealCertificate} <: SumOfSquares.Certificate.AbstractIdealCertificate
+        sparsity::S
+        certificate::C
+    end
+
+Same certificate as `certificate` except that the Sum-of-Squares polynomial `σ`
+is modelled as a sum of Sum-of-Squares polynomials with smaller bases
+using the sparsity reduction `sparsity`.
+"""
+struct Ideal{S <: Pattern, C <: SumOfSquares.Certificate.AbstractIdealCertificate} <: SumOfSquares.Certificate.AbstractIdealCertificate
+    sparsity::S
+    certificate::C
+end
+
+function Ideal(sp::Variable, basis, cone, maxdegree::Nothing, newton_polytope)
+    error("`maxdegree` cannot be `nothing` when `sparsity` is `Sparsity.Variable`.")
+end
+function Ideal(sp::Variable, basis, cone, maxdegree::Integer, newton_polytope)
+    return Ideal(sp, SumOfSquares.Certificate.MaxDegree(cone, basis, maxdegree))
+end
+function Ideal(sp::Union{Monomial, SignSymmetry}, basis, cone, maxdegree, newton_polytope)
+    return Ideal(sp, SumOfSquares.Certificate.Newton(cone, basis, newton_polytope))
+end
+
+function sparsity(poly::MP.AbstractPolynomial, ::Variable, certificate::SumOfSquares.Certificate.MaxDegree)
+    H, cliques = chordal_csp_graph(poly, SemialgebraicSets.FullSpace())
+    return map(cliques) do clique
+        return SumOfSquares.Certificate.maxdegree_gram_basis(certificate.basis, clique, certificate.maxdegree)
+    end
+end
+function sparsity(monos, sp::Union{SignSymmetry, Monomial}, gram_basis::MB.MonomialBasis)
+    return MB.MonomialBasis.(sparsity(monos, sp, gram_basis.monomials))
+end
+function sparsity(poly::MP.AbstractPolynomial, sp::Union{SignSymmetry, Monomial}, certificate::SumOfSquares.Certificate.AbstractIdealCertificate)
+    return sparsity(MP.monomials(poly), sp, SumOfSquares.Certificate.get(certificate, SumOfSquares.Certificate.GramBasis(), poly))
+end
+function SumOfSquares.Certificate.get(certificate::Ideal, ::SumOfSquares.Certificate.GramBasis, poly)
+    return sparsity(poly, certificate.sparsity, certificate.certificate)
+end
+function SumOfSquares.Certificate.get(::Type{Ideal{S, C}}, attr::SumOfSquares.Certificate.GramBasisType) where {S, C}
+    return Vector{<:SumOfSquares.Certificate.get(C, attr)}
+end
+function SumOfSquares.Certificate.get(certificate::Ideal, attr::SumOfSquares.Certificate.ReducedPolynomial, poly, domain)
+    return SumOfSquares.Certificate.get(certificate.certificate, attr, poly, domain)
+end
+function SumOfSquares.Certificate.get(certificate::Ideal, attr::SumOfSquares.Certificate.Cone)
+    return SumOfSquares.Certificate.get(certificate.certificate, attr)
+end
+SumOfSquares.matrix_cone_type(::Type{Ideal{S, C}}) where {S, C} = SumOfSquares.matrix_cone_type(C)
+
+SumOfSquares.Certificate.zero_basis(certificate::Ideal) = SumOfSquares.Certificate.zero_basis(certificate.certificate)
+SumOfSquares.Certificate.zero_basis_type(::Type{Ideal{S, C}}) where {S, C} = SumOfSquares.Certificate.zero_basis_type(C)

src/Certificate/sparse/monomial_sparsity.jl → src/Certificate/Sparsity/monomial.jl

@@ -1,7 +1,7 @@
-const CEG = SumOfSquares.Certificate.ChordalExtensionGraph
+const CEG = ChordalExtensionGraph
 
 """
-    struct MonomialSparsity{C<:CEG.AbstractCompletion} <: Sparsity
+    struct Sparsity.Monomial{C<:CEG.AbstractCompletion} <: Sparsity.Pattern
         completion::C
         k::Int
         use_all_monomials::Bool
@@ -23,11 +23,11 @@ arXiv preprint arXiv:1912.08899 (2020).
 *Chordal-TSSOS: a moment-SOS hierarchy that exploits term sparsity with chordal extension*.
 arXiv preprint arXiv:2003.03210 (2020).
 """
-struct MonomialSparsity{C<:CEG.AbstractCompletion} <: Sparsity
+struct Monomial{C<:CEG.AbstractCompletion} <: Pattern
     completion::C
     k::Int
     use_all_monomials::Bool
-    function MonomialSparsity(
+    function Monomial(
         completion::CEG.AbstractCompletion=CEG.ClusterCompletion(),
         k::Int=0,
         use_all_monomials::Bool=false,
@@ -122,7 +122,7 @@ function monomial_sparsity_iteration(P, completion, use_all_monomials::Bool, mon
 end
 _monovec(cliques::AbstractVector{<:MP.AbstractMonomial}) = MP.monovec(cliques)
 _monovec(cliques) = _monovec.(cliques)
-function sparsity(monos::AbstractVector{<:MP.AbstractMonomial}, sp::MonomialSparsity,
+function sparsity(monos::AbstractVector{<:MP.AbstractMonomial}, sp::Monomial,
                   gram_monos::AbstractVector = SumOfSquares.Certificate.monomials_half_newton_polytope(monos, tuple()),
                   args...)
     P = Set(monos)
@@ -149,21 +149,21 @@ function sparsity(monos::AbstractVector{<:MP.AbstractMonomial}, sp::MonomialSpar
 end
 # This also checks that it is indeed a monomial basis
 _monos(basis::MB.MonomialBasis) = basis.monomials
-function _gram_monos(vars, certificate::MaxDegree{CT, MB.MonomialBasis}) where CT
-    return _monos(maxdegree_gram_basis(MB.MonomialBasis, vars, certificate.maxdegree))
+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
-function sparsity(poly::MP.AbstractPolynomial, domain::BasicSemialgebraicSet, sp::MonomialSparsity, certificate::AbstractPreorderCertificate)
+function sparsity(poly::MP.AbstractPolynomial, domain::SemialgebraicSets.BasicSemialgebraicSet, sp::Monomial, certificate::SumOfSquares.Certificate.AbstractPreorderCertificate)
     gram_monos = _gram_monos(
-        reduce((v, q) -> unique!(sort!([v..., variables(q)...], rev=true)),
-                  domain.p, init = variables(poly)),
-        get(certificate, IdealCertificate())
+        reduce((v, q) -> unique!(sort!([v..., MP.variables(q)...], rev=true)),
+                  domain.p, init = MP.variables(poly)),
+        SumOfSquares.Certificate.get(certificate, SumOfSquares.Certificate.IdealCertificate())
     )
-    processed = get(certificate, PreprocessedDomain(), domain, poly)
+    processed = SumOfSquares.Certificate.get(certificate, SumOfSquares.Certificate.PreprocessedDomain(), domain, poly)
     multiplier_generator_monos = [
-        (_monos(get(certificate, MultiplierBasis(), index, processed)),
-         monomials(get(certificate, Generator(), index, processed)))
-        for index in get(certificate, PreorderIndices(), processed)
+        (_monos(SumOfSquares.Certificate.get(certificate, SumOfSquares.Certificate.MultiplierBasis(), index, processed)),
+         MP.monomials(SumOfSquares.Certificate.get(certificate, SumOfSquares.Certificate.Generator(), index, processed)))
+        for index in SumOfSquares.Certificate.get(certificate, SumOfSquares.Certificate.PreorderIndices(), processed)
     ]
-    cliques, multiplier_cliques = sparsity(monomials(poly), sp, gram_monos, 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/Certificate/Sparsity/preorder.jl

@@ -0,0 +1,47 @@
+"""
+    struct Sparsity.Preorder{S <: Sparsity.Pattern, C <: SumOfSquares.Certificate.AbstractPreorderCertificate} <: SumOfSquares.Certificate.AbstractPreorderCertificate
+        sparsity::S
+        certificate::C
+    end
+
+Same certificate as `C` except that the Sum-of-Squares polynomials `σ_i`
+are modelled as a sum of Sum-of-Squares polynomials with smaller basis
+using the sparsity reduction `sparsity`.
+"""
+struct Preorder{S <: Pattern, C <: SumOfSquares.Certificate.AbstractPreorderCertificate} <: SumOfSquares.Certificate.AbstractPreorderCertificate
+    sparsity::S
+    certificate::C
+end
+
+SumOfSquares.Certificate.get(certificate::Preorder, attr::SumOfSquares.Certificate.Cone) = SumOfSquares.Certificate.get(certificate.certificate, attr)
+
+struct Domain{S, P, B}
+    domain::S
+    processed::P
+    bases::Vector{Vector{B}}
+end
+
+function SumOfSquares.Certificate.get(certificate::Preorder, attr::SumOfSquares.Certificate.PreprocessedDomain, domain::SemialgebraicSets.BasicSemialgebraicSet, p)
+    basis, Preorder_bases = sparsity(p, domain, certificate.sparsity, certificate.certificate)
+    return Domain(domain, SumOfSquares.Certificate.get(certificate.certificate, attr, domain, p), Preorder_bases)
+end
+
+function SumOfSquares.Certificate.get(certificate::Preorder, attr::SumOfSquares.Certificate.PreorderIndices, domain::Domain)
+    return SumOfSquares.Certificate.get(certificate.certificate, attr, domain.processed)
+end
+
+function SumOfSquares.Certificate.get(::Preorder, ::SumOfSquares.Certificate.MultiplierBasis, index::SumOfSquares.Certificate.PreorderIndex, domain::Domain)
+    return domain.bases[index.value]
+end
+function SumOfSquares.Certificate.get(::Type{Preorder{S, C}}, attr::SumOfSquares.Certificate.MultiplierBasisType) where {S, C}
+    return Vector{SumOfSquares.Certificate.get(C, attr)}
+end
+
+function SumOfSquares.Certificate.get(certificate::Preorder, attr::SumOfSquares.Certificate.Generator, index::SumOfSquares.Certificate.PreorderIndex, domain::Domain)
+    return SumOfSquares.Certificate.get(certificate.certificate, attr, index, domain.processed)
+end
+
+SumOfSquares.Certificate.get(certificate::Preorder, attr::SumOfSquares.Certificate.IdealCertificate) = Ideal(certificate.sparsity, SumOfSquares.Certificate.get(certificate.certificate, attr))
+SumOfSquares.Certificate.get(::Type{<:Preorder{S, C}}, attr::SumOfSquares.Certificate.IdealCertificate) where {S, C} = Ideal{S, SumOfSquares.Certificate.get(C, attr)}
+
+SumOfSquares.matrix_cone_type(::Type{Preorder{S, C}}) where {S, C} = SumOfSquares.matrix_cone_type(C)

src/Certificate/sparse/sign.jl → src/Certificate/Sparsity/sign.jl

@@ -1,13 +1,13 @@
 """
-    struct SignSymmetry <: Sparsity end
+    struct SignSymmetry <: Sparsity.Pattern end
 
 Sign symmetry as developed in [Section III.C, L09].
 
 [L09] Lofberg, Johan.
 *Pre-and post-processing sum-of-squares programs in practice*.
 IEEE transactions on automatic control 54, no. 5 (2009): 1007-1011.
 """
-struct SignSymmetry <: Sparsity end
+struct SignSymmetry <: Pattern end
 
 function binary_exponent(exponents, ::Type{T}) where T
     cur = zero(T)

src/Certificate/sparse/variable_sparsity.jl → src/Certificate/Sparsity/variable.jl

@@ -1,11 +1,11 @@
 """
-    struct VariableSparsity <: Sparsity end
+    struct Sparsity.Variable <: Sparsity.Pattern end
 
 Variable or correlative sparsity as developed in [WSMM06].
 
 [WSMM06] Waki, Hayato, Sunyoung Kim, Masakazu Kojima, and Masakazu Muramatsu. "Sums of squares and semidefinite program relaxations for polynomial optimization problems with structured sparsity." SIAM Journal on Optimization 17, no. 1 (2006): 218-242.
 """
-struct VariableSparsity <: Sparsity end
+struct Variable <: Pattern end
 
 const CEG = ChordalExtensionGraph
 
@@ -44,13 +44,13 @@ function chordal_csp_graph(poly::MP.APL, domain::AbstractBasicSemialgebraicSet)
 end
 
 function sparsity(poly::MP.AbstractPolynomial, domain::BasicSemialgebraicSet,
-                  sp::VariableSparsity, certificate::Putinar)
+                  sp::Variable, certificate::SumOfSquares.Certificate.Putinar)
     H, cliques = chordal_csp_graph(poly, domain)
     function bases(q)
         return [
-            maxdegree_gram_basis(certificate.basis, clique,
-                                 multiplier_maxdegree(certificate.maxdegree, q))
-            for clique in cliques if variables(q) ⊆ clique
+            SumOfSquares.Certificate.maxdegree_gram_basis(certificate.basis, clique,
+                                 SumOfSquares.Certificate.multiplier_maxdegree(certificate.maxdegree, q))
+            for clique in cliques if MP.variables(q) ⊆ clique
         ]
     end
     return bases(poly), map(bases, domain.p)

src/Certificate/Symmetry/wedderburn.jl

@@ -56,6 +56,17 @@ function MP.polynomialtype(
     return MP.polynomialtype(PT, V)
 end
 
+"""
+    struct Symmetry.Ideal{C,GT,AT<:SymbolicWedderburn.Action} <: SumOfSquares.Certificate.AbstractIdealCertificate
+        pattern::Symmetry.Pattern{GT,AT}
+        certificate::C
+    end
+
+Same certificate as `certificate` except that the Sum-of-Squares polynomial `σ`
+is modelled as a sum of Sum-of-Squares polynomials with smaller bases
+using the Symbolic Wedderburn decomposition of the symmetry pattern specified
+by `pattern`.
+"""
 struct Ideal{C,GT,AT<:SymbolicWedderburn.Action} <: SumOfSquares.Certificate.AbstractIdealCertificate
     pattern::Pattern{GT,AT}
     certificate::C