Add support for Homotopy Continuation (#23)

* Add support for Homotopy Continuation

* Don't use HomotopyContinuation in tests as it does not support Julia v1.0

* Fix CI on 32 bits

Project.toml

@@ -16,7 +16,7 @@ MultivariateBases = "0.1"
 MultivariatePolynomials = "0.3"
 MutableArithmetics = "0.2.10"
 RowEchelon = "0.2"
-SemialgebraicSets = "0.2"
+SemialgebraicSets = "0.2.3"
 julia = "1"
 
 [extras]

src/extract.jl

@@ -45,7 +45,7 @@ using SemialgebraicSets
 #    r, vals
 #end
 
-function build_system(U::AbstractMatrix, basis::MB.MonomialBasis, ztol)
+function build_system(U::AbstractMatrix, basis::MB.MonomialBasis, ztol, args...)
     # System is
     # y = [U 0] * y
     # where y = x[end:-1:1]
@@ -74,9 +74,9 @@ function build_system(U::AbstractMatrix, basis::MB.MonomialBasis, ztol)
     system = filter(p -> maxdegree(p) > 0, map(equation, 1:length(monos)))
     # Type instability here :(
     if mindegree(monos) == maxdegree(monos) # Homogeneous
-        projectivealgebraicset(system, Buchberger(ztol))
+        projectivealgebraicset(system, Buchberger(ztol), args...)
     else
-        algebraicset(system, Buchberger(ztol))
+        algebraicset(system, Buchberger(ztol), args...)
     end
 end
 
@@ -146,7 +146,7 @@ end
 Computes the `support` field of `ν`.
 The `ranktol` and `dec` parameters are passed as is to the [`lowrankchol`](@ref) function.
 """
-function computesupport!(μ::MomentMatrix, ranktol::Real, dec::LowRankChol=SVDChol())
+function computesupport!(μ::MomentMatrix, ranktol::Real, dec::LowRankChol, args...)
     # We reverse the ordering so that the first columns corresponds to low order monomials
     # so that we have more chance that low order monomials are in β and then more chance
     # v[i] * β to be in μ.x
@@ -159,14 +159,34 @@ 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, μ.basis, √cM) # TODO determine what is better between ranktol and sqrt(ranktol) here
+    μ.support = build_system(W, μ.basis, √cM, args...) # TODO determine what is better between ranktol and sqrt(ranktol) here
 end
 
-"""
-    extractatoms(ν::MomentMatrix, ranktol, [dec])
+function computesupport!(μ::MomentMatrix, ranktol::Real, args...)
+    return computesupport!(μ::MomentMatrix, ranktol::Real, SVDChol(), args...)
+end
 
-Return an `AtomicMeasure` with the atoms of `ν` if it is atomic or `nothing` if `ν` is not atomic.
-The `ranktol` and `dec` parameters are passed as is to the [`lowrankchol`](@ref) function.
+"""
+    extractatoms(ν::MomentMatrix, ranktol, [dec::LowRankChol], [solver::SemialgebraicSets.AbstractAlgebraicSolver])
+
+Return an `AtomicMeasure` with the atoms of `ν` if it is atomic or `nothing` if
+`ν` is not atomic. The `ranktol` and `dec` parameters are passed as is to the
+[`lowrankchol`](@ref) function. By default, `dec` is an instance of
+[`SVDChol`](@ref). The extraction relies on the solution of a system of
+algebraic equations. using `solver`. For instance, given a
+[`MomentMatrix`](@ref), `μ`, the following extract atoms using a `ranktol` of
+`1e-4` for the low-rank decomposition and homotopy continuation to solve the
+obtained system of algebraic equations.
+```julia
+using HomotopyContinuation
+solver = SemialgebraicSetsHCSolver(; compile = false)
+atoms = extractatoms(ν, 1e-4, solver)
+```
+If no solver is given, the default solver of SemialgebraicSets is used which
+currently computes the Gröbner basis, then the multiplication matrices and
+then the Schur decomposition of a random combination of these matrices.
+For floating point arithmetics, homotopy continuation is recommended as it is
+more numerically stable than Gröbner basis computation.
 """
 function extractatoms(ν::MomentMatrix{T}, ranktol, args...) where T
     computesupport!(ν, ranktol, args...)

test/moment_matrix.jl

@@ -1,3 +1,13 @@
+using Test
+
+struct DummySolver <: SemialgebraicSets.AbstractAlgebraicSolver end
+function SemialgebraicSets.solvealgebraicequations(
+    ::SemialgebraicSets.AlgebraicSet,
+    ::DummySolver,
+)
+    error("Dummy solver")
+end
+
 @testset "MomentMatrix" begin
     Mod.@polyvar x y
     @test_throws ArgumentError moment_matrix(measure([1], [x]), [y])
@@ -23,13 +33,17 @@ end
 # [HL05] Henrion, D. & Lasserre, J-B.
 # Detecting Global Optimality and Extracting Solutions of GloptiPoly 2005
 
+Mod.@polyvar x
+const DEFAULT_SOLVER = SemialgebraicSets.defaultalgebraicsolver([1.0x - 1.0x])
+
 @testset "[HL05] Section 2.3" begin
     Mod.@polyvar x y
     η = AtomicMeasure([x, y], WeightedDiracMeasure.([[1, 2], [2, 2], [2, 3]], [0.4132, 0.3391, 0.2477]))
     μ = measure(η, [x^4, x^3*y, x^2*y^2, x*y^3, y^4, x^3, x^2*y, x*y^2, y^3, x^2, x*y, y^2, x, y, 1])
     ν = moment_matrix(μ, [1, x, y, x^2, x*y, y^2])
     for lrc in (SVDChol(), ShiftChol(1e-14))
-        atoms = extractatoms(ν, 1e-4, lrc)
+        @test_throws ErrorException("Dummy solver") extractatoms(ν, 1e-4, lrc, DummySolver())
+        atoms = extractatoms(ν, 1e-4, lrc, DEFAULT_SOLVER)
         @test atoms !== nothing
         @test atoms ≈ η
     end
@@ -69,8 +83,11 @@ end
     ν = moment_matrix(μ, [1, x, y, x^2, x*y, y^2])
     for lrc in (SVDChol(), ShiftChol(1e-16))
         atoms = extractatoms(ν, 1e-16, lrc)
-        @test atoms !== nothing
-        @test atoms ≈ η
+        # Fails on 32-bits in CI
+        if Sys.WORD_SIZE != 32
+            @test atoms !== nothing
+            @test atoms ≈ η
+        end
     end
 end
 
@@ -95,7 +112,8 @@ end
     Mod.@polyvar x[1:2]
     η = AtomicMeasure(x, [WeightedDiracMeasure([1., 0.], 2.53267)])
     ν = moment_matrix([2.53267 -0.0 -5.36283e-19; -0.0 -5.36283e-19 -0.0; -5.36283e-19 -0.0 7.44133e-6], monomials(x, 2))
-    atoms = extractatoms(ν, 1e-5)
+    @test_throws ErrorException("Dummy solver") extractatoms(ν, 1e-5, DummySolver())
+    atoms = extractatoms(ν, 1e-5, DEFAULT_SOLVER)
     @test atoms !== nothing
     @test atoms ≈ η
 end