Add soft_symcover! and soft_cover! refiners

docs/src/index.md

@@ -145,12 +145,29 @@ You can override the default penalty by supplying it as an argument to the solve
 | [`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}`
+The two soft tiers are separated by a different axis than the two hard ones.  For hard
+covers, [`symcover`](@ref) trades optimality for speed while guaranteeing feasibility, and
+[`symcover_min`](@ref) is optimal.  Both soft tiers minimize the same unconstrained
+objective, and differ instead in what they promise about reaching its minimum:
+
+- [`soft_symcover`](@ref) and [`soft_cover`](@ref) are **always native and best-effort.**
+  They require no extension for any penalty, and they own their multistart — but what they
+  return is a coordinate-descent fixed point, which for the non-convex and nonsmooth
+  penalties need not be a minimizer.
+- [`soft_symcover_min`](@ref) and [`soft_cover_min`](@ref) return a **true minimizer of the
+  basin they start in, and may require an extension.**  `AbsLog{2}` is native; the
+  `AbsLinear` penalties need JuMP and Ipopt; `AbsLog{1}` is not implemented.
+
+Both reduce to the same trade: cheap and always available, against best quality and
+possibly an extra dependency.
+
+Under `AbsLog{2}` the objective is convex with a single minimizer, so the tiers coincide —
+[`soft_symcover`](@ref) *is* [`soft_symcover_min`](@ref) there, and likewise for the
+asymmetric pair.  That is the degenerate case of the contract rather than an exception to
+it: with one minimizer there is nothing for a best-effort descent and a minimizer to
+disagree about.  Under `AbsLog{1}` they part company most sharply — 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}`.
+minimizer, and `soft_symcover_min`/`soft_cover_min` do not accept `AbsLog{1}` at all.
 
 **[`symcover`](@ref), [`cover`](@ref), and any native implementation can be recommended for production use,**
 possibly with relaxed convergence bounds.
@@ -208,8 +225,11 @@ a    = symcover_min(AbsLog{1}(), A)   # L1-minimal symmetric hard cover
 a, b = cover_min(AbsLog{1}(), A)      # L1-minimal general hard cover
 ```
 
-The [`soft_symcover_min`](@ref) soft solver is likewise JuMP-backed (HiGHS for
-`AbsLog{2}`, Ipopt for the `AbsLinear` penalties).
+The soft `*_min` solvers divide along the same line, but not at the same place:
+[`soft_symcover_min`](@ref) and [`soft_cover_min`](@ref) solve `AbsLog{2}` natively and
+reach for JuMP with Ipopt only for the `AbsLinear` penalties.  They do not accept
+`AbsLog{1}`; the soft `AbsLog{1}` covers are available through [`soft_symcover`](@ref) and
+[`soft_cover`](@ref), which are native.
 
 ### Uniqueness
 
@@ -250,15 +270,24 @@ At a lower level, this package's interface is organized in three layers:
   basins — which is why the choice is a named part of the start rather than an internal
   detail.  The heuristic [`cover`](@ref) is itself a composition of these: the geometric
   mean, boosted, then tightened.
