Add tests, improve README + docs

README.md

@@ -6,11 +6,16 @@
 [![Coverage](https://codecov.io/gh/HolyLab/MatrixCovers.jl/graph/badge.svg?token=trG4HXo9N4)](https://codecov.io/gh/HolyLab/MatrixCovers.jl)
 [![Aqua QA](https://juliatesting.github.io/Aqua.jl/dev/assets/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
-This package computes **covers** of matrices: non-negative vectors `a` (and `b`)
+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
-building block for scale-invariant numerical analysis.
+building block for scale-invariant numerical analysis. In particular, 
+``\hat A = A ./ (a * b^T)`` is scale-invariant, and because ``|\hat A[i, j]| <= 1``
+for all `i` and `j`, this simple construct finds applications that range from
+[statistical normalization](https://en.wikipedia.org/wiki/Normalization_(statistics))
+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

docs/src/index.md

@@ -60,9 +60,10 @@ julia> round.(typeof.(a), a; digits=6)
 
 `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 names each variable's natural scale outright, 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.
+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.
 
 Normalizing by the cover cancels the units along with the magnitudes, leaving a
 matrix that is all-ones, dimensionless, and scale-invariant:
@@ -80,18 +81,18 @@ It is worth noting that this yields 1 only for entries where the cover bound is
 
 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` naming the entries that conflict.  The requirement is not one
-this package adds: without it the terms in a row of `A*x` carry different units,
-so `A*x` is undefined for every `x`. If a matrix can be used in matrix-vector
-multiplication, it has a cover.
+`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`.
+If a matrix can be used in matrix-vector multiplication, it has a cover.
 
 ## Penalty functions
 
 A cover is valid as long as every constraint is satisfied, but tighter covers
 better capture the scaling of `A`.  Cover quality is measured through the ratios
-`r_{ij} = |A[i,j]| / (a[i] * b[j])`: a hard cover has every `r ≤ 1`, and `r = 1`
-means the constraint is exactly tight.  A **penalty function** `ϕ` combines those
-ratios into a scalar objective
+`r[i, j] = |A[i,j]| / (a[i] * b[j])`: a hard cover has every `0 ≤ r[i, j] ≤ 1`,
+and `r[i, j] == 1` means the constraint is exactly tight.
+
+A **penalty function** `ϕ` combines those ratios into a scalar objective
 
 ```math
 \sum_{i,j} \phi\!\left(\frac{|A_{ij}|}{a_i\, b_j}\right),
@@ -103,9 +104,9 @@ which the solvers minimize.  Two penalty families are provided:
   space, which makes them a favorable (and therefore default) penalty for *hard*
   covers, where `r ≤ 1` and `|log r|` is the log-excess of a constraint.
   `AbsLog{1}` sums the log-excesses (L1), `AbsLog{2}` sums their squares (L2).
-  Their principal disadvantage is the discontinuity at `r = 0`.
+  Their principal disadvantage is the divergence and discontinuity at `r = 0`.
 - [`AbsLinear`](@ref)`{p}` — `ϕ(r) = |1 - r|^p`.  Non-convex, but unlike
-  `AbsLog` these are continuous at `r = 0` (`ϕ(0) = 1`), so zero entries of `A`
+  `AbsLog` these are finite and continuous at `r = 0` (`ϕ(0) = 1`), so zero entries of `A`
   contribute a bounded penalty.  This is the penalty used by default for the
   *soft* covers, where `r > 1` (an uncovered entry) is allowed but penalized.
 

test/gram_covers.jl

@@ -17,6 +17,30 @@
         end
     end
 
+    @testset "one component: the global bound is attained" begin
+        # With a single support component there is nothing to accumulate separately, so the
+        # per-component construction returns exactly the global pair bound `norm(a)*b`. The
+        # agreement is two-sided: an inequality alone would also pass for a construction that
+        # gave up slack it did not have to.
+        #
+        # It is not bitwise, and the deviation is one-sided by design. Each `sqrt` is inflated by
+        # `1 + (n+3)*eps` so the cover holds despite naive summation; `norm` accumulates its own
+        # `n` terms with no such margin. So the tolerance is `(2n+3)*eps` — the declared inflation
+        # plus the reference's own roundoff — and `all(s .<= norm(a) .* b)` is *false* by a few
+        # ulps. That inequality is exposition, not a test.
+        for (seed, m, k) in ((3, 5, 4), (7, 4, 4), (11, 12, 3))
+            rng = StableRNG(seed)
+            J = randn(rng, m, k)
+            a, b = cover(J)
+            s = gramcover(a, b, J)
+            @test ncomponents(support_components(J)) == 1
+            @test s ≈ norm(a) .* b rtol = (2 * length(a) + 3) * eps(Float64)
+            @test !all(s .<= norm(a) .* b)
+            # Coverage, which the margin exists to deliver, is unconditional.
+            @test all(s * s' .>= abs.(J' * J))
+        end
+    end
+
     @testset "block-diagonal: per-component structure" begin
         rng = StableRNG(1)
         B = randn(rng, 4, 3)

test/runtests.jl

@@ -1,5 +1,6 @@
 using MatrixCovers
 using MatrixCovers: foreach_support, foreach_support_sym, unconstrained_min!, tighten_cover!
+using MatrixCovers: ncomponents, support_components
 using JuMP, HiGHS, Ipopt   # triggers MatrixCoversJuMPExt and MatrixCoversIpoptExt extensions
 using SparseArrays  # triggers MatrixCoversSparseArraysExt extension
 using Unitful       # triggers MatrixCoversUnitfulExt extension