From a6fa3b990a2467c5f2467861e44a98ffeecd7359 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 23 Aug 2026 19:08:38 -0500 Subject: [PATCH 01/15] Evaluate the heuristic covers' logarithms once Cache dense log magnitudes and sweep the grid for heuristic covers. Preserve traversal order while storing boost entries contiguously and deriving balance counts from the grid. Assisted-by: Claude Opus (claude-opus-5) --- src/MatrixCovers.jl | 1 + src/dense_heuristic.jl | 314 ++++++++++++++++++++++++++++++++++++++++ src/heuristic_covers.jl | 129 ++++++++++------- src/support.jl | 48 +++++- 4 files changed, 437 insertions(+), 55 deletions(-) create mode 100644 src/dense_heuristic.jl diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index 056a566..50de63a 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -25,6 +25,7 @@ include("penalties.jl") include("support.jl") include("iscover.jl") include("heuristic_covers.jl") +include("dense_heuristic.jl") # full-grid kernels for the heuristic covers include("gram_covers.jl") # symmetric covers of A'*W*A from an asymmetric cover of A include("initializers.jl") # the start menu; consumed by both solver families below include("soft_covers.jl") diff --git a/src/dense_heuristic.jl b/src/dense_heuristic.jl new file mode 100644 index 0000000..18964c6 --- /dev/null +++ b/src/dense_heuristic.jl @@ -0,0 +1,314 @@ +# Dense-grid kernels for the heuristic covers. +# +# `symcover!` and `cover!` sweep the support five or six times over, and every +# sweep of the callback-driven implementations recomputes `log(abs(A[i,j]))`. +# When the support is the full grid those logarithms dominate the run time, so +# these kernels evaluate them once into a log-magnitude grid and run the +# remaining sweeps over that grid. Every sum accumulates in the traversal order +# of `foreach_support`/`foreach_support_sym`, and each residual is formed from +# the same three numbers, so the covers agree with the callback-driven path bit +# for bit. + +# Use the grid only when the support traversal is already dense and its +# allocation is worthwhile. +const DENSE_GRID_MIN = 64 + +_dense_grid_storage(::AbstractMatrix) = false +_dense_grid_storage(::StridedMatrix) = true +_dense_grid_storage(A::Union{Symmetric,Hermitian}) = _dense_grid_storage(parent(A)) + +const DENSE_GRID_FLOAT = Union{Float32,Float64} + +_use_dense_grid(A::AbstractMatrix, ::Type{T}) where {T} = + T <: DENSE_GRID_FLOAT && _dense_grid_storage(A) && minimum(size(A)) >= DENSE_GRID_MIN + +# Slot of `L[1, j]` minus one in a column-packed upper triangle: column `j` +# occupies `_trioff(j)+1 : _trioff(j)+j`. +_trioff(j::Int) = (j * (j - 1)) >> 1 + +# `Lp[_trioff(j)+i] = log(abs(A[i,j]))` for `i <= j`, `-Inf` where the entry is +# zero, alongside the per-row log sums and support counts of +# `unconstrained_min!`. `s[r]` receives the partners of row `r` in increasing +# partner order — the columns `1:r-1` of row `r`'s own column first, then the +# diagonal, then the later columns — which is the order `foreach_support_sym` +# feeds them in. +function _tri_logabs!(Lp::Vector{T}, s::Vector{T}, cnt::Vector{Int}, A::AbstractMatrix) where {T} + ax = axes(A, 1) + or = first(ax) - 1 + n = length(ax) + fill!(s, zero(T)) + fill!(cnt, 0) + for jp in 1:n + j = jp + or + o = _trioff(jp) + sj = zero(T) + cj = 0 + for ip in 1:jp-1 + v = abs(A[ip+or, j]) + if iszero(v) + Lp[o+ip] = T(-Inf) + else + l = log(T(v)) + Lp[o+ip] = l + s[ip] += l + cnt[ip] += 1 + sj += l + cj += 1 + end + end + s[jp] += sj + cnt[jp] += cj + v = abs(A[j, j]) + if iszero(v) + Lp[o+jp] = T(-Inf) + else + l = log(T(v)) + Lp[o+jp] = l + s[jp] += l + cnt[jp] += 1 + end + end + return Lp +end + +# `L[i,j] = log(abs(A[i,j]))`, `-Inf` where the entry is zero, alongside the +# per-row and per-column log sums and support counts of `unconstrained_min!`. +function _grid_logabs!(L::Matrix{T}, sa::Vector{T}, sb::Vector{T}, + na::Vector{Int}, nb::Vector{Int}, A::AbstractMatrix) where {T} + or = first(axes(A, 1)) - 1 + oc = first(axes(A, 2)) - 1 + m, n = size(A) + fill!(sa, zero(T)) + fill!(na, 0) + for jp in 1:n + j = jp + oc + sj = zero(T) + cj = 0 + for ip in 1:m + v = abs(A[ip+or, j]) + if iszero(v) + L[ip, jp] = T(-Inf) + else + l = log(T(v)) + L[ip, jp] = l + sa[ip] += l + na[ip] += 1 + sj += l + cj += 1 + end + end + sb[jp] = sj + nb[jp] = cj + end + return L +end + +# Keep supported scales positive when `exp` underflows. +_uncon_scale(si::T, ni::Int, halfmu::T) where {T} = + iszero(ni) ? zero(T) : max(exp(si / ni - halfmu), floatmin(T)) + +# Greedy boost that updates scales and log scales together. +function _dense_boost!(α::Vector{T}, lα::Vector{T}, entries, zmax::T) where {T} + deficit((i, j, lv)) = lv - lα[i] - lα[j] + function apply!((i, j, lv), z) + h = z / 2 + lα[i] += h; α[i] = exp(lα[i]) + i == j || (lα[j] += h; α[j] = exp(lα[j])) + end + bucket_boost!(deficit, apply!, entries, T, zmax) + return α +end + +# `symcover!` over a packed upper-triangular log-magnitude grid. +function _symcover_dense!(a::AbstractVector, A::AbstractMatrix, ::Type{T}, maxiter::Int) where {T} + ax = axes(A, 1) + or = first(ax) - 1 + n = length(ax) + Lp = Vector{T}(undef, _trioff(n) + n) + α = Vector{T}(undef, n) + lα = Vector{T}(undef, n) + cnt = Vector{Int}(undef, n) + _tri_logabs!(Lp, α, cnt, A) # `α` carries the row log sums here + nztotal = sum(cnt) + halfmu = iszero(nztotal) ? zero(T) : sum(α) / (2 * nztotal) + for ip in 1:n + α[ip] = _uncon_scale(α[ip], cnt[ip], halfmu) + lα[ip] = log(α[ip]) + end + + # Only initially violated entries can require a boost. + nviol = 0 + zmax = zero(T) + for jp in 1:n + o = _trioff(jp) + lj = lα[jp] + for ip in 1:jp + z = Lp[o+ip] - lα[ip] - lj + nviol += ifelse(z > zero(T), 1, 0) + zmax = ifelse(z > zmax, z, zmax) + end + end + # A supported zero scale produces an infinite deficit. + isfinite(zmax) || + throw(ArgumentError("boost_feasible! requires a start with positive scale on every supported row")) + entries = Vector{Tuple{Int,Int,T}}(undef, nviol + 1) + k = 1 + for jp in 1:n + o = _trioff(jp) + lj = lα[jp] + for ip in 1:jp + lv = Lp[o+ip] + entries[k] = (ip, jp, lv) + k += ifelse(lv - lα[ip] - lj > zero(T), 1, 0) + end + end + resize!(entries, nviol) + _dense_boost!(α, lα, entries, zmax) + + lratio = Vector{T}(undef, n) + for _ in 1:maxiter + map!(log, lα, α) + fill!(lratio, T(Inf)) + # `-Inf` slots produce `Inf` or `NaN` ratios, both ignored below. + for jp in 1:n + o = _trioff(jp) + lj = lα[jp] + mj = T(Inf) + for ip in 1:jp-1 + lr = lα[ip] + lj - Lp[o+ip] + lratio[ip] = ifelse(lr < lratio[ip], lr, lratio[ip]) + mj = ifelse(lr < mj, lr, mj) + end + lr = lj + lj - Lp[o+jp] + mj = ifelse(lr < mj, lr, mj) + lratio[jp] = ifelse(mj < lratio[jp], mj, lratio[jp]) + end + for ip in 1:n + lr = lratio[ip] + # Infinite ratios require no update. + isinf(lr) || (α[ip] = _tighten_shrink(α[ip], lr)) + end + end + for ip in 1:n + a[ip+or] = α[ip] + end + return a +end + +# `cover!` over a dense log-magnitude grid, up to but not including the balance +# convention. +function _cover_dense!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, + ::Type{T}, maxiter::Int) where {T} + or = first(axes(A, 1)) - 1 + oc = first(axes(A, 2)) - 1 + m, n = size(A) + L = Matrix{T}(undef, m, n) + α = Vector{T}(undef, m) + β = Vector{T}(undef, n) + lα = Vector{T}(undef, m) + lβ = Vector{T}(undef, n) + na = Vector{Int}(undef, m) + nb = Vector{Int}(undef, n) + _grid_logabs!(L, α, β, na, nb, A) # `α`, `β` carry the log sums here + nztotal = sum(na) + halfmu = iszero(nztotal) ? zero(T) : sum(α) / (2 * nztotal) + for ip in 1:m + α[ip] = _uncon_scale(α[ip], na[ip], halfmu) + lα[ip] = log(α[ip]) + end + for jp in 1:n + β[jp] = _uncon_scale(β[jp], nb[jp], halfmu) + lβ[jp] = log(β[jp]) + end + + # Branchless selection, as in the symmetric kernel. + nviol = 0 + zmax = zero(T) + for jp in 1:n + lj = lβ[jp] + for ip in 1:m + z = L[ip, jp] - lα[ip] - lj + nviol += ifelse(z > zero(T), 1, 0) + zmax = ifelse(z > zmax, z, zmax) + end + end + isfinite(zmax) || + throw(ArgumentError("boost_feasible! requires a start with positive scale on every supported row/column")) + entries = Vector{Tuple{Int,Int,T}}(undef, nviol + 1) + k = 1 + for jp in 1:n + lj = lβ[jp] + for ip in 1:m + lv = L[ip, jp] + entries[k] = (ip, jp, lv) + k += ifelse(lv - lα[ip] - lj > zero(T), 1, 0) + end + end + resize!(entries, nviol) + # Row and column scales require separate updates. + deficit((i, j, lv)) = lv - lα[i] - lβ[j] + function apply!((i, j, lv), z) + h = z / 2 + lα[i] += h; α[i] = exp(lα[i]) + lβ[j] += h; β[j] = exp(lβ[j]) + end + bucket_boost!(deficit, apply!, entries, T, zmax) + + ratioa = Vector{T}(undef, m) + for _ in 1:maxiter + map!(log, lα, α) + map!(log, lβ, β) + fill!(ratioa, T(Inf)) + for jp in 1:n + lj = lβ[jp] + mj = T(Inf) + for ip in 1:m + lr = lα[ip] + lj - L[ip, jp] + ratioa[ip] = ifelse(lr < ratioa[ip], lr, ratioa[ip]) + mj = ifelse(lr < mj, lr, mj) + end + isinf(mj) || (β[jp] = _tighten_shrink(β[jp], mj)) + end + for ip in 1:m + lr = ratioa[ip] + isinf(lr) || (α[ip] = _tighten_shrink(α[ip], lr)) + end + end + + for ip in 1:m + a[ip+or] = α[ip] + end + for jp in 1:n + b[jp+oc] = β[jp] + end + # Balance the factors, then restore coverage lost to rounding. + _balance_cover!(a, b, A) + for ip in 1:m + lα[ip] = log(T(a[ip+or])) + end + for jp in 1:n + lβ[jp] = log(T(b[jp+oc])) + end + t = zero(T) + for jp in 1:n + lj = lβ[jp] + tj = zero(T) + # Off-support entries produce `-Inf` or `NaN`, both ignored by `>`. + for ip in 1:m + u = (L[ip, jp] - lα[ip] - lj) / 2 + tj = ifelse(u > tj, u, tj) + end + t = ifelse(tj > t, tj, t) + end + # A supported row or column with zero scale gives `t = +Inf`. + isfinite(t) || + throw(ArgumentError("inflate_feasible! requires a start with positive scale on every supported row/column")) + iszero(t) && return a, b + for ip in 1:m + iszero(a[ip+or]) || (a[ip+or] = exp(lα[ip] + t)) + end + for jp in 1:n + iszero(b[jp+oc]) || (b[jp+oc] = exp(lβ[jp] + t)) + end + return a, b +end diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index 24d4ee4..d0dd9a4 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -64,9 +64,15 @@ function symcover!(a::AbstractVector, A::AbstractMatrix; kwargs...) axes(A, 2) == ax || throw(ArgumentError("symcover! requires a square matrix")) require_abs_symmetric(A, :symcover!) eachindex(a) == ax || throw(DimensionMismatch("indices of `a` must match the indexing of `A`, got eachindex(a)=$(string(eachindex(a))), axes(A, 1)=$(string(ax))")) + return _symcover!(a, A; kwargs...) +end + +function _symcover!(a::AbstractVector, A::AbstractMatrix; maxiter::Int=3) + T = float(real(eltype(a))) + _use_dense_grid(A, T) && return _symcover_dense!(a, A, T, maxiter) unconstrained_min!(AbsLog{2}(), a, A) boost_feasible!(a, A) - return tighten_cover!(a, A; kwargs...) + return tighten_cover!(a, A; maxiter) end """ @@ -137,9 +143,15 @@ cover!(ϕ::AbstractCoverPenalty, a::AbstractVector, b::AbstractVector, A::Abstra function cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) axes(A, 1) == eachindex(a) || throw(DimensionMismatch("indices of `a` must match row-indexing of `A`, got eachindex(a)=$(string(eachindex(a))), axes(A, 1)=$(string(axes(A, 1)))")) axes(A, 2) == eachindex(b) || throw(DimensionMismatch("indices of `b` must match column-indexing of `A`, got eachindex(b)=$(string(eachindex(b))), axes(A, 2)=$(string(axes(A, 2)))")) + return _cover!(a, b, A; kwargs...) +end + +function _cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=3) + T = float(promote_type(eltype(a), eltype(b))) + _use_dense_grid(A, T) && return _cover_dense!(a, b, A, T, maxiter) unconstrained_min!(AbsLog{2}(), a, b, A) boost_feasible!(a, b, A) - tighten_cover!(a, b, A; kwargs...) + tighten_cover!(a, b, A; maxiter) # Apply the package's balance convention, then restore coverage lost to rounding. _balance_cover!(a, b, A) return inflate_feasible!(a, b, A) @@ -163,28 +175,32 @@ end # exactly, at the cost of balancing only within a factor of `sqrt(2)`. function _balance_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) T = float(promote_type(eltype(a), eltype(b))) - rowcomp, colcomp, ncomp = _support_components(A) + rowcomp, colcomp, ncomp, nzrow, nzcol = _support_components(A) iszero(ncomp) && return a, b Lα = zeros(T, ncomp) Lβ = zeros(T, ncomp) nnz = zeros(Int, ncomp) - or = first(axes(A, 1)) - 1 - oc = first(axes(A, 2)) - 1 - foreach_support(A) do i, j, v - c = rowcomp[i-or] - Lα[c] += log2(T(a[i])) - Lβ[c] += log2(T(b[j])) - nnz[c] += 1 + # Weight each scale by its support count. + for (p, i) in enumerate(eachindex(a)) + c = rowcomp[p] + c == 0 && continue + Lα[c] += nzrow[p] * log2(T(a[i])) + nnz[c] += nzrow[p] + end + for (q, j) in enumerate(eachindex(b)) + c = colcomp[q] + c == 0 && continue + Lβ[c] += nzcol[q] * log2(T(b[j])) end # An integer base-2 exponent makes the rescaling exact. gamma = [exp2(round((Lβ[c] - Lα[c]) / (2 * nnz[c]))) for c in 1:ncomp] - for i in eachindex(a) - c = rowcomp[i-or] + for (p, i) in enumerate(eachindex(a)) + c = rowcomp[p] c == 0 && continue a[i] *= gamma[c] end - for j in eachindex(b) - c = colcomp[j-oc] + for (q, j) in enumerate(eachindex(b)) + c = colcomp[q] c == 0 && continue b[j] /= gamma[c] end @@ -361,45 +377,52 @@ end # to lower buckets. Log-deficit buckets preserve covariance except for ties. const BOOST_BUCKET_WIDTH = log(2) / 4 # quality indistinguishable from exact greedy; only bucket count grows as w shrinks -# Flat linked bucket queue: `head[b]` is the first entry and `nxt[k]` links the -# rest. Deficits are recomputed instead of cached. -function bucket_boost!(deficit::F, apply!::G, entries, ::Type{T}) where {F,G,T} - n = length(entries) - zmax = zero(T) - for entry in entries - z = deficit(entry) - z > zero(T) || continue - zmax = max(zmax, z) - end +# Visit deficit buckets from highest to lowest. Original entries are stored +# contiguously; entries demoted from higher buckets use per-bucket stacks. +function bucket_boost!(deficit::F, apply!::G, entries::AbstractVector, ::Type{T}, zmax::T) where {F,G,T} zmax > zero(T) || return w = T(BOOST_BUCKET_WIDTH) B = max(1, ceil(Int, zmax / w)) bucketof(z) = clamp(ceil(Int, z / w), 1, B) - head = zeros(Int, B) - nxt = zeros(Int, n) - for k in eachindex(entries) - z = deficit(entries[k]) + ptr = zeros(Int, B + 1) + for entry in entries + z = deficit(entry) + z > zero(T) || continue + ptr[bucketof(z)+1] += 1 + end + ptr[1] = 1 + cumsum!(ptr, ptr) + cursor = ptr[1:end-1] # next free slot of each level + sorted = similar(entries, ptr[end] - 1) + for k in reverse(eachindex(entries)) + entry = entries[k] + z = deficit(entry) z > zero(T) || continue b = bucketof(z) - nxt[k] = head[b] - head[b] = k + sorted[cursor[b]] = entry + cursor[b] += 1 + end + # Demoted entries, as a stack per level over one shared array. + dhead = zeros(Int, B) + dentry = similar(entries, 0) + dnext = Int[] + demote!(entry, b2) = (push!(dentry, entry); push!(dnext, dhead[b2]); dhead[b2] = length(dentry)) + function visit!(entry, b) + z = deficit(entry) + z > zero(T) || return + b2 = bucketof(z) + b2 < b ? demote!(entry, b2) : apply!(entry, z) + return end for b in B:-1:1 - k = head[b] - while k != 0 - knext = nxt[k] # save before a possible demotion overwrites nxt[k] - e = entries[k] - z = deficit(e) - if z > zero(T) - b2 = bucketof(z) - if b2 < b - nxt[k] = head[b2] - head[b2] = k # deficit shrank: demote, revisit later - else - apply!(e, z) - end - end - k = knext + e = dhead[b] + while e != 0 + enext = dnext[e] # save before a further demotion appends + visit!(dentry[e], b) + e = enext + end + for s in ptr[b]:ptr[b+1]-1 + visit!(sorted[s], b) end end return @@ -427,9 +450,13 @@ function boost_feasible!(a::AbstractVector{T}, A::AbstractMatrix) where T # `entries` once at its exact size, instead of the repeated grow-and-copy # of building it with `push!`. nviol = Ref(0) + zmax = Ref(zero(T)) foreach_support_sym(A) do i, j, v z = log(T(v)) - la[i] - la[j] - z > zero(T) && (nviol[] += 1) + if z > zero(T) + nviol[] += 1 + zmax[] = max(zmax[], z) + end end entries = Vector{Tuple{IdxT,IdxT,T}}(undef, nviol[]) k = Ref(0) @@ -449,7 +476,7 @@ function boost_feasible!(a::AbstractVector{T}, A::AbstractMatrix) where T la[i] += h; a[i] = exp(la[i]) i == j || (la[j] += h; a[j] = exp(la[j])) end - bucket_boost!(deficit, apply!, entries, T) + bucket_boost!(deficit, apply!, entries, T, zmax[]) return a end @@ -470,9 +497,13 @@ function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix # fill) allocate `entries` once at its exact size, instead of the # repeated grow-and-copy of building it with `push!`. nviol = Ref(0) + zmax = Ref(zero(T)) foreach_support(A) do i, j, v z = log(T(v)) - la[i] - lb[j] - z > zero(T) && (nviol[] += 1) + if z > zero(T) + nviol[] += 1 + zmax[] = max(zmax[], z) + end end entries = Vector{Tuple{IdxA,IdxB,T}}(undef, nviol[]) k = Ref(0) @@ -492,7 +523,7 @@ function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix la[i] += h; a[i] = exp(la[i]) lb[j] += h; b[j] = exp(lb[j]) end - bucket_boost!(deficit, apply!, entries, T) + bucket_boost!(deficit, apply!, entries, T, zmax[]) return a, b end diff --git a/src/support.jl b/src/support.jl index 2e6d796..e9c47a7 100644 --- a/src/support.jl +++ b/src/support.jl @@ -97,6 +97,38 @@ function require_abs_symmetric(A::AbstractMatrix, fname) return nothing end +# Cache-blocked check for dense storage: the transposed reads of a column-major +# sweep miss on every entry once `A` outgrows the cache, while a block of rows and +# its transpose both fit. Only `i < j` needs testing, the diagonal being its own +# partner. +const SYMMETRY_BLOCK = 64 + +function require_abs_symmetric(A::StridedMatrix, fname) + ax = axes(A, 1) + axes(A, 2) == ax || + throw(DimensionMismatch("$fname requires a square matrix, got axes $(string(axes(A)))")) + n = length(ax) + o = first(ax) - 1 + for jb in 1:SYMMETRY_BLOCK:n + jlast = min(jb + SYMMETRY_BLOCK - 1, n) + for ib in 1:SYMMETRY_BLOCK:jlast + for jp in jb:jlast + for ip in ib:min(ib + SYMMETRY_BLOCK - 1, jp - 1) + i, j = ip + o, jp + o + v = abs(A[i, j]) + w = abs(A[j, i]) + m = max(v, w) + abs(v - w) <= ASYMMETRY_ULPS * eps(float(real(typeof(m)))) * m || throw(ArgumentError(""" + $fname requires `abs.(A)` to be symmetric, but abs(A[$(string(i)),$(string(j))]) = $(string(v)) and \ + abs(A[$(string(j)),$(string(i))]) = $(string(w)). Wrap `A` in `Symmetric` (or `Hermitian`) to name the \ + triangle to read; that also skips this check.""")) + end + end + end + end + return nothing +end + # Storage that makes the precondition structural: the wrapper or the type's own # invariant already guarantees `abs(A[i,j]) == abs(A[j,i])`. require_abs_symmetric(::Union{Symmetric,Hermitian,Diagonal,SymTridiagonal}, fname) = nothing @@ -106,7 +138,9 @@ require_abs_symmetric(::Union{Symmetric,Hermitian,Diagonal,SymTridiagonal}, fnam # ncomp)`, where `rowcomp` and `colcomp` are `Vector{Int}` indexed by *position* # within `axes(A, 1)` and `axes(A, 2)` (so offset axes need no special case, as # with `GroupedSupport.ptr`), holding the component id in `1:ncomp` — or 0 for -# rows/columns with empty support, which belong to no component. +# rows/columns with empty support, which belong to no component. `nzrow` and +# `nzcol` count the stored entries of each row and column, in the same position +# space. # # The gauge orbit of an asymmetric cover has one dimension per component: the # rescaling `a -> γ*a`, `b -> b/γ` acts independently on each, because no @@ -118,7 +152,6 @@ function _support_components(A::AbstractMatrix) m = length(axes(A, 1)) n = length(axes(A, 2)) parent = collect(1:(m + n)) - touched = falses(m + n) function find(p) while parent[p] != p parent[p] = parent[parent[p]] # path halving @@ -128,10 +161,13 @@ function _support_components(A::AbstractMatrix) end or = first(axes(A, 1)) - 1 oc = first(axes(A, 2)) - 1 + nzrow = zeros(Int, m) + nzcol = zeros(Int, n) foreach_support(A) do i, j, _ p = i - or q = m + j - oc - touched[p] = touched[q] = true + nzrow[p] += 1 + nzcol[q-m] += 1 rp, rq = find(p), find(q) rp == rq || (parent[rp] = rq) end @@ -140,7 +176,7 @@ function _support_components(A::AbstractMatrix) rowcomp = zeros(Int, m) colcomp = zeros(Int, n) for p in 1:(m + n) - touched[p] || continue + (p <= m ? nzrow[p] : nzcol[p-m]) > 0 || continue r = find(p) if label[r] == 0 ncomp += 1 @@ -152,7 +188,7 @@ function _support_components(A::AbstractMatrix) colcomp[p-m] = label[r] end end - return rowcomp, colcomp, ncomp + return rowcomp, colcomp, ncomp, nzrow, nzcol end """ @@ -183,7 +219,7 @@ read through [`foreach_support`](@ref). See also: [`SupportComponents`](@ref). """ function support_components(A::AbstractMatrix) - rowcomp, colcomp, ncomp = _support_components(A) + rowcomp, colcomp, ncomp, _, _ = _support_components(A) return SupportComponents(rowcomp, colcomp, ncomp, axes(A, 1), axes(A, 2)) end From 43495fe542e23726e85b63aa0b01995e956b1b07 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 23 Aug 2026 19:26:42 -0500 Subject: [PATCH 02/15] Read the native solvers' layout from the matrix Build Woodbury grids directly from `A` and construct edge-list supports only for the other solver paths. Count support first to select the solver. Assisted-by: Claude Opus (claude-opus-5) --- src/minimal_covers.jl | 72 +++++++++++++++++++------------------------ 1 file changed, 31 insertions(+), 41 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 68748da..1bc1f36 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -1038,16 +1038,20 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end n = length(ax) use_lsqr = linsolve === :lsqr - # Build only the support layout needed by the chosen solver. - G = _sym_support(A, T) + # One counting traversal decides the solver; the layout it needs is built after. + o = first(ax) - 1 + nza = zeros(Int, n) # support entries per row, counted in both orientations + foreach_support_sym(A) do i, j, v + nza[i-o] += 1 + i == j || (nza[j-o] += 1) + end hassupp = falses(n) nsupp = 0 # support entries of `A`, counted in both orientations maxzero = 0 # largest number of zeros in any row of `A` - for (ip, i) in enumerate(ax) - ns = length(_slots(G, i)) - hassupp[ip] = ns > 0 - nsupp += ns - maxzero = max(maxzero, n - ns) + for ip in 1:n + hassupp[ip] = nza[ip] > 0 + nsupp += nza[ip] + maxzero = max(maxzero, n - nza[ip]) end # `n*I - L_Z` is positive definite only while no row carries more than # `n ÷ 4` zeros; the total budget bounds cost. @@ -1072,14 +1076,12 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # Woodbury uses a grid; dense and LSQR use an edge list. supp = if use_woodbury C = fill(T(-Inf), n, n) - for (ip, i) in enumerate(ax) - for s in _slots(G, i) - jp = G.idx[s] - first(ax) + 1 - jp >= ip && (C[ip, jp] = log(G.val[s])) - end + foreach_support_sym(A) do i, j, v + C[i-o, j-o] = log(T(v)) end Grid{T}(C) else + G = _sym_support(A, T) edges = Tuple{Int,Int}[] cvals = T[] for (ip, i) in enumerate(ax) @@ -1093,18 +1095,14 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end EdgeList{T}(edges, cvals) end - # Zero set defining the off-diagonal pattern of sparse `C`. + # Zero set defining the off-diagonal pattern of sparse `C`, grouped by row. zedges = Tuple{Int,Int}[] if use_woodbury - mark = falses(n) - for (ip, i) in enumerate(ax) - for s in _slots(G, i) - mark[G.idx[s] - first(ax) + 1] = true - end + Cgrid = supp.C + for ip in 1:n for jp in ip:n - mark[jp] || push!(zedges, (ip, jp)) + isfinite(Cgrid[ip, jp]) || push!(zedges, (ip, jp)) end - fill!(mark, false) end end # The ridge handles singular symmetric support graphs. @@ -1144,18 +1142,15 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), N = m + n use_lsqr = linsolve === :lsqr # Stack row positions before column positions; scatter results back to `A`'s axes. - G = _row_support(A, T) + or = first(axr) - 1 + oc = first(axc) - 1 nzrow = zeros(Int, m) # support entries per row, for the balance convention nzcol = zeros(Int, n) # ditto per column - ne = 0 - for (ip, i) in enumerate(axr) - for s in _slots(G, i) - jp = G.idx[s] - first(axc) + 1 - ne += 1 - nzrow[ip] += 1 - nzcol[jp] += 1 - end + foreach_support(A) do i, j, v + nzrow[i-or] += 1 + nzcol[j-oc] += 1 end + ne = sum(nzrow) hasrow = nzrow .> 0 hascol = nzcol .> 0 # `min(m,n)*I - L_Z` is positive definite only while no row or column @@ -1190,13 +1185,12 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # Woodbury uses a grid; dense and LSQR use an edge list. supp = if use_woodbury C = fill(T(-Inf), m, n) - for (ip, i) in enumerate(axr) - for s in _slots(G, i) - C[ip, G.idx[s]-first(axc)+1] = log(G.val[s]) - end + foreach_support(A) do i, j, v + C[i-or, j-oc] = log(T(v)) end Grid{T}(C) else + G = _row_support(A, T) edges = Tuple{Int,Int}[] cvals = T[] for (ip, i) in enumerate(axr) @@ -1220,15 +1214,11 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), dfull[m+jp] = T(m) Umat[m+jp, 2] = oneunit(T) end - mark = falses(n) - for (ip, i) in enumerate(axr) - for s in _slots(G, i) - mark[G.idx[s] - first(axc) + 1] = true - end + Cgrid = supp.C + for ip in 1:m for jp in 1:n - mark[jp] || push!(zedges, (ip, m + jp)) + isfinite(Cgrid[ip, jp]) || push!(zedges, (ip, m + jp)) end - fill!(mark, false) end end # Pin the global row/column gauge on supported variables. @@ -1255,7 +1245,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end x, stats = _abslog2_continuation(sys, x0; κs, maxiter, linsolve, boost) # Apply the balance convention independently to each support component. - rowcomp, colcomp, ncomp = _support_components(A) + rowcomp, colcomp, ncomp, _, _ = _support_components(A) Lα = zeros(T, ncomp) Lβ = zeros(T, ncomp) nec = zeros(Int, ncomp) From a72038a18fc9f49b20dc4705902080e68bf0239f Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 23 Aug 2026 19:26:42 -0500 Subject: [PATCH 03/15] Rebuild the Woodbury correction in place Assemble the sparse Woodbury correction in reusable CSC storage and multiply its symmetric upper triangle directly. Separate Boolean pattern scans from floating-point loops. Assisted-by: Claude Opus (claude-opus-5) --- src/minimal_covers.jl | 189 ++++++++++++++++++++++++++++-------------- 1 file changed, 129 insertions(+), 60 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 1bc1f36..f619ef5 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -492,21 +492,28 @@ function _fκpat(x, κ, pat, supp::Grid{T}, symmetric::Bool) where {T} pj = view(pat, 1:j-1, j) xi = view(x, 1:j-1) vj = zero(T) - dj = 0 - @simd for i in eachindex(cj, pj, xi) + # Keep the Boolean pattern out of the vectorized floating-point loop. + @simd for i in eachindex(cj, xi) c = cj[i] - fin = isfinite(c) z = xi[i] + xj - c w = ifelse(z < 0, κT, oneunit(T)) - vj += ifelse(fin, w * z^2, zero(T)) - dj += ifelse((fin & (z < 0)) == pj[i], 0, 1) + vj += ifelse(isfinite(c), w * z^2, zero(T)) + end + # Stop comparing after the first pattern change. + if ndiff == 0 + dj = 0 + @simd for i in eachindex(cj, pj, xi) + c = cj[i] + dj += ifelse((isfinite(c) & (xi[i] + xj - c < 0)) == pj[i], 0, 1) + end + ndiff += dj end c = C[j, j] fin = isfinite(c) z = 2xj - c w = ifelse(z < 0, κT, oneunit(T)) v += 2vj + ifelse(fin, w * z^2, zero(T)) - ndiff += dj + ifelse((fin & (z < 0)) == pat[j, j], 0, 1) + ndiff += ifelse((fin & (z < 0)) == pat[j, j], 0, 1) end else xr = view(x, 1:m) @@ -515,17 +522,21 @@ function _fκpat(x, κ, pat, supp::Grid{T}, symmetric::Bool) where {T} cj = view(C, :, j) pj = view(pat, :, j) vj = zero(T) - dj = 0 - @simd for i in eachindex(cj, pj, xr) + @simd for i in eachindex(cj, xr) c = cj[i] - fin = isfinite(c) z = xr[i] + xj - c w = ifelse(z < 0, κT, oneunit(T)) - vj += ifelse(fin, w * z^2, zero(T)) - dj += ifelse((fin & (z < 0)) == pj[i], 0, 1) + vj += ifelse(isfinite(c), w * z^2, zero(T)) + end + if ndiff == 0 + dj = 0 + @simd for i in eachindex(cj, pj, xr) + c = cj[i] + dj += ifelse((isfinite(c) & (xr[i] + xj - c < 0)) == pj[i], 0, 1) + end + ndiff += dj end v += vj - ndiff += dj end end return v, ndiff == 0 @@ -579,7 +590,9 @@ end # `κ === nothing` denotes the unweighted solve. Off-support entries of `C` are # `-Inf`, so products with them go through `ifelse(isfinite(c), ...)`: `0 * -Inf` # is NaN. -function _assemble_woodbury!(f, dg, degV, vedges, vpat, x, κ, supp::Grid{T}, +# `vrow` stores violated off-diagonal rows in compressed-column order; `vcnt` +# stores their column counts. `dg` and `degV` carry diagonal entries. +function _assemble_woodbury!(f, dg, degV, vrow, vcnt, vpat, x, κ, supp::Grid{T}, symmetric::Bool, dκ) where {T} C = supp.C m, n = size(C) @@ -609,18 +622,20 @@ function _assemble_woodbury!(f, dg, degV, vedges, vpat, x, κ, supp::Grid{T}, w = ifelse(viol, κT, oneunit(T)) f[j] += fq + ifelse(fin, w * c, zero(T)) vpat[j, j] = viol + nv = 0 for i in eachindex(vj) vj[i] || continue - push!(vedges, (i, j)) + push!(vrow, i) + nv += 1 degV[i] += 1 degV[j] += 1 dg[i] += dκ dg[j] += dκ end + vcnt[j] = nv # A symmetric diagonal entry sits at both ends of its own residual, so it # lands on `dg[j]` twice while counting once in `degV`. if vpat[j, j] - push!(vedges, (j, j)) degV[j] += 1 dg[j] += 2dκ end @@ -645,19 +660,82 @@ function _assemble_woodbury!(f, dg, degV, vedges, vpat, x, κ, supp::Grid{T}, vj[i] = viol end f[q] += fq + nv = 0 for i in eachindex(vj) vj[i] || continue - push!(vedges, (i, q)) + push!(vrow, i) + nv += 1 degV[i] += 1 degV[q] += 1 dg[i] += dκ dg[q] += dκ end + vcnt[q] = nv end end return f end +# Assemble the upper triangle of the sparse Woodbury correction in CSC order, +# reusing its storage across solves. +function _assemble_C!(colptr::Vector{Int}, rowval::Vector{Int}, nzval::Vector{T}, + zptr, zrow, vptr, vrow, dg, dκ, N::Int) where {T} + colptr[1] = 1 + for q in 1:N + colptr[q+1] = colptr[q] + (zptr[q+1] - zptr[q]) + (vptr[q+1] - vptr[q]) + 1 + end + nnz = colptr[N+1] - 1 + length(rowval) == nnz || resize!(rowval, nnz) + length(nzval) == nnz || resize!(nzval, nnz) + for q in 1:N + s = colptr[q] + zs, ze = zptr[q], zptr[q+1] - 1 + vs, ve = vptr[q], vptr[q+1] - 1 + while zs <= ze && vs <= ve + if zrow[zs] < vrow[vs] + rowval[s] = zrow[zs]; nzval[s] = -oneunit(T); zs += 1 + else + rowval[s] = vrow[vs]; nzval[s] = dκ; vs += 1 + end + s += 1 + end + while zs <= ze + rowval[s] = zrow[zs]; nzval[s] = -oneunit(T); zs += 1; s += 1 + end + while vs <= ve + rowval[s] = vrow[vs]; nzval[s] = dκ; vs += 1; s += 1 + end + # `dg` carries the same diagonal plus the identity the ridge loop added. + rowval[s] = q + nzval[s] = dg[q] - oneunit(T) + end + return SparseMatrixCSC(N, N, colptr, rowval, nzval) +end + +# `y = Symmetric(Cu) * x` for an upper-triangular compressed-column `Cu`. +function _symmul!(y::AbstractVector{T}, Cu::SparseMatrixCSC{T}, x::AbstractVector{T}) where {T} + fill!(y, zero(T)) + rv = rowvals(Cu) + nz = nonzeros(Cu) + for q in axes(Cu, 2) + xq = x[q] + s = zero(T) + for k in nzrange(Cu, q) + p = rv[k] + v = nz[k] + if p == q + s += v * xq + else + y[p] += v * xq + s += v * x[p] + end + end + # Later columns add the remaining terms to `y[q]`. + y[q] += s + end + return y +end + # `AbsLog{2}` penalty continuation. Each stage freezes residual weights, solves # the weighted least-squares problem, and backtracks. `boost=true` applies a final # feasibility shift. The support layout selects the inner solver. @@ -689,9 +767,34 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; cv = zeros(T, ne + 1) # √weight · log|A_ij|, with a trailing 0 gauge target # Violated entries under the current frozen weights. vpat = _violation_pattern(supp) - vedges = Tuple{Int,Int}[] # the violated entries of the current solve + vedges = Tuple{Int,Int}[] # the violated entries of the current solve (LSQR path only) + vrow = Int[] # violated rows, grouped by column (Woodbury path only) + vcnt = zeros(Int, use_woodbury ? N : 0) # violated off-diagonal entries per column + vptr = zeros(Int, use_woodbury ? N + 1 : 0) degV = zeros(Int, use_woodbury ? N : 0) # violated entries per unknown dg = zeros(T, use_woodbury ? N : 0) # diagonal of `B`, for the ridge and the CG preconditioner + # Zero set grouped by column, excluding the diagonal, which `dg` carries. + zptr = zeros(Int, use_woodbury ? N + 1 : 0) + zrow = Int[] + if use_woodbury + for (p, q) in sys.zedges + p == q || (zptr[q+1] += 1) + end + zptr[1] = 1 + cumsum!(zptr, zptr) + resize!(zrow, zptr[end] - 1) + zcursor = zptr[1:end-1] + # `sys.zedges` ordering keeps each compressed column sorted. + for (p, q) in sys.zedges + p == q && continue + zrow[zcursor[q]] = p + zcursor[q] += 1 + end + end + # Storage for the sparse correction, reused across solves. + Ccolptr = zeros(Int, use_woodbury ? N + 1 : 0) + Crowval = Int[] + Cnzval = T[] # Diagonal of the unweighted, gauge-augmented normal matrix. dpart = zeros(T, use_lsqr ? N : 0) if supp isa EdgeList && use_lsqr @@ -719,10 +822,6 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; pg = zeros(T, use_lsqr ? N : 0) # `Rᵀ√W y` before the preconditioner is applied # `K` of the diagonal preconditioner, which is κ-independent and so built once. psqrt = use_lsqr ? sqrt.(dpart) : T[] - # COO triplets of `C`, refilled whenever a Woodbury solve is factorized. - Ci = Int[] - Cj = Int[] - Cv = T[] rhs = zeros(T, use_woodbury ? N : 0, size(U, 2) + 1) dmin = use_woodbury ? minimum(sys.dfull) : oneunit(T) cgx = zeros(T, use_woodbury ? N : 0) @@ -743,8 +842,12 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; fill!(f, zero(T)) copyto!(dg, czero) fill!(degV, 0) - empty!(vedges) - _assemble_woodbury!(f, dg, degV, vedges, vpat, x, κ, supp, symmetric, dκ) + empty!(vrow) + _assemble_woodbury!(f, dg, degV, vrow, vcnt, vpat, x, κ, supp, symmetric, dκ) + vptr[1] = 1 + for q in 1:N + vptr[q+1] = vptr[q] + vcnt[q] + end # Match the ridge used by the dense path. dmax = zero(T) maxdegV = 0 @@ -756,47 +859,13 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; for p in 1:N dg[p] += oneunit(T) + ridge end - # Store full `C` for both matrix-vector products and factorization. - empty!(Ci) - empty!(Cj) - empty!(Cv) - for p in 1:N - push!(Ci, p) - push!(Cj, p) - push!(Cv, czero[p] + ridge) - end - for (p, q) in sys.zedges - p == q && continue - push!(Ci, p) - push!(Cj, q) - push!(Cv, -oneunit(T)) - push!(Ci, q) - push!(Cj, p) - push!(Cv, -oneunit(T)) - end - for (p, q) in vedges - push!(Ci, p) - push!(Cj, p) - push!(Cv, dκ) - push!(Ci, p) - push!(Cj, q) - push!(Cv, dκ) - if q != p - push!(Ci, q) - push!(Cj, q) - push!(Cv, dκ) - push!(Ci, q) - push!(Cj, p) - push!(Cv, dκ) - end - end - C = sparse(Ci, Cj, Cv, N, N) + C = _assemble_C!(Ccolptr, Crowval, Cnzval, zptr, zrow, vptr, vrow, dg, dκ, N) # Use CG while the Gershgorin condition estimate remains small. κest = oneunit(T) + dκ * 2 * maxdegV / dmin if κest <= WOODBURY_CG_KAPPA copyto!(cgx, x) Bmul! = function (yy, xx) - mul!(yy, C, xx) + _symmul!(yy, C, xx) # Indicator columns make the low-rank term a block sum. for k in axes(U, 2) s = zero(T) @@ -817,7 +886,7 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; ok && return copy(cgx) end nchol[] += 1 - return _woodbury_solve!(zeros(T, N), cholesky(Symmetric(C)), U, f, rhs) + return _woodbury_solve!(zeros(T, N), cholesky(Symmetric(C, :U)), U, f, rhs) elseif use_lsqr edges = supp.edges cvals = supp.cvals From 312a7bfcab1a6f76a4b932a17506ff5015dd515e Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 24 Aug 2026 04:16:12 -0500 Subject: [PATCH 04/15] Narrow the boost list and its bucket index Use `Int32` labels for dense boost entries when dimensions fit, and form already-bounded bucket indices without clamping. Assisted-by: Claude Opus (claude-opus-5) --- src/dense_heuristic.jl | 69 ++++++++++++++++++++++++++--------------- src/heuristic_covers.jl | 4 ++- 2 files changed, 47 insertions(+), 26 deletions(-) diff --git a/src/dense_heuristic.jl b/src/dense_heuristic.jl index 18964c6..87aa366 100644 --- a/src/dense_heuristic.jl +++ b/src/dense_heuristic.jl @@ -103,6 +103,48 @@ function _grid_logabs!(L::Matrix{T}, sa::Vector{T}, sb::Vector{T}, return L end +# Use compact boost-list indices when the dimensions fit. +_grid_label(m::Int, n::Int) = max(m, n) <= typemax(Int32) ? Int32 : Int + +# The violated entries of a packed upper triangle, in the traversal order of +# `foreach_support_sym`. Half the entries of a fresh unconstrained start +# violate, so a branch on the test would mispredict on half the grid; the sweep +# selects branchlessly instead, writing every entry to the slot after the last +# one kept and advancing only on a violation. That needs one slot of slack. +function _tri_violated(Lp::Vector{T}, lα::Vector{T}, n::Int, nviol::Int, + ::Type{IT}) where {T,IT} + entries = Vector{Tuple{IT,IT,T}}(undef, nviol + 1) + k = 1 + for jp in 1:n + o = _trioff(jp) + lj = lα[jp] + for ip in 1:jp + lv = Lp[o+ip] + entries[k] = (ip % IT, jp % IT, lv) + k += ifelse(lv - lα[ip] - lj > zero(T), 1, 0) + end + end + resize!(entries, nviol) + return entries +end + +# The violated entries of a full grid, selected as in `_tri_violated`. +function _grid_violated(L::Matrix{T}, lα::Vector{T}, lβ::Vector{T}, m::Int, n::Int, + nviol::Int, ::Type{IT}) where {T,IT} + entries = Vector{Tuple{IT,IT,T}}(undef, nviol + 1) + k = 1 + for jp in 1:n + lj = lβ[jp] + for ip in 1:m + lv = L[ip, jp] + entries[k] = (ip % IT, jp % IT, lv) + k += ifelse(lv - lα[ip] - lj > zero(T), 1, 0) + end + end + resize!(entries, nviol) + return entries +end + # Keep supported scales positive when `exp` underflows. _uncon_scale(si::T, ni::Int, halfmu::T) where {T} = iszero(ni) ? zero(T) : max(exp(si / ni - halfmu), floatmin(T)) @@ -151,19 +193,7 @@ function _symcover_dense!(a::AbstractVector, A::AbstractMatrix, ::Type{T}, maxit # A supported zero scale produces an infinite deficit. isfinite(zmax) || throw(ArgumentError("boost_feasible! requires a start with positive scale on every supported row")) - entries = Vector{Tuple{Int,Int,T}}(undef, nviol + 1) - k = 1 - for jp in 1:n - o = _trioff(jp) - lj = lα[jp] - for ip in 1:jp - lv = Lp[o+ip] - entries[k] = (ip, jp, lv) - k += ifelse(lv - lα[ip] - lj > zero(T), 1, 0) - end - end - resize!(entries, nviol) - _dense_boost!(α, lα, entries, zmax) + _dense_boost!(α, lα, _tri_violated(Lp, lα, n, nviol, _grid_label(n, n)), zmax) lratio = Vector{T}(undef, n) for _ in 1:maxiter @@ -221,7 +251,6 @@ function _cover_dense!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, lβ[jp] = log(β[jp]) end - # Branchless selection, as in the symmetric kernel. nviol = 0 zmax = zero(T) for jp in 1:n @@ -234,17 +263,7 @@ function _cover_dense!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, end isfinite(zmax) || throw(ArgumentError("boost_feasible! requires a start with positive scale on every supported row/column")) - entries = Vector{Tuple{Int,Int,T}}(undef, nviol + 1) - k = 1 - for jp in 1:n - lj = lβ[jp] - for ip in 1:m - lv = L[ip, jp] - entries[k] = (ip, jp, lv) - k += ifelse(lv - lα[ip] - lj > zero(T), 1, 0) - end - end - resize!(entries, nviol) + entries = _grid_violated(L, lα, lβ, m, n, nviol, _grid_label(m, n)) # Row and column scales require separate updates. deficit((i, j, lv)) = lv - lα[i] - lβ[j] function apply!((i, j, lv), z) diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index d0dd9a4..5e71a4e 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -383,7 +383,9 @@ function bucket_boost!(deficit::F, apply!::G, entries::AbstractVector, ::Type{T} zmax > zero(T) || return w = T(BOOST_BUCKET_WIDTH) B = max(1, ceil(Int, zmax / w)) - bucketof(z) = clamp(ceil(Int, z / w), 1, B) + # Deficits only fall, so every `z` reaching `bucketof` lies in `(0, zmax]` + # and `ceil(z / w)` lies in `1:B`, a range the truncation cannot leave. + bucketof(z) = clamp(unsafe_trunc(Int, ceil(z / w)), 1, B) ptr = zeros(Int, B + 1) for entry in entries z = deficit(entry) From 8b694e467c15ca537ae94b0e7b19f15e4f8efabd Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 24 Aug 2026 04:16:12 -0500 Subject: [PATCH 05/15] Read the balance components off the grid Use grid support counts to avoid a second traversal when support is full; retain component traversal otherwise. Assisted-by: Claude Opus (claude-opus-5) --- src/dense_heuristic.jl | 10 +++++++++- src/heuristic_covers.jl | 9 ++++++++- src/support.jl | 19 ++++--------------- 3 files changed, 21 insertions(+), 17 deletions(-) diff --git a/src/dense_heuristic.jl b/src/dense_heuristic.jl index 87aa366..3b315c7 100644 --- a/src/dense_heuristic.jl +++ b/src/dense_heuristic.jl @@ -103,6 +103,14 @@ function _grid_logabs!(L::Matrix{T}, sa::Vector{T}, sb::Vector{T}, return L end +# Derive the balance summary directly for full support; traverse sparse support. +function _grid_components(A::AbstractMatrix, na::Vector{Int}, nb::Vector{Int}, m::Int, n::Int) + for ip in 1:m + na[ip] == n || return _support_components(A) + end + return ones(Int, m), ones(Int, n), 1, na, nb +end + # Use compact boost-list indices when the dimensions fit. _grid_label(m::Int, n::Int) = max(m, n) <= typemax(Int32) ? Int32 : Int @@ -301,7 +309,7 @@ function _cover_dense!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, b[jp+oc] = β[jp] end # Balance the factors, then restore coverage lost to rounding. - _balance_cover!(a, b, A) + _balance_cover!(a, b, _grid_components(A, na, nb, m, n)...) for ip in 1:m lα[ip] = log(T(a[ip+or])) end diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index 5e71a4e..32b1a01 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -174,8 +174,15 @@ end # component. Rounding the shift to a power of two preserves cover products # exactly, at the cost of balancing only within a factor of `sqrt(2)`. function _balance_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) - T = float(promote_type(eltype(a), eltype(b))) rowcomp, colcomp, ncomp, nzrow, nzcol = _support_components(A) + return _balance_cover!(a, b, rowcomp, colcomp, ncomp, nzrow, nzcol) +end + +# Balance from precomputed component labels and support counts. +function _balance_cover!(a::AbstractVector, b::AbstractVector, rowcomp::Vector{Int}, + colcomp::Vector{Int}, ncomp::Int, nzrow::Vector{Int}, + nzcol::Vector{Int}) + T = float(promote_type(eltype(a), eltype(b))) iszero(ncomp) && return a, b Lα = zeros(T, ncomp) Lβ = zeros(T, ncomp) diff --git a/src/support.jl b/src/support.jl index e9c47a7..30915aa 100644 --- a/src/support.jl +++ b/src/support.jl @@ -133,21 +133,10 @@ end # invariant already guarantees `abs(A[i,j]) == abs(A[j,i])`. require_abs_symmetric(::Union{Symmetric,Hermitian,Diagonal,SymTridiagonal}, fname) = nothing -# Connected components of the bipartite support graph of `A`: one vertex per row -# and one per column, one edge per stored nonzero. Returns `(rowcomp, colcomp, -# ncomp)`, where `rowcomp` and `colcomp` are `Vector{Int}` indexed by *position* -# within `axes(A, 1)` and `axes(A, 2)` (so offset axes need no special case, as -# with `GroupedSupport.ptr`), holding the component id in `1:ncomp` — or 0 for -# rows/columns with empty support, which belong to no component. `nzrow` and -# `nzcol` count the stored entries of each row and column, in the same position -# space. -# -# The gauge orbit of an asymmetric cover has one dimension per component: the -# rescaling `a -> γ*a`, `b -> b/γ` acts independently on each, because no -# product `a[i]*b[j]` spans two components. Any convention that pins the split -# between `a` and `b` must therefore be imposed per component; a single global -# constraint leaves `ncomp - 1` directions to the whim of whichever pass ran -# last. +# Connected components of the bipartite support graph. Labels and support +# counts use positions within each axis; unsupported rows and columns have label +# zero. Each component has an independent `a -> γ*a`, `b -> b/γ` gauge, so +# balancing must also be per component. function _support_components(A::AbstractMatrix) m = length(axes(A, 1)) n = length(axes(A, 2)) From f15309acda8aee6561880049a734f693b7f29180 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 24 Aug 2026 04:24:54 -0500 Subject: [PATCH 06/15] Form the deficit bucket index unbounded Deficits passed to `bucketof` are already in range, so compute the index without clamping. Assisted-by: Claude Opus (claude-opus-5) --- src/heuristic_covers.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index 32b1a01..2495a32 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -390,9 +390,9 @@ function bucket_boost!(deficit::F, apply!::G, entries::AbstractVector, ::Type{T} zmax > zero(T) || return w = T(BOOST_BUCKET_WIDTH) B = max(1, ceil(Int, zmax / w)) - # Deficits only fall, so every `z` reaching `bucketof` lies in `(0, zmax]` - # and `ceil(z / w)` lies in `1:B`, a range the truncation cannot leave. - bucketof(z) = clamp(unsafe_trunc(Int, ceil(z / w)), 1, B) + # `z` is a positive log difference no greater than `zmax`; its spacing keeps + # `z / w` from underflowing, so the ceiling remains in `1:B`. + bucketof(z) = unsafe_trunc(Int, ceil(z / w)) ptr = zeros(Int, B + 1) for entry in entries z = deficit(entry) From 5cebd2d6f6bcbb2b3ca7f3baf0fc86bc08dfb5b4 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 24 Aug 2026 05:24:53 -0500 Subject: [PATCH 07/15] Certify hard covers against exact iscover Certify returned hard covers in `iscover` arithmetic by uniformly inflating factors for the worst rounding shortfall. Leave initializer starts unchanged. Assisted-by: Claude Opus (claude-opus-5) --- src/heuristic_covers.jl | 31 +++++++----- src/iscover.jl | 106 ++++++++++++++++++++++++++++++++++++++++ src/minimal_covers.jl | 12 ++++- test/invariants.jl | 33 +++++++++++++ 4 files changed, 169 insertions(+), 13 deletions(-) diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index 2495a32..5166e8f 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -69,10 +69,15 @@ end function _symcover!(a::AbstractVector, A::AbstractMatrix; maxiter::Int=3) T = float(real(eltype(a))) - _use_dense_grid(A, T) && return _symcover_dense!(a, A, T, maxiter) - unconstrained_min!(AbsLog{2}(), a, A) - boost_feasible!(a, A) - return tighten_cover!(a, A; maxiter) + if _use_dense_grid(A, T) + _symcover_dense!(a, A, T, maxiter) + else + unconstrained_min!(AbsLog{2}(), a, A) + boost_feasible!(a, A) + tighten_cover!(a, A; maxiter) + end + # Certify against `A` after log-domain tightening. + return _certify_cover!(a, A, :symcover) end """ @@ -148,13 +153,17 @@ end function _cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=3) T = float(promote_type(eltype(a), eltype(b))) - _use_dense_grid(A, T) && return _cover_dense!(a, b, A, T, maxiter) - unconstrained_min!(AbsLog{2}(), a, b, A) - boost_feasible!(a, b, A) - tighten_cover!(a, b, A; maxiter) - # Apply the package's balance convention, then restore coverage lost to rounding. - _balance_cover!(a, b, A) - return inflate_feasible!(a, b, A) + if _use_dense_grid(A, T) + _cover_dense!(a, b, A, T, maxiter) + else + unconstrained_min!(AbsLog{2}(), a, b, A) + boost_feasible!(a, b, A) + tighten_cover!(a, b, A; maxiter) + # Apply the package's balance convention, then restore coverage lost to rounding. + _balance_cover!(a, b, A) + inflate_feasible!(a, b, A) + end + return _certify_cover!(a, b, A, :cover) end # Adjoint/Transpose wrappers for cover!. diff --git a/src/iscover.jl b/src/iscover.jl index 52c48b7..3a23ada 100644 --- a/src/iscover.jl +++ b/src/iscover.jl @@ -14,6 +14,10 @@ requires `A` to be square. Both tolerances default to zero. A nonzero `atol` breaks scale invariance. +`cover`, `symcover`, and native `AbsLog{2}` minimal covers certify their results +at zero tolerance. Initializers and extension solvers may require a nonzero +`rtol`. + `a` and `b` must be nonnegative; a negative scale raises an `ArgumentError`. Zero is allowed for unsupported rows and columns. @@ -68,3 +72,105 @@ function _require_nonneg(x::AbstractVector, name::String) end return nothing end + +# Log-domain solvers can lose coverage to rounding. Measure the largest +# linear-arithmetic shortfall and apply a uniform inflation without changing the +# balance convention. +const CERTIFY_SWEEPS = 4 + +function _certify_cover!(a::AbstractVector, A::AbstractMatrix, fname::Symbol) + T = scalar_type(eltype(a)) + for _ in 1:CERTIFY_SWEEPS + r = _worst_shortfall(a, A, T, fname) + r > one(T) || return a + _inflate_nonzero!(a, _certify_factor(r)) + end + throw(ArgumentError("$fname could not certify a cover of `A` within $CERTIFY_SWEEPS inflation sweeps")) +end + +function _certify_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, fname::Symbol) + T = scalar_type(promote_type(eltype(a), eltype(b))) + for _ in 1:CERTIFY_SWEEPS + r = _worst_shortfall(a, b, A, T, fname) + r > one(T) || return a, b + s = _certify_factor(r) + _inflate_nonzero!(a, s) + _inflate_nonzero!(b, s) + end + throw(ArgumentError("$fname could not certify a cover of `A` within $CERTIFY_SWEEPS inflation sweeps")) +end + +# How far short of `v` the cover product `p` falls, as the factor `p` must grow +# by. A vanishing or non-finite product cannot be lifted onto a positive entry. +function _shortfall(p, v, i, j, fname::Symbol) + r = v / p + isfinite(r) || + throw(ArgumentError("$fname requires a positive, finite cover product on every supported entry, got $(string(p)) at ($(string(i)), $(string(j)))")) + return r +end + +# Round each factor's share of the required inflation upward. +_certify_factor(r::T) where {T} = max(nextfloat(sqrt(r)), nextfloat(one(T))) + +function _inflate_nonzero!(x::AbstractVector, s) + for i in eachindex(x) + iszero(x[i]) || (x[i] *= s) + end + return x +end + +# Worst factor by which a cover product must grow to reach its entry, or `one` +# when the cover already holds everywhere. +function _worst_shortfall(a::AbstractVector, A::AbstractMatrix, ::Type{T}, fname::Symbol) where {T} + worst = Ref(one(T)) + foreach_support_sym(A) do i, j, v + p = a[i] * a[j] + p >= v && return + worst[] = max(worst[], convert(T, _shortfall(p, v, i, j, fname))) + end + return worst[] +end + +function _worst_shortfall(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, + ::Type{T}, fname::Symbol) where {T} + worst = Ref(one(T)) + foreach_support(A) do i, j, v + p = a[i] * b[j] + p >= v && return + worst[] = max(worst[], convert(T, _shortfall(p, v, i, j, fname))) + end + return worst[] +end + +# Direct dense-storage implementations avoid callback state. +function _worst_shortfall(a::AbstractVector, A::StridedMatrix, ::Type{T}, fname::Symbol) where {T} + worst = one(T) + ax = axes(A, 1) + for j in ax + aj = a[j] + for i in first(ax):j + v = abs(A[i, j]) + iszero(v) && continue + p = a[i] * aj + p >= v && continue + worst = max(worst, convert(T, _shortfall(p, v, i, j, fname))) + end + end + return worst +end + +function _worst_shortfall(a::AbstractVector, b::AbstractVector, A::StridedMatrix, + ::Type{T}, fname::Symbol) where {T} + worst = one(T) + for j in axes(A, 2) + bj = b[j] + for i in axes(A, 1) + v = abs(A[i, j]) + iszero(v) && continue + p = a[i] * bj + p >= v && continue + worst = max(worst, convert(T, _shortfall(p, v, i, j, fname))) + end + end + return worst +end diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index f619ef5..f127ecc 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -1103,7 +1103,10 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), if eps(T) > eps(Float64) a64, stats = _symcover_min_abslog2(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A); κs, maxiter, linsolve, start, boost, fname) - return T.(a64), stats + # Narrowing rounds to nearest and so can round a product below its entry. + a = T.(a64) + boost && _certify_cover!(a, A, fname) + return a, stats end n = length(ax) use_lsqr = linsolve === :lsqr @@ -1187,6 +1190,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for (ip, i) in enumerate(ax) a[i] = hassupp[ip] ? exp(α[ip]) : zero(T) end + boost && _certify_cover!(a, A, fname) return a, stats end @@ -1204,7 +1208,10 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), if eps(T) > eps(Float64) a64, b64, stats = _cover_min_abslog2(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A); κs, maxiter, linsolve, start, boost) - return T.(a64), T.(b64), stats + # Narrowing rounds to nearest and so can round a product below its entry. + a, b = T.(a64), T.(b64) + boost && _certify_cover!(a, b, A, :cover_min) + return a, b, stats end m = length(axr) n = length(axc) @@ -1342,6 +1349,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for (jp, j) in enumerate(axc) b[j] = hascol[jp] ? exp(x[m+jp] - s[colcomp[jp]]) : zero(T) end + boost && _certify_cover!(a, b, A, :cover_min) return a, b, stats end diff --git a/test/invariants.jl b/test/invariants.jl index cf3c0a2..5a40780 100644 --- a/test/invariants.jl +++ b/test/invariants.jl @@ -172,3 +172,36 @@ const GEN_NOTIONS = ( end end end + +@testset "exact coverage" begin + # Native hard-cover results are certified at zero tolerance. + certified_sym = ( + A -> symcover(A), + A -> symcover(A; maxiter=0), + A -> symcover_min(AbsLog{2}(), A), + ) + certified_gen = ( + A -> cover(A), + A -> cover(A; maxiter=0), + A -> cover_min(AbsLog{2}(), A), + ) + # Exercise both flattened-support and dense-grid kernels. + for T in (Float64, Float32), σ in (1, 5), n in (7, 70), seed in 1:3 + rng = StableRNG(97 * seed + 13 * n + 3 * σ + (T === Float32)) + G = randn(rng, n, n) + Msym = T.(exp.((σ / 2) .* (G .+ G'))) + Mgen = T.(exp.(σ .* randn(rng, n, n))) + Mrec = T.(exp.(σ .* randn(rng, n, n ÷ 2 + 1))) + # Zeros exercise the sparse-support traversals and the component split. + drop = rand(rng, n, n) .< 0.3 + Zsym = copy(Msym); Zsym[drop .| drop'] .= 0 + Zgen = copy(Mgen); Zgen[drop] .= 0 + for f in certified_sym, A in (Msym, Zsym, sparse(Zsym)) + @test iscover(f(A), A) + end + for f in certified_gen, A in (Mgen, Zgen, sparse(Zgen), Mrec) + a, b = f(A) + @test iscover(a, b, A) + end + end +end From 095a00348064736f0bddf4b034a8536ebd2e2aa7 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 24 Aug 2026 05:51:29 -0500 Subject: [PATCH 08/15] Report continuation stages that run out of steps Record each continuation stage's exit reason and final decrease; warn when it reaches `maxiter` while still descending. Assisted-by: Claude Opus (claude-opus-5) --- src/minimal_covers.jl | 42 ++++++++++++++++++++++++++++++++++++------ test/minimal_covers.jl | 27 +++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 6 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index f127ecc..f0d7715 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -35,6 +35,9 @@ The native solver accepts `κs` (penalty-continuation schedule), `maxiter` - `:auto` chooses `:woodbury` when supported, `:lsqr` when the stored support fills at most a quarter of the grid, and `:dense` otherwise. +If a stage reaches `maxiter`, the solver warns that the cover may not minimize +the objective. Increase `maxiter` or supply more continuation stages. + The native solver computes in `Float64` for narrower input types, then converts the result to the required element type. @@ -66,8 +69,9 @@ method selects the one with the smallest `AbsLog{2}` objective. # Extended help The native solver accepts the same `κs`, `maxiter`, and `linsolve` keywords as -[`symcover_min`](@ref). For `:woodbury`, an `m × n` matrix may omit at most -`min(m,n) ÷ 4` entries per row or column and `4 * max(m,n)` entries in total. +[`symcover_min`](@ref), and the same warning when a stage runs out of Newton +steps. For `:woodbury`, an `m × n` matrix may omit at most `min(m,n) ÷ 4` entries +per row or column and `4 * max(m,n)` entries in total. `:dense` costs O((m+n)³) per Newton step; sparse matrices default to `:lsqr`. See also: [`symcover_min`](@ref), [`cover`](@ref), [`cover_min!`](@ref). @@ -736,6 +740,16 @@ function _symmul!(y::AbstractVector{T}, Cu::SparseMatrixCSC{T}, x::AbstractVecto return y end +# Warn when a continuation stage reaches `maxiter` while still descending. +function _warn_truncated(fname::Symbol, κs, stats, maxiter::Int) + exits = stats.exits + any(==(:maxiter), exits) || return nothing + stalled = [(k, κs[k], stats.stagedrops[k]) for k in eachindex(exits) if exits[k] === :maxiter] + detail = join(("stage $k (κ = $κ) was still decreasing by $d per step" for (k, κ, d) in stalled), "; ") + @warn "$fname: $(length(stalled)) of $(length(exits)) continuation stages reached maxiter=$maxiter; the result covers `A` but may not minimize the objective ($detail). Increase `maxiter` or supply more `κs` stages." + return nothing +end + # `AbsLog{2}` penalty continuation. Each stage freezes residual weights, solves # the weighted least-squares problem, and backtracks. `boost=true` applies a final # feasibility shift. The support layout selects the inner solver. @@ -1054,8 +1068,13 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; end end x = x0 === nothing ? solve_weighted(zeros(T, N), nothing) : x0 - for κ in κs + # Record each stage's exit reason and final relative decrease. + exits = Vector{Symbol}(undef, length(κs)) + drops = Vector{T}(undef, length(κs)) + for (k, κ) in enumerate(κs) fcur = _fκ(x, κ, supp, symmetric) + exit = :maxiter + drop = zero(T) for _ in 1:maxiter xnew = solve_weighted(x, κ) t = one(T) @@ -1068,11 +1087,20 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; stable = false end x = xt + drop = (fcur - fnew) / max(fcur, one(T)) # For exact inner solves, an unchanged violation pattern ends the stage. - !use_lsqr && stable && break - fcur - fnew <= 5000 * eps(T) * max(fcur, one(T)) && break + if !use_lsqr && stable + exit = :stable + break + end + if fcur - fnew <= 5000 * eps(T) * max(fcur, one(T)) + exit = :decrease + break + end fcur = fnew end + exits[k] = exit + drops[k] = drop end # Hard covers receive a final uniform feasibility shift. if boost @@ -1082,7 +1110,7 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; end end return x, (; nsolves=nsolves[], lsqriters=nlsqr[], cgiters=ncg[], - cholsolves=nchol[], + cholsolves=nchol[], exits=Tuple(exits), stagedrops=Tuple(drops), linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense)) end @@ -1185,6 +1213,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), x0 = start === nothing ? nothing : T[hassupp[ip] ? log(T(start[i])) : zero(T) for (ip, i) in enumerate(ax)] α, stats = _abslog2_continuation(sys, x0; κs, maxiter, linsolve, boost) + _warn_truncated(fname, κs, stats, maxiter) # Dense scale vector matching cover/symcover; `similar(A, …)` is a SparseVector for sparse A. a = similar(Array{T}, ax) for (ip, i) in enumerate(ax) @@ -1320,6 +1349,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), s0 end x, stats = _abslog2_continuation(sys, x0; κs, maxiter, linsolve, boost) + _warn_truncated(:cover_min, κs, stats, maxiter) # Apply the balance convention independently to each support component. rowcomp, colcomp, ncomp, _, _ = _support_components(A) Lα = zeros(T, ncomp) diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 31af75d..53fc813 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -575,3 +575,30 @@ end # Allow iteration growth while rejecting an added dense workspace. @test large < 10 * small end + +@testset "MMC continuation reports how each stage ended" begin + rng = StableRNG(17) + A = (X = exp.(randn(rng, 40, 40)); (X .+ X') ./ 2) + _, s = MatrixCovers._symcover_min_abslog2(A) + @test length(s.exits) == 4 + @test length(s.stagedrops) == length(s.exits) + @test all(in((:stable, :decrease, :maxiter)), s.exits) + # With room to converge, no stage runs out of Newton steps. + @test all(!=(:maxiter), s.exits) + + # A one-step limit truncates every `:lsqr` stage and emits a warning. + @test_logs (:warn, r"reached maxiter=1") match_mode=:any begin + _, sl = MatrixCovers._symcover_min_abslog2(A; maxiter=1, linsolve=:lsqr) + @test all(==(:maxiter), sl.exits) + end + + G = exp.(randn(rng, 40, 30)) + _, _, t = MatrixCovers._cover_min_abslog2(G) + @test length(t.exits) == 4 + @test all(!=(:maxiter), t.exits) + + # A supplied schedule sets how many stages are reported. + _, s8 = MatrixCovers._symcover_min_abslog2(A; κs=10 .^ range(2, 8, length=8)) + @test length(s8.exits) == 8 + @test all(!=(:maxiter), s8.exits) +end From 35c268216217bc93c33250a1306dfebe57e85fd8 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 24 Aug 2026 06:01:27 -0500 Subject: [PATCH 09/15] Set the continuation schedule from the solver path MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Build working-precision schedules ending at `1e8`: eight stages for exact solvers and four for `:lsqr`. An explicit `κs` overrides the default. Assisted-by: Claude Opus (claude-opus-5) --- src/minimal_covers.jl | 39 +++++++++++++++++++++++++++------------ test/minimal_covers.jl | 40 ++++++++++++++++++++++++++++++---------- 2 files changed, 57 insertions(+), 22 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index f0d7715..171245d 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -35,6 +35,9 @@ The native solver accepts `κs` (penalty-continuation schedule), `maxiter` - `:auto` chooses `:woodbury` when supported, `:lsqr` when the stored support fills at most a quarter of the grid, and `:dense` otherwise. +`κs` defaults to a geometric schedule ending at `1e8`: eight stages for exact +solves and four for `:lsqr`. An explicit `κs` overrides this default. + If a stage reaches `maxiter`, the solver warns that the cover may not minimize the objective. Increase `maxiter` or supply more continuation stages. @@ -69,9 +72,10 @@ method selects the one with the smallest `AbsLog{2}` objective. # Extended help The native solver accepts the same `κs`, `maxiter`, and `linsolve` keywords as -[`symcover_min`](@ref), and the same warning when a stage runs out of Newton -steps. For `:woodbury`, an `m × n` matrix may omit at most `min(m,n) ÷ 4` entries -per row or column and `4 * max(m,n)` entries in total. +[`symcover_min`](@ref), with the same solver-dependent default schedule and the +same warning when a stage runs out of Newton steps. For `:woodbury`, an `m × n` +matrix may omit at most `min(m,n) ÷ 4` entries per row or column and +`4 * max(m,n)` entries in total. `:dense` costs O((m+n)³) per Newton step; sparse matrices default to `:lsqr`. See also: [`symcover_min`](@ref), [`cover`](@ref), [`cover_min!`](@ref). @@ -261,9 +265,8 @@ end # - `: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. +# Maximum support density for the `:auto` LSQR path. At higher densities the +# exact dense solve has lower overhead and can stop on a sign-stable step. const AUTO_LSQR_MAX_DENSITY = 1 // 4 # Condition estimate above which Woodbury uses sparse Cholesky instead of CG. @@ -740,6 +743,14 @@ function _symmul!(y::AbstractVector{T}, Cu::SparseMatrixCSC{T}, x::AbstractVecto return y end +# Geometric continuation schedules ending at `1e8`. An exact solve ends a stage +# on the first sign-stable Newton step, so its finer eight-stage schedule costs +# about one extra solve per added stage; `:lsqr` has no such exit, pays a full +# descent per stage, and keeps four. +_kappa_schedule(::Type{T}, use_lsqr::Bool) where {T} = + use_lsqr ? ntuple(k -> T(10)^(2k), 4) : + ntuple(k -> T(10)^(T(2) + T(6) * T(k - 1) / T(7)), 8) + # Warn when a continuation stage reaches `maxiter` while still descending. function _warn_truncated(fname::Symbol, κs, stats, maxiter::Int) exits = stats.exits @@ -1116,7 +1127,7 @@ end # Worker for `symcover_min(::AbsLog{2})`, returning `(a, stats)`. A supplied # `start` replaces the cold initial solve. Narrow types compute in `Float64`. -function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), +function _symcover_min_abslog2(A::AbstractMatrix; κs=nothing, maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, boost::Bool=true, fname=:symcover_min) linsolve in (:auto, :dense, :lsqr, :woodbury) || @@ -1173,6 +1184,8 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), use_lsqr = true linsolve = :lsqr end + # Select the default schedule after selecting the solver. + κsched = κs === nothing ? _kappa_schedule(T, use_lsqr) : κs # Woodbury uses a grid; dense and LSQR use an edge list. supp = if use_woodbury C = fill(T(-Inf), n, n) @@ -1212,8 +1225,8 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), zeros(T, n)) x0 = start === nothing ? nothing : T[hassupp[ip] ? log(T(start[i])) : zero(T) for (ip, i) in enumerate(ax)] - α, stats = _abslog2_continuation(sys, x0; κs, maxiter, linsolve, boost) - _warn_truncated(fname, κs, stats, maxiter) + α, stats = _abslog2_continuation(sys, x0; κs=κsched, maxiter, linsolve, boost) + _warn_truncated(fname, κsched, stats, maxiter) # Dense scale vector matching cover/symcover; `similar(A, …)` is a SparseVector for sparse A. a = similar(Array{T}, ax) for (ip, i) in enumerate(ax) @@ -1224,7 +1237,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end # Worker for `cover_min(::AbsLog{2})`, returning `(a, b, stats)`. -function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), +function _cover_min_abslog2(A::AbstractMatrix; κs=nothing, maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, boost::Bool=true) linsolve in (:auto, :dense, :lsqr, :woodbury) || @@ -1287,6 +1300,8 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), use_lsqr = true linsolve = :lsqr end + # Select the default schedule after selecting the solver. + κsched = κs === nothing ? _kappa_schedule(T, use_lsqr) : κs # Woodbury uses a grid; dense and LSQR use an edge list. supp = if use_woodbury C = fill(T(-Inf), m, n) @@ -1348,8 +1363,8 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end s0 end - x, stats = _abslog2_continuation(sys, x0; κs, maxiter, linsolve, boost) - _warn_truncated(:cover_min, κs, stats, maxiter) + x, stats = _abslog2_continuation(sys, x0; κs=κsched, maxiter, linsolve, boost) + _warn_truncated(:cover_min, κsched, stats, maxiter) # Apply the balance convention independently to each support component. rowcomp, colcomp, ncomp, _, _ = _support_components(A) Lα = zeros(T, ncomp) diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 53fc813..7e876c7 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -330,10 +330,14 @@ end al, sl = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) @test ad ≈ al rtol=1e-6 @test aw ≈ al rtol=1e-6 - # Exact paths save one solve per continuation stage. @test sd.nsolves == sw.nsolves - @test sd.nsolves <= sl.nsolves - length((1e2, 1e4, 1e6, 1e8)) @test sd.nsolves <= 26 + # Exact paths save one solve per continuation stage. The default schedule is + # solver-dependent, so the comparison runs both paths on one schedule. + κs8 = MatrixCovers._kappa_schedule(Float64, false) + _, sd8 = MatrixCovers._symcover_min_abslog2(A; linsolve=:dense, κs=κs8) + _, sl8 = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, κs=κs8) + @test sd8.nsolves <= sl8.nsolves - length(κs8) G = exp.(randn(rng, 60, 45)) gd, hd, td = MatrixCovers._cover_min_abslog2(G; linsolve=:dense) @@ -342,9 +346,11 @@ end @test gd .* hd' ≈ gl .* hl' rtol=1e-6 @test gw .* hw' ≈ gl .* hl' rtol=1e-6 @test td.nsolves == tw.nsolves - @test td.nsolves <= tl.nsolves - length((1e2, 1e4, 1e6, 1e8)) # Leave a small margin in the solve-count bound. - @test td.nsolves <= 24 + @test td.nsolves <= 28 + _, _, td8 = MatrixCovers._cover_min_abslog2(G; linsolve=:dense, κs=κs8) + _, _, tl8 = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr, κs=κs8) + @test td8.nsolves <= tl8.nsolves - length(κs8) end # The Float64 LSQR preconditioner includes κ-weighted rows; other types use the @@ -580,7 +586,7 @@ end rng = StableRNG(17) A = (X = exp.(randn(rng, 40, 40)); (X .+ X') ./ 2) _, s = MatrixCovers._symcover_min_abslog2(A) - @test length(s.exits) == 4 + @test length(s.exits) == 8 @test length(s.stagedrops) == length(s.exits) @test all(in((:stable, :decrease, :maxiter)), s.exits) # With room to converge, no stage runs out of Newton steps. @@ -594,11 +600,25 @@ end G = exp.(randn(rng, 40, 30)) _, _, t = MatrixCovers._cover_min_abslog2(G) - @test length(t.exits) == 4 + @test length(t.exits) == 8 @test all(!=(:maxiter), t.exits) - # A supplied schedule sets how many stages are reported. - _, s8 = MatrixCovers._symcover_min_abslog2(A; κs=10 .^ range(2, 8, length=8)) - @test length(s8.exits) == 8 - @test all(!=(:maxiter), s8.exits) + # Exact solves default to eight stages; `:lsqr` defaults to four. + _, sl = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) + @test length(sl.exits) == 4 + # Check the default `Float64` LSQR schedule exactly. + @test MatrixCovers._kappa_schedule(Float64, true) === (1e2, 1e4, 1e6, 1e8) + @test length(MatrixCovers._kappa_schedule(Float64, false)) == 8 + @test MatrixCovers._kappa_schedule(Float64, false)[1] == 1e2 + @test MatrixCovers._kappa_schedule(Float64, false)[end] == 1e8 + @test issorted(MatrixCovers._kappa_schedule(Float64, false)) + # The schedule uses the working precision. + @test eltype(MatrixCovers._kappa_schedule(BigFloat, false)) === BigFloat + @test eltype(MatrixCovers._kappa_schedule(BigFloat, true)) === BigFloat + + # Explicit schedules override the defaults. + _, s4 = MatrixCovers._symcover_min_abslog2(A; κs=(1e2, 1e4, 1e6, 1e8)) + @test length(s4.exits) == 4 + _, s4l = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, κs=10 .^ range(2, 8, length=8)) + @test length(s4l.exits) == 8 end From 26d4b97befd0e04b1cd632d873aca3757c177267 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 24 Aug 2026 08:10:03 -0500 Subject: [PATCH 10/15] Open the LSQR continuation from the heuristic Start LSQR continuation from the heuristic cover; exact solvers retain the unweighted start. An empty schedule still returns the unweighted fit. Assisted-by: Claude Opus (claude-opus-5) --- src/minimal_covers.jl | 9 +++++++++ test/minimal_covers.jl | 39 +++++++++++++++++++++++++++++++++------ 2 files changed, 42 insertions(+), 6 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 171245d..2ac2f6b 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -1186,6 +1186,11 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=nothing, end # Select the default schedule after selecting the solver. κsched = κs === nothing ? _kappa_schedule(T, use_lsqr) : κs + # Start LSQR continuation from the heuristic cover. An empty schedule returns + # the unweighted fit instead. + if start === nothing && use_lsqr && !isempty(κsched) + start = symcover(A) + end # Woodbury uses a grid; dense and LSQR use an edge list. supp = if use_woodbury C = fill(T(-Inf), n, n) @@ -1302,6 +1307,10 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=nothing, end # Select the default schedule after selecting the solver. κsched = κs === nothing ? _kappa_schedule(T, use_lsqr) : κs + # Start LSQR continuation from the heuristic cover. + if start === nothing && use_lsqr && !isempty(κsched) + start = cover(A) + 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/test/minimal_covers.jl b/test/minimal_covers.jl index 7e876c7..1309927 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -332,11 +332,14 @@ end @test aw ≈ al rtol=1e-6 @test sd.nsolves == sw.nsolves @test sd.nsolves <= 26 - # Exact paths save one solve per continuation stage. The default schedule is - # solver-dependent, so the comparison runs both paths on one schedule. + # Exact paths save one solve per continuation stage. The comparison runs both + # paths on one schedule, because the default is solver-dependent, and from one + # start, because the LSQR path otherwise supplies its own: only a shared + # starting iterate leaves the stage-exit rule as the difference between them. κs8 = MatrixCovers._kappa_schedule(Float64, false) - _, sd8 = MatrixCovers._symcover_min_abslog2(A; linsolve=:dense, κs=κs8) - _, sl8 = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, κs=κs8) + a0 = symcover(A) + _, sd8 = MatrixCovers._symcover_min_abslog2(A; linsolve=:dense, κs=κs8, start=a0) + _, sl8 = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, κs=κs8, start=a0) @test sd8.nsolves <= sl8.nsolves - length(κs8) G = exp.(randn(rng, 60, 45)) @@ -348,11 +351,35 @@ end @test td.nsolves == tw.nsolves # Leave a small margin in the solve-count bound. @test td.nsolves <= 28 - _, _, td8 = MatrixCovers._cover_min_abslog2(G; linsolve=:dense, κs=κs8) - _, _, tl8 = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr, κs=κs8) + g0, h0 = cover(G) + _, _, td8 = MatrixCovers._cover_min_abslog2(G; linsolve=:dense, κs=κs8, start=(g0, h0)) + _, _, tl8 = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr, κs=κs8, start=(g0, h0)) @test td8.nsolves <= tl8.nsolves - length(κs8) end +# LSQR continuation starts from the heuristic cover. +@testset "MMC :lsqr continuation starts from the heuristic cover" begin + rng = StableRNG(17) + A = (X = exp.(randn(rng, 50, 50)); (X .+ X') ./ 2) + al, sl = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) + ah, sh = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, start=symcover(A)) + @test al == ah + @test (sl.nsolves, sl.lsqriters) == (sh.nsolves, sh.lsqriters) + ad, _ = MatrixCovers._symcover_min_abslog2(A; linsolve=:dense) + @test ad ≈ al rtol=1e-6 + + G = exp.(randn(rng, 40, 30)) + gl, hl, tl = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr) + g0, h0 = cover(G) + gh, hh, th = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr, start=(g0, h0)) + @test gl == gh + @test hl == hh + @test (tl.nsolves, tl.lsqriters) == (th.nsolves, th.lsqriters) + + # With no stages the unweighted fit is the answer, not a start. + @test soft_symcover_min(AbsLog{2}(), A; linsolve=:lsqr) != symcover(A) +end + # The Float64 LSQR preconditioner includes κ-weighted rows; other types use the # plain matrix-free iteration. @testset "MMC :lsqr iteration count is bounded across the continuation" begin From 0debd4e41d38eb31966dc67dd72f5058ed514aba Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 24 Aug 2026 11:04:55 -0500 Subject: [PATCH 11/15] Share one log per entry across the cover sweeps Cache one log magnitude per support entry across heuristic sweeps while preserving traversal order and results. The dense-grid path is unchanged. Assisted-by: Claude Opus (claude-opus-5) --- src/heuristic_covers.jl | 106 +++++++++++++++++++++++++++++----------- 1 file changed, 77 insertions(+), 29 deletions(-) diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index 5166e8f..46a9abe 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -72,9 +72,10 @@ function _symcover!(a::AbstractVector, A::AbstractMatrix; maxiter::Int=3) if _use_dense_grid(A, T) _symcover_dense!(a, A, T, maxiter) else - unconstrained_min!(AbsLog{2}(), a, A) - boost_feasible!(a, A) - tighten_cover!(a, A; maxiter) + logs = support_logs_sym(A, eltype(a)) + unconstrained_min!(AbsLog{2}(), a, A; logs) + boost_feasible!(a, A; logs) + tighten_cover!(a, A; maxiter, logs) end # Certify against `A` after log-domain tightening. return _certify_cover!(a, A, :symcover) @@ -156,12 +157,13 @@ function _cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxite if _use_dense_grid(A, T) _cover_dense!(a, b, A, T, maxiter) else - unconstrained_min!(AbsLog{2}(), a, b, A) - boost_feasible!(a, b, A) - tighten_cover!(a, b, A; maxiter) + logs = support_logs(A, T) + unconstrained_min!(AbsLog{2}(), a, b, A; logs) + boost_feasible!(a, b, A; logs) + tighten_cover!(a, b, A; maxiter, logs) # Apply the package's balance convention, then restore coverage lost to rounding. _balance_cover!(a, b, A) - inflate_feasible!(a, b, A) + inflate_feasible!(a, b, A; logs) end return _certify_cover!(a, b, A, :cover) end @@ -179,6 +181,32 @@ end # ============================================================ # Internal helpers # ============================================================ +# Logarithms of the stored magnitudes, in the order the named traversal visits +# them. Every sweep of a heuristic cover reads `log|A_ij|` over the whole +# support; one shared vector replaces a `log` call per entry per sweep. A +# consumer handed one must walk the same traversal, in the same order, and take +# its `k`-th entry on the `k`-th callback. +function support_logs_sym(A::AbstractMatrix, ::Type{T}) where T + logs = T[] + foreach_support_sym(A) do i, j, v + push!(logs, log(T(v))) + end + return logs +end + +function support_logs(A::AbstractMatrix, ::Type{T}) where T + logs = T[] + foreach_support(A) do i, j, v + push!(logs, log(T(v))) + end + return logs +end + +# `log|A_ij|` for the `k`-th entry of a traversal, from the cache when the +# caller supplied one and from `v` when it did not. +@inline _entrylog(::Nothing, k::Int, v, ::Type{T}) where T = log(T(v)) +@inline _entrylog(logs::AbstractVector, k::Int, v, ::Type{T}) where T = logs[k]::T + # Apply the row/column balance convention independently to each support # component. Rounding the shift to a power of two preserves cover products # exactly, at the cost of balancing only within a factor of `sqrt(2)`. @@ -228,13 +256,16 @@ end # ∑_{i,j: A[i,j]≠0} (log(a[i]*a[j]) - log|A[i,j]|)² # Returns row support counts. The Sherman-Morrison approximation is exact on # complete support. -function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, A::AbstractMatrix) where T +function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, A::AbstractMatrix; + logs=nothing) where T ax = eachindex(a) axes(A) == (ax, ax) || throw(DimensionMismatch("`unconstrained_min!(ϕ, a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(string(axes(A))), axes(a)=$(string(axes(a)))")) loga = fill!(similar(a), zero(T)) nza = zeros(Int, ax) + k = Ref(0) foreach_support_sym(A) do i, j, v - lAij = log(T(v)) + k[] += 1 + lAij = _entrylog(logs, k[], v, T) loga[i] += lAij nza[i] += 1 if i != j @@ -253,7 +284,8 @@ function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, A::AbstractMatrix return nza end -function unconstrained_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix) +function unconstrained_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; + logs=nothing) T = float(promote_type(eltype(a), eltype(b))) axes(A, 1) == eachindex(a) || throw(DimensionMismatch("`unconstrained_min!(ϕ, a, b, A)` requires row indices of `A` to match `a`, got axes(A, 1)=$(string(axes(A, 1))), axes(a)=$(string(axes(a)))")) axes(A, 2) == eachindex(b) || throw(DimensionMismatch("`unconstrained_min!(ϕ, a, b, A)` requires column indices of `A` to match `b`, got axes(A, 2)=$(string(axes(A, 2))), axes(b)=$(string(axes(b)))")) @@ -261,8 +293,10 @@ function unconstrained_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A logb = fill!(similar(b, T), zero(T)) nza = zeros(Int, axes(A, 1)) nzb = zeros(Int, axes(A, 2)) + k = Ref(0) foreach_support(A) do i, j, v - lAij = log(T(v)) + k[] += 1 + lAij = _entrylog(logs, k[], v, T) loga[i] += lAij logb[j] += lAij nza[i] += 1 @@ -312,7 +346,7 @@ function _tighten_shrink(x, lr) return y end -function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; maxiter::Int=3) where T +function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; maxiter::Int=3, logs=nothing) where T ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("`tighten_cover!(a, A)` requires a square matrix `A`")) eachindex(a) == ax || throw(DimensionMismatch("indices of `a` must match the indexing of `A`")) @@ -326,8 +360,10 @@ function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; maxiter::Int=3) # in log-ratio (rather than a[i]*a[j]/v) keeps the comparison finite even when # the linear-space product would overflow for extreme dynamic range; a zero # scale gives lr = -Inf, marking the row uncoverable by any finite rescale. + k = Ref(0) foreach_support_sym(A) do i, j, v - lr = la[i] + la[j] - log(T(v)) + k[] += 1 + lr = la[i] + la[j] - _entrylog(logs, k[], v, T) lratio[i] = min(lratio[i], lr) i == j || (lratio[j] = min(lratio[j], lr)) end @@ -342,7 +378,7 @@ function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; maxiter::Int=3) return a end -function tighten_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=3) +function tighten_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=3, logs=nothing) T = float(promote_type(eltype(a), eltype(b))) eachindex(a) == axes(A, 1) || throw(DimensionMismatch("indices of a must match row-indexing of A")) eachindex(b) == axes(A, 2) || throw(DimensionMismatch("indices of b must match column-indexing of A")) @@ -358,8 +394,10 @@ function tighten_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; # finite even when the linear-space product would overflow for extreme # dynamic range; a zero scale gives lr = -Inf, marking the row or column # uncoverable by any finite rescale. + k = Ref(0) foreach_support(A) do i, j, v - lr = la[i] + lb[j] - log(T(v)) + k[] += 1 + lr = la[i] + lb[j] - _entrylog(logs, k[], v, T) lratioa[i] = min(lratioa[i], lr) lratiob[j] = min(lratiob[j], lr) end @@ -452,7 +490,7 @@ end # (the diagonal included, so no separate clamp step is needed). Requires a # start with strictly positive scale on every supported row (the geometric-mean # init from `unconstrained_min!` guarantees this). -function boost_feasible!(a::AbstractVector{T}, A::AbstractMatrix) where T +function boost_feasible!(a::AbstractVector{T}, A::AbstractMatrix; logs=nothing) where T IdxT = eltype(eachindex(a)) # `la` caches log.(a) and is updated alongside `a`, so deficits cost no log # calls; log(0) = -Inf on unsupported rows is never read (every entry's @@ -469,23 +507,27 @@ function boost_feasible!(a::AbstractVector{T}, A::AbstractMatrix) where T # of building it with `push!`. nviol = Ref(0) zmax = Ref(zero(T)) + k = Ref(0) foreach_support_sym(A) do i, j, v - z = log(T(v)) - la[i] - la[j] + k[] += 1 + z = _entrylog(logs, k[], v, T) - la[i] - la[j] if z > zero(T) nviol[] += 1 zmax[] = max(zmax[], z) end end entries = Vector{Tuple{IdxT,IdxT,T}}(undef, nviol[]) - k = Ref(0) + k[] = 0 + nfill = Ref(0) foreach_support_sym(A) do i, j, v - lv = log(T(v)) + k[] += 1 + lv = _entrylog(logs, k[], v, T) z = lv - la[i] - la[j] if z > zero(T) (iszero(a[i]) || iszero(a[j])) && throw(ArgumentError("boost_feasible! requires a start with positive scale on every supported row")) - k[] += 1 - entries[k[]] = (i, j, lv) + nfill[] += 1 + entries[nfill[]] = (i, j, lv) end end deficit((i, j, lv)) = lv - la[i] - la[j] @@ -503,7 +545,7 @@ end # for every entry visited by `foreach_support`. The # diagonal is treated as an ordinary entry. Requires a start with strictly # positive scale on every supported row of `a` and column of `b`. -function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) +function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; logs=nothing) T = float(promote_type(eltype(a), eltype(b))) IdxA, IdxB = eltype(eachindex(a)), eltype(eachindex(b)) # `la`/`lb` cache log.(a)/log.(b) and are updated alongside `a`/`b`; see @@ -516,23 +558,27 @@ function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix # repeated grow-and-copy of building it with `push!`. nviol = Ref(0) zmax = Ref(zero(T)) + k = Ref(0) foreach_support(A) do i, j, v - z = log(T(v)) - la[i] - lb[j] + k[] += 1 + z = _entrylog(logs, k[], v, T) - la[i] - lb[j] if z > zero(T) nviol[] += 1 zmax[] = max(zmax[], z) end end entries = Vector{Tuple{IdxA,IdxB,T}}(undef, nviol[]) - k = Ref(0) + k[] = 0 + nfill = Ref(0) foreach_support(A) do i, j, v - lv = log(T(v)) + k[] += 1 + lv = _entrylog(logs, k[], v, T) z = lv - la[i] - lb[j] if z > zero(T) (iszero(a[i]) || iszero(b[j])) && throw(ArgumentError("boost_feasible! requires a start with positive scale on every supported row/column")) - k[] += 1 - entries[k[]] = (i, j, lv) + nfill[] += 1 + entries[nfill[]] = (i, j, lv) end end deficit((i, j, lv)) = lv - la[i] - lb[j] @@ -661,12 +707,14 @@ end # Requires a start with strictly positive scale on every supported row and column. # The shift is covariant under an independent row/column rescaling `D_r*A*D_c`, # and is accumulated in the log domain for the reasons given in the symmetric method. -function inflate_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) +function inflate_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; logs=nothing) T = float(promote_type(eltype(a), eltype(b))) la, lb = map(log, a), map(log, b) tref = Ref(zero(T)) + k = Ref(0) foreach_support(A) do i, j, v - tref[] = max(tref[], (log(T(v)) - la[i] - lb[j]) / 2) + k[] += 1 + tref[] = max(tref[], (_entrylog(logs, k[], v, T) - la[i] - lb[j]) / 2) end t = tref[] # A supported row or column with zero scale gives la (or lb) = -Inf, hence t = +Inf. From 732cba09c7d30ba51fb913b33b272e1cb78699d9 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 24 Aug 2026 13:14:35 -0500 Subject: [PATCH 12/15] Match sparse symmetry partners by cursor Match `SparseMatrixCSC` symmetry partners with per-column cursors in O(nnz), preserving generic-check behavior for zeros and roundoff. Assisted-by: Claude Opus (claude-opus-5) --- src/sparse_support.jl | 36 ++++++++++++++++++++++++ src/support.jl | 10 +++++++ test/support.jl | 64 ++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 109 insertions(+), 1 deletion(-) diff --git a/src/sparse_support.jl b/src/sparse_support.jl index 4121d1d..2af33f0 100644 --- a/src/sparse_support.jl +++ b/src/sparse_support.jl @@ -30,6 +30,42 @@ function foreach_support_sym(f, A::SparseMatrixCSC) return nothing end +# Match symmetric partners in O(nnz) with one cursor per column. `tr[c]` points +# to the first unpaired entry in column `c`; absent entries are zeros. +function require_abs_symmetric(A::SparseMatrixCSC, fname) + ax = axes(A, 1) + axes(A, 2) == ax || + throw(DimensionMismatch("$fname requires a square matrix, got axes $(string(axes(A)))")) + rv, nzs = rowvals(A), nonzeros(A) + tr = [first(nzrange(A, j)) for j in axes(A, 2)] + for col in axes(A, 2) + for p in tr[col]:last(nzrange(A, col)) + v = abs(nzs[p]) + iszero(v) && continue + row = rv[p] + row == col && continue + row < col && _abs_asymmetry_error(fname, row, col, v, zero(v)) + off, stop = tr[row], last(nzrange(A, row)) + 1 + w = zero(v) + while off < stop + r2 = rv[off] + r2 > col && break + if r2 == col + w = abs(nzs[off]) + tr[row] = off + 1 + break + end + u = abs(nzs[off]) + iszero(u) || _abs_asymmetry_error(fname, r2, row, u, zero(u)) + off += 1 + tr[row] = off + end + _abs_symmetric(v, w) || _abs_asymmetry_error(fname, row, col, v, w) + end + end + return nothing +end + # Emitted pairs are canonical (row <= col) regardless of uplo: for uplo='L' # the stored (i, j) with i >= j is reported as (j, i). Complex `Hermitian` is # admitted alongside the real case because only `abs` of a stored value is ever diff --git a/src/support.jl b/src/support.jl index 30915aa..67d2505 100644 --- a/src/support.jl +++ b/src/support.jl @@ -97,6 +97,16 @@ function require_abs_symmetric(A::AbstractMatrix, fname) return nothing end +# Shared predicate and error for storage-specific symmetry checks. +_abs_symmetric(v, w) = (m = max(v, w); abs(v - w) <= ASYMMETRY_ULPS * eps(float(real(typeof(m)))) * m) + +@noinline function _abs_asymmetry_error(fname, i, j, v, w) + throw(ArgumentError(""" + $fname requires `abs.(A)` to be symmetric, but abs(A[$(string(i)),$(string(j))]) = $(string(v)) and \ + abs(A[$(string(j)),$(string(i))]) = $(string(w)). Wrap `A` in `Symmetric` (or `Hermitian`) to name the \ + triangle to read; that also skips this check.""")) +end + # Cache-blocked check for dense storage: the transposed reads of a column-major # sweep miss on every entry once `A` outgrows the cache, while a block of rows and # its transpose both fit. Only `i < j` needs testing, the diagonal being its own diff --git a/test/support.jl b/test/support.jl index af2a12b..f473e37 100644 --- a/test/support.jl +++ b/test/support.jl @@ -137,7 +137,7 @@ end @test symcover(Symmetric(M, :L)) isa AbstractVector @test symcover(Diagonal([1.0, 2.0])) isa AbstractVector - # Banded storage earns no exemption. A `Bidiagonal` reads one of its + # Banded storage is still checked. A `Bidiagonal` reads one of its # off-diagonals as a structural zero, so any nonzero band makes it asymmetric; # a `Tridiagonal` stores both bands and qualifies only when they agree. for A in (Bidiagonal([3.0, 2.0, 1.0], [6.0, 0.5], :U), @@ -154,6 +154,68 @@ end @test symcover(Bidiagonal([3.0, 2.0, 1.0], [0.0, 0.0], :U)) isa AbstractVector end +# Wrapper that uses the generic `AbstractMatrix` symmetry check. +struct Unstructured{M} <: AbstractMatrix{Float64} + A::M +end +Base.size(W::Unstructured) = size(W.A) +Base.getindex(W::Unstructured, i::Int, j::Int) = W.A[i, j] + +# The sparse and generic symmetry checks must agree. +@testset "the symmetry precondition on compressed columns" begin + outcome(A) = try + MatrixCovers.require_abs_symmetric(A, :f) + "accepted" + catch e + sprint(showerror, e) + end + + rng = StableRNG(31) + accepted = rejected = 0 + for _ in 1:300 + n = rand(rng, 4:12) + S = sprandn(rng, n, n, 0.25) + A = S + S' + defect = rand(rng, 1:4) + if !iszero(nnz(A)) + k = rand(rng, 1:nnz(A)) + if defect == 2 + nonzeros(A)[k] *= 1 + 1e-3 # value asymmetry + elseif defect == 3 + i = rowvals(A)[k] + j = findfirst(c -> k in nzrange(A, c), 1:n) + A[i, j] = 0 # structural asymmetry + dropzeros!(A) + elseif defect == 4 + nonzeros(A)[k] = 0.0 # explicit zero on one side + end + end + both = outcome(A) + both == "accepted" ? (accepted += 1) : (rejected += 1) + @test both == outcome(Unstructured(A)) + end + @test accepted > 20 && rejected > 20 + + # An absent partner reads as zero. + lone = sparse([2], [1], [3.0], 2, 2) + @test_throws "abs(A[2,1]) = 3.0 and abs(A[1,2]) = 0.0" MatrixCovers.require_abs_symmetric(lone, :f) + @test_throws "abs(A[1,2]) = 3.0 and abs(A[2,1]) = 0.0" MatrixCovers.require_abs_symmetric(sparse([1], [2], [3.0], 2, 2), :f) + + # Stored and implicit zeros behave alike. + ez = SparseMatrixCSC(2, 2, [1, 2, 2], [1], [0.0]) + @test outcome(ez) == "accepted" == outcome(Unstructured(ez)) + @test MatrixCovers.require_abs_symmetric(sparse([1, 2, 2], [1, 1, 2], [1.0, 0.0, 1.0], 2, 2), :f) === nothing + mixed = sparse([1, 2, 1, 2], [1, 1, 2, 2], [1.0, 0.0, 4.0, 1.0], 2, 2) + @test outcome(mixed) == outcome(Unstructured(mixed)) != "accepted" + + # Sparse storage uses the same roundoff allowance and wrapper exemptions. + B = sprandn(rng, 20, 20, 0.3); S = B + B' + d = exp.(randn(rng, 20)) + @test MatrixCovers.require_abs_symmetric((d .* S) .* d', :f) === nothing + @test symcover(S) isa AbstractVector + @test_throws DimensionMismatch MatrixCovers.require_abs_symmetric(sprandn(rng, 3, 4, 0.5), :f) +end + # Grouped support must preserve entries and symmetric full-grid multiplicity. @testset "grouped support reproduces the matrix" begin function regroup(S, groups_are_rows::Bool, sz) From 4936758d2fb1647ec6269a66b271dedea9764a4d Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 24 Aug 2026 15:23:05 -0500 Subject: [PATCH 13/15] Factor the LSQR preconditioner Right-precondition `Float64` LSQR with a reusable sparse Cholesky factor of the weighted normal matrix. A ridge handles singular support graphs. Assisted-by: Claude Opus (claude-opus-5) --- src/MatrixCovers.jl | 4 +- src/minimal_covers.jl | 182 +++++++++++++++++++---------------------- test/minimal_covers.jl | 16 +++- 3 files changed, 98 insertions(+), 104 deletions(-) diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index 50de63a..f3837ab 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -1,11 +1,11 @@ module MatrixCovers using LinearAlgebra: LinearAlgebra, Adjoint, Bidiagonal, Diagonal, Hermitian, - SymTridiagonal, Symmetric, Transpose, Tridiagonal, cholesky, + SymTridiagonal, Symmetric, Transpose, Tridiagonal, cholesky, cholesky!, dot, mul!, norm using PrecompileTools: PrecompileTools, @compile_workload using Random: Random, AbstractRNG, MersenneTwister -using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse +using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse, spzeros export AbsLog, AbsLinear export cover_objective, iscover diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 2ac2f6b..fc79a04 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -257,8 +257,8 @@ end # - `:woodbury` represents them as sparse `C + U*U'`, using conjugate gradients # while the condition estimate is small and sparse Cholesky otherwise. It is # restricted to nearly dense `Float64` problems. -# - `:lsqr` applies the weighted residual operator `M` matrix-free. In `Float64`, its -# right preconditioner includes violated rows once diagonal scaling is inadequate. +# - `:lsqr` applies the weighted residual operator `M` matrix-free, with sparse +# Cholesky right-preconditioning for `Float64`. # 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. @@ -272,9 +272,6 @@ const AUTO_LSQR_MAX_DENSITY = 1 // 4 # Condition estimate above which Woodbury uses sparse Cholesky instead of CG. const WOODBURY_CG_KAPPA = 1000 -# Condition estimate above which LSQR includes the weighted rows in its preconditioner. -const LSQR_PRECOND_KAPPA = 1000 - # Solve `(C + U*U')x = f` from a sparse factorization of `C` using the # Woodbury identity. `rhs` stores the combined `[f U]` solve. function _woodbury_solve!(x, F, U, f, rhs) @@ -761,6 +758,50 @@ function _warn_truncated(fname::Symbol, κs, stats, maxiter::Int) return nothing end +# Position of `S[i,j]` in `nonzeros(S)`; the entry must be stored. +function _nzindex(S::SparseMatrixCSC, i::Int, j::Int) + r = nzrange(S, j) + rv = rowvals(S) + k = searchsortedfirst(view(rv, r), i) + k <= length(r) && rv[r[k]] == i || + throw(ArgumentError("the preconditioner pattern is missing entry ($i, $j)")) + return r[k] +end + +# Unweighted normal-matrix pattern for the LSQR preconditioner. The ridge makes +# bipartite support components positive definite. `N == 0` disables it. +function _precond_pattern(::Type{T}, supp::EdgeList, v0, N::Int, mult) where {T} + N == 0 && return spzeros(T, 0, 0) + Mi, Mj, Mv = collect(1:N), collect(1:N), zeros(T, N) + for (p, q) in supp.edges + if p == q + Mv[p] += 4 * oneunit(T) + else + w = mult(p, q) * oneunit(T) + Mv[p] += w + Mv[q] += w + push!(Mi, p, q) + push!(Mj, q, p) + push!(Mv, w, w) + end + end + dmax = zero(T) + for p in 1:N + Mv[p] += v0[p]^2 + dmax = max(dmax, Mv[p]) + end + ρ = _precond_ridge(dmax) + for p in 1:N + Mv[p] += ρ + end + return sparse(Mi, Mj, Mv, N, N) +end + +_precond_pattern(::Type{T}, ::Grid, v0, N::Int, mult) where {T} = spzeros(T, 0, 0) + +# Scale-relative ridge for a positive-definite preconditioner. +_precond_ridge(dmax::T) where {T} = (dmax > 0 ? dmax : oneunit(T)) * sqrt(eps(T)) + # `AbsLog{2}` penalty continuation. Each stage freezes residual weights, solves # the weighted least-squares problem, and backtracks. `boost=true` applies a final # feasibility shift. The support layout selects the inner solver. @@ -792,7 +833,6 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; cv = zeros(T, ne + 1) # √weight · log|A_ij|, with a trailing 0 gauge target # Violated entries under the current frozen weights. vpat = _violation_pattern(supp) - vedges = Tuple{Int,Int}[] # the violated entries of the current solve (LSQR path only) vrow = Int[] # violated rows, grouped by column (Woodbury path only) vcnt = zeros(Int, use_woodbury ? N : 0) # violated off-diagonal entries per column vptr = zeros(Int, use_woodbury ? N + 1 : 0) @@ -820,33 +860,27 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; Ccolptr = zeros(Int, use_woodbury ? N + 1 : 0) Crowval = Int[] Cnzval = T[] - # Diagonal of the unweighted, gauge-augmented normal matrix. - dpart = zeros(T, use_lsqr ? N : 0) - if supp isa EdgeList && use_lsqr - for (p, q) in supp.edges - if p == q - dpart[p] += 4 * oneunit(T) - else - w = mult(p, q) * oneunit(T) - dpart[p] += w - dpart[q] += w - end - end - for p in 1:N - dpart[p] += v0[p]^2 - end - for p in 1:N - dpart[p] > 0 || (dpart[p] = oneunit(T)) + px = zeros(T, use_precond ? N : 0) # scale vector recovered from the LSQR variable + pg = zeros(T, use_precond ? N : 0) # `Rᵀ√W y` before the preconditioner is applied + mdiag = zeros(T, use_precond ? N : 0) # weighted degrees, the preconditioner's diagonal + # Right-preconditioner: a Cholesky factor of the weighted, gauge-augmented + # normal matrix. Its sparsity pattern is the support's and does not depend on + # the weights, so one symbolic analysis serves the whole continuation and each + # solve refactors in place. The gauge is a dense rank-1 term that would fill + # the factor, so it enters as its diagonal alone and LSQR absorbs the rest. + Msp = _precond_pattern(T, supp, v0, use_precond ? N : 0, mult) + # Positions of the entries each solve overwrites: the diagonal, and both + # copies of each off-diagonal support entry. + dpos = use_precond ? [_nzindex(Msp, p, p) for p in 1:N] : Int[] + epos = zeros(Int, use_precond ? 2 * ne : 0) + if use_precond + for (e, (p, q)) in enumerate(supp.edges) + p == q && continue + epos[2 * e - 1] = _nzindex(Msp, p, q) + epos[2 * e] = _nzindex(Msp, q, p) end end - mdiag = zeros(T, use_lsqr ? N : 0) # the violated rows' diagonal, per unit of κ−1 - Mi = Int[] # COO triplets of the preconditioner - Mj = Int[] - Mv = T[] - px = zeros(T, use_lsqr ? N : 0) # scale vector recovered from the LSQR variable - pg = zeros(T, use_lsqr ? N : 0) # `Rᵀ√W y` before the preconditioner is applied - # `K` of the diagonal preconditioner, which is κ-independent and so built once. - psqrt = use_lsqr ? sqrt.(dpart) : T[] + MF = use_precond ? cholesky(Symmetric(Msp)) : nothing rhs = zeros(T, use_woodbury ? N : 0, size(U, 2) + 1) dmin = use_woodbury ? minimum(sys.dfull) : oneunit(T) cgx = zeros(T, use_woodbury ? N : 0) @@ -915,9 +949,6 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; elseif use_lsqr edges = supp.edges cvals = supp.cvals - dκ = κ === nothing ? zero(T) : T(κ) - oneunit(T) - empty!(vedges) - fill!(mdiag, zero(T)) for (e, (p, q)) in enumerate(edges) c = cvals[e] viol = κ !== nothing && (x[p] + x[q] - c) < 0 @@ -925,80 +956,35 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; sw = sqrt(mult(p, q) * (viol ? T(κ) : oneunit(T))) ws[e] = sw cv[e] = sw * c - if viol && use_precond - push!(vedges, (p, q)) + end + g = ne + 1 # index of the appended gauge row + if use_precond + # Refill the preconditioner with the weights this solve freezes and + # refactor it in place: the pattern, and so the symbolic analysis, + # is the same on every solve. + nzv = nonzeros(Msp) + fill!(mdiag, zero(T)) + for (e, (p, q)) in enumerate(edges) + w = ws[e]^2 if p == q - mdiag[p] += 4 * oneunit(T) + mdiag[p] += 4 * w else - w = mult(p, q) * oneunit(T) mdiag[p] += w mdiag[q] += w + nzv[epos[2 * e - 1]] = w + nzv[epos[2 * e]] = w end end - end - g = ne + 1 # index of the appended gauge row - if use_precond - # Include violated rows once diagonal scaling is inadequate. - κest = oneunit(T) + dmax = zero(T) for p in 1:N - κest = max(κest, oneunit(T) + dκ * 2 * mdiag[p] / dpart[p]) + mdiag[p] += v0[p]^2 + dmax = max(dmax, mdiag[p]) end - if κest <= LSQR_PRECOND_KAPPA - # Diagonal `K` needs only elementwise scaling. - Dmul! = function (y, yv) - @. px = yv / psqrt - for (e, (p, q)) in enumerate(edges) - y[e] = ws[e] * (px[p] + px[q]) - end - y[g] = dot(v0, px) - return y - end - Dtmul! = function (z, y) - fill!(pg, zero(T)) - for (e, (p, q)) in enumerate(edges) - t = ws[e] * y[e] - pg[p] += t - pg[q] += t - end - @. pg += v0 * y[g] - @. z = pg / psqrt - return z - end - soly, it = _lsqr(Dmul!, Dtmul!, cv, psqrt .* x) - nlsqr[] += it - return soly ./ psqrt - end - empty!(Mi) - empty!(Mj) - empty!(Mv) + ρ = _precond_ridge(dmax) for p in 1:N - push!(Mi, p) - push!(Mj, p) - push!(Mv, dpart[p]) - end - for (p, q) in vedges - if p == q - push!(Mi, p) - push!(Mj, p) - push!(Mv, 4 * dκ) - else - w = mult(p, q) * dκ - push!(Mi, p) - push!(Mj, p) - push!(Mv, w) - push!(Mi, q) - push!(Mj, q) - push!(Mv, w) - push!(Mi, p) - push!(Mj, q) - push!(Mv, w) - push!(Mi, q) - push!(Mj, p) - push!(Mv, w) - end + nzv[dpos[p]] = mdiag[p] + ρ end - Msp = sparse(Mi, Mj, Mv, N, N) - MF = cholesky(Symmetric(Msp)) + cholesky!(MF, Symmetric(Msp)) Kc = MF.PtL Uc = MF.UP # CHOLMOD factor-component solves allocate their result. diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 1309927..00a480a 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -380,8 +380,7 @@ end @test soft_symcover_min(AbsLog{2}(), A; linsolve=:lsqr) != symcover(A) end -# The Float64 LSQR preconditioner includes κ-weighted rows; other types use the -# plain matrix-free iteration. +# Bound iterations for Cholesky-preconditioned `Float64` LSQR. @testset "MMC :lsqr iteration count is bounded across the continuation" begin rng = StableRNG(5) A = (X = exp.(randn(rng, 120, 120)); (X .+ X') ./ 2) @@ -389,14 +388,23 @@ end al, sl = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) @test al ≈ ad rtol=1e-6 # Bound the preconditioned iteration count with margin. - @test sl.lsqriters <= 60 * sl.nsolves + @test sl.lsqriters <= 12 * sl.nsolves G = exp.(randn(rng, 120, 90)) gd, hd, _ = MatrixCovers._cover_min_abslog2(G; linsolve=:dense) gl, hl, tl = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr) @test gl .* hl' ≈ gd .* hd' rtol=1e-6 # Apply the same iteration bound to asymmetric problems. - @test tl.lsqriters <= 60 * tl.nsolves + @test tl.lsqriters <= 12 * tl.nsolves + + # The ridge handles bipartite support. + rngb = StableRNG(11) + Tsp = sparse(Matrix(SymTridiagonal(exp.(randn(rngb, 40)), exp.(randn(rngb, 39))))) + at, st = MatrixCovers._symcover_min_abslog2(Tsp; linsolve=:lsqr) + atd, _ = MatrixCovers._symcover_min_abslog2(Matrix(Tsp); linsolve=:dense) + @test at ≈ atd rtol=1e-6 + @test iscover(at, at, Tsp) + @test st.lsqriters <= 12 * st.nsolves # A working type CHOLMOD cannot factor keeps the plain matrix-free iteration. A32 = Float32.([4.0 1.0 0.5; 1.0 3.0 1.0; 0.5 1.0 2.5]) From db8c0631f0f44ae16537edf3ad10abae2f9f0d8b Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 25 Aug 2026 03:25:03 -0500 Subject: [PATCH 14/15] Flatten heuristic support sweeps, vectorize logs Flatten support into compact row, column, and log-magnitude arrays, using a SIMD logarithm approximation. Preserve traversal order and certify returned covers in linear arithmetic. Assisted-by: Claude Fable 5 (claude-fable-5) --- docs/src/index.md | 8 +- src/MatrixCovers.jl | 3 +- src/dense_heuristic.jl | 70 ++++---- src/fastlog.jl | 31 ++++ src/heuristic_covers.jl | 342 ++++++++++++++++++++++++++++++++-------- src/sparse_support.jl | 6 + 6 files changed, 353 insertions(+), 107 deletions(-) create mode 100644 src/fastlog.jl diff --git a/docs/src/index.md b/docs/src/index.md index 46bb182..1d7647c 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -332,7 +332,9 @@ julia> round.(a; digits=6) 1000.0 2.0 -julia> mag = sum(abs(bi / ai) for (bi, ai) in zip(b, a)) +julia> mag = sum(abs(bi / ai) for (bi, ai) in zip(b, a)); + +julia> round(mag; digits=6) 4.5 ``` @@ -341,7 +343,7 @@ Here `mag` is within a factor of 1.5 of the scaled solution norm: ```jldoctest roundoff julia> x = A \ b; -julia> sum(abs.(x .* a)) +julia> round(sum(abs.(x .* a)); digits=6) 3.0 ``` @@ -354,7 +356,7 @@ julia> Ad, bd = d .* A .* d', d .* b; julia> ad = symcover(Ad); -julia> sum(abs(bi / ai) for (bi, ai) in zip(bd, ad)) +julia> round(sum(abs(bi / ai) for (bi, ai) in zip(bd, ad)); digits=6) 4.5 ``` diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index f3837ab..beb8848 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -5,7 +5,7 @@ using LinearAlgebra: LinearAlgebra, Adjoint, Bidiagonal, Diagonal, Hermitian, dot, mul!, norm using PrecompileTools: PrecompileTools, @compile_workload using Random: Random, AbstractRNG, MersenneTwister -using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse, spzeros +using SparseArrays: SparseArrays, SparseMatrixCSC, nnz, nonzeros, nzrange, rowvals, sparse, spzeros export AbsLog, AbsLinear export cover_objective, iscover @@ -24,6 +24,7 @@ end include("penalties.jl") include("support.jl") include("iscover.jl") +include("fastlog.jl") include("heuristic_covers.jl") include("dense_heuristic.jl") # full-grid kernels for the heuristic covers include("gram_covers.jl") # symmetric covers of A'*W*A from an asymmetric cover of A diff --git a/src/dense_heuristic.jl b/src/dense_heuristic.jl index 3b315c7..2a092a9 100644 --- a/src/dense_heuristic.jl +++ b/src/dense_heuristic.jl @@ -1,13 +1,7 @@ # Dense-grid kernels for the heuristic covers. # -# `symcover!` and `cover!` sweep the support five or six times over, and every -# sweep of the callback-driven implementations recomputes `log(abs(A[i,j]))`. -# When the support is the full grid those logarithms dominate the run time, so -# these kernels evaluate them once into a log-magnitude grid and run the -# remaining sweeps over that grid. Every sum accumulates in the traversal order -# of `foreach_support`/`foreach_support_sym`, and each residual is formed from -# the same three numbers, so the covers agree with the callback-driven path bit -# for bit. +# Full-grid storage needs only a log-magnitude array because its indices are +# implicit. Its accumulation order matches `FlatSupport`. # Use the grid only when the support traversal is already dense and its # allocation is worthwhile. @@ -26,30 +20,30 @@ _use_dense_grid(A::AbstractMatrix, ::Type{T}) where {T} = # occupies `_trioff(j)+1 : _trioff(j)+j`. _trioff(j::Int) = (j * (j - 1)) >> 1 -# `Lp[_trioff(j)+i] = log(abs(A[i,j]))` for `i <= j`, `-Inf` where the entry is -# zero, alongside the per-row log sums and support counts of -# `unconstrained_min!`. `s[r]` receives the partners of row `r` in increasing -# partner order — the columns `1:r-1` of row `r`'s own column first, then the -# diagonal, then the later columns — which is the order `foreach_support_sym` -# feeds them in. +# Pack the upper triangle as log magnitudes, with `-Inf` for zeros, and compute +# the row sums and support counts used by `unconstrained_min!`. function _tri_logabs!(Lp::Vector{T}, s::Vector{T}, cnt::Vector{Int}, A::AbstractMatrix) where {T} ax = axes(A, 1) or = first(ax) - 1 n = length(ax) + for jp in 1:n + j = jp + or + o = _trioff(jp) + for ip in 1:jp + Lp[o+ip] = abs(A[ip+or, j]) + end + end + _fastlog!(Lp) fill!(s, zero(T)) fill!(cnt, 0) + ninf = T(-Inf) for jp in 1:n - j = jp + or o = _trioff(jp) sj = zero(T) cj = 0 for ip in 1:jp-1 - v = abs(A[ip+or, j]) - if iszero(v) - Lp[o+ip] = T(-Inf) - else - l = log(T(v)) - Lp[o+ip] = l + l = Lp[o+ip] + if l != ninf s[ip] += l cnt[ip] += 1 sj += l @@ -58,12 +52,8 @@ function _tri_logabs!(Lp::Vector{T}, s::Vector{T}, cnt::Vector{Int}, A::Abstract end s[jp] += sj cnt[jp] += cj - v = abs(A[j, j]) - if iszero(v) - Lp[o+jp] = T(-Inf) - else - l = log(T(v)) - Lp[o+jp] = l + l = Lp[o+jp] + if l != ninf s[jp] += l cnt[jp] += 1 end @@ -71,26 +61,28 @@ function _tri_logabs!(Lp::Vector{T}, s::Vector{T}, cnt::Vector{Int}, A::Abstract return Lp end -# `L[i,j] = log(abs(A[i,j]))`, `-Inf` where the entry is zero, alongside the -# per-row and per-column log sums and support counts of `unconstrained_min!`. +# Fill a log-magnitude grid and the row and column summaries. function _grid_logabs!(L::Matrix{T}, sa::Vector{T}, sb::Vector{T}, na::Vector{Int}, nb::Vector{Int}, A::AbstractMatrix) where {T} or = first(axes(A, 1)) - 1 oc = first(axes(A, 2)) - 1 m, n = size(A) + for jp in 1:n + j = jp + oc + for ip in 1:m + L[ip, jp] = abs(A[ip+or, j]) + end + end + _fastlog!(L) fill!(sa, zero(T)) fill!(na, 0) + ninf = T(-Inf) for jp in 1:n - j = jp + oc sj = zero(T) cj = 0 for ip in 1:m - v = abs(A[ip+or, j]) - if iszero(v) - L[ip, jp] = T(-Inf) - else - l = log(T(v)) - L[ip, jp] = l + l = L[ip, jp] + if l != ninf sa[ip] += l na[ip] += 1 sj += l @@ -114,11 +106,7 @@ end # Use compact boost-list indices when the dimensions fit. _grid_label(m::Int, n::Int) = max(m, n) <= typemax(Int32) ? Int32 : Int -# The violated entries of a packed upper triangle, in the traversal order of -# `foreach_support_sym`. Half the entries of a fresh unconstrained start -# violate, so a branch on the test would mispredict on half the grid; the sweep -# selects branchlessly instead, writing every entry to the slot after the last -# one kept and advancing only on a violation. That needs one slot of slack. +# Select violated upper-triangle entries in traversal order without branching. function _tri_violated(Lp::Vector{T}, lα::Vector{T}, n::Int, nviol::Int, ::Type{IT}) where {T,IT} entries = Vector{Tuple{IT,IT,T}}(undef, nviol + 1) diff --git a/src/fastlog.jl b/src/fastlog.jl new file mode 100644 index 0000000..02ded3d --- /dev/null +++ b/src/fastlog.jl @@ -0,0 +1,31 @@ +# SIMD-friendly logarithm using exponent reduction and an `atanh` series. +# Returned covers are certified separately. +# +# Callers pass nonnegative values. Zero, infinity, and NaN match `Base.log`. +@inline function _fastlog(x::Float64) + # Scale subnormals into the normal range so exponent extraction sees them. + sub = x < floatmin(Float64) + xs = ifelse(sub, x * 0x1p54, x) + bits = reinterpret(UInt64, xs) + e = Float64(Int64(bits >> 52) - 1023) - ifelse(sub, 54.0, 0.0) + m = reinterpret(Float64, (bits & 0x000f_ffff_ffff_ffff) | 0x3ff0_0000_0000_0000) + # Reduce the mantissa from [1, 2) to [√2/2, √2), centering log(m) on zero. + big = m > 1.4142135623730951 + m = ifelse(big, 0.5 * m, m) + e = ifelse(big, e + 1.0, e) + f = m - 1.0 + s = f / (2.0 + f) + # Six terms bound the truncation error below 6e-11. + p = @evalpoly(s * s, 2.0, 2 / 3, 2 / 5, 2 / 7, 2 / 9, 2 / 11) + r = e * 0.6931471805599453 + s * p + return ifelse(iszero(x), -Inf, ifelse(x < Inf, r, x)) +end +_fastlog(x::Float32) = Float32(_fastlog(Float64(x))) +_fastlog(x::Real) = log(x) # element types the kernel does not cover + +function _fastlog!(x::AbstractArray) + @simd for k in eachindex(x) + x[k] = _fastlog(x[k]) + end + return x +end diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index 46a9abe..125da45 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -72,10 +72,10 @@ function _symcover!(a::AbstractVector, A::AbstractMatrix; maxiter::Int=3) if _use_dense_grid(A, T) _symcover_dense!(a, A, T, maxiter) else - logs = support_logs_sym(A, eltype(a)) - unconstrained_min!(AbsLog{2}(), a, A; logs) - boost_feasible!(a, A; logs) - tighten_cover!(a, A; maxiter, logs) + sup = flat_support_sym(A, T) + unconstrained_min!(AbsLog{2}(), a, sup) + boost_feasible!(a, sup) + tighten_cover!(a, sup; maxiter) end # Certify against `A` after log-domain tightening. return _certify_cover!(a, A, :symcover) @@ -157,13 +157,13 @@ function _cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxite if _use_dense_grid(A, T) _cover_dense!(a, b, A, T, maxiter) else - logs = support_logs(A, T) - unconstrained_min!(AbsLog{2}(), a, b, A; logs) - boost_feasible!(a, b, A; logs) - tighten_cover!(a, b, A; maxiter, logs) + sup = flat_support(A, T) + unconstrained_min!(AbsLog{2}(), a, b, sup) + boost_feasible!(a, b, sup) + tighten_cover!(a, b, sup; maxiter) # Apply the package's balance convention, then restore coverage lost to rounding. _balance_cover!(a, b, A) - inflate_feasible!(a, b, A; logs) + inflate_feasible!(a, b, sup) end return _certify_cover!(a, b, A, :cover) end @@ -181,31 +181,72 @@ end # ============================================================ # Internal helpers # ============================================================ -# Logarithms of the stored magnitudes, in the order the named traversal visits -# them. Every sweep of a heuristic cover reads `log|A_ij|` over the whole -# support; one shared vector replaces a `log` call per entry per sweep. A -# consumer handed one must walk the same traversal, in the same order, and take -# its `k`-th entry on the `k`-th callback. -function support_logs_sym(A::AbstractMatrix, ::Type{T}) where T - logs = T[] +# Matrix support flattened in traversal order as row, column, and `log|A_ij|` +# arrays. +struct FlatSupport{Ti<:Integer,Tj<:Integer,T} + is::Vector{Ti} + js::Vector{Tj} + lv::Vector{T} +end + +# Use `Int32` when it contains the axis; otherwise preserve the axis index type. +function _flat_index_type(ax) + I = eltype(ax) + I <: Integer || return I + isempty(ax) && return Int32 + return (typemin(Int32) <= first(ax) && last(ax) <= typemax(Int32)) ? Int32 : I +end + +# Storage-specific upper bounds for `sizehint!`; zero means unknown. +_support_sizehint(::AbstractMatrix) = 0 +_support_sizehint_sym(::AbstractMatrix) = 0 + +# The outer methods select concrete index types for the traversal. +flat_support_sym(A::AbstractMatrix, ::Type{T}) where T = + _flat_support_sym(A, T, _flat_index_type(axes(A, 1))) + +function _flat_support_sym(A::AbstractMatrix, ::Type{T}, ::Type{Ti}) where {T,Ti} + is, js, lv = Ti[], Ti[], T[] + hint = _support_sizehint_sym(A) + if hint > 0 + sizehint!(is, hint); sizehint!(js, hint); sizehint!(lv, hint) + end foreach_support_sym(A) do i, j, v - push!(logs, log(T(v))) + push!(is, i); push!(js, j); push!(lv, T(v)) end - return logs + _fastlog!(lv) # one vectorized pass over the collected magnitudes + return FlatSupport(is, js, lv) end -function support_logs(A::AbstractMatrix, ::Type{T}) where T - logs = T[] +flat_support(A::AbstractMatrix, ::Type{T}) where T = + _flat_support(A, T, _flat_index_type(axes(A, 1)), _flat_index_type(axes(A, 2))) + +function _flat_support(A::AbstractMatrix, ::Type{T}, ::Type{Ti}, ::Type{Tj}) where {T,Ti,Tj} + is, js, lv = Ti[], Tj[], T[] + hint = _support_sizehint(A) + if hint > 0 + sizehint!(is, hint); sizehint!(js, hint); sizehint!(lv, hint) + end foreach_support(A) do i, j, v - push!(logs, log(T(v))) + push!(is, i); push!(js, j); push!(lv, T(v)) end - return logs + _fastlog!(lv) # one vectorized pass over the collected magnitudes + return FlatSupport(is, js, lv) end -# `log|A_ij|` for the `k`-th entry of a traversal, from the cache when the -# caller supplied one and from `v` when it did not. -@inline _entrylog(::Nothing, k::Int, v, ::Type{T}) where T = log(T(v)) -@inline _entrylog(logs::AbstractVector, k::Int, v, ::Type{T}) where T = logs[k]::T +# Select violated entries in traversal order without branching. +function _flat_violated(sup::FlatSupport{Ti,Tj,T}, la, lb, nviol::Int) where {Ti,Tj,T} + is, js, lv = sup.is, sup.js, sup.lv + entries = Vector{Tuple{Ti,Tj,T}}(undef, nviol + 1) + k = 1 + for p in eachindex(is, js, lv) + i, j, lvp = is[p], js[p], lv[p] + entries[k] = (i, j, lvp) + k += ifelse(lvp - la[i] - lb[j] > zero(T), 1, 0) + end + resize!(entries, nviol) + return entries +end # Apply the row/column balance convention independently to each support # component. Rounding the shift to a power of two preserves cover products @@ -256,16 +297,13 @@ end # ∑_{i,j: A[i,j]≠0} (log(a[i]*a[j]) - log|A[i,j]|)² # Returns row support counts. The Sherman-Morrison approximation is exact on # complete support. -function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, A::AbstractMatrix; - logs=nothing) where T +function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, A::AbstractMatrix) where T ax = eachindex(a) axes(A) == (ax, ax) || throw(DimensionMismatch("`unconstrained_min!(ϕ, a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(string(axes(A))), axes(a)=$(string(axes(a)))")) loga = fill!(similar(a), zero(T)) nza = zeros(Int, ax) - k = Ref(0) foreach_support_sym(A) do i, j, v - k[] += 1 - lAij = _entrylog(logs, k[], v, T) + lAij = log(T(v)) loga[i] += lAij nza[i] += 1 if i != j @@ -284,8 +322,33 @@ function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, A::AbstractMatrix return nza end -function unconstrained_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; - logs=nothing) +# The symmetric objective over a flattened support: `sup` must have been built +# by `flat_support_sym` over a matrix whose axes match `eachindex(a)`. +function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, sup::FlatSupport) where T + is, js, lv = sup.is, sup.js, sup.lv + loga = fill!(similar(a), zero(T)) + nza = zeros(Int, eachindex(a)) + for k in eachindex(is, js, lv) + i, j, lAij = is[k], js[k], lv[k] + loga[i] += lAij + nza[i] += 1 + if i != j + loga[j] += lAij + nza[j] += 1 + end + end + nztotal = sum(nza) + halfmu = iszero(nztotal) ? zero(T) : sum(loga) / (2 * nztotal) + for i in eachindex(a) + # exp can underflow for extreme dynamic range; a zero scale on a + # supported row would make the boost's log-deficits infinite, so + # clamp to the smallest normal positive value. + a[i] = iszero(nza[i]) ? zero(T) : max(exp(loga[i] / nza[i] - halfmu), floatmin(T)) + end + return nza +end + +function unconstrained_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix) T = float(promote_type(eltype(a), eltype(b))) axes(A, 1) == eachindex(a) || throw(DimensionMismatch("`unconstrained_min!(ϕ, a, b, A)` requires row indices of `A` to match `a`, got axes(A, 1)=$(string(axes(A, 1))), axes(a)=$(string(axes(a)))")) axes(A, 2) == eachindex(b) || throw(DimensionMismatch("`unconstrained_min!(ϕ, a, b, A)` requires column indices of `A` to match `b`, got axes(A, 2)=$(string(axes(A, 2))), axes(b)=$(string(axes(b)))")) @@ -293,10 +356,8 @@ function unconstrained_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A logb = fill!(similar(b, T), zero(T)) nza = zeros(Int, axes(A, 1)) nzb = zeros(Int, axes(A, 2)) - k = Ref(0) foreach_support(A) do i, j, v - k[] += 1 - lAij = _entrylog(logs, k[], v, T) + lAij = log(T(v)) loga[i] += lAij logb[j] += lAij nza[i] += 1 @@ -318,6 +379,39 @@ function unconstrained_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A return nza, nzb end +# The asymmetric objective over a flattened support: `sup` must have been +# built by `flat_support` over a matrix whose row and column axes match +# `eachindex(a)` and `eachindex(b)`. +function unconstrained_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, sup::FlatSupport) + T = float(promote_type(eltype(a), eltype(b))) + is, js, lv = sup.is, sup.js, sup.lv + loga = fill!(similar(a, T), zero(T)) + logb = fill!(similar(b, T), zero(T)) + nza = zeros(Int, eachindex(a)) + nzb = zeros(Int, eachindex(b)) + for k in eachindex(is, js, lv) + i, j, lAij = is[k], js[k], lv[k] + loga[i] += lAij + logb[j] += lAij + nza[i] += 1 + nzb[j] += 1 + end + # Each stored entry contributes lAij to loga exactly once and increments + # nza exactly once, so these sums equal the per-entry running totals. + nztotal = sum(nza) + halfmu = iszero(nztotal) ? zero(T) : sum(loga) / (2 * nztotal) + for i in eachindex(a) + # exp can underflow for extreme dynamic range; a zero scale on a + # supported row would make the boost's log-deficits infinite, so + # clamp to the smallest normal positive value. + a[i] = iszero(nza[i]) ? zero(T) : max(exp(loga[i] / nza[i] - halfmu), floatmin(T)) + end + for j in eachindex(b) + b[j] = iszero(nzb[j]) ? zero(T) : max(exp(logb[j] / nzb[j] - halfmu), floatmin(T)) + end + return nza, nzb +end + # Feasible cover starting from the diagonal alone, resolved by # `boost_feasible_seq!`. Unlike `boost_feasible!`, a zero entry of `a` going # into that call means "not yet resolved", not "permanently unsupported" — @@ -346,7 +440,7 @@ function _tighten_shrink(x, lr) return y end -function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; maxiter::Int=3, logs=nothing) where T +function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; maxiter::Int=3) where T ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("`tighten_cover!(a, A)` requires a square matrix `A`")) eachindex(a) == ax || throw(DimensionMismatch("indices of `a` must match the indexing of `A`")) @@ -360,10 +454,8 @@ function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; maxiter::Int=3, # in log-ratio (rather than a[i]*a[j]/v) keeps the comparison finite even when # the linear-space product would overflow for extreme dynamic range; a zero # scale gives lr = -Inf, marking the row uncoverable by any finite rescale. - k = Ref(0) foreach_support_sym(A) do i, j, v - k[] += 1 - lr = la[i] + la[j] - _entrylog(logs, k[], v, T) + lr = la[i] + la[j] - log(T(v)) lratio[i] = min(lratio[i], lr) i == j || (lratio[j] = min(lratio[j], lr)) end @@ -378,7 +470,33 @@ function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; maxiter::Int=3, return a end -function tighten_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=3, logs=nothing) +# Symmetric tightening over a flattened support; the log-ratio convention is +# that of the matrix method. The `ifelse` minimum rejects NaN (an infinite +# entry against a zero scale), leaving such rows at the +Inf no-op, and lets +# the loop run branch-free; a diagonal entry updates its row twice, which the +# minimum absorbs. +function tighten_cover!(a::AbstractVector{T}, sup::FlatSupport; maxiter::Int=3) where T + is, js, lv = sup.is, sup.js, sup.lv + lratio = similar(a) + la = similar(a) + for _ in 1:maxiter + map!(log, la, a) # log(0) = -Inf marks zero scales; see the matrix method + fill!(lratio, T(Inf)) + for k in eachindex(is, js, lv) + i, j = is[k], js[k] + lr = la[i] + la[j] - lv[k] + lratio[i] = ifelse(lr < lratio[i], lr, lratio[i]) + lratio[j] = ifelse(lr < lratio[j], lr, lratio[j]) + end + for i in eachindex(a) + lr = lratio[i] + isinf(lr) || (a[i] = _tighten_shrink(a[i], lr)) + end + end + return a +end + +function tighten_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=3) T = float(promote_type(eltype(a), eltype(b))) eachindex(a) == axes(A, 1) || throw(DimensionMismatch("indices of a must match row-indexing of A")) eachindex(b) == axes(A, 2) || throw(DimensionMismatch("indices of b must match column-indexing of A")) @@ -394,10 +512,8 @@ function tighten_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; # finite even when the linear-space product would overflow for extreme # dynamic range; a zero scale gives lr = -Inf, marking the row or column # uncoverable by any finite rescale. - k = Ref(0) foreach_support(A) do i, j, v - k[] += 1 - lr = la[i] + lb[j] - _entrylog(logs, k[], v, T) + lr = la[i] + lb[j] - log(T(v)) lratioa[i] = min(lratioa[i], lr) lratiob[j] = min(lratiob[j], lr) end @@ -417,6 +533,37 @@ function tighten_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; return a, b end +# Asymmetric tightening over a flattened support; see the symmetric flat +# method for the `ifelse`-minimum convention. +function tighten_cover!(a::AbstractVector, b::AbstractVector, sup::FlatSupport; maxiter::Int=3) + T = float(promote_type(eltype(a), eltype(b))) + is, js, lv = sup.is, sup.js, sup.lv + lratioa = fill(T(Inf), eachindex(a)) + lratiob = fill(T(Inf), eachindex(b)) + la, lb = similar(a, T), similar(b, T) + for _ in 1:maxiter + map!(log, la, a) # log(0) = -Inf marks zero scales; see the matrix method + map!(log, lb, b) + fill!(lratioa, T(Inf)) + fill!(lratiob, T(Inf)) + for k in eachindex(is, js, lv) + i, j = is[k], js[k] + lr = la[i] + lb[j] - lv[k] + lratioa[i] = ifelse(lr < lratioa[i], lr, lratioa[i]) + lratiob[j] = ifelse(lr < lratiob[j], lr, lratiob[j]) + end + for i in eachindex(a) + lr = lratioa[i] + isinf(lr) || (a[i] = _tighten_shrink(a[i], lr)) + end + for j in eachindex(b) + lr = lratiob[j] + isinf(lr) || (b[j] = _tighten_shrink(b[j], lr)) + end + end + return a, b +end + # Adjoint/Transpose wrappers for tighten_cover!. function tighten_cover!(a::AbstractVector, b::AbstractVector, A::Adjoint; kwargs...) tighten_cover!(b, a, parent(A); kwargs...) @@ -490,7 +637,7 @@ end # (the diagonal included, so no separate clamp step is needed). Requires a # start with strictly positive scale on every supported row (the geometric-mean # init from `unconstrained_min!` guarantees this). -function boost_feasible!(a::AbstractVector{T}, A::AbstractMatrix; logs=nothing) where T +function boost_feasible!(a::AbstractVector{T}, A::AbstractMatrix) where T IdxT = eltype(eachindex(a)) # `la` caches log.(a) and is updated alongside `a`, so deficits cost no log # calls; log(0) = -Inf on unsupported rows is never read (every entry's @@ -507,21 +654,17 @@ function boost_feasible!(a::AbstractVector{T}, A::AbstractMatrix; logs=nothing) # of building it with `push!`. nviol = Ref(0) zmax = Ref(zero(T)) - k = Ref(0) foreach_support_sym(A) do i, j, v - k[] += 1 - z = _entrylog(logs, k[], v, T) - la[i] - la[j] + z = log(T(v)) - la[i] - la[j] if z > zero(T) nviol[] += 1 zmax[] = max(zmax[], z) end end entries = Vector{Tuple{IdxT,IdxT,T}}(undef, nviol[]) - k[] = 0 nfill = Ref(0) foreach_support_sym(A) do i, j, v - k[] += 1 - lv = _entrylog(logs, k[], v, T) + lv = log(T(v)) z = lv - la[i] - la[j] if z > zero(T) (iszero(a[i]) || iszero(a[j])) && @@ -540,12 +683,42 @@ function boost_feasible!(a::AbstractVector{T}, A::AbstractMatrix; logs=nothing) return a end +# Symmetric boost over a flattened support. As in the matrix method, only +# entries already violated at the start are stored; the count pass runs +# branchlessly (about half a fresh start's entries violate, so a data-dependent +# branch would mispredict constantly), and `_flat_violated` selects them the +# same way. A zero scale on a supported row makes some deficit +Inf, which the +# `isfinite` check below turns into the matrix method's error. +function boost_feasible!(a::AbstractVector{T}, sup::FlatSupport) where T + is, js, lv = sup.is, sup.js, sup.lv + # `la` caches log.(a) and is updated alongside `a`; see the matrix method. + la = map(log, a) + nviol = 0 + zmax = zero(T) + for k in eachindex(is, js, lv) + z = lv[k] - la[is[k]] - la[js[k]] + nviol += ifelse(z > zero(T), 1, 0) + zmax = ifelse(z > zmax, z, zmax) + end + isfinite(zmax) || + throw(ArgumentError("boost_feasible! requires a start with positive scale on every supported row")) + entries = _flat_violated(sup, la, la, nviol) + deficit((i, j, lvk)) = lvk - la[i] - la[j] + function apply!((i, j, lvk), z) + h = z / 2 + la[i] += h; a[i] = exp(la[i]) + i == j || (la[j] += h; a[j] = exp(la[j])) + end + bucket_boost!(deficit, apply!, entries, T, zmax) + return a +end + # Asymmetric feasibility boost: scale `a`, `b` in place so that # `a[i]*b[j] >= |A[i,j]|`, up to the round-off of the log-domain updates, # for every entry visited by `foreach_support`. The # diagonal is treated as an ordinary entry. Requires a start with strictly # positive scale on every supported row of `a` and column of `b`. -function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; logs=nothing) +function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) T = float(promote_type(eltype(a), eltype(b))) IdxA, IdxB = eltype(eachindex(a)), eltype(eachindex(b)) # `la`/`lb` cache log.(a)/log.(b) and are updated alongside `a`/`b`; see @@ -558,21 +731,17 @@ function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix # repeated grow-and-copy of building it with `push!`. nviol = Ref(0) zmax = Ref(zero(T)) - k = Ref(0) foreach_support(A) do i, j, v - k[] += 1 - z = _entrylog(logs, k[], v, T) - la[i] - lb[j] + z = log(T(v)) - la[i] - lb[j] if z > zero(T) nviol[] += 1 zmax[] = max(zmax[], z) end end entries = Vector{Tuple{IdxA,IdxB,T}}(undef, nviol[]) - k[] = 0 nfill = Ref(0) foreach_support(A) do i, j, v - k[] += 1 - lv = _entrylog(logs, k[], v, T) + lv = log(T(v)) z = lv - la[i] - lb[j] if z > zero(T) (iszero(a[i]) || iszero(b[j])) && @@ -591,6 +760,33 @@ function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix return a, b end +# Asymmetric boost over a flattened support; see the symmetric flat method. +function boost_feasible!(a::AbstractVector, b::AbstractVector, sup::FlatSupport) + T = float(promote_type(eltype(a), eltype(b))) + is, js, lv = sup.is, sup.js, sup.lv + # `la`/`lb` cache log.(a)/log.(b) and are updated alongside `a`/`b`; see + # the matrix methods. + la, lb = map(log, a), map(log, b) + nviol = 0 + zmax = zero(T) + for k in eachindex(is, js, lv) + z = lv[k] - la[is[k]] - lb[js[k]] + nviol += ifelse(z > zero(T), 1, 0) + zmax = ifelse(z > zmax, z, zmax) + end + isfinite(zmax) || + throw(ArgumentError("boost_feasible! requires a start with positive scale on every supported row/column")) + entries = _flat_violated(sup, la, lb, nviol) + deficit((i, j, lvk)) = lvk - la[i] - lb[j] + function apply!((i, j, lvk), z) + h = z / 2 + la[i] += h; a[i] = exp(la[i]) + lb[j] += h; b[j] = exp(lb[j]) + end + bucket_boost!(deficit, apply!, entries, T, zmax) + return a, b +end + # Sequential nearest-neighbor feasibility propagation by diagonal offset. function boost_feasible_seq!(a::AbstractVector{T}, A::AbstractMatrix) where T ax = eachindex(a) @@ -707,14 +903,12 @@ end # Requires a start with strictly positive scale on every supported row and column. # The shift is covariant under an independent row/column rescaling `D_r*A*D_c`, # and is accumulated in the log domain for the reasons given in the symmetric method. -function inflate_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; logs=nothing) +function inflate_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) T = float(promote_type(eltype(a), eltype(b))) la, lb = map(log, a), map(log, b) tref = Ref(zero(T)) - k = Ref(0) foreach_support(A) do i, j, v - k[] += 1 - tref[] = max(tref[], (_entrylog(logs, k[], v, T) - la[i] - lb[j]) / 2) + tref[] = max(tref[], (log(T(v)) - la[i] - lb[j]) / 2) end t = tref[] # A supported row or column with zero scale gives la (or lb) = -Inf, hence t = +Inf. @@ -729,3 +923,27 @@ function inflate_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatr end return a, b end + +# Asymmetric uniform inflation over a flattened support; the shift convention +# is that of the matrix method. +function inflate_feasible!(a::AbstractVector, b::AbstractVector, sup::FlatSupport) + T = float(promote_type(eltype(a), eltype(b))) + is, js, lv = sup.is, sup.js, sup.lv + la, lb = map(log, a), map(log, b) + t = zero(T) + for k in eachindex(is, js, lv) + u = (lv[k] - la[is[k]] - lb[js[k]]) / 2 + t = ifelse(u > t, u, t) + end + # A supported row or column with zero scale gives la (or lb) = -Inf, hence t = +Inf. + isfinite(t) || + throw(ArgumentError("inflate_feasible! requires a start with positive scale on every supported row/column")) + iszero(t) && return a, b + for i in eachindex(a) + iszero(a[i]) || (a[i] = exp(la[i] + t)) + end + for j in eachindex(b) + iszero(b[j]) || (b[j] = exp(lb[j] + t)) + end + return a, b +end diff --git a/src/sparse_support.jl b/src/sparse_support.jl index 2af33f0..f5f1cc8 100644 --- a/src/sparse_support.jl +++ b/src/sparse_support.jl @@ -30,6 +30,12 @@ function foreach_support_sym(f, A::SparseMatrixCSC) return nothing end +# Upper bounds for sparse support traversals, used by `sizehint!`. +_support_sizehint(A::SparseMatrixCSC) = nnz(A) +_support_sizehint_sym(A::SparseMatrixCSC) = (nnz(A) + size(A, 1) + 1) >> 1 +_support_sizehint_sym(S::Union{Symmetric{<:Any,<:SparseMatrixCSC},Hermitian{<:Any,<:SparseMatrixCSC}}) = + nnz(parent(S)) + # Match symmetric partners in O(nnz) with one cursor per column. `tr[c]` points # to the first unpaired entry in column `c`; absent entries are zeros. function require_abs_symmetric(A::SparseMatrixCSC, fname) From dfd77f75fc31a2732b3f8c24365c1c6b450fd4dd Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 24 Aug 2026 17:18:58 -0500 Subject: [PATCH 15/15] Cap the LSQR preconditioner's factor Use CHOLMOD's symbolic fill estimate to cap Cholesky preconditioner storage at `fillbudget`. Fall back to diagonal preconditioning and report the selected regime in the statistics. Assisted-by: Claude Opus (claude-opus-5) --- src/minimal_covers.jl | 120 ++++++++++++++++++++++++++++++----------- test/minimal_covers.jl | 37 +++++++++++++ test/runtests.jl | 3 +- 3 files changed, 128 insertions(+), 32 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index fc79a04..25c31ed 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -25,7 +25,7 @@ method selects the one with the smallest `AbsLog{2}` objective. # Extended help The native solver accepts `κs` (penalty-continuation schedule), `maxiter` -(Newton steps per stage), and `linsolve`: +(Newton steps per stage), `fillbudget` (see below), and `linsolve`: - `:dense` factorizes dense normal equations at O(n³) per Newton step. - `:woodbury` handles nearly dense `Float64` support as a sparse correction. It @@ -38,6 +38,10 @@ The native solver accepts `κs` (penalty-continuation schedule), `maxiter` `κs` defaults to a geometric schedule ending at `1e8`: eight stages for exact solves and four for `:lsqr`. An explicit `κs` overrides this default. +For `Float64`, `:lsqr` uses a Cholesky preconditioner when its predicted storage +does not exceed `fillbudget` bytes (default `2^30`). Otherwise it uses a diagonal +preconditioner. The returned statistics identify the choice as `precond`. + If a stage reaches `maxiter`, the solver warns that the cover may not minimize the objective. Increase `maxiter` or supply more continuation stages. @@ -71,11 +75,11 @@ method selects the one with the smallest `AbsLog{2}` objective. # Extended help -The native solver accepts the same `κs`, `maxiter`, and `linsolve` keywords as -[`symcover_min`](@ref), with the same solver-dependent default schedule and the -same warning when a stage runs out of Newton steps. For `:woodbury`, an `m × n` -matrix may omit at most `min(m,n) ÷ 4` entries per row or column and -`4 * max(m,n)` entries in total. +The native solver accepts the same `κs`, `maxiter`, `fillbudget`, and `linsolve` +keywords as [`symcover_min`](@ref), with the same solver-dependent default +schedule, the same preconditioner budget, and the same warning when a stage runs +out of Newton steps. For `:woodbury`, an `m × n` matrix may omit at most +`min(m,n) ÷ 4` entries per row or column and `4 * max(m,n)` entries in total. `:dense` costs O((m+n)³) per Newton step; sparse matrices default to `:lsqr`. See also: [`symcover_min`](@ref), [`cover`](@ref), [`cover_min!`](@ref). @@ -802,11 +806,30 @@ _precond_pattern(::Type{T}, ::Grid, v0, N::Int, mult) where {T} = spzeros(T, 0, # Scale-relative ridge for a positive-definite preconditioner. _precond_ridge(dmax::T) where {T} = (dmax > 0 ? dmax : oneunit(T)) * sqrt(eps(T)) +# Default storage limit, in bytes, for the LSQR Cholesky preconditioner. +# Tripping this switches to diagonal preconditioning, reducing memory +# consumption but increasing the number of iterations for convergence. +const LSQR_FILL_BUDGET = 1 << 30 + +# Return CHOLMOD's symbolic factorization and its predicted number of values. +function _precond_analysis(M::SparseMatrixCSC) + F = SparseArrays.CHOLMOD.symbolic(SparseArrays.CHOLMOD.Sparse(Symmetric(M))) + s = unsafe_load(pointer(F)) + Int(s.n) == size(M, 1) || + error("CHOLMOD analyzed a matrix of order $(Int(s.n)), but `M` has order $(size(M, 1))") + s.is_super == 0 || return F, Int(s.xsize) + counts = unsafe_wrap(Array, convert(Ptr{_factor_index(F)}, s.ColCount), Int(s.n)) + return F, sum(Int, counts) +end + +_factor_index(::SparseArrays.CHOLMOD.Factor{<:Any,Ti}) where {Ti} = Ti + # `AbsLog{2}` penalty continuation. Each stage freezes residual weights, solves # the weighted least-squares problem, and backtracks. `boost=true` applies a final # feasibility shift. The support layout selects the inner solver. function _abslog2_continuation(sys::SupportSystem{T}, x0; - κs, maxiter::Int, linsolve::Symbol, boost::Bool) where {T} + κs, maxiter::Int, linsolve::Symbol, boost::Bool, + fillbudget::Real=LSQR_FILL_BUDGET) where {T} N = sys.N supp = sys.supp v0 = sys.v0 @@ -863,24 +886,23 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; px = zeros(T, use_precond ? N : 0) # scale vector recovered from the LSQR variable pg = zeros(T, use_precond ? N : 0) # `Rᵀ√W y` before the preconditioner is applied mdiag = zeros(T, use_precond ? N : 0) # weighted degrees, the preconditioner's diagonal - # Right-preconditioner: a Cholesky factor of the weighted, gauge-augmented - # normal matrix. Its sparsity pattern is the support's and does not depend on - # the weights, so one symbolic analysis serves the whole continuation and each - # solve refactors in place. The gauge is a dense rank-1 term that would fill - # the factor, so it enters as its diagonal alone and LSQR absorbs the rest. + # The normal-matrix pattern is constant, so one symbolic analysis serves all + # stages. Use its diagonal if the predicted Cholesky factor exceeds the budget. Msp = _precond_pattern(T, supp, v0, use_precond ? N : 0, mult) - # Positions of the entries each solve overwrites: the diagonal, and both - # copies of each off-diagonal support entry. - dpos = use_precond ? [_nzindex(Msp, p, p) for p in 1:N] : Int[] - epos = zeros(Int, use_precond ? 2 * ne : 0) - if use_precond + MF, fill_entries = use_precond ? _precond_analysis(Msp) : (nothing, 0) + use_factor = use_precond && sizeof(T) * fill_entries <= fillbudget + # Positions of the entries each factored solve overwrites: the diagonal, and + # both copies of each off-diagonal support entry. + dpos = use_factor ? [_nzindex(Msp, p, p) for p in 1:N] : Int[] + epos = zeros(Int, use_factor ? 2 * ne : 0) + if use_factor for (e, (p, q)) in enumerate(supp.edges) p == q && continue epos[2 * e - 1] = _nzindex(Msp, p, q) epos[2 * e] = _nzindex(Msp, q, p) end end - MF = use_precond ? cholesky(Symmetric(Msp)) : nothing + psqrt = zeros(T, use_factor ? 0 : (use_precond ? N : 0)) # `K` of the diagonal preconditioner rhs = zeros(T, use_woodbury ? N : 0, size(U, 2) + 1) dmin = use_woodbury ? minimum(sys.dfull) : oneunit(T) cgx = zeros(T, use_woodbury ? N : 0) @@ -959,11 +981,9 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; end g = ne + 1 # index of the appended gauge row if use_precond - # Refill the preconditioner with the weights this solve freezes and - # refactor it in place: the pattern, and so the symbolic analysis, - # is the same on every solve. - nzv = nonzeros(Msp) + # Weighted degrees, the diagonal of `RᵀWR`. fill!(mdiag, zero(T)) + nzv = nonzeros(Msp) for (e, (p, q)) in enumerate(edges) w = ws[e]^2 if p == q @@ -971,10 +991,17 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; else mdiag[p] += w mdiag[q] += w - nzv[epos[2 * e - 1]] = w - nzv[epos[2 * e]] = w + if use_factor + nzv[epos[2 * e - 1]] = w + nzv[epos[2 * e]] = w + end end end + end + if use_factor + # Refill the preconditioner with the weights this solve freezes + # and refactor it in place: the pattern, and so the symbolic + # analysis, is the same on every solve. dmax = zero(T) for p in 1:N mdiag[p] += v0[p]^2 @@ -1011,6 +1038,35 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; soly, it = _lsqr(Pmul!, Ptmul!, cv, Kc \ px) nlsqr[] += it return (Uc \ soly)::Vector{T} + elseif use_precond + # An unknown outside the support gets an identity row. + for p in 1:N + d = mdiag[p] + v0[p]^2 + psqrt[p] = sqrt(d > 0 ? d : oneunit(T)) + end + # Diagonal `K` needs only elementwise scaling. + Dmul! = function (y, yv) + @. px = yv / psqrt + for (e, (p, q)) in enumerate(edges) + y[e] = ws[e] * (px[p] + px[q]) + end + y[g] = dot(v0, px) + return y + end + Dtmul! = function (z, y) + fill!(pg, zero(T)) + for (e, (p, q)) in enumerate(edges) + t = ws[e] * y[e] + pg[p] += t + pg[q] += t + end + @. pg += v0 * y[g] + @. z = pg / psqrt + return z + end + soly, it = _lsqr(Dmul!, Dtmul!, cv, psqrt .* x) + nlsqr[] += it + return soly ./ psqrt end Amul! = function (y, xx) for (e, (p, q)) in enumerate(edges) @@ -1108,14 +1164,16 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; end return x, (; nsolves=nsolves[], lsqriters=nlsqr[], cgiters=ncg[], cholsolves=nchol[], exits=Tuple(exits), stagedrops=Tuple(drops), - linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense)) + linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense), + precond=(!use_precond ? :none : use_factor ? :factor : :diagonal)) end # Worker for `symcover_min(::AbsLog{2})`, returning `(a, stats)`. A supplied # `start` replaces the cold initial solve. Narrow types compute in `Float64`. function _symcover_min_abslog2(A::AbstractMatrix; κs=nothing, maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, - boost::Bool=true, fname=:symcover_min) + boost::Bool=true, fillbudget::Real=LSQR_FILL_BUDGET, + fname=:symcover_min) linsolve in (:auto, :dense, :lsqr, :woodbury) || throw(ArgumentError("linsolve must be :auto, :dense, :lsqr, or :woodbury; got :$linsolve")) # Shared symmetry check for native symmetric minimal covers. @@ -1127,7 +1185,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=nothing, # Continuation tolerances require at least Float64 resolution. if eps(T) > eps(Float64) a64, stats = _symcover_min_abslog2(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A); - κs, maxiter, linsolve, start, boost, fname) + κs, maxiter, linsolve, start, boost, fillbudget, fname) # Narrowing rounds to nearest and so can round a product below its entry. a = T.(a64) boost && _certify_cover!(a, A, fname) @@ -1216,7 +1274,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=nothing, zeros(T, n)) x0 = start === nothing ? nothing : T[hassupp[ip] ? log(T(start[i])) : zero(T) for (ip, i) in enumerate(ax)] - α, stats = _abslog2_continuation(sys, x0; κs=κsched, maxiter, linsolve, boost) + α, stats = _abslog2_continuation(sys, x0; κs=κsched, maxiter, linsolve, boost, fillbudget) _warn_truncated(fname, κsched, stats, maxiter) # Dense scale vector matching cover/symcover; `similar(A, …)` is a SparseVector for sparse A. a = similar(Array{T}, ax) @@ -1230,7 +1288,7 @@ end # Worker for `cover_min(::AbsLog{2})`, returning `(a, b, stats)`. function _cover_min_abslog2(A::AbstractMatrix; κs=nothing, maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, - boost::Bool=true) + boost::Bool=true, fillbudget::Real=LSQR_FILL_BUDGET) linsolve in (:auto, :dense, :lsqr, :woodbury) || throw(ArgumentError("linsolve must be :auto, :dense, :lsqr, or :woodbury; got :$linsolve")) axr = axes(A, 1) @@ -1240,7 +1298,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=nothing, # Continuation tolerances require at least Float64 resolution. if eps(T) > eps(Float64) a64, b64, stats = _cover_min_abslog2(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A); - κs, maxiter, linsolve, start, boost) + κs, maxiter, linsolve, start, boost, fillbudget) # Narrowing rounds to nearest and so can round a product below its entry. a, b = T.(a64), T.(b64) boost && _certify_cover!(a, b, A, :cover_min) @@ -1358,7 +1416,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=nothing, end s0 end - x, stats = _abslog2_continuation(sys, x0; κs=κsched, maxiter, linsolve, boost) + x, stats = _abslog2_continuation(sys, x0; κs=κsched, maxiter, linsolve, boost, fillbudget) _warn_truncated(:cover_min, κsched, stats, maxiter) # Apply the balance convention independently to each support component. rowcomp, colcomp, ncomp, _, _ = _support_components(A) diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 00a480a..10fab1e 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -357,6 +357,43 @@ end @test td8.nsolves <= tl8.nsolves - length(κs8) end +# The LSQR preconditioner runs in two regimes: a Cholesky factor of the weighted +# normal matrix while the fill budget allows it, and that matrix's diagonal +# beyond it. Both must reach the same cover. +@testset "MMC :lsqr preconditioner regimes" begin + rng = StableRNG(23) + n = 200 + Ssp = sprandn(rng, n, n, 8 / (2n)) + Ssp = Ssp + Ssp' + A = SparseMatrixCSC(size(Ssp)..., Ssp.colptr, Ssp.rowval, exp.(Ssp.nzval)) + Gsp = sprandn(rng, n, n, 8 / n) + G = SparseMatrixCSC(size(Gsp)..., Gsp.colptr, Gsp.rowval, exp.(Gsp.nzval)) + + af, sf = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) + ad, sd = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, fillbudget=0) + @test sf.precond === :factor + @test sd.precond === :diagonal + @test af ≈ ad rtol=1e-6 + @test iscover(af, af, A) && iscover(ad, ad, A) + @test cover_objective(AbsLog{2}(), af, af, A) ≈ cover_objective(AbsLog{2}(), ad, ad, A) rtol=1e-8 + # The factored regime is the one that converges in a few iterations per solve. + @test sf.lsqriters < sd.lsqriters + + gf, hf, tf = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr) + gd, hd, td = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr, fillbudget=0) + @test (tf.precond, td.precond) === (:factor, :diagonal) + @test gf .* hf' ≈ gd .* hd' rtol=1e-6 + @test iscover(gf, hf, G) && iscover(gd, hd, G) + @test cover_objective(AbsLog{2}(), gf, hf, G) ≈ cover_objective(AbsLog{2}(), gd, hd, G) rtol=1e-8 + + # The budget is a keyword of the public solvers, and the paths that never + # precondition say so. + @test symcover_min(AbsLog{2}(), A; linsolve=:lsqr, fillbudget=0) ≈ ad rtol=1e-6 + @test MatrixCovers._symcover_min_abslog2(A; linsolve=:dense)[2].precond === :none + Abig = BigFloat.([4.0 1.0 0.5; 1.0 3.0 1.0; 0.5 1.0 2.5]) + @test MatrixCovers._symcover_min_abslog2(Abig; linsolve=:lsqr)[2].precond === :none +end + # LSQR continuation starts from the heuristic cover. @testset "MMC :lsqr continuation starts from the heuristic cover" begin rng = StableRNG(17) diff --git a/test/runtests.jl b/test/runtests.jl index 3b03953..9cc0c66 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,7 +44,8 @@ include("helpers.jl") # isbalanced, covaries, PENALTIES :_edge_list, :_sym_edge_list, :_degrees, :_balance_cover!, :inflate_feasible!) # External non-public names with no usable public equivalent. - foreign = (:FreeUnits, :Unit, :Units, :Optimizer, :Experimental, :register_error_hint) + foreign = (:FreeUnits, :Unit, :Units, :Optimizer, :Experimental, :register_error_hint, + :CHOLMOD, :symbolic) test_explicit_imports( MatrixCovers; all_explicit_imports_are_public = VERSION >= v"1.11" ?