From 13e61405f7393241f1c2eede1369247b88702f01 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 23 Aug 2026 15:57:56 -0500 Subject: [PATCH 1/4] Default soft AbsLog{2} minimizers to sparse LSQR soft_symcover_min and soft_cover_min used :auto, which could select dense normal equations for sparse matrices. Give their sparse-storage methods the :lsqr default used by the hard minimal-cover solvers. The other soft AbsLog entry points delegate to these methods. Add Unitful disambiguation methods for allocating sparse minimal-cover calls. Without them, sparse quantity matrices lost unit-stripping dispatch and could fail at log(::Quantity). Assisted-by: Claude Fable 5 --- ext/MatrixCoversUnitfulExt.jl | 14 +++++++++++++- src/sparse_support.jl | 13 +++++++++++++ test/soft_covers.jl | 20 ++++++++++++++++++++ test/unitful.jl | 7 ++++++- 4 files changed, 52 insertions(+), 2 deletions(-) diff --git a/ext/MatrixCoversUnitfulExt.jl b/ext/MatrixCoversUnitfulExt.jl index accc42b..8fc724a 100644 --- a/ext/MatrixCoversUnitfulExt.jl +++ b/ext/MatrixCoversUnitfulExt.jl @@ -274,12 +274,24 @@ MC.soft_symcover_min!(a::QVector, A::QMatrix; kwargs...) = symstart!(MC.soft_sym MC.soft_cover_min(A::QMatrix; kwargs...) = asym(MC.soft_cover_min, A; kwargs...) MC.soft_cover_min!(a::QVector, b::QVector, A::QMatrix; kwargs...) = asymstart!(MC.soft_cover_min!, a, b, A; kwargs...) -# Resolve the overlap between sparse refiner and unitful matrix methods. +# Disambiguate sparse unitful matrices without losing sparse storage. const QSparse = SparseMatrixCSC{<:Quantity} const QSparseSym = Union{QSparse, Symmetric{<:Quantity,<:SparseMatrixCSC}, Hermitian{<:Quantity,<:SparseMatrixCSC}} +MC.symcover_min(ϕ::AbsLog{2}, A::QSparseSym; kwargs...) = + sym(MC.symcover_min, A, ϕ; kwargs...) + +MC.cover_min(ϕ::AbsLog{2}, A::QSparse; kwargs...) = + asym(MC.cover_min, A, ϕ; kwargs...) + +MC.soft_symcover_min(ϕ::AbsLog{2}, A::QSparseSym; kwargs...) = + sym(MC.soft_symcover_min, A, ϕ; kwargs...) + +MC.soft_cover_min(ϕ::AbsLog{2}, A::QSparse; kwargs...) = + asym(MC.soft_cover_min, A, ϕ; kwargs...) + MC.symcover_min!(ϕ::AbsLog{2}, a::QVector, A::QSparseSym; kwargs...) = symstart!(MC.symcover_min!, a, A, ϕ; kwargs...) diff --git a/src/sparse_support.jl b/src/sparse_support.jl index e863ef4..6836ebb 100644 --- a/src/sparse_support.jl +++ b/src/sparse_support.jl @@ -117,3 +117,16 @@ function cover_min!(ϕ::AbsLog{2}, a::AbstractVector, b::AbstractVector, b .= bnew return a, b end + +# Soft minimizers use the same sparse default. +function soft_symcover_min(ϕ::AbsLog{2}, + S::Union{SparseMatrixCSC,Symmetric{<:Any,<:SparseMatrixCSC},Hermitian{<:Any,<:SparseMatrixCSC}}; + linsolve::Symbol=:lsqr, kwargs...) + a, _ = _soft_symcover_min_abslog2(S; linsolve, kwargs...) + return a +end + +function soft_cover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) + a, b, _ = _soft_cover_min_abslog2(A; linsolve, kwargs...) + return a, b +end diff --git a/test/soft_covers.jl b/test/soft_covers.jl index 5663114..bdc3777 100644 --- a/test/soft_covers.jl +++ b/test/soft_covers.jl @@ -411,6 +411,26 @@ end @test MatrixCovers._soft_cover_min_abslog2(Y)[3].linsolve === :woodbury end +@testset "sparse soft AbsLog{2} defaults to matrix-free LSQR" begin + rng = StableRNG(19) + S0 = sprand(rng, 60, 60, 0.1) + Ssp = S0 + S0' + I # supported diagonal keeps every component non-bipartite + a = soft_symcover_min(AbsLog{2}(), Ssp) + @test a isa Vector{Float64} + @test a == soft_symcover_min(AbsLog{2}(), Ssp; linsolve=:lsqr) + @test a ≈ soft_symcover_min(AbsLog{2}(), Matrix(Ssp); linsolve=:dense) rtol=1e-6 + U = Symmetric(sparse(triu(Ssp))) + @test soft_symcover_min(AbsLog{2}(), U) == soft_symcover_min(AbsLog{2}(), U; linsolve=:lsqr) + # Asymmetric form; the cover products are gauge-free. + G = sprand(rng, 50, 40, 0.12) + ga, gb = soft_cover_min(AbsLog{2}(), G) + @test (ga, gb) == soft_cover_min(AbsLog{2}(), G; linsolve=:lsqr) + gad, gbd = soft_cover_min(AbsLog{2}(), Matrix(G); linsolve=:dense) + @test ga .* gb' ≈ gad .* gbd' rtol=1e-6 + a0 = copy(a) + @test soft_symcover_min!(AbsLog{2}(), a0, Ssp) == a +end + @testset "soft_symcover!/soft_cover! refiners" begin Asym = [4.0 1.0 0.5; 1.0 3.0 1.0; 0.5 1.0 2.5] Agen = [1.0 2.0 0.5; 0.25 3.0 1.0] diff --git a/test/unitful.jl b/test/unitful.jl index 3cde2c9..2dfb6dd 100644 --- a/test/unitful.jl +++ b/test/unitful.jl @@ -93,7 +93,12 @@ @test all(unit.(v) == fill(u"m^-1", 3) for v in cover(S)) @test unit.(symcover(Symmetric(S))) == fill(u"m^-1", 3) - # The refiners are where the sparse methods and MatrixCoversUnitfulExt overlap. + @test unit.(symcover_min(AbsLog{2}(), S)) == fill(u"m^-1", 3) + @test unit.(symcover_min(AbsLog{2}(), Symmetric(S))) == fill(u"m^-1", 3) + @test all(unit.(v) == fill(u"m^-1", 3) for v in cover_min(AbsLog{2}(), S)) + @test unit.(soft_symcover_min(AbsLog{2}(), S)) == fill(u"m^-1", 3) + @test unit.(soft_symcover_min(AbsLog{2}(), Symmetric(S))) == fill(u"m^-1", 3) + @test all(unit.(v) == fill(u"m^-1", 3) for v in soft_cover_min(AbsLog{2}(), S)) a = initialize_symcover(S) @test unit.(symcover_min!(AbsLog{2}(), a, S)) == fill(u"m^-1", 3) @test unit.(symcover_min!(AbsLog{2}(), a, Symmetric(S))) == fill(u"m^-1", 3) From 5235045d74e1258030afd178f90f5705ddd10ecc Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 23 Aug 2026 15:59:35 -0500 Subject: [PATCH 2/4] Traverse gramcover weights by stored support Replace three full row-pair scans of W with foreach_support. Sparse weights now take O(nnz(W)) traversal instead of O(m^2 log m); dense behavior is unchanged. Count only stored nonzero summands when inflating for roundoff, since zero terms add no error. Assisted-by: Claude Fable 5 --- src/gram_covers.jl | 56 +++++++++++++++++++-------------------------- test/gram_covers.jl | 16 +++++++++++++ 2 files changed, 40 insertions(+), 32 deletions(-) diff --git a/src/gram_covers.jl b/src/gram_covers.jl index 0bca308..dbbf117 100644 --- a/src/gram_covers.jl +++ b/src/gram_covers.jl @@ -181,7 +181,9 @@ function gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc: throw(DimensionMismatch("`W` couples support rows, so it must be square on the row axis: axes(W) must be $(string((sc.rowax, sc.rowax))), got $(string(axes(W)))")) ncomp = ncomponents(sc) - # Merge support components coupled by nonzero entries of `W`. + # Merge components coupled by `W`. Entries involving unsupported rows are + # ignored: their `a`-scales lie outside the cover, so summing them would + # corrupt the bound. parent = collect(1:ncomp) function find(p) while parent[p] != p @@ -190,33 +192,28 @@ function gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc: end return p end - merged = false - for i in sc.rowax + mergedref = Ref(false) + foreach_support(W) do i, ip, _ ci = rowcomponent(sc, i) - iszero(ci) && continue - for ip in sc.rowax - cip = rowcomponent(sc, ip) - (iszero(cip) || ci == cip) && continue - iszero(abs(W[i, ip])) && continue - ri, rip = find(ci), find(cip) - ri == rip && continue - parent[ri] = rip - merged = true - end + cip = rowcomponent(sc, ip) + (iszero(ci) || iszero(cip) || ci == cip) && return + ri, rip = find(ci), find(cip) + ri == rip && return + parent[ri] = rip + mergedref[] = true end # Avoid group bookkeeping when all components remain independent. - if !merged + if !mergedref[] m = zeros(typeof(_gc_term(a, W)), ncomp) n = zeros(Int, ncomp) - for i in sc.rowax + # Without merges, each supported entry stays within one component. + foreach_support(W) do i, ip, v ci = rowcomponent(sc, i) - iszero(ci) && continue - for ip in sc.rowax - rowcomponent(sc, ip) == ci || continue - m[ci] += a[i] * abs(W[i, ip]) * a[ip] - n[ci] += 1 - end + iszero(ci) && return + rowcomponent(sc, ip) == ci || return + m[ci] += a[i] * v * a[ip] + n[ci] += 1 end return _write_gramcover!(s, b, sc.colcomp, first(sc.colax) - 1, m, n) end @@ -244,19 +241,14 @@ function gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc: M[r] = zeros(T, k, k) nterm[r] = zeros(Int, k, k) end - for i in sc.rowax + # Every supported entry now has endpoints with the same root. + foreach_support(W) do i, ip, v ci = rowcomponent(sc, i) - iszero(ci) && continue + cip = rowcomponent(sc, ip) + (iszero(ci) || iszero(cip)) && return r = find(ci) - p = local_idx[ci] - for ip in sc.rowax - cip = rowcomponent(sc, ip) - iszero(cip) && continue - find(cip) == r || continue - q = local_idx[cip] - M[r][p, q] += a[i] * abs(W[i, ip]) * a[ip] - nterm[r][p, q] += 1 - end + M[r][local_idx[ci], local_idx[cip]] += a[i] * v * a[ip] + nterm[r][local_idx[ci], local_idx[cip]] += 1 end sq = _gc_group_scales(members, M, nterm, degenerate) diff --git a/test/gram_covers.jl b/test/gram_covers.jl index fd1f061..a3f140f 100644 --- a/test/gram_covers.jl +++ b/test/gram_covers.jl @@ -171,6 +171,22 @@ @test all(sw * sw' .>= abs.(Matrix(Js)' * Diagonal(w) * Matrix(Js))) end + @testset "sparse W matches its dense reading" begin + rng = StableRNG(21) + B = randn(rng, 4, 3); C = randn(rng, 3, 2) + J = [B zeros(4, 2); zeros(3, 3) C] # two support components + a, b = cover(J) + m = size(J, 1) + Wsp = sparse(1.0I, m, m) + Wsp[1, 5] = Wsp[5, 1] = 0.5 + s = gramcover(a, b, J, Wsp) + @test s == gramcover(a, b, J, Matrix(Wsp)) + @test all(s * s' .>= abs.(J' * Matrix(Wsp) * J)) + # Uncoupled weights. + Dsp = sparse(2.0I, m, m) + @test gramcover(a, b, J, Dsp) == gramcover(a, b, J, Matrix(Dsp)) + end + @testset "empty column" begin J = [1.0 0.0; 2.0 0.0; 0.0 0.0] a, b = cover(J) From a25d7a4cbd7108d5e6376ee13777839081bb47c9 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 23 Aug 2026 16:01:59 -0500 Subject: [PATCH 3/4] Resolve deferred diagfeasible pairs by worklist boost_feasible_seq! repeatedly rescanned deferred pairs, taking Theta(D^2) when scales propagated against scan order. Queue pairs incident to each newly scaled vertex instead. Each vertex resolves once, reducing deferred resolution to O(D) and leaving the sort as the O(nnz log nnz) term. Pairs unreachable from a scaled endpoint retain the ordered equal-split fallback. Worklist order can change heuristic outputs when multiple neighbors could supply a scale, but every result remains feasible. Assisted-by: Claude Fable 5 --- src/heuristic_covers.jl | 54 +++++++++++++++++++++++++---------------- test/initializers.jl | 22 +++++++++++++++++ 2 files changed, 55 insertions(+), 21 deletions(-) diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index 74b74ac..24d4ee4 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -496,9 +496,7 @@ function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix return a, b end -# Sequential nearest-neighbor feasibility propagation in increasing diagonal -# offset. Zero scales are unresolved; deferred pairs are revisited, then split -# equally if neither endpoint acquires a scale. The method costs O(n²). +# Sequential nearest-neighbor feasibility propagation by diagonal offset. function boost_feasible_seq!(a::AbstractVector{T}, A::AbstractMatrix) where T ax = eachindex(a) axes(A) == (ax, ax) || throw(DimensionMismatch("`boost_feasible_seq!(a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(string(axes(A))), axes(a)=$(string(axes(a))))")) @@ -530,31 +528,46 @@ function boost_feasible_seq!(a::AbstractVector{T}, A::AbstractMatrix) where T end end - # Resolve deferred constraints: re-scan until no more progress, then equal-split. - while !isempty(deferred) - changed = false - filter!(deferred) do (k, l, v) + # Each vertex is assigned at most once, so worklist resolution is linear in + # the number of deferred pairs. + if !isempty(deferred) + o = first(ax) - 1 + inc = [Int[] for _ in eachindex(ax)] # deferred pairs touching each vertex + for (e, (k, l, _)) in enumerate(deferred) + push!(inc[k-o], e) + push!(inc[l-o], e) + end + done = falses(length(deferred)) + queue = collect(eachindex(deferred)) + qi = firstindex(queue) + while qi <= lastindex(queue) + e = queue[qi] + qi += 1 + done[e] && continue + k, l, v = deferred[e] ak, al = a[k], a[l] - if !iszero(ak) && !iszero(al) + if iszero(ak) && iszero(al) + continue + elseif iszero(al) + a[l] = v / ak + done[e] = true + append!(queue, inc[l-o]) + elseif iszero(ak) + a[k] = v / al + done[e] = true + append!(queue, inc[k-o]) + else aprod = ak * al if aprod < v s = sqrt(v / aprod) a[k] *= s; a[l] *= s end - elseif !iszero(ak) - a[l] = v / ak - elseif !iszero(al) - a[k] = v / al - else - return true # still unresolvable; keep in list + done[e] = true end - changed = true - return false # resolved; drop from list end - changed && continue - # No progress: all remaining have both indices zero. - # Process in order so earlier equal-splits can inform later ones in the same pass. - for (k, l, v) in deferred + # Split components that no diagonal scale reached, preserving order. + for (e, (k, l, v)) in enumerate(deferred) + done[e] && continue ak, al = a[k], a[l] if iszero(ak) && iszero(al) a[k] = a[l] = sqrt(v) @@ -570,7 +583,6 @@ function boost_feasible_seq!(a::AbstractVector{T}, A::AbstractMatrix) where T end end end - break end return a diff --git a/test/initializers.jl b/test/initializers.jl index 470d2ab..e0feef0 100644 --- a/test/initializers.jl +++ b/test/initializers.jl @@ -88,6 +88,28 @@ @test ah[3] < ag[3] end + @testset ":diagfeasible propagates from sparse anchors" begin + # Propagate the diagonal scale backward through a path. + n = 6 + P = SymTridiagonal([zeros(n - 1); 4.0], fill(2.0, n - 1)) + a = initialize_symcover(P; strategy=:diagfeasible, feasible=:none) + @test iscover(a, P) + @test a[n] == 2.0 + @test all(a[i] * a[i+1] == 2.0 for i in 1:n-1) + + # A component with no diagonal at all falls back to equal splits. + Q = [0.0 9.0; 9.0 0.0] + aq = initialize_symcover(Q; strategy=:diagfeasible, feasible=:none) + @test iscover(aq, Q) + @test aq == [3.0, 3.0] + + # Guard against quadratic rescans. + n = 100_000 + P = SymTridiagonal([zeros(n - 1); 4.0], fill(2.0, n - 1)) + a = initialize_symcover(P; strategy=:diagfeasible, feasible=:none) + @test iscover(a, P) + end + @testset "no penalty argument" begin # Initializers depend on `A`, not on a penalty. @test_throws MethodError initialize_symcover(AbsLog{2}(), Asyms[1]) From fe38d5c5a5de00856d09cb77871329634564cfa8 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 23 Aug 2026 16:06:39 -0500 Subject: [PATCH 4/4] Choose linsolve=:auto by measured support density Outside sparse storage, :auto used dense normal equations whenever Woodbury was unavailable. Structured and wrapped sparse-support matrices could therefore incur O(N^2) memory and O(N^3) factorization costs. Select Woodbury when eligible, LSQR at support density <= 1/4, and dense solves otherwise. The choice depends on foreach_support rather than storage type. Results near the threshold may differ at solver tolerance. Assisted-by: Claude Fable 5 --- src/minimal_covers.jl | 19 +++++++++++++++++-- src/sparse_support.jl | 3 +-- test/minimal_covers.jl | 27 ++++++++++++++++++++++++++- 3 files changed, 44 insertions(+), 5 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index c687b6f..68748da 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -32,7 +32,8 @@ The native solver accepts `κs` (penalty-continuation schedule), `maxiter` requires at most `n ÷ 4` missing entries per row and `4n` in total. - `:lsqr` is matrix-free with O(nnz) work per iteration and is the sparse-matrix default. -- `:auto` chooses `:woodbury` when supported and `:dense` otherwise. +- `:auto` chooses `:woodbury` when supported, `:lsqr` when the stored support + fills at most a quarter of the grid, and `:dense` otherwise. The native solver computes in `Float64` for narrower input types, then converts the result to the required element type. @@ -253,7 +254,13 @@ end # It is not interchangeable with CG on the normal equations: LSQR's accuracy # tracks the condition number of `M` (≈ √κ), CG's that of `MᵀM` (≈ κ), and at # κ = 1e8 the latter exhausts double precision. -# - `:auto` selects `:woodbury` when supported and `:dense` otherwise. +# - `:auto` selects `:woodbury` when supported, `:lsqr` when the stored support +# fills at most `AUTO_LSQR_MAX_DENSITY` of the grid, and `:dense` otherwise. + +# Maximum support density for the `:auto` LSQR path. At or above it the exact +# dense solve is the better bargain: it terminates a stage on a sign-stable +# Newton step and has smaller constants. +const AUTO_LSQR_MAX_DENSITY = 1 // 4 # Condition estimate above which Woodbury uses sparse Cholesky instead of CG. const WOODBURY_CG_KAPPA = 1000 @@ -1058,6 +1065,10 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end use_woodbury = ok end + if linsolve === :auto && !use_woodbury && nsupp <= AUTO_LSQR_MAX_DENSITY * (n * n) + use_lsqr = true + linsolve = :lsqr + end # Woodbury uses a grid; dense and LSQR use an edge list. supp = if use_woodbury C = fill(T(-Inf), n, n) @@ -1172,6 +1183,10 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end use_woodbury = ok end + if linsolve === :auto && !use_woodbury && ne <= AUTO_LSQR_MAX_DENSITY * (m * n) + use_lsqr = true + linsolve = :lsqr + end # Woodbury uses a grid; dense and LSQR use an edge list. supp = if use_woodbury C = fill(T(-Inf), m, n) diff --git a/src/sparse_support.jl b/src/sparse_support.jl index 6836ebb..4121d1d 100644 --- a/src/sparse_support.jl +++ b/src/sparse_support.jl @@ -77,8 +77,7 @@ end # Native minimal-cover (MMC) solvers # ============================================================ -# Sparse `AbsLog{2}` solvers default to matrix-free LSQR. Use `linsolve=:dense` -# or `:auto` to request factorization-based paths. +# Sparse `AbsLog{2}` solvers default to LSQR; `:auto` uses support density. function symcover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) a, _ = _symcover_min_abslog2(A; linsolve, kwargs...) return a diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index a95b131..31af75d 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -295,6 +295,31 @@ end @test MatrixCovers._cover_min_abslog2(Gbig)[3].linsolve === :dense end +@testset "MMC :auto measures the stored support" begin + T3 = SymTridiagonal(fill(4.0, 30), fill(1.0, 29)) + aT, sT = MatrixCovers._symcover_min_abslog2(T3) + @test sT.linsolve === :lsqr + @test aT ≈ MatrixCovers._symcover_min_abslog2(Matrix(T3); linsolve=:dense)[1] rtol=1e-6 + @test iscover(aT, T3; rtol=1e-8) + + # Wrapped sparse support selects LSQR. + rng = StableRNG(3) + S0 = sprand(rng, 40, 40, 0.05) + Sw = Symmetric(S0 + S0' + I) + aw, bw, stw = MatrixCovers._cover_min_abslog2(Sw) + @test stw.linsolve === :lsqr + @test iscover(aw, bw, Sw; rtol=1e-8) + + X8 = exp.(randn(rng, 8, 8)) + @test MatrixCovers._symcover_min_abslog2((X8 .+ X8') ./ 2)[2].linsolve === :woodbury + + # The threshold includes equality. + D4 = Matrix(Diagonal([4.0, 9.0, 1.0, 16.0])) + @test MatrixCovers._symcover_min_abslog2(D4)[2].linsolve === :lsqr + D4[1, 2] = D4[2, 1] = 1.0 + @test MatrixCovers._symcover_min_abslog2(D4)[2].linsolve === :dense +end + # Exact inner solves stop a stage when the violated set is unchanged; LSQR uses # the decrease test. @testset "MMC exact paths stop on a sign-stable Newton step" begin @@ -353,7 +378,7 @@ end singletons(vals) = Matrix(sparse(1:length(vals), 1:length(vals), float.(vals))) # k singleton components block2(k) = cat(([2.0+i i; i 3.0+i] for i in 1:k)...; dims = (1, 2)) # k dense 2×2 components for M in (singletons([4.0, 9.0, 1.0]), singletons(1.0:6.0), block2(3), block2(6)) - ad, bd = cover_min(AbsLog{2}(), M) # :auto = dense + ridge + ad, bd = cover_min(AbsLog{2}(), M; linsolve = :dense) # ridge lifts the unpinned gauges al, bl = cover_min(AbsLog{2}(), M; linsolve = :lsqr) @test iscover(ad, bd, M; atol=1e-8) @test iscover(al, bl, M; atol=1e-8)