Skip to content

Reduce optimality gap for symmetric #10

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
timholy merged 3 commits into from
Mar 10, 2026
Merged

Reduce optimality gap for symmetric #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b47893b
Reduce optimality gap for symmetric
timholy Mar 10, 2026
a80c50e
fix docs
timholy Mar 10, 2026
ff1e78b
Improve coverage
timholy Mar 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/src/index.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down
97 changes: 85 additions & 12 deletions src/covers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
Expand Down Expand Up @@ -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...)

Expand Down
53 changes: 29 additions & 24 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Expand Down Expand Up @@ -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])
Expand Down Expand Up @@ -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
Loading