Merge pull request #13 from JuliaAlgebra/bl/basis

Add support for arbitrary polynomial basis

.travis.yml

@@ -5,11 +5,13 @@ os:
   - osx
 julia:
   - 1.0
-  - 1.1
+  - 1.3
 notifications:
   email: false
 git:
   depth: 99999999
+before_script:
+  - julia -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/JuliaAlgebra/MultivariateBases.jl", rev="master"))'
 after_success:
   - julia -e 'using Pkg; Pkg.add("Coverage"); cd(Pkg.dir("MultivariateMoments")); using Coverage; Coveralls.submit(process_folder()); Codecov.submit(process_folder())'
 jobs:

Project.toml

@@ -5,14 +5,16 @@ version = "0.2.4"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 RowEchelon = "af85af4c-bcd5-5d23-b03a-a909639aa875"
 SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204"
 
 [compat]
+MultivariateBases = "0.1"
 MultivariatePolynomials = "0.3"
-SemialgebraicSets = "0.2"
 RowEchelon = "0.1"
+SemialgebraicSets = "0.2"
 julia = "1"
 
 [extras]

src/MultivariateMoments.jl

@@ -5,6 +5,9 @@ using LinearAlgebra
 using MultivariatePolynomials
 const MP = MultivariatePolynomials
 
+import MultivariateBases
+const MB = MultivariateBases
+
 export AbstractMeasureLike, AbstractMomentLike
 abstract type AbstractMeasureLike{T} end
 abstract type AbstractMomentLike{T} <: AbstractMeasureLike{T} end

src/extract.jl

@@ -45,13 +45,13 @@ using SemialgebraicSets
 #    r, vals
 #end
 
-function build_system(U::AbstractMatrix, mv::AbstractVector, ztol)
+function build_system(U::AbstractMatrix, basis::MB.MonomialBasis, ztol)
     # System is
     # y = [U 0] * y
     # where y = x[end:-1:1]
     # which is
     # y = U * β
-    m = length(mv)
+    m = length(basis)
     r = size(U, 1)
     pivots = [findfirst(j -> U[i, j] != 0, 1:m) for i in 1:r]
     if any(pivots .== 0)
@@ -60,19 +60,20 @@ function build_system(U::AbstractMatrix, mv::AbstractVector, ztol)
         r = length(pivots)
         U = U[keep, :]
     end
-    β = monovec(mv[m + 1 .- pivots]) # monovec makes sure it stays sorted, TypedPolynomials wouldn't let it sorted
+    monos = basis.monomials
+    β = monovec(monos[m + 1 .- pivots]) # monovec makes sure it stays sorted, TypedPolynomials wouldn't let it sorted
     function equation(i)
         if iszero(r) # sum throws ArgumentError: reducing over an empty collection is not allowed, if r is zero
             z = zero(eltype(β)) * zero(eltype(U))
             s = z + z # For type stability
         else
             s = sum(j -> β[r+1-j] * U[j, i], 1:r)
         end
-        s - mv[m+1-i]
+        s - monos[m+1-i]
     end
-    system = filter(p -> maxdegree(p) > 0, map(equation, 1:length(mv)))
+    system = filter(p -> maxdegree(p) > 0, map(equation, 1:length(monos)))
     # Type instability here :(
-    if mindegree(mv) == maxdegree(mv) # Homogeneous
+    if mindegree(monos) == maxdegree(monos) # Homogeneous
         projectivealgebraicset(system, Buchberger(ztol))
     else
         algebraicset(system, Buchberger(ztol))
