Add Cluster Completion for Term Sparsity (#178)

Author: blegat

Project.toml

@@ -5,6 +5,7 @@ version = "0.4.4"
 
 [deps]
 ComplexOptInterface = "25c3070e-1275-41b5-afef-2f982c87090a"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
@@ -19,6 +20,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 ComplexOptInterface = "0.0.2"
+DataStructures = "0.18"
 DynamicPolynomials = "0.3.6"
 JuMP = "0.21"
 MathOptInterface = "0.9.18"

docs/src/constraints.md

@@ -285,6 +285,11 @@ Approximations of the cone of copositive matrices:
 SumOfSquares.CopositiveInner
 ```
 
+Sparsity certificates:
+```@docs
+SumOfSquares.Certificate.MonomialSparsity
+```
+
 Attributes
 ```@docs
 SumOfSquares.PolyJuMP.MomentsAttribute
@@ -320,5 +325,5 @@ SumOfSquares.Certificate.ChordalExtensionGraph.LabelledGraph
 SumOfSquares.Certificate.ChordalExtensionGraph.add_node!
 SumOfSquares.Certificate.ChordalExtensionGraph.add_edge!
 SumOfSquares.Certificate.ChordalExtensionGraph.add_clique!
-SumOfSquares.Certificate.ChordalExtensionGraph.chordal_extension
+SumOfSquares.Certificate.ChordalExtensionGraph.completion
 ```

examples/Symmetry reduction.jl

@@ -1,8 +1,8 @@
-# # Lyapunov Function Search
+# # Symmetry reduction
 
-#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Lyapunov Function Search.ipynb)
-#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Lyapunov Function Search.ipynb)
-# **Adapted from**: SOSTOOLS' SOSDEMO2 (See Section 4.2 of [SOSTOOLS User's Manual](http://sysos.eng.ox.ac.uk/sostools/sostools.pdf))
+#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Symmetry reduction.ipynb)
+#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Symmetry reduction.ipynb)
+# **Adapted from**: https://github.com/kalmarek/SymbolicWedderburn.jl/blob/tw/ex_sos/examples/ex_C4.jl
 
 using Pkg
 pkg"add https://github.com/kalmarek/SymbolicWedderburn.jl#tw/ex_sos"

examples/Term sparsity.jl

@@ -0,0 +1,76 @@
+# # Term sparsity
+
+#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Term sparsity.ipynb)
+#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Term sparsity.ipynb)
+# **Adapted from**: Example 3.5 of [WML20b]
+#
+# [WML20a] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre.
+# *TSSOS: A Moment-SOS hierarchy that exploits term sparsity*.
+# arXiv preprint arXiv:1912.08899 (2020).
+#
+# [WML20b] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre.
+# *Chordal-TSSOS: a moment-SOS hierarchy that exploits term sparsity with chordal extension*.
+# arXiv preprint arXiv:2003.03210 (2020).
+
+using Test #src
+using DynamicPolynomials
+@polyvar x[1:3]
+
+# We would like to find the minimum value of the polynomial
+
+poly = x[1]^2 - 2x[1]*x[2] + 3x[2]^2 - 2x[1]^2*x[2] + 2x[1]^2*x[2]^2 - 2x[2]*x[3] + 6x[3]^2 + 18x[2]^2*x[3] - 54x[2]*x[3]^2 + 142x[2]^2*x[3]^2
+
+# The minimum value of the polynomial can be found to be zero.
+
+import CSDP
+solver = CSDP.Optimizer
+using SumOfSquares
+function sos_min(sparsity)
+    model = Model(solver)
+    @variable(model, t)
+    @objective(model, Max, t)
+    con_ref = @constraint(model, poly - t in SOSCone(), sparsity = sparsity)
+    optimize!(model)
+    return value(t), moment_matrix(con_ref)
+end
+
+bound, ν = sos_min(NoSparsity())
+@test bound ≈ 0 atol=1e-6 #src
+bound
+
+# We find the corresponding minimizer `(0, 0, 0)` by matching the moments
+# of the moment matrix with a dirac measure centered at this minimizer.
+
+extractatoms(ν, 1e-6)
+
+# We can see below that the basis contained 6 monomials hence we needed to use 6x6 PSD matrix variables.
+
+@test ν.basis.monomials == [x[1]*x[2], x[2]*x[3], x[1], x[2], x[3], 1] #src
+ν.basis
+
+# Using the monomial/term sparsity method of [WML20a] based on cluster completion, we find the same bound.
+
+bound, ν = sos_min(MonomialSparsity())
+@test bound ≈ 0 atol=1e-6 #src
+bound
+
+# Which is not suprising as no sparsity reduction could be performed.
+
+@test length(ν.sub_moment_matrices) == 1 #src
+@test ν.sub_moment_matrices[1].basis.monomials == [x[1]*x[2], x[2]*x[3], x[1], x[2], x[3], 1] #src
+[sub.basis for sub in ν.sub_moment_matrices]
+
+# Using the monomial/term sparsity method of [WML20b] based on chordal completion, the lower bound is smaller than 0.
+
+bound, ν = sos_min(MonomialSparsity(ChordalCompletion()))
+@test bound ≈ -0.00355 rtol=1e-3 #src
+bound
+
+# However, this bound was obtained with an SDP with 4 matrices of size 3x3.
+
+@test length(ν.sub_moment_matrices) == 4                                  #src
+@test ν.sub_moment_matrices[1].basis.monomials == [x[2]*x[3], x[2], x[3]] #src
+@test ν.sub_moment_matrices[2].basis.monomials == [x[2]*x[3], x[2], 1]    #src
+@test ν.sub_moment_matrices[3].basis.monomials == [x[1], x[2], 1]         #src
+@test ν.sub_moment_matrices[4].basis.monomials == [x[1]*x[2], x[1], 1]    #src
+[sub.basis for sub in ν.sub_moment_matrices]

src/Certificate/Certificate.jl

@@ -175,6 +175,9 @@ 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")
 
 end

src/Certificate/sparse/ChordalExtensionGraph.jl

@@ -1,7 +1,21 @@
 module ChordalExtensionGraph
 
+using DataStructures
+
+export AbstractCompletion, ChordalCompletion, ClusterCompletion
 export LabelledGraph, add_node!, add_edge!, add_clique!, chordal_extension
 
+abstract type AbstractGreedyAlgorithm end
+
+abstract type AbstractCompletion end
+
+struct ChordalCompletion{A<:AbstractGreedyAlgorithm} <: AbstractCompletion
+    algo::A
+end
+ChordalCompletion() = ChordalCompletion(GreedyFillIn())
+
+struct ClusterCompletion <: AbstractCompletion end
+
 # With a `Vector{Vector{Int}}` with unsorted neighbors, computing fill-in
 # would be inefficient.
 # With a `Vector{Vector{Int}}` with sorted neighbors, adding an edge would be
@@ -213,8 +227,6 @@ function add_clique!(G::LabelledGraph{T}, x::Vector{T}) where T
     add_clique!(G.graph, add_node!.(G, x))
 end
 
-abstract type AbstractGreedyAlgorithm end
-
 struct GreedyFillIn <: AbstractGreedyAlgorithm end
 cache(G::Graph, ::GreedyFillIn) = FillInCache(G)
 heuristic_value(G::FillInCache, node::Int, ::GreedyFillIn) = fill_in(G, node)
@@ -262,10 +274,22 @@ function _greedy_triangulation!(G, algo::AbstractGreedyAlgorithm)
     return candidate_cliques
 end
 
-function chordal_extension(G::Graph, algo::AbstractGreedyAlgorithm)
+"""
+    completion(G::Graph, comp::ChordalCompletion)
+
+Return a chordal extension of `G` and the corresponding maximal cliques.
+
+The algoritm is Algorithm 3 in [BA10] with the GreedyFillIn heuristic of Table I.
+
+[BA10] Bodlaender, Hans L., and Arie MCA Koster.
+Treewidth computations I. Upper bounds.
+Information and Computation 208, no. 3 (2010): 259-275.
+Utrecht University, Utrecht, The Netherlands www.cs.uu.nl
+"""
+function completion(G::Graph, comp::ChordalCompletion)
     H = copy(G)
 
-    candidate_cliques = _greedy_triangulation!(cache(H, algo), algo)
+    candidate_cliques = _greedy_triangulation!(cache(H, comp.algo), comp.algo)
     enable_all_nodes!(H)
 
     sort!.(candidate_cliques)
@@ -291,19 +315,41 @@ function chordal_extension(G::Graph, algo::AbstractGreedyAlgorithm)
 end
 
 """
-    chordal_extension(G::LabelledGraph{T})
-
-Return a chordal extension of G and the corresponding maximal cliques.
+    completion(G::Graph, comp::ChordalCompletion)
 
-The algoritm is Algorithm 3 in [BA10] with the GreedyFillIn heuristic of Table I.
-
-[BA10] Bodlaender, Hans L., and Arie MCA Koster.
-Treewidth computations I. Upper bounds.
-Information and Computation 208, no. 3 (2010): 259-275.
-Utrecht University, Utrecht, The Netherlands www.cs.uu.nl
+Return a cluster completion of `G` and the corresponding maximal cliques.
 """
-function chordal_extension(G::LabelledGraph{T}, algo::AbstractGreedyAlgorithm) where T
-    H, cliques = chordal_extension(G.graph, algo)
+function completion(G::Graph, ::ClusterCompletion)
+    H = copy(G)
+    union_find = IntDisjointSets(num_nodes(G))
+    for from in 1:num_nodes(G)
+        for to in neighbors(G, from)
+            union!(union_find, from, to)
+        end
+    end
+    cliques = [Int[] for i in 1:num_groups(union_find)]
+    ids = zeros(Int, num_nodes(G))
+    k = 0
+    for node in 1:num_nodes(G)
+        root = find_root!(union_find, node)
+        if iszero(ids[root])
+            k += 1
+            ids[root] = k
+        end
+        push!(cliques[ids[root]], node)
+    end
+    @assert k == length(cliques)
+    for clique in cliques
+        add_clique!(H, clique)
+    end
+    return H, cliques
+end
+
+function completion(G::LabelledGraph{T}, comp::AbstractCompletion) where T
+    H, cliques = completion(G.graph, comp)
     return LabelledGraph{T}(G.n2int, G.int2n, H), [[G.int2n[i] for i in clique] for clique in cliques]
 end
+
+chordal_extension(G, algo) = completion(G, ChordalCompletion(algo))
+
 end #module

src/Certificate/sparse/monomial_sparsity.jl

@@ -1,13 +1,41 @@
-struct MonomialSparsity <: Sparsity
+const CEG = SumOfSquares.Certificate.ChordalExtensionGraph
+
+"""
+    struct MonomialSparsity{C<:CEG.AbstractCompletion} <: Sparsity
+        k::Int
+        completion::C
+        use_all_monomials::Bool
+    end
+
+Monomial or term sparsity as developed in [].
+
+# [WML20a] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre.
+# *TSSOS: A Moment-SOS hierarchy that exploits term sparsity*.
+# arXiv preprint arXiv:1912.08899 (2020).
+#
+# [WML20b] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre.
+# *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
+    completion::C
     k::Int
     use_all_monomials::Bool
-    function MonomialSparsity(k::Int=0, use_all_monomials::Bool=false)
-        return new(k, use_all_monomials)
+    function MonomialSparsity(
+        completion::CEG.AbstractCompletion=CEG.ClusterCompletion(),
+        k::Int=0,
+        use_all_monomials::Bool=false,
+    )
+        return new{typeof(completion)}(completion, k, use_all_monomials)
     end
 end
 
+# Note: For some implementation of MultivariatePolynomials such as
+# https://github.com/blegat/CondensedMatterSOS.jl,
+# the product of monomials may be a term so wrap any multiplication of monomials by
+# `MP.monomial`.
+
 const MP = SumOfSquares.MP
-const CEG = SumOfSquares.Certificate.ChordalExtensionGraph
 function monomial_sparsity_graph(monos, P, use_all_monomials::Bool)
     g = CEG.LabelledGraph{eltype(monos)}()
     if use_all_monomials
@@ -20,7 +48,7 @@ function monomial_sparsity_graph(monos, P, use_all_monomials::Bool)
     end
     for a in monos
         for b in monos
-            if (a * b) in P
+            if MP.monomial(a * b) in P
                 if a != b
                     CEG.add_edge!(g, a, b)
                 elseif squares !== nothing
@@ -44,21 +72,21 @@ function _add_monos(add, monos, H, squares)
         end
         for bi in neighbors
             b = H.int2n[bi]
-            add(a * b)
+            add(MP.monomial(a * b))
         end
     end
 end
-function chordal_with_squares(g, squares)
-    H, cliques = CEG.chordal_extension(g, CEG.GreedyFillIn())
+function completion_with_squares(g, squares, completion)
+    H, cliques = CEG.completion(g, completion)
     if squares !== nothing
         cliques = filter(monos -> length(monos) > 1 || (monos[1] in squares), collect(cliques))
     end
     return H, cliques
 end
-function monomial_sparsity_iteration(P, use_all_monomials::Bool, monos)
+function monomial_sparsity_iteration(P, completion, use_all_monomials::Bool, monos)
     P_next = Set{eltype(P)}()
     g, squares = monomial_sparsity_graph(monos, P, use_all_monomials)
-    H, cliques = chordal_with_squares(g, squares)
+    H, cliques = completion_with_squares(g, squares, completion)
     _add_monos(mono -> push!(P_next, mono), monos, H, squares)
     return P_next, cliques
 end
@@ -68,18 +96,18 @@ struct GeneratorP{PT, GT}
 end
 function Base.in(mono, g::GeneratorP)
     return any(g.generator_monos) do g_mono
-        (mono * g_mono) in g.P
+        MP.monomial(mono * g_mono) in g.P
     end
 end
-function monomial_sparsity_iteration(P, use_all_monomials::Bool, monos, multiplier_generator_monos)
-    P_next, cliques = monomial_sparsity_iteration(P, use_all_monomials, monos)
+function monomial_sparsity_iteration(P, completion, use_all_monomials::Bool, monos, multiplier_generator_monos)
+    P_next, cliques = monomial_sparsity_iteration(P, completion, use_all_monomials, monos)
     multiplier_cliques = map(multiplier_generator_monos) do m
         multiplier_monos, generator_monos = m
         g, squares = monomial_sparsity_graph(multiplier_monos, GeneratorP(P, generator_monos), use_all_monomials)
-        H, _cliques = chordal_with_squares(g, squares)
+        H, _cliques = completion_with_squares(g, squares, completion)
         _add_monos(multiplier_monos, H, squares) do mono
             for b in generator_monos
-                push!(P_next, mono * b)
+                push!(P_next, MP.monomial(mono * b))
             end
         end
         return _cliques
@@ -101,7 +129,7 @@ function sparsity(monos::AbstractVector{<:MP.AbstractMonomial}, sp::MonomialSpar
     iter = 0
     while iter < sp.k || iszero(sp.k)
         P_prev = P
-        P, cliques = monomial_sparsity_iteration(P_prev, sp.use_all_monomials, gram_monos, args...)
+        P, cliques = monomial_sparsity_iteration(P_prev, sp.completion, sp.use_all_monomials, gram_monos, args...)
         if iszero(iter)
             # If gram_monos + gram_monos !⊆ monos, then it's possible that P_prev !⊆ P
             P == P_prev && break

src/SumOfSquares.jl

@@ -27,8 +27,8 @@ function matrix_cone_type end
 
 export Certificate
 include("Certificate/Certificate.jl")
-using .Certificate: Sparsity, NoSparsity, VariableSparsity, MonomialSparsity, SignSymmetry
-export Sparsity, NoSparsity, VariableSparsity, MonomialSparsity, SignSymmetry
+using .Certificate: Sparsity, NoSparsity, VariableSparsity, MonomialSparsity, SignSymmetry, ChordalCompletion, ClusterCompletion
+export Sparsity, NoSparsity, VariableSparsity, MonomialSparsity, SignSymmetry, ChordalCompletion, ClusterCompletion
 include("rand.jl")
 
 # MOI extension

test/certificate.jl

@@ -135,7 +135,7 @@ end
         certificate_api(preorder)
         sparsities = Sparsity[VariableSparsity()]
         if certificate isa Certificate.MaxDegree
-            push!(sparsities, MonomialSparsity(1))
+            push!(sparsities, MonomialSparsity(ChordalCompletion(), 1))
         end
         @testset "$(typeof(sparsity))" for sparsity in sparsities
             certificate_api(Certificate.SparsePreorder(sparsity, preorder))
@@ -152,7 +152,7 @@ end
         if certificate isa Certificate.MaxDegree
             _test(Certificate.SparseIdeal(VariableSparsity(), certificate))
         end
-        @testset "$(typeof(sparsity))" for sparsity in [SignSymmetry(), MonomialSparsity(1)]
+        @testset "$(typeof(sparsity))" for sparsity in [SignSymmetry(), MonomialSparsity(ChordalCompletion(), 1)]
             _test(Certificate.SparseIdeal(sparsity, certificate))
             _test(Certificate.SparseIdeal(sparsity, Certificate.Remainder(certificate)))
         end

test/runtests.jl

@@ -33,8 +33,9 @@ include("constraint.jl")
 
 include("Mock/mock_tests.jl")
 
-# Tests needed a solver
-# FIXME these tests should be moved to `Tests` and tested in `Mock`
+# Tests needing a solver
+# FIXME these tests should be converted to Literate and moved to `examples` or
+# converted to be used with `MockOptimizer` and moved to `test/Tests`
 include("solvers.jl")
 include("sospoly.jl")
 include("domain.jl")