From 8d6fa294823cd6c8da27b7fdae3e9d459cbe20bb Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 8 Mar 2026 10:20:51 -0500 Subject: [PATCH 1/2] Ensure initializer produce feasible solutions --- src/covers.jl | 24 +++++++++++++++++++++++- test/runtests.jl | 13 +++++++++++-- 2 files changed, 34 insertions(+), 3 deletions(-) diff --git a/src/covers.jl b/src/covers.jl index 4abab26..433aa1b 100644 --- a/src/covers.jl +++ b/src/covers.jl @@ -81,7 +81,7 @@ function symcover(A::AbstractMatrix; kwargs...) else aprod = ai * aj aprod >= Aij && continue - s = sqrt(Aij / sqrt(aprod)) + s = sqrt(Aij / aprod) a[i] = s * ai a[j] = s * aj end @@ -146,6 +146,28 @@ function cover(A::AbstractMatrix; kwargs...) for j in axes(A, 2) b[j] = iszero(nzb[j]) ? zero(T) : exp(logb[j] / nzb[j] - halfmu) end + # Now we have sums of (log(a[i]) + log(b[j]) - log(A[i, j])) to be zero across rows or columns. + # Now it needs to be boosted to cover A. + for j in axes(A, 2) + for i in axes(A, 1) + Aij, ai, bj = abs(A[i, j]), a[i], b[j] + if iszero(bj) + if !iszero(ai) + b[j] = Aij / ai + else + a[i] = b[j] = sqrt(Aij) + end + elseif iszero(ai) + a[i] = Aij / bj + else + aprod = ai * bj + aprod >= Aij && continue + s = sqrt(Aij / aprod) + a[i] = s * ai + b[j] = s * bj + end + end + end return tighten_cover!(a, b, A; kwargs...) end diff --git a/test/runtests.jl b/test/runtests.jl index cb4a973..523391c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -134,16 +134,25 @@ using Test sym_ratios = Float64[] for (_, A) in symmetric_matrices Af = Float64.(A) + 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(symcover(Af; iter=10), Af) + lfast = cover_lobjective(a, Af) iszero(lopt) || push!(sym_ratios, lfast / lopt) end - @test median(sym_ratios) < 1.01 + @test median(sym_ratios) < 1.1 # cover cover_lobjective should be close to optimal (cover_lmin) gen_ratios = Float64[] for (_, A) in general_matrices Af = Float64.(A) + 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) From 675cf7d4160d7e83e9f4e445ed4ea59ad13501c6 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 8 Mar 2026 10:27:49 -0500 Subject: [PATCH 2/2] fix docs --- docs/src/index.md | 6 +++--- src/covers.jl | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 5f653ea..bca7adc 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -79,15 +79,15 @@ julia> using ScaleInvariantAnalysis, JuMP, HiGHS julia> A = [1 2 3; 6 5 4]; julia> a, b = cover(A) -([1.2544610775677627, 3.4759059767492304], [1.7261686708831454, 1.6581941714076147, 2.391465190627206]) +([1.2674308473260654, 3.4759059767492304], [1.7261686708831454, 1.61137045961268, 2.366993044495631]) julia> aq, bq = cover_qmin(A) ([1.1986299970143055, 3.25358233351279], [1.8441211516912772, 1.6685716234216104, 2.5028574351324164]) julia> a * b' 2×3 Matrix{Float64}: - 2.16541 2.08014 3.0 - 6.0 5.76373 8.31251 + 2.1878 2.0423 3.0 + 6.0 5.60097 8.22745 julia> aq * bq' 2×3 Matrix{Float64}: diff --git a/src/covers.jl b/src/covers.jl index 433aa1b..cf050ee 100644 --- a/src/covers.jl +++ b/src/covers.jl @@ -108,12 +108,12 @@ See also: [`cover_lmin`](@ref), [`cover_qmin`](@ref), [`symcover`](@ref). julia> A = [1 2 3; 6 5 4]; julia> a, b = cover(A) -([1.2544610775677627, 3.4759059767492304], [1.7261686708831454, 1.6581941714076147, 2.391465190627206]) +([1.2674308473260654, 3.4759059767492304], [1.7261686708831454, 1.61137045961268, 2.366993044495631]) julia> a * b' 2×3 Matrix{Float64}: - 2.16541 2.08014 3.0 - 6.0 5.76373 8.31251 + 2.1878 2.0423 3.0 + 6.0 5.60097 8.22745 ``` """ function cover(A::AbstractMatrix; kwargs...)