diff --git a/docs/src/index.md b/docs/src/index.md index 1d8ca9b..08147e9 100644 --- a/docs/src/index.md +++ b/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 diff --git a/src/ScaleInvariantAnalysis.jl b/src/ScaleInvariantAnalysis.jl index cb9d706..7579d19 100644 --- a/src/ScaleInvariantAnalysis.jl +++ b/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") diff --git a/src/covers.jl b/src/covers.jl index bf8dffc..c95ed6f 100644 --- a/src/covers.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 523391c..50384b6 100644 --- a/test/runtests.jl +++ b/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)