Add MacaulyNullspaceSolver (#57)

* Add MacaulyNullspaceSolver

* Fix format

* Add doc

Author: blegat

docs/src/atoms.md

@@ -46,6 +46,13 @@ FlatExtension
 The second approach is to first obtain the image space of the moment matrix,
 represented as a [`MacaulayNullspace`](@ref)
 and to then compute the multiplication matrices from this image space.
+This is implemented by the [`ImageSpaceSolver`](@ref).
+
+```@docs
+MacaulayNullspace
+ImageSpaceSolver
+```
+
 This image space is obtained from a low rank LDLT decomposition of the moment matrix.
 This decomposition can either be obtained by a cholesky or SVD decomposition from which we remove the rows corresponding to the negligeable eigen/singular values.
 
@@ -55,7 +62,6 @@ ShiftCholeskyLDLT
 SVDLDLT
 low_rank_ldlt
 LowRankLDLT
-MacaulayNullspace
 ```
 
 The choice of cutoff between the significant and neglibeable eigen/singular values is

src/echelon.jl

@@ -78,16 +78,42 @@ Foundations of Computational Mathematics 8 (2008): 607-647.
 *Numerical polynomial algebra.*
 Society for Industrial and Applied Mathematics, 2004.
 """
-struct Echelon end
+struct Echelon{S<:Union{Nothing,SemialgebraicSets.AbstractAlgebraicSolver}} <:
+       MacaulayNullspaceSolver
+    solver::S
+end
+Echelon() = Echelon(nothing)
+
+# TODO remove in next breaking release
+function compute_support!(
+    ν::MomentMatrix,
+    rank_check::RankCheck,
+    ldlt::LowRankLDLTAlgorithm,
+    ::Echelon,
+    arg::SemialgebraicSets.AbstractAlgebraicSolver,
+)
+    return compute_support!(ν, rank_check, ldlt, Echelon(arg))
+end
 
 import RowEchelon
 
-function solve(null::MacaulayNullspace, ::Echelon, args...)
+function solve(null::MacaulayNullspace, solver::Echelon)
+    # Ideally, we should determine the pivots with a greedy sieve algorithm [LLR08, Algorithm 1]
+    # so that we have more chance that low order monomials are in β and then more chance
+    # so that the pivots form an order ideal. We just use `rref` which does not implement the sieve
+    # v[i] * β to be in μ.x
+    # but maybe it's sufficient due to the shift structure of the matrix.
+    #
+    # [LLR08] Lasserre, Jean Bernard and Laurent, Monique, and Rostalski, Philipp.
+    # "Semidefinite characterization and computation of zero-dimensional real radical ideals."
+    # Foundations of Computational Mathematics 8 (2008): 607-647.
+
     # If M is multiplied by λ, W is multiplied by √λ
     # so we take √||M|| = √nM
     Z = Matrix(null.matrix')
     RowEchelon.rref!(Z, null.accuracy / sqrt(size(Z, 2)))
     #r, vals = solve_system(U', μ.x)
     # TODO determine what is better between rank_check and sqrt(rank_check) here
+    args = isnothing(solver.solver) ? tuple() : (solver.solver,)
     return build_system(Z, null.basis, √null.accuracy, args...)
 end

src/extract.jl

@@ -2,45 +2,12 @@ using SemialgebraicSets
 import CommonSolve: solve
 
 """
-    compute_support!(ν::MomentMatrix, rank_check, [dec])
+    function compute_support!(ν::MomentMatrix, rank_check, solver) end
 
-Computes the `support` field of `ν`.
-The `rank_check` and `dec` parameters are passed as is to the [`low_rank_ldlt`](@ref) function.
+Computes the `support` field of `ν` wth `solver` using a low-rank decomposition
+with the rank evaluated with `rank_check`.
 """
-function compute_support!(
-    μ::MomentMatrix,
-    rank_check::RankCheck,
-    dec::LowRankLDLTAlgorithm,
-    nullspace_solver = Echelon(),
-    args...,
-)
-    # Ideally, we should determine the pivots with a greedy sieve algorithm [LLR08, Algorithm 1]
-    # so that we have more chance that low order monomials are in β and then more chance
-    # so that the pivots form an order ideal. We just use `rref` which does not implement the sieve
-    # v[i] * β to be in μ.x
-    # but maybe it's sufficient due to the shift structure of the matrix.
-    #
-    # [LLR08] Lasserre, Jean Bernard and Laurent, Monique, and Rostalski, Philipp.
-    # "Semidefinite characterization and computation of zero-dimensional real radical ideals."
-    # Foundations of Computational Mathematics 8 (2008): 607-647.
-    M = value_matrix(μ)
-    chol = low_rank_ldlt(M, dec, rank_check)
-    @assert size(chol.L, 1) == LinearAlgebra.checksquare(M)
-    return μ.support = solve(
-        MacaulayNullspace(chol.L, μ.basis, accuracy(chol)),
-        nullspace_solver,
-        args...,
-    )
-end
-
-function compute_support!(μ::MomentMatrix, rank_check::RankCheck, args...)
-    return compute_support!(
-        μ::MomentMatrix,
-        rank_check::RankCheck,
-        SVDLDLT(),
-        args...,
-    )
-end
+function compute_support! end
 
 # Determines weight
 

src/null.jl

@@ -27,3 +27,59 @@ function Base.getindex(
         null.accuracy,
     )
 end
+
+abstract type MacaulayNullspaceSolver end
+
+"""
+    struct ImageSpaceSolver{A<:LowRankLDLTAlgorithm,S<:MacaulayNullspaceSolver}
+        ldlt::A
+        null::S
+    end
+
+Computes the image space of the moment matrix using `ldlt` and then solve it
+with `null`.
+"""
+struct ImageSpaceSolver{A<:LowRankLDLTAlgorithm,S<:MacaulayNullspaceSolver}
+    ldlt::A
+    null::S
+end
+
+function compute_support!(
+    ν::MomentMatrix,
+    rank_check::RankCheck,
+    solver::ImageSpaceSolver,
+)
+    M = value_matrix(ν)
+    chol = low_rank_ldlt(M, solver.ldlt, rank_check)
+    @assert size(chol.L, 1) == LinearAlgebra.checksquare(M)
+    ν.support =
+        solve(MacaulayNullspace(chol.L, ν.basis, accuracy(chol)), solver.null)
+    return
+end
+
+function compute_support!(
+    ν::MomentMatrix,
+    rank_check::RankCheck,
+    ldlt::LowRankLDLTAlgorithm,
+    null::MacaulayNullspaceSolver,
+)
+    return compute_support!(ν, rank_check, ImageSpaceSolver(ldlt, null))
+end
+
+function compute_support!(
+    ν::MomentMatrix,
+    rank_check::RankCheck,
+    ldlt::LowRankLDLTAlgorithm,
+    args...,
+)
+    return compute_support!(ν, rank_check, ldlt, Echelon(args...))
+end
+
+function compute_support!(μ::MomentMatrix, rank_check::RankCheck, args...)
+    return compute_support!(
+        μ::MomentMatrix,
+        rank_check::RankCheck,
+        SVDLDLT(),
+        args...,
+    )
+end

src/shift.jl

@@ -78,16 +78,20 @@ function gap_zone_standard_monomials(monos, maxdegree)
 end
 
 """
-    struct ShiftNullspace end
+    struct ShiftNullspace{C<:RankCheck} <: MacaulayNullspaceSolver
+        check::C
+    end
 
 From the [`MacaulayNullspace`](@ref), computes multiplication matrices
 by exploiting the shift property of the rows [DBD12].
+The rank check `check` is used to determine the standard monomials among
+the row indices of the null space.
 
 [DBD12] Dreesen, Philippe, Batselier, Kim, and De Moor, Bart.
 *Back to the roots: Polynomial system solving, linear algebra, systems theory.*
 IFAC Proceedings Volumes 45.16 (2012): 1203-1208.
 """
-struct ShiftNullspace{C<:RankCheck}
+struct ShiftNullspace{C<:RankCheck} <: MacaulayNullspaceSolver
     check::C
 end
 # Because the matrix is orthogonal, we know the SVD of the whole matrix is

test/extract.jl

@@ -372,11 +372,12 @@ end
 function test_extract()
     default_solver = SemialgebraicSets.default_algebraic_solver([1.0x - 1.0x])
     for solver in [
-        SVDLDLT(),
-        ShiftCholeskyLDLT(1e-15),
         FlatExtension(),
         FlatExtension(NewtonTypeDiagonalization()),
+        Echelon(),
+        ImageSpaceSolver(ShiftCholeskyLDLT(1e-15), Echelon()),
         ShiftNullspace(),
+        ImageSpaceSolver(ShiftCholeskyLDLT(1e-15), ShiftNullspace()),
     ]
         atoms_1(1e-10, solver)
         atoms_2(1e-10, solver)