Add partial multiplication matrices (#73)

* Add partial multiplication matrices

* Fixes

* Fixes

* Fixes

* fix format

Author: blegat

src/border.jl

@@ -62,20 +62,40 @@ function BorderBasis{StaircaseDependence}(b::BorderBasis{LinearDependence})
     return BorderBasis(d, b.matrix[rows, cols])
 end
 
+struct StaircaseSolver{
+    T,
+    R<:RankCheck,
+    M<:SemialgebraicSets.AbstractMultiplicationMatricesSolver,
+}
+    max_partial_iterations::Int
+    max_iterations::Int
+    rank_check::R
+    solver::M
+end
+function StaircaseSolver{T}(;
+    max_partial_iterations::Int = 0,
+    max_iterations::Int = -1,
+    rank_check::RankCheck = LeadingRelativeRankTol(Base.rtoldefault(T)),
+    solver = SS.ReorderedSchurMultiplicationMatricesSolver{T}(),
+) where {T}
+    return StaircaseSolver{T,typeof(rank_check),typeof(solver)}(
+        max_partial_iterations,
+        max_iterations,
+        rank_check,
+        solver,
+    )
+end
+
 function solve(
     b::BorderBasis{LinearDependence,T},
-    solver::SemialgebraicSets.AbstractMultiplicationMatricesSolver = MultivariateMoments.SemialgebraicSets.ReorderedSchurMultiplicationMatricesSolver{
-        T,
-    }(),
+    solver::StaircaseSolver = StaircaseSolver{T}(),
 ) where {T}
     return solve(BorderBasis{StaircaseDependence}(b), solver)
 end
 
 function solve(
     b::BorderBasis{<:StaircaseDependence,T},
-    solver::SemialgebraicSets.AbstractMultiplicationMatricesSolver = MultivariateMoments.SemialgebraicSets.ReorderedSchurMultiplicationMatricesSolver{
-        T,
-    }(),
+    solver::StaircaseSolver{T} = StaircaseSolver{T}(),
 ) where {T}
     d = b.dependence
     dependent = dependent_basis(d)
@@ -180,34 +200,70 @@ function solve(
         end
     end
     @assert o <= length(vars)
-    if o < length(vars)
+    partial = o < length(vars)
+    if partial
         # 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
+        # We now try to build new relation by comparing partial multiplication matrices
+        # We store them in a vector and reshape in a matrix after as it's easy to append to a vector in-place.
+        # a matrix after
+        if solver.max_partial_iterations == 0
+            return
+        else
+            com_fix = partial_commutation_fix(
+                known_border_coefficients,
+                border_coefficients,
+                T,
+                standard,
+                vars,
+                solver.rank_check,
+            )
+        end
+    else
+        if solver.max_iterations == 0
+            Uperp = nothing
+        else
+            Uperp = commutation_fix(mult, solver.solver.ε)
+        end
+        com_fix = if isnothing(Uperp)
+            nothing
+        else
+            Uperp, standard
+        end
     end
-    Uperp = commutation_fix(mult, solver.ε)
-    if isnothing(Uperp)
+    if isnothing(com_fix)
         # The matrices commute, let's simultaneously diagonalize them
-        sols = SS.solve(SS.MultiplicationMatrices(mult), solver)
+        sols = SS.solve(SS.MultiplicationMatrices(mult), solver.solver)
         return ZeroDimensionalVariety(sols)
     else
+        Uperp, Ubasis = com_fix
         # The matrices don't commute, let's find the updated staircase and start again
-        new_basis, I1, I2 = MB.merge_bases(standard, dependent)
+        new_basis, I1, I2 = MB.merge_bases(Ubasis, dependent)
         new_matrix = Matrix{T}(undef, length(new_basis), size(Uperp, 2))
+        I_nontrivial_standard = [
+            _index(Ubasis, std) for
+            std in standard_basis(b.dependence, trivial = false).monomials
+        ]
+        Uperpstd = Uperp[I_nontrivial_standard, :]
         for i in axes(new_matrix, 1)
             if iszero(I1[i])
                 @assert !iszero(I2[i])
-                new_matrix[i, :] = b.matrix[:, I2[i]]' * Uperp
+                new_matrix[i, :] = b.matrix[:, I2[i]]' * Uperpstd
             else
                 @assert iszero(I2[i])
                 new_matrix[i, :] = Uperp[I1[i], :]
             end
         end
         null = MacaulayNullspace(new_matrix, new_basis)
-        return solve(null, ShiftNullspace{StaircaseDependence}())
+        new_solver = StaircaseSolver{T}(;
+            max_partial_iterations = solver.max_partial_iterations - partial,
+            max_iterations = solver.max_iterations - !partial,
+            solver.rank_check,
+            solver.solver,
+        )
+        return solve(null, ShiftNullspace{StaircaseDependence}(new_solver))
     end
 end
 
@@ -249,6 +305,117 @@ function commutation_fix(matrices, ε)
     end
 end
 
+function partial_commutation_fix(
+    known_border_coefficients,
+    border_coefficients,
+    ::Type{T},
+    standard,
+    vars,
+    rank_check::RankCheck,
+) where {T}
+    function shifted_border_coefficients(mono, shift)
+        coef = border_coefficients(mono)
+        ret = zero(coef)
+        unknown = zero(MP.polynomial_type(mono, T))
+        for i in eachindex(coef)
+            if iszero(coef)
+                continue
+            end
+            shifted = shift * standard.monomials[i]
+            j = _index(standard, shifted)
+            if !isnothing(j)
+                ret[j] += coef[i]
+            elseif known_border_coefficients(shifted)
+                ret .+= coef[i] .* border_coefficients(shifted)
+            else
+                MA.operate!(MA.add_mul, unknown, coef[i], shifted)
+            end
+        end
+        return ret, unknown
+    end
+    new_relations = T[]
+    unknowns = MP.polynomial_type(prod(vars), T)[]
+    for std in standard.monomials
+        for x in vars
+            mono_x = x * std
+            if !known_border_coefficients(mono_x)
+                # FIXME what do we do if one of the two monomials is unknown
+                # but the other one is known ?
+                continue
+            end
+            for y in vars
+                mono_y = y * std
+                if !known_border_coefficients(mono_y)
+                    # FIXME what do we do if one of the two monomials is unknown
+                    # but the other one is known ?
+                    continue
+                end
+                if isnothing(_index(standard, mono_x))
+                    if isnothing(_index(standard, mono_y))
+                        coef_xy, unknowns_xy =
+                            shifted_border_coefficients(mono_y, x)
+                    else
+                        mono_xy = x * mono_y
+                        if known_border_coefficients(mono_xy)
+                            coef_xy = border_coefficients(mono_xy)
+                            unknowns_yx = zero(PT)
+                        else
+                            coef_xy = zeros(length(standard.monomials))
+                            unknowns_xy = mono_xy
+                        end
+                    end
+                    coef_yx, unknowns_yx =
+                        shifted_border_coefficients(mono_x, y)
+                else
+                    if !isnothing(_index(standard, mono_y))
+                        # Let `f` be `known_border_coefficients`.
+                        # They are both standard so we'll get
+                        # `f(y * mono_x) - f(x * mono_y)`
+                        # which will give a zero column, let's just ignore it
+                        continue
+                    end
+                    mono_yx = y * mono_x
+                    if known_border_coefficients(mono_yx)
+                        coef_yx = border_coefficients(mono_yx)
+                        unknowns_yx = zero(PT)
+                    else
+                        coef_yx = zeros(length(standard.monomials))
+                        unknowns_yx = mono_yx
+                    end
+                    coef_xy, unknowns_xy =
+                        shifted_border_coefficients(mono_y, x)
+                end
+                append!(new_relations, coef_xy - coef_yx)
+                push!(unknowns, unknowns_xy - unknowns_yx)
+            end
+        end
+    end
+    standard_part = reshape(
+        new_relations,
+        length(standard.monomials),
+        div(length(new_relations), length(standard.monomials)),
+    )
+    unknown_monos = MP.merge_monomial_vectors(MP.monomials.(unknowns))
+    unknown_part = Matrix{T}(undef, length(unknown_monos), length(unknowns))
+    for i in eachindex(unknowns)
+        unknown_part[:, i] = MP.coefficients(unknowns[i], unknown_monos)
+    end
+    basis, I1, I2 = MB.merge_bases(standard, MB.MonomialBasis(unknown_monos))
+    M = Matrix{T}(undef, length(basis.monomials), size(standard_part, 2))
+    for i in eachindex(basis.monomials)
+        if iszero(I1[i])
+            @assert !iszero(I2[i])
+            M[i, :] = unknown_part[I2[i], :]
+        else
+            @assert iszero(I2[i])
+            M[i, :] = standard_part[I1[i], :]
+        end
+    end
+    F = LinearAlgebra.svd(M, full = true)
+    r = rank_from_singular_values(F.S, rank_check)
+    return F.U[:, (r+1):end], basis
+end
+
 """
     Base.@kwdef struct AlgebraicBorderSolver{
         D,

src/shift.jl

@@ -103,7 +103,8 @@ 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{D,C<:RankCheck} <: MacaulayNullspaceSolver
+struct ShiftNullspace{D,S,C<:RankCheck} <: MacaulayNullspaceSolver
+    solver::S
     check::C
 end
 # Because the matrix is orthogonal, we know the SVD of the whole matrix is
@@ -112,14 +113,30 @@ end
 # constant monomial) should be a standard monomial, `LeadingRelativeRankTol`
 # ensures that we will take it.
 function ShiftNullspace{D}(check::RankCheck) where {D}
-    return ShiftNullspace{D,typeof(check)}(check)
+    return ShiftNullspace{D,Nothing,typeof(check)}(nothing, check)
+end
+function ShiftNullspace{D}(solver::StaircaseSolver) where {D}
+    check = solver.rank_check
+    return ShiftNullspace{D,typeof(solver),typeof(check)}(solver, check)
+end
+function ShiftNullspace{D}() where {D}
+    return ShiftNullspace{D}(LeadingRelativeRankTol(Base.rtoldefault(Float64)))
 end
-ShiftNullspace{D}() where {D} = ShiftNullspace{D}(LeadingRelativeRankTol(1e-8))
 ShiftNullspace(args...) = ShiftNullspace{StaircaseDependence}(args...)
 
+function _solver(
+    shift::ShiftNullspace{StaircaseDependence,Nothing},
+    ::Type{T},
+) where {T}
+    return StaircaseSolver{T}(rank_check = shift.check)
+end
+function _solver(shift::ShiftNullspace, ::Type)
+    return shift.solver
+end
+
 function border_basis_and_solver(
-    null::MacaulayNullspace,
+    null::MacaulayNullspace{T},
     shift::ShiftNullspace{D},
-) where {D}
-    return BorderBasis{D}(null, shift.check), nothing
+) where {T,D}
+    return BorderBasis{D}(null, shift.check), _solver(shift, T)
 end

test/null.jl

@@ -9,6 +9,24 @@ b(x) = MB.MonomialBasis(x)
 
 include("utils.jl")
 
+function test_partial_commutative_fix(x, y)
+    matrix = Float64[
+        1 0 0 0 # 1
+        0 1 0 0 # y
+        0 0 1 0 # x
+        0 0 0 1 # y^2
+        0 0 1 0 # x * y
+        1 0 0 0 # x^2
+    ]
+    basis = MB.MonomialBasis(MP.monomials((x, y), 0:2))
+    null = MM.MacaulayNullspace(matrix, basis, 1e-8)
+    D = MM.StaircaseDependence
+    solver = MM.StaircaseSolver{Float64}(max_partial_iterations = 1)
+    shift_solver = MM.ShiftNullspace{D}(solver)
+    sol = MM.solve(null, shift_solver)
+    return testelements(sol, [[1, 1], [-1, 1]], 1e-8)
+end
+
 function test_dependent_border(x, y)
     matrix = Float64[
         1 0 0 0 0 # 1