Add symdiagcover

docs/src/index.md

@@ -5,15 +5,20 @@ CurrentModule = ScaleInvariantAnalysis
 # ScaleInvariantAnalysis
 
 This package computes **covers** of matrices.  Given a matrix `A`, a cover is a
-pair of non-negative vectors `a`, `b` satisfying
+matrix `C` that can be defined as `C = a * b'`, where `a` and `b` are non-negative vectors.
+`C` must satisfy
 
 ```math
-a_i \cdot b_j \;\geq\; |A_{ij}| \quad \text{for all } i, j.
+C_{ij} \;\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])`.
 
+For symmetric matrices, this package also supports **diagonalized covers**,
+where `C = Diagonal(d) + a * a'` is a cover for `A`. Here, both `a` and `d` are
+nonnegative vectors.
+
 ## Why covers?
 
 Covers are the natural **scale-covariant** representation of a matrix.  If you
@@ -101,14 +106,16 @@ julia> aq * bq'
 |---|---|---|---|
 | [`symcover`](@ref) | yes | (fast heuristic) | — |
 | [`cover`](@ref) | no | (fast heuristic) | — |
+| [`symdiagcover`](@ref) | yes | (fast heuristic) | — |
 | [`symcover_lmin`](@ref) | yes | `cover_lobjective` | JuMP + HiGHS |
 | [`cover_lmin`](@ref) | no | `cover_lobjective` | JuMP + HiGHS |
 | [`symcover_qmin`](@ref) | yes | `cover_qobjective` | JuMP + HiGHS |
 | [`cover_qmin`](@ref) | no | `cover_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 `cover_lobjective`-optimal
-cover (see the quality tests involving `test/testmatrices.jl`).
+**`symcover`, `cover`, and `symdiagcover` are recommended for production use.**
+They run in O(mn) time for an ``m\times n`` matrix `A` and often land within a few
+percent of the `cover_lobjective`-optimal cover (see the quality tests involving
+`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
@@ -127,7 +134,6 @@ a, b   = cover_qmin(A)        # globally quadratic-minimal general cover
 
 ```@index
 Modules = [ScaleInvariantAnalysis]
-Private = false
 ```
 
 ## Reference documentation

src/ScaleInvariantAnalysis.jl

@@ -5,7 +5,7 @@ using SparseArrays
 using LoopVectorization
 using PrecompileTools
 
-export cover_lobjective, cover_qobjective, symcover, cover, symcover_lmin, cover_lmin, symcover_qmin, cover_qmin
+export cover_lobjective, cover_qobjective, cover, symcover, symdiagcover, cover_lmin, symcover_lmin, cover_qmin, symcover_qmin
 export dotabs
 
 include("covers.jl")

src/covers.jl

@@ -57,12 +57,16 @@ julia> a * a'
  1.0  0.25
 ```
 """
-function symcover(A::AbstractMatrix; kwargs...)
+function symcover(A::AbstractMatrix; exclude_diagonal::Bool=false, kwargs...)
     ax = axes(A, 1)
     axes(A, 2) == ax || throw(ArgumentError("symcover requires a square matrix"))
     a = similar(A, float(eltype(A)), ax)
-    for j in ax
-        a[j] = sqrt(abs(A[j, j]))
+    if exclude_diagonal
+        fill!(a, zero(eltype(a)))
+    else
+        for j in ax
+            a[j] = sqrt(abs(A[j, j]))
+        end
     end
     # Iterate over the diagonals of A, and update a[i] and a[j] to satisfy |A[i, j]| ≤ a[i] * a[j] whenever this constraint is violated
     # Iterating over the diagonals gives a more "balanced" result and typically results in lower loss than iterating in a triangular pattern.
@@ -87,7 +91,53 @@ function symcover(A::AbstractMatrix; kwargs...)
             end
         end
     end
-    return tighten_cover!(a, A; kwargs...)
+    return tighten_cover!(a, A; exclude_diagonal, kwargs...)
+end
+
+"""
+    d, a = symdiagcover(A; kwargs...)
+
+Given a square matrix `A` assumed to be symmetric, return vectors `d` and `a`
+representing a symmetric diagonalized cover `Diagonal(d) + a * a'` of `A` with
+the diagonal as tight as possible given `A` and `a`. In particular,
+
+    abs(A[i, j]) ≤ a[i] * a[j] for all i ≠ j, and
+    abs(A[i, i]) ≤ a[i]^2 + d[i] for all i.
+
+# Examples
+
+```jldoctest; setup=:(using LinearAlgebra), filter = r"(\\d+\\.\\d{6})\\d+" => s"\\1"
+julia> A = [4 1e-8; 1e-8 1];
+
+julia> a = symcover(A)
+2-element Vector{Float64}:
+ 2.0
+ 1.0
+
+julia> a * a'
+2×2 Matrix{Float64}:
+ 4.0  2.0
+ 2.0  1.0
+
+julia> d, a = symdiagcover(A)
+([3.99999999, 0.99999999], [0.0001, 0.0001])
+
+julia> Diagonal(d) + a * a'
+2×2 Matrix{Float64}:
+ 4.0     1.0e-8
+ 1.0e-8  1.0
+```
+For this case, one sees much tighter covering with `symdiagcover`.
+"""
+function symdiagcover(A::AbstractMatrix; kwargs...)
+    ax = axes(A, 1)
+    axes(A, 2) == ax || throw(ArgumentError("symcover requires a square matrix"))
+    a = symcover(A; exclude_diagonal=true, kwargs...)
+    d = map(ax) do i
+        Aii, ai = abs(A[i, i]), a[i]
+        max(zero(Aii), Aii - ai^2)
+    end
+    return d, a
 end
 
 """
