Close some solver discrepancies

docs/src/index.md

@@ -140,11 +140,18 @@ You can override the default penalty by supplying it as an argument to the solve
 | [`cover`](@ref) | no | hard (`r ≤ 1`) | heuristic | — |
 | [`symcover_min`](@ref) | yes | hard (`r ≤ 1`) | `AbsLog{2}` (or `AbsLog{1}`, `AbsLinear`) | native for `AbsLog{2}`; else JuMP |
 | [`cover_min`](@ref) | no | hard (`r ≤ 1`) | `AbsLog{2}` (or `AbsLog{1}`, `AbsLinear`) | native for `AbsLog{2}`; else JuMP |
-| [`soft_symcover`](@ref) | yes | soft (penalized) | `AbsLinear{2}` (or `AbsLog`, `AbsLinear{1}`) | — |
-| [`soft_cover`](@ref) | no | soft (penalized) | `AbsLinear{2}` (or `AbsLinear{1}`) | — |
-| [`soft_symcover_min`](@ref) | yes | soft (penalized) | `AbsLog{2}`, `AbsLinear` | JuMP |
+| [`soft_symcover`](@ref) | yes | soft (penalized) | `AbsLinear{2}` (or `AbsLog`, `AbsLinear{1}`) | native for `AbsLog`; else — |
+| [`soft_cover`](@ref) | no | soft (penalized) | `AbsLinear{2}` (or `AbsLog`, `AbsLinear{1}`) | native for `AbsLog`; else — |
+| [`soft_symcover_min`](@ref) | yes | soft (penalized) | `AbsLog{2}`, `AbsLinear` | native for `AbsLog{2}`; else JuMP |
 | [`soft_cover_min`](@ref) | no | soft (penalized) | `AbsLog{2}`, `AbsLinear` | native for `AbsLog{2}`; else JuMP |
 
+Under `AbsLog{2}` the soft objective is convex with a single minimizer, so
+[`soft_symcover`](@ref) and [`soft_symcover_min`](@ref) are the same function, as are
+[`soft_cover`](@ref) and [`soft_cover_min`](@ref): there is nothing for a heuristic and a
+minimizer to disagree about. Under `AbsLog{1}` they part company — the soft `AbsLog{1}`
+covers are coordinate descents that reach a deterministic fixed point rather than a
+minimizer, and `soft_symcover_min`/`soft_cover_min` do not yet accept `AbsLog{1}`.
+
 **[`symcover`](@ref), [`cover`](@ref), and any native implementation can be recommended for production use,**
 possibly with relaxed convergence bounds.
 The heuristic solvers are particularly fast: they run in ``O(mn)`` time for an

ext/SIAJuMP.jl

@@ -207,42 +207,5 @@ function _cover_min_abslog1(A, start)
     return a, b
 end
 
-# Soft (unconstrained) symmetric cover: minimize ∑ (log r_ij)² with no constraints.
-# The objective is quadratic in α = log a, solved as a QP. Convex with a unique minimizer,
-# so no multistart is needed and the start — when the refiner supplies one — is a hint the
-# result does not record.
-ScaleInvariantAnalysis.soft_symcover_min(::AbsLog{2}, A) = _soft_symcover_min_abslog2(A, nothing)
-
-function ScaleInvariantAnalysis.soft_symcover_min!(::AbsLog{2}, a::AbstractVector, A)
-    ScaleInvariantAnalysis._prepare_soft_symcover_start!(a, A)
-    a .= _soft_symcover_min_abslog2(A, a)
-    return a
-end
-
-function _soft_symcover_min_abslog2(A, start)
-    axr = axes(A, 1)
-    axes(A, 2) == axr || throw(ArgumentError("soft_symcover_min requires a square matrix"))
-    T = float(real(eltype(A)))
-    pr = collect(axr)
-    n = length(pr)
-    Apos = [A[pr[i], pr[j]] for i in 1:n, j in 1:n]
-    logA = log.(abs.(Apos))
-    supported = [any(!iszero, @view Apos[i, :]) || any(!iszero, @view Apos[:, i]) for i in 1:n]
-    model = JuMP.Model(HiGHS.Optimizer)
-    JuMP.set_silent(model)
-    if start === nothing
-        @variable(model, α[1:n])
-    else
-        α0 = [supported[k] ? log(T(start[pr[k]])) : zero(T) for k in 1:n]
-        @variable(model, α[k=1:n], start = α0[k])
-    end
-    @objective(model, Min, sum(abs2, α[i] + α[j] - logA[i, j] for i in 1:n, j in 1:n if Apos[i, j] != 0))
-    JuMP.optimize!(model)
-    a = similar(Array{T}, axr)
-    for (i, k) in pairs(pr)
-        a[k] = supported[i] ? exp(JuMP.value(α[i])) : zero(T)
-    end
-    return a
-end
 
 end

