diff --git a/Project.toml b/Project.toml index 2464435..08a2ea8 100644 --- a/Project.toml +++ b/Project.toml @@ -6,14 +6,15 @@ authors = ["Tim Holy and contributors"] [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [weakdeps] HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [extensions] SIAJuMP = ["HiGHS", "JuMP"] +SIASparseArrays = "SparseArrays" [compat] Graphs = "1" @@ -36,8 +37,9 @@ NautyGraphs = "7509a0a4-015a-4167-b44b-0799a1a2605e" ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Graphs", "HiGHS", "JuMP", "NautyGraphs", "ParametricOptInterface", "ProgressMeter", "Random", "Statistics", "Test"] +test = ["Graphs", "HiGHS", "JuMP", "NautyGraphs", "ParametricOptInterface", "ProgressMeter", "Random", "SparseArrays", "Statistics", "Test"] diff --git a/ext/SIASparseArrays.jl b/ext/SIASparseArrays.jl new file mode 100644 index 0000000..07fcea5 --- /dev/null +++ b/ext/SIASparseArrays.jl @@ -0,0 +1,430 @@ +module SIASparseArrays + +using LinearAlgebra +using SparseArrays +using ScaleInvariantAnalysis + +# ============================================================ +# Private helpers (operate on a plain SparseMatrixCSC parent) +# ============================================================ + +# Quadratic initialisation of the symmetric cover. +# `uplo` controls which triangle is treated as canonical ('U' → i ≤ j, 'L' → i ≥ j). +# For a fully-stored SparseMatrixCSC, use 'U' to avoid double-counting. +function _sparse_symcover_init_quadratic!(a::AbstractVector{T}, P::SparseMatrixCSC, uplo::Char; exclude_diagonal::Bool=false) where T + rv = rowvals(P) + nzs = nonzeros(P) + loga = zeros(T, size(P, 1)) + nza = zeros(Int, size(P, 1)) + for j in axes(P, 2) + for k in nzrange(P, j) + i = rv[k] + (uplo == 'U' ? i > j : i < j) && continue # canonical triangle only + exclude_diagonal && i == j && continue + Aij = abs(nzs[k]) + iszero(Aij) && continue + lAij = log(Aij) + loga[i] += lAij + nza[i] += 1 + if i != j + loga[j] += lAij + nza[j] += 1 + end + end + end + nztotal = sum(nza) + halfmu = iszero(nztotal) ? zero(T) : sum(loga) / (2 * nztotal) + for i in eachindex(a) + a[i] = iszero(nza[i]) ? zero(T) : exp(loga[i] / nza[i] - halfmu) + end + if !exclude_diagonal + # Enforce a[i]^2 >= |A[i,i]| for stored diagonal entries + for j in axes(P, 2) + for k in nzrange(P, j) + i = rv[k] + i != j && continue + Aii = abs(nzs[k]) + if a[i]^2 < Aii + a[i] = sqrt(Aii) + end + break + end + end + end + return a +end + +# Fast (diagonal-based) initialisation of the symmetric cover. +function _sparse_symcover_init_fast!(a::AbstractVector{T}, P::SparseMatrixCSC; exclude_diagonal::Bool=false) where T + fill!(a, zero(T)) + exclude_diagonal && return a + rv = rowvals(P) + nzs = nonzeros(P) + for j in axes(P, 2) + for k in nzrange(P, j) + i = rv[k] + if i == j + a[j] = sqrt(abs(nzs[k])) + break + end + end + end + return a +end + +# Boost: for each stored off-diagonal entry in the canonical triangle, ensure a[i]*a[j] >= |A[i,j]|. +function _sparse_symcover_boost!(a::AbstractVector, P::SparseMatrixCSC, uplo::Char) + rv = rowvals(P) + nzs = nonzeros(P) + for j in axes(P, 2) + for k in nzrange(P, j) + i = rv[k] + i == j && continue # skip diagonal + (uplo == 'U' ? i > j : i < j) && continue # canonical off-diagonal only + Aij = abs(nzs[k]) + ai, aj = a[i], a[j] + if iszero(aj) + if !iszero(ai) + a[j] = Aij / ai + else + a[i] = a[j] = sqrt(Aij) + end + elseif iszero(ai) + a[i] = Aij / aj + else + aprod = ai * aj + aprod >= Aij && continue + s = sqrt(Aij / aprod) + a[i] = s * ai + a[j] = s * aj + end + end + end +end + +# Tighten the symmetric cover, iterating only stored nonzeros. +# Works for fully-stored symmetric matrices and for one-triangle-stored Symmetric wrappers: +# for each column j, updating both aratio[i] (directly) and aratio[j] (via aratioj) correctly +# propagates the constraint a[i]*a[j] >= |A[i,j]| to both sides. +function _tighten_cover_sym_sparse!(a::AbstractVector{T}, P::SparseMatrixCSC; iter::Int, exclude_diagonal::Bool) where T + rv = rowvals(P) + nzs = nonzeros(P) + aratio = similar(a) + for _ in 1:iter + fill!(aratio, typemax(T)) + for j in axes(P, 2) + aratioj = aratio[j] + aj = a[j] + for k in nzrange(P, j) + i = rv[k] + (exclude_diagonal && i == j) && continue + Aij = T(abs(nzs[k])) + 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 + a ./= sqrt.(aratio) + end + return a +end + +# Tighten the asymmetric cover, iterating only stored nonzeros. +function _tighten_cover_asym_sparse!(a::AbstractVector{T}, b::AbstractVector{T}, P::SparseMatrixCSC; iter::Int) where T + rv = rowvals(P) + nzs = nonzeros(P) + aratio = fill(typemax(T), eachindex(a)) + bratio = fill(typemax(T), eachindex(b)) + for _ in 1:iter + fill!(aratio, typemax(T)) + fill!(bratio, typemax(T)) + for j in eachindex(b) + bratioj = bratio[j] + bj = b[j] + for k in nzrange(P, j) + i = rv[k] + Aij = T(abs(nzs[k])) + r = ifelse(iszero(Aij), typemax(T), a[i] * bj / Aij) + aratio[i] = min(aratio[i], r) + bratioj = min(bratioj, r) + end + bratio[j] = bratioj + end + a ./= sqrt.(aratio) + b ./= sqrt.(bratio) + end + return a, b +end + +# Accumulate objective sums over the canonical triangle of a symmetric matrix. +# For each off-diagonal stored entry (i,j), both (i,j) and (j,i) contribute to the sum. +function _cover_lobjective_sym_sparse(a, b, P::SparseMatrixCSC, uplo::Char) + rv = rowvals(P) + nzs = nonzeros(P) + s = zero(float(promote_type(eltype(a), eltype(b), eltype(P)))) + for j in axes(P, 2) + bj = b[j] + for k in nzrange(P, j) + i = rv[k] + (uplo == 'U' ? i > j : i < j) && continue + Aij = abs(nzs[k]) + iszero(Aij) && continue + s += log(a[i] * bj / Aij) + if i != j + s += log(a[j] * b[i] / Aij) + end + end + end + return s +end + +function _cover_qobjective_sym_sparse(a, b, P::SparseMatrixCSC, uplo::Char) + rv = rowvals(P) + nzs = nonzeros(P) + s = zero(float(promote_type(eltype(a), eltype(b), eltype(P)))) + for j in axes(P, 2) + bj = b[j] + for k in nzrange(P, j) + i = rv[k] + (uplo == 'U' ? i > j : i < j) && continue + Aij = abs(nzs[k]) + iszero(Aij) && continue + s += log(a[i] * bj / Aij)^2 + if i != j + s += log(a[j] * b[i] / Aij)^2 + end + end + end + return s +end + +# Extract the diagonal correction vector for symdiagcover. +function _symdiagcover_d!(d::AbstractVector{T}, a::AbstractVector, P::SparseMatrixCSC) where T + rv = rowvals(P) + nzs = nonzeros(P) + for j in axes(P, 2) + for k in nzrange(P, j) + i = rv[k] + if i == j + Aii = abs(nzs[k]) + d[i] = max(zero(T), Aii - a[i]^2) + break + end + end + end +end + +# ============================================================ +# SparseMatrixCSC methods +# ============================================================ + +function ScaleInvariantAnalysis.cover_lobjective(a, b, A::SparseMatrixCSC) + rv = rowvals(A) + nzs = nonzeros(A) + s = zero(float(promote_type(eltype(a), eltype(b), eltype(A)))) + for j in axes(A, 2) + bj = b[j] + for k in nzrange(A, j) + Aij = abs(nzs[k]) + iszero(Aij) && continue + s += log(a[rv[k]] * bj / Aij) + end + end + return s +end + +function ScaleInvariantAnalysis.cover_qobjective(a, b, A::SparseMatrixCSC) + rv = rowvals(A) + nzs = nonzeros(A) + s = zero(float(promote_type(eltype(a), eltype(b), eltype(A)))) + for j in axes(A, 2) + bj = b[j] + for k in nzrange(A, j) + Aij = abs(nzs[k]) + iszero(Aij) && continue + s += log(a[rv[k]] * bj / Aij)^2 + end + end + return s +end + +function ScaleInvariantAnalysis.tighten_cover!(a::AbstractVector{T}, A::SparseMatrixCSC; 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`")) + return _tighten_cover_sym_sparse!(a, A; iter, exclude_diagonal) +end + +function ScaleInvariantAnalysis.tighten_cover!(a::AbstractVector{T}, b::AbstractVector{T}, A::SparseMatrixCSC; iter::Int=3) where T + eachindex(a) == axes(A, 1) || throw(DimensionMismatch("indices of a must match row-indexing of A")) + eachindex(b) == axes(A, 2) || throw(DimensionMismatch("indices of b must match column-indexing of A")) + return _tighten_cover_asym_sparse!(a, b, A; iter) +end + +function ScaleInvariantAnalysis.symcover(A::SparseMatrixCSC; 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")) + T = float(eltype(A)) + a = zeros(T, size(A, 1)) + if prioritize == :quality + _sparse_symcover_init_quadratic!(a, A, 'U'; exclude_diagonal) + else + _sparse_symcover_init_fast!(a, A; exclude_diagonal) + end + _sparse_symcover_boost!(a, A, 'U') + return _tighten_cover_sym_sparse!(a, A; iter=get(kwargs, :iter, 3), exclude_diagonal) +end + +function ScaleInvariantAnalysis.symdiagcover(A::SparseMatrixCSC; kwargs...) + axes(A, 1) == axes(A, 2) || throw(ArgumentError("symcover requires a square matrix")) + a = ScaleInvariantAnalysis.symcover(A; exclude_diagonal=true, kwargs...) + T = float(eltype(A)) + d = zeros(T, size(A, 1)) + _symdiagcover_d!(d, a, A) + return d, a +end + +function ScaleInvariantAnalysis.cover(A::SparseMatrixCSC; kwargs...) + T = float(eltype(A)) + a = zeros(T, size(A, 1)) + b = zeros(T, size(A, 2)) + loga = zeros(T, size(A, 1)) + logb = zeros(T, size(A, 2)) + nza = zeros(Int, size(A, 1)) + nzb = zeros(Int, size(A, 2)) + logmu = zero(T) + nztotal = 0 + rv = rowvals(A) + nzs = nonzeros(A) + for j in axes(A, 2) + for k in nzrange(A, j) + Aij = abs(nzs[k]) + iszero(Aij) && continue + i = rv[k] + lAij = log(Aij) + loga[i] += lAij + logb[j] += lAij + nza[i] += 1 + nzb[j] += 1 + logmu += lAij + nztotal += 1 + end + end + halfmu = iszero(nztotal) ? zero(T) : T(logmu / (2 * nztotal)) + for i in axes(A, 1) + a[i] = iszero(nza[i]) ? zero(T) : exp(loga[i] / nza[i] - halfmu) + end + for j in axes(A, 2) + b[j] = iszero(nzb[j]) ? zero(T) : exp(logb[j] / nzb[j] - halfmu) + end + for j in axes(A, 2) + bj = b[j] + for k in nzrange(A, j) + i = rv[k] + Aij, ai = abs(nzs[k]), a[i] + aprod = ai * bj + aprod >= Aij && continue + s = sqrt(Aij / aprod) + a[i] = s * ai + b[j] = bj = s * bj + end + end + return _tighten_cover_asym_sparse!(a, b, A; iter=get(kwargs, :iter, 3)) +end + +# ============================================================ +# Symmetric{<:Any, <:SparseMatrixCSC} methods +# ============================================================ + +function ScaleInvariantAnalysis.cover_lobjective(a, b, S::Symmetric{<:Any, <:SparseMatrixCSC}) + return _cover_lobjective_sym_sparse(a, b, parent(S), S.uplo) +end + +function ScaleInvariantAnalysis.cover_qobjective(a, b, S::Symmetric{<:Any, <:SparseMatrixCSC}) + return _cover_qobjective_sym_sparse(a, b, parent(S), S.uplo) +end + +function ScaleInvariantAnalysis.tighten_cover!(a::AbstractVector{T}, S::Symmetric{<:Any, <:SparseMatrixCSC}; iter::Int=3, exclude_diagonal::Bool=false) where T + P = parent(S) + ax = axes(P, 1) + axes(P, 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`")) + return _tighten_cover_sym_sparse!(a, P; iter, exclude_diagonal) +end + +function ScaleInvariantAnalysis.symcover(S::Symmetric{<:Any, <:SparseMatrixCSC}; exclude_diagonal::Bool=false, prioritize::Symbol=:quality, kwargs...) + prioritize in (:quality, :speed) || throw(ArgumentError("prioritize must be :quality or :speed")) + P = parent(S) + axes(P, 1) == axes(P, 2) || throw(ArgumentError("symcover requires a square matrix")) + T = float(eltype(P)) + a = zeros(T, size(P, 1)) + uplo = S.uplo + if prioritize == :quality + _sparse_symcover_init_quadratic!(a, P, uplo; exclude_diagonal) + else + _sparse_symcover_init_fast!(a, P; exclude_diagonal) + end + _sparse_symcover_boost!(a, P, uplo) + return _tighten_cover_sym_sparse!(a, P; iter=get(kwargs, :iter, 3), exclude_diagonal) +end + +function ScaleInvariantAnalysis.symdiagcover(S::Symmetric{<:Any, <:SparseMatrixCSC}; kwargs...) + P = parent(S) + axes(P, 1) == axes(P, 2) || throw(ArgumentError("symcover requires a square matrix")) + a = ScaleInvariantAnalysis.symcover(S; exclude_diagonal=true, kwargs...) + T = float(eltype(P)) + d = zeros(T, size(P, 1)) + _symdiagcover_d!(d, a, P) + return d, a +end + +# ============================================================ +# Hermitian{<:Real, <:SparseMatrixCSC} methods +# (Real-valued Hermitian is equivalent to Symmetric) +# ============================================================ + +function ScaleInvariantAnalysis.cover_lobjective(a, b, H::Hermitian{<:Real, <:SparseMatrixCSC}) + return _cover_lobjective_sym_sparse(a, b, parent(H), H.uplo) +end + +function ScaleInvariantAnalysis.cover_qobjective(a, b, H::Hermitian{<:Real, <:SparseMatrixCSC}) + return _cover_qobjective_sym_sparse(a, b, parent(H), H.uplo) +end + +function ScaleInvariantAnalysis.tighten_cover!(a::AbstractVector{T}, H::Hermitian{<:Real, <:SparseMatrixCSC}; iter::Int=3, exclude_diagonal::Bool=false) where T + P = parent(H) + ax = axes(P, 1) + axes(P, 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`")) + return _tighten_cover_sym_sparse!(a, P; iter, exclude_diagonal) +end + +function ScaleInvariantAnalysis.symcover(H::Hermitian{<:Real, <:SparseMatrixCSC}; exclude_diagonal::Bool=false, prioritize::Symbol=:quality, kwargs...) + prioritize in (:quality, :speed) || throw(ArgumentError("prioritize must be :quality or :speed")) + P = parent(H) + axes(P, 1) == axes(P, 2) || throw(ArgumentError("symcover requires a square matrix")) + T = float(eltype(P)) + a = zeros(T, size(P, 1)) + uplo = H.uplo + if prioritize == :quality + _sparse_symcover_init_quadratic!(a, P, uplo; exclude_diagonal) + else + _sparse_symcover_init_fast!(a, P; exclude_diagonal) + end + _sparse_symcover_boost!(a, P, uplo) + return _tighten_cover_sym_sparse!(a, P; iter=get(kwargs, :iter, 3), exclude_diagonal) +end + +function ScaleInvariantAnalysis.symdiagcover(H::Hermitian{<:Real, <:SparseMatrixCSC}; kwargs...) + P = parent(H) + axes(P, 1) == axes(P, 2) || throw(ArgumentError("symcover requires a square matrix")) + a = ScaleInvariantAnalysis.symcover(H; exclude_diagonal=true, kwargs...) + T = float(eltype(P)) + d = zeros(T, size(P, 1)) + _symdiagcover_d!(d, a, P) + return d, a +end + +end # module SIASparseArrays diff --git a/src/ScaleInvariantAnalysis.jl b/src/ScaleInvariantAnalysis.jl index b64c038..837b395 100644 --- a/src/ScaleInvariantAnalysis.jl +++ b/src/ScaleInvariantAnalysis.jl @@ -1,13 +1,13 @@ module ScaleInvariantAnalysis using LinearAlgebra -using SparseArrays using PrecompileTools export cover_lobjective, cover_qobjective, cover, symcover, symdiagcover, cover_lmin, symcover_lmin, cover_qmin, symcover_qmin export dotabs include("covers.jl") +include("structured.jl") """ diff --git a/src/covers.jl b/src/covers.jl index 59b7bad..d3120be 100644 --- a/src/covers.jl +++ b/src/covers.jl @@ -404,3 +404,35 @@ function cover_qmin end Similar to [`cover_qmin`](@ref), but returns a linear-minimal cover of `A`. """ function cover_lmin end + +# ============================================================ +# Adjoint and Transpose wrappers +# +# If B = adjoint(P) or transpose(P), then |B[i,j]| = |P[j,i]|, so a cover +# (a, b) of B satisfies a[i]*b[j] >= |P[j,i]|, which is exactly a cover +# (b, a) of P. Therefore: unwrap the parent, compute its cover, and swap. +# ============================================================ + +cover_lobjective(a, b, A::Adjoint) = cover_lobjective(b, a, parent(A)) +cover_lobjective(a, b, A::Transpose) = cover_lobjective(b, a, parent(A)) + +cover_qobjective(a, b, A::Adjoint) = cover_qobjective(b, a, parent(A)) +cover_qobjective(a, b, A::Transpose) = cover_qobjective(b, a, parent(A)) + +function tighten_cover!(a::AbstractVector{T}, b::AbstractVector{T}, A::Adjoint; kwargs...) where T + tighten_cover!(b, a, parent(A); kwargs...) + return a, b +end +function tighten_cover!(a::AbstractVector{T}, b::AbstractVector{T}, A::Transpose; kwargs...) where T + tighten_cover!(b, a, parent(A); kwargs...) + return a, b +end + +function cover(A::Adjoint; kwargs...) + a, b = cover(parent(A); kwargs...) + return b, a +end +function cover(A::Transpose; kwargs...) + a, b = cover(parent(A); kwargs...) + return b, a +end diff --git a/src/structured.jl b/src/structured.jl new file mode 100644 index 0000000..ba629e9 --- /dev/null +++ b/src/structured.jl @@ -0,0 +1,330 @@ +# Cover methods for structured LinearAlgebra matrix types. +# +# Diagonal is handled individually. +# SymTridiagonal, Bidiagonal, and Tridiagonal share a single implementation +# via the PlusMinus1Banded union: all are square with entries only on the main +# diagonal and the ±1 off-diagonals. Structural zeros (e.g. the sub-diagonal +# of an upper Bidiagonal) are returned as exact zero by getindex, so they are +# skipped cheaply by the iszero checks. + +const PlusMinus1Banded = Union{SymTridiagonal, Bidiagonal, Tridiagonal} + +# ============================================================ +# Diagonal +# ============================================================ + +function cover_lobjective(a, b, D::Diagonal) + T = float(promote_type(eltype(a), eltype(b), eltype(D))) + s = zero(T) + for i in eachindex(D.diag) + Dii = abs(D.diag[i]) + iszero(Dii) && continue + s += log(T(a[i]) * T(b[i]) / T(Dii)) + end + return s +end + +function cover_qobjective(a, b, D::Diagonal) + T = float(promote_type(eltype(a), eltype(b), eltype(D))) + s = zero(T) + for i in eachindex(D.diag) + Dii = abs(D.diag[i]) + iszero(Dii) && continue + s += log(T(a[i]) * T(b[i]) / T(Dii))^2 + end + return s +end + +function tighten_cover!(a::AbstractVector{T}, D::Diagonal; iter::Int=3, exclude_diagonal::Bool=false) where T + exclude_diagonal && return a + for i in eachindex(a, D.diag) + Dii = T(abs(D.diag[i])) + iszero(Dii) || (a[i] = sqrt(Dii)) + end + return a +end + +function tighten_cover!(a::AbstractVector{T}, b::AbstractVector{T}, D::Diagonal; iter::Int=3) where T + for i in eachindex(a, b, D.diag) + Dii = T(abs(D.diag[i])) + iszero(Dii) && continue + aprod = a[i] * b[i] + iszero(aprod) && continue + s = sqrt(aprod / Dii) + a[i] /= s + b[i] /= s + end + return a, b +end + +function symcover(D::Diagonal; exclude_diagonal::Bool=false, kwargs...) + T = float(eltype(D)) + a = zeros(T, length(D.diag)) + exclude_diagonal && return a + for i in eachindex(a, D.diag) + a[i] = sqrt(T(abs(D.diag[i]))) + end + return a +end + +function symdiagcover(D::Diagonal; kwargs...) + T = float(eltype(D)) + return T.(abs.(D.diag)), zeros(T, length(D.diag)) +end + +function cover(D::Diagonal; kwargs...) + a = symcover(D) + return tighten_cover!(a, copy(a), D; kwargs...) +end + +# ============================================================ +# PlusMinus1Banded (SymTridiagonal, Bidiagonal, Tridiagonal) +# ============================================================ + +function cover_lobjective(a, b, A::PlusMinus1Banded) + T = float(promote_type(eltype(a), eltype(b), eltype(A))) + s = zero(T) + n = size(A, 1) + for i in 1:n + Aii = abs(A[i, i]) + iszero(Aii) || (s += log(T(a[i]) * T(b[i]) / T(Aii))) + end + for i in 1:n-1 + Aij = abs(A[i, i+1]) + iszero(Aij) || (s += log(T(a[i]) * T(b[i+1]) / T(Aij))) + Aij = abs(A[i+1, i]) + iszero(Aij) || (s += log(T(a[i+1]) * T(b[i]) / T(Aij))) + end + return s +end + +function cover_qobjective(a, b, A::PlusMinus1Banded) + T = float(promote_type(eltype(a), eltype(b), eltype(A))) + s = zero(T) + n = size(A, 1) + for i in 1:n + Aii = abs(A[i, i]) + iszero(Aii) || (s += log(T(a[i]) * T(b[i]) / T(Aii))^2) + end + for i in 1:n-1 + Aij = abs(A[i, i+1]) + iszero(Aij) || (s += log(T(a[i]) * T(b[i+1]) / T(Aij))^2) + Aij = abs(A[i+1, i]) + iszero(Aij) || (s += log(T(a[i+1]) * T(b[i]) / T(Aij))^2) + end + return s +end + +# Symmetric tighten: both A[i,i+1] and A[i+1,i] are checked independently. +# For SymTridiagonal they are equal (no-op redundancy); for a Tridiagonal that +# the caller asserts is symmetric, using both is consistent with the general +# tighten_cover!(a, A::AbstractMatrix) which iterates over every entry. +function tighten_cover!(a::AbstractVector{T}, A::PlusMinus1Banded; iter::Int=3, exclude_diagonal::Bool=false) where T + n = size(A, 1) + aratio = similar(a) + for _ in 1:iter + fill!(aratio, typemax(T)) + if !exclude_diagonal + for i in 1:n + Aii = T(abs(A[i, i])) + iszero(Aii) && continue + aratio[i] = min(aratio[i], a[i]^2 / Aii) + end + end + for i in 1:n-1 + Aij = T(abs(A[i, i+1])) + if !iszero(Aij) + r = a[i] * a[i+1] / Aij + aratio[i] = min(aratio[i], r) + aratio[i+1] = min(aratio[i+1], r) + end + Aij = T(abs(A[i+1, i])) + if !iszero(Aij) + r = a[i] * a[i+1] / Aij + aratio[i] = min(aratio[i], r) + aratio[i+1] = min(aratio[i+1], r) + end + end + a ./= sqrt.(aratio) + end + return a +end + +# Asymmetric tighten for all ±1-banded types. +function tighten_cover!(a::AbstractVector{T}, b::AbstractVector{T}, A::PlusMinus1Banded; iter::Int=3) where T + n = size(A, 1) + aratio = similar(a) + bratio = similar(b) + for _ in 1:iter + fill!(aratio, typemax(T)) + fill!(bratio, typemax(T)) + for i in 1:n + Aii = T(abs(A[i, i])) + iszero(Aii) && continue + r = a[i] * b[i] / Aii + aratio[i] = min(aratio[i], r) + bratio[i] = min(bratio[i], r) + end + for i in 1:n-1 + Aij = T(abs(A[i, i+1])) + if !iszero(Aij) + r = a[i] * b[i+1] / Aij + aratio[i] = min(aratio[i], r) + bratio[i+1] = min(bratio[i+1], r) + end + Aij = T(abs(A[i+1, i])) + if !iszero(Aij) + r = a[i+1] * b[i] / Aij + aratio[i+1] = min(aratio[i+1], r) + bratio[i] = min(bratio[i], r) + end + end + a ./= sqrt.(aratio) + b ./= sqrt.(bratio) + end + return a, b +end + +# symcover and symdiagcover apply to any PlusMinus1Banded matrix; the caller +# asserts that A is symmetric in value (or accepts that only the upper triangle +# is used for initialization). +function symcover(A::PlusMinus1Banded; exclude_diagonal::Bool=false, prioritize::Symbol=:quality, kwargs...) + prioritize in (:quality, :speed) || throw(ArgumentError("prioritize must be :quality or :speed")) + n = size(A, 1) + T = float(eltype(A)) + a = zeros(T, n) + if prioritize == :quality + loga = zeros(T, n) + nza = zeros(Int, n) + if !exclude_diagonal + for i in 1:n + Aii = abs(A[i, i]) + iszero(Aii) && continue + loga[i] += log(Aii) + nza[i] += 1 + end + end + for i in 1:n-1 + Aij = abs(A[i, i+1]) # upper triangle; caller asserts A[i+1,i] matches + iszero(Aij) && continue + lAij = log(Aij) + loga[i] += lAij; nza[i] += 1 + loga[i+1] += lAij; nza[i+1] += 1 + end + nztotal = sum(nza) + halfmu = iszero(nztotal) ? zero(T) : sum(loga) / (2 * nztotal) + for i in 1:n + a[i] = iszero(nza[i]) ? zero(T) : exp(loga[i] / nza[i] - halfmu) + end + if !exclude_diagonal + for i in 1:n + Aii = T(abs(A[i, i])) + a[i]^2 < Aii && (a[i] = sqrt(Aii)) + end + end + else # :speed + if !exclude_diagonal + for i in 1:n + a[i] = sqrt(T(abs(A[i, i]))) + end + end + end + # Boost from off-diagonal entries (upper triangle) + for i in 1:n-1 + Aij = T(abs(A[i, i+1])) + iszero(Aij) && continue + ai, aj = a[i], a[i+1] + if iszero(aj) + iszero(ai) ? (a[i] = a[i+1] = sqrt(Aij)) : (a[i+1] = Aij / ai) + elseif iszero(ai) + a[i] = Aij / aj + else + aprod = ai * aj + if aprod < Aij + s = sqrt(Aij / aprod) + a[i] *= s; a[i+1] *= s + end + end + end + return tighten_cover!(a, A; iter=get(kwargs, :iter, 3), exclude_diagonal) +end + +function symdiagcover(A::PlusMinus1Banded; kwargs...) + a = symcover(A; exclude_diagonal=true, kwargs...) + T = float(eltype(A)) + d = map(1:size(A, 1)) do i + Aii = T(abs(A[i, i])) + max(zero(T), Aii - a[i]^2) + end + return d, a +end + +function cover(A::PlusMinus1Banded; kwargs...) + T = float(eltype(A)) + n = size(A, 1) + a = zeros(T, n) + b = zeros(T, n) + loga = zeros(T, n) + logb = zeros(T, n) + nza = zeros(Int, n) + nzb = zeros(Int, n) + logmu = zero(T) + nztotal = 0 + for i in 1:n + Aii = abs(A[i, i]) + iszero(Aii) && continue + lAii = log(T(Aii)) + loga[i] += lAii; logb[i] += lAii + nza[i] += 1; nzb[i] += 1 + logmu += lAii; nztotal += 1 + end + for i in 1:n-1 + Aij = abs(A[i, i+1]) + if !iszero(Aij) + lAij = log(T(Aij)) + loga[i] += lAij; logb[i+1] += lAij + nza[i] += 1; nzb[i+1] += 1 + logmu += lAij; nztotal += 1 + end + Aij = abs(A[i+1, i]) + if !iszero(Aij) + lAij = log(T(Aij)) + loga[i+1] += lAij; logb[i] += lAij + nza[i+1] += 1; nzb[i] += 1 + logmu += lAij; nztotal += 1 + end + end + halfmu = iszero(nztotal) ? zero(T) : T(logmu / (2 * nztotal)) + for i in 1:n + a[i] = iszero(nza[i]) ? zero(T) : exp(loga[i] / nza[i] - halfmu) + b[i] = iszero(nzb[i]) ? zero(T) : exp(logb[i] / nzb[i] - halfmu) + end + # Boost diagonal + for i in 1:n + Aii = T(abs(A[i, i])) + aprod = a[i] * b[i] + aprod >= Aii && continue + s = sqrt(Aii / aprod) + a[i] *= s; b[i] *= s + end + # Boost off-diagonal + for i in 1:n-1 + Aij = T(abs(A[i, i+1])) + if !iszero(Aij) + aprod = a[i] * b[i+1] + if aprod < Aij + s = sqrt(Aij / aprod) + a[i] *= s; b[i+1] *= s + end + end + Aij = T(abs(A[i+1, i])) + if !iszero(Aij) + aprod = a[i+1] * b[i] + if aprod < Aij + s = sqrt(Aij / aprod) + a[i+1] *= s; b[i] *= s + end + end + end + return tighten_cover!(a, b, A; kwargs...) +end diff --git a/test/runtests.jl b/test/runtests.jl index f5abd59..96a6057 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ using ScaleInvariantAnalysis using ScaleInvariantAnalysis: divmag, dotabs using JuMP, HiGHS # triggers the SIAJuMP extension (symcover_lmin, cover_lmin, etc.) +using SparseArrays # triggers the SIASparseArrays extension using LinearAlgebra using Statistics: median using Test @@ -160,6 +161,124 @@ using Test end end + @testset "SparseMatrixCSC" begin + for Adense in ([2.0 1.0; 1.0 3.0], [1.0 -0.2; -0.2 0.0], [0.0 12.0 9.0; 12.0 7.0 12.0; 9.0 12.0 0.0], + [100.0 1.0; 1.0 0.01]) + for A in (sparse(Adense), Symmetric(sparse(tril(Adense)), :L), Symmetric(sparse(triu(Adense)), :U), + Hermitian(sparse(tril(Adense)), :L), Hermitian(sparse(triu(Adense)), :U)) + a = symcover(A) + # Cover property + @test all(a[i] * a[j] >= abs(Adense[i, j]) - 1e-12 for i in axes(Adense, 1), j in axes(Adense, 2)) + # Objectives match dense + @test cover_lobjective(a, A) ≈ cover_lobjective(a, Adense) + @test cover_qobjective(a, A) ≈ cover_qobjective(a, Adense) + end + end + for Adense in ([2.0 1.0; 1.0 3.0], [0.0 1.0; -2.0 0.0], [1.0 2.0 3.0; 4.0 5.0 6.0]) + A = sparse(Adense) + a, b = cover(A) + @test all(a[i] * b[j] >= abs(Adense[i, j]) - 1e-12 for i in axes(Adense, 1), j in axes(Adense, 2)) + @test cover_lobjective(a, b, A) ≈ cover_lobjective(a, b, Adense) + @test cover_qobjective(a, b, A) ≈ cover_qobjective(a, b, Adense) + end + for Adense in ([2.0 1.0; 1.0 3.0], [4.0 1e-8; 1e-8 1.0], [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0]) + for A in (sparse(Adense), Symmetric(sparse(tril(Adense)), :L), Symmetric(sparse(triu(Adense)), :U), + Hermitian(sparse(tril(Adense)), :L), Hermitian(sparse(triu(Adense)), :U)) + d, a = symdiagcover(A) + # Off-diagonal cover property + @test all(a[i] * a[j] >= abs(Adense[i, j]) - 1e-12 for i in axes(Adense, 1), j in axes(Adense, 2) if i != j) + # Diagonal cover property + @test all(a[i]^2 + d[i] >= abs(Adense[i, i]) - 1e-12 for i in axes(Adense, 1)) + # d is non-negative + @test all(d[i] >= 0 for i in axes(Adense, 1)) + end + end + # Zero-diagonal sparse matrix + A0 = sparse([0.0 1.0; 1.0 0.0]) + @test symcover(A0) == [1.0, 1.0] + # Diagonal sparse matrix: a should be zero, d covers diagonal + Adiag = sparse(Diagonal([4.0, 9.0])) + d, a = symdiagcover(Adiag) + @test all(iszero, a) + @test d ≈ [4.0, 9.0] + end + + @testset "structured matrices" begin + @testset "Diagonal" begin + for dv in ([4.0, 9.0, 1.0], [4.0, 0.0, 1.0], [0.25, 100.0]) + D = Diagonal(dv) + Ddense = Matrix(D) + a = symcover(D) + @test all(a[i] * a[j] >= abs(Ddense[i, j]) - 1e-12 for i in axes(Ddense, 1), j in axes(Ddense, 2)) + @test cover_lobjective(a, D) ≈ cover_lobjective(a, Ddense) + @test cover_qobjective(a, D) ≈ cover_qobjective(a, Ddense) + d, a2 = symdiagcover(D) + @test all(iszero, a2) + @test d ≈ abs.(dv) + a3, b3 = cover(D) + @test all(a3[i] * b3[j] >= abs(Ddense[i, j]) - 1e-12 for i in axes(Ddense, 1), j in axes(Ddense, 2)) + @test cover_lobjective(a3, b3, D) ≈ cover_lobjective(a3, b3, Ddense) + @test cover_qobjective(a3, b3, D) ≈ cover_qobjective(a3, b3, Ddense) + end + @test all(iszero, symcover(Diagonal([1.0, 2.0]); exclude_diagonal=true)) + end + + @testset "PlusMinus1Banded" begin + asym_cases = [ + Bidiagonal([3.0, 2.0, 1.0], [6.0, 0.5], :U), + Bidiagonal([3.0, 2.0, 1.0], [6.0, 0.5], :L), + Tridiagonal([1.0, 0.5], [3.0, 2.0, 1.0], [4.0, 0.5]), + ] + for A in asym_cases + Adense = Matrix(A) + n = size(A, 1) + a, b = cover(A) + @test all(a[i] * b[j] >= abs(Adense[i, j]) - 1e-12 for i in 1:n, j in 1:n) + @test cover_lobjective(a, b, A) ≈ cover_lobjective(a, b, Adense) + @test cover_qobjective(a, b, A) ≈ cover_qobjective(a, b, Adense) + end + sym_cases = [ + SymTridiagonal([4.0, 3.0, 1.0], [2.0, 0.5]), + SymTridiagonal([0.0, 3.0, 0.0], [2.0, 0.5]), + # Tridiagonal that the caller asserts is symmetric + Tridiagonal([2.0, 0.5], [4.0, 3.0, 1.0], [2.0, 0.5]), + ] + for A in sym_cases + Adense = Matrix(A) + n = size(A, 1) + a, b = cover(A) + @test all(a[i] * b[j] >= abs(Adense[i, j]) - 1e-12 for i in 1:n, j in 1:n) + for prioritize in (:speed, :quality) + a = symcover(A; prioritize) + @test all(a[i] * a[j] >= abs(Adense[i, j]) - 1e-12 for i in 1:n, j in 1:n) + @test cover_lobjective(a, A) ≈ cover_lobjective(a, Adense) + @test cover_qobjective(a, A) ≈ cover_qobjective(a, Adense) + end + d, a = symdiagcover(A) + @test all(a[i] * a[j] >= abs(Adense[i, j]) - 1e-12 for i in 1:n, j in 1:n if i != j) + @test all(a[i]^2 + d[i] >= abs(Adense[i, i]) - 1e-12 for i in 1:n) + @test all(d[i] >= 0 for i in 1:n) + end + end + + @testset "Adjoint and Transpose" begin + for Adense in ([1.0 2.0; 3.0 4.0], [1.0 2.0 3.0; 4.0 5.0 6.0]) + for wrapper in (adjoint, transpose) + A = wrapper(Adense) + Adense_wrap = Matrix(A) + a, b = cover(A) + @test all(a[i] * b[j] >= abs(Adense_wrap[i, j]) - 1e-12 + for i in axes(Adense_wrap, 1), j in axes(Adense_wrap, 2)) + @test cover_lobjective(a, b, A) ≈ cover_lobjective(a, b, Adense_wrap) + @test cover_qobjective(a, b, A) ≈ cover_qobjective(a, b, Adense_wrap) + # Objectives are the same as computing cover on the parent and swapping + a0, b0 = cover(Adense) + @test cover_lobjective(a, b, A) ≈ cover_lobjective(a0, b0, Adense) + end + end + end + end + @testset "quality vs optimal (testmatrices)" begin if !isdefined(@__MODULE__, :symmetric_matrices) || !isdefined(@__MODULE__, :general_matrices) include("testmatrices.jl") # defines symmetric_matrices and general_matrices