diff --git a/README.md b/README.md index 0ab7b0c..554c233 100644 --- a/README.md +++ b/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 diff --git a/docs/src/index.md b/docs/src/index.md index 65440c0..0cb24cd 100644 --- a/docs/src/index.md +++ b/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. diff --git a/test/gram_covers.jl b/test/gram_covers.jl index dd38a20..864fe01 100644 --- a/test/gram_covers.jl +++ b/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) diff --git a/test/runtests.jl b/test/runtests.jl index e2bbdf0..00412e0 100644 --- a/test/runtests.jl +++ b/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