diff --git a/docs/src/index.md b/docs/src/index.md index fd70c3d..c43bc3f 100644 --- a/docs/src/index.md +++ b/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: diff --git a/src/ScaleInvariantAnalysis.jl b/src/ScaleInvariantAnalysis.jl index 0ab7c16..947b0a8 100644 --- a/src/ScaleInvariantAnalysis.jl +++ b/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! diff --git a/src/soft_covers.jl b/src/soft_covers.jl index 513e4e7..19b4540 100644 --- a/src/soft_covers.jl +++ b/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 diff --git a/test/soft_covers.jl b/test/soft_covers.jl index 1b65fc1..29cbaff 100644 --- a/test/soft_covers.jl +++ b/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