Improve asymptotics

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...)
 

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)

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

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)

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
@@ -117,3 +116,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

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)

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])

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)

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]