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 056a566..beb8848 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, nnz, nonzeros, nzrange, rowvals, sparse, spzeros export AbsLog, AbsLinear export cover_objective, iscover @@ -24,7 +24,9 @@ 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 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..2a092a9 --- /dev/null +++ b/src/dense_heuristic.jl @@ -0,0 +1,329 @@ +# Dense-grid kernels for the heuristic covers. +# +# 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. +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 + +# 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 + o = _trioff(jp) + sj = zero(T) + cj = 0 + for ip in 1:jp-1 + l = Lp[o+ip] + if l != ninf + s[ip] += l + cnt[ip] += 1 + sj += l + cj += 1 + end + end + s[jp] += sj + cnt[jp] += cj + l = Lp[o+jp] + if l != ninf + s[jp] += l + cnt[jp] += 1 + end + end + return Lp +end + +# 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 + sj = zero(T) + cj = 0 + for ip in 1:m + l = L[ip, jp] + if l != ninf + sa[ip] += l + na[ip] += 1 + sj += l + cj += 1 + end + end + sb[jp] = sj + nb[jp] = cj + end + 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 + +# 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) + 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)) + +# 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")) + _dense_boost!(α, lα, _tri_violated(Lp, lα, n, nviol, _grid_label(n, n)), 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 + + 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 = _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) + 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, _grid_components(A, na, nb, m, n)...) + 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/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 24d4ee4..125da45 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -64,9 +64,21 @@ 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))")) - unconstrained_min!(AbsLog{2}(), a, A) - boost_feasible!(a, A) - return tighten_cover!(a, A; kwargs...) + return _symcover!(a, A; kwargs...) +end + +function _symcover!(a::AbstractVector, A::AbstractMatrix; maxiter::Int=3) + T = float(real(eltype(a))) + if _use_dense_grid(A, T) + _symcover_dense!(a, A, T, maxiter) + else + 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) end """ @@ -137,12 +149,23 @@ 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)))")) - unconstrained_min!(AbsLog{2}(), a, b, A) - boost_feasible!(a, b, A) - tighten_cover!(a, b, A; kwargs...) - # Apply the package's balance convention, then restore coverage lost to rounding. - _balance_cover!(a, b, A) - return inflate_feasible!(a, b, A) + 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))) + if _use_dense_grid(A, T) + _cover_dense!(a, b, A, T, maxiter) + else + 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, sup) + end + return _certify_cover!(a, b, A, :cover) end # Adjoint/Transpose wrappers for cover!. @@ -158,33 +181,111 @@ end # ============================================================ # Internal helpers # ============================================================ +# 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!(is, i); push!(js, j); push!(lv, T(v)) + end + _fastlog!(lv) # one vectorized pass over the collected magnitudes + return FlatSupport(is, js, lv) +end + +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!(is, i); push!(js, j); push!(lv, T(v)) + end + _fastlog!(lv) # one vectorized pass over the collected magnitudes + return FlatSupport(is, js, lv) +end + +# 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 # exactly, at the cost of balancing only within a factor of `sqrt(2)`. function _balance_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) + 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))) - rowcomp, colcomp, ncomp = _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 @@ -221,6 +322,32 @@ function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, A::AbstractMatrix return nza end +# 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)))")) @@ -252,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" — @@ -310,6 +470,32 @@ function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; maxiter::Int=3) return a end +# 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")) @@ -347,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...) @@ -361,45 +578,54 @@ 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) +# 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)) + # `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) z > zero(T) || continue - zmax = max(zmax, z) + ptr[bucketof(z)+1] += 1 end - 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[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,20 +653,24 @@ 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) + nfill = Ref(0) foreach_support_sym(A) do i, j, v lv = log(T(v)) 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] @@ -449,7 +679,37 @@ 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 + +# 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 @@ -470,20 +730,24 @@ 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) + nfill = Ref(0) foreach_support(A) do i, j, v lv = log(T(v)) 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] @@ -492,7 +756,34 @@ 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 + +# 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 @@ -632,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/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 68748da..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 @@ -35,6 +35,16 @@ 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. + +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. + The native solver computes in `Float64` for narrower input types, then converts the result to the required element type. @@ -65,8 +75,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). For `:woodbury`, an `m × n` matrix may omit at most +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`. @@ -249,25 +261,21 @@ 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. # - `: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. 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) @@ -492,21 +500,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 +530,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 +598,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 +630,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,24 +668,168 @@ 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 + +# 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 + 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 + +# 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)) + +# 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 @@ -689,40 +856,53 @@ 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 + 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 - # 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 + # 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 - for p in 1:N - dpart[p] += v0[p]^2 + 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 - for p in 1:N - dpart[p] > 0 || (dpart[p] = oneunit(T)) + end + # Storage for the sparse correction, reused across solves. + Ccolptr = zeros(Int, use_woodbury ? N + 1 : 0) + Crowval = Int[] + Cnzval = 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 + # 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) + 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 - 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[] - # COO triplets of `C`, refilled whenever a Woodbury solve is factorized. - Ci = Int[] - Cj = Int[] - Cv = T[] + 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) @@ -743,8 +923,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 +940,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,13 +967,10 @@ 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 - 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 @@ -831,80 +978,40 @@ 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 + # 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 - mdiag[p] += 4 * oneunit(T) + mdiag[p] += 4 * w else - w = mult(p, q) * oneunit(T) mdiag[p] += w mdiag[q] += w + if use_factor + nzv[epos[2 * e - 1]] = w + nzv[epos[2 * e]] = w + end 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) + 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 - κ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]) + nzv[dpos[p]] = mdiag[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 - 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. @@ -931,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) @@ -985,8 +1121,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) @@ -999,11 +1140,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 @@ -1013,15 +1163,17 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; end end return x, (; nsolves=nsolves[], lsqriters=nlsqr[], cgiters=ncg[], - cholsolves=nchol[], - linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense)) + cholsolves=nchol[], exits=Tuple(exits), stagedrops=Tuple(drops), + 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=(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) + 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. @@ -1033,21 +1185,28 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # 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) - return T.(a64), stats + κ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) + return a, stats 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. @@ -1069,17 +1228,22 @@ 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 + # 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) - 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 +1257,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. @@ -1114,19 +1274,21 @@ 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) + α, 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) 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 # 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) + 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) @@ -1136,26 +1298,26 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # 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) - return T.(a64), T.(b64), stats + κ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) + return a, b, stats end m = length(axr) n = length(axc) 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 @@ -1187,16 +1349,21 @@ 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 + # 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) - 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 +1387,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. @@ -1253,9 +1416,10 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end s0 end - x, stats = _abslog2_continuation(sys, x0; κs, 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) + rowcomp, colcomp, ncomp, _, _ = _support_components(A) Lα = zeros(T, ncomp) Lβ = zeros(T, ncomp) nec = zeros(Int, ncomp) @@ -1283,6 +1447,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/src/sparse_support.jl b/src/sparse_support.jl index 4121d1d..f5f1cc8 100644 --- a/src/sparse_support.jl +++ b/src/sparse_support.jl @@ -30,6 +30,48 @@ 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) + 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 2e6d796..67d2505 100644 --- a/src/support.jl +++ b/src/support.jl @@ -97,28 +97,60 @@ 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 +# 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 -# 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. -# -# 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)) 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 +160,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 +175,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 +187,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 +218,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 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 diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 31af75d..10fab1e 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -330,10 +330,17 @@ 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 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) + 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)) gd, hd, td = MatrixCovers._cover_min_abslog2(G; linsolve=:dense) @@ -342,13 +349,75 @@ 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 + 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 -# The Float64 LSQR preconditioner includes κ-weighted rows; other types use the -# plain matrix-free iteration. +# 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) + 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 + +# 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) @@ -356,14 +425,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]) @@ -575,3 +653,44 @@ 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) == 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. + @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) == 8 + @test all(!=(:maxiter), t.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 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" ? 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)