Switch to inequality bounds

Project.toml

@@ -5,15 +5,34 @@ authors = ["Tim Holy <tim.holy@gmail.com> and contributors"]
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[weakdeps]
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+
+[extensions]
+SIAJuMP = ["HiGHS", "JuMP"]
+
 [compat]
+Graphs = "1"
+HiGHS = "1"
+JuMP = "1"
 LinearAlgebra = "1.10.0"
+NautyGraphs = "0.7"
+ParametricOptInterface = "0.15"
 SparseArrays = "1.10.0"
 julia = "1.10.10"
 
 [extras]
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+NautyGraphs = "7509a0a4-015a-4167-b44b-0799a1a2605e"
+ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["Graphs", "HiGHS", "JuMP", "NautyGraphs", "ParametricOptInterface", "Statistics", "Test"]

README.md

@@ -5,19 +5,16 @@
 [![Build Status](https://github.com/HolyLab/ScaleInvariantAnalysis.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/HolyLab/ScaleInvariantAnalysis.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/HolyLab/ScaleInvariantAnalysis.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/HolyLab/ScaleInvariantAnalysis.jl)
 
-This package provides tools for numerical analysis under conditions where the
-underlying mathematics is scale-invariant. We work on computers with finite
-precision, so operations like matrix-multiplication and matrix-division are
-expected to have some error. However, naive estimates of the error, based on
-quantities like `norm(x)`, may not be scale-invariant.
+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.
 
-For example, if `H` is a diagonal nonnegative (Hessian) matrix (i.e., a rank-2
-covariant tensor), with a change-of-scale in the variables all such `H` are
-equivalent to the identity matrix. Therefore we might think that its [condition
-number](https://en.wikipedia.org/wiki/Condition_number) should in fact be 1.
-This package provides tool to calculate condition number, as well as several
-other quantities, under conditions of scale-invariance.
+Fast O(n²) heuristics (`symcover`, `cover`) are provided for everyday use.
+Globally optimal covers minimizing a log-domain linear or quadratic objective
+(`symcover_lmin`, `cover_lmin`, `symcover_qmin`, `cover_qmin`) are available
+as an extension when JuMP and HiGHS are loaded.
 
-See the
-[documentation](https://HolyLab.github.io/ScaleInvariantAnalysis.jl/dev/) for
-more information.
+See the [documentation](https://HolyLab.github.io/ScaleInvariantAnalysis.jl/dev/)
+for motivation, examples, and a full API reference.

docs/src/index.md

@@ -4,58 +4,114 @@ CurrentModule = ScaleInvariantAnalysis
 
 # ScaleInvariantAnalysis
 
-This small package provides a number of tools for numerical analysis under
-conditions where the underlying problems are scale-invariant. At present it is
-oriented toward the types of problems that appear in mathematical optimization.
-Under scaling transformations ``x \rightarrow s \odot x`` (``\odot`` is the
-[Hadamard product](https://en.wikipedia.org/wiki/Hadamard_product_(matrices))),
-a Hessian matrix ``H`` transforms as ``H \rightarrow H \oslash (s \otimes s)``.
-Therefore, operations like `norm(H)` are non-sensical. Nevertheless, we work on
-computers with finite precision, so operations like ``Hx`` and ``H^{-1} g`` are
-expected to have some error. This package provides tools for calculating and
-estimating errors in a scale-covariant manner.
+This package computes **covers** of matrices.  Given a matrix `A`, a cover is a
+pair of non-negative vectors `a`, `b` satisfying
 
-## Example
+```math
+a_i \cdot b_j \;\geq\; |A_{ij}| \quad \text{for all } i, j.
+```
+
+For a symmetric matrix the cover is symmetric (`b = a`), so a single vector
+suffices: `a[i] * a[j] >= abs(A[i, j])`.
+
+## Why covers?
 
-Suppose you have a diagonal Hessian matrix and you estimate its condition number
-using tools that are not scale-invariant:
+Covers are the natural **scale-covariant** representation of a matrix.  If you
+rescale rows by a positive diagonal factor `D_r` and columns by `D_c`, the
+optimal cover transforms as `a → D_r * a`, `b → D_c * b` — exactly mirroring
+how the matrix entries change.  Scalar summaries like `norm(A)` or
+`maximum(abs, A)` do not have this property and therefore implicitly encode an
+arbitrary choice of units.
 
-```jldoctest example
-julia> using LinearAlgebra
+A concrete example: a 3×3 matrix whose rows and columns correspond to physical
+variables at very different scales (position in metres, velocity in m/s, force
+in N):
+
+```jldoctest coverones
+julia> using ScaleInvariantAnalysis
 
-julia> H = [1 0; 0 1e-8];
+julia> A = [1e6 1e3 1.0; 1e3 1.0 1e-3; 1.0 1e-3 1e-6];
 
-julia> cond(H)
-1.0e8
+julia> a = symcover(A)
+3-element Vector{Float64}:
+ 1000.0
+    1.0
+    0.001
 ```
 
-You might declare this matrix to be "poorly scaled." However, the operations
-`H * x` and `H \ g` both have coordinatewise relative errors of size `eps()`: there
-are no delicate cancelations and thus operations involving `H` reach the full
-machine precision. This does not seem entirely consistent with common
-expectations of working with matrices with large condition numbers.
+The cover `a` captures the natural per-variable scale.  The normalised matrix
+`A ./ (a .* a')` is all-ones and is scale-invariant.
 
-Under a coordinate transformation `x → [x[1], x[2]/10^4]`, `H` becomes the
-identity matrix which has a condition number of 1, and this better reflects our
-actual experience with operations involving `H`. This package provides a
-scale-invariant analog of the condition number:
+## Measuring cover quality
 
-```jldoctest example; filter = r"1\.0\d*" => "1.0"
+A cover is valid as long as every constraint is satisfied, but tighter covers
+better capture the scaling of `A`.  The *log-excess* of an entry is
+`log(a[i] * b[j] / abs(A[i, j])) >= 0`; it is zero when the constraint is
+exactly tight.  Two summary statistics aggregate these excesses:
+
+- [`lobjective`](@ref) — sum of log-excesses (L1 in log space).
+- [`qobjective`](@ref) — sum of squared log-excesses (L2 in log space).
+
+Both equal zero if and only if every constraint is tight.
+
+```jldoctest quality; filter = r"(\d+\.\d{6})\d+" => s"\1"
 julia> using ScaleInvariantAnalysis
 
-julia> condscale(H)
-1.0
+julia> A = [4.0 2.0; 2.0 9.0];
+
+julia> a = symcover(A)
+2-element Vector{Float64}:
+ 2.0
+ 3.0
+
+julia> lobjective(a, A)
+2.1972245773362196
+
+julia> qobjective(a, A)
+2.413897921625164
 ```
 
-(You may have some roundoff error in the last few digits.) This version of the
-condition number matches our actual experience using `H`. In contrast,
+## Choosing a cover algorithm
+
+| Function | Symmetric | Minimizes | Requires |
+|---|---|---|---|
+| [`symcover`](@ref) | yes | (fast heuristic) | — |
+| [`cover`](@ref) | no | (fast heuristic) | — |
+| [`symcover_lmin`](@ref) | yes | `lobjective` | JuMP + HiGHS |
+| [`cover_lmin`](@ref) | no | `lobjective` | JuMP + HiGHS |
+| [`symcover_qmin`](@ref) | yes | `qobjective` | JuMP + HiGHS |
+| [`cover_qmin`](@ref) | no | `qobjective` | JuMP + HiGHS |
+
+**`symcover` and `cover` are recommended for production use.**  They run in
+O(n²) time and often land within a few percent of the `lobjective`-optimal
+cover (see the quality tests in `test/testmatrices.jl`).
+
+The `*_lmin` and `*_qmin` variants solve a convex program (via
+[JuMP](https://jump.dev/) and [HiGHS](https://highs.dev/)) and return a
+global optimum of the respective objective.  They are loaded on demand as a
+package extension — simply load both libraries before calling them:
+
+```julia
+using JuMP, HiGHS
+using ScaleInvariantAnalysis
+
+a      = symcover_lmin(A)     # globally l-minimal symmetric cover
+a, b   = cover_qmin(A)        # globally q-minimal general cover
+```
+
+## Scale-invariant magnitude estimation
+
+[`divmag`](@ref) combines `symcover` with a right-hand side vector to produce a
+scale-covariant estimate of the magnitude of `A \ b` without solving the system:
 
-```jldoctest example; filter = r"(19999\.0\d*|19998\.9\d+)" => "19999.0"
-julia> condscale([1 0.9999; 0.9999 1])
-19999.0
+```julia
+a, mag = divmag(A, b)
 ```
 
-remains poorly-conditioned under all scale-transformations of the matrix.
+`a` is `symcover(A)` and `mag` estimates `dotabs(A \ b, a)`.  Both transform
+covariantly when `A` and `b` are rescaled together, so `mag` serves as a
+reliable unit for assessing roundoff in the solution.  Pass `cond=true` to fold
+in the scale-invariant condition number for ill-conditioned systems.
 
 ## Index of available tools
 

ext/SIAJuMP.jl

@@ -0,0 +1,91 @@
+module SIAJuMP
+
+using JuMP: JuMP, @variable, @objective, @constraint
+using HiGHS: HiGHS
+using ScaleInvariantAnalysis
+
+# Reference implementation for a symmetric matrix cover
+function ScaleInvariantAnalysis.symcover_qmin(A)
+    ax = axes(A, 1)
+    axes(A, 2) == ax || throw(ArgumentError("symcover_qmin requires a square matrix"))
+    logA = log.(abs.(A))
+    n = size(logA, 1)
+    model = JuMP.Model(HiGHS.Optimizer)
+    JuMP.set_silent(model)
+    @variable(model, α[1:n])
+    @objective(model, Min, sum(abs2, α[i] + α[j] - logA[i, j] for i in 1:n, j in 1:n if A[i, j] != 0))
+    for i in 1:n
+        for j in i:n
+            if A[i, j] != 0
+                @constraint(model, α[i] + α[j] - logA[i, j] >= 0)
+            end
+        end
+    end
+    JuMP.optimize!(model)
+    return exp.(JuMP.value.(α))
+end
+
+function ScaleInvariantAnalysis.symcover_lmin(A)
+    ax = axes(A, 1)
+    axes(A, 2) == ax || throw(ArgumentError("symcover_lmin_ref requires a square matrix"))
+    logA = log.(abs.(A))
+    n = size(logA, 1)
+    model = JuMP.Model(HiGHS.Optimizer)
+    JuMP.set_silent(model)
+    @variable(model, α[1:n])
+    @objective(model, Min, sum(α[i] + α[j] - logA[i, j] for i in 1:n, j in 1:n if A[i, j] != 0))
+    for i in 1:n
+        for j in i:n
+            if A[i, j] != 0
+                @constraint(model, α[i] + α[j] - logA[i, j] >= 0)
+            end
+        end
+    end
+    JuMP.optimize!(model)
+    return exp.(JuMP.value.(α))
+end
+
+# Reference implementation for a general matrix cover
+function ScaleInvariantAnalysis.cover_qmin(A)
+    logA = log.(abs.(A))
+    m, n = size(logA)
+    model = JuMP.Model(HiGHS.Optimizer)
+    JuMP.set_silent(model)
+    @variable(model, α[1:m])
+    @variable(model, β[1:n])
+    @objective(model, Min, sum(abs2, α[i] + β[j] - logA[i, j] for i in 1:m, j in 1:n if A[i, j] != 0))
+    for i in 1:m
+        for j in 1:n
+            if A[i, j] != 0
+                @constraint(model, α[i] + β[j] - logA[i, j] >= 0)
+            end
+        end
+    end
+    nza, nzb = sum(!iszero, A; dims=2)[:], sum(!iszero, A; dims=1)[:]
+    @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n))
+    JuMP.optimize!(model)
+    return exp.(JuMP.value.(α)), exp.(JuMP.value.(β))
+end
+
+function ScaleInvariantAnalysis.cover_lmin(A)
+    logA = log.(abs.(A))
+    m, n = size(logA)
+    model = JuMP.Model(HiGHS.Optimizer)
+    JuMP.set_silent(model)
+    @variable(model, α[1:m])
+    @variable(model, β[1:n])
+    @objective(model, Min, sum(α[i] + β[j] - logA[i, j] for i in 1:m, j in 1:n if A[i, j] != 0))
+    for i in 1:m
+        for j in 1:n
+            if A[i, j] != 0
+                @constraint(model, α[i] + β[j] - logA[i, j] >= 0)
+            end
+        end
+    end
+    nza, nzb = sum(!iszero, A; dims=2)[:], sum(!iszero, A; dims=1)[:]
+    @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n))
+    JuMP.optimize!(model)
+    return exp.(JuMP.value.(α)), exp.(JuMP.value.(β))
+end
+
+end