ext/SIASparseArrays.jl

@@ -71,27 +71,24 @@ end
 # is the intended path when nnz ≪ n²; pass `linsolve=:auto`/`:dense` to force the
 # dense factorization. Only AbsLog{2} is native; other penalties dispatch to the
 # JuMP extension.
-# The worker allocates its scale vectors with `similar(A, ...)`, which is a
-# `SparseVector` for a sparse `A`; the scales are dense objects, so return plain
-# `Vector`s, matching `cover`/`symcover` on the same input.
 function ScaleInvariantAnalysis.symcover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...)
     a, _ = _symcover_min_abslog2(A; linsolve, kwargs...)
-    return Vector(a)
+    return a
 end
 
 function ScaleInvariantAnalysis.cover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...)
     a, b, _ = _cover_min_abslog2(A; linsolve, kwargs...)
-    return Vector(a), Vector(b)
+    return a, b
 end
 
 function ScaleInvariantAnalysis.symcover_min(ϕ::AbsLog{2}, S::Symmetric{<:Any, <:SparseMatrixCSC}; linsolve::Symbol=:lsqr, kwargs...)
     a, _ = _symcover_min_abslog2(S; linsolve, kwargs...)
-    return Vector(a)
+    return a
 end
 
 function ScaleInvariantAnalysis.symcover_min(ϕ::AbsLog{2}, H::Hermitian{<:Any, <:SparseMatrixCSC}; linsolve::Symbol=:lsqr, kwargs...)
     a, _ = _symcover_min_abslog2(H; linsolve, kwargs...)
-    return Vector(a)
+    return a
 end
 
 # The refiners take the same sparse `linsolve` default as the solvers above.

src/ScaleInvariantAnalysis.jl

@@ -75,11 +75,8 @@ function __init__()
             printstyled(io, "\nAbsLog{2} is solved natively; other penalties require loading JuMP plus HiGHS (for AbsLog{1}) or Ipopt (for AbsLinear)."; color=:yellow)
             return true
         end
-        if exc.f === soft_symcover_min || exc.f === soft_symcover_min!
-            printstyled(io, "\nThis method requires loading JuMP plus HiGHS (for AbsLog{2}) or Ipopt (for AbsLinear)."; color=:yellow)
-            return true
-        end
-        if exc.f === soft_cover_min || exc.f === soft_cover_min!
+        if exc.f === soft_symcover_min || exc.f === soft_symcover_min! ||
+           exc.f === soft_cover_min || exc.f === soft_cover_min!
             printstyled(io, "\nAbsLog{2} is solved natively; AbsLinear penalties require loading JuMP plus Ipopt. AbsLog{1} is not yet supported."; color=:yellow)
             return true
         end

src/initializers.jl

@@ -26,9 +26,10 @@ starting point does not depend on the objective it will be refined against.
 
 `strategy` names the point:
 
-- `:geomean` — the AbsLog{2} unconstrained minimum, the geometric mean of the
-  nonzero entries of each row. This is the minimizer of the soft AbsLog{2}
-  objective, and is *not* a cover.
+- `:geomean` — the geometric mean of the nonzero entries of each row, and *not* a
+  cover. It minimizes the soft AbsLog{2} objective exactly when every entry of `A`
+  is nonzero; on a sparse support it approximates that minimum, which
+  [`soft_symcover_min`](@ref)`(AbsLog{2}(), A)` returns exactly.
 - `:leaveout` — the geometric mean recomputed with the most-underweighted
   support entry dropped, which lands in the basin that treats that entry as
   effectively zero. Raises an `ArgumentError` when no entry can be dropped

src/minimal_covers.jl

@@ -380,7 +380,8 @@ end
 # indexed like `axes(A, 1)` and supplies the first iterate in place of the cold
 # unweighted solve; the objective is convex, so it changes the path but not the result.
 function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
-                               maxiter::Int=40, linsolve::Symbol=:auto, start=nothing)
+                               maxiter::Int=40, linsolve::Symbol=:auto, start=nothing,
+                               boost::Bool=true)
     linsolve in (:auto, :dense, :lsqr) ||
         throw(ArgumentError("linsolve must be :auto, :dense, or :lsqr; got :$linsolve"))
     ax = axes(A, 1)
@@ -503,12 +504,17 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
         end
     end
     # Uniform boost to exact feasibility: α_i + α_j ≥ log|A_ij| for all support.
