Remove dotabs and divmag

docs/src/index.md

@@ -280,6 +280,104 @@ perturbations of a base point up to a user-controllable number of `starts`.
 For the convex `AbsLog` penalties the start cannot change the result, and the refiners
 accept one only so that the two families share an interface.
 
+## Worked example: roundoff in `A \ b`
+
+Because a cover names each variable's natural scale, it also says how to measure a
+solution in units that do not depend on how the problem was parameterized.
+
+Solving `x = A \ b` is *contravariant*: rescaling `A → D*A*D` and `b → D*b` sends
+`x → x ./ d`, while the cover is covariant, `a → d .* a`.  The products `x .* a` are
+therefore unchanged, and `∑ᵢ |xᵢ * aᵢ|` is a measure of the solution's size that is the
+same in every frame.
+
+That quantity can be estimated from the magnitudes of `A` and `b` alone, without
+forming `x` at all:
+
+```jldoctest roundoff
+julia> using ScaleInvariantAnalysis, LinearAlgebra
+
+julia> A = [1e6 1e3; 1e3 4.0];
+
+julia> b = [1.5e3, 6.0];
+
+julia> a = symcover(A);
+
+julia> round.(a; digits=6)
+2-element Vector{Float64}:
+ 1000.0
+    2.0
+
+julia> mag = sum(abs(bi / ai) for (bi, ai) in zip(b, a))
+4.5
+```
+
+The cover reports natural scales of 1000 and 2, and `mag` estimates the size of the
+solution measured against them — here within a factor of 1.5 of the truth:
+
+```jldoctest roundoff
+julia> x = A \ b;
+
+julia> sum(abs.(x .* a))
+3.0
+```
+
+Both numbers are scale-invariant, so the estimate is unchanged by any diagonal
+rescaling of the problem:
+
+```jldoctest roundoff
+julia> d = [0.05, 3.0];
+
+julia> Ad, bd = d .* A .* d', d .* b;
+
+julia> ad = symcover(Ad);
+
+julia> sum(abs(bi / ai) for (bi, ai) in zip(bd, ad))
+4.5
+```
+
+This makes `eps(mag)` a scale-invariant estimate of the roundoff floor of the sum.
+For a well-conditioned `A`, the error of the `Float64` solve meets that floor:
+
+```jldoctest roundoff
+julia> xbig = big.(A) \ big.(b);
+
+julia> abs(sum(abs.(x .* a)) - sum(abs.(Float64.(xbig) .* a))) <= 2 * eps(mag)
+true
+```
+
+The estimate is built from magnitudes only, so it knows nothing about the conditioning
+of `A` or about cancellation during the solve.  When `A` is ill-conditioned the true
+error sits far above the floor:
+
+```jldoctest roundoff
+julia> Aill = [1.0 -0.9999; -0.9999 1.0];
+
+julia> bill = [0.75, 7.0];
+
+julia> aill = symcover(Aill);
+
+julia> magill = sum(abs(bi / ai) for (bi, ai) in zip(bill, aill));
+
+julia> xill = Aill \ bill;
+
+julia> xbigill = big.(Aill) \ big.(bill);
+
+julia> err = abs(sum(abs.(xill .* aill)) - sum(abs.(Float64.(xbigill) .* aill)));
+
+julia> err > 1e6 * eps(magill)
+true
+```
+
+Folding in the condition number of the *normalized* matrix `A ./ (a .* a')` — itself
+scale-invariant, since normalizing cancels the frame — restores a usable bound:
+
+```jldoctest roundoff
+julia> κ = cond(Aill ./ (aill .* aill'));
+
+julia> err <= 1e3 * eps(κ * magill)
+true
+```
+
 ## Index of available tools
 
 ```@index

src/ScaleInvariantAnalysis.jl

@@ -11,7 +11,6 @@ export cover, cover!, symcover, symcover!, soft_symcover, soft_cover
 export initialize_cover, initialize_cover!, initialize_symcover, initialize_symcover!
 export symcover_min, symcover_min!, cover_min, cover_min!
 export soft_symcover_min, soft_symcover_min!, soft_cover_min, soft_cover_min!
-export dotabs, divmag
 
 include("penalties.jl")
 include("support.jl")
@@ -21,42 +20,6 @@ include("soft_covers.jl")
 include("minimal_covers.jl")
 
 
-"""
-    dotabs(x, y)
-
-Compute the sum of absolute values of elementwise products of `x` and `y`:
-
-    ∑_i |x[i] * y[i]|
-"""
-function dotabs(x::AbstractVector, y::AbstractVector)
-    s = zero(eltype(x)) * zero(eltype(y))
-    for i in eachindex(x, y)
-        s += abs(x[i] * y[i])
-    end
-    return s
-end
-
-"""
-    a, mag = divmag(A, b; use_cond::Bool=false)
-
-Given a symmetric matrix `A` and vector `b`, for `x = A \\ b` return a pair
-where `mag` is a naive estimate of the magnitude of `sum(abs.(x .* a))`. `a` and
-`mag` are scale-covariant in circumstances where `A \\ b` is contravariant. With
-`use_cond=false`, the estimate is based only on the magnitudes of the numbers
-in `A` and `b`, and does not account for the conditioning of `A` or
-cancellation in the solution process.
-
-This can be used to form scale-invariant estimates of roundoff errors in
-computations involving `A`, `b`, and `x`.
-"""
-function divmag(A, b; use_cond::Bool=false)
-    a = symcover(A)
-    κ = use_cond ? LinearAlgebra.cond(A ./ (a .* a')) : 1
-    return a, κ * sum(abs ∘ splat(ratio_nz), zip(b, a))
-end
-ratio_nz(n, d) = iszero(d) ? zero(n) / oneunit(d) : n / d
-
-
 # True only when a MethodError's argument types are consistent with the calling
 # convention of the `*_min` solvers — a penalty, the scale vectors the mutating
 # forms refine in place, and the matrix — i.e. the failure could plausibly be fixed

test/penalties.jl

@@ -1,4 +1,4 @@
-# Penalty objectives (cover_objective) and the scalar utilities dotabs/divmag.
+# Penalty objectives (cover_objective).
 
 @testset "cover_objective" begin
     A = [4.0 1.5; 1.5 1.0]
@@ -35,34 +35,3 @@ end
     end
 end
 
-@testset "dotabs" begin
-    @test dotabs([1.0, -2.0, 3.0], [4.0, 5.0, -6.0]) ≈ 4.0 + 10.0 + 18.0
-    @test dotabs([0.0, 1.0], [1.0, 0.0]) == 0.0
-    @test dotabs(big.([1.0, 2.0]), big.([3.0, 4.0])) ≈ 3.0 + 8.0
-    @test dotabs([1.0 + 2.0im, -1.0im], [3.0, 1.0 + 1.0im]) ≈ abs((1.0 + 2.0im) * 3.0) + abs(-1.0im * (1.0 + 1.0im))
-end
-
-@testset "divmag" begin
-    A = [1.0 -0.2; -0.2 0]
-    b = [0.75, 7.0]
-    a, mag = divmag(A, b)
-    @test abs(dotabs(A \ b, a) - dotabs(big.(A) \ big.(b), a)) <= 2 * eps(mag)
-    # Uniform scaling covariance
-    asc, magsc = divmag(1000 * A, 1000 * b)
-    @test asc ≈ sqrt(1000) .* a
-    @test magsc ≈ sqrt(1000) * mag
-    # Diagonal scaling covariance, from mild to strong anisotropy
-    for d in ([1.0, 1.0], [0.5, 2.0], [0.05, 3.0], [10.0, 0.01])
-        Ad = A .* d .* d'
-        bd = b .* d
-        a_d, mag_d = divmag(Ad, bd)
-        @test abs(dotabs(Ad \ bd, a_d) - dotabs(big.(Ad) \ big.(bd), a_d)) <= 100 * eps(mag_d)
-    end
-    # Ill-conditioned matrix
-    A_ill = [1.0 -0.9999; -0.9999 1]
-    a_ill, mag_ill = divmag(A_ill, b)
-    @test abs(dotabs(A_ill \ b, a_ill) - dotabs(big.(A_ill) \ big.(b), a_ill)) > 10^6 * eps(mag_ill)
-    # With use_cond=true, accuracy is recovered
-    a_ill2, mag_ill2 = divmag(A_ill, b; use_cond=true)
-    @test abs(dotabs(A_ill \ b, a_ill2) - dotabs(big.(A_ill) \ big.(b), a_ill2)) <= 10^3 * eps(mag_ill2)
-end

test/runtests.jl

@@ -1,5 +1,5 @@
 using ScaleInvariantAnalysis
-using ScaleInvariantAnalysis: divmag, dotabs, foreach_support, foreach_support_sym, unconstrained_min!, tighten_cover!
+using ScaleInvariantAnalysis: foreach_support, foreach_support_sym, unconstrained_min!, tighten_cover!
 using JuMP, HiGHS, Ipopt   # triggers SIAJuMP and SIAIpopt extensions
 using SparseArrays  # triggers SIASparseArrays extension
 using Unitful       # triggers SIAUnitful extension
@@ -16,7 +16,7 @@ include("helpers.jl")               # iscover, covaries, PENALTIES
 
 @testset "ScaleInvariantAnalysis.jl" begin
 
-    include("penalties.jl")         # cover_objective, dotabs, divmag
+    include("penalties.jl")         # cover_objective
     include("support.jl")           # foreach_support(_sym) traversal
     include("heuristic_covers.jl")  # symcover/cover and their internals
     include("soft_covers.jl")       # soft_symcover/soft_cover multistart descent