Make the symmetric heuristic start scale-covariant

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

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
 

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)