From 6852f07ea04b7403c7d379edc1b3f3cc6dcf5509 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Thu, 3 Sep 2026 05:05:02 -0500 Subject: [PATCH 1/3] Make the symmetric heuristic start scale-covariant The unconstrained AbsLog{2} start for symcover now fits the diagonally normalized entries |A_ij|/sqrt(|A_ii A_jj|) and rescales by sqrt|A_ii|. The start, and hence the whole heuristic, is scale-covariant whenever every component of the support contains a nonzero diagonal entry. Assisted-by: Claude Fable 5.1 --- docs/src/index.md | 24 +++++-- src/dense_heuristic.jl | 49 ++++++++++--- src/heuristic_covers.jl | 144 ++++++++++++++++++++++++++++----------- src/initializers.jl | 7 +- test/heuristic_covers.jl | 43 ++++++++++++ 5 files changed, 211 insertions(+), 56 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 1d7647c..465953d 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -136,11 +136,10 @@ per iteration. ### Covariance of the heuristics -The heuristic solvers are exactly covariant when every row and column has the -same nonzero pattern, including dense matrices without zeros. On irregular -sparse support they may be only approximately covariant: +[`symcover`](@ref) is exactly covariant whenever every connected component of +the support contains a nonzero diagonal entry: -```jldoctest +```jldoctest covariance julia> using MatrixCovers, LinearAlgebra julia> A = [1.0 1 0; 1 1 1; 0 1 1]; @@ -152,7 +151,22 @@ julia> a1 = symcover(A); a2 = symcover(D * A * D); julia> P1 = (d .* a1) * (d .* a1)'; P2 = a2 * a2'; julia> round.(extrema(P2 ./ P1); digits=3) -(1.0, 1.077) +(1.0, 1.0) +``` + +[`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: + +```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> Q1 = (d .* r1) * (e .* c1)'; Q2 = r2 * c2'; + +julia> round.(extrema(Q2 ./ Q1); digits=3) +(0.872, 1.205) ``` Use [`symcover_min`](@ref) or [`cover_min`](@ref) when exact covariance is diff --git a/src/dense_heuristic.jl b/src/dense_heuristic.jl index 2a092a9..b04e50b 100644 --- a/src/dense_heuristic.jl +++ b/src/dense_heuristic.jl @@ -20,9 +20,8 @@ _use_dense_grid(A::AbstractMatrix, ::Type{T}) where {T} = # 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} +# Pack the upper triangle as log magnitudes, with `-Inf` for zeros. +function _tri_logabs!(Lp::Vector{T}, A::AbstractMatrix) where {T} ax = axes(A, 1) or = first(ax) - 1 n = length(ax) @@ -34,16 +33,24 @@ function _tri_logabs!(Lp::Vector{T}, s::Vector{T}, cnt::Vector{Int}, A::Abstract end end _fastlog!(Lp) + return Lp +end + +# Row sums of the diagonally normalized log magnitudes `Lp[i,j] - ρ[i] - ρ[j]` +# and the row support counts; see `_sym_unconstrained!`. +function _tri_normsums!(s::Vector{T}, cnt::Vector{Int}, Lp::Vector{T}, ρ::Vector{T}, n::Int) where {T} fill!(s, zero(T)) fill!(cnt, 0) ninf = T(-Inf) for jp in 1:n o = _trioff(jp) + rj = ρ[jp] sj = zero(T) cj = 0 for ip in 1:jp-1 l = Lp[o+ip] if l != ninf + l -= ρ[ip] + rj s[ip] += l cnt[ip] += 1 sj += l @@ -54,11 +61,11 @@ function _tri_logabs!(Lp::Vector{T}, s::Vector{T}, cnt::Vector{Int}, A::Abstract cnt[jp] += cj l = Lp[o+jp] if l != ninf - s[jp] += l + s[jp] += l - 2 * rj cnt[jp] += 1 end end - return Lp + return s, cnt end # Fill a log-magnitude grid and the row and column summaries. @@ -141,9 +148,8 @@ function _grid_violated(L::Matrix{T}, lα::Vector{T}, lβ::Vector{T}, m::Int, n: 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)) +# The asymmetric start has no reference shift; see `_uncon_scale` in heuristic_covers.jl. +_uncon_scale(si::T, ni::Int, halfmu::T) where {T} = _uncon_scale(si, ni, halfmu, zero(T)) # Greedy boost that updates scales and log scales together. function _dense_boost!(α::Vector{T}, lα::Vector{T}, entries, zmax::T) where {T} @@ -166,11 +172,34 @@ function _symcover_dense!(a::AbstractVector, A::AbstractMatrix, ::Type{T}, maxit α = 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 + _tri_logabs!(Lp, A) + # Covariant reference from the diagonal; see `_sym_reference!`. + ρ = Vector{T}(undef, n) + for jp in 1:n + l = Lp[_trioff(jp)+jp] + ρ[jp] = ifelse(l == T(-Inf), T(NaN), l / 2) + end + _tri_normsums!(α, cnt, Lp, ρ, n) # `α` carries the row log sums here + nmissing = count(ip -> cnt[ip] > 0 && isnan(ρ[ip]), 1:n) + if nmissing > 0 + # Rare path: some supported row has a zero diagonal entry. + function foreach_entries(f) + for jp in 1:n + o = _trioff(jp) + for ip in 1:jp + l = Lp[o+ip] + l == T(-Inf) || f(ip, jp, l) + end + end + end + _sym_reference!(ρ, foreach_entries, nmissing) + _tri_normsums!(α, cnt, Lp, ρ, n) + end + # Unsupported rows may keep NaN references; `_uncon_scale` never reads them. nztotal = sum(cnt) halfmu = iszero(nztotal) ? zero(T) : sum(α) / (2 * nztotal) for ip in 1:n - α[ip] = _uncon_scale(α[ip], cnt[ip], halfmu) + α[ip] = _uncon_scale(α[ip], cnt[ip], halfmu, ρ[ip]) lα[ip] = log(α[ip]) end diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index 125da45..e5d72fc 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -13,8 +13,10 @@ Given a square matrix `A` assumed to be symmetric, return a vector `a` representing a symmetric hard cover of `A`: `a[i] * a[j] >= abs(A[i, j])` for all `i`, `j`. -The method initializes from per-row geometric means, covers the most-violated -entries first, then applies `maxiter` tightening iterations. +The method initializes from per-row geometric means of the diagonally +normalized entries `abs(A[i, j]) / sqrt(abs(A[i, i] * A[j, j]))`, rescaled by +`sqrt(abs(A[i, i]))`, covers the most-violated entries first, then applies +`maxiter` tightening iterations. `ϕ` is accepted for API compatibility but is currently ignored. For a cover that provably minimizes a given `ϕ`, use [`symcover_min`](@ref). @@ -36,6 +38,12 @@ julia> a * a' # covers |A|: a[i]*a[j] >= abs(A[i, j]) 4.0 4.0 4.0 4.0 ``` + +# Extended help + +The result is scale-covariant whenever every connected component of the support has a +nonzero diagonal entry (rows with a zero diagonal take their reference from +neighbors that have one). """ symcover(ϕ::AbstractCoverPenalty, A::AbstractMatrix; kwargs...) = symcover(A; kwargs...) @@ -293,59 +301,119 @@ function _balance_cover!(a::AbstractVector, b::AbstractVector, rowcomp::Vector{I end -# Analytical minimizer of the unconstrained `AbsLog{2}` symmetric objective +# Symmetric unconstrained `AbsLog{2}` start. The objective # ∑_{i,j: A[i,j]≠0} (log(a[i]*a[j]) - log|A[i,j]|)² -# Returns row support counts. The Sherman-Morrison approximation is exact on -# complete support. -function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, A::AbstractMatrix) where T +# is stationary where (diag(n) + S) α = 𝔞, with `n` the row support counts, `S` +# the support indicator, and 𝔞 the row sums of log|A_ij|. The start is one +# Jacobi-type sweep of that system, α = diag(n)⁻¹(𝔞 - S ρ) + c e, from the +# reference point ρ of `_sym_reference!`, with the constant `c` fixed by the +# balance nᵀα = eᵀ𝔞/2 that every exact solution satisfies. Equivalently, it is +# the per-row geometric mean of the diagonally normalized entries +# |A_ij| / sqrt(|A_ii A_jj|), scaled back by sqrt|A_ii|. On complete support the +# sweep is exact (Sherman–Morrison), and because ρ co-varies with `A`, the +# start is exactly scale-covariant on any support whose components each hold a +# nonzero diagonal entry. +# +# Returns the row support counts. +function _sym_unconstrained!(a::AbstractVector{T}, foreach_entries::F) where {T,F} ax = eachindex(a) - axes(A) == (ax, ax) || throw(DimensionMismatch("`unconstrained_min!(ϕ, a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(string(axes(A))), axes(a)=$(string(axes(a)))")) - loga = fill!(similar(a), zero(T)) - nza = zeros(Int, ax) - foreach_support_sym(A) do i, j, v - lAij = log(T(v)) - loga[i] += lAij - nza[i] += 1 - if i != j - loga[j] += lAij - nza[j] += 1 + nza = zeros(Int, ax) + ρ = fill!(similar(a, T), T(NaN)) # NaN marks a row without a reference + foreach_entries() do i, j, lv + nza[i] += 1 + if i == j + ρ[i] = lv / 2 + else + nza[j] += 1 end end + nmissing = count(i -> !iszero(nza[i]) && isnan(ρ[i]), ax) + _sym_reference!(ρ, foreach_entries, nmissing) + loga = fill!(similar(a, T), zero(T)) + foreach_entries() do i, j, lv + l = lv - ρ[i] - ρ[j] + loga[i] += l + i == j || (loga[j] += l) + end nztotal = sum(nza) halfmu = iszero(nztotal) ? zero(T) : sum(loga) / (2 * nztotal) for i in ax - # 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)) + a[i] = _uncon_scale(loga[i], nza[i], halfmu, ρ[i]) end return nza end +# Keep supported scales positive when `exp` underflows. +_uncon_scale(si::T, ni::Int, halfmu::T, ρi::T) where {T} = + iszero(ni) ? zero(T) : max(exp(si / ni - halfmu + ρi), floatmin(T)) + +# Covariant reference log-scales. The caller is expected to initialize +# {log|A_ii|/2 if A_ii ≠ 0 +# ρ[i] = { +# {NaN otherwise +# `nmissing` counts supported rows (those with some A_ij ≠ 0) for which +# A_ii == 0. +# +# Each pass assigns every such row the mean of `log|A_ik| - ρ[k]` over its +# already-referenced neighbors `k`, so the reference spreads outward by graph +# distance from the diagonal. A pass costs one traversal, so the total is +# proportional to the largest graph distance from the diagonal. +# +# Rows in components with no nonzero diagonal entry never receive a reference +# and are set to zero, on which the start is not covariant. +function _sym_reference!(ρ::AbstractVector{T}, foreach_entries::F, nmissing::Int) where {T,F} + if nmissing > 0 + acc = similar(ρ, T) + cnt = zeros(Int, eachindex(ρ)) + while nmissing > 0 + fill!(acc, zero(T)) + fill!(cnt, 0) + foreach_entries() do i, j, lv + i == j && return + ri, rj = ρ[i], ρ[j] + # Only entries joining a referenced row to an unreferenced one contribute. + if isnan(ri) && !isnan(rj) + acc[i] += lv - rj + cnt[i] += 1 + elseif isnan(rj) && !isnan(ri) + acc[j] += lv - ri + cnt[j] += 1 + end + end + nnew = 0 + for i in eachindex(ρ) + if cnt[i] > 0 + ρ[i] = acc[i] / cnt[i] + nnew += 1 + end + end + nnew == 0 && break + nmissing -= nnew + end + end + for i in eachindex(ρ) + isnan(ρ[i]) && (ρ[i] = zero(T)) + end + return ρ +end + +function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, A::AbstractMatrix) where T + ax = eachindex(a) + axes(A) == (ax, ax) || throw(DimensionMismatch("`unconstrained_min!(ϕ, a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(string(axes(A))), axes(a)=$(string(axes(a)))")) + foreach_entries(f) = foreach_support_sym((i, j, v) -> f(i, j, log(T(v))), A) + return _sym_unconstrained!(a, foreach_entries) +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 + function foreach_entries(f) + for k in eachindex(is, js, lv) + f(is[k], js[k], lv[k]) 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 + return _sym_unconstrained!(a, foreach_entries) end function unconstrained_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix) diff --git a/src/initializers.jl b/src/initializers.jl index 59cf9d8..b75fdf2 100644 --- a/src/initializers.jl +++ b/src/initializers.jl @@ -16,7 +16,8 @@ Build a symmetric starting point for [`symcover_min`](@ref) or `strategy` names the point: -- `:geomean` — geometric means of the nonzero entries in each row. +- `:geomean` — the unconstrained `AbsLog{2}` start of [`symcover`](@ref): + geometric means of the diagonally normalized entries in each row. - `:leaveout` — the geometric mean recomputed with the most-underweighted support entry omitted. It fails if removing that entry empties a row. - `:diagfeasible` — a cover grown from the diagonal by nearest-neighbor @@ -189,8 +190,8 @@ function _leaveout_logmean_init!(a::AbstractVector{T}, A::AbstractMatrix) where # Account for both endpoints of an off-diagonal entry. nza[ibest] > 1 || return false ibest == jbest || nza[jbest] > 1 || return false - # Minimize the reduced-support `AbsLog{2}` objective by Gauss-Seidel. Starting - # from a covariant point preserves covariance at every sweep. + # Minimize the reduced-support `AbsLog{2}` objective by Gauss-Seidel. Each + # sweep preserves covariance, so the result is covariant whenever the start is. α = similar(a) for i in ax α[i] = iszero(nza[i]) ? zero(T) : log(a[i]) diff --git a/test/heuristic_covers.jl b/test/heuristic_covers.jl index 1a40f03..59ad8f7 100644 --- a/test/heuristic_covers.jl +++ b/test/heuristic_covers.jl @@ -241,3 +241,46 @@ end a, b = tighten_cover!(zeros(3), zeros(3), Diagonal([1.0, 2.0, 3.0])) @test all(iszero, a) && all(iszero, b) end + +@testset "symcover scale covariance on irregular support" begin + # The unconstrained start references each row to its diagonal entry, so + # symcover co-varies with D*A*D on any support whose components contain a + # nonzero diagonal entry, not only on complete support. + rng = StableRNG(7) + start(A) = (a = zeros(size(A, 1)); unconstrained_min!(AbsLog{2}(), a, A); a) + for n in (12, 40, 100) # n = 100 uses the dense-grid kernel for Matrix storage + S = sprandn(rng, n, n, 0.2) + A = Matrix(S + S') + Diagonal(randn(rng, n) .+ 0.1) + d = exp.(2 .* randn(rng, n)) + @test covaries(start, A, d; rtol=1e-9) + @test covaries(A -> symcover(A; maxiter=0), A, d; rtol=1e-9) + @test covaries(symcover, A, d; rtol=1e-9) + @test covaries(symcover, sparse(A), d; rtol=1e-9) + @test symcover(sparse(A)) ≈ symcover(A) rtol=1e-9 + end + # A power-of-two rescaling is exact in floating point. + n = 30 + S = sprandn(rng, n, n, 0.2) + A = Matrix(S + S') + Diagonal(randn(rng, n) .+ 0.1) + d = exp2.(rand(rng, -20:20, n)) + @test covaries(symcover, A, d; rtol=4eps()) + + # Rows with a zero diagonal take their reference from neighbors, layer by + # layer: a path whose only diagonal entry sits at one end propagates through + # every row. + for n in (10, 70) + A = Matrix(SymTridiagonal([1.0; zeros(n - 1)], exp.(randn(rng, n - 1)))) + d = exp.(2 .* randn(rng, n)) + @test covaries(start, A, d; rtol=1e-9) + @test covaries(symcover, A, d; rtol=1e-9) + @test iscover(symcover(A), A; rtol=8eps()) + end + + # On complete support the start is the exact minimizer of the unconstrained + # objective: (diag(n) + S) α = 𝔞 with S the all-ones support. + n = 8 + B = randn(rng, n, n); A = B + B' + L = log.(abs.(A)) + α = (n * I + ones(n, n)) \ vec(sum(L, dims=2)) + @test start(A) ≈ exp.(α) rtol=1e-10 +end From 71a3c4c46bbe2b041e227c07c5be0e9a4254f35d Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 4 Sep 2026 04:28:25 -0500 Subject: [PATCH 2/3] Build the leave-one-out start with one reduced-support sweep The :leaveout initializer now reruns the unconstrained AbsLog{2} kernel with the dropped entry excluded from the support iterator, instead of polishing the full-support start with Gauss-Seidel sweeps. Excluding the entry also removes it from the diagonal reference, so a dropped diagonal entry's row takes its reference from its neighbors, and the start equals the :geomean start of the matrix with that entry zeroed. Assisted-by: Claude Fable 5.1 --- src/initializers.jl | 68 +++++++++++++------------------------------- test/initializers.jl | 17 +++++++++++ 2 files changed, 37 insertions(+), 48 deletions(-) diff --git a/src/initializers.jl b/src/initializers.jl index b75fdf2..dad1bb7 100644 --- a/src/initializers.jl +++ b/src/initializers.jl @@ -158,67 +158,39 @@ function _reject_kwargs(strategy::Symbol, kwargs) throw(ArgumentError("strategy=:$strategy accepts no further keyword arguments, got $(string(join(keys(kwargs), ", ")))")) end -# Recompute the geometric mean after dropping the most negative log-residual. -# Residual ties use raw magnitude; this is the strategy's only covariance -# exception. Return `false` if no entry can be removed without emptying a row. +# Recompute the unconstrained start with the most negative log-residual entry +# removed from the support. Residual ties use raw magnitude; this is the +# strategy's only covariance exception. Return `false` if no entry can be +# removed without emptying a row. function _leaveout_logmean_init!(a::AbstractVector{T}, A::AbstractMatrix) where T ax = eachindex(a) axes(A) == (ax, ax) || throw(DimensionMismatch("`_leaveout_logmean_init!(a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(string(axes(A))), axes(a)=$(string(axes(a))))")) nza = unconstrained_min!(AbsLog{2}(), a, A) sum(nza) == 0 && return false - # Pair scans use the upper triangle; Gauss-Seidel uses complete rows. - S = _sym_support(A, T) # Use a roundoff-tolerant set for tied residuals. - zmin = T(Inf) - for i in ax, s in _slots(S, i) - j = S.idx[s] - j < i && continue - zmin = min(zmin, log(S.val[s]) - log(a[i]) - log(a[j])) + zmin = Ref(T(Inf)) + foreach_support_sym(A) do i, j, v + zmin[] = min(zmin[], log(T(v)) - log(a[i]) - log(a[j])) end - ztol = 64 * eps(T) * max(one(T), abs(zmin)) - ibest = jbest = first(ax) - 1 - Abest = T(Inf) - for i in ax, s in _slots(S, i) - j = S.idx[s] - j < i && continue - Aij = S.val[s] + ztol = 64 * eps(T) * max(one(T), abs(zmin[])) + best = Ref((first(ax) - 1, first(ax) - 1, T(Inf))) + foreach_support_sym(A) do i, j, v + Aij = T(v) z = log(Aij) - log(a[i]) - log(a[j]) - if z <= zmin + ztol && Aij < Abest - ibest, jbest, Abest = i, j, Aij + if z <= zmin[] + ztol && Aij < best[][3] + best[] = (i, j, Aij) end end + ibest, jbest, _ = best[] # Account for both endpoints of an off-diagonal entry. nza[ibest] > 1 || return false ibest == jbest || nza[jbest] > 1 || return false - # Minimize the reduced-support `AbsLog{2}` objective by Gauss-Seidel. Each - # sweep preserves covariance, so the result is covariant whenever the start is. - α = similar(a) - for i in ax - α[i] = iszero(nza[i]) ? zero(T) : log(a[i]) - end - for _ in 1:8 - for i in ax - iszero(nza[i]) && continue - num = zero(T) # Σ_j W[i,j] (log|A[i,j]| - α[j]), α[i]-coefficient split out - den = zero(T) - for s in _slots(S, i) - j = S.idx[s] - (min(i, j) == ibest && max(i, j) == jbest) && continue - lAij = log(S.val[s]) - if j == i - num += lAij - den += 2 - else - num += lAij - α[j] - den += 1 - end - end - iszero(den) && continue # row's only support was the dropped entry (guarded above) - α[i] = num / den - end - end - for i in ax - a[i] = iszero(nza[i]) ? zero(T) : exp(α[i]) + # One sweep of the same kernel on the reduced support. Omitting the entry + # from the iterator also omits it from the diagonal reference, so a row + # whose diagonal entry is dropped takes its reference from its neighbors. + foreach_entries(f) = foreach_support_sym(A) do i, j, v + (i == ibest && j == jbest) || f(i, j, log(T(v))) end + _sym_unconstrained!(a, foreach_entries) return true end diff --git a/test/initializers.jl b/test/initializers.jl index e0feef0..cc249c1 100644 --- a/test/initializers.jl +++ b/test/initializers.jl @@ -88,6 +88,23 @@ @test ah[3] < ag[3] end + @testset ":leaveout is :geomean on the reduced support" begin + # Dropping an entry removes it from the row means and from the diagonal + # reference alike, so the start equals the :geomean start of the matrix + # with that entry zeroed. The most-underweighted entry is the one scaled + # by 1e-8, whether off-diagonal or diagonal. + rng = StableRNG(11) + B = Matrix(Symmetric(exp.(randn(rng, 5, 5)))) + for (i, j) in ((2, 4), (3, 3)) + A = copy(B) + A[i, j] = A[j, i] = 1e-8 * B[i, j] + Z = copy(A) + Z[i, j] = Z[j, i] = 0.0 + @test initialize_symcover(A; strategy=:leaveout, feasible=:none) == + initialize_symcover(Z; strategy=:geomean, feasible=:none) + end + end + @testset ":diagfeasible propagates from sparse anchors" begin # Propagate the diagonal scale backward through a path. n = 6 From ff07b717ab11022b2b2dcbe27834be07b2a50565 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 4 Sep 2026 05:07:00 -0500 Subject: [PATCH 3/3] Widen the symmetric solve-count bound The count of inner solves on the dense exact path varies by a few across machines with BLAS roundoff; on one CI runner the symmetric 60x60 case took 35. Give the symmetric bound the same kind of margin the nonsymmetric bound already has. Assisted-by: Claude Fable 5.1 --- test/minimal_covers.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 375b070..b432c5c 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -336,9 +336,11 @@ end @test ad ≈ al rtol=1e-6 @test aw ≈ al rtol=1e-6 # The dense and Woodbury trajectories agree to roundoff, which adaptive - # escalation can turn into slightly different outer-iteration counts. + # escalation can turn into slightly different outer-iteration counts. The + # same roundoff varies the count by a few solves across machines, so the + # absolute bounds below carry a margin over observed counts. @test abs(sd.nsolves - sw.nsolves) <= 10 - @test sd.nsolves <= 34 + @test sd.nsolves <= 40 # Once the active set settles, each outer pass of an exact path ends on its # first, sign-stable Newton step. @test count(==(:stable), sd.exits) >= max(1, sd.nouter ÷ 2) @@ -350,7 +352,6 @@ end @test gd .* hd' ≈ gl .* hl' rtol=1e-6 @test gw .* hw' ≈ gl .* hl' rtol=1e-6 @test abs(td.nsolves - tw.nsolves) <= 10 - # Leave a small margin in the solve-count bound. @test td.nsolves <= 48 @test count(==(:stable), td.exits) >= max(1, td.nouter ÷ 2) end