Fix and document XORSpace (#314)

* Fix and document XORSpace

* Add tests

* Fix

docs/src/constraints.md

@@ -295,6 +295,7 @@ Types of sparsity
 SumOfSquares.Certificate.Sparsity.Variable
 SumOfSquares.Certificate.Sparsity.Monomial
 SumOfSquares.Certificate.Sparsity.SignSymmetry
+SumOfSquares.Certificate.Sparsity.XORSpace
 ```
 
 Ideal certificates:

src/Certificate/Sparsity/sign.jl

@@ -2,8 +2,40 @@
     struct SignSymmetry <: Sparsity.Pattern end
 
 Sign symmetry as developed in [Section III.C, L09].
+Let `n` be the number of variables.
+A *sign-symmetry* is a binary vector `r` of `{0, 1}^n`
+such that `dot(r, e)` is even for all exponent `e`.
 
-[L09] Lofberg, Johan.
+Let `o(e)` be the binary vector of `{0, 1}^n` for which the `i`th bit is `i`
+iff the `i`th entry of `e` is odd.
+Let `O` be the set of `o(e)` for exponents of `e`.
+The sign symmetry of `r` is then equivalent to its orthogonality
+with all elements of `O`.
+
+Since we are only interested in the orthogonal subspace, say `R`, of `O`,
+even if `O` is not a linear subspace (i.e., it is not invariant under `xor`),
+we compute its span.
+We start by creating row echelon form of the span of `O` using
+[`XORSpace`](@ref).
+We then compute a basis for `R`.
+The set `R` of sign symmetries will be the span of this basis.
+
+Let `a`, `b` be exponents of monomials of the gram basis. For any `r in R`,
+```
+⟨r, o(a + b)⟩ = ⟨r, xor(o(a), o(b)⟩ = xor(⟨r, o(a)⟩, ⟨r, o(b)⟩)
+```
+For `o(a, b)` to be sign symmetric, this scalar product should be zero for all
+sign symmetry `r`.
+This is equivalent to saying that `⟨r, o(a)⟩` and `⟨r, o(b)⟩` are equal
+for all `r in R`.
+In other words, the projection of `o(a)` and `o(b)` to `R` have the same
+coordinates in the basis.
+If we order the monomials by grouping them by equal coordinates of projection,
+we see that the product that are sign symmetric form a block diagonal
+structure.
+This means that we can group the monomials by these coordinates.
+
+[L09] Löfberg, Johan.
 *Pre-and post-processing sum-of-squares programs in practice*.
 IEEE transactions on automatic control 54, no. 5 (2009): 1007-1011.
 """
@@ -24,7 +56,9 @@ function buckets_sign_symmetry(monos, r, ::Type{T}, ::Type{U}) where {T,U}
     # We store vector of indices instead of vector of monos in the bucket
     # as `monos` is already sorted and the buckets content will be sorted
     # so it will remain sorted and `DynamicPolynomials` can keep the
-    # `MonomialVector` struct in `monos[I]`.
+    # `MonomialVector` struct in `monos[I]`. With MultivariateBases, the
+    # sorted state will be ensured by the basis struct so it will also
+    # serve other MultivariatePolynomials implementations
     buckets = Dict{U,Vector{Int}}()
     for (idx, mono) in enumerate(monos)
         exp = binary_exponent(MP.exponents(mono), T)
@@ -51,11 +85,8 @@ function sign_symmetry(
 end
 function sparsity(
     monos::AbstractVector{<:MP.AbstractMonomial},
-    sp::SignSymmetry,
-    gram_monos::AbstractVector = SumOfSquares.Certificate.monomials_half_newton_polytope(
-        monos,
-        tuple(),
-    ),
+    ::SignSymmetry,
+    gram_monos::AbstractVector{<:MP.AbstractMonomial},
 )
     n = MP.nvariables(monos)
     return sign_symmetry(monos, n, appropriate_type(n), gram_monos)

src/Certificate/Sparsity/xor_space.jl

@@ -1,3 +1,44 @@
+"""
+    mutable struct XORSpace{T<:Integer}
+        dimension::Int
+        basis::Vector{T}
+    end
+
+Basis for a linear subspace of the Hamming space (i.e. set of binary string
+`{0, 1}^n` of length `n`) represented in the bits of an integer of type `T`.
+This is used to implement `Certificate.Sparsity.SignSymmetry`.
+
+Consider the scalar product `⟨a, b⟩` which returns the
+the xor of the bits of `a & b`.
+It is a scalar product since `⟨a, b⟩ = ⟨b, a⟩` and
+`⟨a, xor(b, c)⟩ = xor(⟨a, b⟩, ⟨a, c⟩)`.
+
+We have two options here to compute the orthogonal space.
+
+The first one is to build an orthogonal basis
+with some kind of Gram-Schmidt process and then to obtain the orthogonal
+space by removing the projection from each vector of the basis.
+
+The second option, which is the one we use here is to compute the row echelon
+form and then read out the orthogonal subspace directly from it.
+For instance, if the row echelon form is
+```
+1 a 0 c e
+  b 1 d f
+```
+then the orthogonal basis is
+```
+a 1 b 0 0
+c 0 d 1 0
+e 0 f 0 1
+```
+
+The `basis` vector has `dimension` nonzero elements. Any element added with
+`push!` can be obtained as the `xor` of some of the elements of `basis`.
+Moreover, the `basis` is kept in row echelon form in the sense that the first,
+second, ..., `i - 1`th bits of `basis[i]` are zero and `basis[i]` is zero if its
+`i`th bit is not one.
+"""
 mutable struct XORSpace{T<:Integer}
     dimension::Int
     basis::Vector{T}
@@ -19,18 +60,29 @@ end
 function XORSpace(n)
     return XORSpace{appropriate_type(n)}(n)
 end
+b(x) = bitstring(x)[end-5:end]
+e_i(T, i) = (one(T) << (i - 1))
+has_bit(x, i) = !iszero(x & e_i(typeof(x), i))
 function Base.push!(xs::XORSpace{T}, x::T) where {T}
-    for i in eachindex(xs.basis)
-        if !iszero(x & (one(T) << (i - 1)))
-            if iszero(xs.basis[i])
-                xs.dimension += 1
-                xs.basis[i] = x
-                break
-            else
-                x = xor(x, xs.basis[i])
-            end
+    I = eachindex(xs.basis)
+    for i in I
+        if !iszero(xs.basis) && has_bit(x, i)
+            x = xor(x, xs.basis[i])
+        end
+    end
+    i = findfirst(Base.Fix1(has_bit, x), I)
+    if isnothing(i)
+        return
+    end
+    xs.dimension += 1
+    # Clear `i`th column elements of the basis already added
+    for j in I
+        if has_bit(xs.basis[j], i)
+            xs.basis[j] = xor(xs.basis[j], x)
         end
     end
+    # `x` will be the row used as pivot for the `i`th column
+    xs.basis[i] = x
     return xs
 end
 function orthogonal_complement(xs::XORSpace{T}) where {T}

src/Certificate/newton_polytope.jl

@@ -213,7 +213,11 @@ end
 # 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 [Theorem 2, L09] or [Section 2.4, BKP16].
+# 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*.

test/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

test/sparsity.jl

@@ -1,18 +1,31 @@
 using Test
+import Combinatorics
 using SumOfSquares
 import MultivariateBases as MB
 
+function _xor_complement_test(
+    exps,
+    expected,
+    n = maximum(ndigits(exp, base = 2) for exp in exps; init = 1),
+)
+    for e in Combinatorics.permutations(exps)
+        @test Certificate.Sparsity.xor_complement(e, n) == expected
+    end
+end
+
 function xor_complement_test()
-    @test Certificate.Sparsity.xor_complement([1], 1) == Int[]
-    @test Certificate.Sparsity.xor_complement(Int[], 1) == [1]
-    @test Certificate.Sparsity.xor_complement([1], 2) == [2]
-    @test Certificate.Sparsity.xor_complement([2], 2) == [1]
-    @test Certificate.Sparsity.xor_complement([1, 2], 2) == Int[]
-    @test Certificate.Sparsity.xor_complement([1, 3], 2) == Int[]
-    @test Certificate.Sparsity.xor_complement(Int[], 2) == [1, 2]
-    @test Certificate.Sparsity.xor_complement([7], 3) == [3, 5]
-    @test Certificate.Sparsity.xor_complement([5, 6, 3], 3) == [7]
-    @test Certificate.Sparsity.xor_complement([3], 3) == [3, 4]
+    _xor_complement_test([1], Int[])
+    _xor_complement_test(Int[], [1])
+    _xor_complement_test([1], [2], 2)
+    _xor_complement_test([2], [1])
+    _xor_complement_test([1, 2], Int[])
+    _xor_complement_test([1, 3], Int[])
+    _xor_complement_test(Int[], [1, 2], 2)
+    _xor_complement_test([7], [3, 5])
+    _xor_complement_test([5, 6, 3], [7])
+    _xor_complement_test([3], [3, 4], 3)
+    _xor_complement_test([32, 3, 24, 14, 21, 56], [7, 27])
+    return
 end
 
 function set_monos(bases::Vector{<:MB.MonomialBasis})