Reduce optimality gap for symmetric (#10)

It turns out that initializing with the solution to the unconstrained
problem produces better results (in the sense of a closer match to
quadratic optimality) than the previous diagonal initialization. Keep
both as options, but default to the higher-quality one. The reason to
keep the old strategy is that it is faster (about 2x on a 5x5 matrix),
largely by avoiding the need to take logarithms.

Author: timholy

docs/src/index.md

@@ -32,7 +32,7 @@ A concrete example: a 3×3 matrix whose rows and columns correspond to physical
 variables with different units (position in meters, velocity in m/s, force
 in N):
 
-```jldoctest coverones; filter = r"(\d+\.\d{6})\d+" => s"\1"
+```jldoctest coverones; filter = r"(\d+\.\d{2})\d+" => s"\1"
 julia> using ScaleInvariantAnalysis
 
 julia> A = [1e6 1e3 1.0; 1e3 1.0 1e-3; 1.0 1e-3 1e-6];

src/covers.jl

@@ -29,24 +29,29 @@ cover_qobjective(a, b, A) = sum(log(a[i] * b[j] / abs(A[i, j]))^2 for i in axes(
 cover_qobjective(a, A) = cover_qobjective(a, a, A)
 
 """
-    a = symcover(A; iter=3)
+    a = symcover(A; prioritize::Symbol=:quality, iter=3)
 
 Given a square matrix `A` assumed to be symmetric, return a vector `a`
 representing the symmetric cover of `A`, so that `a[i] * a[j] >= abs(A[i, j])`
 for all `i`, `j`.
 
-`a` may not be minimal, but it is tightened iteratively, with `iter` specifying
-the number of iterations (more iterations make tighter covers).
-`symcover` is fast and generally recommended for production use.
+`prioritize=:quality` yields a cover that is typically closer to being
+quadratically optimal, though there are exceptions.
+`prioritize=:speed` is often about twice as fast (with default `iter=3`). In
+either case, after initialization `a` is tightened iteratively, with `iter`
+specifying the number of iterations (more iterations make tighter covers).
+
+Regardless of which `prioritize` option is chosen, `symcover` is fast and
+generally recommended for production use.
 
 See also: [`symcover_lmin`](@ref), [`symcover_qmin`](@ref), [`cover`](@ref).
 
 # Examples
 
-```jldoctest
+```jldoctest; filter = r"(\\d+\\.\\d{6})\\d+" => s"\\1"
 julia> A = [4 -1; -1 0];
 
-julia> a = symcover(A)
+julia> a = symcover(A; prioritize=:speed)
 2-element Vector{Float64}:
  2.0
  0.5
@@ -55,18 +60,43 @@ julia> a * a'
 2×2 Matrix{Float64}:
  4.0  1.0
  1.0  0.25
+
+julia> A = [0 12 9; 12 7 12; 9 12 0];
+
+julia> a = symcover(A; prioritize=:quality)
+3-element Vector{Float64}:
+ 3.4021999694928753
+ 3.54528705924512
+ 3.3847752803845172
+
+julia> a * a'
+3×3 Matrix{Float64}:
+ 11.575   12.0618  11.5157
+ 12.0618  12.5691  12.0
+ 11.5157  12.0     11.4567
+
+julia> a = symcover(A; prioritize=:speed)
+3-element Vector{Float64}:
+ 4.535573676110727
+ 2.6457513110645907
+ 4.535573676110727
+
+julia> a * a'
+3×3 Matrix{Float64}:
+ 20.5714  12.0  20.5714
+ 12.0      7.0  12.0
+ 20.5714  12.0  20.5714
 ```
 """
-function symcover(A::AbstractMatrix; exclude_diagonal::Bool=false, kwargs...)
+function symcover(A::AbstractMatrix; exclude_diagonal::Bool=false, prioritize::Symbol=:quality, kwargs...)
+    prioritize in (:quality, :speed) || throw(ArgumentError("prioritize must be :quality or :speed"))
     ax = axes(A, 1)
     axes(A, 2) == ax || throw(ArgumentError("symcover requires a square matrix"))
     a = similar(A, float(eltype(A)), ax)
-    if exclude_diagonal
-        fill!(a, zero(eltype(a)))
+    if prioritize == :quality
+        _symcover_init_quadratic!(a, A; exclude_diagonal)
     else
-        for j in ax
-            a[j] = sqrt(abs(A[j, j]))
-        end
+        _symcover_init_fast!(a, A; exclude_diagonal)
     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.
@@ -94,6 +124,49 @@ function symcover(A::AbstractMatrix; exclude_diagonal::Bool=false, kwargs...)
     return tighten_cover!(a, A; exclude_diagonal, kwargs...)
 end
 
+function _symcover_init_quadratic!(a::AbstractVector{T}, A::AbstractMatrix; exclude_diagonal::Bool=false) where T
+    ax = eachindex(a)
+    loga = fill!(similar(a), zero(T))
+    nza  = fill(0, ax)
+    for j in ax
+        for i in first(ax):j - exclude_diagonal
+            Aij = abs(A[i, j])
+            iszero(Aij) && continue
+            lAij = log(Aij)
+            loga[i] += lAij
+            nza[i] += 1
+            if j != i
+                loga[j] += lAij
+                nza[j] += 1
+            end
+        end
+    end
+    nztotal = sum(nza)
+    halfmu = iszero(nztotal) ? zero(T) : sum(loga) / (2 * nztotal)
+    for i in ax
+        a[i] = ai = iszero(nza[i]) ? zero(T) : exp(loga[i] / nza[i] - halfmu)
+        if !exclude_diagonal
+            # The rest of the algorithm will ensure the initialization is a valid cover, but we have to do the diagonal here.
+            Aii = abs(A[i, i])
+            if ai^2 < Aii
+                a[i] = sqrt(Aii)
+            end
+        end
+    end
+    return a
+end
+
+function _symcover_init_fast!(a::AbstractVector{T}, A::AbstractMatrix; exclude_diagonal::Bool=false) where T
+    if exclude_diagonal
+        fill!(a, zero(T))
+    else
+        for j in eachindex(a)
+            a[j] = sqrt(abs(A[j, j]))
+        end
+    end
+    return a
+end
+
 """
     d, a = symdiagcover(A; kwargs...)
 

test/runtests.jl

@@ -10,8 +10,10 @@ using Test
         # Cover property: a[i]*a[j] >= abs(A[i,j]) for all i, j
         for A in ([2.0 1.0; 1.0 3.0], [1.0 -0.2; -0.2 0.0], [1.0 0.0; 0.0 0.0],
                   [100.0 1.0; 1.0 0.01])
-            a = symcover(A)
-            @test all(a[i] * a[j] >= abs(A[i, j]) - 1e-12 for i in axes(A, 1), j in axes(A, 2))
+            for prioritize in (:speed, :quality)
+                a = symcover(A; prioritize)
+                @test all(a[i] * a[j] >= abs(A[i, j]) - 1e-12 for i in axes(A, 1), j in axes(A, 2))
+            end
         end
         # Zero diagonal gives zero scale
         a = symcover([1.0 0; 0 0])
@@ -54,15 +56,17 @@ using Test
     @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))
+            for prioritize in (:speed, :quality)
+                d, a = symdiagcover(A; prioritize)
+                # 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
         end
         # Non-square input is rejected
         @test_throws ArgumentError symdiagcover([1.0 2.0; 3.0 4.0; 5.0 6.0])
@@ -161,35 +165,36 @@ using Test
             include("testmatrices.jl")   # defines symmetric_matrices and general_matrices
         end
 
-        # symcover cover_lobjective should be close to optimal (symcover_lmin)
+        # symcover cover_qobjective should be close to optimal (symcover_qmin)
         sym_ratios = Float64[]
         for (_, A) in symmetric_matrices
             Af = Float64.(A)
+            # Initialization should give a valid cover
             a0 = symcover(Af; iter=0)
             @test all(a0[i] * a0[j] >= abs(Af[i, j]) - 1e-12 for i in axes(Af, 1), j in axes(Af, 2))
             a0 = symcover(Af / 100; iter=0)
             @test all(a0[i] * a0[j] >= abs(Af[i, j])/100 - 1e-12 for i in axes(Af, 1), j in axes(Af, 2))
-            a = symcover(Af; iter=10)
-            lopt  = cover_lobjective(symcover_lmin(Af), Af)
-            lfast = cover_lobjective(a, Af)
-            iszero(lopt) || push!(sym_ratios, lfast / lopt)
+            # Covers are nearly quadratically optimal
+            qopt  = cover_qobjective(symcover_qmin(Af), Af)
+            qfast = cover_qobjective(symcover(Af; iter=10), Af)
+            iszero(qopt) || push!(sym_ratios, qfast / qopt)
         end
-        @test median(sym_ratios) < 1.1
+        @test median(sym_ratios) < 1.02
 
-        # cover cover_lobjective should be close to optimal (cover_lmin)
+        # cover cover_qobjective should be close to optimal (cover_qmin)
         gen_ratios = Float64[]
         for (_, A) in general_matrices
             Af = Float64.(A)
+            # Initialization should give a valid cover
             a0, b0  = cover(Af; iter=0)
             @test all(a0[i] * b0[j] >= abs(Af[i, j]) - 1e-12 for i in axes(Af, 1), j in axes(Af, 2))
             a0, b0  = cover(Af / 100; iter=0)
             @test all(a0[i] * b0[j] >= abs(Af[i, j])/100 - 1e-12 for i in axes(Af, 1), j in axes(Af, 2))
-            al, bl = cover_lmin(Af)
-            a,  b  = cover(Af; iter=10)
-            lopt  = cover_lobjective(al, bl, Af)
-            lfast = cover_lobjective(a,  b,  Af)
-            iszero(lopt) || push!(gen_ratios, lfast / lopt)
+            # Covers are nearly quadratically optimal
+            qopt  = cover_qobjective(cover_qmin(Af)..., Af)
+            qfast = cover_qobjective(cover(Af; iter=10)..., Af)
+            iszero(qopt) || push!(gen_ratios, qfast / qopt)
         end
-        @test median(gen_ratios) < 1.1
+        @test median(gen_ratios) < 1.02
     end
 end