@@ -4,58 +4,114 @@ CurrentModule = ScaleInvariantAnalysis
44
55# ScaleInvariantAnalysis
66
7- This small package provides a number of tools for numerical analysis under
8- conditions where the underlying problems are scale-invariant. At present it is
9- oriented toward the types of problems that appear in mathematical optimization.
10- Under scaling transformations `` x \rightarrow s \odot x `` (`` \odot `` is the
11- [ Hadamard product] ( https://en.wikipedia.org/wiki/Hadamard_product_(matrices) ) ),
12- a Hessian matrix `` H `` transforms as `` H \rightarrow H \oslash (s \otimes s) `` .
13- Therefore, operations like ` norm(H) ` are non-sensical. Nevertheless, we work on
14- computers with finite precision, so operations like `` Hx `` and `` H^{-1} g `` are
15- expected to have some error. This package provides tools for calculating and
16- estimating errors in a scale-covariant manner.
7+ This package computes ** covers** of matrices. Given a matrix ` A ` , a cover is a
8+ pair of non-negative vectors ` a ` , ` b ` satisfying
179
18- ## Example
10+ ``` math
11+ a_i \cdot b_j \;\geq\; |A_{ij}| \quad \text{for all } i, j.
12+ ```
13+
14+ For a symmetric matrix the cover is symmetric (` b = a ` ), so a single vector
15+ suffices: ` a[i] * a[j] >= abs(A[i, j]) ` .
16+
17+ ## Why covers?
1918
20- Suppose you have a diagonal Hessian matrix and you estimate its condition number
21- using tools that are not scale-invariant:
19+ Covers are the natural ** scale-covariant** representation of a matrix. If you
20+ rescale rows by a positive diagonal factor ` D_r ` and columns by ` D_c ` , the
21+ optimal cover transforms as ` a → D_r * a ` , ` b → D_c * b ` — exactly mirroring
22+ how the matrix entries change. Scalar summaries like ` norm(A) ` or
23+ ` maximum(abs, A) ` do not have this property and therefore implicitly encode an
24+ arbitrary choice of units.
2225
23- ``` jldoctest example
24- julia> using LinearAlgebra
26+ A concrete example: a 3×3 matrix whose rows and columns correspond to physical
27+ variables at very different scales (position in metres, velocity in m/s, force
28+ in N):
29+
30+ ``` jldoctest coverones
31+ julia> using ScaleInvariantAnalysis
2532
26- julia> H = [1 0; 0 1e-8 ];
33+ julia> A = [1e6 1e3 1.0; 1e3 1.0 1e-3; 1. 0 1e-3 1e-6 ];
2734
28- julia> cond(H)
29- 1.0e8
35+ julia> a = symcover(A)
36+ 3-element Vector{Float64}:
37+ 1000.0
38+ 1.0
39+ 0.001
3040```
3141
32- You might declare this matrix to be "poorly scaled." However, the operations
33- ` H * x ` and ` H \ g ` both have coordinatewise relative errors of size ` eps() ` : there
34- are no delicate cancelations and thus operations involving ` H ` reach the full
35- machine precision. This does not seem entirely consistent with common
36- expectations of working with matrices with large condition numbers.
42+ The cover ` a ` captures the natural per-variable scale. The normalised matrix
43+ ` A ./ (a .* a') ` is all-ones and is scale-invariant.
3744
38- Under a coordinate transformation ` x → [x[1], x[2]/10^4] ` , ` H ` becomes the
39- identity matrix which has a condition number of 1, and this better reflects our
40- actual experience with operations involving ` H ` . This package provides a
41- scale-invariant analog of the condition number:
45+ ## Measuring cover quality
4246
43- ``` jldoctest example; filter = r"1\. 0\d*" => "1.0"
47+ A cover is valid as long as every constraint is satisfied, but tighter covers
48+ better capture the scaling of ` A ` . The * log-excess* of an entry is
49+ ` log(a[i] * b[j] / abs(A[i, j])) >= 0 ` ; it is zero when the constraint is
50+ exactly tight. Two summary statistics aggregate these excesses:
51+
52+ - [ ` lobjective ` ] ( @ref ) — sum of log-excesses (L1 in log space).
53+ - [ ` qobjective ` ] ( @ref ) — sum of squared log-excesses (L2 in log space).
54+
55+ Both equal zero if and only if every constraint is tight.
56+
57+ ``` jldoctest quality; filter = r"(\d+\. \d{6})\d+" => s"\1"
4458julia> using ScaleInvariantAnalysis
4559
46- julia> condscale(H)
47- 1.0
60+ julia> A = [4.0 2.0; 2.0 9.0];
61+
62+ julia> a = symcover(A)
63+ 2-element Vector{Float64}:
64+ 2.0
65+ 3.0
66+
67+ julia> lobjective(a, A)
68+ 2.1972245773362196
69+
70+ julia> qobjective(a, A)
71+ 2.413897921625164
4872```
4973
50- (You may have some roundoff error in the last few digits.) This version of the
51- condition number matches our actual experience using ` H ` . In contrast,
74+ ## Choosing a cover algorithm
75+
76+ | Function | Symmetric | Minimizes | Requires |
77+ | ---| ---| ---| ---|
78+ | [ ` symcover ` ] ( @ref ) | yes | (fast heuristic) | — |
79+ | [ ` cover ` ] ( @ref ) | no | (fast heuristic) | — |
80+ | [ ` symcover_lmin ` ] ( @ref ) | yes | ` lobjective ` | JuMP + HiGHS |
81+ | [ ` cover_lmin ` ] ( @ref ) | no | ` lobjective ` | JuMP + HiGHS |
82+ | [ ` symcover_qmin ` ] ( @ref ) | yes | ` qobjective ` | JuMP + HiGHS |
83+ | [ ` cover_qmin ` ] ( @ref ) | no | ` qobjective ` | JuMP + HiGHS |
84+
85+ ** ` symcover ` and ` cover ` are recommended for production use.** They run in
86+ O(n²) time and often land within a few percent of the ` lobjective ` -optimal
87+ cover (see the quality tests in ` test/testmatrices.jl ` ).
88+
89+ The ` *_lmin ` and ` *_qmin ` variants solve a convex program (via
90+ [ JuMP] ( https://jump.dev/ ) and [ HiGHS] ( https://highs.dev/ ) ) and return a
91+ global optimum of the respective objective. They are loaded on demand as a
92+ package extension — simply load both libraries before calling them:
93+
94+ ``` julia
95+ using JuMP, HiGHS
96+ using ScaleInvariantAnalysis
97+
98+ a = symcover_lmin (A) # globally l-minimal symmetric cover
99+ a, b = cover_qmin (A) # globally q-minimal general cover
100+ ```
101+
102+ ## Scale-invariant magnitude estimation
103+
104+ [ ` divmag ` ] ( @ref ) combines ` symcover ` with a right-hand side vector to produce a
105+ scale-covariant estimate of the magnitude of ` A \ b ` without solving the system:
52106
53- ``` jldoctest example; filter = r"(19999\. 0\d*|19998\. 9\d+)" => "19999.0"
54- julia> condscale([1 0.9999; 0.9999 1])
55- 19999.0
107+ ``` julia
108+ a, mag = divmag (A, b)
56109```
57110
58- remains poorly-conditioned under all scale-transformations of the matrix.
111+ ` a ` is ` symcover(A) ` and ` mag ` estimates ` dotabs(A \ b, a) ` . Both transform
112+ covariantly when ` A ` and ` b ` are rescaled together, so ` mag ` serves as a
113+ reliable unit for assessing roundoff in the solution. Pass ` cond=true ` to fold
114+ in the scale-invariant condition number for ill-conditioned systems.
59115
60116## Index of available tools
61117
0 commit comments