From 34fcc135249bca0c72e04f23b7bd7f680e4e0c54 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 10 Oct 2026 11:20:46 -0500 Subject: [PATCH] Raise the LSQR penalty cap to 1e6 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit On sparse nonsymmetric lognormal matrices `cover_transversal` was 1.5 to 29 times slower than `cover_min`, the ratio growing with size. Its multiplier iteration contracted by less than a factor of two per update once the penalty reached the LSQR cap of 1e5, stalled with a KKT residual of 1e-6 to 1e-4, and left the active-set finish tens to hundreds of steps. The cause is the shape of the active set. At the minimizer the active constraints of the transversal-tight problem form one large tree: on a 10000x10000 matrix with six entries per row, 2792 rows with diameter 163 and Laplacian Fiedler value 2.6e-5, against trees of at most 41 nodes and Fiedler value 1.3e-2 for `cover_min` on the same matrix. The method of multipliers contracts at a rate near 1/(1 + κ s), with s bounded by that eigenvalue, so κ = 1e5 barely contracts while 1e6 gains a factor of two to four per update. Neither problem is degenerate: both have a clean split between active and slack constraints and one positive multiplier per active constraint. With the cap at 1e6, `cover_transversal` certifies with 0 to 2 finish steps instead of 10 to 119 on four such matrices (three at n = 10000, one at n = 31623), in equal or shorter time; at n = 31623 the time drops from 14.2 s to 8.8 s. `cover_min` is unchanged within three solves. At 1e7 the iteration levels off near 1e-7, as does the exact dense path, so 1e6 is the usable ceiling in Float64. The drain test moves to a banded matrix whose slack multipliers still need a penalty above the new cap. The cross-check of the symmetric finish against the asymmetric solve loosens to 1e-9: at the higher cap the asymmetric multiplier iteration can meet its own KKT tolerance and skip the finish, which leaves its scaled entries about 1e-10 from the certified minimizer. Assisted-by: Claude Fable 5.1 --- src/minimal_covers.jl | 4 ++-- test/minimal_covers.jl | 20 +++++++++++--------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 61008ff..b056414 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -1806,9 +1806,9 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; # Absolute thresholds preserve covariance under rescaling because the # residuals are logarithmic. vtol = tol === nothing ? 1000 * eps(T) : T(tol) - # Precision bounds the usable penalty: the cap is `1e5` (LSQR) or `1e8` + # Precision bounds the usable penalty: the cap is `1e6` (LSQR) or `1e8` # (exact solves) in `Float64` and grows as `1/sqrt(eps(T))`. - κcap = T(use_lsqr ? 1e5 : 1e8) * sqrt(T(eps(Float64)) / eps(T)) + κcap = T(use_lsqr ? 1e6 : 1e8) * sqrt(T(eps(Float64)) / eps(T)) # Slack below the inner solver's resolution is indistinguishable from an # active entry. zdrain = sqrt(eps(T)) diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 596e7ce..85834c2 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -582,7 +582,9 @@ end # cover. With either preconditioner it certifies the minimizer, which is also # the minimizer of the asymmetric problem. (Whether the full multiplier # iteration converges on its own here depends on rounding, so the test - # does not rely on it.) + # does not rely on it. When the asymmetric solve does stop at the multiplier + # iteration's KKT tolerance, its scaled entries are within about 1e-10 of the + # certified minimizer rather than at roundoff.) an, sn = @test_logs (:warn, r"Increase `maxouter` or `κ`") MatrixCovers._symcover_min_abslog2(A; flopbudget=Inf, maxouter=3, polish=false) @test !sn.converged @@ -592,10 +594,10 @@ end ap, sp = MatrixCovers._symcover_min_abslog2(A; flopbudget=budget, maxouter=3) @test sp.converged && sp.polish.certified @test cover_objective(AbsLog{2}(), ap, ap, A) <= cover_objective(AbsLog{2}(), an, an, A) - @test ap[Ai] .* ap[Aj] ≈ ga[Ai] .* gb[Aj] rtol=1e-12 + @test ap[Ai] .* ap[Aj] ≈ ga[Ai] .* gb[Aj] rtol=1e-9 @test cover_objective(AbsLog{2}(), ap, ap, A) ≈ cover_objective(AbsLog{2}(), ga, gb, A) rtol=1e-12 end - @test af[Ai] .* af[Aj] ≈ ga[Ai] .* gb[Aj] rtol=1e-12 + @test af[Ai] .* af[Aj] ≈ ga[Ai] .* gb[Aj] rtol=1e-9 # The forest factor needs no budget: any budget below the full factor's # gives the forest regime for every solve, and the cover is unchanged. @@ -884,7 +886,7 @@ end # The LSQR path caps escalation where its accuracy floor makes larger # weights useless. _, sl = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) - @test sl.κs[end] <= 1e5 + @test sl.κs[end] <= 1e6 @test sl.kkt[end] < 1e-8 # Stopping far short of convergence without the active-set finish emits a warning. @@ -930,14 +932,14 @@ end end @testset "MMC :lsqr drains slack multipliers at the penalty cap" begin - # Draining these multipliers needs a penalty above the `:lsqr` cap of 1e5. + # Draining these multipliers needs a penalty above the `:lsqr` cap of 1e6. # With the penalty held at the cap they are zeroed directly; left to drain by - # `2(κ-1)z` per update they hold the KKT residual at 2.6e-6. + # `2(κ-1)z` per update they hold the KKT residual at 6.2e-8. # Without the active-set finish (`polish=false`), slack multipliers are zeroed. # Whether the residual then reaches the tolerance depends on rounding (with the # forest preconditioner it can stall just above it); the solve warns exactly # when it does not. The default solve's finish certifies the result. - A = banded_sparse_sym(StableRNG(5), 10_000, 3, 1.0) + A = banded_sparse_sym(StableRNG(1), 10_000, 3, 1.0) for fillbudget in (MatrixCovers.LSQR_FILL_BUDGET, 0) (a, s), logs = capture_logs() do MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, fillbudget, polish=false) @@ -946,11 +948,11 @@ end @test all(l -> l.level == Base.CoreLogging.Warn && occursin("try a larger `κ`", l.message), logs) @test s.precond === (fillbudget == 0 ? :forest : :factor) @test s.nzeroed > 0 - @test maximum(s.κs) == 1e5 + @test maximum(s.κs) == 1e6 @test s.kkt[end] < 1e-9 _, sp = @test_logs symcover_min(AbsLog{2}(), A; fillbudget) @test MatrixCovers.converged(sp) - @test cover_objective(AbsLog{2}(), a, A) ≈ 161350.449205 rtol=1e-9 + @test cover_objective(AbsLog{2}(), a, A) ≈ 163679.149704 rtol=1e-9 @test iscover(a, A) end