From b0c9d8957e3fb600d973ed26c6848e6c6db47339 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 21 Aug 2026 15:08:27 -0500 Subject: [PATCH 1/2] Use minimal covers for coupled Gram components Symmetrize component block sums and solve them with `symcover_min`. Reject loopless bipartite couplings by default, with `degenerate=:uniform` available as a gauge-dependent fallback. Assisted-by: Claude Opus 5 --- src/gram_covers.jl | 306 +++++++++++++++++++++----------------------- test/gram_covers.jl | 181 ++++++++++++++++++++++++-- test/unitful.jl | 23 ++++ 3 files changed, 344 insertions(+), 166 deletions(-) diff --git a/src/gram_covers.jl b/src/gram_covers.jl index 1900fa4..ac33887 100644 --- a/src/gram_covers.jl +++ b/src/gram_covers.jl @@ -1,6 +1,9 @@ # Symmetric covers of a (weighted) Gram matrix `A'*W*A`, computed directly from # an asymmetric cover of `A` — without ever forming the Gram matrix. +# Default for the general-`W` methods. +const GRAMCOVER_DEGENERATE = :error + # ============================================================ # Public interface # ============================================================ @@ -8,31 +11,26 @@ """ s = gramcover(a, b, A) s = gramcover(a, b, A, w::AbstractVector) - s = gramcover(a, b, A, W::AbstractMatrix) + s = gramcover(a, b, A, W::AbstractMatrix; degenerate=:$(GRAMCOVER_DEGENERATE)) s = gramcover(a, b, sc::SupportComponents) s = gramcover(a, b, sc::SupportComponents, w::AbstractVector) - s = gramcover(a, b, sc::SupportComponents, W::AbstractMatrix) - -Given an asymmetric cover `a[i]*b[j] >= abs(A[i,j])` — from [`cover`](@ref), -[`cover_min`](@ref), or any other solver producing such a pair — return a -symmetric cover `s` of a (weighted) Gram matrix of `A`, without forming that -Gram matrix: `s[j]*s[k] >= abs(G[j,k])` for every `j`, `k`, where `G = A'*A` -for the two-argument form, `G = A'*Diagonal(w)*A` for the vector-weighted -form, and `G = A'*W*A` for the general form. Only `abs.(W)` enters the bound, -so `W` need not be positive definite, positive semidefinite, or even -symmetric; passing `W::Diagonal` is equivalent to passing `W.diag` as `w`. - -`(a, b)` covering `A` is a precondition, not verified here — use -[`iscover`](@ref)`(a, b, A)` to check it beforehand. `gramcover` composes with -any asymmetric cover this package produces: [`cover`](@ref), [`cover_min`](@ref), -[`soft_cover`](@ref), and their mutating and `_min` forms. - -Only the connected components of `A`'s support enter the two- and -vector-weighted forms, so a caller holding an -[`support_components`](@ref)`(A)` result may pass it in place of `A`; the -matrix forms are exactly that call followed by the `sc` form. `sc.rowax` and -`sc.colax` then play the roles of `axes(A, 1)` and `axes(A, 2)` in the axis -requirements on `a`, `b`, and `w`/`W`. + s = gramcover(a, b, sc::SupportComponents, W::AbstractMatrix; degenerate=:$(GRAMCOVER_DEGENERATE)) + +Given an asymmetric cover `a[i]*b[j] >= abs(A[i,j])`, return a symmetric cover +`s` of a weighted Gram matrix without forming it: `s[j]*s[k] >= abs(G[j,k])`, +where `G` is `A'*A`, `A'*Diagonal(w)*A`, or `A'*W*A`. Only `abs.(W)` enters the +bound, so `W` need not be symmetric or positive semidefinite. Passing +`W::Diagonal` is equivalent to passing `W.diag`. + +`(a, b)` must cover `A`; use [`iscover`](@ref)`(a, b, A)` to check it. + +The methods accepting [`SupportComponents`](@ref) reuse a previous +[`support_components`](@ref)`(A)` computation. For these methods, `sc.rowax` and +`sc.colax` replace `axes(A, 1)` and `axes(A, 2)` in the axis requirements. + +Equivalent componentwise rescalings of `(a, b)` produce the same `s`. Some +coupling patterns make this impossible; the general form then throws an +`ArgumentError`. Set `degenerate=:uniform` to return a gauge-dependent cover. # Extended help @@ -51,45 +49,37 @@ row, so for `W` diagonal the sum needed is exactly the one over the rows of chain of such couplings form a group. Within a group, writing `M[p,q] = Σ_{i ∈ rows(p), i' ∈ rows(q)} a[i]*abs(W[i,i'])*a[i']` for the block sum over components `p` and `q`, `abs(G[j,k]) <= M[p,q]*b[j]*b[k]` for `j ∈ p` -and `k ∈ q`. Take `M` symmetric, as it is whenever `abs.(W)` is; the general -case needs one substitution, made in the remark below. Then - - s[j] = sqrt(Σ_q M[p,q] * sqrt(M[p,p]/M[q,q])) * b[j], j ∈ p - -meets every one of those bounds: the `q` term of the sum for `s[j]` and the `p` -term of the one for `s[k]` already multiply to `M[p,q]*M[q,p] = M[p,q]^2`, and -no term is negative. For a component that no `W` entry couples to another, the -group is a single component and this reduces to the diagonal-`W` formula above. -`G[j,k]` is exactly zero across distinct groups, and columns with no support get -`s[j] = 0`. - -Rescaling `a -> γ*a`, `b -> b/γ` within any support component — independently -per component — leaves `s` unchanged. Unlike `b` alone, `s` is therefore safe to -use as an absolute scale, e.g. a Levenberg-Marquardt damping term -`λ*Diagonal(s.^2)`: a caller should maintain the dimensionless `λ` against `s`, -not against a quantity that depends on which gauge the cover solver happened to -return. What pins the relative scale of coupled components is the ratio -`M[p,p]/M[q,q]`, extended to a component whose own block vanishes (`abs.(W)` zero -throughout it, though `W` couples it to a sibling) by propagating along the -coupling. A group in which *every* component's own block vanishes is the one -exception, and it is a genuine degeneracy rather than a shortcoming of the -formula: the gauge acts on the surviving off-diagonal data as -`M[p,q] -> γ[p]*γ[q]*M[p,q]`, which fixes each *product* `s[j]*s[k]` across two -components but nothing about how it divides between them. Such a group falls back -to the uniform total `sqrt(Σ_{p,q} M[p,q])`, which covers, but there `s` depends -on the factorization and not on the products alone. - -With more than one component and no coupling between them, `s[j] <= norm(a)*b[j]`, -strictly tighter whenever another component carries weight — the naive global -bound obtained by ignoring the block structure entirely. Coupled components admit -no such uniform comparison: fixing the gauge redistributes tightness among them, -resulting in an `s` that depends only on the products at the cost of individual entries -that a gauge-dependent global sum can beat. - -When `abs.(W)` is not symmetric, neither is `G`, and since `s[j]*s[k]` is a -single number bounding both `abs(G[j,k])` and `abs(G[k,j])`, `M[p,q]` is -replaced throughout by `max(M[p,q], M[q,p])`; nothing else changes, the gauge -included, since that replacement is itself symmetric. +and `k ∈ q`. When `abs.(W)` is not symmetric neither is `G`, and since +`s[j]*s[k]` is a single number bounding both `abs(G[j,k])` and `abs(G[k,j])`, the +block sums enter only through their symmetrization +`Ms[p,q] = max(M[p,q], M[q,p])`. What remains is to divide each bound between its +two components: any `σ` with `σ[p]*σ[q] >= Ms[p,q]` for every `p`, `q` yields + + s[j] = σ[p]*b[j], j ∈ p + +This is a symmetric cover of the `k×k` matrix `Ms`, where `k` is the number of +components in the group. The implementation computes +[`symcover_min`](@ref)`(AbsLog{2}(), Ms)`. For a component that no `W` entry +couples to another, this reduces to the diagonal-`W` formula. Entries of `G` +across distinct groups vanish, and unsupported columns get `s[j] = 0`. + +Rescaling `a -> γ*a`, `b -> b/γ` within each support component leaves `s` +unchanged. Under this rescaling, `Ms[p,q] -> γ[p]*γ[q]*Ms[p,q]` and the minimal +cover changes as `σ[p] -> γ[p]*σ[p]`; therefore `σ[p]*b[j]` is invariant. This +makes `s` suitable as an absolute scale, such as in a Levenberg-Marquardt term +`λ*Diagonal(s.^2)`. + +The invariant cover does not exist when a connected part of `Ms`'s support graph +has an edge, no loop, and is bipartite. Scaling one color class by `t` and the +other by `1/t` leaves `Ms` unchanged but changes the individual `σ` values. +`gramcover` throws an `ArgumentError` in this case. With +`degenerate=:uniform`, it instead uses +`σ[p] = sqrt(Σ_{p,q} Ms[p,q])`, which depends on the gauge of `(a, b)`. A loop +or odd cycle removes this degeneracy. + +For uncoupled components, `s[j] <= norm(a)*b[j]`, with a strict inequality when +another component carries weight. There is no uniform comparison for coupled +components because changing the gauge redistributes tightness among them. When a positive-semidefinite `W` is available only as an operator — `W[i,i]` readable, `W[i,i']` for `i != i'` not — `abs(W[i,i']) <= sqrt(W[i,i]*W[i',i'])` @@ -98,8 +88,11 @@ yields the looser diagonal-only bound from `diag(W)`. The methods here always compute the tighter entrywise form above, which requires `W`'s entries. -The accumulation feeding each `sqrt` is inflated to guarantee coverage despite -naive-summation roundoff, without ever forming `G` to check it. +The minimal-cover solve has size `k`, the number of support components coupled +by `W`. Constructing `Ms` also costs `O(k^2)` and usually dominates the solve. + +Roundoff margins on the block sums and a final feasibility check preserve the +cover in floating-point arithmetic. See also: [`gramcover!`](@ref), [`symcover`](@ref), [`cover`](@ref), [`iscover`](@ref). @@ -129,11 +122,12 @@ gramcover(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) = gramcover(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, w::AbstractVector) = gramcover(a, b, support_components(A), w) -gramcover(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, W::Diagonal) = - gramcover(a, b, A, W.diag) +# Accept and validate `degenerate` consistently for every matrix weight. +gramcover(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, W::Diagonal; degenerate::Symbol=GRAMCOVER_DEGENERATE) = + (_gc_check_degenerate(degenerate); gramcover(a, b, A, W.diag)) -gramcover(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, W::AbstractMatrix) = - gramcover(a, b, support_components(A), W) +gramcover(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, W::AbstractMatrix; kwargs...) = + gramcover(a, b, support_components(A), W; kwargs...) function gramcover(a::AbstractVector, b::AbstractVector, sc::SupportComponents) T = _gc_eltype(a, b) @@ -147,28 +141,29 @@ function gramcover(a::AbstractVector, b::AbstractVector, sc::SupportComponents, return gramcover!(s, a, b, sc, w) end -gramcover(a::AbstractVector, b::AbstractVector, sc::SupportComponents, W::Diagonal) = - gramcover(a, b, sc, W.diag) +gramcover(a::AbstractVector, b::AbstractVector, sc::SupportComponents, W::Diagonal; degenerate::Symbol=GRAMCOVER_DEGENERATE) = + (_gc_check_degenerate(degenerate); gramcover(a, b, sc, W.diag)) -function gramcover(a::AbstractVector, b::AbstractVector, sc::SupportComponents, W::AbstractMatrix) +function gramcover(a::AbstractVector, b::AbstractVector, sc::SupportComponents, W::AbstractMatrix; kwargs...) T = _gc_eltype(a, b, W) s = similar(Array{T}, sc.colax) - return gramcover!(s, a, b, sc, W) + return gramcover!(s, a, b, sc, W; kwargs...) end """ s = gramcover!(s, a, b, A) s = gramcover!(s, a, b, A, w::AbstractVector) - s = gramcover!(s, a, b, A, W::AbstractMatrix) + s = gramcover!(s, a, b, A, W::AbstractMatrix; degenerate=:$(GRAMCOVER_DEGENERATE)) s = gramcover!(s, a, b, sc::SupportComponents) s = gramcover!(s, a, b, sc::SupportComponents, w::AbstractVector) - s = gramcover!(s, a, b, sc::SupportComponents, W::AbstractMatrix) + s = gramcover!(s, a, b, sc::SupportComponents, W::AbstractMatrix; degenerate=:$(GRAMCOVER_DEGENERATE)) Mutating counterpart of [`gramcover`](@ref): writes the symmetric cover of the (weighted) Gram matrix into `s` and returns it, rather than allocating a new vector. `eachindex(s)` must match `axes(A, 2)` — `sc.colax` for the [`SupportComponents`](@ref) forms — in addition to the axis requirements -[`gramcover`](@ref) places on `a`, `b`, and `w`/`W`. +[`gramcover`](@ref) places on `a`, `b`, and `w`/`W`, and shares its +`degenerate` keyword. See also: [`gramcover`](@ref). """ @@ -178,11 +173,11 @@ gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, A::AbstractM gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, A::AbstractMatrix, w::AbstractVector) = gramcover!(s, a, b, support_components(A), w) -gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, A::AbstractMatrix, W::Diagonal) = - gramcover!(s, a, b, A, W.diag) +gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, A::AbstractMatrix, W::Diagonal; degenerate::Symbol=GRAMCOVER_DEGENERATE) = + (_gc_check_degenerate(degenerate); gramcover!(s, a, b, A, W.diag)) -gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, A::AbstractMatrix, W::AbstractMatrix) = - gramcover!(s, a, b, support_components(A), W) +gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, A::AbstractMatrix, W::AbstractMatrix; kwargs...) = + gramcover!(s, a, b, support_components(A), W; kwargs...) function gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc::SupportComponents) _check_gramcover_ab(a, b, sc) @@ -216,10 +211,11 @@ function gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc: return _write_gramcover!(s, b, sc.colcomp, oc, m, n) end -gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc::SupportComponents, W::Diagonal) = - gramcover!(s, a, b, sc, W.diag) +gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc::SupportComponents, W::Diagonal; degenerate::Symbol=GRAMCOVER_DEGENERATE) = + (_gc_check_degenerate(degenerate); gramcover!(s, a, b, sc, W.diag)) -function gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc::SupportComponents, W::AbstractMatrix) +function gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc::SupportComponents, W::AbstractMatrix; degenerate::Symbol=GRAMCOVER_DEGENERATE) + _gc_check_degenerate(degenerate) _check_gramcover_ab(a, b, sc) _check_gramcover_s(s, sc) axes(W) == (sc.rowax, sc.rowax) || @@ -284,7 +280,7 @@ function gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc: # Block sums `M[r][p,q] = Σ_{i ∈ comp p, i' ∈ comp q} a[i]*abs(W[i,i'])*a[i']` # over the components merged into root `r`; `nterm` counts the terms of each, - # for the roundoff inflation in `_gc_group_scales!`. + # for the roundoff inflation in `_gc_group_scales`. T = typeof(_gc_term(a, W)) M = Vector{Matrix{T}}(undef, ncomp) nterm = Vector{Matrix{Int}}(undef, ncomp) @@ -309,7 +305,7 @@ function gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc: end end - sq = _gc_group_scales!(members, M, nterm) + sq = _gc_group_scales(members, M, nterm, degenerate) oc = first(sc.colax) - 1 return _write_gramcover_sq!(s, b, sc.colcomp, oc, sq) end @@ -326,26 +322,27 @@ function _check_gramcover_ab(a::AbstractVector, b::AbstractVector, sc::SupportCo return nothing end +function _gc_check_degenerate(degenerate::Symbol) + degenerate in (:error, :uniform) || + throw(ArgumentError("`degenerate` must be :error or :uniform, got :$degenerate")) + return degenerate +end + function _check_gramcover_s(s::AbstractVector, sc::SupportComponents) eachindex(s) == sc.colax || throw(DimensionMismatch("`s` holds one Gram scale per support column: eachindex(s) must be $(string(sc.colax)), got $(string(eachindex(s)))")) return nothing end -# The per-term type of the sum accumulated into `m[c]`, computed from `zero` -# values so it tracks the units and precision of `a` (and `w`/`W`) without -# hardcoding a numeric type. +# Infer the accumulator type while preserving units and precision. _gc_term(a::AbstractVector) = zero(eltype(a)) * zero(eltype(a)) _gc_term(a::AbstractVector, w::AbstractVector) = abs(zero(eltype(w))) * zero(eltype(a)) * zero(eltype(a)) _gc_term(a::AbstractVector, W::AbstractMatrix) = zero(eltype(a)) * abs(zero(eltype(W))) * zero(eltype(a)) -# The element type `gramcover` allocates `s` at: the type of `sqrt` of a -# `_gc_term` times a `b`-scale, matching what `_write_gramcover!` actually computes. +# Element type of the output scale. _gc_eltype(a, b, args...) = typeof(sqrt(_gc_term(a, args...)) * zero(eltype(b))) -# Turn per-component sums `m[c]`, each accumulated from `n[c]` nonnegative terms, -# into per-component scales: the uncoupled case, where component `c`'s cover is -# just `sqrt(m[c])`. +# Convert uncoupled component sums to scales. # # Naive summation of `n` nonnegative terms computes `fl(Σ) >= Σ/(1+γ)` with # `γ = n*ulp/(1-n*ulp)`, `ulp = eps(scalarT)/2`; inflating `sqrt(m[c])` by @@ -358,38 +355,9 @@ function _write_gramcover!(s::AbstractVector, b::AbstractVector, colcomp::Vector return _write_gramcover_sq!(s, b, colcomp, oc, sq) end -# Per-component scales `sq[c]` from the block sums `M[r][p,q]` of each merged -# group: `s[j] = sq[c]*b[j]` covers `A'*W*A` iff `sq[p]*sq[q] >= max(M[p,q], M[q,p])` -# for every pair of components `p`, `q` in one group (the `max` because `W` need -# not be symmetric, so the `(j,k)` and `(k,j)` entries of the Gram matrix are -# bounded by different block sums). -# -# Writing `d[p] = M[p,p]`, the choice -# -# sq[p] = sqrt(Σ_q max(M[p,q], M[q,p]) * sqrt(d[p]/d[q])) -# -# satisfies that: the `q` term of `sq[p]^2` and the `p` term of `sq[q]^2` alone -# multiply to `max(M[p,q], M[q,p])^2`, and every term is nonnegative. -# -# The proof needs only that `d` be strictly positive, so `d` is purely a choice of -# gauge, and it is the choice that makes `s` independent of the gauge `a -> γ*a`, -# `b -> b/γ` applied independently per component. Under that gauge -# `M[p,q] -> γ[p]*γ[q]*M[p,q]`, so any `d` with `d[p] -> γ[p]^2*d[p]` makes each -# term `M[p,q]*sqrt(d[p]/d[q])` scale by `γ[p]^2`, leaving `sq[p]*b[j]` unchanged. -# A uniform group total `sqrt(Σ_{p,q} M[p,q])`, which mixes blocks that scale by -# different powers of `γ`, does not have that property. -# -# `d[p] = M[p,p]` transforms that way, and is the choice wherever it is nonzero. -# Where it is zero — `abs.(W)` vanishing throughout a component even though `W` -# couples it to a sibling — `d[p] = M[p,q]^2/d[q]` for an already-assigned -# neighbor `q` transforms the same way, and `M[p,q] > 0` for every coupled pair, -# so propagating outward from the nonzero diagonals covers the whole group. -# -# A group whose diagonal vanishes entirely leaves no pivot to propagate from, and -# no choice of `d` would help: the gauge sends `M[p,q] -> γ[p]*γ[q]*M[p,q]`, so -# with only off-diagonal blocks surviving it fixes the products `sq[p]*sq[q]` but -# not the split between them. That group falls back to the uniform total. -function _gc_group_scales!(members::Vector{Vector{Int}}, M::Vector{<:Matrix}, nterm::Vector{Matrix{Int}}) +# Compute a minimal symmetric cover of each symmetrized block matrix. Minimality +# makes the result invariant under componentwise rescaling of `(a, b)`. +function _gc_group_scales(members::Vector{Vector{Int}}, M::Vector{<:Matrix}, nterm::Vector{Matrix{Int}}, degenerate::Symbol) T = eltype(eltype(M)) scalarT = scalar_type(T) sq = [sqrt(zero(T)) for _ in 1:length(members)] @@ -402,54 +370,76 @@ function _gc_group_scales!(members::Vector{Vector{Int}}, M::Vector{<:Matrix}, nt sq[mems[1]] = sqrt(Mr[1, 1]) * (1 + (nr[1, 1] + 3) * eps(scalarT)) continue end - Ms = [max(Mr[p, q], Mr[q, p]) for p in 1:k, q in 1:k] # `W` need not be symmetric - d = _gc_gauge(Ms) - if d === nothing - v = sqrt(sum(Mr)) * (1 + (sum(nr) + 3) * eps(scalarT)) + # Account for summation and multiplication roundoff in each block. + Ms = [max(Mr[p, q], Mr[q, p]) * (1 + (max(nr[p, q], nr[q, p]) + 1) * eps(scalarT)) + for p in 1:k, q in 1:k] + if _loopless_bipartite(Ms) + degenerate === :uniform || throw(ArgumentError( + "gramcover: $k coupled support components have a loopless bipartite support graph, so no gauge-invariant cover exists; pass `degenerate=:uniform` to allow a gauge-dependent cover")) + v = sqrt(sum(Ms)) * (1 + (k * k + 3) * eps(scalarT)) for p in 1:k sq[mems[p]] = v end continue end + σ = _gc_inflate_to_cover(symcover_min(AbsLog{2}(), Ms), Ms, scalarT) for p in 1:k - acc = zero(T) - nt = 0 - for q in 1:k - acc += Ms[p, q] * sqrt(d[p] / d[q]) - nt += max(nr[p, q], nr[q, p]) - end - # `nt` terms of naive summation, plus a division and a `sqrt` per term, - # plus the outer `sqrt` and the multiply by `b[j]`. - sq[mems[p]] = sqrt(acc) * (1 + (nt + 3k + 3) * eps(scalarT)) + # The extra factor is for the multiply by `b[j]` in `_write_gramcover_sq!`. + sq[mems[p]] = σ[p] * (1 + 3 * eps(scalarT)) end end return sq end -# Strictly positive gauge `d` for one group's symmetrized block sums `Ms`: the -# diagonal where it is nonzero, propagated as `Ms[p,q]^2/d[q]` from an assigned -# neighbor `q` where it is not. Returns `nothing` if the diagonal is entirely zero. -# The propagation takes the lowest-indexed assigned neighbor, so `d` depends on -# `Ms` and the component ordering alone — never on the gauge it is fixing. -function _gc_gauge(Ms::Matrix) +# Uniformly inflate `σ` until it covers `Ms` in floating-point arithmetic. The +# retry guards against a downward-rounded inflation factor. +function _gc_inflate_to_cover(σ::AbstractVector, Ms::AbstractMatrix, ::Type{scalarT}) where {scalarT} + for _ in 1:8 + ρ = one(scalarT) + for q in axes(Ms, 2), p in axes(Ms, 1) + iszero(Ms[p, q]) && continue + pq = σ[p] * σ[q] + iszero(pq) && + throw(ErrorException("gramcover: zero scale on a coupled component; please report this with the inputs")) + ρ = max(ρ, scalarT(Ms[p, q] / pq)) + end + ρ <= 1 && return σ + σ = σ .* (sqrt(ρ) * (1 + 4 * eps(scalarT))) + end + throw(ErrorException("gramcover: group cover did not reach feasibility; please report this with the inputs")) +end + +# True when a nontrivial connected component of `Ms`'s support graph is +# loopless and bipartite. +# +# Zero-weight components are isolated and do not create a degeneracy. +function _loopless_bipartite(Ms::AbstractMatrix) k = size(Ms, 1) - d = [Ms[p, p] for p in 1:k] - any(!iszero, d) || return nothing - while any(iszero, d) - progress = false - for p in 1:k - iszero(d[p]) || continue - q = findfirst(q -> !iszero(d[q]) && !iszero(Ms[p, q]), 1:k) - q === nothing && continue - d[p] = Ms[p, q] * Ms[p, q] / d[q] - progress = true + color = zeros(Int8, k) + stack = Int[] + for root in 1:k + iszero(color[root]) || continue + color[root] = 1 + push!(stack, root) + hasedge = hasloop = false + bipartite = true + while !isempty(stack) + p = pop!(stack) + iszero(Ms[p, p]) || (hasloop = true) + for q in 1:k + (q == p || iszero(Ms[p, q])) && continue + hasedge = true + if iszero(color[q]) + color[q] = -color[p] + push!(stack, q) + elseif color[q] == color[p] + bipartite = false + end + end end - # Every pair the union-find merged has `Ms[p,q] > 0`, so the group is - # connected and a sweep that assigns nothing cannot happen; guard anyway - # rather than spin. - progress || error("gramcover: group with $(count(iszero, d)) unreachable component(s); please report this with the inputs") + hasedge && !hasloop && bipartite && return true end - return d + return false end # Write `s[j] = sq[c]*b[j]` for `j` in component `c` (`colcomp[j - oc]`), and diff --git a/test/gram_covers.jl b/test/gram_covers.jl index 864fe01..823aeeb 100644 --- a/test/gram_covers.jl +++ b/test/gram_covers.jl @@ -3,6 +3,34 @@ @testset "gramcover" begin + # Explicit symmetrized block sums covered by `gramcover`. + function blocksums(a, J, W) + sc = MatrixCovers.support_components(J) + k = MatrixCovers.ncomponents(sc) + M = zeros(k, k) + for i in axes(W, 1), ip in axes(W, 2) + ci, cip = MatrixCovers.rowcomponent(sc, i), MatrixCovers.rowcomponent(sc, ip) + (iszero(ci) || iszero(cip)) && continue + M[ci, cip] += a[i] * abs(W[i, ip]) * a[ip] + end + return [max(M[p, q], M[q, p]) for p in 1:k, q in 1:k] + end + + # Recover `σ[p]` from any supported column in component `p`. + function groupscales(s, b, J) + sc = MatrixCovers.support_components(J) + return [(j = findfirst(==(p), sc.colcomp); s[j] / b[j]) + for p in 1:MatrixCovers.ncomponents(sc)] + end + + # Three square blocks on the diagonal, one support component each. + function threeblocks(rng) + b1, b2, b3 = randn(rng, 2, 2), randn(rng, 2, 2), randn(rng, 2, 2) + return [b1 zeros(2, 2) zeros(2, 2) + zeros(2, 2) b2 zeros(2, 2) + zeros(2, 2) zeros(2, 2) b3] + end + @testset "random dense J: exact coverage of A'A" begin rng = StableRNG(3) J = randn(rng, 8, 5) @@ -87,8 +115,7 @@ @test isapprox(gramcover(a, b, J, W), gramcover(a2, b2, J, W); rtol=1e-12) end - # `W` with a vanishing diagonal block on one component: the gauge is - # propagated across the coupling rather than read off that block. + # A coupling fixes the scale of a component whose diagonal block vanishes. Wz = zeros(m, m) Wz[1:4, 1:4] .= 1.0 Wz[1, 5] = Wz[5, 1] = 2.0 @@ -135,18 +162,16 @@ Gc = J' * Wc * J @test all(sc * sc' .>= abs.(Gc) .- 1e-9 * maximum(abs, Gc)) - # Coupled components whose own blocks all vanish: the surviving data fixes - # the products `s[j]*s[k]` across the two but not the split between them, - # so `gramcover` covers, but `s` is not fixed by the products alone. Documented, - # not a bug. + # A single loopless edge admits no gauge-invariant cover. Wo = zeros(m, m) Wo[1, 5] = Wo[5, 1] = 2.0 - so = gramcover(a, b, J, Wo) + @test_throws "no gauge-invariant cover exists" gramcover(a, b, J, Wo) + so = gramcover(a, b, J, Wo; degenerate=:uniform) Go = J' * Wo * J @test all(so * so' .>= abs.(Go) .- 1e-9 * maximum(abs, Go)) a2, b2 = copy(a), copy(b) a2[1:4] .*= 8; b2[1:3] ./= 8 - @test !isapprox(so, gramcover(a2, b2, J, Wo); rtol=1e-6) + @test !isapprox(so, gramcover(a2, b2, J, Wo; degenerate=:uniform); rtol=1e-6) end @testset "sparse J" begin @@ -238,6 +263,16 @@ @test sbuf == gramcover(a, b, J, w) @test gramcover!(sbuf, a, b, sc, Diagonal(w)) === sbuf @test sbuf == gramcover(a, b, sc, w) + + # Matrix-weight methods accept and validate `degenerate` consistently. + @test gramcover(a, b, J, Diagonal(w); degenerate=:uniform) == gramcover(a, b, J, w) + @test gramcover(a, b, sc, Diagonal(w); degenerate=:uniform) == gramcover(a, b, sc, w) + @test gramcover!(sbuf, a, b, J, Diagonal(w); degenerate=:uniform) === sbuf + @test sbuf == gramcover(a, b, J, w) + @test_throws "`degenerate` must be :error or :uniform" gramcover(a, b, J, Diagonal(w); degenerate=:nonsense) + @test_throws "`degenerate` must be :error or :uniform" gramcover(a, b, sc, Diagonal(w); degenerate=:nonsense) + @test_throws "`degenerate` must be :error or :uniform" gramcover!(sbuf, a, b, J, Diagonal(w); degenerate=:nonsense) + @test_throws "`degenerate` must be :error or :uniform" gramcover!(sbuf, a, b, sc, Diagonal(w); degenerate=:nonsense) end @testset "general W leaves an uncoupled component as its own group" begin @@ -259,4 +294,134 @@ @test s[5:6] ≈ gramcover(a, b, J)[5:6] end + @testset "two coupled components, both diagonal blocks nonzero" begin + # For a 2x2 `Ms` with positive diagonal, the minimal cover is + # `σ[p] = sqrt(Ms[p,p]) * sqrt(max(1, κ))`, `κ = Ms[1,2]/sqrt(Ms[1,1]*Ms[2,2])`: + # below `κ = 1` the diagonal constraints bind alone, above it the coupling does. + rng = StableRNG(9) + B = randn(rng, 4, 3); C = randn(rng, 3, 2) + J = [B zeros(4, 2); zeros(3, 3) C] + a, b = cover(J) + m = size(J, 1) + κs = Float64[] + for c in (0.05, 40.0) + W = Matrix{Float64}(I, m, m) + W[1, 5] = W[5, 1] = c + s = gramcover(a, b, J, W) + Ms = blocksums(a, J, W) + κ = Ms[1, 2] / sqrt(Ms[1, 1] * Ms[2, 2]) + push!(κs, κ) + σ = groupscales(s, b, J) + @test σ ≈ sqrt.([Ms[1, 1], Ms[2, 2]]) .* sqrt(max(1, κ)) rtol = 1e-6 + G = J' * W * J + @test all(s * s' .>= abs.(G)) + end + @test κs[1] < 1 < κs[2] + end + + @testset "two coupled components, one diagonal block zero" begin + # For `Ms = [d e; e 0]`, the minimal cover saturates its two constraints. + rng = StableRNG(11) + B = randn(rng, 4, 3); C = randn(rng, 3, 2) + J = [B zeros(4, 2); zeros(3, 3) C] + a, b = cover(J) + m = size(J, 1) + W = zeros(m, m) + W[1:4, 1:4] .= 1.0 + W[1, 5] = W[5, 1] = 2.0 + s = gramcover(a, b, J, W) + Ms = blocksums(a, J, W) + @test iszero(Ms[2, 2]) + σ = groupscales(s, b, J) + @test σ[1] ≈ sqrt(Ms[1, 1]) rtol = 1e-6 + @test σ[2] ≈ Ms[1, 2] / sqrt(Ms[1, 1]) rtol = 1e-6 + G = J' * W * J + @test all(s * s' .>= abs.(G)) + end + + @testset "three coupled components in a triangle, every diagonal block zero" begin + # An odd cycle fixes the gauge without a diagonal block. + rng = StableRNG(13) + J = threeblocks(rng) + a, b = cover(J) + m = size(J, 1) + W = zeros(m, m) + W[1, 3] = W[3, 1] = 1.5 + W[1, 5] = W[5, 1] = 0.75 + W[3, 5] = W[5, 3] = 2.25 + s = gramcover(a, b, J, W) + Ms = blocksums(a, J, W) + @test all(iszero, [Ms[p, p] for p in 1:3]) + σ = groupscales(s, b, J) + @test σ ≈ [sqrt(Ms[1, 2] * Ms[1, 3] / Ms[2, 3]), + sqrt(Ms[1, 2] * Ms[2, 3] / Ms[1, 3]), + sqrt(Ms[1, 3] * Ms[2, 3] / Ms[1, 2])] rtol = 1e-6 + G = J' * W * J + @test all(s * s' .>= abs.(G)) + + # The result remains gauge-invariant. + a2, b2 = copy(a), copy(b) + for (p, (rows, cols)) in enumerate(((1:2, 1:2), (3:4, 3:4), (5:6, 5:6))) + γ = (2.0, 0.3, 7.0)[p] + a2[rows] .*= γ; b2[cols] ./= γ + end + @test isapprox(s, gramcover(a2, b2, J, W); rtol=1e-6) + end + + @testset "loopless bipartite coupling graphs are refused" begin + # An even cycle without loops is bipartite. + rng = StableRNG(17) + blocks = [randn(rng, 2, 2) for _ in 1:4] + J = zeros(8, 8) + for p in 1:4 + J[2p-1:2p, 2p-1:2p] .= blocks[p] + end + a, b = cover(J) + @test MatrixCovers.ncomponents(MatrixCovers.support_components(J)) == 4 + W = zeros(8, 8) + for (i, ip) in ((1, 3), (3, 5), (5, 7), (7, 1)) + W[i, ip] = W[ip, i] = 1.0 + 0.5 * i + end + @test_throws "no gauge-invariant cover exists" gramcover(a, b, J, W) + s = gramcover(a, b, J, W; degenerate=:uniform) + G = J' * W * J + @test all(s * s' .>= abs.(G)) + + @test_throws "`degenerate` must be :error or :uniform" gramcover(a, b, J, W; degenerate=:nonsense) + + # One loop removes the obstruction. + W[1, 1] = 1.0 + s1 = gramcover(a, b, J, W) + G1 = J' * W * J + @test all(s1 * s1' .>= abs.(G1)) + end + + @testset "generic W over three coupled components" begin + # Compare the minimal cover with a valid diagonal-normalized cover. + rng = StableRNG(19) + J = threeblocks(rng) + a, b = cover(J) + m = size(J, 1) + W = abs.(randn(rng, m, m)) .+ 0.1 + s = gramcover(a, b, J, W) + G = J' * W * J + @test all(s * s' .>= abs.(G)) + + Ms = blocksums(a, J, W) + σ = groupscales(s, b, J) + @test σ ≈ symcover_min(AbsLog{2}(), Ms) rtol = 1e-6 + + # This diagonal-normalized formula covers `Ms` but need not be minimal. + d = [Ms[p, p] for p in 1:3] + σh = [sqrt(sum(Ms[p, q] * sqrt(d[p] / d[q]) for q in 1:3)) for p in 1:3] + @test all(σh * σh' .>= Ms) + @test cover_objective(AbsLog{2}(), σ, Ms) <= cover_objective(AbsLog{2}(), σh, Ms) + + a2, b2 = copy(a), copy(b) + for (p, γ) in enumerate((5.0, 0.2, 1.7)) + a2[2p-1:2p] .*= γ; b2[2p-1:2p] ./= γ + end + @test isapprox(s, gramcover(a2, b2, J, W); rtol=1e-6) + end + end diff --git a/test/unitful.jl b/test/unitful.jl index bbc5c85..38d160e 100644 --- a/test/unitful.jl +++ b/test/unitful.jl @@ -294,6 +294,29 @@ @test all(ustrip.(sw * sw') .>= ustrip.(abs.(Gw))) end + @testset "gramcover with a general W carries units" begin + # Coupled block sums have units `a^2*W`. + J = [4.0 1.0 0.0 0.0 + 1.0 3.0 0.0 0.0 + 0.0 0.0 2.0 1.0 + 0.0 0.0 0.5 4.0] .* u"N/m" + a, b = cover(J) + W = [1.0 0.0 0.3 0.0 + 0.0 2.0 0.0 0.0 + 0.3 0.0 1.5 0.0 + 0.0 0.0 0.0 0.5] .* u"s" + s = gramcover(a, b, J, W) + @test unit.(s) == unit(a[1]) .* unit.(b) .* unit(sqrt(1.0u"s")) + G = J' * (W * J) + @test all(ustrip.(s * s') .>= ustrip.(abs.(G))) + + # Units do not alter gauge invariance. + a2, b2 = copy(a), copy(b) + a2[1:2] .*= 3; b2[1:2] ./= 3 + a2[3:4] .*= 0.4; b2[3:4] ./= 0.4 + @test isapprox(ustrip.(s), ustrip.(gramcover(a2, b2, J, W)); rtol=1e-6) + end + @testset "unit types the cover cannot use" begin # An affine unit measures from a shifted origin, so no product a[i]*b[j] # scales it: refuse it rather than silently reducing it to its atom. From c04cb4cca7265c54eaf0c7c49fb1ae7fe743d38b Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 21 Aug 2026 16:02:39 -0500 Subject: [PATCH 2/2] Tighten the text & correct stale comments --- README.md | 3 +- docs/src/index.md | 180 ++++++++++-------------- ext/MatrixCoversIpoptExt.jl | 7 +- ext/MatrixCoversJuMPExt.jl | 15 +- src/heuristic_covers.jl | 39 ++---- src/initializers.jl | 30 ++-- src/minimal_covers.jl | 70 +++------ src/penalties.jl | 10 +- src/soft_covers.jl | 273 +++++++++++------------------------- test/heuristic_covers.jl | 8 +- test/initializers.jl | 6 +- test/minimal_covers.jl | 4 +- 12 files changed, 216 insertions(+), 429 deletions(-) diff --git a/README.md b/README.md index 77f4fb3..c57edc4 100644 --- a/README.md +++ b/README.md @@ -8,8 +8,7 @@ This package computes **covers** of matrices: non-negative vectors `a` and `b` such that `a[i] * b[j] >= abs(A[i, j])` for all `i`, `j`. Covers are the -natural scale-covariant representation of a matrix — under row/column diagonal -scaling they transform exactly as the matrix entries do — making them a useful +natural scale-covariant representation of a matrix, making them a useful building block for scale-invariant numerical analysis. In particular, $`\hat A = A ./ (a b^T)`$ is scale-invariant, and because $`|\hat A[i, j]| \le 1`$ for all `i` and `j`, this simple construct finds applications that range from diff --git a/docs/src/index.md b/docs/src/index.md index 0cb24cd..07ae7d3 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -15,30 +15,26 @@ C_{ij} \;\geq\; |A_{ij}| \quad \text{for all } i, j. For a symmetric matrix the cover is symmetric (`b = a`), so a single vector suffices: `a[i] * a[j] >= abs(A[i, j])`. -An *optimal cover* is one for which `C` is "as tight as possible" in bounding -`A` (equivalently, no larger than strictly necessary), by criteria that will be -described below. +A *minimal cover* chooses `C` "as tight as possible" in bounding `A`, by +criteria that will be described below. -This package also supports *soft covers*: these satisfy -`C[i, j] ⪆ abs(A[i, j])`, meaning that `C` matches or exceeds `A` at most indexes -but not necessarily all; this intuitive notion will be made concrete -below. +The package also supports *soft covers*, which penalize uncovered entries instead +of requiring every inequality to hold. ## Why covers? Covers provide a natural **scale-covariant** "summary" of a matrix. If you rescale rows by a positive diagonal factor `D_r` and columns by `D_c`, the -optimal cover transforms as `a → D_r * a`, `b → D_c * b`, exactly mirroring -how the matrix entries change. Scalar summaries like `norm(A)` or +optimal cover transforms as `a → D_r * a`, `b → D_c * b`, so the product `a * b'` +transforms identically to `A`. Scalar summaries like `norm(A)` or `maximum(abs, A)` do not have this property and therefore implicitly encode an arbitrary choice of units. -Most users will employ matrices that store pure numbers, and this package works -well with such matrices. But to emphasize the scale-covariance, we'll start with -an example of a 3×3 matrix whose rows and columns correspond to physical -variables with different units — position in meters, velocity in m/s, force in +While most users will employ matrices that store pure numbers, we'll start with +an example of a 3×3 matrix whose rows and columns correspond to *physical +variables with different units* — position in meters, velocity in m/s, force in Newtons. Loading [Unitful](https://github.com/PainterQubits/Unitful.jl) lets -the matrix carry those units itself: +the matrix express those units directly: ```jldoctest coverunits julia> using MatrixCovers, Unitful @@ -58,12 +54,9 @@ julia> round.(typeof.(a), a; digits=6) 0.001 N^-1 ``` -`A[i,j]` has units `1/(u[i]*u[j])` (modeling a Hessian matrix for functions of -parameter vectors with units `u[i]`), so `a[i]` comes back with gradient-like -units of `1/u[i]`: the cover provides each variable's natural scale inferred -from `A`, here 1 mm, 1 m/s, and 1 kN. Had we expressed the original matrix in -those units directly, we would have gotten the equivalent cover stated in those -units. +`A[i,j]` has units `1/(u[i]*u[j])`, as in a Hessian whose parameters have units +`u[i]`. The cover has units `1/u[i]` and identifies scales of 1 mm, 1 m/s, and +1 kN. Had we expressed `A` in those units, we would have gotten the equivalent cover. Normalizing by the cover cancels the units along with the magnitudes, leaving a matrix that is all-ones, dimensionless, and scale-invariant: @@ -76,13 +69,13 @@ julia> round.(A ./ (a .* a'); digits=6) 1.0 1.0 1.0 ``` -It is worth noting that this yields 1 only for entries where the cover bound is -*tight*; had `A` been, say, diagonal, then `A ./ (a .* a')` would also be diagonal. +An entry is 1 only where the cover bound is tight, and this is not guaranteed for all matrices. +For example, given diagonal `A`, the normalized matrix is also diagonal. A cover exists only when the units of `A` factor as `unit(A[i,j]) == unit(a[i])*unit(b[j])`, and a matrix that fails this is rejected with a -`DimensionMismatch`. This requirement is not exhorbitant: without it, the terms -in a row of `A*x` carry different units, so `A*x` is undefined for every `x`. +`DimensionMismatch`. Without it, the terms in a row of `A*x` have incommensurate +units and cannot be added, so `A*x` is undefined for every `x`. If a matrix can be used in matrix-vector multiplication, it has a cover. ## Penalty functions @@ -146,38 +139,26 @@ You can override the default penalty by supplying it as an argument to the solve | [`soft_symcover_min`](@ref) | yes | soft (penalized) | `AbsLog{2}`, `AbsLinear` | native for `AbsLog{2}`; else JuMP | | [`soft_cover_min`](@ref) | no | soft (penalized) | `AbsLog{2}`, `AbsLinear` | native for `AbsLog{2}`; else JuMP | -The two soft tiers are separated by a different axis than the two hard ones. For hard -covers, [`symcover`](@ref) trades optimality for speed while guaranteeing feasibility, and -[`symcover_min`](@ref) is optimal. Both soft tiers minimize the same unconstrained -objective, and differ instead in what they promise about reaching its minimum: - -- [`soft_symcover`](@ref) and [`soft_cover`](@ref) are **always native and best-effort.** - They require no extension for any penalty, and they own their multistart — but what they - return is a coordinate-descent fixed point, which for the non-convex and nonsmooth - penalties need not be a minimizer. -- [`soft_symcover_min`](@ref) and [`soft_cover_min`](@ref) return a **true minimizer of the - basin they start in, and may require an extension.** `AbsLog{2}` is native; the - `AbsLinear` penalties need JuMP and Ipopt; `AbsLog{1}` is not implemented. - -Both reduce to the same trade: cheap and always available, against best quality and -possibly an extra dependency. - -Under `AbsLog{2}` the objective is convex with a single minimizer, so the tiers coincide — -[`soft_symcover`](@ref) *is* [`soft_symcover_min`](@ref) there, and likewise for the -asymmetric pair. That is the degenerate case of the contract rather than an exception to -it: with one minimizer there is nothing for a best-effort descent and a minimizer to -disagree about. Under `AbsLog{1}` they part company most sharply — the soft `AbsLog{1}` -covers are coordinate descents that reach a deterministic fixed point rather than a -minimizer, and `soft_symcover_min`/`soft_cover_min` do not accept `AbsLog{1}` at all. - -**[`symcover`](@ref), [`cover`](@ref), and any native implementation can be recommended for production use,** -possibly with relaxed convergence bounds. -The heuristic solvers are particularly fast: they run in ``O(mn)`` time for an -``m\times n`` matrix and often land within a few percent of the -objective-minimal cover (see the quality tests involving -`test/testmatrices.jl`). Native solvers (both hard and soft) are intermediate, still roughly ``O(mn)`` but -requiring many iterations (and for non-convex cases, multiple start points by default) for convergence; -still, they are much faster than their JuMP-counterparts, which are provided mainly as a reference. +For hard covers, [`symcover`](@ref) trades optimality for speed while +[`symcover_min`](@ref) minimizes the selected objective. The soft solvers differ +in their convergence guarantees: + +- [`soft_symcover`](@ref) and [`soft_cover`](@ref) use native coordinate descent + and multistart. For nonconvex or nonsmooth penalties they may stop at a fixed + point that is not a local minimum. +- [`soft_symcover_min`](@ref) and [`soft_cover_min`](@ref) find a local minimum. + `AbsLog{2}` is native; `AbsLinear` requires JuMP and Ipopt; `AbsLog{1}` is not + implemented. + +Under `AbsLog{2}` the objective is convex with a unique minimizer, so both soft +tiers return the same result. Under `AbsLog{1}`, only the native coordinate-descent +solvers are available. + +The heuristic solvers run in ``O(mn)`` time for an ``m\times n`` matrix and +often come within a few percent of the minimum objective. Native iterative +solvers are also roughly ``O(mn)`` per iteration. JuMP-based methods are slower +and primarily useful when the native methods do not support the selected +penalty. ### Covariance of the heuristics @@ -201,10 +182,8 @@ julia> round.(extrema(P2 ./ P1); digits=3) # not for the heuristic sol (1.0, 1.077) ``` -The departure is small (bounded by the heuristic's own suboptimality) and both covers are valid, so it matters -mostly when the covariance itself is what you are relying on. When it is, use +Both covers are valid, but the difference matters when covariance is required. Use [`symcover_min`](@ref) or [`cover_min`](@ref), whose minimizer is scale-covariant -by construction. ### Objective-minimal covers @@ -239,9 +218,8 @@ julia> round(cover_objective(AbsLog{2}(), aq, bq, A); digits=6) 1.141281 ``` -The native solver is near-exact (relative objective excess typically a few -``\times 10^{-7}``, growing slowly with problem size) and orders of magnitude -faster than a general-purpose convex solver. The other penalties — `AbsLog{1}` +The native solver typically has relative objective error of a few +``\times 10^{-7}``, growing slowly with problem size. The other penalties — `AbsLog{1}` (a linear program) and the non-convex `AbsLinear` variants — are solved through [JuMP](https://jump.dev/) and are loaded on demand as a package extension: @@ -269,9 +247,8 @@ julia> round.(a * b'; digits=6) # tight on four of the six entries The solver returns values good to roughly solver tolerance, so these examples round before displaying. -The soft `*_min` solvers divide along the same line, but not at the same place: -[`soft_symcover_min`](@ref) and [`soft_cover_min`](@ref) solve `AbsLog{2}` natively and -reach for JuMP with Ipopt only for the `AbsLinear` penalties. They do not accept +[`soft_symcover_min`](@ref) and [`soft_cover_min`](@ref) solve `AbsLog{2}` +natively and use JuMP with Ipopt for `AbsLinear`. They do not accept `AbsLog{1}`; the soft `AbsLog{1}` covers are available through [`soft_symcover`](@ref) and [`soft_cover`](@ref), which are native. @@ -288,11 +265,9 @@ support graph of `A` (rows and columns as vertices, stored nonzeros as edges), so the sums are taken within each component separately. This convention is not scale-invariant but has no impact on the cover itself. -Other penalties may be more degenerate. The `AbsLog{1}()` penalty is identical -over a whole face of the feasible polytope, and its members are genuinely -different covers — the products `a[i]*b[j]` differ — that merely happen to score -the same objective. To make the result deterministic, we select the one that -additionally minimizes the `AbsLog{2}` objective. +The `AbsLog{1}()` minimum can be a face of the feasible polytope containing +different covers with the same objective. The implementation returns the member +that also minimizes the `AbsLog{2}` objective. `AbsLinear` penalties typically have isolated minima, so are not as degenerate as `AbsLog{1}()`, but these minima occur in separate basins. There is no @@ -300,41 +275,34 @@ guarantee of global optimality. ### Starting points: initialize and refine -For objectives with multiple minima, the solver starts from a specified point and descends. -At a lower level, this package's interface is organized in three layers: +For objectives with multiple minima, the result can depend on the starting +point. The interface has three layers: -- **Initializers** name the starting points. [`initialize_symcover`](@ref) and +- **Initializers** construct starting points. [`initialize_symcover`](@ref) and [`initialize_cover`](@ref) take a `strategy` — `:geomean`, `:leaveout`, `:diagfeasible`, or `:hardcover` — and return that point. Each is a property of `A` alone; no objective - is involved, so an initializer takes no penalty. A second keyword, `feasible`, says how - the point is brought up to covering `A`: `:inflate` (the default) scales it bodily by one + is involved, so an initializer takes no penalty. A second keyword, `feasible`, + controls how the point is made to cover `A`: `:inflate` (the default) applies one common factor, `:boost` raises only the rows touching a violated entry, and `:none` - leaves it as it is. The hard-cover solvers need a cover, so they take one of the first - two; the soft covers want `:none`, since forcing the geometric mean to cover `A` would - spoil the very property that makes it the soft `AbsLog{2}` optimum. - - The two feasible routes land on the boundary at different points, hence in different - basins — which is why the choice is a named part of the start rather than an internal - detail. The heuristic [`cover`](@ref) is itself a composition of these: the geometric - mean, boosted, then tightened. + leaves it unchanged. Hard-cover solvers use one of the first two; soft-cover + solvers use `:none`. + + The two feasible routes reach different points on the boundary and can enter + different basins. The heuristic [`cover`](@ref) computes the geometric mean, + boosts it, and then tightens it. - **Refiners** improve a starting point in place, and are the `!`-suffixed forms of the solvers: [`symcover_min!`](@ref), [`cover_min!`](@ref), [`soft_symcover!`](@ref), [`soft_cover!`](@ref), [`soft_symcover_min!`](@ref), and [`soft_cover_min!`](@ref) - validate the start, then optimize from it. Which basin they reach is the caller's - choice, by construction, and supplying the start is the caller's job. The hard refiners - require a start that covers `A`; the soft ones do not, since their objective constrains - nothing — build theirs with `feasible=:none`. + validate a caller-provided start, then optimize from it. Hard refiners require + a cover; soft refiners accept starts built with `feasible=:none`. - **Solvers** bundle the two. [`symcover_min`](@ref), [`cover_min`](@ref), [`soft_symcover`](@ref), [`soft_cover`](@ref), [`soft_symcover_min`](@ref), and - [`soft_cover_min`](@ref) refine every start on a menu (the `strategies` keyword, or the - multistart's own list) and return the best cover by [`cover_objective`](@ref), so their - result depends on `A` and not on an initialization the caller never chose. + [`soft_cover_min`](@ref) refine a set of starts (the `strategies` keyword or the + multistart list) and return the best result by [`cover_objective`](@ref). -That is the rule for the whole grid: **the plain form owns the menu, so its result is a -property of `A`; the `!` form refines the one start you give it, so its result is a -property of `A` and that start.** [`symcover!`](@ref) and [`cover!`](@ref) are the -exception that proves it — they are initializers, not refiners, and construct their cover -from scratch rather than reading the vector passed in. +In general, a plain form chooses among several starts, while a `!` form refines +the supplied start. [`symcover!`](@ref) and [`cover!`](@ref) are initializers and +overwrite their vector arguments. For finer control, you can run these manually: @@ -343,19 +311,19 @@ julia> using MatrixCovers, JuMP, Ipopt # Ipopt for the AbsLinear penalties julia> S = [4 1 0; 1 1 5; 0 5 2]; -julia> round.(symcover_min(AbsLinear{2}(), S); digits=6) # multistart over the whole menu +julia> round.(symcover_min(AbsLinear{2}(), S); digits=6) # use all default starts 3-element Vector{Float64}: 2.0 1.0 5.0 -julia> round.(symcover_min(AbsLinear{2}(), S; strategies=(:geomean,)); digits=6) # or commit to one start +julia> round.(symcover_min(AbsLinear{2}(), S; strategies=(:geomean,)); digits=6) # use one strategy 3-element Vector{Float64}: 2.0 1.0 5.0 -julia> a0 = initialize_symcover(S; strategy=:geomean); # or drive it yourself +julia> a0 = initialize_symcover(S; strategy=:geomean); # construct a start explicitly julia> symcover_min!(AbsLinear{2}(), a0, S); @@ -383,13 +351,11 @@ rescaling: if `a*b'` covers `A`, then `a*(D*b)'` covers `A*D` — but the balanced representative of the rescaled problem is `(γ*a, D*b/γ)` for a per-component constant `γ ≠ 1` that depends on `D`. -This has important implications for applications where you might estimate covers -by composition. Let's take the example of the -[Levenberg-Marquardt algorithm](https://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm), -where you form products `J'*J` of the Jacobian matrix `J`. Suppose `a*b'` is -a cover of `J`: then `(a'*a) * b * b'` is a cover of `J'*J` (note -`a'*a` is a scalar). The *tightness* of this cover for `J'*J` depends on the convention -used to balance `a` and `b`. +This matters when covers are composed. For example, the +[Levenberg-Marquardt algorithm](https://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm) +uses products `J'*J` of a Jacobian `J`. If `a*b'` covers `J`, then +`(a'*a) * b * b'` covers `J'*J`, but its tightness and invariance depends on the balance +convention for `a` and `b`. To do better, this package provides the Gram cover `s = ` [`gramcover`](@ref)`(a, b, J[, W])`, a symmetric cover of `J'*W*J` built from the asymmetric cover of `J`. @@ -418,8 +384,8 @@ true ## Worked example: roundoff in `A \ b` -Because a cover names each variable's natural scale, it also says how to measure a -solution in units that do not depend on how the problem was parameterized. +A cover provides scales for each variable with which to measure a solution independently +of the parameterization. Solving `x = A \ b` is *contravariant*: rescaling `A → D*A*D` and `b → D*b` sends `x → x ./ d`, while the cover is covariant, `a → d .* a`. The products `x .* a` are diff --git a/ext/MatrixCoversIpoptExt.jl b/ext/MatrixCoversIpoptExt.jl index 922fcf9..bf47724 100644 --- a/ext/MatrixCoversIpoptExt.jl +++ b/ext/MatrixCoversIpoptExt.jl @@ -18,11 +18,8 @@ using MatrixCovers: _edge_list, _sym_edge_list, _degrees # hard-cover models below sum over its `ei <= ej` half instead, per the objective each # one is defined by. -# The AbsLinear objectives are non-convex, so Ipopt returns a local minimum of -# whichever basin it descends into from the start it is given. That makes the start a -# genuine input rather than a hint, and it is why the hard-cover entry points here are the -# `*_min!` refiners, which take the start from the caller. The non-mutating `symcover_min` -# and `cover_min` are multistart drivers over these kernels and live in the main package. +# Ipopt returns a local minimum selected by the start. These kernels therefore +# implement the mutating refiners; the main package supplies multistart drivers. check_solved(model, fname) = MatrixCovers.check_solved(JuMP.termination_status(model), "Ipopt", fname) diff --git a/ext/MatrixCoversJuMPExt.jl b/ext/MatrixCoversJuMPExt.jl index 054f87d..05cffa3 100644 --- a/ext/MatrixCoversJuMPExt.jl +++ b/ext/MatrixCoversJuMPExt.jl @@ -62,19 +62,10 @@ end # bounds how far the reported AbsLog{1} objective can drift above its true optimum. const LEX_L1_SLACK = 1e-9 -# Second stage of the lexicographic AbsLog{1} solve, run on the just-optimized `model`. -# The AbsLog{1} optimum is a whole face of the feasible polytope rather than a point: its -# members are genuinely different covers — the products a[i]*a[j] differ — that happen to -# score the same objective, so the solver would otherwise return whichever vertex it landed -# on. Pinning the AbsLog{1} objective at its optimum and minimizing the AbsLog{2} objective -# over what remains selects one canonical member: the strictly convex quadratic has a unique -# minimizer over the face, and both objectives are functions of the residuals alone (which a -# rescaling A -> D*A*D leaves invariant), so the choice is scale-covariant. This is what makes -# the result independent of the start. +# Select a unique point on the optimal AbsLog{1} face by minimizing AbsLog{2} +# over it. Both objectives depend only on scale-invariant residuals. # -# `lin` is the AbsLog{1} objective and `residuals` the expressions α[i]+α[j]-log|A[i,j]| over -# the support, one per stored entry — the same convention `cover_objective` sums over, so the -# quadratic minimized here is the AbsLog{2} objective it reports. +# `residuals` uses the same support weighting as `cover_objective`. function _minimize_l2_over_l1_face!(model, lin, residuals, fname) isempty(residuals) && return nothing linopt = JuMP.value(lin) diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index dc5ee0a..b0c8a83 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -13,14 +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 initialization is the AbsLog{2} unconstrained minimum (geometric mean of -nonzero entries per row). It is then boosted to feasibility by a greedy -max-deficit rule (the most-violated entries are covered first), and `maxiter` -iterations of the tightening algorithm (Algorithm 1 of the manuscript) are -applied. - -`ϕ` names the penalty the caller would like the cover to do well on. Currently, -the heuristic covers ignore `ϕ`, although this behavior may change in future versions. +The method initializes from per-row geometric means, 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). See also: [`symcover!`](@ref), [`symcover_min`](@ref), [`soft_symcover`](@ref), [`cover`](@ref). @@ -77,11 +73,10 @@ end a, b = cover(ϕ, A; maxiter=3) a, b = cover(A; maxiter=3) -Given a matrix `A`, return vectors `a` and `b` such that `a[i] * b[j] >= abs(A[i, j])` -for all `i`, `j`. The initialization is the AbsLog{2} -unconstrained minimum (geometric mean of nonzero entries per row/column). It is -then boosted to feasibility by a greedy max-deficit rule (the most-violated -entries are covered first), and `maxiter` tightening iterations are applied. +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. Only the products `a[i] * b[j]` are determined by the problem: `a -> c*a`, `b -> b/c` leaves every one of them unchanged. The split is fixed by the balance convention @@ -89,8 +84,7 @@ leaves every one of them unchanged. The split is fixed by the balance convention column `j`), imposed within each connected component of the support (the gauge acts independently on each), as it is throughout the package; see [`cover_min`](@ref). -`ϕ` names the penalty the caller would like the cover to do well on. Currently, -the heuristic covers ignore `ϕ`, although this behavior may change in future versions. +`ϕ` is accepted for API compatibility but is currently ignored. For a cover that provably minimizes a given `ϕ`, use [`cover_min`](@ref). See also: [`cover!`](@ref), [`cover_min`](@ref), [`symcover`](@ref). @@ -149,12 +143,7 @@ function cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs. unconstrained_min!(AbsLog{2}(), a, b, A) boost_feasible!(a, b, A) tighten_cover!(a, b, A; kwargs...) - # The boost and the tightening raise and shrink rows and columns independently, so they - # leave the gauge wherever they happen to land it. Pin it, so the split reported between - # `a` and `b` is the package's convention rather than a residue of the passes above. - # The shift is exact in the log-scales but rounds in the products, which the tightening - # has driven onto the coverage boundary; the uniform inflation restores exact coverage, - # and preserves the balance just set because ∑ nzaᵢ = ∑ nzbⱼ. + # Apply the package's balance convention, then restore coverage lost to rounding. _balance_cover!(a, b, A) return inflate_feasible!(a, b, A) end @@ -237,7 +226,6 @@ end # ∑_{i,j: A[i,j]≠0} (log(a[i]*a[j]) - log|A[i,j]|)² # Fills `a` in-place and returns nza[i] = number of nonzero entries in row i. # For efficiency, uses a Sherman-Morrison approximation for the pattern of nonzeros. (It's exact when there are no zeros.) -# This is the "rank-1 solution" described in manuscript section 5.2. 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)))")) @@ -655,10 +643,9 @@ end # `foreach_support_sym`. Requires a start with strictly positive scale on every # supported row (the geometric-mean init from `unconstrained_min!` guarantees this). # -# Unlike `boost_feasible!`, which raises only the rows touching violated entries, -# this leaves the shape of the starting point untouched and moves it bodily to the -# feasibility boundary. The two reach different basins of the non-convex AbsLinear -# objective, which is why `soft_symcover` offers both as starts. The shift depends +# Unlike `boost_feasible!`, this preserves the shape of the starting point. The +# two methods can reach different basins of the nonconvex AbsLinear objective. +# The shift depends # on `A` only through the log-deficits at the starting point, which are invariant # under a diagonal rescaling `D*A*D`, so the result is scale-covariant. Growing the # log-scales directly (rather than multiplying by `exp(t)`) stays finite even when diff --git a/src/initializers.jl b/src/initializers.jl index 01a57fa..471cfa3 100644 --- a/src/initializers.jl +++ b/src/initializers.jl @@ -1,8 +1,4 @@ -# Named starting covers, shared by every algorithm that needs a starting point: -# the soft-cover multistarts and the `*_min` solvers alike. Each strategy is a -# deterministic, scale-covariant point built from the primitives in -# `heuristic_covers.jl`; nothing here calls a solver, so the dependency runs one -# way and the start menu has a single definition. +# Starting covers shared by the soft-cover multistarts and `*_min` solvers. # Starting covers the non-convex AbsLinear solvers refine, in the order they are tried. # `:leaveout` and `:diagfeasible` have no asymmetric formulation, so the two menus differ. @@ -21,39 +17,33 @@ const COVER_MIN_STRATEGIES = (:hardcover, :geomean) Build a starting point for the symmetric cover of `A`, as consumed by [`symcover_min`](@ref) and by the [`soft_symcover`](@ref) multistart. -No penalty is taken: every strategy below is a property of `A` alone, so the -starting point does not depend on the objective it will be refined against. +The strategies depend only on `A`, so this function takes no penalty. `strategy` names the point: -- `:geomean` — the geometric mean of the nonzero entries of each row, and *not* a - cover. It minimizes the soft AbsLog{2} objective exactly when every entry of `A` - is nonzero; on a sparse support it approximates that minimum, which - [`soft_symcover_min`](@ref)`(AbsLog{2}(), A)` returns exactly. +- `:geomean` — the geometric mean of each row's nonzero entries. It is not + generally a cover. - `:leaveout` — the geometric mean recomputed with the most-underweighted support entry dropped, which lands in the basin that treats that entry as effectively zero. Raises an `ArgumentError` when no entry can be dropped (empty support, or dropping it would empty a row). Not a cover. - `:diagfeasible` — a cover grown from the diagonal by nearest-neighbor - propagation. Feasible by construction. + propagation. - `:hardcover` — the tightened hard cover of [`symcover`](@ref), which is `:geomean` boosted to feasibility and then tightened. Forwards `maxiter` to the - tightening pass. Feasible by construction, so `feasible` has no effect on it. + tightening pass. `feasible` has no effect on it. `feasible` names how the point is brought up to covering `A` — that is, to `a[i]*a[j] >= abs(A[i,j])`, up to the roundoff of the log-domain arithmetic: - `:inflate` (the default) multiplies every scale by the smallest common factor - that achieves coverage, moving the point bodily and leaving its shape intact. + that achieves coverage. - `:boost` raises only the rows that touch a violated entry, so it changes the shape of the point. This is the route [`symcover`](@ref) itself takes. -- `:none` returns the strategy's own point, with no coverage guarantee. This is - what the soft covers want: forcing the geometric mean to cover `A` would - destroy the very property that makes it the soft AbsLog{2} optimum. +- `:none` returns the strategy's point without a coverage guarantee. -The two feasible routes land on the boundary at different points, and so in -different basins of the non-convex `AbsLinear` objectives — which is exactly why a -menu of starts is worth having, and why the choice is exposed rather than fixed. +The two feasible routes reach different points on the boundary and can enter +different basins of the nonconvex `AbsLinear` objectives. Under every setting the result is strictly positive on every row that carries support and exactly zero on every row that carries none. diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 9b56887..50cf6fa 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -27,18 +27,12 @@ Supported ϕ values: `nnz ≪ n²`. - `AbsLog{1}()`: requires JuMP and HiGHS. - `AbsLinear{1}()`, `AbsLinear{2}()`: requires JuMP and Ipopt. These objectives are - non-convex, so the solver returns the minimum of the basin it starts in. Rather than - commit to one start, these methods refine each of `strategies` — the - [`initialize_symcover`](@ref) menu, by default `$(SYMCOVER_MIN_STRATEGIES)` — and return - the best cover found, at a cost of one solve per start. A strategy that `A` admits no - start for is skipped. The result is the best *local* minimum on that menu: the multistart - is a hedge against a poor basin, not a certificate of global optimality. - -The `AbsLog` penalties are convex in the log-scales, so for them the minimum value is -unique and no such hedge is needed. `AbsLog{2}` has a unique minimizer too. `AbsLog{1}` -does not: its optimum is a whole face of the feasible polytope, whose members are -genuinely different covers that happen to score alike. The one returned is the member -of that face minimizing the `AbsLog{2}` objective. + nonconvex. Each strategy in `strategies` is refined, and the best local + minimum is returned. A strategy that cannot produce a start is skipped. + +The `AbsLog` penalties are convex in the log-scales. `AbsLog{2}` has a unique +minimizer. When the `AbsLog{1}` optimum is a face, the method returns the member +that minimizes the `AbsLog{2}` objective. !!! note Even the native solver is more expensive than the [`symcover`](@ref) heuristic. @@ -52,17 +46,12 @@ symcover_min(A::AbstractMatrix; kwargs...) = symcover_min(AbsLog{2}(), A; kwargs a, b = cover_min(ϕ, A) a, b = cover_min(A) -Return the ϕ-minimal asymmetric hard cover of `A`: the vectors `a`, `b` minimizing +Return the ϕ-minimal asymmetric hard cover of `A`: vectors `a`, `b` minimizing `∑_{i,j} ϕ(|A[i,j]|/(a[i]*b[j]))` subject to `a[i]*b[j] >= |A[i,j]|` for every nonzero -entry of `A`. Only the products `a[i]*b[j]` are determined by the problem: the gauge -`a -> γ*a`, `b -> b/γ` leaves every one of them unchanged, and acts independently on -each connected component of the bipartite support graph of `A` (rows and columns as -vertices, stored nonzeros as edges), since no product spans two components. The split -is pinned by imposing, within each component, the balance convention +entry of `A`. The split between `a` and `b` is set within each support component +by the balance convention `∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]` (`nzaᵢ`, `nzbⱼ` = nonzero counts of row `i`, -column `j`, summed over that component's rows/columns) — so the result is a -deterministic function of the support and the products, and block-diagonal assembly -commutes with the split. The no-ϕ form defaults to `AbsLog{2}()`, matching +column `j`). The no-ϕ form defaults to `AbsLog{2}()`, matching [`cover`](@ref). Supported ϕ values: @@ -77,18 +66,12 @@ Supported ϕ values: reweighted normal equations is the wrong solve when `nnz ≪ n²`. - `AbsLog{1}()`: requires JuMP and HiGHS. - `AbsLinear{1}()`, `AbsLinear{2}()`: requires JuMP and Ipopt. These objectives are - non-convex, so the solver returns the minimum of the basin it starts in. Rather than - commit to one start, these methods refine each of `strategies` — the - [`initialize_cover`](@ref) menu, by default `$(COVER_MIN_STRATEGIES)` — and return the - best cover found, at a cost of one solve per start. The result is the best *local* - minimum on that menu: the multistart is a hedge against a poor basin, not a certificate - of global optimality. - -The `AbsLog` penalties are convex in the log-scales, so for them the minimum value is -unique and no such hedge is needed. `AbsLog{2}` has a unique minimizer too. `AbsLog{1}` -does not: its optimum is a whole face of the feasible polytope, whose members are -genuinely different covers that happen to score alike. The one returned is the member -of that face minimizing the `AbsLog{2}` objective. + nonconvex. Each strategy in `strategies` is refined, and the best local + minimum is returned. + +The `AbsLog` penalties are convex in the log-scales. `AbsLog{2}` has a unique +minimizer. When the `AbsLog{1}` optimum is a face, the method returns the member +that minimizes the `AbsLog{2}` objective. !!! note Even the native solver is more expensive than the [`cover`](@ref) heuristic. @@ -114,12 +97,9 @@ cover `A` — `a[i]*a[j] >= abs(A[i,j])` — to within the roundoff of the log-d arithmetic; otherwise an `ArgumentError` is raised. Scales on rows carrying no support are inert: whatever they hold on input, they are zero on output. -How much the start matters depends on ϕ. Under the `AbsLog` penalties the result is -start-independent: they are convex in the log-scales, `AbsLog{2}` has a unique -minimizer, and `AbsLog{1}` — whose optimum is a whole face of equally-scoring covers -— is pinned to the member of that face minimizing the `AbsLog{2}` objective. The -`AbsLinear` penalties are non-convex, and the identified local minima depend on -the start(s). +The `AbsLog` result is independent of the start. For `AbsLog{1}`, ties are broken +by the `AbsLog{2}` objective. Local minima under `AbsLinear` can depend on the +start. See also: [`initialize_symcover`](@ref), [`symcover_min`](@ref), [`cover_min!`](@ref). """ @@ -198,16 +178,10 @@ function cover_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::Abstra return a, b end -# The AbsLinear objectives are non-convex, so a refinement reports the minimum of whichever -# basin its start lies in. The drivers below therefore refine every start on a menu and keep -# the best, which is what makes the result of the non-mutating entry point a property of `A` -# rather than of an initialization the caller never chose. The kernels they call — the -# `*_min!` refiners for the AbsLinear penalties — live in the MatrixCoversIpoptExt extension, but the -# menu and the selection are native, so the two families need only one description. +# Multistart drivers for nonconvex AbsLinear objectives. The `*_min!` kernels live +# in MatrixCoversIpoptExt; the main package owns the starts and selection. # -# The winner is picked by `_multistart_select`, the same scale-covariant rule the soft-cover -# multistarts use: a later start replaces the incumbent only on a genuine relative -# improvement, never on the roundoff by which two starts reaching the same basin differ. +# Use the same roundoff-tolerant selection rule as the soft-cover multistarts. function symcover_min(ϕ::AbsLinear, A::AbstractMatrix; strategies=SYMCOVER_MIN_STRATEGIES) ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("symcover_min requires a square matrix")) diff --git a/src/penalties.jl b/src/penalties.jl index ce799ab..c30d6b4 100644 --- a/src/penalties.jl +++ b/src/penalties.jl @@ -24,13 +24,9 @@ A subtype must be callable on a nonnegative real: whenever `A[i,j]` is zero, and `cover_objective` passes `typemax` for an entry left uncovered by a zero scale. Penalties are conventionally singleton structs. -That call is the whole contract, and it buys exactly one thing: -[`cover_objective`](@ref) works for any subtype. **The solvers do not.** Every -solver in this package dispatches on a concrete built-in penalty — `AbsLog{2}` -is solved natively, the `AbsLinear` penalties through JuMP — so a custom subtype -passed to [`symcover_min`](@ref), [`soft_symcover`](@ref), or any other solver -raises a `MethodError`. Scoring covers with your own penalty is supported; -minimizing it is not. +[`cover_objective`](@ref) works for any subtype, but solvers support only +specific built-in penalties: `AbsLog{2}` natively and `AbsLinear` through JuMP. +Passing a custom subtype to a solver raises a `MethodError`. """ abstract type AbstractCoverPenalty<:Function end diff --git a/src/soft_covers.jl b/src/soft_covers.jl index f7b1ff7..05899e3 100644 --- a/src/soft_covers.jl +++ b/src/soft_covers.jl @@ -16,42 +16,25 @@ minimizing the soft-cover objective `∑_{i,j} ϕ(|A[i,j]| / (a[i]*a[j]))`. Unlike [`symcover`](@ref), there is no hard coverage constraint: `a[i]*a[j]` may be less than `|A[i,j]|`, with violations penalized by `ϕ`. -Supported penalty functions: -- `AbsLog{2}()`: convex, and returns its exact unconstrained minimum from a single linear - solve. Identical to [`soft_symcover_min`](@ref)`(AbsLog{2}(), A)` — with one minimizer - there is nothing for a heuristic and a minimizer to disagree about. -- `AbsLog{1}()`: initializes from the AbsLog{2} minimum, then refines by coordinate descent - with a log-space weighted-median step, reaching a deterministic and scale-covariant fixed - point. That point is not in general a minimizer: each step minimizes exactly over one - coordinate, but the objective's nonsmoothness couples `a[i]` with `a[j]`, so the descent - can settle where no single-coordinate move improves and the objective still sits - materially above its minimum. [`soft_symcover_min`](@ref) does not yet offer an exact - `AbsLog{1}` alternative. -- `AbsLinear{2}()` (default): non-convex; refined by coordinate descent from `starts` - scale-covariant starting points, keeping the lowest-objective result (see below). -- `AbsLinear{1}()`: initializes from the `AbsLinear{2}()` result, coordinate descent uses a - weighted-median step. - -For the `AbsLinear` penalties the objective is non-convex, so `starts` starting points are -tried and the best kept, taken in this order: the geometric-mean minimum, the tightened hard -cover, the geometric-mean minimum inflated uniformly until it covers `A`, a leave-one-out -geometric mean that drops the support entry with the most negative log-residual (this start -keeps the result continuous as an entry `|A[i,j]|` approaches zero), and — only when `A` has a -zero entry — a greedy feasible cover. Any slots left over are multiplicative log-normal -perturbations `a .* exp.(σ .* ξ)` of the geometric-mean point with spread `σ`, `ξ` drawn from -`rng`; at the default `starts=5` there is at most one such perturbation. Every start co-varies -with a diagonal rescaling of `A` and the objective is scale-invariant, so the selection is -scale-covariant. The default `rng` is a fresh `MersenneTwister(0)` per call, making repeated -calls (and the two frames of a covariance check) agree; pass your own `rng` for reproducibility -you control, since default RNG streams are not stable across Julia versions. `sigma` is -accepted as an ASCII alias for `σ`. +Supported penalties are: + +- `AbsLog{2}()`: the convex minimum, computed by one linear solve. +- `AbsLog{1}()`: weighted-median coordinate descent to a fixed point, which need + not be a local minimum. +- `AbsLinear{2}()` (default): multistart coordinate descent. +- `AbsLinear{1}()`: weighted-median descent initialized from the `AbsLinear{2}` + result. + +For `AbsLinear`, `starts` controls the number of deterministic and perturbed +starting points. Perturbations have the form `a .* exp.(σ .* ξ)`, where `ξ` is +drawn from `rng`. The default RNG is reset for each call. Pass an explicit `rng` +to control reproducibility. `sigma` is an ASCII alias for `σ`. See also: [`symcover`](@ref), [`cover_objective`](@ref), [`soft_symcover_min`](@ref). # Examples -The multistart converges to the covariant minimizer to within its objective -tolerance; round to compare against exact values. +Round the multistart result when comparing it with exact values. ```jldoctest julia> A = [4 -1; -1 0]; @@ -111,19 +94,15 @@ Refine the starting point `a` into a symmetric soft cover of `A`, in place, and The no-ϕ form defaults to `AbsLinear{2}()`, matching [`soft_symcover`](@ref), whose supported ϕ values these methods share. -This is the refiner half of [`soft_symcover`](@ref): the non-mutating form owns a -multistart menu, so its result is a property of `A`, while this one descends from the -single start you hand it, so its result is a property of `A` *and* that start. Building -the start is the caller's job — see [`initialize_symcover`](@ref), and pass -`feasible=:none`, since a soft cover is under no obligation to cover `A`. +Unlike [`soft_symcover`](@ref), this method refines one caller-provided start. +Build it with [`initialize_symcover`](@ref) and `feasible=:none`. `a` must be finite and strictly positive on every row of `A` that carries support; scales on rows carrying no support are inert, and are zero on output. Unlike [`symcover_min!`](@ref), `a` need *not* cover `A` — the soft objective imposes no coverage constraint. -`maxiter` bounds the descent sweeps; its default matches the corresponding -[`soft_symcover`](@ref) method. Under `AbsLog{2}` the objective is convex with a unique -minimizer, so the start is honored but not visible in the result. +`maxiter` bounds the descent sweeps. Under `AbsLog{2}`, the unique minimizer is +independent of the start. See also: [`soft_symcover`](@ref), [`soft_symcover_min!`](@ref), [`initialize_symcover`](@ref), [`soft_cover!`](@ref). """ @@ -166,44 +145,24 @@ objective `∑_{i,j} ϕ(|A[i,j]| / (a[i]*b[j]))`. This is the asymmetric analog Unlike [`cover`](@ref), there is no hard coverage constraint: `a[i]*b[j]` may be less than `|A[i,j]|`, with violations penalized by `ϕ`. -Supported penalty functions: -- `AbsLog{2}()`: convex, and returns its exact unconstrained minimum from a single linear - solve. Identical to [`soft_cover_min`](@ref)`(AbsLog{2}(), A)` — with one minimizer - there is nothing for a heuristic and a minimizer to disagree about. -- `AbsLog{1}()`: initializes from the `AbsLog{2}()` minimum, then refines by alternating - weighted-median row and column updates, reaching a deterministic and scale-covariant fixed - point. As in [`soft_symcover`](@ref), that point is not in general a minimizer: each - half-sweep minimizes exactly, but the objective's nonsmoothness couples `a[i]` with `b[j]`, - so the descent can settle where no such sweep improves and the objective still sits - materially above its minimum. [`soft_cover_min`](@ref) does not yet offer an exact - `AbsLog{1}` alternative. -- `AbsLinear{2}()` (default): in the inverse-scale variables `u = 1 ./ a`, `v = 1 ./ b`, the - objective `∑_{i,j∈S} (1 - |A[i,j]| u[i] v[j])²` (sum over the nonzero support `S`) is - biconvex, so alternating least squares with the closed-form half-sweeps - - u[i] = ∑_j |A[i,j]| v[j] / ∑_j (|A[i,j]| v[j])² (dually for v[j]) - - is monotone, stopping when the relative objective decrease falls to rounding level - for the element type, or after `maxiter` sweeps. -- `AbsLinear{1}()`: initializes from the `AbsLinear{2}()` result, then refines by alternating - weighted-median updates — each row/column block is minimized exactly, so the descent is - monotone. Its flat basins are broken by a deterministic lower-median tie-break, giving a - scale-covariant representative. +Supported penalties are: + +- `AbsLog{2}()`: the convex minimum, computed by one linear solve. +- `AbsLog{1}()`: alternating weighted-median updates to a fixed point, which + need not be a local minimum. +- `AbsLinear{2}()` (default): alternating least squares. +- `AbsLinear{1}()`: alternating weighted-median updates initialized from the + `AbsLinear{2}` result. Rows or columns of `A` that are entirely zero receive scale `0`. As with [`cover`](@ref), only the products `a[i] * b[j]` are determined by the problem; the split is fixed by the balance convention `∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]`, imposed within each connected component of the support (the gauge acts independently on each). -The objective is non-convex, so `starts` starting points are tried and the lowest-objective -result kept: the geometric mean boosted until it covers `A`, the tightened hard cover -[`cover`](@ref), and — for the remaining starts — multiplicative log-normal perturbations -`a .* exp.(σ .* ξ)`, with spread `σ`, of that boosted point, `ξ` drawn from `rng`. Every start co-varies with an -independent row/column rescaling of `A` and the objective is scale-invariant, so the selection -is scale-covariant. The default `rng` is a fresh `MersenneTwister(0)` per call, making repeated -calls (and the two frames of a covariance check) agree; pass your own `rng` for reproducibility -you control, since default RNG streams are not stable across Julia versions. `sigma` is accepted -as an ASCII alias for `σ`. +For `AbsLinear`, `starts` controls the number of deterministic and perturbed +starting points. Perturbations have spread `σ` and use `rng`. The default RNG is +reset for each call; pass one explicitly to control reproducibility. `sigma` is +an ASCII alias for `σ`. See also: [`cover`](@ref), [`soft_symcover`](@ref), [`cover_objective`](@ref). @@ -253,18 +212,13 @@ end a, b = soft_cover!(ϕ, a, b, A; maxiter=...) a, b = soft_cover!(a, b, A; maxiter=...) -Refine the starting point `(a, b)` into a soft cover of `A`, in place, and return it. This -is the asymmetric counterpart of [`soft_symcover!`](@ref) and the refiner half of -[`soft_cover`](@ref), carrying the same contract on the start: finite and strictly -positive on every supported row and column, inert (and zero on output) elsewhere, and -under no obligation to cover `A`. Build one with [`initialize_cover`](@ref) and -`feasible=:none`. The no-ϕ form defaults to `AbsLinear{2}()`, matching -[`soft_cover`](@ref). +Refine the starting point `(a, b)` into a soft cover of `A` in place. Scales must +be finite and positive on supported rows and columns; unsupported scales are +zeroed. The start need not cover `A`. Build one with [`initialize_cover`](@ref) +and `feasible=:none`. The no-ϕ form uses `AbsLinear{2}()`. -The product `a[i]*b[j]` is unchanged by `a -> c*a`, `b -> b/c`, so the start is read only -up to that gauge: `(a, b)` and `(2a, b/2)` give the same result. The result itself is -pinned to the balance convention of [`cover_min`](@ref), as every asymmetric cover in this -package is. +The result uses the balance convention of [`cover_min`](@ref), so equivalent +rescalings `(c*a, b/c)` give the same result. See also: [`soft_cover`](@ref), [`soft_cover_min!`](@ref), [`initialize_cover`](@ref), [`soft_symcover!`](@ref). """ @@ -307,21 +261,11 @@ with no coverage constraints. The no-ϕ form defaults to `AbsLinear{2}()`, match [`soft_symcover`](@ref). Supported ϕ values and required extensions: -- `AbsLog{2}()`: solved natively (no external solver). In log space the objective is a - linear least-squares, so one solve settles it, and being convex it has a unique - minimizer that no start can influence. `linsolve` selects the inner solve, exactly as - in [`symcover_min`](@ref). -- `AbsLinear{1}()`, `AbsLinear{2}()`: requires JuMP and Ipopt. These objectives are - non-convex, so the solver returns the minimum of the basin it starts in. Rather than - commit to one start, these methods refine each of `strategies` — the - [`initialize_symcover`](@ref) menu, by default `$(SYMCOVER_MIN_STRATEGIES)`, without - forcing feasibility — and return the best cover found, at a cost of one solve per start. -- `AbsLog{1}()`: not yet implemented. The objective is an LP in log space, but its optimum - is a face, and the lexicographic AbsLog{2} selection that [`symcover_min`](@ref) uses to - pin one member of the corresponding hard face does not carry over: the hard face is bounded - by the coverage constraints, while this one is a level set of an unconstrained piecewise- - linear objective, across which the quadratic pulls far enough to cost most of the exactly - tight residuals that make `AbsLog{1}` worth choosing. +- `AbsLog{2}()`: solved natively as linear least squares; `linsolve` has the same + meaning as in [`symcover_min`](@ref). +- `AbsLinear{1}()`, `AbsLinear{2}()`: require JuMP and Ipopt. Each strategy in + `strategies` is refined, and the best local minimum is returned. +- `AbsLog{1}()`: not implemented. See also: [`soft_symcover_min!`](@ref), [`soft_symcover`](@ref), [`symcover_min`](@ref). """ @@ -357,22 +301,14 @@ end a = soft_symcover_min!(ϕ, a, A) a = soft_symcover_min!(a, A) -Refine the starting point `a` into the ϕ-minimal symmetric soft cover of `A`, in place. -This is the soft counterpart of [`symcover_min!`](@ref), and the second half of the -initialize/refine pair whose first half is [`initialize_symcover`](@ref). The no-ϕ form -defaults to `AbsLinear{2}()`, matching [`soft_symcover_min`](@ref), whose supported ϕ -values these methods share. +Refine `a` into a local minimum of the symmetric soft-cover objective, in place. +The no-ϕ form uses `AbsLinear{2}()`. -`a` must be strictly positive on every row of `A` that carries support; scales on rows -carrying no support are inert, and are zero on output. Unlike [`symcover_min!`](@ref), `a` -need *not* cover `A` — the soft objective imposes no coverage constraint, and the natural -starts do not satisfy one. Pass `feasible=:none` when building a start with -[`initialize_symcover`](@ref). +`a` must be positive on supported rows; unsupported scales are zeroed. It need +not cover `A`. Use `feasible=:none` with [`initialize_symcover`](@ref). -The `AbsLinear` penalties are non-convex, so the start selects the local minimum the solver -descends into; that is why [`soft_symcover_min`](@ref) tries several rather than committing -to one. Under `AbsLog{2}` the objective is convex with a unique minimizer, so the start is -honored but not visible in the result. +For `AbsLinear`, the result can depend on the start. The `AbsLog{2}` result is +unique. See also: [`initialize_symcover`](@ref), [`soft_symcover_min`](@ref), [`symcover_min!`](@ref). """ @@ -386,10 +322,7 @@ function soft_symcover_min!(::AbsLog{2}, a::AbstractVector, A::AbstractMatrix; k return a end -# Shared prologue of the symmetric soft refiners (`soft_symcover!`, `soft_symcover_min!`). -# The soft objective constrains nothing, so — unlike `_prepare_symcover_start!` — this -# checks positivity only, and moves the start nowhere. `fname` names the caller so the -# error reports the function the user actually called. +# Validate a symmetric soft-cover start and clear unsupported scales. function _prepare_soft_symcover_start!(a::AbstractVector, A::AbstractMatrix, fname::Symbol=:soft_symcover_min!) ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("$fname requires a square matrix")) @@ -427,25 +360,12 @@ the objective depends on `a` and `b` only through the products `a[i]*b[j]`, so w convention the split between them would be arbitrary. Supported ϕ values and required extensions: -- `AbsLog{2}()`: solved natively (no external solver) — the same analytic geometric-mean - minimum [`cover`](@ref) computes as its initial point. Convex, so the minimizer is unique. -- `AbsLinear{1}()`, `AbsLinear{2}()`: requires JuMP and Ipopt. These objectives are - non-convex, so the solver returns the minimum of the basin it starts in. Rather than - commit to one start, these methods refine each of `strategies` — the - [`initialize_cover`](@ref) menu, by default `$(COVER_MIN_STRATEGIES)`, taken raw - (`feasible=:none`, since the soft objective constrains nothing) — and return the best - cover found, at a cost of one solve per start. The result is the best *local* minimum on - that menu: the multistart is a hedge against a poor basin, not a certificate of global - optimality. -- `AbsLog{1}()`: not yet implemented. The objective is an LP in log space, but its optimum - is a face, and the lexicographic AbsLog{2} selection that [`symcover_min`](@ref) uses to - pin one member of the corresponding hard face does not carry over: the hard face is bounded - by the coverage constraints, while this one is a level set of an unconstrained piecewise- - linear objective, across which the quadratic pulls far enough to cost most of the exactly - tight residuals that make `AbsLog{1}` worth choosing. - -Every start on the menu co-varies with a rescaling of `A` and the objective is -scale-invariant, so the selection — and hence the result — is scale-covariant. +- `AbsLog{2}()`: solved natively; the minimizer is unique. +- `AbsLinear{1}()`, `AbsLinear{2}()`: require JuMP and Ipopt. Each strategy in + `strategies` is refined, and the best local minimum is returned. +- `AbsLog{1}()`: not implemented. + +The starting points and objective are scale-covariant, as is the selected result. See also: [`soft_cover_min!`](@ref), [`soft_symcover_min`](@ref), [`soft_cover`](@ref). """ @@ -473,19 +393,14 @@ end a, b = soft_cover_min!(ϕ, a, b, A) a, b = soft_cover_min!(a, b, A) -Refine the starting point `(a, b)` into the ϕ-minimal asymmetric soft cover of `A`, in -place. This is the asymmetric counterpart of [`soft_symcover_min!`](@ref). The no-ϕ form -defaults to `AbsLinear{2}()`, matching [`soft_cover_min`](@ref), whose supported ϕ values -these methods share. +Refine `(a, b)` into a local minimum of the asymmetric soft-cover objective, in +place. The no-ϕ form uses `AbsLinear{2}()`. -`a` and `b` must be strictly positive on every supported row and column; scales on -unsupported rows and columns are inert, and are zero on output. As with -[`soft_symcover_min!`](@ref) — and unlike [`cover_min!`](@ref) — the start need *not* cover -`A`. Build one with `feasible=:none`. +`a` and `b` must be positive on supported rows and columns; unsupported scales +are zeroed. The start need not cover `A`. Build one with `feasible=:none`. -The product `a[i]*b[j]` is unchanged by `a -> c*a`, `b -> b/c`, so the start is read only -up to that gauge, and the result is pinned to the balance convention of -[`soft_cover_min`](@ref). +The result uses the balance convention of [`soft_cover_min`](@ref), so equivalent +rescalings `(c*a, b/c)` give the same result. See also: [`initialize_cover`](@ref), [`soft_cover_min`](@ref), [`soft_symcover_min!`](@ref). """ @@ -501,11 +416,7 @@ function soft_cover_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::A return a, b end -# Shared prologue of the asymmetric soft refiners (`soft_cover!`, `soft_cover_min!`); the -# counterpart of `_prepare_soft_symcover_start!`. Positivity only — the soft objective -# constrains nothing — plus the balance pin, so the refiners read the start only up to the -# row/column gauge. `fname` names the caller so the error reports the function the user -# actually called. +# Validate an asymmetric soft-cover start, clear unsupported scales, and balance it. function _prepare_soft_cover_start!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, fname::Symbol=:soft_cover_min!) 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)))")) @@ -554,27 +465,12 @@ function _resolve_alias(primary, alias, default, primary_name::Symbol, alias_nam return primary end -# Relative-objective margin a later start must beat the incumbent by to replace it. -# -# Two candidates that reach the same point differ, between one frame and a rescaled one, only -# by roundoff in evaluating the objective — a sum of O(n²) terms, so of relative size O(n·eps). -# The margin must exceed that, or the selection could flip with the frame and forfeit -# covariance; it must also stay below any genuine basin gap, which is orders of magnitude -# larger. Candidates stopped short of convergence may differ by far more than roundoff, but -# such differences are deterministic and co-vary with the frame, so switching on them is safe. +# Relative improvement required to replace the incumbent. This prevents +# roundoff-equivalent candidates from changing the selection after rescaling. _multistart_switchtol(::Type{T}) where {T} = 5_000_000 * eps(T) -# Labeled candidate starts for the symmetric AbsLinear{2} multistart, in selection order. -# The deterministic starts are the `initialize_symcover` menu: the geometric mean (raw — it is -# the exact soft AbsLog{2} optimum, so forcing it to feasibility would spoil it — and also the -# perturbation base), the tightened hard cover, the uniformly inflated geometric mean, the -# leave-one-out geometric mean (when a support entry can be dropped), and — only when `A` has -# a zero entry — the greedy feasible cover. The feasible start is gated because on a fully -# dense `A` it never uniquely wins, and "which entries are zero" is invariant under a diagonal -# rescaling `D*A*D`, so the gate keeps the selection scale-covariant. Remaining slots, up to -# `starts` total, are multiplicative log-normal perturbations `a_g .* exp.(σ .* ξ)` of the -# geometric-mean point, `ξ` drawn from `rng` (drawn for every index so the stream is -# frame-independent). `starts` below the number of deterministic starts truncates the list. +# Symmetric AbsLinear{2} starts in selection order, followed by log-normal +# perturbations of the geometric-mean point. `starts` truncates or extends the list. function _soft_symcover_abslinear2_inits(A::AbstractMatrix, starts::Int, σ::Real, rng) ax = axes(A, 1) T = float(real(eltype(A))) @@ -608,13 +504,7 @@ function _soft_symcover_abslinear2_inits(A::AbstractMatrix, starts::Int, σ::Rea return labels, inits end -# Index of the multistart winner among candidate objectives `objs`: the earliest candidate not -# beaten by a strict relative improvement. Switching only on a genuine improvement is what keeps -# the selection scale-covariant — candidates landing in the same basin converge to the same -# objective only to the descent tolerance (~1e-14), and switching on that noise would forfeit -# covariance, since the incumbent (ordered geometric-mean first) is the most covariant start. -# Real basin improvements are far larger. This is the single source of the selection rule, so a -# caller that captures `objs` (below) recovers the winner exactly, without re-deriving it. +# Return the first candidate not beaten by more than the roundoff margin. function _multistart_select(objs) besti = firstindex(objs) Ebest = objs[besti] @@ -662,14 +552,16 @@ end # Coordinate-descent iteration for AbsLinear{2} soft cover. # Each coordinate a[k] is updated to the exact minimizer of -# (1 - d/x²)² + ∑_{j≠k} (1 - c_j/x)² -# where d = |A[k,k]| and c_j = |A[k,j]|/a[j]. +# ½(1 - d/x²)² + ∑_{j≠k} (1 - c_j/x)² +# where d = |A[k,k]| and c_j = |A[k,j]|/a[j]. This is half the part of the +# `cover_objective` sum that depends on a[k]: that sum runs over the full grid, +# so each off-diagonal pair contributes twice and the diagonal once. # Closed form when d=0 (x = s2/s1); Newton on a cubic otherwise. # # `iter` bounds the sweeps; the descent exits early once every coordinate's -# stationarity residual r_k = ∑_j (1 - ρ)ρ (ρ = |A[k,j]|/(a[k]a[j])), the -# gradient of the objective in log a[k], has magnitude below `tol` at the start -# of a sweep. The residual is available for free from the sums already formed +# stationarity residual r_k = ∑_j (1 - ρ)ρ (ρ = |A[k,j]|/(a[k]a[j])), half the +# gradient of the objective above in log a[k], has magnitude below `tol` at the +# start of a sweep. The residual is available for free from the sums already formed # for the update (r_k = s1/a[k] - s2/a[k]² + d/a[k]² - d²/a[k]⁴), it is the exact # quantity optimality demands be zero, and it is scale-invariant (each ρ is), so # covariant restarts of a rescaled problem exit on the same sweep and the @@ -753,16 +645,20 @@ function _weighted_self_median!(c::AbstractVector{T}) where T return wm end -# AbsLinear{1} coordinate objective at candidate `x`: |1 - d/x²| + ∑ᵢ |1 - cᵢ/x|. +# AbsLinear{1} coordinate objective at candidate `x`: |1 - d/x²| + 2∑ᵢ |1 - cᵢ/x|, +# the part of the full-grid `cover_objective` sum that depends on a[k] (each +# off-diagonal pair appears twice there, the diagonal once). # This is a top-level function because a closure in `_abslinear1_iter!` would # capture and box the reassigned `d`, allocating in the inner loop and causing # juliac's trim verifier to report a dynamic call. -_abslinear1_obj(x, d, c) = abs(1 - d/x^2) + sum(abs(1 - ci/x) for ci in c) +_abslinear1_obj(x, d, c) = abs(1 - d/x^2) + 2 * sum(abs(1 - ci/x) for ci in c) # Coordinate-descent iteration for AbsLinear{1} soft cover. -# Each coordinate a[k] is updated to minimize ∑_j |1 - |A[k,j]|/(a[k]*a[j])|. -# For the off-diagonal sum, the minimizer is the weighted median of c_j with weights c_j, -# where c_j = |A[k,j]|/a[j]. When A[k,k] ≠ 0 we also compare against sqrt(|A[k,k]|). +# Each coordinate a[k] is updated to reduce |1 - d/x²| + 2∑_{j≠k} |1 - c_j/x|, +# where d = |A[k,k]| and c_j = |A[k,j]|/a[j]. The off-diagonal sum is minimized +# by the weighted median of the c_j with weights c_j (the factor 2 does not move +# it). When d ≠ 0 the update takes the better of that median and sqrt(d) under +# `_abslinear1_obj`; this is not an exact coordinate minimization. # # `iter` bounds the sweeps; the descent exits early once the largest relative # coordinate movement in a sweep drops to `tol`. The median update reaches an @@ -941,9 +837,10 @@ function _soft_cover_abslinear2_inits(A::AbstractMatrix, starts::Int, σ::Real, ag, bg = initialize_cover(A; strategy=:geomean, feasible=:boost) labels = ["boost"] inits = [(copy(ag), copy(bg))] - # The tightened hard cover `cover(A)` is exactly this point tightened, so tighten a copy + # `cover(A)` is this point tightened, then balanced and re-inflated. Tighten a copy # (at `tighten_cover!`'s own default `maxiter`) rather than recomputing the shared - # geometric-mean and boost passes. + # geometric-mean and boost passes; the balance is gauge-only and the re-inflation + # only recovers roundoff, so this start differs from `cover(A)` negligibly. length(inits) < starts && (push!(labels, "hardcover"); push!(inits, tighten_cover!(copy(ag), copy(bg), A))) k = 0 while length(inits) < starts diff --git a/test/heuristic_covers.jl b/test/heuristic_covers.jl index e3bec33..6f08e44 100644 --- a/test/heuristic_covers.jl +++ b/test/heuristic_covers.jl @@ -33,10 +33,7 @@ end @testset "symcover does not currently consult ϕ" begin - # A regression check on the heuristic as it stands, not an API guarantee: the docstring - # says the heuristic covers ignore `ϕ` *currently*, and that this may change. Should a - # penalty-tuned heuristic land, this test records what changes — update it; do not read a - # failure here as a broken promise to callers. + # The current heuristic accepts but ignores the penalty. rng = StableRNG(1) for n in (2, 5, 40) B = randn(rng, n, n); A = (B + B') / 2 @@ -46,8 +43,7 @@ end end end - # Ignored is not the same as unchecked: the slot takes a penalty, so a wrong - # first argument fails rather than being silently dropped. + # The first argument must still be a penalty. A = [4.0 1.5; 1.5 1.0] a = symcover(A) b = copy(a) diff --git a/test/initializers.jl b/test/initializers.jl index 26fb4c1..e88e625 100644 --- a/test/initializers.jl +++ b/test/initializers.jl @@ -80,11 +80,7 @@ end @testset "the strategies are distinct starts" begin - # :hardcover raises only the rows touching violated entries; an inflated :geomean moves - # the whole point bodily to the boundary, so it lifts the already-slack rows too. - # Reaching the boundary by different routes is what gives the AbsLinear multistart - # different basins to choose between. Here rows 1-2 are the violated block and row 3 is - # slack, so only the inflation disturbs row 3. + # :hardcover changes only violated rows; inflation also changes slack rows. A = [1.0 0.1 0.0; 0.1 1.0 0.0; 0.0 0.0 100.0] ah = initialize_symcover(A; strategy=:hardcover) ag = initialize_symcover(A; strategy=:geomean) diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 0c344cc..4e53aa2 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -212,9 +212,7 @@ end end @testset "soft_cover_min native AbsLog{2}" begin - # `A` has no zero entry, the case in which the geometric mean coincides with the - # minimum; the two compute it differently, so they agree to roundoff, not bitwise. - # On a sparse support they part company -- see the oracle in `test/soft_covers.jl`. + # On dense support, the geometric mean equals the minimum up to roundoff. A = [1.0 2.0 3.0; 6.0 5.0 4.0] a, b = soft_cover_min(AbsLog{2}(), A) a_ref, b_ref = similar(a), similar(b)