From 1fad3d2ba844faabd9ce6779273a8a2a97b6cdc6 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 5 Sep 2026 22:23:53 -0500 Subject: [PATCH 1/2] Refine asymmetric starts with conjugate gradients Add `cgiter` (default 4) to `cover` and `cover!` to refine the geometric-mean start toward a scale-covariant least-squares log fit. `cgiter=0` disables refinement. Assisted-by: Claude Fable 5.1 --- src/dense_heuristic.jl | 17 +++++- src/heuristic_covers.jl | 121 +++++++++++++++++++++++++++++++++++++-- test/heuristic_covers.jl | 55 ++++++++++++++++++ 3 files changed, 186 insertions(+), 7 deletions(-) diff --git a/src/dense_heuristic.jl b/src/dense_heuristic.jl index b04e50b..29cd9c7 100644 --- a/src/dense_heuristic.jl +++ b/src/dense_heuristic.jl @@ -253,7 +253,7 @@ 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} + ::Type{T}, maxiter::Int, cgiter::Int) where {T} or = first(axes(A, 1)) - 1 oc = first(axes(A, 2)) - 1 m, n = size(A) @@ -275,6 +275,21 @@ function _cover_dense!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, β[jp] = _uncon_scale(β[jp], nb[jp], halfmu) lβ[jp] = log(β[jp]) end + if cgiter > 0 + function foreach_grid(f) + for jp in 1:n, ip in 1:m + l = L[ip, jp] + isfinite(l) && f(ip, jp, l) # -Inf marks a zero entry + end + end + _cg_refine!(lα, lβ, foreach_grid, cgiter) + for ip in 1:m + iszero(na[ip]) || (α[ip] = max(exp(lα[ip]), floatmin(T)); lα[ip] = log(α[ip])) + end + for jp in 1:n + iszero(nb[jp]) || (β[jp] = max(exp(lβ[jp]), floatmin(T)); lβ[jp] = log(β[jp])) + end + end nviol = 0 zmax = zero(T) diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index e5d72fc..7b41536 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -90,13 +90,16 @@ function _symcover!(a::AbstractVector, A::AbstractMatrix; maxiter::Int=3) end """ - a, b = cover(ϕ, A; maxiter=3) - a, b = cover(A; maxiter=3) + a, b = cover(ϕ, A; maxiter=3, cgiter=4) + a, b = cover(A; maxiter=3, cgiter=4) Given a matrix `A`, return vectors `a` and `b` such that `a[i] * b[j] >= abs(A[i, j])` for all `i`, `j`. The method initializes from row -and column geometric means, covers the most-violated entries first, then applies -`maxiter` tightening iterations. +and column geometric means, refines that start with up to `cgiter` conjugate-gradient +iterations, covers the most-violated entries first, then applies `maxiter` +tightening iterations. + +`cgiter=0` disables refinement. The factors use the per-component balance convention described by [`cover_min`](@ref). @@ -119,6 +122,13 @@ julia> a * b' 2.16541 2.03444 3.0 6.0 5.63709 8.31251 ``` + +# Extended help + +Conjugate-gradient refinement fits `log(a[i]*b[j])` to `log(abs(A[i, j]))` +in least squares over the nonzero entries. The exact fit is scale-covariant; +a finite number of iterations need not achieve this. Each iteration makes +one pass over the support, stopping early when the residual is small. """ cover(ϕ::AbstractCoverPenalty, A::AbstractMatrix; kwargs...) = cover(A; kwargs...) @@ -160,13 +170,15 @@ function cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs. return _cover!(a, b, A; kwargs...) end -function _cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=3) +function _cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=3, cgiter::Int=4) + cgiter >= 0 || throw(ArgumentError("cgiter must be nonnegative, got $cgiter")) T = float(promote_type(eltype(a), eltype(b))) if _use_dense_grid(A, T) - _cover_dense!(a, b, A, T, maxiter) + _cover_dense!(a, b, A, T, maxiter, cgiter) else sup = flat_support(A, T) unconstrained_min!(AbsLog{2}(), a, b, sup) + cg_refine_start!(a, b, sup, cgiter) boost_feasible!(a, b, sup) tighten_cover!(a, b, sup; maxiter) # Apply the package's balance convention, then restore coverage lost to rounding. @@ -480,6 +492,103 @@ function unconstrained_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, s return nza, nzb end +# Refine the log scales by CG on the normal equations for +# ∑_{ij ∈ support} (lα[i] + lβ[j] - log|A_ij|)². +# `foreach_entries(f)` calls `f(i, j, log|A_ij|)` over the support. +# Updates preserve the initial gauge and leave unsupported log scales at -Inf. +function _cg_refine!(lα::AbstractVector{T}, lβ::AbstractVector{T}, foreach_entries::F, maxiter::Int) where {T,F} + maxiter > 0 || return lα, lβ + axa, axb = eachindex(lα), eachindex(lβ) + na = zeros(Int, axa) + nb = zeros(Int, axb) + fa = fill!(similar(lα, T), zero(T)) + fb = fill!(similar(lβ, T), zero(T)) + ra = fill!(similar(lα, T), zero(T)) + rb = fill!(similar(lβ, T), zero(T)) + # Accumulate support counts, right-hand sides, and off-diagonal residuals. + foreach_entries() do i, j, lv + na[i] += 1 + nb[j] += 1 + fa[i] += lv + fb[j] += lv + ra[i] += lv - lβ[j] + rb[j] += lv - lα[i] + end + rr = zero(T) # squared residual norm + rref = zero(T) # squared right-hand-side norm + for i in axa + if iszero(na[i]) + ra[i] = zero(T) + else + ra[i] -= na[i] * lα[i] + rr += ra[i]^2 + rref += fa[i]^2 + end + end + for j in axb + if iszero(nb[j]) + rb[j] = zero(T) + else + rb[j] -= nb[j] * lβ[j] + rr += rb[j]^2 + rref += fb[j]^2 + end + end + tol = eps(T) * rref + pa = copy(ra) + pb = copy(rb) + Apa = similar(ra) + Apb = similar(rb) + for _ in 1:maxiter + rr <= tol && break + for i in axa + Apa[i] = na[i] * pa[i] + end + for j in axb + Apb[j] = nb[j] * pb[j] + end + foreach_entries() do i, j, _ + Apa[i] += pb[j] + Apb[j] += pa[i] + end + pAp = LinearAlgebra.dot(pa, Apa) + LinearAlgebra.dot(pb, Apb) + pAp > 0 || break + γ = rr / pAp + lα .+= γ .* pa + lβ .+= γ .* pb + ra .-= γ .* Apa + rb .-= γ .* Apb + rr2 = LinearAlgebra.dot(ra, ra) + LinearAlgebra.dot(rb, rb) + δ = rr2 / rr + pa .= ra .+ δ .* pa + pb .= rb .+ δ .* pb + rr = rr2 + end + return lα, lβ +end + +# Apply log-space refinement to `a`, `b`, preserving zeros and clamping underflow. +function cg_refine_start!(a::AbstractVector, b::AbstractVector, sup::FlatSupport, maxiter::Int) + maxiter > 0 || return a, b + T = float(promote_type(eltype(a), eltype(b))) + lα = map(x -> log(T(x)), a) + lβ = map(x -> log(T(x)), b) + is, js, lv = sup.is, sup.js, sup.lv + function foreach_entries(f) + for k in eachindex(is, js, lv) + f(is[k], js[k], lv[k]) + end + end + _cg_refine!(lα, lβ, foreach_entries, maxiter) + for i in eachindex(a) + iszero(a[i]) || (a[i] = max(exp(lα[i]), floatmin(T))) + end + for j in eachindex(b) + iszero(b[j]) || (b[j] = max(exp(lβ[j]), floatmin(T))) + end + return a, b +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" — diff --git a/test/heuristic_covers.jl b/test/heuristic_covers.jl index 59ad8f7..281ef9f 100644 --- a/test/heuristic_covers.jl +++ b/test/heuristic_covers.jl @@ -284,3 +284,58 @@ end α = (n * I + ones(n, n)) \ vec(sum(L, dims=2)) @test start(A) ≈ exp.(α) rtol=1e-10 end + +@testset "cover: conjugate-gradient refinement of the start" begin + # Centered log-product deviation under row and column scaling. + function covdev(coverfn, A, dr, dc) + a, b = coverfn(A) + aD, bD = coverfn(dr .* A .* dc') + d = log.((aD .* bD') ./ ((dr .* a) .* (dc .* b)'))[A .!= 0] + return maximum(abs, d .- sum(d) / length(d)) + end + rng = StableRNG(20260905) + # Refinement reduces scaling dependence on banded support. + n = 12 + A = zeros(n, n) + for i in 1:n, j in max(1, i - 1):min(n, i + 1) + A[i, j] = exp(2 * randn(rng)) + end + dr, dc = exp.(4 .* randn(rng, n)), exp.(4 .* randn(rng, n)) + dev0 = covdev(A -> cover(A; cgiter=0), A, dr, dc) + dev4 = covdev(A -> cover(A), A, dr, dc) + devx = covdev(A -> cover(A; cgiter=4n), A, dr, dc) + @test dev0 > 0.1 + @test dev4 < dev0 / 3 + @test devx < 1e-6 + @test covaries(A -> cover(A; cgiter=4n), A, dr, dc; rtol=1e-6) + for k in (0, 4, 4n) + a, b = cover(A; cgiter=k) + @test iscover(a, b, A; rtol=8eps()) + end + @test_throws ArgumentError cover(A; cgiter=-1) + + # The dense-grid and flattened-support paths refine identically. + m = 2 * MatrixCovers.DENSE_GRID_MIN + B = zeros(m, m) + for i in 1:m, j in max(1, i - 2):min(m, i + 1) + B[i, j] = exp(2 * randn(rng)) * (rand(rng) < 0.7) + end + for i in 1:m + B[i, i] = exp(2 * randn(rng)) + end + for k in (0, 4, 40) + ad, bd = cover(B; cgiter=k) + as, bs = cover(sparse(B); cgiter=k) + @test ad .* bd' ≈ as .* bs' rtol = 1e-10 + end + + # Fully populated support needs no refinement. + C = exp.(2 .* randn(rng, 9, 7)) + @test cover(C; cgiter=4) == cover(C; cgiter=0) + @test cover(sparse(C); cgiter=4) == cover(sparse(C); cgiter=0) + + # Empty rows and columns keep zero scales through the refinement. + Z = [0.0 2.0 0.0; 0.0 0.0 0.0; 1.0 0.0 0.0] + a, b = cover(Z) + @test a[2] == 0 && b[3] == 0 && iscover(a, b, Z; rtol=8eps()) +end From 49e4ca93688133c0b2218ffb692d07ca9e0a273e Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 6 Sep 2026 10:26:27 -0500 Subject: [PATCH 2/2] Show the effect of CG refinement on covariance in the docs The covariance example for `cover` now compares `cgiter=0` with the default refinement: the unrefined start is only approximately covariant, while the refined start reaches the covariant fit on this small example. The prose notes that larger problems need not converge within the default iteration count. Assisted-by: Claude Fable 5.1 --- docs/src/index.md | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 465953d..d1c1528 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -156,17 +156,27 @@ julia> round.(extrema(P2 ./ P1); digits=3) [`cover`](@ref) is exactly covariant when every row and column has the same nonzero pattern, including dense matrices without zeros. On irregular sparse -support it may be only approximately covariant: +support its geometric-mean start is only approximately covariant; the +conjugate-gradient refinement (`cgiter`, default 4) moves it toward a +covariant least-squares fit, reaching it here but not necessarily on larger +problems: ```jldoctest covariance julia> e = [2.0, 0.3, 5.0]; E = Diagonal(e); -julia> r1, c1 = cover(A); r2, c2 = cover(D * A * E); +julia> function covariance_spread(cgiter) + r1, c1 = cover(A; cgiter) + r2, c2 = cover(D * A * E; cgiter) + Q1 = (d .* r1) * (e .* c1)' + Q2 = r2 * c2' + return round.(extrema(Q2 ./ Q1); digits=3) + end; -julia> Q1 = (d .* r1) * (e .* c1)'; Q2 = r2 * c2'; - -julia> round.(extrema(Q2 ./ Q1); digits=3) +julia> covariance_spread(0) (0.872, 1.205) + +julia> covariance_spread(4) +(1.0, 1.0) ``` Use [`symcover_min`](@ref) or [`cover_min`](@ref) when exact covariance is