From 2541c71ec45530a17560c28f0543ea161fee0867 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 26 May 2021 08:40:51 -0400 Subject: [PATCH 1/3] Add support for Homotopy Continuation --- Project.toml | 5 +++-- src/extract.jl | 38 +++++++++++++++++++++++++++++--------- test/moment_matrix.jl | 32 +++++++++++++++++++++++--------- 3 files changed, 55 insertions(+), 20 deletions(-) diff --git a/Project.toml b/Project.toml index 21ce4d2..4c43559 100644 --- a/Project.toml +++ b/Project.toml @@ -16,12 +16,13 @@ MultivariateBases = "0.1" MultivariatePolynomials = "0.3" MutableArithmetics = "0.2.10" RowEchelon = "0.2" -SemialgebraicSets = "0.2" +SemialgebraicSets = "0.2.3" julia = "1" [extras] DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" +HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["DynamicPolynomials", "Test"] +test = ["DynamicPolynomials", "HomotopyContinuation", "Test"] diff --git a/src/extract.jl b/src/extract.jl index aed46e8..743e80e 100644 --- a/src/extract.jl +++ b/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...) diff --git a/test/moment_matrix.jl b/test/moment_matrix.jl index cd39ac3..7d389ae 100644 --- a/test/moment_matrix.jl +++ b/test/moment_matrix.jl @@ -1,3 +1,8 @@ +using Test + +import HomotopyContinuation +const HC = HomotopyContinuation + @testset "MomentMatrix" begin Mod.@polyvar x y @test_throws ArgumentError moment_matrix(measure([1], [x]), [y]) @@ -23,15 +28,20 @@ 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 atoms !== nothing - @test atoms ≈ η + for solver in (HC.SemialgebraicSetsHCSolver(; compile = false), DEFAULT_SOLVER) + atoms = extractatoms(ν, 1e-4, lrc, solver) + @test atoms !== nothing + @test atoms ≈ η + end end end @@ -68,9 +78,11 @@ end [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-16)) - atoms = extractatoms(ν, 1e-16, lrc) - @test atoms !== nothing - @test atoms ≈ η + for solver in (HC.SemialgebraicSetsHCSolver(; compile = false), DEFAULT_SOLVER) + atoms = extractatoms(ν, 1e-16, lrc, solver) + @test atoms !== nothing + @test atoms ≈ η + end end end @@ -95,9 +107,11 @@ 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 atoms !== nothing - @test atoms ≈ η + for solver in (HC.SemialgebraicSetsHCSolver(; compile = false), DEFAULT_SOLVER) + atoms = extractatoms(ν, 1e-5, solver) + @test atoms !== nothing + @test atoms ≈ η + end end @testset "[LJP17] Example ?" begin From 755dde0b6bd13e0a5cf6d8e58bae1a00eccedd73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 26 May 2021 20:07:31 -0400 Subject: [PATCH 2/3] Don't use HomotopyContinuation in tests as it does not support Julia v1.0 --- Project.toml | 3 +-- test/moment_matrix.jl | 35 ++++++++++++++++++----------------- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/Project.toml b/Project.toml index 4c43559..0e86f26 100644 --- a/Project.toml +++ b/Project.toml @@ -21,8 +21,7 @@ julia = "1" [extras] DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" -HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["DynamicPolynomials", "HomotopyContinuation", "Test"] +test = ["DynamicPolynomials", "Test"] diff --git a/test/moment_matrix.jl b/test/moment_matrix.jl index 7d389ae..d4bcbcf 100644 --- a/test/moment_matrix.jl +++ b/test/moment_matrix.jl @@ -1,7 +1,12 @@ using Test -import HomotopyContinuation -const HC = HomotopyContinuation +struct DummySolver <: SemialgebraicSets.AbstractAlgebraicSolver end +function SemialgebraicSets.solvealgebraicequations( + ::SemialgebraicSets.AlgebraicSet, + ::DummySolver, +) + error("Dummy solver") +end @testset "MomentMatrix" begin Mod.@polyvar x y @@ -37,11 +42,10 @@ const DEFAULT_SOLVER = SemialgebraicSets.defaultalgebraicsolver([1.0x - 1.0x]) μ = 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)) - for solver in (HC.SemialgebraicSetsHCSolver(; compile = false), DEFAULT_SOLVER) - atoms = extractatoms(ν, 1e-4, lrc, solver) - @test atoms !== nothing - @test atoms ≈ η - end + @test_throws ErrorException("Dummy solver") extractatoms(ν, 1e-4, lrc, DummySolver()) + atoms = extractatoms(ν, 1e-4, lrc, DEFAULT_SOLVER) + @test atoms !== nothing + @test atoms ≈ η end end @@ -78,11 +82,9 @@ end [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-16)) - for solver in (HC.SemialgebraicSetsHCSolver(; compile = false), DEFAULT_SOLVER) - atoms = extractatoms(ν, 1e-16, lrc, solver) - @test atoms !== nothing - @test atoms ≈ η - end + atoms = extractatoms(ν, 1e-16, lrc) + @test atoms !== nothing + @test atoms ≈ η end end @@ -107,11 +109,10 @@ 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)) - for solver in (HC.SemialgebraicSetsHCSolver(; compile = false), DEFAULT_SOLVER) - atoms = extractatoms(ν, 1e-5, solver) - @test atoms !== nothing - @test atoms ≈ η - end + @test_throws ErrorException("Dummy solver") extractatoms(ν, 1e-5, DummySolver()) + atoms = extractatoms(ν, 1e-5, DEFAULT_SOLVER) + @test atoms !== nothing + @test atoms ≈ η end @testset "[LJP17] Example ?" begin From c5d29c2965695bfda15fee4713809fe72930ff2d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 27 May 2021 10:26:52 -0400 Subject: [PATCH 3/3] Fix CI on 32 bits --- test/moment_matrix.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/test/moment_matrix.jl b/test/moment_matrix.jl index d4bcbcf..e636bea 100644 --- a/test/moment_matrix.jl +++ b/test/moment_matrix.jl @@ -83,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