diff --git a/docs/src/index.md b/docs/src/index.md index a910ad4..ddc1508 100644 --- a/docs/src/index.md +++ b/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 diff --git a/ext/SIAJuMP.jl b/ext/SIAJuMP.jl index ce4dda8..837807c 100644 --- a/ext/SIAJuMP.jl +++ b/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 diff --git a/ext/SIASparseArrays.jl b/ext/SIASparseArrays.jl index b9a476e..fbf453a 100644 --- a/ext/SIASparseArrays.jl +++ b/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. diff --git a/src/ScaleInvariantAnalysis.jl b/src/ScaleInvariantAnalysis.jl index 8c03236..8f51630 100644 --- a/src/ScaleInvariantAnalysis.jl +++ b/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 diff --git a/src/initializers.jl b/src/initializers.jl index ad1541d..d61ff9d 100644 --- a/src/initializers.jl +++ b/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 diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 1f787b7..30f2692 100644 --- a/src/minimal_covers.jl +++ b/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 diff --git a/src/soft_covers.jl b/src/soft_covers.jl index f9ac5a4..513e4e7 100644 --- a/src/soft_covers.jl +++ b/src/soft_covers.jl @@ -17,11 +17,16 @@ Unlike [`symcover`](@ref), there is no hard coverage constraint: `a[i]*a[j]` may less than `|A[i,j]|`, with violations penalized by `ϕ`. Supported penalty functions: -- `AbsLog{2}()`: returns the analytical unconstrained minimum (no iterations needed). +- `AbsLog{2}()`: convex, and returns its exact unconstrained minimum from a single linear + solve. Identical to [`soft_symcover_min`](@ref)`(AbsLog{2}(), A)` — with one minimizer + there is nothing for a heuristic and a minimizer to disagree about. - `AbsLog{1}()`: initializes from the AbsLog{2} minimum, then refines by coordinate descent - with a log-space weighted-median step. The AbsLog{1} objective has a flat basin of equally - good minima; this returns the deterministic, scale-covariant representative reached by - coordinate descent from the AbsLog{2} minimum. + with a log-space weighted-median step, reaching a deterministic and scale-covariant fixed + point. That point is not in general a minimizer: each step minimizes exactly over one + coordinate, but the objective's nonsmoothness couples `a[i]` with `a[j]`, so the descent + can settle where no single-coordinate move improves and the objective still sits + materially above its minimum. [`soft_symcover_min`](@ref) does not yet offer an exact + `AbsLog{1}` alternative. - `AbsLinear{2}()` (default): non-convex; refined by coordinate descent from `starts` scale-covariant starting points, keeping the lowest-objective result (see below). - `AbsLinear{1}()`: initializes from the `AbsLinear{2}()` result, coordinate descent uses a @@ -64,23 +69,12 @@ julia> round.(soft_symcover([0 1; 1 0]); digits=4) """ soft_symcover(A::AbstractMatrix; kwargs...) = soft_symcover(AbsLinear{2}(), A; kwargs...) -function soft_symcover(ϕ::AbsLog{2}, A::AbstractMatrix) - ax = axes(A, 1) - axes(A, 2) == ax || throw(ArgumentError("soft_symcover requires a square matrix")) - T = float(real(eltype(A))) - # Dense scale vector matching cover/symcover; `similar(A, …)` is a SparseVector for sparse A. - a = similar(Array{T}, ax) - unconstrained_min!(ϕ, a, A) # analytical minimum; no iterations needed - return a -end +# The soft AbsLog{2} objective is convex with one minimizer, so the heuristic and the +# minimizer coincide: both are this solve. +soft_symcover(::AbsLog{2}, A::AbstractMatrix; kwargs...) = soft_symcover_min(AbsLog{2}(), A; kwargs...) function soft_symcover(::AbsLog{1}, A::AbstractMatrix; maxiter::Int=20) - ax = axes(A, 1) - axes(A, 2) == ax || throw(ArgumentError("soft_symcover requires a square matrix")) - T = float(real(eltype(A))) - # Dense scale vector matching cover/symcover; `similar(A, …)` is a SparseVector for sparse A. - a = similar(Array{T}, ax) - unconstrained_min!(AbsLog{2}(), a, A) # convex AbsLog{2} minimum: a good start + a = soft_symcover_min(AbsLog{2}(), A) # convex AbsLog{2} minimum: a good start _abslog1_iter!(a, A, maxiter) return a end @@ -118,6 +112,16 @@ Unlike [`cover`](@ref), there is no hard coverage constraint: `a[i]*b[j]` may be `|A[i,j]|`, with violations penalized by `ϕ`. Supported penalty functions: +- `AbsLog{2}()`: convex, and returns its exact unconstrained minimum from a single linear + solve. Identical to [`soft_cover_min`](@ref)`(AbsLog{2}(), A)` — with one minimizer + there is nothing for a heuristic and a minimizer to disagree about. +- `AbsLog{1}()`: initializes from the `AbsLog{2}()` minimum, then refines by alternating + weighted-median row and column updates, reaching a deterministic and scale-covariant fixed + point. As in [`soft_symcover`](@ref), that point is not in general a minimizer: each + half-sweep minimizes exactly, but the objective's nonsmoothness couples `a[i]` with `b[j]`, + so the descent can settle where no such sweep improves and the objective still sits + materially above its minimum. [`soft_cover_min`](@ref) does not yet offer an exact + `AbsLog{1}` alternative. - `AbsLinear{2}()` (default): in the inverse-scale variables `u = 1 ./ a`, `v = 1 ./ b`, the objective `∑_{i,j∈S} (1 - |A[i,j]| u[i] v[j])²` (sum over the nonzero support `S`) is biconvex, so alternating least squares with the closed-form half-sweeps @@ -162,6 +166,16 @@ julia> a * b' """ soft_cover(A::AbstractMatrix; kwargs...) = soft_cover(AbsLinear{2}(), A; kwargs...) +# The soft AbsLog{2} objective is convex with one minimizer, so the heuristic and the +# minimizer coincide: both are this solve. +soft_cover(::AbsLog{2}, A::AbstractMatrix; kwargs...) = soft_cover_min(AbsLog{2}(), A; kwargs...) + +function soft_cover(::AbsLog{1}, A::AbstractMatrix; maxiter::Int=20) + a, b = soft_cover_min(AbsLog{2}(), A) # convex AbsLog{2} minimum: a good start + _abslog1_iter_asym!(a, b, A, maxiter) + return _balance_cover!(a, b, A) +end + # Sole owner of the starts/σ/rng defaults for the AbsLinear{2} soft-cover family; # every other method in that family (the no-ϕ wrapper, the AbsLinear{1} method) # forwards them via `kwargs...` rather than restating the default. @@ -188,23 +202,38 @@ with no coverage constraints. The no-ϕ form defaults to `AbsLinear{2}()`, match [`soft_symcover`](@ref). Supported ϕ values and required extensions: -- `AbsLog{2}()`: requires JuMP and HiGHS. Convex, so the minimizer is unique. +- `AbsLog{2}()`: solved natively (no external solver). In log space the objective is a + linear least-squares, so one solve settles it, and being convex it has a unique + minimizer that no start can influence. `linsolve` selects the inner solve, exactly as + in [`symcover_min`](@ref). - `AbsLinear{1}()`, `AbsLinear{2}()`: requires JuMP and Ipopt. These objectives are non-convex, so the solver returns the minimum of the basin it starts in. Rather than commit to one start, these methods refine each of `strategies` — the [`initialize_symcover`](@ref) menu, by default `$(SYMCOVER_MIN_STRATEGIES)`, without forcing feasibility — and return the best cover found, at a cost of one solve per start. -- `AbsLog{1}()`: not yet implemented. +- `AbsLog{1}()`: not yet implemented. The objective is an LP in log space, but its optimum + is a face, and the lexicographic AbsLog{2} selection that [`symcover_min`](@ref) uses to + pin one member of the corresponding hard face does not carry over: the hard face is bounded + by the coverage constraints, while this one is a level set of an unconstrained piecewise- + linear objective, across which the quadratic pulls far enough to cost most of the exactly + tight residuals that make `AbsLog{1}` worth choosing. See also: [`soft_symcover_min!`](@ref), [`soft_symcover`](@ref), [`symcover_min`](@ref). """ function soft_symcover_min end soft_symcover_min(A::AbstractMatrix; kwargs...) = soft_symcover_min(AbsLinear{2}(), A; kwargs...) +function soft_symcover_min(::AbsLog{2}, A::AbstractMatrix; kwargs...) + ax = axes(A, 1) + axes(A, 2) == ax || throw(ArgumentError("soft_symcover_min requires a square matrix")) + a, _ = _soft_symcover_min_abslog2(A; kwargs...) + return a +end + # Multistart driver for the non-convex soft AbsLinear covers, the unconstrained counterpart # of `symcover_min(::AbsLinear)`. The kernels (`soft_symcover_min!`) live in SIAIpopt; the # menu and the selection are native. Starts are taken raw: a soft cover is under no -# obligation to cover `A`, and the raw geometric mean is the exact soft AbsLog{2} optimum. +# obligation to cover `A`. function soft_symcover_min(ϕ::AbsLinear, A::AbstractMatrix; strategies=SYMCOVER_MIN_STRATEGIES) ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("soft_symcover_min requires a square matrix")) @@ -244,6 +273,12 @@ function soft_symcover_min! end soft_symcover_min!(a::AbstractVector, A::AbstractMatrix; kwargs...) = soft_symcover_min!(AbsLinear{2}(), a, A; kwargs...) +function soft_symcover_min!(::AbsLog{2}, a::AbstractVector, A::AbstractMatrix; kwargs...) + _prepare_soft_symcover_start!(a, A) + a .= soft_symcover_min(AbsLog{2}(), A; kwargs...) # convex: the start is not read + return a +end + # Shared prologue of the `soft_symcover_min!` kernels. The soft objective constrains # nothing, so — unlike `_prepare_symcover_start!` — this checks positivity only, and moves # the start nowhere. @@ -292,7 +327,12 @@ Supported ϕ values and required extensions: cover found, at a cost of one solve per start. The result is the best *local* minimum on that menu: the multistart is a hedge against a poor basin, not a certificate of global optimality. -- `AbsLog{1}()`: not yet implemented. +- `AbsLog{1}()`: not yet implemented. The objective is an LP in log space, but its optimum + is a face, and the lexicographic AbsLog{2} selection that [`symcover_min`](@ref) uses to + pin one member of the corresponding hard face does not carry over: the hard face is bounded + by the coverage constraints, while this one is a level set of an unconstrained piecewise- + linear objective, across which the quadratic pulls far enough to cost most of the exactly + tight residuals that make `AbsLog{1}` worth choosing. Every start on the menu co-varies with a rescaling of `A` and the objective is scale-invariant, so the selection — and hence the result — is scale-covariant. @@ -302,11 +342,8 @@ See also: [`soft_cover_min!`](@ref), [`soft_symcover_min`](@ref), [`soft_cover`] function soft_cover_min end soft_cover_min(A::AbstractMatrix; kwargs...) = soft_cover_min(AbsLinear{2}(), A; kwargs...) -function soft_cover_min(::AbsLog{2}, A::AbstractMatrix) - T = float(real(eltype(A))) - a = similar(Array{T}, axes(A, 1)) - b = similar(Array{T}, axes(A, 2)) - unconstrained_min!(AbsLog{2}(), a, b, A) +function soft_cover_min(::AbsLog{2}, A::AbstractMatrix; kwargs...) + a, b, _ = _soft_cover_min_abslog2(A; kwargs...) return a, b end @@ -346,9 +383,11 @@ function soft_cover_min! end soft_cover_min!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) = soft_cover_min!(AbsLinear{2}(), a, b, A; kwargs...) -function soft_cover_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix) +function soft_cover_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) _prepare_soft_cover_start!(a, b, A) - unconstrained_min!(AbsLog{2}(), a, b, A) # analytic minimum; the start is not needed + anew, bnew = soft_cover_min(AbsLog{2}(), A; kwargs...) # convex: the start is not read + a .= anew + b .= bnew return a, b end @@ -695,6 +734,73 @@ function _abslog1_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; tol: return a end +# Alternating weighted-median descent for the asymmetric AbsLog{1} soft cover, working in +# log space (α = log a, β = log b). Updating α[i] with β fixed minimizes +# ∑_{j: A[i,j]≠0} |α[i] + β[j] - log|A[i,j]||, whose minimizer is the median of the points +# log|A[i,j]| - β[j], one per nonzero entry of row i; the β-update is dual. Row and column +# scales are distinct variables, so no term is self-coupled — the symmetric solver's +# double-weighted diagonal has no counterpart here — and each half-sweep is an exact block +# minimization. The AbsLog{1} minimum is a flat basin; the lower median is chosen for a +# deterministic, scale-covariant result. +# +# Each half-sweep minimizes exactly, but the objective's nonsmoothness couples α[i] with +# β[j], so a fixed point of the sweeps need not minimize it. See `soft_cover`. +# +# `iter` bounds the sweeps; the descent exits early once the largest relative coordinate +# movement in a sweep drops to `tol` (scale-invariant, so covariant restarts of a rescaled +# problem exit on the same sweep). Rows/columns with empty support keep scale 0. +function _abslog1_iter_asym!(a::AbstractVector{T}, b::AbstractVector{T}, A::AbstractMatrix, + iter::Int; tol::Real=1e-12) where T + axr, axc = axes(A, 1), axes(A, 2) + eachindex(a) == axr || throw(DimensionMismatch("row indices of `A` must match `a`, got $(axr) vs $(eachindex(a))")) + eachindex(b) == axc || throw(DimensionMismatch("column indices of `A` must match `b`, got $(axc) vs $(eachindex(b))")) + bufc = Vector{T}(undef, length(axc)) # log-points for an a-row update + bufr = Vector{T}(undef, length(axr)) # log-points for a b-column update + for _ in 1:iter + maxrel = zero(T) + for i in axr + iszero(a[i]) && continue # unsupported rows/columns stay at zero + nc = 0 + for j in axc + Aij = T(abs(A[i, j])) + iszero(Aij) && continue + bj = b[j] + iszero(bj) && continue + bufc[nc += 1] = log(Aij) - log(bj) + end + nc == 0 && continue + c = view(bufc, 1:nc) + sort!(c) + x = exp(c[(nc + 1) ÷ 2]) # lower median + ai = a[i] + den = max(abs(x), abs(ai)) + iszero(den) || (maxrel = max(maxrel, abs(x - ai) / den)) + a[i] = x + end + for j in axc + iszero(b[j]) && continue + nr = 0 + for i in axr + Aij = T(abs(A[i, j])) + iszero(Aij) && continue + ai = a[i] + iszero(ai) && continue + bufr[nr += 1] = log(Aij) - log(ai) + end + nr == 0 && continue + c = view(bufr, 1:nr) + sort!(c) + x = exp(c[(nr + 1) ÷ 2]) + bj = b[j] + den = max(abs(x), abs(bj)) + iszero(den) || (maxrel = max(maxrel, abs(x - bj) / den)) + b[j] = x + end + maxrel <= T(tol) && break + end + return a, b +end + # Labeled `(a, b)` candidate starts for the asymmetric AbsLinear{2} multistart, in selection # order. Deterministic starts: the boosted geometric mean (also the perturbation base) and the # tightened hard cover, obtained by tightening a copy of the former so the shared passes run diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 9e391b9..2b9f84c 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -192,12 +192,14 @@ end end @testset "soft_cover_min native AbsLog{2}" begin - # Matches the analytic asymmetric minimizer directly. + # `A` has no zero entry, the case in which the geometric mean coincides with the + # minimum; the two compute it differently, so they agree to roundoff, not bitwise. + # On a sparse support they part company -- see the oracle in `test/soft_covers.jl`. A = [1.0 2.0 3.0; 6.0 5.0 4.0] a, b = soft_cover_min(AbsLog{2}(), A) a_ref, b_ref = similar(a), similar(b) ScaleInvariantAnalysis.unconstrained_min!(AbsLog{2}(), a_ref, b_ref, A) - @test a == a_ref && b == b_ref + @test a ≈ a_ref && b ≈ b_ref # It's the unconstrained minimum: any perturbation can only raise the objective. obj0 = cover_objective(AbsLog{2}(), a, b, A) diff --git a/test/soft_covers.jl b/test/soft_covers.jl index b00f0e5..1b65fc1 100644 --- a/test/soft_covers.jl +++ b/test/soft_covers.jl @@ -167,8 +167,51 @@ end @test isfinite(cover_objective(AbsLinear{1}(), az, bz, Az1)) end - # AbsLog penalties are unsupported for soft_cover. - @test_throws MethodError soft_cover(AbsLog{2}(), A) + # AbsLog{2} is the convex soft cover, and identical to its minimizer. + @test soft_cover(AbsLog{2}(), A) == soft_cover_min(AbsLog{2}(), A) +end + +@testset "soft_cover AbsLog{1}" begin + # Each half-sweep is an exact block minimization, so the descent never worsens the + # AbsLog{2} start it refines. + for M in ([1.0 2.0 3.0; 6.0 5.0 4.0], + [2.0 1.0 0.5; 0.1 4.0 3.0; 1.0 2.0 0.2; 5.0 0.3 1.0], + [1.0 2.0 0.0 4.0; 0.0 5.0 6.0 1.0; 3.0 0.0 2.0 8.0]) + a0, b0 = soft_cover_min(AbsLog{2}(), M) + a, b = soft_cover(AbsLog{1}(), M) + @test cover_objective(AbsLog{1}(), a, b, M) <= cover_objective(AbsLog{1}(), a0, b0, M) + 1e-12 + @test isbalanced(a, b, M) + @test a isa Vector{Float64} && b isa Vector{Float64} + end + + # A rank-1 matrix is exactly coverable, so the L1 objective reaches 0. + A1 = [2.0, 0.5, 3.0] * [1.0, 4.0, 0.25, 2.0]' + a, b = soft_cover(AbsLog{1}(), A1) + @test cover_objective(AbsLog{1}(), a, b, A1) ≈ 0 atol=1e-10 + + # Deterministic, and scale-covariant in the product under row/column rescaling. + @test soft_cover(AbsLog{1}(), A1) == soft_cover(AbsLog{1}(), A1) + Ac = [2.0 1.0 0.5; 0.1 4.0 3.0; 1.0 2.0 0.2; 5.0 0.3 1.0] + dr = [3.0, 0.5, 2.0, 0.25]; dc = [4.0, 0.1, 1.5] + @test covaries(A -> soft_cover(AbsLog{1}(), A), Ac, dr, dc; rtol=1e-8) + + # An entirely-zero row keeps scale 0 and the rest stays finite. + Az = [0.0 0.0 0.0; 1.0 2.0 3.0; 4.0 0.0 5.0] + az, bz = soft_cover(AbsLog{1}(), Az) + @test az[1] == 0 + @test all(isfinite, az) && all(isfinite, bz) + + # On a symmetric matrix this does not reduce to `soft_symcover`: freeing `a` from `b` + # relaxes the problem, so the two descents are minimizing over different sets and their + # fixed points differ. Only exact coverability forces them to agree. + S = [4.0 1.0 2.0; 1.0 9.0 3.0; 2.0 3.0 16.0] + S1 = (v = [2.0, 0.5, 3.0]; v * v') + a, b = soft_cover(AbsLog{1}(), S1) + @test a .* b' ≈ soft_symcover(AbsLog{1}(), S1) .* soft_symcover(AbsLog{1}(), S1)' rtol=1e-8 + + # The minimizers stay unimplemented; the heuristic is not a stand-in for one. + @test_throws MethodError soft_cover_min(AbsLog{1}(), A1) + @test_throws MethodError soft_symcover_min(AbsLog{1}(), S) end @testset "AbsLinear soft-cover multistart" begin @@ -313,3 +356,87 @@ end # than merely documenting. @test covaries(soft_cover, A, dr, dc; rtol=1e-9) end + +@testset "soft AbsLog{2} is the exact unconstrained minimum" begin + # The soft AbsLog{2} objective is a linear least-squares in log space, so an + # oracle needs no solver: `pinv(M) * z` settles it directly. `M` always carries + # the (e; −e) gauge null direction in the asymmetric case and can be singular in + # the symmetric one (bipartite support), so compare objectives — which the gauge + # cannot move — rather than the scale vectors. + function exact_sym(A) + n = size(A, 1) + S = [(i, j) for i in 1:n, j in 1:n if !iszero(A[i, j])] + M = zeros(length(S), n) + z = zeros(length(S)) + for (k, (i, j)) in enumerate(S) + M[k, i] += 1 + M[k, j] += 1 + z[k] = log(abs(A[i, j])) + end + return exp.(pinv(M) * z) + end + function exact_asym(A) + m, n = size(A) + S = [(i, j) for i in 1:m, j in 1:n if !iszero(A[i, j])] + M = zeros(length(S), m + n) + z = zeros(length(S)) + for (k, (i, j)) in enumerate(S) + M[k, i] = 1 + M[k, m+j] = 1 + z[k] = log(abs(A[i, j])) + end + x = pinv(M) * z + return exp.(x[1:m]), exp.(x[m+1:end]) + end + + # A zero entry is what separates the exact minimum from the geometric mean; a + # fully supported `A` cannot tell them apart. + sym_zeros = [4.0 1.0 0.0 2.0; 1.0 9.0 3.0 0.0; 0.0 3.0 1.0 5.0; 2.0 0.0 5.0 16.0] + asym_zeros = [1.0 2.0 0.0 4.0; 0.0 5.0 6.0 1.0; 3.0 0.0 2.0 8.0] + rng = StableRNG(11) + sym_dense = (M = rand(rng, 5, 5); (M + M') ./ 2) + asym_dense = rand(rng, 4, 6) + + for A in (sym_zeros, sym_dense, float.(last(symmetric_matrices[1]))) + a = soft_symcover_min(AbsLog{2}(), A) + @test cover_objective(AbsLog{2}(), a, A) ≈ cover_objective(AbsLog{2}(), exact_sym(A), A) rtol=1e-8 + # Stationarity of the convex objective: ∂/∂α_k ∑ (α_i + α_j − log|A_ij|)² = 0. + α = log.(a) + g = [sum(2 * (α[i] + α[j] - log(abs(A[i, j]))) * ((i == k) + (j == k)) + for i in axes(A, 1), j in axes(A, 2) if !iszero(A[i, j])) for k in axes(A, 1)] + @test maximum(abs, g) < 1e-8 * max(1, maximum(abs, α)) + end + for A in (asym_zeros, asym_dense, float.(last(general_matrices[1]))) + a, b = soft_cover_min(AbsLog{2}(), A) + ae, be = exact_asym(A) + @test cover_objective(AbsLog{2}(), a, b, A) ≈ cover_objective(AbsLog{2}(), ae, be, A) rtol=1e-8 + @test isbalanced(a, b, A) + end + + # The geometric mean is the minimum only on a full support. Were the solvers to + # fall back to it, the sparse cases above would silently regress. + @test cover_objective(AbsLog{2}(), initialize_symcover(sym_dense; strategy=:geomean, feasible=:none), sym_dense) ≈ + cover_objective(AbsLog{2}(), soft_symcover_min(AbsLog{2}(), sym_dense), sym_dense) rtol=1e-8 + @test cover_objective(AbsLog{2}(), initialize_symcover(sym_zeros; strategy=:geomean, feasible=:none), sym_zeros) > + cover_objective(AbsLog{2}(), soft_symcover_min(AbsLog{2}(), sym_zeros), sym_zeros) * (1 + 1e-6) + + # Convex with one minimizer, so the heuristic and the minimizer are the same + # function, and a refiner cannot be steered by its start. + @test soft_symcover(AbsLog{2}(), sym_zeros) == soft_symcover_min(AbsLog{2}(), sym_zeros) + @test soft_cover(AbsLog{2}(), asym_zeros) == soft_cover_min(AbsLog{2}(), asym_zeros) + for strategy in (:geomean, :leaveout, :diagfeasible, :hardcover) + a0 = initialize_symcover(sym_zeros; strategy, feasible=:none) + @test soft_symcover_min!(AbsLog{2}(), a0, sym_zeros) ≈ soft_symcover_min(AbsLog{2}(), sym_zeros) + end + a0, b0 = initialize_cover(asym_zeros; strategy=:geomean, feasible=:none) + ar, br = soft_cover_min!(AbsLog{2}(), a0, b0, asym_zeros) + @test (ar, br) == soft_cover_min(AbsLog{2}(), asym_zeros) + + # A bipartite support graph makes the signless Laplacian singular; the balanced + # representative is the one reported. + @test soft_symcover_min(AbsLog{2}(), [0 1; 1 0]) ≈ [1.0, 1.0] + + # The `:lsqr` and dense paths solve the same problem. + @test soft_symcover_min(AbsLog{2}(), sym_zeros; linsolve=:lsqr) ≈ + soft_symcover_min(AbsLog{2}(), sym_zeros; linsolve=:dense) rtol=1e-6 +end