@@ -158,7 +159,7 @@ function computesupport!(μ::MomentMatrix, ranktol::Real, dec::LowRankChol=SVDCh
     # so we take √||M|| = √nM
     rref!(W, √(nM) * cM / sqrt(m))
     #r, vals = solve_system(U', μ.x)
-    μ.support = build_system(W, μ.x, √cM) # TODO determine what is better between ranktol and sqrt(ranktol) here
+    μ.support = build_system(W, μ.basis, √cM) # TODO determine what is better between ranktol and sqrt(ranktol) here
 end
 
 """

src/moment_matrix.jl

@@ -3,38 +3,38 @@ export getmat, moment_matrix, symmetric_setindex!
 
 using SemialgebraicSets
 
-abstract type AbstractMomentMatrix{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMeasureLike{T} end
+abstract type AbstractMomentMatrix{T, B<:MB.AbstractPolynomialBasis} <: AbstractMeasureLike{T} end
 
 """
-    mutable struct MomentMatrix{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMeasureLike{T}
+    mutable struct MomentMatrix{T, B<:MultivariateBases.AbstractPolynomialBasis} <: AbstractMeasureLike{T}
         Q::SymMatrix{T}
-        x::MVT
+        basis::B
         support::Union{Nothing, AlgebraicSet}
     end
 
 Measure ``\\nu`` represented by the moments of the monomial matrix ``x x^\\top`` in the symmetric matrix `Q`.
 The set of points that are zeros of all the polynomials ``p`` such that ``\\mathbb{E}_{\\nu}[p] = 0`` is stored in the field `support` when it is computed.
 """
-mutable struct MomentMatrix{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMomentMatrix{T, MT, MVT}
+mutable struct MomentMatrix{T, B<:MB.AbstractPolynomialBasis} <: AbstractMomentMatrix{T, B}
     Q::SymMatrix{T}
-    x::MVT
+    basis::B
     support::Union{Nothing, AlgebraicSet}
 end
-MomentMatrix{T, MT, MVT}(Q::SymMatrix{T}, x::MVT) where {T, MT, MVT} = MomentMatrix{T, MT, MVT}(Q, x, nothing)
+MomentMatrix{T, B}(Q::SymMatrix{T}, basis::MB.AbstractPolynomialBasis) where {T, B} = MomentMatrix{T, B}(Q, basis, nothing)
+function MomentMatrix(Q::SymMatrix{T}, basis::MB.AbstractPolynomialBasis) where T
+    return MomentMatrix{T, typeof(basis)}(Q, basis)
+end
 
-MP.variables(μ::MomentMatrix) = variables(μ.x)
-MP.nvariables(μ::MomentMatrix) = nvariables(μ.x)
+MP.variables(μ::MomentMatrix) = variables(μ.basis)
+MP.nvariables(μ::MomentMatrix) = nvariables(μ.basis)
 
-function MomentMatrix{T}(f::Function, x::AbstractVector{MT}, σ) where {T, MT<:AbstractMonomial}
-    MomentMatrix{T, MT, monovectype(x)}(trimat(T, f, length(x), σ), x)
+function MomentMatrix{T}(f::Function, basis::MB.AbstractPolynomialBasis, σ=1:length(basis)) where T
+    return MomentMatrix(trimat(T, f, length(basis), σ), basis)
 end
-MomentMatrix{T}(f::Function, x::AbstractVector, σ) where T = MomentMatrix{T}(f, monomial.(x), σ)
-function MomentMatrix{T}(f::Function, x::AbstractVector) where T
-    σ, X = sortmonovec(x)
-    MomentMatrix{T}(f, X, σ)
+function MomentMatrix{T}(f::Function, monos::AbstractVector) where T
+    σ, sorted_monos = sortmonovec(monos)
+    return MomentMatrix{T}(f, MB.MonomialBasis(sorted_monos), σ)
 end
-MomentMatrix(f::Function, x) = MomentMatrix{Base.promote_op(f, Int, Int)}(f, x)
-moment_matrix(f::Function, x) = MomentMatrix(f, x)
 
 """
     moment_matrix(μ::Measure, x)
@@ -51,27 +51,28 @@ function moment_matrix(μ::Measure{T}, X) where T
         end
         throw(ArgumentError("μ does not have the moment $(x)"))
     end
-    MomentMatrix{T}(getmom, X)
+    return MomentMatrix{T}(getmom, X)
 end
 
-function MomentMatrix(Q::AbstractMatrix{T}, x, σ) where T
-    MomentMatrix{T}((i,j) -> Q[σ[i], σ[j]], x)
+function MomentMatrix(Q::AbstractMatrix{T}, basis::MB.AbstractPolynomialBasis, σ) where T
+    return MomentMatrix{T}((i, j) -> Q[σ[i], σ[j]], basis)
 end
-function MomentMatrix(Q::AbstractMatrix, x)
-    σ, X = sortmonovec(x)
-    MomentMatrix(Q, X, σ)
+function MomentMatrix(Q::AbstractMatrix, monos::AbstractVector)
+    σ, sorted_monos = MP.sortmonovec(monos)
+    return MomentMatrix(Q, MB.MonomialBasis(sorted_monos), σ)
 end
-moment_matrix(Q::AbstractMatrix, x) = MomentMatrix(Q, x)
+moment_matrix(Q::AbstractMatrix, monos) = MomentMatrix(Q, monos)
 
 function getmat(μ::MomentMatrix)
     Matrix(μ.Q)
 end
 
-function measure(ν::MomentMatrix)
-    n = length(ν.x)
-    measure(ν.Q.Q, [ν.x[i] * ν.x[j] for i in 1:n for j in 1:i])
+function measure(ν::MomentMatrix{T, <:MB.MonomialBasis}) where T
+    n = length(ν.basis)
+    monos = ν.basis.monomials
+    measure(ν.Q.Q, [monos[i] * monos[j] for i in 1:n for j in 1:i])
 end
 
-struct SparseMomentMatrix{T, MT <: MP.AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMomentMatrix{T, MT, MVT}
-    sub_moment_matrices::Vector{MomentMatrix{T, MT, MVT}}
+struct SparseMomentMatrix{T, B <: MB.AbstractPolynomialBasis} <: AbstractMomentMatrix{T, B}
+    sub_moment_matrices::Vector{MomentMatrix{T, B}}
 end

src/symmatrix.jl

@@ -16,6 +16,9 @@ struct SymMatrix{T} <: AbstractMatrix{T}
     n::Int
 end
 
+Base.copy(Q::SymMatrix) = SymMatrix(copy(Q.Q), Q.n)
+Base.map(f::Function, Q::SymMatrix) = SymMatrix(map(f, Q.Q), Q.n)
+
 # i <= j
 #function trimap(i, j, n) # It was used when Q was the lower triangular part
 #    div(n*(n+1), 2) - div((n-i+1)*(n-i+2), 2) + j-i+1

test/moment_matrix.jl

@@ -4,8 +4,15 @@
     @test Q.Q == [4, 2, 3]
     symmetric_setindex!(Q, 5, 1, 2)
     @test Q.Q == [4, 5, 3]
-    symmetric_setindex!(Q, 6, 2, 2)
-    @test Q.Q == [4, 5, 6]
+    P = copy(Q)
+    @test P.n == 2
+    symmetric_setindex!(P, 6, 2, 2)
+    @test Q.Q == [4, 5, 3]
+    @test P.n == 2
+    @test P.Q == [4, 5, 6]
+    R = map(i -> i - 1, P)
+    @test R.n == 2
+    @test R.Q == [3, 4, 5]
 end
 
 @testset "MomentMatrix" begin
@@ -14,7 +21,7 @@ end
     μ = measure(1:3, [x^2, x*y, y^2])
     X = [x, y]
     ν1 = moment_matrix(μ, X)
-    ν2 = moment_matrix((i, j) -> i + j - 1, X)
+    ν2 = MomentMatrix{Int}((i, j) -> i + j - 1, X)
     for ν in (ν1, ν2)
         @test ν.Q[1:4] == [1, 2, 2, 3]
         @test ν.Q[1, 1] == 1
@@ -27,7 +34,7 @@ end
     end
     @test_throws ArgumentError moment_matrix(measure([1], [x]), [y])
     sparse_ν = SparseMomentMatrix([ν1, ν2])
-    @test sparse_ν isa SparseMomentMatrix{Int,eltype(ν1.x),typeof(ν1.x)}
+    @test sparse_ν isa SparseMomentMatrix{Int,typeof(ν1.basis)}
 end
 
 # [HL05] Henrion, D. & Lasserre, J-B.
@@ -65,7 +72,7 @@ end
     # β will be [z, y, x, y*z, x*z, x*y]
     x = monovec([z, y, x, z^2, y*z, y^2, x*z, x*y, x^2])
     # x*β contains x*y*z, x^2*z, x^2*y which are not present so it show fail
-    V = MultivariateMoments.build_system(U', x, sqrt(eps(Float64)))
+    V = MultivariateMoments.build_system(U', MB.MonomialBasis(x), sqrt(eps(Float64)))
     @test iszerodimensional(V)
     testelements(V, [[2.0, 2.0, 2.0], [2.0, 2.0, 0.0], [2.0, 0.0, 2.0], [0.0, 2.0, 2.0], [2.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0], [0.0, 0.0, 0.0]], 1e-11)
 end

test/runtests.jl

@@ -1,5 +1,7 @@
-using MultivariateMoments
 using MultivariatePolynomials
+import MultivariateBases
+const MB = MultivariateBases
+using MultivariateMoments
 using SemialgebraicSets
 
 using Test