Tighten the prose

README.md

@@ -16,12 +16,10 @@ for all `i` and `j`, this simple construct finds applications that range from
 of data to the design of well-behaved numerical algorithms (thanks, e.g., 
 to [bounds on $`\hat A`$'s eigenvalues](https://en.wikipedia.org/wiki/Gershgorin_circle_theorem)).
 
-Fast O(mn) heuristics (`symcover`, `cover`) are provided for everyday use, along
-with *soft* covers (`soft_symcover`, `soft_cover`) that penalize under-coverage
-rather than forbid it.  Objective-minimal hard covers (`symcover_min`,
-`cover_min`) minimize a penalty subject to the coverage constraint: the default
-squared-log-excess penalty is solved natively, with no external solver, while
-the other penalties are available when JuMP and HiGHS (or Ipopt) are loaded.
+The package provides O(mn) heuristics (`symcover`, `cover`), *soft* covers that
+penalize violations (`soft_symcover`, `soft_cover`), and objective-minimal hard
+covers (`symcover_min`, `cover_min`). The default squared-log penalty uses a
+built-in solver; other penalties use JuMP with HiGHS or Ipopt.
 
 ## Example
 
@@ -44,7 +42,7 @@ julia> iscover(a, A)
 true
 ```
 
-Covers are scale-covariant: rescaling the matrix rescales the cover the same way.
+Covers are scale-covariant:
 
 ```julia
 julia> D = Diagonal([10.0, 0.5]);
@@ -53,9 +51,8 @@ julia> symcover(D * A * D) ≈ D * a
 true
 ```
 
-Non-symmetric matrices get separate row and column scales from `cover`, and
-`cover_min`/`symcover_min` trade the fast heuristic for a cover that minimizes a
-penalty subject to the same constraint:
+For nonsymmetric matrices, `cover` returns separate row and column scales.
+`cover_min` minimizes a penalty subject to the coverage constraint:
 
 ```julia
 julia> M = [1.0 2.0 3.0; 6.0 5.0 4.0];
@@ -65,11 +62,11 @@ julia> a, b = cover(M);
 julia> iscover(a, b, M)
 true
 
-julia> aq, bq = cover_min(AbsLog{2}(), M);   # minimal, solved natively
+julia> aq, bq = cover_min(AbsLog{2}(), M);
 
 julia> cover_objective(AbsLog{2}(), aq, bq, M) <= cover_objective(AbsLog{2}(), a, b, M)
 true
 ```
 
-See the [documentation](https://HolyLab.github.io/MatrixCovers.jl/dev/)
-for motivation, examples, and a full API reference.
+See the [documentation](https://HolyLab.github.io/MatrixCovers.jl/dev/) for the
+algorithm guide and API reference.

ext/MatrixCoversIpoptExt.jl

@@ -6,38 +6,23 @@ using MatrixCovers
 using MatrixCovers: AbsLinear
 using MatrixCovers: _edge_list, _sym_edge_list, _degrees
 
-# The models are built over 1-based positions 1:n; `pr`/`pc` map each position to the
-# corresponding axis index of `A`, and results are scattered back onto vectors
-# whose axes match `A`'s so offset axes are honored. A row/column of `A` with no
-# nonzero entry appears in no constraint or objective term; its scale is set to
-# exactly 0, matching the native solvers.
-#
-# `A` is read through the support hook and gathered into a flat edge list in position
-# space — `ei`/`ej` the endpoints, `elog` the log-magnitude — so a model costs O(nnz)
-# to build rather than O(length(A)). The symmetric list is the full-grid reading; the
-# hard-cover models below sum over its `ei <= ej` half instead, per the objective each
-# one is defined by.
+# Models use 1-based positions and scatter results back to `A`'s axes. Support is
+# gathered as an O(nnz) edge list; unsupported scales are zero.
 
-# Ipopt returns a local minimum selected by the start. These kernels therefore
-# implement the mutating refiners; the main package supplies multistart drivers.
+# Ipopt kernels refine one start; the main package supplies multistart drivers.
 
 check_solved(model, fname) =
     MatrixCovers.check_solved(JuMP.termination_status(model), "Ipopt", fname)
 
-# `set_silent` alone still lets Ipopt print its startup banner, once per session, from
-# its C++ core; `sb` ("suppress banner") is the option that covers it. Every model here
-# is solved for a caller who asked for a cover, not for solver output.
+# Suppress both solver output and Ipopt's startup banner.
 function _ipopt_model()
     model = JuMP.Model(Ipopt.Optimizer)
     JuMP.set_silent(model)
     JuMP.set_attribute(model, "sb", "yes")
     return model
 end
 
-# The `i ≤ j` half of a symmetric gather, paired with the multiplicity `w` each entry
-# stands for in the full grid: an off-diagonal pair is two entries of `A`, a diagonal
-# entry one. Carrying the weight is equivalent to summing over both orientations and
-# costs half the terms — and, for the AbsLinear{1} models, half the auxiliary variables.
+# Symmetric triangle with full-grid multiplicities.
 function _triangle(fi, fj, flog)
     keep = [e for e in eachindex(fi) if fi[e] <= fj[e]]
     return fi[keep], fj[keep], flog[keep], [fi[e] == fj[e] ? 1 : 2 for e in keep]
@@ -112,16 +97,8 @@ end
 
 # ============================================================
 # Hard cover: cover_min!(::AbsLinear{p}, a, b, A)
-# The bipartite analog of symcover_min!: row scales α = log a, column scales
-# β = log b, residuals over every stored (i, j) rather than over i ≤ j. The product
-# 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. 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.
+# Asymmetric hard-cover model in row and column log scales. The model pins the
+# global gauge; post-processing balances components and restores feasibility.
 # ============================================================
 
 function MatrixCovers.cover_min!(::AbsLinear{2}, a::AbstractVector, b::AbstractVector, A)
@@ -196,9 +173,7 @@ end
 
 # ============================================================
 # Soft cover: soft_symcover_min!(::AbsLinear{p}, a, A)
-# Same objective, no coverage constraints — so the start need not cover `A`, and the
-# raw geometric mean (the exact soft AbsLog{2} optimum) is a natural one. The multistart
-# driver over these kernels is native; see soft_symcover_min.
+# Same objective without coverage constraints; starts need not cover `A`.
 # ============================================================
 
 function MatrixCovers.soft_symcover_min!(::AbsLinear{2}, a::AbstractVector, A)
@@ -255,14 +230,8 @@ 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 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.
+# Asymmetric soft-cover model. Post-processing balances component gauges; zero
+# entries contribute the constant `ϕ(0) = 1`.
 # ============================================================
 
 function MatrixCovers.soft_cover_min!(::AbsLinear{2}, a::AbstractVector, b::AbstractVector, A)

ext/MatrixCoversJuMPExt.jl

@@ -10,19 +10,10 @@ using LinearAlgebra: dot
 check_solved(model, fname) =
     MatrixCovers.check_solved(JuMP.termination_status(model), "HiGHS", fname)
 
-# The models are built over 1-based positions 1:n; `pr`/`pc` map each position to
-# the corresponding axis index of `A`, and results are scattered back onto vectors
-# whose axes match `A`'s so offset axes are honored. A row/column of `A` with no
-# nonzero entry carries no constraint or objective term; its scale is set to exactly
-# 0, matching the native solvers.
-#
-# `A` is read through the support hook and gathered into a flat edge list in position
-# space — `ei`/`ej` the endpoints, `elog` the log-magnitude — so a model costs O(nnz)
-# to build rather than O(length(A)). The symmetric list is the full-grid reading, whose
-# `ei <= ej` half is the constraint set.
+# Models use 1-based positions and scatter results back to `A`'s axes. Support is
+# gathered as an O(nnz) edge list; unsupported scales are zero.
 
-# Exact reference for the native `symcover_min(::AbsLog{2})`: same QP, solved by
-# HiGHS. Not exported; used by the test suite to cross-check the native solver.
+# HiGHS reference for tests of native symmetric `AbsLog{2}`.
 function MatrixCovers.symcover_min_jump(::AbsLog{2}, A)
     axr = axes(A, 1)
     axes(A, 2) == axr || throw(ArgumentError("symcover_min_jump requires a square matrix"))
@@ -56,16 +47,10 @@ function MatrixCovers.symcover_min!(::AbsLog{1}, a::AbstractVector, A)
     return a
 end
 
-# Relative slack allowed on the AbsLog{1} optimum while the AbsLog{2} objective is
-# minimized over it. The incumbent attains the bound exactly, so the face is never empty;
-# the slack only has to absorb the rounding of re-evaluating the objective row, and it
-# bounds how far the reported AbsLog{1} objective can drift above its true optimum.
+# Slack for re-evaluating the `AbsLog{1}` optimum during tie-breaking.
 const LEX_L1_SLACK = 1e-9
 
-# Select a unique point on the optimal AbsLog{1} face by minimizing AbsLog{2}
-# over it. Both objectives depend only on scale-invariant residuals.
-#
-# `residuals` uses the same support weighting as `cover_objective`.
+# Break `AbsLog{1}` ties with `AbsLog{2}` using full-grid support weights.
 function _minimize_l2_over_l1_face!(model, lin, residuals, fname)
     isempty(residuals) && return nothing
     linopt = JuMP.value(lin)
@@ -77,11 +62,8 @@ function _minimize_l2_over_l1_face!(model, lin, residuals, fname)
     return nothing
 end
 
-# The AbsLog{1} hard cover is an LP: the coverage constraint forces every residual
-# α[i]+α[j]-log|A[i,j]| to be nonnegative, so |·| drops away and the objective is
-# linear in α. Its optimum is a face, not a point, so a second stage picks the canonical
-# member of that face. `start`, when given, is a cover of `A` supplying the initial point;
-# it is a hint to the solver, and the canonical selection keeps it out of the result.
+# Symmetric `AbsLog{1}` LP. A second stage selects the canonical point on the
+# optimal face; `start` is only a solver hint.
 function _symcover_min_abslog1(A, start)
     axr = axes(A, 1)
     axes(A, 2) == axr || throw(ArgumentError("symcover_min requires a square matrix"))
@@ -136,10 +118,7 @@ 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.
+    # The post-solve balance handles component gauges not pinned here.
     @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")
@@ -165,11 +144,8 @@ function MatrixCovers.cover_min!(::AbsLog{1}, a::AbstractVector, b::AbstractVect
     return a, b
 end
 
-# Asymmetric counterpart of `_symcover_min_abslog1`, on the bipartite support: the
-# 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.
+# Asymmetric `AbsLog{1}` LP. Balance globally in the model and per component
+# after solving.
 function _cover_min_abslog1(A, start)
     axr = axes(A, 1)
     axc = axes(A, 2)
@@ -199,11 +175,7 @@ function _cover_min_abslog1(A, start)
         @constraint(model, α[ei[e]] + β[ej[e]] - elog[e] >= 0)
     end
     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. 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.
+    # Pin the global row/column gauge; post-processing handles components.
     @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")