From 27b57fce8cbb91abaff550e9bf7839398df7957b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 11 Jul 2023 10:12:42 +0200 Subject: [PATCH] Add MacaulyNullspaceSolver (#57) * Add MacaulyNullspaceSolver * Fix format * Add doc --- docs/src/atoms.md | 8 ++++++- src/echelon.jl | 30 +++++++++++++++++++++++-- src/extract.jl | 41 ++++------------------------------ src/null.jl | 56 +++++++++++++++++++++++++++++++++++++++++++++++ src/shift.jl | 8 +++++-- test/extract.jl | 5 +++-- 6 files changed, 104 insertions(+), 44 deletions(-) diff --git a/docs/src/atoms.md b/docs/src/atoms.md index bb58ec9..2b8d5dd 100644 --- a/docs/src/atoms.md +++ b/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 diff --git a/src/echelon.jl b/src/echelon.jl index ee0a551..dfebb51 100644 --- a/src/echelon.jl +++ b/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 diff --git a/src/extract.jl b/src/extract.jl index e7825ab..5752087 100644 --- a/src/extract.jl +++ b/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 diff --git a/src/null.jl b/src/null.jl index 99b96ee..d7125f4 100644 --- a/src/null.jl +++ b/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 diff --git a/src/shift.jl b/src/shift.jl index b511d7e..0d89cdb 100644 --- a/src/shift.jl +++ b/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 diff --git a/test/extract.jl b/test/extract.jl index 24e2b5b..b060d12 100644 --- a/test/extract.jl +++ b/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)