Minimal gram covers

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

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)

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)

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,20 +73,18 @@ 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
 `∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]` (`nzaᵢ`, `nzbⱼ` = nonzero counts of row `i`,
 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

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.

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

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