Fix tolerance for commutation fix (#88)

* Fix tolerance for commutation fix

* Fix

* Fix

* Fix format

src/border.jl

@@ -71,18 +71,21 @@ struct StaircaseSolver{
     max_iterations::Int
     rank_check::R
     solver::M
+    print_level::Int
 end
 function StaircaseSolver{T}(;
     max_partial_iterations::Int = 0,
     max_iterations::Int = -1,
     rank_check::RankCheck = LeadingRelativeRankTol(Base.rtoldefault(T)),
     solver = SS.ReorderedSchurMultiplicationMatricesSolver{T}(),
+    print_level::Int = 1,
 ) where {T}
     return StaircaseSolver{T,typeof(rank_check),typeof(solver)}(
         max_partial_iterations,
         max_iterations,
         rank_check,
         solver,
+        print_level,
     )
 end
 
@@ -225,7 +228,11 @@ function solve(
         if solver.max_iterations == 0
             Uperp = nothing
         else
-            Uperp = commutation_fix(mult, solver.solver.ε)
+            Uperp = commutation_fix(
+                mult,
+                √solver.solver.ε;
+                print_level = solver.print_level,
+            )
         end
         com_fix = if isnothing(Uperp)
             nothing
@@ -267,7 +274,7 @@ function solve(
     end
 end
 
-function commutation_fix(matrices, ε)
+function commutation_fix(matrices, ε; print_level)
     if isempty(first(matrices))
         return
     end
@@ -301,6 +308,11 @@ function commutation_fix(matrices, ε)
     if isnothing(leading_F)
         return
     else
+        if print_level >= 1
+            @info(
+                "Found multiplication matrices with commutation error of `$leading_ratio` which is larger than the tolerance of `$ε`. Adding this to the equations and continuing."
+            )
+        end
         return leading_F.U[:, 2:end]
     end
 end

src/shift.jl

@@ -106,37 +106,53 @@ IFAC Proceedings Volumes 45.16 (2012): 1203-1208.
 struct ShiftNullspace{D,S,C<:RankCheck} <: MacaulayNullspaceSolver
     solver::S
     check::C
+    print_level::Int
 end
 # Because the matrix is orthogonal, we know the SVD of the whole matrix is
 # `ones(...)` so an `AbsoluteRankTol` would be fine here.
 # 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.
-function ShiftNullspace{D}(check::RankCheck) where {D}
-    return ShiftNullspace{D,Nothing,typeof(check)}(nothing, check)
+function ShiftNullspace{D}(check::RankCheck; print_level = 1) where {D}
+    return ShiftNullspace{D,Nothing,typeof(check)}(nothing, check, print_level)
 end
-function ShiftNullspace{D}(solver::StaircaseSolver) where {D}
+function ShiftNullspace{D}(solver::StaircaseSolver; print_level = 1) where {D}
     check = solver.rank_check
-    return ShiftNullspace{D,typeof(solver),typeof(check)}(solver, check)
+    return ShiftNullspace{D,typeof(solver),typeof(check)}(
+        solver,
+        check,
+        print_level,
+    )
 end
-function ShiftNullspace{D}() where {D}
-    return ShiftNullspace{D}(LeadingRelativeRankTol(Base.rtoldefault(Float64)))
+function ShiftNullspace{D}(; kws...) where {D}
+    return ShiftNullspace{D}(
+        LeadingRelativeRankTol(Base.rtoldefault(Float64));
+        kws...,
+    )
 end
 ShiftNullspace(args...) = ShiftNullspace{StaircaseDependence}(args...)
 
 function _solver(
+    null::MacaulayNullspace,
     shift::ShiftNullspace{StaircaseDependence,Nothing},
     ::Type{T},
 ) where {T}
-    return StaircaseSolver{T}(rank_check = shift.check)
+    return StaircaseSolver{T}(;
+        rank_check = shift.check,
+        solver = SS.ReorderedSchurMultiplicationMatricesSolver(
+            convert(T, null.accuracy),
+        ),
+        print_level = shift.print_level,
+    )
 end
-function _solver(shift::ShiftNullspace, ::Type)
+
+function _solver(::MacaulayNullspace, shift::ShiftNullspace, ::Type)
     return shift.solver
 end
 
 function border_basis_and_solver(
     null::MacaulayNullspace{T},
     shift::ShiftNullspace{D},
 ) where {T,D}
-    return BorderBasis{D}(null, shift.check), _solver(shift, T)
+    return BorderBasis{D}(null, shift.check), _solver(null, shift, T)
 end

test/extract.jl

@@ -362,6 +362,29 @@ function large_norm(rank_check)
     @test nothing === atomic_measure(ν, rank_check)
 end
 
+function no_com(solver, T)
+    Mod.@polyvar x[1:2]
+    Q = T[
+        1.0000000000000002 2.9999999777389235 1.408602019734018 8.999999911981313 4.225805987406138 3.043010098670093
+        2.9999999777389235 8.999999911981313 4.225805987406138 26.999999824988187 12.677417999642149 9.129030026074997
+        1.408602019734018 4.225805987406138 3.043010098670093 12.677417999642149 9.129030026074997 9.580642414414385
+        8.999999911981313 26.999999824988187 12.677417999642149 80.99999955990646 38.03225428611012 27.387090350285558
+        4.225805987406138 12.677417999642149 9.129030026074997 38.03225428611012 27.387090350285558 28.74192618075039
+        3.043010098670093 9.129030026074997 9.580642414414385 27.387090350285558 28.74192618075039 35.73117167739159
+    ]
+    ν = moment_matrix(Q, monomials(x, 0:2))
+    η = AtomicMeasure(
+        x,
+        [
+            WeightedDiracMeasure(T[1, 3], T(0.863799337)),
+            WeightedDiracMeasure(T[4, 3], T(0.136200673)),
+        ],
+    )
+    rank_check = LeadingRelativeRankTol(1e-7)
+    atoms = atomic_measure(ν, rank_check, solver)
+    @test atoms ≈ η rtol = 1e-4
+end
+
 _short(x::String) = x[1:min(10, length(x))]
 _short(x) = _short(string(x))
 
@@ -443,6 +466,12 @@ function test_extract()
     jcg14_6_1(6e-3)
     jcg14_6_1(8e-4, false)
     large_norm(1e-2)
+    @testset "No comm fix $(_short(solver)) $T" for solver in
+                                                    [ShiftNullspace()],
+        T in [Float64]
+        #, BigFloat]
+        no_com(solver, T)
+    end
     return
 end