@@ -161,22 +211,41 @@ function cover(A::AbstractMatrix; kwargs...)
     return tighten_cover!(a, b, A; kwargs...)
 end
 
-function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; iter::Int=3) where T
+function tighten_cover!(a::AbstractVector{T}, A::AbstractMatrix; iter::Int=3, exclude_diagonal::Bool=false) where T
     ax = axes(A, 1)
     axes(A, 2) == ax || throw(ArgumentError("`tighten_cover!(a, A)` requires a square matrix `A`"))
     eachindex(a) == ax || throw(DimensionMismatch("indices of `a` must match the indexing of `A`"))
     aratio = similar(a)
     for _ in 1:iter
         fill!(aratio, typemax(T))
-        @turbo for j in eachindex(a)
-            aratioj, aj = aratio[j], a[j]
-            for i in eachindex(a)
-                Aij = T(abs(A[i, j]))
-                r = ifelse(iszero(Aij), typemax(T), a[i] * aj / Aij)
-                aratio[i] = min(aratio[i], r)
-                aratioj = min(aratioj, r)
+        if exclude_diagonal
+            for j in eachindex(a)
+                aratioj, aj = aratio[j], a[j]
+                for i in first(ax):j-1
+                    Aij = T(abs(A[i, j]))
+                    r = ifelse(iszero(Aij), typemax(T), a[i] * aj / Aij)
+                    aratio[i] = min(aratio[i], r)
+                    aratioj = min(aratioj, r)
+                end
+                for i in j+1:last(ax)
+                    Aij = T(abs(A[i, j]))
+                    r = ifelse(iszero(Aij), typemax(T), a[i] * aj / Aij)
+                    aratio[i] = min(aratio[i], r)
+                    aratioj = min(aratioj, r)
+                end
+                aratio[j] = aratioj
+            end
+        else
+            @turbo for j in eachindex(a)
+                aratioj, aj = aratio[j], a[j]
+                for i in eachindex(a)
+                    Aij = T(abs(A[i, j]))
+                    r = ifelse(iszero(Aij), typemax(T), a[i] * aj / Aij)
+                    aratio[i] = min(aratio[i], r)
+                    aratioj = min(aratioj, r)
+                end
+                aratio[j] = aratioj
             end
-            aratio[j] = aratioj
         end
         a ./= sqrt.(aratio)
     end

test/runtests.jl

@@ -51,6 +51,37 @@ using Test
         @test b / c ≈ dc .* b0
     end
 
+    @testset "symdiagcover" begin
+        for A in ([2.0 1.0; 1.0 3.0], [1.0 -0.2; -0.2 0.0], [4.0 1e-8; 1e-8 1.0],
+                  [100.0 1.0; 1.0 0.01], [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0])
+            d, a = symdiagcover(A)
+            # Off-diagonal cover property
+            @test all(a[i] * a[j] >= abs(A[i, j]) - 1e-12 for i in axes(A, 1), j in axes(A, 2) if i != j)
+            # Diagonal cover property
+            @test all(a[i]^2 + d[i] >= abs(A[i, i]) - 1e-12 for i in axes(A, 1))
+            # d is non-negative
+            @test all(d[i] >= 0 for i in axes(A, 1))
+            # d is as tight as possible: d[i] == max(0, |A[i,i]| - a[i]^2)
+            @test all(d[i] ≈ max(0.0, abs(A[i, i]) - a[i]^2) for i in axes(A, 1))
+        end
+        # Non-square input is rejected
+        @test_throws ArgumentError symdiagcover([1.0 2.0; 3.0 4.0; 5.0 6.0])
+        # For a diagonal matrix, a should be all zeros and d should cover the diagonal
+        A_diag = [4.0 0.0; 0.0 9.0]
+        d, a = symdiagcover(A_diag)
+        @test all(iszero, a)
+        @test d ≈ [4.0, 9.0]
+        # symdiagcover gives a tighter diagonal cover than symcover when off-diagonal entries are tiny
+        A_tiny = [4.0 1e-8; 1e-8 1.0]
+        d2, a2 = symdiagcover(A_tiny)
+        a_sym = symcover(A_tiny)
+        # The Diagonal(d)+a*a' cover is valid
+        cover_mat = Diagonal(d2) + a2 * a2'
+        @test all(cover_mat[i, j] >= abs(A_tiny[i, j]) - 1e-12 for i in axes(A_tiny, 1), j in axes(A_tiny, 2))
+        # symdiagcover uses a smaller a[i] than symcover (the diagonal slack goes to d instead)
+        @test any(a_sym[i] > a2[i] + 1e-12 for i in axes(A_tiny, 1))
+    end
+
     @testset "cover_lobjective and cover_qobjective" begin
         A = [4.0 2.0; 2.0 1.0]
         a = symcover(A)