+    # `boost=false` leaves the iterate untouched, for the soft objective, which
+    # imposes no coverage constraint and whose optimum the boost would move off.
     γ = zero(T)
-    for jp in 1:n, ip in 1:n
-        S[ip, jp] || continue
-        γ = max(γ, (C[ip, jp] - α[ip] - α[jp]) / 2)
+    if boost
+        for jp in 1:n, ip in 1:n
+            S[ip, jp] || continue
+            γ = max(γ, (C[ip, jp] - α[ip] - α[jp]) / 2)
+        end
     end
-    a = similar(A, T, ax)
+    # Dense scale vector matching cover/symcover; `similar(A, …)` is a SparseVector for sparse A.
+    a = similar(Array{T}, ax)
     for (ip, i) in enumerate(ax)
         a[i] = hassupp[ip] ? exp(α[ip] + γ) : zero(T)
     end
@@ -520,7 +526,8 @@ end
 # `start`, when given, is a positive cover `(a, b)` indexed like the rows and columns
 # of `A`, supplying the first iterate in place of the cold unweighted solve.
 function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
-                            maxiter::Int=40, linsolve::Symbol=:auto, start=nothing)
+                            maxiter::Int=40, linsolve::Symbol=:auto, start=nothing,
+                            boost::Bool=true)
     linsolve in (:auto, :dense, :lsqr) ||
         throw(ArgumentError("linsolve must be :auto, :dense, or :lsqr; got :$linsolve"))
     axr = axes(A, 1)
@@ -684,13 +691,19 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
         end
     end
     # Uniform boost to exact feasibility: α_i + β_j ≥ log|A_ij| on the support.
-    γ = zero(T)
-    for jp in 1:n, ip in 1:m
-        S[ip, jp] || continue
-        γ = max(γ, (C[ip, jp] - x[ip] - x[m+jp]) / 2)
-    end
-    for p in 1:N
-        x[p] += γ
+    # `boost=false` leaves the iterate untouched, for the soft objective, which
+    # imposes no coverage constraint and whose optimum the boost would move off.
+    # The balance shift below still applies: the gauge is a convention, not a
+    # constraint, and every cover this package returns satisfies it.
+    if boost
+        γ = zero(T)
+        for jp in 1:n, ip in 1:m
+            S[ip, jp] || continue
+            γ = max(γ, (C[ip, jp] - x[ip] - x[m+jp]) / 2)
+        end
+        for p in 1:N
+            x[p] += γ
+        end
     end
     # Shift along the (e; -e) gauge to the balance convention ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ.
     nnz = count(S)
@@ -703,8 +716,9 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
         Lβ += count(@view S[:, jp]) * x[m+jp]
     end
     s = nnz > 0 ? (Lβ - Lα) / (2 * nnz) : zero(T)
-    a = similar(A, T, axr)
-    b = similar(A, T, axc)
+    # Dense scale vectors matching cover/symcover; `similar(A, …)` is a SparseVector for sparse A.
+    a = similar(Array{T}, axr)
+    b = similar(Array{T}, axc)
     for (ip, i) in enumerate(axr)
         a[i] = hasrow[ip] ? exp(x[ip] + s) : zero(T)
     end
@@ -714,6 +728,21 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8),
     return a, b, (; nsolves=nsolves[], lsqriters=nlsqr[], linsolve=(use_lsqr ? :lsqr : :dense))
 end
 
+# Workers for the soft (unconstrained) AbsLog{2} covers. The soft objective
+# `∑_{i,j∈S} (log a_i + log a_j - log|A_ij|)²` is the hard workers' reweighted
+# least-squares problem with every weight held at 1, which is the cold solve they
+# already take as their first iterate: `κs=()` runs no penalty continuation, and
+# `boost=false` keeps the unconstrained minimizer where it is. It is convex, so one
+# linear solve settles it — no iteration and no multistart, unlike the non-convex
+# `AbsLinear` soft covers.
+#
+# Both paths inherit the hard workers' handling of a singular signless Laplacian (the
+# `[0 1; 1 0]` support graph among them) and of support-free rows and columns.
+_soft_symcover_min_abslog2(A::AbstractMatrix; kwargs...) =
+    _symcover_min_abslog2(A; κs=(), boost=false, kwargs...)
+_soft_cover_min_abslog2(A::AbstractMatrix; kwargs...) =
+    _cover_min_abslog2(A; κs=(), boost=false, kwargs...)
+
 # Internal exact reference implemented by the SIAJuMP extension; used only to
 # cross-check the native `symcover_min(::AbsLog{2})` in the test suite.
 function symcover_min_jump end