diff --git a/docs/src/index.md b/docs/src/index.md index 6eb8249..f8b5c44 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -250,12 +250,15 @@ reach for JuMP with Ipopt only for the `AbsLinear` penalties. They do not accep ### Uniqueness The `AbsLog{2}()` penalty generally has a unique minimum, with one exception: -row/column scaling `a → γ*a`, `b → b/γ` does not affect `C` and thus invisible -to the objective function. For non-symmetric (i.e., not `symcover`) problems, -the scaling of each is pinned separately by the balance convention -`∑ n_i log a[i] = ∑ m_j log b[j]`, where `n_i`, `m_j` are the the nonzero counts -of row `i` and column `j`, respectively. This convention is not scale-invariant -but has no impact on the cover itself. +row/column scaling `a → γ*a`, `b → b/γ` does not affect `C` and is thus +invisible to the objective function. For non-symmetric (i.e., not `symcover`) problems, +the scaling of each is pinned by the balance convention +`∑ n_i log a[i] = ∑ m_j log b[j]`, where `n_i`, `m_j` are the nonzero counts +of row `i` and column `j`, respectively. The gauge freedom, and hence this +convention, acts independently on each connected component of the bipartite +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 diff --git a/ext/MatrixCoversIpoptExt.jl b/ext/MatrixCoversIpoptExt.jl index 16ebb6f..922fcf9 100644 --- a/ext/MatrixCoversIpoptExt.jl +++ b/ext/MatrixCoversIpoptExt.jl @@ -120,8 +120,11 @@ end # a[i]*b[j] is invariant under (α, β) → (α + s, β - s), so — unlike the symmetric # problem, which has no such freedom — the model is degenerate along that direction # until the balance constraint ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ pins it, exactly as -# cover_min(::AbsLog{1}) does. Without it the split Ipopt reports between `a` and `b` -# would be arbitrary. +# cover_min(::AbsLog{1}) does. That constraint pins only the global gauge direction; +# a support with more than one connected component carries one such gauge per +# component (see MatrixCovers._support_components), so each kernel below finishes +# with a post-solve balance shift (`_balance_cover!`, then `inflate_feasible!` to +# restore exact coverage) that pins the rest. # ============================================================ function MatrixCovers.cover_min!(::AbsLinear{2}, a::AbstractVector, b::AbstractVector, A) @@ -153,7 +156,8 @@ function MatrixCovers.cover_min!(::AbsLinear{2}, a::AbstractVector, b::AbstractV for (j, k) in pairs(pc) b[k] = nzb[j] > 0 ? exp(JuMP.value(β[j])) : zero(T) end - return a, b + MatrixCovers._balance_cover!(a, b, A) + return MatrixCovers.inflate_feasible!(a, b, A) end function MatrixCovers.cover_min!(::AbsLinear{1}, a::AbstractVector, b::AbstractVector, A) @@ -189,7 +193,8 @@ function MatrixCovers.cover_min!(::AbsLinear{1}, a::AbstractVector, b::AbstractV for (j, k) in pairs(pc) b[k] = nzb[j] > 0 ? exp(JuMP.value(β[j])) : zero(T) end - return a, b + MatrixCovers._balance_cover!(a, b, A) + return MatrixCovers.inflate_feasible!(a, b, A) end # ============================================================ @@ -254,8 +259,12 @@ end # ============================================================ # Soft cover: soft_cover_min!(::AbsLinear{p}, a, b, A) # The bipartite analog of soft_symcover_min!, and the unconstrained analog of cover_min!: -# no coverage constraints, but the same row/column gauge, pinned by the same balance -# constraint ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ. A zero entry of `A` contributes ϕ(0) = 1 whatever the +# no coverage constraints, but the same row/column gauge, pinned in the model by the same +# balance constraint ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ. As in cover_min!, that constraint pins only the +# global gauge direction, so each kernel below finishes with a post-solve `_balance_cover!` +# that pins the rest (one per connected component of the support); unlike the hard-cover +# kernels, no `inflate_feasible!` follows, since the soft objective imposes no coverage +# constraint for it to restore. A zero entry of `A` contributes ϕ(0) = 1 whatever the # scales, so the count of zeros enters the objective as a constant, matching cover_objective. # ============================================================ @@ -286,7 +295,7 @@ function MatrixCovers.soft_cover_min!(::AbsLinear{2}, a::AbstractVector, b::Abst for (j, k) in pairs(pc) b[k] = nzb[j] > 0 ? exp(JuMP.value(β[j])) : zero(T) end - return a, b + return MatrixCovers._balance_cover!(a, b, A) end function MatrixCovers.soft_cover_min!(::AbsLinear{1}, a::AbstractVector, b::AbstractVector, A) @@ -320,7 +329,7 @@ function MatrixCovers.soft_cover_min!(::AbsLinear{1}, a::AbstractVector, b::Abst for (j, k) in pairs(pc) b[k] = nzb[j] > 0 ? exp(JuMP.value(β[j])) : zero(T) end - return a, b + return MatrixCovers._balance_cover!(a, b, A) end end # module MatrixCoversIpoptExt diff --git a/ext/MatrixCoversJuMPExt.jl b/ext/MatrixCoversJuMPExt.jl index 5f1c9b5..054f87d 100644 --- a/ext/MatrixCoversJuMPExt.jl +++ b/ext/MatrixCoversJuMPExt.jl @@ -145,6 +145,10 @@ function MatrixCovers.cover_min_jump(::AbsLog{2}, A) @constraint(model, α[ei[e]] + β[ej[e]] - elog[e] >= 0) end nza, nzb = _degrees(ei, m), _degrees(ej, n) + # Pins only the global gauge direction; a disconnected support carries one (e; -e) + # gauge per component (see MatrixCovers._support_components), so the remaining + # directions are left to whichever vertex HiGHS returns. The post-solve balance + # shift below fixes all of them, matching the native solver. @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n)) JuMP.optimize!(model) check_solved(model, "cover_min_jump") @@ -156,7 +160,8 @@ function MatrixCovers.cover_min_jump(::AbsLog{2}, A) for (j, k) in pairs(pc) b[k] = nzb[j] > 0 ? exp(JuMP.value(β[j])) : zero(T) end - return a, b + MatrixCovers._balance_cover!(a, b, A) + return MatrixCovers.inflate_feasible!(a, b, A) end MatrixCovers.cover_min(::AbsLog{1}, A) = _cover_min_abslog1(A, nothing) @@ -170,8 +175,10 @@ function MatrixCovers.cover_min!(::AbsLog{1}, a::AbstractVector, b::AbstractVect end # Asymmetric counterpart of `_symcover_min_abslog1`, on the bipartite support: the -# same LP over row scales α and column scales β, with the row/column gauge pinned by -# the balance constraint so the split between `a` and `b` is deterministic. +# same LP over row scales α and column scales β. The balance constraint below pins +# the global row/column gauge; a support with more than one connected component +# carries additional per-component gauges that the post-solve balance shift pins, +# so the split between `a` and `b` is deterministic. function _cover_min_abslog1(A, start) axr = axes(A, 1) axc = axes(A, 2) @@ -203,7 +210,9 @@ function _cover_min_abslog1(A, start) nza, nzb = rowcount, colcount # Gauge pin: the products a[i]*b[j] are unchanged by a -> c*a, b -> b/c, so without this # the split between `a` and `b` would be arbitrary. It is orthogonal to the AbsLog{1} - # degeneracy the second stage resolves, and stays in force there. + # degeneracy the second stage resolves, and stays in force there. It pins only the + # global gauge direction, one of possibly several (one per connected component of the + # support); the post-solve balance shift below pins the rest. @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n)) JuMP.optimize!(model) check_solved(model, "cover_min") @@ -217,7 +226,8 @@ function _cover_min_abslog1(A, start) for (j, k) in pairs(pc) b[k] = nzb[j] > 0 ? exp(JuMP.value(β[j])) : zero(T) end - return a, b + MatrixCovers._balance_cover!(a, b, A) + return MatrixCovers.inflate_feasible!(a, b, A) end diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index 9967fa7..44fb463 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -86,7 +86,8 @@ entries are covered first), and `maxiter` tightening iterations are applied. 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 `∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]` (`nzaᵢ`, `nzbⱼ` = nonzero counts of row `i`, -column `j`), as it is throughout the package; see [`cover_min`](@ref). +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. @@ -179,30 +180,49 @@ end # the objective cannot see it, so without a convention the split between `a` and `b` would # be an artifact of whichever pass last touched them. # +# The gauge acts independently on each connected component of the bipartite +# support graph (`_support_components`), so the convention is imposed per +# component: within every component, the row-side and column-side weighted log +# sums agree. This makes the split a well-defined function of the support and +# the products — block-diagonal assembly commutes with balancing — rather than +# pinning only the global scalar and leaving the per-component splits to +# whichever solver internals ran last. Rows and columns with empty support +# belong to no component and are left untouched. +# # Two points differing only by the gauge land on the same point here, so a refiner given # either start cannot tell them apart. The uniform inflation to feasibility raises every -# supported scale of `a` and `b` alike, and `∑ nzaᵢ = ∑ nzbⱼ = nnz`, so it preserves the -# balance it finds. +# supported scale of `a` and `b` alike, and within each component +# `∑ nzaᵢ = ∑ nzbⱼ = nnz`, so it preserves the balance it finds. function _balance_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) T = float(promote_type(eltype(a), eltype(b))) - Lα = Ref(zero(T)) - Lβ = Ref(zero(T)) - nnz = Ref(0) + rowcomp, colcomp, ncomp = _support_components(A) + iszero(ncomp) && return a, b + Lα = zeros(T, ncomp) + Lβ = zeros(T, ncomp) + nnz = zeros(Int, ncomp) + or = first(axes(A, 1)) - 1 + oc = first(axes(A, 2)) - 1 foreach_support(A) do i, j, v - Lα[] += log(T(a[i])) - Lβ[] += log(T(b[j])) - nnz[] += 1 + c = rowcomp[i-or] + Lα[c] += log(T(a[i])) + Lβ[c] += log(T(b[j])) + nnz[c] += 1 + end + s = Lβ + for c in 1:ncomp + s[c] = (Lβ[c] - Lα[c]) / (2 * nnz[c]) end - iszero(nnz[]) && return a, b - s = (Lβ[] - Lα[]) / (2 * nnz[]) - iszero(s) && return a, b # Grow the log-scales directly: `exp(s)` alone can overflow where the shifted # scale is perfectly representable. for i in eachindex(a) - iszero(a[i]) || (a[i] = exp(log(T(a[i])) + s)) + c = rowcomp[i-or] + c == 0 && continue + iszero(s[c]) || iszero(a[i]) || (a[i] = exp(log(T(a[i])) + s[c])) end for j in eachindex(b) - iszero(b[j]) || (b[j] = exp(log(T(b[j])) - s)) + c = colcomp[j-oc] + c == 0 && continue + iszero(s[c]) || iszero(b[j]) || (b[j] = exp(log(T(b[j])) - s[c])) end return a, b end diff --git a/src/initializers.jl b/src/initializers.jl index e56e855..14655b8 100644 --- a/src/initializers.jl +++ b/src/initializers.jl @@ -107,8 +107,8 @@ formulation and raise an `ArgumentError`, as does any unrecognized `strategy` or Under every `feasible` setting the result is strictly positive on every supported row and column and exactly zero on the unsupported ones, and the split between `a` and `b` is fixed by the balance convention -`∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]` that every asymmetric cover in the package -uses (see [`cover_min`](@ref)). +`∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]`, imposed within each connected component of the +support, that every asymmetric cover in the package uses (see [`cover_min`](@ref)). See also: [`initialize_cover!`](@ref), [`initialize_symcover`](@ref), [`cover`](@ref), [`cover_min`](@ref). """ diff --git a/src/iscover.jl b/src/iscover.jl index 2a68554..252015d 100644 --- a/src/iscover.jl +++ b/src/iscover.jl @@ -5,33 +5,23 @@ iscover(a, b, A; rtol=0, atol=0) iscover(a, A; rtol=0, atol=0) -Test whether `a` and `b` cover `A`, that is, whether `a[i]*b[j] >= abs(A[i,j])` for every -entry. The two-argument form tests the symmetric cover `a*a'`, and requires `A` to be -square. - -This is the inequality the whole package is organized around. [`cover`](@ref), -[`symcover`](@ref), and the `*_min` solvers all satisfy it by construction — up to the -roundoff noted below — so the predicate earns its keep mainly on covers that carry no such -guarantee: those from [`soft_cover`](@ref) and [`soft_symcover`](@ref), which penalize -under-coverage rather than forbidding it, and covers a caller has adjusted by hand. - -`rtol` and `atol` supply the slack the producing algorithm warrants, testing -`a[i]*b[j] >= abs(A[i,j])*(1 - rtol) - atol`. Neither defaults to any slack at all, so the -bare call is the exact inequality. Use `rtol` for roundoff that scales with entry magnitude -(the log-domain arithmetic behind the heuristics warrants a few multiples of `eps`), and -`atol` for the convergence tolerance of an iterative solver. `atol` is subtracted only when -it is nonzero: `abs(A[i,j]) - atol` is undefined when the two carry different units, so an -`rtol`-only check never forms it. That makes `rtol` the only one of the two meaningful for a -dimensional `A`, whose entries need not share units — no single scalar `atol` is -commensurate with all of them. - -`a` and `b` must be nonnegative; a negative scale raises an `ArgumentError`. Zero is -allowed, and is what every solver here returns for a row or column with no support. The -requirement is not cosmetic: only nonnegativity makes it sound to skip the zero entries of -`A`, which is what lets this run in time proportional to the support rather than to -`length(A)`. - -`eachindex(a)` must match `axes(A, 1)` and `eachindex(b)` must match `axes(A, 2)`. +Test whether `a` and `b` cover `A`, that is, whether `a[i]*b[j] >= abs(A[i,j])` +for every entry. The two-argument form tests the symmetric cover `a*a'`, and +requires `A` to be square. + +`rtol` and `atol` allow for small violations, testing + + a[i]*b[j] >= abs(A[i,j])*(1 - rtol) - atol + +The default for both tolerances is zero (no slack, test that the cover condition +holds); note that `atol != 0` breaks scale-invariance. + +`a` and `b` must be nonnegative; a negative scale raises an `ArgumentError`. +Zero is allowed, and is what every solver here returns for a row or column with +no support. + +`eachindex(a)` must match `axes(A, 1)` and `eachindex(b)` must match `axes(A, +2)`. See also: [`cover_objective`](@ref), [`cover`](@ref), [`symcover`](@ref). diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index fcb6401..3ec638b 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -54,10 +54,16 @@ symcover_min(A::AbstractMatrix; kwargs...) = symcover_min(AbsLog{2}(), A; kwargs Return the ϕ-minimal asymmetric hard cover of `A`: the 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`. The row/column scales are pinned to the balance convention +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 `∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]` (`nzaᵢ`, `nzbⱼ` = nonzero counts of row `i`, -column `j`) so the result is deterministic. The no-ϕ form defaults to `AbsLog{2}()`, -matching [`cover`](@ref). +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 +[`cover`](@ref). Supported ϕ values: - `AbsLog{2}()`: solved natively (no external solver). Accepts keyword arguments @@ -162,11 +168,13 @@ end # Asymmetric AbsLog{2} hard cover via the same one-sided quadratic penalty as # `symcover_min`, on stacked log-scales x = (α; β) (α = log a over rows, β = log b # over columns) with residuals z_ij = α_i + β_j - log|A_ij|. The row and column -# scales share a gauge freedom (α_i, β_j) → (α_i + s, β_j - s) that leaves every -# residual unchanged; during the solve it is fixed by adding v0*v0ᵀ, -# v0 = [ones(m); -ones(n)], to the normal equations, and afterwards the result is -# shifted along that gauge to the balance convention ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ -# (nzaᵢ, nzbⱼ = nonzero counts of row i, column j) so it is deterministic. +# scales share a gauge freedom (α_i, β_j) → (α_i + s, β_j - s), one dimension per +# connected component of the support, that leaves every residual unchanged; during +# the solve the global one is fixed by adding v0*v0ᵀ, v0 = [ones(m); -ones(n)], to +# the normal equations, and afterwards the result is shifted, per component, to the +# balance convention ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ (nzaᵢ, nzbⱼ = nonzero counts of row i, +# column j, summed within the component) so it is deterministic — see +# `_cover_min_abslog2` for how the remaining per-component gauges are lifted and shifted. function cover_min(::AbsLog{2}, A::AbstractMatrix; kwargs...) a, b, _ = _cover_min_abslog2(A; kwargs...) return a, b @@ -272,7 +280,8 @@ end # Shared prologue of the `cover_min!` kernels; the asymmetric counterpart of # `_prepare_symcover_start!`. The start is additionally pinned to the balance -# convention, so the refiners read it only up to the row/column gauge. +# convention (imposed within each connected component of the support), so the +# refiners read it only up to the per-component row/column gauge. function _prepare_cover_start!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) axes(A, 1) == eachindex(a) || throw(DimensionMismatch("indices of `a` must match row-indexing of `A`, got eachindex(a)=$(eachindex(a)), axes(A, 1)=$(axes(A, 1))")) axes(A, 2) == eachindex(b) || throw(DimensionMismatch("indices of `b` must match column-indexing of `A`, got eachindex(b)=$(eachindex(b)), axes(A, 2)=$(axes(A, 2))")) @@ -582,12 +591,14 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end # Each Newton step solves the reweighted least-squares problem for the stacked # scales x = (α; β), residuals z_ij = α_i + β_j - log|A_ij|. Row and column scales - # share the (e; −e) gauge; both paths pin it. The dense path adds the rank-1 term - # v0*v0ᵀ to the normal equations `B x = f` and factorizes (support-free variables - # get an identity row). The LSQR path appends one gauge row `v0ᵀ x = 0` to the - # least-squares system so `√W R` has full column rank, applies it matrix-free, and - # warm-starts from the incoming iterate. After the solve a closed-form shift moves - # the result to the balance convention, so the pinned gauge is not observable. + # share the global (e; −e) gauge; both paths pin it. The dense path adds the rank-1 + # term v0*v0ᵀ to the normal equations `B x = f` and factorizes (support-free + # variables get an identity row; a support with more than one connected component + # carries additional per-component gauges, lifted by the ridge below). The LSQR + # path appends one gauge row `v0ᵀ x = 0` to the least-squares system so `√W R` has + # full column rank, applies it matrix-free, and warm-starts from the incoming + # iterate. After the solve a closed-form shift, applied within each component, + # moves the result to the balance convention, so the pinned gauge is not observable. f = zeros(T, N) ws = zeros(T, ne) # √weight per support entry (LSQR path) cv = zeros(T, ne + 1) # √weight · log|A_ij|, with a trailing 0 gauge target @@ -702,24 +713,37 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), x[p] += γ end end - # Shift along the (e; -e) gauge to the balance convention ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ. - Lα = zero(T) - Lβ = zero(T) + # Shift along the (e; -e) gauges to the balance convention ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ, + # imposed within each connected component of the support: the gauge acts + # independently on each component, so a single global shift would leave the + # per-component splits wherever the ridge (or LSQR's gauge row) put them. + rowcomp, colcomp, ncomp = _support_components(A) + Lα = zeros(T, ncomp) + Lβ = zeros(T, ncomp) + nec = zeros(Int, ncomp) for ip in 1:m - Lα += nzrow[ip] * x[ip] + c = rowcomp[ip] + c == 0 && continue + Lα[c] += nzrow[ip] * x[ip] + nec[c] += nzrow[ip] end for jp in 1:n - Lβ += nzcol[jp] * x[m+jp] + c = colcomp[jp] + c == 0 && continue + Lβ[c] += nzcol[jp] * x[m+jp] + end + s = Lβ + for c in 1:ncomp + s[c] = (Lβ[c] - Lα[c]) / (2 * nec[c]) end - s = ne > 0 ? (Lβ - Lα) / (2 * ne) : zero(T) # Dense scale vectors matching cover/symcover; `similar(A, …)` is a SparseVector for sparse A. a = similar(Array{T}, axr) b = similar(Array{T}, axc) for (ip, i) in enumerate(axr) - a[i] = hasrow[ip] ? exp(x[ip] + s) : zero(T) + a[i] = hasrow[ip] ? exp(x[ip] + s[rowcomp[ip]]) : zero(T) end for (jp, j) in enumerate(axc) - b[j] = hascol[jp] ? exp(x[m+jp] - s) : zero(T) + b[j] = hascol[jp] ? exp(x[m+jp] - s[colcomp[jp]]) : zero(T) end return a, b, (; nsolves=nsolves[], lsqriters=nlsqr[], linsolve=(use_lsqr ? :lsqr : :dense)) end diff --git a/src/soft_covers.jl b/src/soft_covers.jl index f3c72ea..4ae7bbd 100644 --- a/src/soft_covers.jl +++ b/src/soft_covers.jl @@ -192,7 +192,8 @@ Supported penalty functions: 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]`. +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 @@ -421,8 +422,9 @@ matching [`soft_cover`](@ref). The row/column scales are pinned to the balance convention `∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]` (`nzaᵢ`, `nzbⱼ` = nonzero counts of row `i`, column -`j`), as in [`cover_min`](@ref): the objective depends on `a` and `b` only through the -products `a[i]*b[j]`, so without a convention the split between them would be arbitrary. +`j`), imposed within each connected component of the support, as in [`cover_min`](@ref): +the objective depends on `a` and `b` only through the products `a[i]*b[j]`, so without a +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 diff --git a/src/support.jl b/src/support.jl index 0c92486..2d93123 100644 --- a/src/support.jl +++ b/src/support.jl @@ -147,6 +147,60 @@ end # invariant already guarantees `abs(A[i,j]) == abs(A[j,i])`. require_abs_symmetric(::Union{Symmetric,Hermitian,Diagonal,SymTridiagonal}, fname) = nothing +# Connected components of the bipartite support graph of `A`: one vertex per row +# and one per column, one edge per stored nonzero. Returns `(rowcomp, colcomp, +# ncomp)`, where `rowcomp` and `colcomp` are `Vector{Int}` indexed by *position* +# within `axes(A, 1)` and `axes(A, 2)` (so offset axes need no special case, as +# with `GroupedSupport.ptr`), holding the component id in `1:ncomp` — or 0 for +# rows/columns with empty support, which belong to no component. +# +# The gauge orbit of an asymmetric cover has one dimension per component: the +# rescaling `a -> γ*a`, `b -> b/γ` acts independently on each, because no +# product `a[i]*b[j]` spans two components. Any convention that pins the split +# between `a` and `b` must therefore be imposed per component; a single global +# constraint leaves `ncomp - 1` directions to the whim of whichever pass ran +# last. +function _support_components(A::AbstractMatrix) + m = length(axes(A, 1)) + n = length(axes(A, 2)) + parent = collect(1:(m + n)) + touched = falses(m + n) + function find(p) + while parent[p] != p + parent[p] = parent[parent[p]] # path halving + p = parent[p] + end + return p + end + or = first(axes(A, 1)) - 1 + oc = first(axes(A, 2)) - 1 + foreach_support(A) do i, j, _ + p = i - or + q = m + j - oc + touched[p] = touched[q] = true + rp, rq = find(p), find(q) + rp == rq || (parent[rp] = rq) + end + label = zeros(Int, m + n) + ncomp = 0 + rowcomp = zeros(Int, m) + colcomp = zeros(Int, n) + for p in 1:(m + n) + touched[p] || continue + r = find(p) + if label[r] == 0 + ncomp += 1 + label[r] = ncomp + end + if p <= m + rowcomp[p] = label[r] + else + colcomp[p-m] = label[r] + end + end + return rowcomp, colcomp, ncomp +end + # Number of entries `foreach_support` reports, i.e. the size of the stored support. function _nsupport(A::AbstractMatrix) n = Ref(0) diff --git a/test/extensions.jl b/test/extensions.jl index a0349f5..5743382 100644 --- a/test/extensions.jl +++ b/test/extensions.jl @@ -525,6 +525,27 @@ end # the latter convention puts the optimum near [2.0, 3.177, 1.574] instead. Most # matrices do not discriminate, because the binding constraints have zero residual # at the optimum and a zero residual is weight-independent. +@testset "cover_min balances a block-diagonal support (JuMP/HiGHS)" begin + # Two connected components: the model constraint pins only the global gauge + # direction, so the per-component balance must come from the post-solve + # `_balance_cover!` shift `cover_min!` applies after the JuMP solve. + B = [2.0 1.0; 0.5 3.0] + C = [1.5 4.0 2.0; 3.0 0.5 1.0] + A = [B zeros(2, 3); zeros(2, 2) C] + a, b = cover_min(AbsLog{1}(), A) + @test iscover(a, b, A; atol=1e-8) + @test isbalanced(a, b, A) +end + +@testset "cover_min balances a block-diagonal support (Ipopt)" begin + B = [2.0 1.0; 0.5 3.0] + C = [1.5 4.0 2.0; 3.0 0.5 1.0] + A = [B zeros(2, 3); zeros(2, 2) C] + a, b = cover_min(AbsLinear{2}(), A) + @test iscover(a, b, A; atol=1e-8) + @test isbalanced(a, b, A) +end + @testset "the Ipopt sym objective uses the full-grid weighting" begin A = Float64[4 1 0; 1 1 5; 0 5 2] a0 = fill(sqrt(maximum(A)), 3) diff --git a/test/helpers.jl b/test/helpers.jl index 436e01b..810cf12 100644 --- a/test/helpers.jl +++ b/test/helpers.jl @@ -10,18 +10,28 @@ The balance convention `∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]` (`nzaᵢ`, `n nonzero counts of row `i` and column `j`), which fixes the row/column gauge `a → c*a`, `b → b/c`. The gauge leaves every product `a[i]*b[j]` unchanged, so no objective and no coverage constraint can see it; without a convention the split -between `a` and `b` would be an artifact of whichever pass last touched them. Every -asymmetric cover the package returns satisfies this. +between `a` and `b` would be an artifact of whichever pass last touched them. It is +imposed within each connected component of the bipartite support graph of `A` +separately, since the gauge acts independently on each; this checks the sum on +every component and returns `true` only if all of them balance. Every asymmetric +cover the package returns satisfies this. Note it is *not* scale-invariant: rescaling `A` moves the balance point, which is why [`covaries`](@ref) compares outer products rather than the vectors themselves. """ function isbalanced(a, b, A; atol=1e-8) - nza = vec(count(!iszero, A, dims=2)) - nzb = vec(count(!iszero, A, dims=1)) - La = sum(nza[i] * log(a[i]) for i in axes(A, 1) if nza[i] > 0; init=0.0) - Lb = sum(nzb[j] * log(b[j]) for j in axes(A, 2) if nzb[j] > 0; init=0.0) - return isapprox(La, Lb; atol=atol * max(1, abs(La), abs(Lb))) + rowcomp, colcomp, ncomp = MatrixCovers._support_components(A) + iszero(ncomp) && return true + La = zeros(ncomp) + Lb = zeros(ncomp) + or = first(axes(A, 1)) - 1 + oc = first(axes(A, 2)) - 1 + MatrixCovers.foreach_support(A) do i, j, v + c = rowcomp[i-or] + La[c] += log(a[i]) + Lb[c] += log(b[j]) + end + return all(isapprox(La[c], Lb[c]; atol=atol * max(1, abs(La[c]), abs(Lb[c]))) for c in 1:ncomp) end """ diff --git a/test/invariants.jl b/test/invariants.jl index 8d6e1a9..8a649d4 100644 --- a/test/invariants.jl +++ b/test/invariants.jl @@ -65,6 +65,8 @@ const GEN_NOTIONS = ( Azgen = [1.0 0.0 2.0; 0.0 0.0 0.0; 3.0 0.0 4.0] # row/column 2 unsupported Gc = [1.0+1.0im 2.0 0.5; 0.25im 3.0 1.0-2.0im] dr, dc = [2.0, 0.5], [3.0, 0.25, 1.5] + # Two connected components, so the balance convention pins two independent gauges. + Ablk = [1.0 2.0 0.0 0.0 0.0; 0.25 3.0 0.0 0.0 0.0; 0.0 0.0 1.5 2.5 0.75] @testset "gen: $(nt.name)" for nt in GEN_NOTIONS a, b = nt.f(Agen) @@ -90,11 +92,17 @@ const GEN_NOTIONS = ( # balance convention does, and every asymmetric cover reports its result in it. @test isbalanced(a, b, Agen) @test isbalanced(az, bz, Azgen) + # The balance convention is imposed per connected component: on a support with + # more than one component, `isbalanced` must hold within each rather than only + # in a single global sum. + ablk, bblk = nt.f(Ablk) + @test isbalanced(ablk, bblk, Ablk) end @testset "gen: initialize_cover($strategy, $feasible) is balanced" for strategy in (:hardcover, :geomean), feasible in (:inflate, :boost, :none) @test isbalanced(initialize_cover(Agen; strategy, feasible)..., Agen) + @test isbalanced(initialize_cover(Ablk; strategy, feasible)..., Ablk) end # The sym objective is summed over the full grid: each off-diagonal pair counts diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 31e837d..75810ca 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -165,6 +165,26 @@ end @test cover_objective(AbsLog{2}(), a, b, Matrix(D)) ≈ 0.0 atol = 1e-10 end +@testset "MCM block-diagonal assembly consistency" begin + # The balance convention is imposed per connected component, so the (a, b) split + # is a function of the support and the products alone: covering a block-diagonal + # matrix must reproduce the factors of covering each block separately. + rng = StableRNG(23) + for _ in 1:5 + mB, nB = rand(rng, 2:5), rand(rng, 2:5) + mC, nC = rand(rng, 2:5), rand(rng, 2:5) + B = rand(rng, mB, nB) .+ 0.1 + C = rand(rng, mC, nC) .+ 0.1 + A = [B zeros(mB, nC); zeros(mC, nB) C] + aA, bA = cover_min(AbsLog{2}(), A) + aB, bB = cover_min(AbsLog{2}(), B) + aC, bC = cover_min(AbsLog{2}(), C) + @test aA ≈ vcat(aB, aC) rtol=1e-5 + @test bA ≈ vcat(bB, bC) rtol=1e-5 + @test isbalanced(aA, bA, A) + end +end + @testset "symcover_min(AbsLog{2}) on complex Hermitian/Symmetric input" begin # The cover problem only depends on abs.(A), so a complex Hermitian (real # diagonal, conjugate-symmetric off-diagonals) or complex Symmetric matrix diff --git a/test/runtests.jl b/test/runtests.jl index c6d33d2..97b1303 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -45,7 +45,8 @@ include("helpers.jl") # isbalanced, covaries, PENALTIES :foreach_support, :foreach_support_sym, :cover_min_jump, :symcover_min_jump, :check_solved, :require_abs_symmetric, - :_edge_list, :_sym_edge_list, :_degrees) + :_edge_list, :_sym_edge_list, :_degrees, + :_balance_cover!, :inflate_feasible!) # Non-public names owned by other packages, each with no public equivalent: # `FreeUnits`/`Unit` are Unitful's unit representation and `Units` their # abstract supertype, needed to reject the unit types this package cannot diff --git a/test/support.jl b/test/support.jl index 2a66ad0..b437532 100644 --- a/test/support.jl +++ b/test/support.jl @@ -197,3 +197,62 @@ end @test SO.ax == -1:0 @test sort([(SO.idx[s], SO.val[s]) for s in MatrixCovers._slots(SO, -1)]) == [(-1, 2.0), (0, 1.0)] end + +# The gauge freedom of an asymmetric cover (`a -> γ*a`, `b -> b/γ`) acts independently +# on each connected component of the bipartite support graph, so the balance +# convention that pins it must be imposed per component; these tests check the +# labeling `_support_components` builds directly. +@testset "_support_components" begin + rng = StableRNG(11) + + @testset "dense: one component" begin + A = randn(rng, 4, 5) + rowcomp, colcomp, ncomp = MatrixCovers._support_components(A) + @test ncomp == 1 + @test all(==(1), rowcomp) + @test all(==(1), colcomp) + end + + @testset "block-diagonal: two components, correctly partitioned" begin + B = randn(rng, 3, 2) + C = randn(rng, 2, 4) + A = [B zeros(3, 4); zeros(2, 2) C] + rowcomp, colcomp, ncomp = MatrixCovers._support_components(A) + @test ncomp == 2 + # Rows/columns of B share one label, rows/columns of C the other, and the + # two labels differ. + @test allequal(rowcomp[1:3]) + @test allequal(colcomp[1:2]) + @test allequal(rowcomp[4:5]) + @test allequal(colcomp[3:6]) + @test rowcomp[1] == colcomp[1] + @test rowcomp[4] == colcomp[3] + @test rowcomp[1] != rowcomp[4] + end + + @testset "empty row and column labeled 0" begin + A = randn(rng, 4, 4) + A[2, :] .= 0 # empty row + A[:, 3] .= 0 # empty column + rowcomp, colcomp, ncomp = MatrixCovers._support_components(A) + @test rowcomp[2] == 0 + @test colcomp[3] == 0 + @test all(!=(0), rowcomp[[1, 3, 4]]) + @test all(!=(0), colcomp[[1, 2, 4]]) + # The rest of the support is still one connected component (row 2/col 3 + # aside, every other row touches every other column through a nonzero). + @test ncomp == 1 + end + + @testset "OffsetArray: position-indexed results match the parent" begin + B = randn(rng, 3, 2) + C = randn(rng, 2, 4) + A = [B zeros(3, 4); zeros(2, 2) C] + O = OffsetArray(A, -1:3, 10:15) + rowcomp, colcomp, ncomp = MatrixCovers._support_components(A) + orowcomp, ocolcomp, oncomp = MatrixCovers._support_components(O) + @test oncomp == ncomp + @test orowcomp == rowcomp + @test ocolcomp == colcomp + end +end