Merge pull request #570 from JuliaHomotopyContinuation/scaling_bug

Fixing a bug in row scaling

src/endgame_tracker.jl

@@ -40,6 +40,7 @@ These parameters control the behaviour during the endgame.
 * `val_finite_tol = 1e-3`: Tolerance on the valuation which has to be satisfied before the endgame is started.
 * `sing_cond = 1e14`: value for the condition number above which a solution is considered singular.
 * `sing_accuracy = 1e-12`: value for the accuracy number above which a solution is considered singular.
+* `scaling_threshold = -30.0`: Row scaling of matrices is only applied to rows with `e < scaling_threshold`, where is the norm of the row is estimated to be `2^e`. See [`skeel_row_scaling`](@skeel_row_scaling) for details.
 * `refine_steps = 3`: number of steps for refining solutions at the end.
 """
 Base.@kwdef mutable struct EndgameOptions
@@ -63,6 +64,7 @@ Base.@kwdef mutable struct EndgameOptions
     # singular solutions parameters
     sing_cond::Float64 = 1e14
     sing_accuracy::Float64 = 1e-12
+    scaling_threshold::Float64 = -30.0
 
     # refinement parameters
     refine_steps::Int = 3
@@ -386,7 +388,7 @@ function step!(endgame_tracker::EndgameTracker, debug::Bool = false)
     update!(state.val, tracker.predictor, t)
     debug && println(state.val)
     if check_finite!(state, options)
-        switch_to_singular!(state, tracker; debug = debug)
+        switch_to_singular!(state, tracker, options; debug = debug)
         return state.code
     elseif check_at_infinity!(state, tracker, options; debug = debug)
         return state.code
@@ -442,6 +444,7 @@ function check_at_infinity!(state, tracker, options; debug::Bool = false)
                         state.row_scaling,
                         workspace(tracker.state.jacobian),
                         state.col_scaling,
+                        options.scaling_threshold,
                     )
                 end
                 κ = LA.cond(tracker.state.jacobian, state.row_scaling, state.col_scaling)
@@ -494,7 +497,7 @@ function check_at_infinity!(state, tracker, options; debug::Bool = false)
     false
 end
 
-function switch_to_singular!(state, tracker; debug::Bool)
+function switch_to_singular!(state, tracker, options; debug::Bool)
     state.singular_endgame = true
     t = real(tracker.state.t)
     state.singular_start = t
@@ -506,6 +509,7 @@ function switch_to_singular!(state, tracker; debug::Bool)
             state.row_scaling,
             workspace(tracker.state.jacobian),
             state.col_scaling,
+            options.scaling_threshold,
         )
     end
     add_sample!(state, tracker, t)
@@ -703,6 +707,7 @@ function tracking_stopped!(endgame_tracker::EndgameTracker)
             state.row_scaling,
             workspace(tracker.state.jacobian),
             state.col_scaling,
+            options.scaling_threshold,
         )
         state.cond =
             LA.cond(tracker, state.solution, 0.0, state.row_scaling, state.col_scaling)

src/linear_algebra.jl

@@ -418,19 +418,21 @@ end
 
 """
     skeel_row_scaling!(W::MatrixWorkspace, c)
-    skeel_row_scaling!(d, A, c)
+    skeel_row_scaling!(d, A, c;
+    scaling_threshold = -30.0)
 
 Compute optimal scaling factors `d` for the matrix `A` following Skeel [^S79]
 if `c` is approximately of the order of the solution of the linear system
 of interest.
-The scaling factors are rounded to powers of the base radix 2.
+The scaling factors are rounded to powers `e` of the base radix 2. Row scaling is only applied to rows with `e - m ≥ scaling_threshold` (this is an inequality of norms at log-scale) to prevent zero rows from being scaled, where `2^m` is approximately the maximum norm of a row of `A`. 
 
 [^S79]: Skeel, Robert D. "Scaling for numerical stability in Gaussian elimination." Journal of the ACM (JACM) 26.3 (1979): 494-526.
 """
 function skeel_row_scaling!(
     d::AbstractVector{<:Real},
     A::AbstractMatrix{<:Complex},
-    c::AbstractVector{<:Real},
+    c::AbstractVector{<:Real};
+    scaling_threshold::Float64 = -30.0,
 )
     n = length(c)
     @inbounds d .= zero(eltype(d))
@@ -440,8 +442,16 @@ function skeel_row_scaling!(
             d[i] += fast_abs(A[i, j]) * cj
         end
     end
+
+    m = maximum(d)
+    s = scaling_threshold + m
     @inbounds for i = 1:n
-        d[i] = exp2(-last(frexp(d[i])))
+        e = last(frexp(d[i]))
+        if e < s
+            d[i] = 1.0
+        else
+            d[i] = exp2(-e)
+        end
     end
 
     d
@@ -456,10 +466,11 @@ function row_scaling!(
     d::AbstractVector{<:Real},
     WS::MatrixWorkspace,
     c::AbstractVector{<:Real},
+    scaling_threshold::Float64,
 )
     m, n = size(WS)
     if m == n
-        skeel_row_scaling!(d, WS.A, c)
+        skeel_row_scaling!(d, WS.A, c; scaling_threshold = scaling_threshold)
     else
         d .= 1.0
     end
@@ -754,6 +765,7 @@ function LA.cond(
             rmax = max(rmax, rᵢ)
             rmin = min(rmin, rᵢ)
         end
+        #@show (rmax, rmin)
         rmax / rmin
     else
         inverse_inf_norm_est(WS, d_l, d_r) * inf_norm(WS, d_l, d_r)

src/numerical_irreducible_decomposition.jl

@@ -449,7 +449,7 @@ function regeneration!(
     endgame_options = EndgameOptions(;
         max_endgame_steps = 100,
         max_endgame_extended_steps = 100,
-        sing_accuracy = 1e-10,
+        sing_cond = 1e12,
     ),
     threading::Bool = true,
     seed = rand(UInt32),

test/nid_test.jl

@@ -28,6 +28,9 @@
         N = nid(F; seed = nothing)
         @test isnothing(seed(N))
 
+        # bad seed
+        N = nid(F; seed = 0xc770fa47)
+
         # progress
         N = nid(F; show_progress = false)
         @test isa(N, NumericalIrreducibleDecomposition)