Pin the asymmetric-cover gauge per support component

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

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

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
 
 

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

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).
 """

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).