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/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/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/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 e863ef4..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 @@ -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 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) 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]) 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) 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)