From 453722facf416b26da88adc3e80a8ce23522a1a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 20 Jul 2023 14:57:53 +0200 Subject: [PATCH] Staircase gap (#58) * Staircase gap * Rename * Complete * Fixes * Fix format * Fixes * Fix format * Fixes * Fix homogeneous check * Fix * Fixes * Fixes * Fixes * Fix format --- Project.toml | 2 + docs/src/atoms.md | 32 ++++- src/MultivariateMoments.jl | 5 + src/border.jl | 284 +++++++++++++++++++++++++++++++++++++ src/dependence.jl | 266 ++++++++++++++++++++++++++++++++++ src/echelon.jl | 41 +++--- src/measure.jl | 4 + src/moment_matrix.jl | 3 +- src/null.jl | 10 +- src/shift.jl | 171 +++++++++++----------- test/Project.toml | 2 + test/commutativetests.jl | 1 + test/dependence.jl | 105 ++++++++++++++ test/extract.jl | 71 ++++++---- 14 files changed, 854 insertions(+), 143 deletions(-) create mode 100644 src/border.jl create mode 100644 src/dependence.jl create mode 100644 test/dependence.jl diff --git a/Project.toml b/Project.toml index 5d30b6c..28b934c 100644 --- a/Project.toml +++ b/Project.toml @@ -10,8 +10,10 @@ MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378" MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" +RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" RowEchelon = "af85af4c-bcd5-5d23-b03a-a909639aa875" SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] CommonSolve = "0.2" diff --git a/docs/src/atoms.md b/docs/src/atoms.md index 2b8d5dd..e9f7154 100644 --- a/docs/src/atoms.md +++ b/docs/src/atoms.md @@ -64,7 +64,7 @@ low_rank_ldlt LowRankLDLT ``` -The choice of cutoff between the significant and neglibeable eigen/singular values is +The choice of cutoff between the significant and negligeable eigen/singular values is parametrized by the following interface: ```@docs RankCheck @@ -93,12 +93,32 @@ ShiftNullspace Echelon ``` -The [`Echelon`](@ref) uses the RowEchelon package to determine the standard -monomials (which is not numerically stable) while the [`ShiftNullspace`](@ref) -uses the following function internally which is based on SVD so it should have -better numerical behavior. +The [`Echelon`](@ref) uses the RowEchelon package to determine the independent +rows (which is not numerically stable) while the [`ShiftNullspace`](@ref) uses +[`RankCheck`](@ref)s with the singular values so it should have better numerical +behavior. They can either simply distinguish the dependency of rows with +[`AnyDependence`](@ref) or use a sieve with [`StaircaseDependence`](@ref) to +save some the computation of the singular values for some submatrices: + +```@docs +AnyDependence +StaircaseDependence +``` + +The relationship between the dependent and the independent rows are +then stored in a [`BorderBasis`](@ref). +By default, calling `solve` with only a [`BorderBasis`](@ref)` as argument +or providing a `SemialgebraicSets.AbstractMultiplicationMatricesSolver` +as second argument will try to construct the moment matrices from these +and find the solutions from these moment matrices. +Alternatively, give an [`AlgebraicBorderSolver`](@ref) or a +[`AlgebraicFallbackBorderSolver`](@ref) as second argument +to solve the system formed by these relations with these instead. + ```@docs -standard_monomials_and_border +BorderBasis +AlgebraicBorderSolver +AlgebraicFallbackBorderSolver ``` Once the center of the atoms are determined, a linear system is solved to determine diff --git a/src/MultivariateMoments.jl b/src/MultivariateMoments.jl index 8bfd868..87335da 100644 --- a/src/MultivariateMoments.jl +++ b/src/MultivariateMoments.jl @@ -1,6 +1,7 @@ module MultivariateMoments using LinearAlgebra +import SparseArrays import MutableArithmetics as MA @@ -9,6 +10,8 @@ const _APL = MP.AbstractPolynomialLike import MultivariateBases as MB +import SemialgebraicSets as SS + abstract type AbstractMeasureLike{T} end abstract type AbstractMomentLike{T} <: AbstractMeasureLike{T} end @@ -24,6 +27,8 @@ include("moment_matrix.jl") include("atomic.jl") include("rank.jl") +include("dependence.jl") +include("border.jl") include("null.jl") include("extract.jl") include("echelon.jl") diff --git a/src/border.jl b/src/border.jl new file mode 100644 index 0000000..66a4ee8 --- /dev/null +++ b/src/border.jl @@ -0,0 +1,284 @@ +""" + struct BorderBasis{D<:AbstractDependence,T,MT<:AbstractMatrix{T},BT} + dependence::D + matrix::MT + end + +This matrix with rows indexed by `standard` and columns indexed by `border` +a `standard`-border basis of the ideal `border .- matrix' * standard` +[LLR08, Section 2.5]. +For solving this with a multiplication matrix solver, it is necessary for the +basis `border` to be a superset of the set of *corners* of `standard` and it is +sufficient for it to be the *border* of `standard`. + +[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. +""" +struct BorderBasis{D<:AbstractDependence,T,MT<:AbstractMatrix{T}} + dependence::D + matrix::MT +end + +function Base.show(io::IO, b::BorderBasis) + println(io, "BorderBasis with independent rows and dependent columns in:") + println(io, b.dependence) + print(io, "And entries in a ", summary(b.matrix)) + isempty(b.matrix) && return + println(io, ":") + Base.print_matrix(io, b.matrix) + return +end + +BorderBasis{D}(b::BorderBasis{<:D}) where {D} = b + +function BorderBasis{AnyDependence}(b::BorderBasis{<:StaircaseDependence}) + return BorderBasis(convert(AnyDependence, b.dependence), b.matrix) +end + +function BorderBasis{StaircaseDependence}(b::BorderBasis{<:AnyDependence}) + d = StaircaseDependence(_ -> true, b.dependent) + dependent = convert(AnyDependence, d).dependent + rows = _indices(b.dependence.independent, d.standard) + cols = _indices(b.dependence.dependent, dependent) + return BorderBasis(d, b.matrix[rows, cols]) +end + +function solve( + b::BorderBasis{<:StaircaseDependence,T}, + solver::SemialgebraicSets.AbstractMultiplicationMatricesSolver = MultivariateMoments.SemialgebraicSets.ReorderedSchurMultiplicationMatricesSolver{ + T, + }(), +) where {T} + d = b.dependence + dependent = convert(AnyDependence, d).dependent + vars = MP.variables(d.standard) + m = length(d.standard) + # If a monomial `border` is not in `dependent` and it is not a corner + # then it can be divided by another monomial `x^α in dependent`. + # So there exists a monomial `x^β` such that `border = x^(α + β)`. + # In particular, there exists a variable `v` different from `shift` + # that divides `border` and such that `shift * border / v` is in the border. + # If the multiplication matrix for `v` happens to have been computed + # before the multiplication matrix for `shift` then we can obtained + # the column for `shift * border` as the the multiplication matrix for `v` + # multiplied by the column for `shift * border / v`. + # Let `V` be the set of variables such that `shift * border / v` is in the + # border. We need an order such that at least one of the variables of `V` + # is before `shift`. So we need some kind of topological sort on the + # directed forward hypergraph with the variables as nodes and the forward hyperedges + # `shift => V`. Fortunatly, the Kahn's algorithm generalizes itself quite naturally + # to hypergraphs. + order = zeros(Int, length(vars)) + mult = Matrix{T}[zeros(T, m, m) for _ in eachindex(vars)] + completed_border = Dict{eltype(dependent.monomials),Vector{T}}() + function known_border_coefficients(border) + return !isnothing(_index(d.standard, border)) || + !isnothing(_index(dependent, border)) || + haskey(completed_border, border) + end + function border_coefficients(border) + k = _index(d.standard, border) + if !isnothing(k) + return SparseArrays.sparsevec([k], [one(T)], m) + end + k = _index(dependent, border) + if !isnothing(k) + return b.matrix[:, k] + end + if haskey(completed_border, border) + return completed_border[border] + else + return + end + end + function try_add_to_border(border, o) + if known_border_coefficients(border) + return + end + for i in 1:o + v = vars[order[i]] + if MP.divides(v, border) + other_border = MP.div_multiple(border, v) + a = border_coefficients(other_border) + if !isnothing(a) + completed_border[border] = mult[order[i]] * a + return + end + end + end + end + # Variation of Kahn's algorithm + prev_o = -1 + o = 0 + while o > prev_o + prev_o = o + # We iterate over `standard` monomials first` + # and then over variables so that we encounter + # the monomials `shift * std` in increasing + # monomial order so that we know that if `try_add_to_border` + # fails, it will fail again if we run this for loop again with the same `o`. + # That allows to limit the number of iteration of the outer loop by `length(vars)` + for std in d.standard.monomials + for (k, shift) in enumerate(vars) + if k in view(order, 1:o) + continue + end + try_add_to_border(shift * std, o) + end + end + for (k, shift) in enumerate(vars) + if k in view(order, 1:o) + continue + end + if all(d.standard.monomials) do std + return known_border_coefficients(shift * std) + end + o += 1 + order[o] = k + for (col, std) in enumerate(d.standard.monomials) + mult[k][:, col] = border_coefficients(shift * std) + end + end + end + end + @assert o <= length(vars) + if o < length(vars) + # Several things could have gone wrong here: + # 1) We could be missing corners, + # 2) We have all corners but we could not complete the border because + # there is not topological order working + # In any case, we don't have all multiplication matrices so we abort + return + end + sols = SS.solve(SS.MultiplicationMatrices(mult), solver) + return ZeroDimensionalVariety(sols) +end + +""" + Base.@kwdef struct AlgebraicBorderSolver{ + D<:AbstractDependence, + A<:Union{Nothing,SS.AbstractGröbnerBasisAlgorithm}, + S<:Union{Nothing,SS.AbstractAlgebraicSolver}, + } + algorithm::A = nothing + solver::S = nothing + end + +Solve a border basis using `algorithm` and `solver by first converting it to a +`BorderBasis{D}`. +""" +struct AlgebraicBorderSolver{ + D<:AbstractDependence, + A<:Union{Nothing,SS.AbstractGröbnerBasisAlgorithm}, + S<:Union{Nothing,SS.AbstractAlgebraicSolver}, +} + algorithm::A + solver::S +end +function AlgebraicBorderSolver{D}(; + algorithm::Union{Nothing,SS.AbstractGröbnerBasisAlgorithm} = nothing, + solver::Union{Nothing,SS.AbstractAlgebraicSolver} = nothing, +) where {D} + return AlgebraicBorderSolver{D,typeof(algorithm),typeof(solver)}( + algorithm, + solver, + ) +end +#function AlgebraicBorderSolver{D}(solver::SS.AbstractAlgebraicSolver) where {D} +# return AlgebraicBorderSolver{D}(nothing, solver) +#end +#function AlgebraicBorderSolver{D}(algo::SS.AbstractGröbnerBasisAlgorithm) where {D} +# return AlgebraicBorderSolver{D}(algo) +#end + +_some_args(::Nothing) = tuple() +_some_args(arg) = (arg,) + +function solve(b::BorderBasis{E}, solver::AlgebraicBorderSolver{D}) where {D,E} + if !(E <: D) + solve(BorderBasis{D}(b), solver) + end + # Form the system described in [HL05, (8)], note that `b.matrix` + # is the transpose to the matrix in [HL05, (8)]. + # [HL05] Henrion, D. & Lasserre, J-B. + # *Detecting Global Optimality and Extracting Solutions of GloptiPoly* + # 2005 + d = convert(AnyDependence, b.dependence) + system = [ + d.dependent.monomials[col] - + MP.polynomial(b.matrix[:, col], d.independent) for + col in eachindex(d.dependent.monomials) + ] + filter!(!MP.isconstant, system) + V = if MP.mindegree(d) == MP.maxdegree(d) + # Homogeneous + projective_algebraic_set( + system, + _some_args(solver.algorithm)..., + _some_args(solver.solver)..., + ) + else + algebraic_set( + system, + _some_args(solver.algorithm)..., + _some_args(solver.solver)..., + ) + end + return V +end + +""" + struct AlgebraicFallbackBorderSolver{M<:Union{Nothing,SS.AbstractMultiplicationMatricesSolver},S<:Union{Nothing,SS.AbstractAlgebraicSolver}} + multiplication_matrices_solver::M + algebraic_solver::A + end + +Solve with `multiplication_matrices_solver` and if it fails, falls back to +solving the algebraic system formed by the border basis with `algebraic_solver`. +""" +struct AlgebraicFallbackBorderSolver{ + M<:Union{Nothing,SS.AbstractMultiplicationMatricesSolver}, + S<:AlgebraicBorderSolver, +} + multiplication_matrices_solver::M + algebraic_solver::S +end +function AlgebraicFallbackBorderSolver(; + multiplication_matrices_solver::Union{ + Nothing, + SS.AbstractMultiplicationMatricesSolver, + } = nothing, + algorithm::Union{Nothing,SS.AbstractGröbnerBasisAlgorithm} = nothing, + solver::Union{Nothing,SS.AbstractAlgebraicSolver} = nothing, + algebraic_solver::AlgebraicBorderSolver = AlgebraicBorderSolver{ + StaircaseDependence, + }(; + algorithm, + solver, + ), +) + return AlgebraicFallbackBorderSolver( + multiplication_matrices_solver, + algebraic_solver, + ) +end + +function solve(b::BorderBasis, solver::AlgebraicFallbackBorderSolver) + sol = solve(b, _some_args(solver.multiplication_matrices_solver)...) + if isnothing(sol) + sol = solve(b, solver.algebraic_solver) + end + return sol +end + +function solve( + b::BorderBasis{AnyDependence}, + solver::AlgebraicFallbackBorderSolver{ + M, + <:AlgebraicBorderSolver{StaircaseDependence}, + }, +) where {M} + # We will need to convert in both cases anyway + return solve(BorderBasis{StaircaseDependence}(b), solver) +end diff --git a/src/dependence.jl b/src/dependence.jl new file mode 100644 index 0000000..5a08402 --- /dev/null +++ b/src/dependence.jl @@ -0,0 +1,266 @@ +import MultivariatePolynomials as MP +import MultivariateBases as MB +import RecipesBase + +abstract type AbstractDependence end + +""" + struct AnyDependence{B<:MB.AbstractPolynomialBasis} <: AbstractDependence + independent::B + dependent::B + end + +The independent and dependent can be arbitrary disjoint bases. + +!!! tip + If the number of variables is 2 or 3, it can be visualized with Plots.jl. +""" +struct AnyDependence{B<:MB.AbstractPolynomialBasis} <: AbstractDependence + independent::B + dependent::B +end + +function Base.isempty(d::AnyDependence) + return isempty(d.independent.monomials) && isempty(d.dependent.monomials) +end + +function Base.show(io::IO, d::AnyDependence) + println(io, "AnyDependence") + println(io, "with independent basis:") + println(io, d.independent) + println(io, "and dependent basis:") + println(io, d.dependent) + return +end + +""" + struct StaircaseDependence{B<:MB.AbstractPolynomialBasis} <: AbstractDependence + standard::B + corners::B + dependent_border::B + independent_border::B + end + +No independent is a multiple of a dependent and each dependent can be obtained +by multiplying an independent with a variable. + +!!! tip + If the number of variables is 2 or 3, it can be visualized with Plots.jl. +""" +struct StaircaseDependence{B<:MB.AbstractPolynomialBasis} <: AbstractDependence + standard::B + corners::B + dependent_border::B + independent_border::B +end + +function Base.isempty(d::StaircaseDependence) + return isempty(d.standard.monomials) && + isempty(d.corners.monomials) && + isempty(d.dependent_border.monomials) && + isempty(d.independent_border.monomials) +end + +function Base.show(io::IO, d::StaircaseDependence) + println(io, "StaircaseDependence") + println(io, "with standard basis:") + println(io, d.standard) + println(io, "and corners basis:") + println(io, d.corners) + println(io, "and dependent border basis:") + println(io, d.dependent_border) + println(io, "and independent border basis:") + print(io, d.independent_border) + return +end + +function Base.convert(::Type{AnyDependence}, d::StaircaseDependence) + dependent = MB.merge_bases(d.corners, d.dependent_border)[1] + return AnyDependence(d.standard, dependent) +end + +""" + function StaircaseDependence( + is_dependent::Function, + basis::MB.AbstractPolynomialBasis, + ) + +Computes the set of standard monomials using the *greedy sieve* algorithm +presented in [LLR08, Algorithm 1]. A monomial outside of `basis` is assumed +to be independent. Otherwise, if its index in `basis` is `i`, `is_dependent` +returns whether it is dependent. + +[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. +""" +function StaircaseDependence( + is_dependent::Function, + basis::MB.MonomialBasis{M}, +) where {M} + standard = M[] + corners = M[] + vars = MP.variables(basis) + for mono in MP.monomials(vars, 0:MP.maxdegree(basis.monomials)) + # This sieve of [LLR08, Algorithm 1] is a performance improvement but not only. + # It also ensures that the standard monomials have the "staircase structure". + if !any(Base.Fix2(MP.divides, mono), corners) + i = _index(basis, mono) + if isnothing(i) || is_dependent(i) + push!(standard, mono) + else + push!(corners, mono) + end + end + end + dependent_border = M[] + independent_border = M[] + for mono in standard + for shift in vars + border = shift * mono + if isnothing(_monomial_index(standard, border)) && + isnothing(_monomial_index(corners, border)) + i = _index(basis, mono) + if isnothing(i) || is_dependent(i) + push!(independent_border, border) + else + push!(dependent_border, border) + end + end + end + end + return StaircaseDependence( + MB.MonomialBasis(standard), + MB.MonomialBasis(corners), + MB.MonomialBasis(dependent_border), + MB.MonomialBasis(independent_border), + ) +end + +function _exponents(monos, i) + return [MP.exponents(mono)[i] for mono in monos] +end + +function _split_exponents(monos) + N = MP.nvariables(monos) + return ntuple(Base.Fix1(_exponents, monos), Val(N)) +end + +function _split_exponents(basis::MB.AbstractPolynomialBasis) + return _split_exponents(basis.monomials) +end + +function _deg(deg, b, args...) + if isempty(b.monomials) + return nothing + else + return deg(b.monomials, args...) + end +end + +__combine(::Function, ::Nothing) = nothing +__combine(::Function, a) = a +__combine(f::Function, ::Nothing, args...) = __combine(f, args...) +function __combine(f::Function, a, args...) + d = __combine(f, args...) + if isnothing(d) + return a + else + f(a, d) + end +end +function _combine(f, args...) + d = __combine(f, args...) + if isnothing(d) + error("Cannot compute `$(f)degree` as all bases are empty") + end + return d +end + +function _combine_deg(combine, deg, d::AnyDependence, args...) + return _combine( + combine, + _deg(deg, d.independent, args...), + _deg(deg, d.dependent, args...), + ) +end + +function _combine_deg(combine, deg, d::StaircaseDependence, args...) + return _combine( + combine, + _deg(deg, d.standard, args...), + _deg(deg, d.corners, args...), + _deg(deg, d.dependent_border, args...), + _deg(deg, d.independent_border, args...), + ) +end + +function MP.mindegree(d::AbstractDependence, args...) + return _combine_deg(min, MP.mindegree, d, args...) +end + +function MP.maxdegree(d::AbstractDependence, args...) + return _combine_deg(max, MP.maxdegree, d, args...) +end + +function _ticks(d::AbstractDependence, v) + return MP.mindegree(d, v):MP.maxdegree(d, v) +end + +RecipesBase.@recipe function f(m::AnyDependence) + vars = MP.variables(m.independent.monomials) + t = _ticks.(Ref(m), vars) + aspect_ratio --> :equal # defaults to `:auto` + xticks --> t[1] + yticks --> t[2] + if length(t) >= 3 + zticks --> t[3] + end + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :circle + label := "Independent" + _split_exponents(m.independent) + end + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :rect + label := "Dependent" + _split_exponents(m.dependent) + end +end + +RecipesBase.@recipe function f(m::StaircaseDependence) + vars = MP.variables(m.standard.monomials) + t = _ticks.(Ref(m), vars) + aspect_ratio --> :equal # defaults to `:auto` + xticks --> t[1] + yticks --> t[2] + if length(t) >= 3 + zticks --> t[3] + end + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :circle + label := "Standard" + _split_exponents(m.standard) + end + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :rect + label := "Corners" + _split_exponents(m.corners) + end + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :rect + label := "Dependent border" + _split_exponents(m.dependent_border) + end + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :circle + label := "Independent border" + _split_exponents(m.independent_border) + end +end diff --git a/src/echelon.jl b/src/echelon.jl index dfebb51..1c3cc1d 100644 --- a/src/echelon.jl +++ b/src/echelon.jl @@ -1,28 +1,28 @@ -function build_system(U::AbstractMatrix, basis::MB.MonomialBasis, ztol, args...) +function build_system(U::AbstractMatrix, basis::MB.MonomialBasis) # System is - # y = [U 0] * y - # where y = x[end:-1:1] + # `basis = [U' 0] * basis` # which is - # y = U * β + # `y = U' * β` + # where `β` contains the pivots m = length(basis) r = size(U, 1) pivots = [findfirst(j -> U[i, j] != 0, 1:m) for i in 1:r] if any(isnothing, pivots) + # Remove zero rows of `U` keep = map(!isnothing, pivots) pivots = pivots[keep] r = length(pivots) U = U[keep, :] end - monos = basis.monomials - β = monos[pivots] - system = [MA.operate(dot, β, U[:, i]) - monos[i] for i in eachindex(monos)] - filter!(!MP.isconstant, system) - # Type instability here :( - if MP.mindegree(monos) == MP.maxdegree(monos) # Homogeneous - projective_algebraic_set(system, Buchberger(ztol), args...) - else - algebraic_set(system, Buchberger(ztol), args...) - end + set_pivots = Set(pivots) + independent = MP.monomial_vector(basis.monomials[pivots]) + non_pivots = setdiff(eachindex(basis.monomials), set_pivots) + dependent = MP.monomial_vector(basis.monomials[non_pivots]) + d = AnyDependence( + MB.MonomialBasis(independent), + MB.MonomialBasis(dependent), + ) + return BorderBasis(d, U[:, non_pivots]) end """ @@ -84,6 +84,10 @@ struct Echelon{S<:Union{Nothing,SemialgebraicSets.AbstractAlgebraicSolver}} <: end Echelon() = Echelon(nothing) +function border_basis_solver(solver::Echelon) + return AlgebraicBorderSolver{AnyDependence}(solver.solver) +end + # TODO remove in next breaking release function compute_support!( ν::MomentMatrix, @@ -97,7 +101,7 @@ end import RowEchelon -function solve(null::MacaulayNullspace, solver::Echelon) +function border_basis_and_solver(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 @@ -114,6 +118,9 @@ function solve(null::MacaulayNullspace, solver::Echelon) 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...) + solver = AlgebraicBorderSolver{AnyDependence}( + algorithm = SS.Buchberger(√null.accuracy), + solver = solver.solver, + ) + return build_system(Z, null.basis), solver end diff --git a/src/measure.jl b/src/measure.jl index dbd6529..b92b395 100644 --- a/src/measure.jl +++ b/src/measure.jl @@ -90,6 +90,10 @@ function _monomial_index(monos::AbstractVector, mono) end end +function _index(basis::MB.MonomialBasis, mono) + return _monomial_index(basis.monomials, mono) +end + function moment_value(μ, mono) i = _monomial_index(μ.x, mono) if isnothing(i) diff --git a/src/moment_matrix.jl b/src/moment_matrix.jl index aa69c63..5bbca8d 100644 --- a/src/moment_matrix.jl +++ b/src/moment_matrix.jl @@ -65,7 +65,8 @@ function show_basis_indexed_matrix(io::IO, A, pre = "") print(io, pre, "And entries in a ", summary(A.Q)) isempty(A.Q) && return println(io, ":") - return Base.print_matrix(io, A.Q, pre * " ") + Base.print_matrix(io, A.Q, pre * " ") + return end function Base.show(io::IO, M::MomentMatrix) diff --git a/src/null.jl b/src/null.jl index d7125f4..85aeaf7 100644 --- a/src/null.jl +++ b/src/null.jl @@ -2,6 +2,7 @@ struct MacaulayNullspace{T,MT<:AbstractMatrix{T},BT} matrix::MT basis::BT + accuracy::T end This matrix with rows indexed by `basis` can either be interpreted as @@ -10,7 +11,7 @@ the image space of a moment matrix with rows and columns indexed by `basis`). The value of `matrix[i, j]` should be interpreted as the value of the `i`th element of `basis` for the `j`th generator of the null space (resp. image) space. """ -struct MacaulayNullspace{T,MT<:AbstractMatrix{T},BT} +struct MacaulayNullspace{T,MT<:AbstractMatrix{T},BT<:MB.AbstractPolynomialBasis} matrix::MT basis::BT accuracy::T @@ -20,7 +21,7 @@ function Base.getindex( null::MacaulayNullspace{T,MT,<:MB.MonomialBasis}, monos, ) where {T,MT} - I = _monomial_index.(Ref(null.basis.monomials), monos) + I = _index.(Ref(null.basis), monos) return MacaulayNullspace( null.matrix[I, :], MB.MonomialBasis(null.basis.monomials[I]), @@ -30,6 +31,11 @@ end abstract type MacaulayNullspaceSolver end +function solve(null::MacaulayNullspace, solver::MacaulayNullspaceSolver) + border, solver = border_basis_and_solver(null, solver) + return solve(border, _some_args(solver)...) +end + """ struct ImageSpaceSolver{A<:LowRankLDLTAlgorithm,S<:MacaulayNullspaceSolver} ldlt::A diff --git a/src/shift.jl b/src/shift.jl index 0d89cdb..574a905 100644 --- a/src/shift.jl +++ b/src/shift.jl @@ -1,84 +1,88 @@ # Inspired from macaulaylab.net -""" - function standard_monomials_and_border( - null::MacaulayNullspace, - rank_check, - ) +mutable struct RankDependence{T,MT<:AbstractMatrix{T},C} + matrix::MT + check::C + independent_rows::Vector{Int} + old_rank::Int +end -Computes the set of standard monomials using the *greedy sieve* algorithm -presented in [LLR08, Algorithm 1]. +RankDependence(matrix, check) = RankDependence(matrix, check, Int[], 0) -[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. -""" -function standard_monomials_and_border( - null::MacaulayNullspace{T,MT,<:MB.MonomialBasis}, - rank_check, -) where {T,MT} - monos = eltype(null.basis.monomials)[] - border = eltype(monos)[] - rows = Int[] - old_rank = 0 - for (k, mono) in enumerate(null.basis.monomials) - if old_rank == size(null.matrix, 2) - break - end - # This sieve of [LLR08, Algorithm 1] is a performance improvement but not only. - # It also ensures that the standard monomials have the "staircase structure". - if !any(Base.Fix2(MP.divides, mono), border) - new_rank = - LinearAlgebra.rank(null.matrix[vcat(rows, k), :], rank_check) - if new_rank < old_rank - @warn( - "After adding rows, the rank dropped from `$old_rank` to `$new_rank`. Correcting the rank to `$old_rank` and continuing." - ) - new_rank = old_rank - elseif new_rank > old_rank - push!(rows, k) - push!(monos, mono) - else - push!(border, mono) - end - old_rank = new_rank +function is_dependent!(r::RankDependence, row) + if r.old_rank == size(r.matrix, 2) + return false + end + rows = vcat(r.independent_rows, row) + new_rank = LinearAlgebra.rank(r.matrix[rows, :], r.check) + if new_rank < r.old_rank + @warn( + "After adding rows, the rank dropped from `$old_rank` to `$new_rank`. Correcting the rank to `$old_rank` and continuing." + ) + new_rank = r.old_rank + end + independent = new_rank > r.old_rank + r.old_rank = new_rank + if independent + r.independent_rows = rows + end + return independent +end + +function AnyDependence(null::MacaulayNullspace, rank_check::RankCheck) + r = RankDependence(null.matrix, rank_check) + return AnyDependence(Base.Fix1(is_dependent!, r), null.basis) +end + +function StaircaseDependence(null::MacaulayNullspace, rank_check::RankCheck) + r = RankDependence(null.matrix, rank_check) + return StaircaseDependence(Base.Fix1(is_dependent!, r), null.basis) +end + +function _indices_or_ignore(in::MB.MonomialBasis, from::MB.MonomialBasis) + indices = Int[] + for mono in from.monomials + i = _index(in, mono) + if !isnothing(i) + push!(indices, i) end end - return monos, border + return indices +end + +function _indices(in::MB.MonomialBasis, from::MB.MonomialBasis) + return Int[_index(in, mono) for mono in from.monomials] end -function shift_nullspace(null::MacaulayNullspace, monos) - S = null[monos] - Sx = [null[monos.*shift] for shift in MP.variables(monos)] - pS = LinearAlgebra.pinv(S.matrix) - mult = SemialgebraicSets.MultiplicationMatrices([pS * S.matrix for S in Sx]) - return MultivariateMoments.SemialgebraicSets.solve( - mult, - MultivariateMoments.SemialgebraicSets.ReorderedSchurMultiplicationMatricesSolver{ - Float64, - }(), +function BorderBasis(d::AnyDependence, null::MacaulayNullspace) + indep_rows = _indices_or_ignore(null.basis, d.independent) + dep_rows = _indices(null.basis, d.dependent) + @assert length(indep_rows) == size(null.matrix, 2) + return BorderBasis( + d, + (null.matrix[dep_rows, :] / null.matrix[indep_rows, :])', ) end -function gap_zone_standard_monomials(monos, maxdegree) - num = zeros(Int, maxdegree + 1) - for mono in monos - num[MP.maxdegree(mono)+1] += 1 - end - i = findfirst(iszero, num) - if isnothing(i) - return +function BorderBasis(d::StaircaseDependence, null::MacaulayNullspace) + indep_rows = _indices_or_ignore(null.basis, d.standard) + dependent = convert(AnyDependence, d).dependent + dep_rows = _indices(null.basis, dependent) + if length(indep_rows) < size(null.matrix, 2) + error("Column compression not supported yet") end - gap_size = something( - findfirst(!iszero, @view(num[(i+1):end])), - length(num) - i + 1, + return BorderBasis( + d, + (null.matrix[dep_rows, :] / null.matrix[indep_rows, :])', ) - num_affine = sum(view(num, 1:(i-1))) - return num_affine, gap_size +end + +function BorderBasis{D}(null::MacaulayNullspace, check::RankCheck) where {D} + return BorderBasis(D(null, check), null) end """ - struct ShiftNullspace{C<:RankCheck} <: MacaulayNullspaceSolver + struct ShiftNullspace{D,C<:RankCheck} <: MacaulayNullspaceSolver check::C end @@ -91,7 +95,7 @@ the row indices of the null space. *Back to the roots: Polynomial system solving, linear algebra, systems theory.* IFAC Proceedings Volumes 45.16 (2012): 1203-1208. """ -struct ShiftNullspace{C<:RankCheck} <: MacaulayNullspaceSolver +struct ShiftNullspace{D,C<:RankCheck} <: MacaulayNullspaceSolver check::C end # Because the matrix is orthogonal, we know the SVD of the whole matrix is @@ -99,30 +103,15 @@ end # However, since we also know that the first row (which correspond to the # constant monomial) should be a standard monomial, `LeadingRelativeRankTol` # ensures that we will take it. -ShiftNullspace() = ShiftNullspace(LeadingRelativeRankTol(1e-8)) - -function solve(null::MacaulayNullspace, shift::ShiftNullspace) - Z = null.matrix - d = MP.maxdegree(null.basis.monomials) - monos, _ = standard_monomials_and_border(null, shift.check) - gap_zone = gap_zone_standard_monomials(monos, d) - if isnothing(gap_zone) - return - end - num_affine, gap_size = gap_zone - if gap_size < 1 - return - end - if num_affine == length(monos) - affine_monos = monos - affine_null = null - else - affine_monos = monos[1:num_affine] - @warn("Column compression not supported yet") - return - end +function ShiftNullspace{D}(check::RankCheck) where {D} + return ShiftNullspace{D,typeof(check)}(check) +end +ShiftNullspace{D}() where {D} = ShiftNullspace{D}(LeadingRelativeRankTol(1e-8)) +ShiftNullspace(args...) = ShiftNullspace{StaircaseDependence}(args...) - # Solve the system: - sols = shift_nullspace(affine_null, affine_monos) - return ZeroDimensionalVariety(sols) +function border_basis_and_solver( + null::MacaulayNullspace, + shift::ShiftNullspace{D}, +) where {D} + return BorderBasis{D}(null, shift.check), nothing end diff --git a/test/Project.toml b/test/Project.toml index 25f3e75..ba987f1 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,5 +5,7 @@ MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378" MultivariateMoments = "f4abf1af-0426-5881-a0da-e2f168889b5e" MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d" diff --git a/test/commutativetests.jl b/test/commutativetests.jl index 622f16f..a5d9ec0 100644 --- a/test/commutativetests.jl +++ b/test/commutativetests.jl @@ -5,3 +5,4 @@ include("rank.jl") include("extract.jl") include("atomic.jl") include("hermitian_poly.jl") +include("dependence.jl") diff --git a/test/dependence.jl b/test/dependence.jl new file mode 100644 index 0000000..2b16fd2 --- /dev/null +++ b/test/dependence.jl @@ -0,0 +1,105 @@ +module TestDependence + +using Test +import RecipesBase as RB +import MultivariatePolynomials as MP +import MultivariateBases as MB +import MultivariateMoments as MM + +b(x) = MB.MonomialBasis(x) + +function test_degree_error(x, y, z) + M = typeof(x * y * z) + m = b(M[]) + a = MM.AnyDependence(m, m) + @test isempty(a) + s = MM.StaircaseDependence(m, m, m, m) + @test isempty(s) + for (f, name) in [(MP.mindegree, "min"), (MP.maxdegree, "max")] + err = ErrorException( + "Cannot compute `$(name)degree` as all bases are empty", + ) + @test_throws err f(a) + @test_throws err f(a, x) + @test_throws err f(s) + @test_throws err f(s, x) + end +end + +function _test_recipe(dep, ticks, args, names, indep) + d = Dict{Symbol,Any}() + r = RB.apply_recipe(d, dep) + @test length(d) == 1 + length(ticks) + @test d[:aspect_ratio] == :equal + @test d[:xticks] == ticks[1] + @test d[:yticks] == ticks[2] + if length(ticks) >= 3 + @test d[:zticks] == ticks[3] + end + @test length(r) == length(names) + for i in eachindex(names) + @test r[i].args == args[i] + @test r[i].plotattributes[:seriestype] == :scatter + shape = indep[i] ? :circle : :rect + @test r[i].plotattributes[:markershape] == shape + @test r[i].plotattributes[:label] == names[i] + end +end + +function test_recipe(x, y, z) + a = [x^0 * y^0] + A = ([0], [0]) + c = [x, y^2] + C = ([1, 0], [0, 2]) + d = [x * z^1, x * y^2] + D = ([1, 1], [0, 2], [1, 0]) + e = [x * y^2, x * y^3, x * z^4] + E = ([1, 1, 1], [2, 3, 0], [0, 0, 4]) + _test_recipe( + MM.AnyDependence(b(a), b(c)), + [0:1, 0:2], + [A, C], + ["Independent", "Dependent"], + [true, false], + ) + _test_recipe( + MM.AnyDependence(b(d), b(e)), + [1:1, 0:3, 0:4], + [D, E], + ["Independent", "Dependent"], + [true, false], + ) + return _test_recipe( + MM.StaircaseDependence(b(a .* z^0), b(c .* z^0), b(d), b(e)), + [0:1, 0:3, 0:4], + [(A..., [0]), (C..., [0, 0]), D, E], + ["Standard", "Corners", "Dependent border", "Independent border"], + [true, false, false, true], + ) +end + +function runtests(args...) + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)(args...) + end + end + end +end + +end + +using Test + +import DynamicPolynomials +@testset "DynamicPolynomials" begin + DynamicPolynomials.@polyvar x y z + TestDependence.runtests(x, y, z) +end + +import TypedPolynomials +@testset "TypedPolynomials" begin + TypedPolynomials.@polyvar x y z + TestDependence.runtests(x, y, z) +end diff --git a/test/extract.jl b/test/extract.jl index b060d12..d2f2e42 100644 --- a/test/extract.jl +++ b/test/extract.jl @@ -32,7 +32,7 @@ function testelements(X, Y, atol) end function _atoms(atoms, rank_check, solver) - Mod.@polyvar x[1:length(first(atoms))] + Mod.@polyvar x[1:2] η = AtomicMeasure(x, WeightedDiracMeasure.(atoms, ones(length(atoms)))) monos = monomials(x, 0:(length(atoms)+2)) μ = measure(η, monos) @@ -138,16 +138,13 @@ function hl05_3_3_1() 0 2 0 0 0 0 # y^2 = 2y 0 0 0 0 1 0 0 0 0 0 0 1 - 0 0 2 0 0 0 - ] # x^2 = 2x + 0 0 2 0 0 0 # x^2 = 2x + ] # β will be [z, y, x, y*z, x*z, x*y] x = monomial_vector([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', - MB.MonomialBasis(x), - sqrt(eps(Float64)), - ) + # x*β contains x*y*z, x^2*z, x^2*y which are not present so it should fail + b = MultivariateMoments.build_system(U', MB.MonomialBasis(x)) + V = solve(b, AlgebraicBorderSolver{AnyDependence}()) @test is_zero_dimensional(V) return testelements( V, @@ -369,9 +366,12 @@ function large_norm(rank_check) @test nothing === atomic_measure(ν, rank_check) end +_short(x::String) = x[1:min(10, length(x))] +_short(x) = _short(string(x)) + function test_extract() default_solver = SemialgebraicSets.default_algebraic_solver([1.0x - 1.0x]) - for solver in [ + @testset "$(_short(solver))" for solver in [ FlatExtension(), FlatExtension(NewtonTypeDiagonalization()), Echelon(), @@ -379,10 +379,15 @@ function test_extract() ShiftNullspace(), ImageSpaceSolver(ShiftCholeskyLDLT(1e-15), ShiftNullspace()), ] - atoms_1(1e-10, solver) - atoms_2(1e-10, solver) + @testset "Atom 1" begin + atoms_1(1e-10, solver) + end + @testset "Atom 2" begin + atoms_2(1e-10, solver) + end end - for lrc in (SVDLDLT(), ShiftCholeskyLDLT(1e-14)) + @testset "hl05_2_3 $(_short(lrc))" for lrc in + (SVDLDLT(), ShiftCholeskyLDLT(1e-14)) perturb = !(lrc isa ShiftCholeskyLDLT) # the shift `1e-14` is too small compared to the noise of `1e-6`. We want high noise so that the default rtol of `Base.rtoldefault` does not work so that it tests that `rtol` is passed around. hl05_2_3(1e-4, lrc, default_solver, perturb) @test_throws ErrorException("Dummy solver") hl05_2_3( @@ -391,23 +396,36 @@ function test_extract() DummySolver(), ) end - hl05_3_3_1() + @testset "hl05_3_3_1" begin + hl05_3_3_1() + end # Fails on 32-bits in CI if Sys.WORD_SIZE != 32 - for lrc in (SVDLDLT(), ShiftCholeskyLDLT(1e-16)) + @testset "hl05_4 $(_short(lrc))" for lrc in ( + SVDLDLT(), + ShiftCholeskyLDLT(1e-16), + ) hl05_4(1e-16, lrc) end end - # All singular values will be at least 1e-6 > 1e-12 it won't eliminate any row - lpj20_3_8_0(1e-12, ShiftCholeskyLDLT(1e-6), false) - # The following tests that the method does not error if ranktol eliminates everything - # In particular, this tests that the function equation(i) do not call sum when r equal to 0 - # this that throws an ArgumentError as details in src/extract.jl - lpj20_3_8_0(1.0, SVDLDLT(), false) - lpj20_3_8_0(LeadingRelativeRankTol(1e-5), SVDLDLT()) - lpj20_3_8_0(AbsoluteRankTol(1e-5), SVDLDLT()) - lpj20_3_8_0(DifferentialRankTol(1e-2), SVDLDLT()) - lpj20_3_8_0(LargestDifferentialRank(), SVDLDLT()) + @testset "lpj20_3_8_0 $(_short(check)) $(_short(lrc)) $ok" for ( + check, + lrc, + ok, + ) in [ + # All singular values will be at least 1e-6 > 1e-12 it won't eliminate any row + (1e-12, ShiftCholeskyLDLT(1e-6), false) + # The following tests that the method does not error if ranktol eliminates everything + # In particular, this tests that the function equation(i) do not call sum when r equal to 0 + # this that throws an ArgumentError as details in src/extract.jl + (1.0, SVDLDLT(), false) + (LeadingRelativeRankTol(1e-5), SVDLDLT(), true) + (AbsoluteRankTol(1e-5), SVDLDLT(), true) + (DifferentialRankTol(1e-2), SVDLDLT(), true) + (LargestDifferentialRank(), SVDLDLT(), true) + ] + lpj20_3_8_0(check, lrc, ok) + end @test_throws ErrorException("Dummy solver") lpj20_3_8(1e-5, DummySolver()) lpj20_3_8(1e-5, default_solver) for ranktol in [1e-3, 1e-4, 1e-5, 1e-6, 1e-7] @@ -428,7 +446,8 @@ function test_extract() lpj20_3_9(FixedRank(7), 0) jcg14_6_1(6e-3) jcg14_6_1(8e-4, false) - return large_norm(1e-2) + large_norm(1e-2) + return end @testset "Atom extraction" begin