-- **Refiners** improve a starting point in place.  [`symcover_min!`](@ref),
-  [`cover_min!`](@ref), and [`soft_symcover_min!`](@ref) validate the start, then optimize
-  from it.  Which basin they reach is the caller's choice, by construction.  The hard
-  refiners require a start that covers `A`; the soft one does not, since its objective
-  constrains nothing.
-- **Solvers** bundle the two.  [`symcover_min`](@ref), [`cover_min`](@ref), and
-  [`soft_symcover_min`](@ref) refine every start on a menu (the `strategies` keyword) and
-  return the best cover by [`cover_objective`](@ref), so their result depends on `A` and not
-  on an initialization the caller never chose.
+- **Refiners** improve a starting point in place, and are the `!`-suffixed forms of the
+  solvers: [`symcover_min!`](@ref), [`cover_min!`](@ref), [`soft_symcover!`](@ref),
+  [`soft_cover!`](@ref), [`soft_symcover_min!`](@ref), and [`soft_cover_min!`](@ref)
+  validate the start, then optimize from it.  Which basin they reach is the caller's
+  choice, by construction, and supplying the start is the caller's job.  The hard refiners
+  require a start that covers `A`; the soft ones do not, since their objective constrains
+  nothing — build theirs with `feasible=:none`.
+- **Solvers** bundle the two.  [`symcover_min`](@ref), [`cover_min`](@ref),
+  [`soft_symcover`](@ref), [`soft_cover`](@ref), [`soft_symcover_min`](@ref), and
+  [`soft_cover_min`](@ref) refine every start on a menu (the `strategies` keyword, or the
+  multistart's own list) and return the best cover by [`cover_objective`](@ref), so their
+  result depends on `A` and not on an initialization the caller never chose.
+
+That is the rule for the whole grid: **the plain form owns the menu, so its result is a
+property of `A`; the `!` form refines the one start you give it, so its result is a
+property of `A` and that start.**  [`symcover!`](@ref) and [`cover!`](@ref) are the
+exception that proves it — they are initializers, not refiners, and construct their cover
+from scratch rather than reading the vector passed in.
 
 For finer control, you can run these manually:
 

src/ScaleInvariantAnalysis.jl

@@ -7,7 +7,8 @@ using Random: Random, AbstractRNG, MersenneTwister
 
 export AbsLog, AbsLinear
 export cover_objective
-export cover, cover!, symcover, symcover!, soft_symcover, soft_cover
+export cover, cover!, symcover, symcover!
+export soft_symcover, soft_symcover!, soft_cover, soft_cover!
 export initialize_cover, initialize_cover!, initialize_symcover, initialize_symcover!
 export symcover_min, symcover_min!, cover_min, cover_min!
 export soft_symcover_min, soft_symcover_min!, soft_cover_min, soft_cover_min!

src/soft_covers.jl

@@ -100,6 +100,58 @@ function soft_symcover(::AbsLinear{1}, A::AbstractMatrix; maxiter::Int=20, kwarg
     return a
 end
 
+"""
+    a = soft_symcover!(ϕ, a, A; maxiter=...)
+    a = soft_symcover!(a, A; maxiter=...)
+
+Refine the starting point `a` into a symmetric soft cover of `A`, in place, and return it.
+The no-ϕ form defaults to `AbsLinear{2}()`, matching [`soft_symcover`](@ref), whose
+supported ϕ values these methods share.
+
+This is the refiner half of [`soft_symcover`](@ref): the non-mutating form owns a
+multistart menu, so its result is a property of `A`, while this one descends from the
+single start you hand it, so its result is a property of `A` *and* that start. Building
+the start is the caller's job — see [`initialize_symcover`](@ref), and pass
+`feasible=:none`, since a soft cover is under no obligation to cover `A`.
+
+`a` must be finite and strictly positive on every row of `A` that carries support; scales
+on rows carrying no support are inert, and are zero on output. Unlike [`symcover_min!`](@ref),
+`a` need *not* cover `A` — the soft objective imposes no coverage constraint.
+
+`maxiter` bounds the descent sweeps; its default matches the corresponding
+[`soft_symcover`](@ref) method. Under `AbsLog{2}` the objective is convex with a unique
+minimizer, so the start is honored but not visible in the result.
+
+See also: [`soft_symcover`](@ref), [`soft_symcover_min!`](@ref), [`initialize_symcover`](@ref), [`soft_cover!`](@ref).
+"""
+function soft_symcover! end
+soft_symcover!(a::AbstractVector, A::AbstractMatrix; kwargs...) =
+    soft_symcover!(AbsLinear{2}(), a, A; kwargs...)
+
+function soft_symcover!(::AbsLog{2}, a::AbstractVector, A::AbstractMatrix; kwargs...)
+    _prepare_soft_symcover_start!(a, A, :soft_symcover!)
+    a .= soft_symcover_min(AbsLog{2}(), A; kwargs...)   # convex: the start is not read
+    return a
+end
+
+function soft_symcover!(::AbsLog{1}, a::AbstractVector, A::AbstractMatrix; maxiter::Int=20)
+    _prepare_soft_symcover_start!(a, A, :soft_symcover!)
+    _abslog1_iter!(a, A, maxiter)
+    return a
+end
+
+function soft_symcover!(::AbsLinear{2}, a::AbstractVector, A::AbstractMatrix; maxiter::Int=32)
+    _prepare_soft_symcover_start!(a, A, :soft_symcover!)
+    _abslinear2_iter!(a, A, maxiter)
+    return a
+end
+
+function soft_symcover!(::AbsLinear{1}, a::AbstractVector, A::AbstractMatrix; maxiter::Int=20)
+    _prepare_soft_symcover_start!(a, A, :soft_symcover!)
+    _abslinear1_iter!(a, A, maxiter)
+    return a
+end
+
 """
     a, b = soft_cover(ϕ, A; maxiter=200, starts=4, σ=2.0, rng=MersenneTwister(0))
     a, b = soft_cover(A; maxiter=200, starts=4, σ=2.0, rng=MersenneTwister(0))
@@ -193,6 +245,55 @@ function soft_cover(ϕ::AbsLinear{1}, A::AbstractMatrix; maxiter::Int=100, kwarg
     return _balance_cover!(a, b, A)
 end
 
+"""
+    a, b = soft_cover!(ϕ, a, b, A; maxiter=...)
+    a, b = soft_cover!(a, b, A; maxiter=...)
+
+Refine the starting point `(a, b)` into a soft cover of `A`, in place, and return it. This
+is the asymmetric counterpart of [`soft_symcover!`](@ref) and the refiner half of
+[`soft_cover`](@ref), carrying the same contract on the start: finite and strictly
+positive on every supported row and column, inert (and zero on output) elsewhere, and
+under no obligation to cover `A`. Build one with [`initialize_cover`](@ref) and
+`feasible=:none`. The no-ϕ form defaults to `AbsLinear{2}()`, matching
+[`soft_cover`](@ref).
+
+The product `a[i]*b[j]` is unchanged by `a -> c*a`, `b -> b/c`, so the start is read only
+up to that gauge: `(a, b)` and `(2a, b/2)` give the same result. The result itself is
+pinned to the balance convention of [`cover_min`](@ref), as every asymmetric cover in this
+package is.
+
+See also: [`soft_cover`](@ref), [`soft_cover_min!`](@ref), [`initialize_cover`](@ref), [`soft_symcover!`](@ref).
+"""
+function soft_cover! end
+soft_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) =
+    soft_cover!(AbsLinear{2}(), a, b, A; kwargs...)
+
+function soft_cover!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...)
+    _prepare_soft_cover_start!(a, b, A, :soft_cover!)
+    anew, bnew = soft_cover_min(AbsLog{2}(), A; kwargs...)   # convex: the start is not read
+    a .= anew
+    b .= bnew
+    return a, b
+end
+
+function soft_cover!(::AbsLog{1}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=20)
+    _prepare_soft_cover_start!(a, b, A, :soft_cover!)
+    _abslog1_iter_asym!(a, b, A, maxiter)
+    return _balance_cover!(a, b, A)
+end
+
+function soft_cover!(::AbsLinear{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=200)
+    _prepare_soft_cover_start!(a, b, A, :soft_cover!)
+    _mscm_als!(a, b, A, maxiter)
+    return _balance_cover!(a, b, A)
+end
+
+function soft_cover!(::AbsLinear{1}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; maxiter::Int=100)
+    _prepare_soft_cover_start!(a, b, A, :soft_cover!)
+    _abslinear1_iter_asym!(a, b, A, maxiter)
+    return _balance_cover!(a, b, A)
+end
+
 """
     a = soft_symcover_min(ϕ, A)
     a = soft_symcover_min(A)
@@ -279,12 +380,13 @@ function soft_symcover_min!(::AbsLog{2}, a::AbstractVector, A::AbstractMatrix; k
     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.
-function _prepare_soft_symcover_start!(a::AbstractVector, A::AbstractMatrix)
+# Shared prologue of the symmetric soft refiners (`soft_symcover!`, `soft_symcover_min!`).
+# The soft objective constrains nothing, so — unlike `_prepare_symcover_start!` — this
+# checks positivity only, and moves the start nowhere. `fname` names the caller so the
+# error reports the function the user actually called.
+function _prepare_soft_symcover_start!(a::AbstractVector, A::AbstractMatrix, fname::Symbol=:soft_symcover_min!)
     ax = axes(A, 1)
-    axes(A, 2) == ax || throw(ArgumentError("soft_symcover_min! requires a square matrix"))
+    axes(A, 2) == ax || throw(ArgumentError("$fname requires a square matrix"))
     eachindex(a) == ax || throw(DimensionMismatch("indices of `a` must match the indexing of `A`, got eachindex(a)=$(eachindex(a)), axes(A, 1)=$ax"))
     supp = fill!(similar(a, Bool), false)
     foreach_support_sym(A) do i, j, v
@@ -297,7 +399,7 @@ function _prepare_soft_symcover_start!(a::AbstractVector, A::AbstractMatrix)
     for i in ax
         supp[i] || continue
         (isfinite(a[i]) && a[i] > zero(a[i])) ||
-            throw(ArgumentError("soft_symcover_min! requires a start with finite positive scale on every supported row, got a[$i] = $(a[i])"))
+            throw(ArgumentError("$fname requires a start with finite positive scale on every supported row, got a[$i] = $(a[i])"))
     end
     return a
 end
@@ -391,10 +493,13 @@ function soft_cover_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::A
     return a, b
 end
 
-# Shared prologue of the `soft_cover_min!` kernels; the asymmetric counterpart of
-# `_prepare_soft_symcover_start!`. Positivity only — the soft objective constrains nothing —
-# plus the balance pin, so the refiners read the start only up to the row/column gauge.
-function _prepare_soft_cover_start!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix)
+# Shared prologue of the asymmetric soft refiners (`soft_cover!`, `soft_cover_min!`); the
+# counterpart of `_prepare_soft_symcover_start!`. Positivity only — the soft objective
+# constrains nothing — plus the balance pin, so the refiners read the start only up to the
+# row/column gauge. `fname` names the caller so the error reports the function the user
+# actually called.
+function _prepare_soft_cover_start!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix,
+                                    fname::Symbol=:soft_cover_min!)
     axes(A, 1) == eachindex(a) || throw(DimensionMismatch("indices of `a` must match row-indexing of `A`, got eachindex(a)=$(eachindex(a)), axes(A, 1)=$(axes(A, 1))"))
     axes(A, 2) == eachindex(b) || throw(DimensionMismatch("indices of `b` must match column-indexing of `A`, got eachindex(b)=$(eachindex(b)), axes(A, 2)=$(axes(A, 2))"))
     suppa = fill!(similar(a, Bool), false)
@@ -412,12 +517,12 @@ function _prepare_soft_cover_start!(a::AbstractVector, b::AbstractVector, A::Abs
     for i in eachindex(a)
         suppa[i] || continue
         (isfinite(a[i]) && a[i] > zero(a[i])) ||
-            throw(ArgumentError("soft_cover_min! requires a start with finite positive scale on every supported row, got a[$i] = $(a[i])"))
+            throw(ArgumentError("$fname requires a start with finite positive scale on every supported row, got a[$i] = $(a[i])"))
     end
     for j in eachindex(b)
         suppb[j] || continue
         (isfinite(b[j]) && b[j] > zero(b[j])) ||
-            throw(ArgumentError("soft_cover_min! requires a start with finite positive scale on every supported column, got b[$j] = $(b[j])"))
+            throw(ArgumentError("$fname requires a start with finite positive scale on every supported column, got b[$j] = $(b[j])"))
     end
     return _balance_cover!(a, b, A)
 end

test/soft_covers.jl

@@ -440,3 +440,89 @@ end
     @test soft_symcover_min(AbsLog{2}(), sym_zeros; linsolve=:lsqr) ≈
           soft_symcover_min(AbsLog{2}(), sym_zeros; linsolve=:dense) rtol=1e-6
 end
+
+@testset "soft_symcover!/soft_cover! refiners" begin
+    Asym = [4.0 1.0 0.5; 1.0 3.0 1.0; 0.5 1.0 2.5]
+    Agen = [1.0 2.0 0.5; 0.25 3.0 1.0]
+
+    @testset "refining the multistart's own start reproduces it: $ϕ" for ϕ in PENALTIES
+        a = initialize_symcover(Asym; strategy=:geomean, feasible=:none)
+        soft_symcover!(ϕ, a, Asym)
+        @test cover_objective(ϕ, a, Asym) ≈ cover_objective(ϕ, soft_symcover(ϕ, Asym), Asym) rtol=1e-6
+    end
+
+    @testset "no-ϕ form defaults to AbsLinear{2}" begin
+        a1, a2 = initialize_symcover(Asym; strategy=:geomean, feasible=:none), initialize_symcover(Asym; strategy=:geomean, feasible=:none)
+        @test soft_symcover!(a1, Asym) == soft_symcover!(AbsLinear{2}(), a2, Asym)
+        b1, c1 = initialize_cover(Agen; strategy=:geomean, feasible=:none)
+        b2, c2 = initialize_cover(Agen; strategy=:geomean, feasible=:none)
+        @test soft_cover!(b1, c1, Agen) == soft_cover!(AbsLinear{2}(), b2, c2, Agen)
+    end
+
+    # The refiner descends from the start it is handed; the multistart owns a menu.
+    # Under the convex AbsLog{2} the start is honored but cannot be seen in the result,
+    # while a non-convex AbsLinear{2} objective with two basins reports whichever the
+    # start lies in.
+    @testset "start-dependence" begin
+        Abasins = [0.021778451276962405 1.5690256886348526
+                   1.5690256886348526  0.20473123461805692]
+        starts() = (initialize_symcover(Abasins; strategy=:geomean, feasible=:none),
+                    initialize_symcover(Abasins; strategy=:hardcover))
+        o(ϕ, a) = cover_objective(ϕ, a, Abasins)
+
+        s1, s2 = starts()
+        @test o(AbsLog{2}(), soft_symcover!(AbsLog{2}(), s1, Abasins)) ≈
+              o(AbsLog{2}(), soft_symcover!(AbsLog{2}(), s2, Abasins))
+
+        s1, s2 = starts()
+        @test !isapprox(o(AbsLinear{2}(), soft_symcover!(AbsLinear{2}(), s1, Abasins)),
+                        o(AbsLinear{2}(), soft_symcover!(AbsLinear{2}(), s2, Abasins)); rtol=1e-6)
+    end
+
+    # Unlike symcover_min!, a soft refiner imposes no coverage constraint on its start.
+    @testset "start need not cover A: $ϕ" for ϕ in PENALTIES
+        a = fill(0.01, 3)
+        @test !iscover(a, Asym)
+        @test soft_symcover!(ϕ, a, Asym) === a
+        b, c = fill(0.01, 2), fill(0.01, 3)
+        @test !iscover(b, c, Agen)
+        @test soft_cover!(ϕ, b, c, Agen) === (b, c)
+    end
+
+    @testset "invalid starts throw, naming the function called" begin
+        @test_throws "soft_symcover! requires a start with finite positive scale" soft_symcover!([1.0, -1.0, 1.0], Asym)
+        @test_throws "soft_symcover! requires a start with finite positive scale" soft_symcover!([1.0, 0.0, 1.0], Asym)
+        @test_throws "soft_symcover! requires a start with finite positive scale" soft_symcover!([1.0, Inf, 1.0], Asym)
+        @test_throws "soft_symcover! requires a square matrix" soft_symcover!([1.0, 1.0], Agen)
+        @test_throws DimensionMismatch soft_symcover!([1.0, 1.0], Asym)
+        @test_throws "soft_cover! requires a start with finite positive scale" soft_cover!([1.0, -1.0], fill(1.0, 3), Agen)
+        @test_throws "soft_cover! requires a start with finite positive scale" soft_cover!(fill(1.0, 2), [1.0, -1.0, 1.0], Agen)
+        @test_throws DimensionMismatch soft_cover!(fill(1.0, 3), fill(1.0, 3), Agen)
+    end
+
+    # Scales on unsupported rows are inert and come back zero, matching every other
+    # cover in the package.
+    @testset "unsupported rows are zeroed: $ϕ" for ϕ in PENALTIES
+        Az = [1.0 0.0 2.0; 0.0 0.0 0.0; 2.0 0.0 3.0]
+        a = initialize_symcover(Az; strategy=:geomean, feasible=:none)
+        @test soft_symcover!(ϕ, a, Az)[2] == 0
+    end
+
+    @testset "asymmetric refiners pin the balance convention: $ϕ" for ϕ in PENALTIES
+        a, b = initialize_cover(Agen; strategy=:geomean, feasible=:none)
+        a .*= 7           # move the gauge; the objective cannot see it
+        b ./= 7
+        soft_cover!(ϕ, a, b, Agen)
+        @test isbalanced(a, b, Agen)
+    end
+
+    @testset "offset axes propagate: $ϕ" for ϕ in PENALTIES
+        Ao = OffsetArray(Asym, -1:1, -1:1)
+        ao = initialize_symcover(Ao; strategy=:geomean, feasible=:none)
+        soft_symcover!(ϕ, ao, Ao)
+        @test axes(ao, 1) == axes(Ao, 1)
+        aref = initialize_symcover(Asym; strategy=:geomean, feasible=:none)
+        soft_symcover!(ϕ, aref, Asym)
+        @test collect(ao) ≈ aref rtol=1e-6
+    end
+end