From a56247e4a50f52a7a48f2de67f1a98d33ea3af84 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 9 Oct 2026 04:22:19 -0500 Subject: [PATCH 1/7] Make AbsLog{2} stats describe the returned cover When the active-set finish certifies a minimizer, `stats.converged` was set from the finish but the only residual reported was the multiplier iteration's `kkt` trace, which can sit far above `vtol`. The finish now reports its own KKT check (`primal`, `ptol`, `dual`, `dtol`) in `stats.polish`. A finish that fails to certify returns the multiplier iterate, and now always warns; previously it warned only when that iterate's residual exceeded `vwarn`. With `maxouter = 0` (the soft AbsLog{2} covers) nothing is iterated, so `converged` is now `false` rather than `true`; no warning is issued. Assisted-by: Claude Opus 5.5 --- src/active_set.jl | 37 +++++++++++++++++++-------------- src/minimal_covers.jl | 47 +++++++++++++++++++++++++----------------- test/active_set.jl | 4 +++- test/minimal_covers.jl | 11 ++++++++++ 4 files changed, 64 insertions(+), 35 deletions(-) diff --git a/src/active_set.jl b/src/active_set.jl index e22f54e..e829ee9 100644 --- a/src/active_set.jl +++ b/src/active_set.jl @@ -530,8 +530,8 @@ function _most_violated(u::Vector{T}, inW::Matrix{Bool}, s::_DiffGridQP, de::T) end """ - u, certified, nsteps = _polish_difference_qp(supp, u0; τ, maxsteps, seed, budgets) - u, certified, nsteps = _polish_difference_qp(edges, cvals, u0; τ, maxsteps, seed, budgets) + u, certified, nsteps, kkt = _polish_difference_qp(supp, u0; τ, maxsteps, seed, budgets) + u, certified, nsteps, kkt = _polish_difference_qp(edges, cvals, u0; τ, maxsteps, seed, budgets) Minimize `½ ∑ₑ (u[p] - u[q] - c[e])²` subject to every term being nonnegative, starting from the near-optimal `u0`. The constraints are those of `supp`, an @@ -546,8 +546,11 @@ an augmented-Lagrangian iteration, the candidates are the constraints with positive weight or residual below `ptol`, largest weight first and then tightest first. `certified` reports that the returned `u` satisfies the KKT conditions to tolerances proportional to the largest magnitude in `c` and `u0`; -when it is false, `u` is `u0`. Each step solves a graph-Laplacian system, and -`maxsteps` bounds their number. `budgets`, a `NamedTuple` with fields +when it is false, `u` is `u0`. `kkt = (; primal, ptol, dual, dtol)` gives the +largest constraint violation of the returned `u` and the most negative active +multiplier (as nonnegative magnitudes), with the tolerances they were certified +against; `primal` and `dual` are `NaN` when `certified` is false. Each step +solves a graph-Laplacian system, and `maxsteps` bounds their number. `budgets`, a `NamedTuple` with fields `fillbudget` and/or `flopbudget`, sets the limits above which these systems are solved by conjugate gradients rather than a sparse Cholesky factorization. The constant on each connected component of the graph is taken from `u0`. @@ -599,6 +602,7 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ u = u0 u, nsteps = _dual_feasible!(inW, W, λ, s, nV, dtol, nsteps, maxsteps; budgets) certified = false + nokkt = (; primal=T(NaN), ptol, dual=T(NaN), dtol) while nsteps <= maxsteps e, de = _most_violated(u, inW, s, -ptol) if e == 0 @@ -619,7 +623,7 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ sv[q] += one(T) _route!(r, sv, F, s) t, wb = _ratio_test(λ, r, W) - wb == 0 && return u0, false, nsteps # the constraints are infeasible + wb == 0 && return u0, false, nsteps, nokkt # the constraints are infeasible for w in W λ[w] += t * r[w] end @@ -634,7 +638,7 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ y = _tree_solve(K, rhs) z = [y[F.tree[v]] for v in 1:nV] a = z[p] - z[q] - a > 0 || return u0, false, nsteps + a > 0 || return u0, false, nsteps, nokkt # Rates of the multipliers in `W`: route `L z - gₑ` onto the forest. sv = _difference_sums!(zeros(T, nV), z, s, false) sv[p] -= one(T) @@ -661,10 +665,12 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ # Recompute from `W` alone, so rounding does not accumulate. u, nsteps = _dual_feasible!(inW, W, λ, s, nV, dtol, nsteps, maxsteps; budgets) end - certified || return u0, false, nsteps + certified || return u0, false, nsteps, nokkt # Restore each component's constant from `u0`. (`u` is reassigned above, so # capturing it in a closure would box it; `ustar` is bound once.) ustar = u + kkt = (; primal=max(zero(T), -minimum(e -> _resid(ustar, s, e), 1:_nconstraints(s))), ptol, + dual=max(zero(T), -minimum(k -> λ[k], W; init=zero(T))), dtol) uf .= 1:nV _foreach_constraint(s) do p, q rp, rq = find(p), find(q) @@ -676,7 +682,7 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ rv = find(v) rv == v && (shift[v] = u0[v] - ustar[v]) end - return [ustar[v] + shift[find(v)] for v in 1:nV], true, nsteps + return [ustar[v] + shift[find(v)] for v in 1:nV], true, nsteps, kkt end function _deactivate!(inW, W::Vector{Int}, λ, w::Int) @@ -722,7 +728,7 @@ function _ratio_test(λ::Vector{T}, r::Vector{T}, W::Vector{Int}) where {T} end # Solve the `AbsLog{2}` problem of `_cover_min_abslog2` exactly from its -# augmented-Lagrangian iterate `x`, returning `(x, certified, nsteps)`. In the +# augmented-Lagrangian iterate `x`, returning `(x, certified, nsteps, kkt)`. In the # stacked layout `(α; β)` an entry's residual is `α[i] + β[j] - c`, which is a # difference in `(α; -β)`; the transversal layout is a difference already. # `multipliers`, the iteration's multipliers in the layout of `supp`, seed the @@ -731,9 +737,9 @@ function _polish_cover(x::Vector{T}, supp::EdgeList{T}, m::Int; multipliers=noth budgets::NamedTuple=(;)) where {T} supp.qsign < 0 && return _polish_difference_qp(supp, x; seed=multipliers, budgets) flip(u) = [k <= m ? u[k] : -u[k] for k in eachindex(u)] - u, certified, nsteps = _polish_difference_qp(EdgeList{T}(supp.edges, supp.cvals, -oneunit(T)), flip(x); - seed=multipliers, budgets) - return flip(u), certified, nsteps + u, certified, nsteps, kkt = _polish_difference_qp(EdgeList{T}(supp.edges, supp.cvals, -oneunit(T)), flip(x); + seed=multipliers, budgets) + return flip(u), certified, nsteps, kkt end function _polish_cover(x::Vector{T}, supp::Grid{T}, m::Int; multipliers=nothing, @@ -755,7 +761,8 @@ _polish_cover(x::Vector{T}, supp::DiffGrid{T}, m::Int; multipliers=nothing, budg _polish_difference_qp(supp, x; seed=multipliers === nothing ? nothing : vec(multipliers), budgets) # Solve the `AbsLog{2}` problem of `_symcover_min_abslog2` exactly from its -# augmented-Lagrangian iterate `x`, returning `(x, certified, nsteps)`. The +# augmented-Lagrangian iterate `x`, returning `(x, certified, nsteps, kkt)`, with +# `kkt` describing the double-cover solution. The # symmetric problem is the asymmetric one on the bipartite double cover: with # unknowns `u = (α; β)` and start `(x; -x)`, each stored entry `(p, q)` gives the # difference constraints `α[p] - β[q] ≥ c` and, off the diagonal, `α[q] - β[p] ≥ c`. @@ -808,8 +815,8 @@ function _polish_symcover(x::Vector{T}, supp::Grid{T}; multipliers=nothing, end # Symmetric scales `(α - β)/2` from the double-cover solution `u = (α; β)`. -function _unfold_symcover(x::Vector{T}, (u, certified, nsteps)) where {T} +function _unfold_symcover(x::Vector{T}, (u, certified, nsteps, kkt)) where {T} n = length(x) xnew = (u[1:n] .- u[n+1:2n]) ./ 2 - return xnew, certified, nsteps + return xnew, certified, nsteps, kkt end diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 4fe35b4..a5489b0 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -1171,13 +1171,16 @@ end const AL_PENALTY = 1e2 const AL_MAXOUTER = 32 -# Warn when the outer iteration ends with a KKT residual well above the -# inner solver's accuracy floor and no active-set finish certified the result. +# Warn when the returned scales may not minimize the objective: an active-set +# finish ran and did not certify a minimizer, or none ran and the outer +# iteration ended with a KKT residual well above the inner solver's accuracy +# floor. With `maxouter = 0` there is no iteration and no warning. function _warn_unconverged(fname::Symbol, stats, maxouter::Int) - stats.converged && return nothing + (stats.converged || maxouter == 0) && return nothing v = isempty(stats.kkt) ? oftype(stats.vwarn, NaN) : stats.kkt[end] - v > stats.vwarn || return nothing - advice = haskey(stats, :polish) && !stats.polish.certified ? + uncertified = haskey(stats, :polish) && !stats.polish.certified + uncertified || v > stats.vwarn || return nothing + advice = uncertified ? "The active-set finish that follows the multiplier iteration did not certify a minimizer, so the result is the multiplier iterate." : stats.nouter < maxouter ? "The residual stopped contracting before the update limit, so a larger `maxouter` will not help; try a larger `κ`." : "Increase `maxouter` or `κ`." @@ -1343,11 +1346,12 @@ _dense_factor_type(::Type{T}) where {T} = LinearAlgebra.LU{T,Matrix{T},Vector{In # log scales and a statistics tuple: the final multipliers (`multipliers`, laid # out by `_multiplier_storage`), inner solve and iteration counts, the # LSQR iterations of each solve (`lsqrtrace`), the per-update exits, drops, KKT -# residuals (`kkt`) and penalty weights (`κs`), the convergence flag with its -# tolerances, the `linsolve` and `precond` choices, the full Cholesky -# preconditioner's predicted entry count `fill_entries` and flop count -# `factor_flops` (`0` and `0.0` without a factor), the number of numeric -# factorizations `nrefactor` (full or heavy-forest), the numbers of LSQR solves +# residuals (`kkt`) and penalty weights (`κs`), the convergence flag +# (`kkt[end] <= vtol`; always `false` for `maxouter = 0`) with its tolerances, +# the `linsolve` and `precond` choices, the full Cholesky preconditioner's +# predicted entry count `fill_entries` and flop count `factor_flops` (`0` and +# `0.0` without a factor), the number of numeric factorizations `nrefactor` +# (full or heavy-forest), the numbers of LSQR solves # preconditioned by the heavy-forest factor (`nforest`) and by the diagonal # (`ndiagonal`) when the full factor is over budget, and the number `nzeroed` # of positive multipliers zeroed directly on entries with slack while the @@ -1835,7 +1839,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; drops = T[] viols = T[] κtrace = T[] - converged = maxouter == 0 + converged = false vprev = T(Inf) nstall = 0 vbest = T(Inf) # residual that the next tenfold contraction is measured from @@ -1938,7 +1942,11 @@ end # With `polish`, the active-set solver of `_polish_symcover` finishes a # multiplier iteration that stops before converging, and the iteration then # stops at its first stall at the penalty cap rather than zeroing slack -# multipliers (`zeroslack`, which defaults to `!polish`). +# multipliers (`zeroslack`, which defaults to `!polish`). When the finish runs, +# `stats.polish` holds `certified`, the step count `nsteps`, and the KKT check of +# `_polish_difference_qp` (`primal`, `ptol`, `dual`, `dtol`), and +# `stats.converged == stats.polish.certified`. A certified finish replaces the +# multiplier iterate, so `stats.kkt` then describes the iterate, not the result. function _symcover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::Int=AL_MAXOUTER, maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, boost::Bool=true, fillbudget::Real=LSQR_FILL_BUDGET, @@ -2045,9 +2053,9 @@ function _symcover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter: T[hassupp[ip] ? log(T(start[i])) : zero(T) for (ip, i) in enumerate(ax)] α, stats = _abslog2_auglag(sys, x0; κ, maxouter, maxiter, linsolve, boost, fillbudget, flopbudget, zeroslack) if polish && !stats.converged && maxouter > 0 - α, certified, nsteps = _polish_symcover(α, supp; multipliers=stats.multipliers, - budgets=(; fillbudget, flopbudget)) - stats = merge(stats, (; converged=certified, polish=(; certified, nsteps))) + α, certified, nsteps, kkt = _polish_symcover(α, supp; multipliers=stats.multipliers, + budgets=(; fillbudget, flopbudget)) + stats = merge(stats, (; converged=certified, polish=(; certified, nsteps, kkt...))) end _warn_unconverged(fname, stats, maxouter) _check_representable(α, ip -> hassupp[ip], nothing, nothing, T, fname; gauge=false) @@ -2068,7 +2076,8 @@ end # transversal. With `polish`, the active-set solver of `_polish_cover` finishes # a multiplier iteration that stops before converging, and the iteration then # stops at its first stall at the penalty cap rather than zeroing slack -# multipliers (`zeroslack`, which defaults to `!polish`). +# multipliers (`zeroslack`, which defaults to `!polish`). `stats` is as for +# `_symcover_min_abslog2`. function _cover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::Int=AL_MAXOUTER, maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, boost::Bool=true, fillbudget::Real=LSQR_FILL_BUDGET, @@ -2304,9 +2313,9 @@ function _cover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::In x, stats = _abslog2_auglag(sys, x0; κ, maxouter, maxiter, linsolve, boost=boost && !tt, fillbudget, flopbudget, zeroslack) if polish && !stats.converged && maxouter > 0 - x, certified, nsteps = _polish_cover(x, supp, m; multipliers=stats.multipliers, - budgets=(; fillbudget, flopbudget)) - stats = merge(stats, (; converged=certified, polish=(; certified, nsteps))) + x, certified, nsteps, kkt = _polish_cover(x, supp, m; multipliers=stats.multipliers, + budgets=(; fillbudget, flopbudget)) + stats = merge(stats, (; converged=certified, polish=(; certified, nsteps, kkt...))) end _warn_unconverged(fname, stats, maxouter) if tt diff --git a/test/active_set.jl b/test/active_set.jl index 7426f6f..708de4d 100644 --- a/test/active_set.jl +++ b/test/active_set.jl @@ -37,6 +37,7 @@ using MatrixCovers: _cover_min_abslog2, _polish_difference_qp @test sr.converged a, b, s = _cover_min_abslog2(A; maxouter=2) @test s.polish.certified + @test s.converged && s.polish.primal <= s.polish.ptol && s.polish.dual <= s.polish.dtol @test logproduct_gap(a, b, ar, br, A) <= 1e-8 at, bt, _ = _cover_min_abslog2(A; transversal=true, fname=:cover_transversal) a, b, s = _cover_min_abslog2(A; transversal=true, maxouter=1, fname=:cover_transversal) @@ -84,9 +85,10 @@ end @test_throws "seed must have axes 1:3" _polish_difference_qp([(1, 2), (2, 3), (1, 3)], [0.0, 0.0, 0.0], zeros(3); seed=ones(2)) # Infeasible constraints (a cycle with positive total cost) are not certified. edges = [(1, 2), (2, 3), (3, 1)] - u, certified, _ = _polish_difference_qp(edges, [1.0, 1.0, 1.0], zeros(3)) + u, certified, _, kkt = _polish_difference_qp(edges, [1.0, 1.0, 1.0], zeros(3)) @test !certified @test u == zeros(3) + @test isnan(kkt.primal) && isnan(kkt.dual) end @testset "tree Laplacian assembly" begin diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index b8b6c4e..aa9f817 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -85,6 +85,17 @@ end stats = (; converged=false, kkt=(1.0,), vtol=1e-13, vwarn=1e-10, nouter=3, polish=(; certified=false, nsteps=5)) @test_logs (:warn, r"active-set finish .* did not certify") MatrixCovers._warn_unconverged(:symcover_min, stats, 32) + # An uncertified finish warns even when the residual is below `vwarn`. + stats = merge(stats, (; kkt=(1e-12,))) + @test_logs (:warn, r"active-set finish .* did not certify") MatrixCovers._warn_unconverged(:symcover_min, stats, 32) + # A certified finish reports a KKT check that meets its tolerances. + _, s = MatrixCovers._symcover_min_abslog2(A; maxouter=3) + @test s.converged && s.polish.certified + @test s.polish.primal <= s.polish.ptol && s.polish.dual <= s.polish.dtol + # With no outer iteration nothing is checked, so the result is not + # reported as converged, and there is no warning. + _, s0 = @test_logs MatrixCovers._symcover_min_abslog2(A; maxouter=0, boost=false, fname=:soft_symcover) + @test !s0.converged && s0.nouter == 0 && !haskey(s0, :polish) end @testset "cover_min native AbsLog{2}" begin From 99cb3b15c58c6e65c647e74a9014ecc1ff80505a Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 9 Oct 2026 05:07:27 -0500 Subject: [PATCH 2/7] Return solver statistics from the AbsLog{2} covers `symcover_min`, `cover_min`, `cover_transversal`, their mutating forms, and the soft `AbsLog{2}` covers now return a statistics object as their last output: `a, stats = symcover_min(A)`, `a, b, stats = cover_min(A)`. The new public (unexported) `AbstractCoverStats` has one subtype per solver family; this commit adds `AugmentedLagrangianStats` (with `ActiveSetStats` for the finish) and `LeastSquaresStats`. The accessors `converged`, `iterations`, `residual`, and `tolerance` describe the returned cover; `LeastSquaresStats` implements none of them, since the soft cover has no stopping criterion. The augmented Lagrangian gains a `tol` keyword (default `1000eps`), and any result that does not meet it warns; the separate warning threshold is gone. The public methods declare their keywords instead of forwarding them to the workers, so internal options are no longer reachable, and the sparse methods delegate to the generic ones. Assisted-by: Claude Opus 5.5 --- ext/MatrixCoversUnitfulExt.jl | 31 ++++-- src/MatrixCovers.jl | 3 +- src/minimal_covers.jl | 80 ++++++++------ src/soft_covers.jl | 23 ++-- src/sparse_support.jl | 67 ++++-------- src/stats.jl | 191 ++++++++++++++++++++++++++++++++++ src/transversal.jl | 14 ++- test/element_types.jl | 28 ++--- test/extensions.jl | 118 ++++++++++----------- test/helpers.jl | 7 ++ test/heuristic_covers.jl | 4 +- test/invariants.jl | 14 +-- test/minimal_covers.jl | 141 ++++++++++++++----------- test/powermean.jl | 8 +- test/soft_covers.jl | 72 ++++++------- test/storage_types.jl | 20 ++-- test/support.jl | 2 +- test/transversal.jl | 10 +- test/unitful.jl | 48 ++++----- 19 files changed, 550 insertions(+), 331 deletions(-) create mode 100644 src/stats.jl diff --git a/ext/MatrixCoversUnitfulExt.jl b/ext/MatrixCoversUnitfulExt.jl index e4b6568..6653fbf 100644 --- a/ext/MatrixCoversUnitfulExt.jl +++ b/ext/MatrixCoversUnitfulExt.jl @@ -191,47 +191,56 @@ reattach(a, ua) = a .* ua # Entry points # ============================================================ -sym(f, A::QMatrix, ϕ...; kwargs...) = reattach(f(ϕ..., strip_matrix(A); kwargs...), factor_units_sym(A)) +# The optimizing solvers return statistics after the scale vectors; split them off +# so they pass through unchanged. +symsplit(x::AbstractVector) = (x, ()) +symsplit(r::Tuple) = (first(r), Base.tail(r)) +asymsplit(r::Tuple) = (r[1], r[2], Base.tail(Base.tail(r))) + +function sym(f, A::QMatrix, ϕ...; kwargs...) + a, rest = symsplit(f(ϕ..., strip_matrix(A); kwargs...)) + return isempty(rest) ? reattach(a, factor_units_sym(A)) : (reattach(a, factor_units_sym(A)), rest...) +end function asym(f, A::QMatrix, ϕ...; kwargs...) ua, ub = factor_units(A) - a, b = f(ϕ..., strip_matrix(A); kwargs...) - return reattach(a, ua), reattach(b, ub) + a, b, rest = asymsplit(f(ϕ..., strip_matrix(A); kwargs...)) + return (reattach(a, ua), reattach(b, ub), rest...) end # Allocating scratch from `A` permits uninitialized destination vectors. function sym!(f, a::QVector, A::QMatrix, ϕ...; kwargs...) ua = factor_units_sym(A) An = strip_matrix(A) - x = f(ϕ..., similar(a, float(real(eltype(An)))), An; kwargs...) + x, rest = symsplit(f(ϕ..., similar(a, float(real(eltype(An)))), An; kwargs...)) a .= reattach(x, ua) - return a + return isempty(rest) ? a : (a, rest...) end function asym!(f, a::QVector, b::QVector, A::QMatrix, ϕ...; kwargs...) ua, ub = factor_units(A) An = strip_matrix(A) T = float(real(eltype(An))) - x, y = f(ϕ..., similar(a, T), similar(b, T), An; kwargs...) + x, y, rest = asymsplit(f(ϕ..., similar(a, T), similar(b, T), An; kwargs...)) a .= reattach(x, ua) b .= reattach(y, ub) - return a, b + return (a, b, rest...) end # `a` is a start, and is read. function symstart!(f, a::QVector, A::QMatrix, ϕ...; kwargs...) ua = factor_units_sym(A) - x = f(ϕ..., strip_start(a, ua), strip_matrix(A); kwargs...) + x, rest = symsplit(f(ϕ..., strip_start(a, ua), strip_matrix(A); kwargs...)) a .= reattach(x, ua) - return a + return isempty(rest) ? a : (a, rest...) end function asymstart!(f, a::QVector, b::QVector, A::QMatrix, ϕ...; kwargs...) ua, ub = factor_units(A) - x, y = f(ϕ..., strip_start(a, ua), strip_start(b, ub), strip_matrix(A); kwargs...) + x, y, rest = asymsplit(f(ϕ..., strip_start(a, ua), strip_start(b, ub), strip_matrix(A); kwargs...)) a .= reattach(x, ua) b .= reattach(y, ub) - return a, b + return (a, b, rest...) end # Mirror the core penalty dispatch while specializing on unitful matrices. diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index 5641ade..a07dec4 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -18,10 +18,11 @@ export cover_transversal, cover_transversal! # `public` is parsed as a keyword only from Julia 1.11; this package supports 1.10. @static if VERSION >= v"1.11" - eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym, scalar_type, SupportComponents, support_components, ncomponents, rowcomponent, colcomponent, SolverFailure")) + eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym, scalar_type, SupportComponents, support_components, ncomponents, rowcomponent, colcomponent, SolverFailure, AbstractCoverStats, AugmentedLagrangianStats, ActiveSetStats, LeastSquaresStats, converged, iterations, residual, tolerance")) end include("penalties.jl") +include("stats.jl") # statistics returned by the optimizing solvers include("support.jl") include("iscover.jl") include("fastlog.jl") diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index a5489b0..39faffe 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -172,33 +172,45 @@ cover_min!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) = # # Each `κ` stage freezes the weights, solves the normal equations, and uses a # backtracking line search. A final uniform shift restores feasibility. -function symcover_min(::AbsLog{2}, A::AbstractMatrix; kwargs...) - a, _ = _symcover_min_abslog2(A; kwargs...) - return a +function symcover_min(::AbsLog{2}, A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::Int=AL_MAXOUTER, + maxiter::Int=AL_MAXITER, linsolve::Symbol=:auto, + fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET, + tol::Union{Real,Nothing}=nothing) + a, nt = _symcover_min_abslog2(A; κ, maxouter, maxiter, linsolve, fillbudget, flopbudget, tol) + return a, AugmentedLagrangianStats(A, nt) end # Asymmetric counterpart on stacked log scales `(α; β)`. The solve pins the # global gauge; `_cover_min_abslog2` regularizes and balances the remaining # component gauges. -function cover_min(::AbsLog{2}, A::AbstractMatrix; kwargs...) - a, b, _ = _cover_min_abslog2(A; kwargs...) - return a, b +function cover_min(::AbsLog{2}, A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::Int=AL_MAXOUTER, + maxiter::Int=AL_MAXITER, linsolve::Symbol=:auto, + fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET, + tol::Union{Real,Nothing}=nothing) + a, b, nt = _cover_min_abslog2(A; κ, maxouter, maxiter, linsolve, fillbudget, flopbudget, tol) + return a, b, AugmentedLagrangianStats(A, nt) end # The convex objective has the same result from every valid start. -function symcover_min!(::AbsLog{2}, a::AbstractVector, A::AbstractMatrix; kwargs...) +function symcover_min!(::AbsLog{2}, a::AbstractVector, A::AbstractMatrix; κ::Real=AL_PENALTY, + maxouter::Int=AL_MAXOUTER, maxiter::Int=AL_MAXITER, linsolve::Symbol=:auto, + fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET, + tol::Union{Real,Nothing}=nothing) _prepare_symcover_start!(a, A) - anew, _ = _symcover_min_abslog2(A; start=a, kwargs...) + anew, nt = _symcover_min_abslog2(A; start=a, κ, maxouter, maxiter, linsolve, fillbudget, flopbudget, tol) a .= anew - return a + return a, AugmentedLagrangianStats(A, nt) end -function cover_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) +function cover_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; κ::Real=AL_PENALTY, + maxouter::Int=AL_MAXOUTER, maxiter::Int=AL_MAXITER, linsolve::Symbol=:auto, + fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET, + tol::Union{Real,Nothing}=nothing) _prepare_cover_start!(a, b, A) - anew, bnew, _ = _cover_min_abslog2(A; start=(a, b), kwargs...) + anew, bnew, nt = _cover_min_abslog2(A; start=(a, b), κ, maxouter, maxiter, linsolve, fillbudget, flopbudget, tol) a .= anew b .= bnew - return a, b + return a, b, AugmentedLagrangianStats(A, nt) end # Multistart drivers for the `AbsLinear` kernels in MatrixCoversIpoptExt. @@ -1170,21 +1182,19 @@ end # Augmented-Lagrangian defaults. const AL_PENALTY = 1e2 const AL_MAXOUTER = 32 +const AL_MAXITER = 40 -# Warn when the returned scales may not minimize the objective: an active-set -# finish ran and did not certify a minimizer, or none ran and the outer -# iteration ended with a KKT residual well above the inner solver's accuracy -# floor. With `maxouter = 0` there is no iteration and no warning. +# Warn when the returned scales may not minimize the objective. With +# `maxouter = 0` there is no iteration and no warning. function _warn_unconverged(fname::Symbol, stats, maxouter::Int) (stats.converged || maxouter == 0) && return nothing - v = isempty(stats.kkt) ? oftype(stats.vwarn, NaN) : stats.kkt[end] - uncertified = haskey(stats, :polish) && !stats.polish.certified - uncertified || v > stats.vwarn || return nothing + v = isempty(stats.kkt) ? oftype(stats.tol, NaN) : stats.kkt[end] + uncertified = stats.polish !== nothing && !stats.polish.certified advice = uncertified ? "The active-set finish that follows the multiplier iteration did not certify a minimizer, so the result is the multiplier iterate." : stats.nouter < maxouter ? "The residual stopped contracting before the update limit, so a larger `maxouter` will not help; try a larger `κ`." : "Increase `maxouter` or `κ`." - @warn "$fname: the multiplier iteration ended after $(stats.nouter) of maxouter=$maxouter updates with KKT residual $v (tolerance $(stats.vtol)); the result covers `A` but may not minimize the objective. $advice" + @warn "$fname: the multiplier iteration ended after $(stats.nouter) of maxouter=$maxouter updates with KKT residual $v (tolerance $(stats.tol)); the result covers `A` but may not minimize the objective. $advice" return nothing end @@ -1347,7 +1357,7 @@ _dense_factor_type(::Type{T}) where {T} = LinearAlgebra.LU{T,Matrix{T},Vector{In # out by `_multiplier_storage`), inner solve and iteration counts, the # LSQR iterations of each solve (`lsqrtrace`), the per-update exits, drops, KKT # residuals (`kkt`) and penalty weights (`κs`), the convergence flag -# (`kkt[end] <= vtol`; always `false` for `maxouter = 0`) with its tolerances, +# (`kkt[end] <= tol`; always `false` for `maxouter = 0`) with its tolerance, # the `linsolve` and `precond` choices, the full Cholesky preconditioner's # predicted entry count `fill_entries` and flop count `factor_flops` (`0` and # `0.0` without a factor), the number of numeric factorizations `nrefactor` @@ -1363,8 +1373,11 @@ _dense_factor_type(::Type{T}) where {T} = LinearAlgebra.LU{T,Matrix{T},Vector{In function _abslog2_auglag(sys::SupportSystem{T}, x0; κ::Real, maxouter::Int, maxiter::Int, linsolve::Symbol, boost::Bool, fillbudget::Real=LSQR_FILL_BUDGET, - flopbudget::Real=LSQR_FLOP_BUDGET, zeroslack::Bool=true) where {T} + flopbudget::Real=LSQR_FLOP_BUDGET, zeroslack::Bool=true, + tol::Union{Real,Nothing}=nothing) where {T} κ > 1 || throw(ArgumentError("κ must exceed 1; got $κ")) + tol === nothing || (isfinite(tol) && tol >= 0) || + throw(ArgumentError("tol must be finite and nonnegative; got $tol")) maxouter >= 0 || throw(ArgumentError("maxouter must be nonnegative; got $maxouter")) N = sys.N supp = sys.supp @@ -1826,8 +1839,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; x = x0 === nothing ? solve_weighted(zeros(T, N), nothing) : x0 # Absolute thresholds preserve covariance under rescaling because the # residuals are logarithmic. - vtol = 1000 * eps(T) - vwarn = max(sqrt(eps(T)), T(1e-6)) + vtol = tol === nothing ? 1000 * eps(T) : T(tol) # Precision bounds the usable penalty: the cap is `1e5` (LSQR) or `1e8` # (exact solves) in `Float64` and grows as `1/sqrt(eps(T))`. κcap = T(use_lsqr ? 1e5 : 1e8) * sqrt(T(eps(Float64)) / eps(T)) @@ -1931,7 +1943,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; return x, (; multipliers=λ, nsolves=nsolves[], lsqriters=nlsqr[], lsqrtrace=Tuple(lsqrtrace), cgiters=ncg[], cholsolves=nchol[], nouter=length(exits), exits=Tuple(exits), drops=Tuple(drops), kkt=Tuple(viols), κs=Tuple(κtrace), - converged, vtol, vwarn, + converged, tol=vtol, polish=nothing, linsolve=(use_lsqr ? :lsqr : use_cg ? :woodbury : :dense), precond=(!use_precond ? :none : use_factor ? :factor : nforest[] > 0 ? :forest : :diagonal), nrefactor=nrefactor[], nforest=nforest[], ndiagonal=ndiagonal[], fill_entries, factor_flops=flops, nzeroed) @@ -1948,10 +1960,11 @@ end # `stats.converged == stats.polish.certified`. A certified finish replaces the # multiplier iterate, so `stats.kkt` then describes the iterate, not the result. function _symcover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::Int=AL_MAXOUTER, - maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, + maxiter::Int=AL_MAXITER, linsolve::Symbol=:auto, start=nothing, boost::Bool=true, fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET, polish::Bool=true, - zeroslack::Bool=!polish, fname=:symcover_min) + zeroslack::Bool=!polish, tol::Union{Real,Nothing}=nothing, + fname=:symcover_min) linsolve in (:auto, :dense, :lsqr, :woodbury) || throw(ArgumentError("linsolve must be :auto, :dense, :lsqr, or :woodbury; got :$linsolve")) # Shared symmetry check for native symmetric minimal covers. @@ -1964,7 +1977,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter: if eps(T) > eps(Float64) a64, stats = _symcover_min_abslog2(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A); κ, maxouter, maxiter, linsolve, start, boost, fillbudget, flopbudget, - polish, zeroslack, fname) + polish, zeroslack, tol, fname) # Narrowing rounds to nearest and so can round a product below its entry. a = T.(a64) boost && _certify_cover!(a, A, fname) @@ -2051,7 +2064,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter: zeros(T, n)) x0 = start === nothing ? nothing : T[hassupp[ip] ? log(T(start[i])) : zero(T) for (ip, i) in enumerate(ax)] - α, stats = _abslog2_auglag(sys, x0; κ, maxouter, maxiter, linsolve, boost, fillbudget, flopbudget, zeroslack) + α, stats = _abslog2_auglag(sys, x0; κ, maxouter, maxiter, linsolve, boost, fillbudget, flopbudget, zeroslack, tol) if polish && !stats.converged && maxouter > 0 α, certified, nsteps, kkt = _polish_symcover(α, supp; multipliers=stats.multipliers, budgets=(; fillbudget, flopbudget)) @@ -2079,11 +2092,12 @@ end # multipliers (`zeroslack`, which defaults to `!polish`). `stats` is as for # `_symcover_min_abslog2`. function _cover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::Int=AL_MAXOUTER, - maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, + maxiter::Int=AL_MAXITER, linsolve::Symbol=:auto, start=nothing, boost::Bool=true, fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET, transversal::Union{Nothing,Bool,AbstractVector{<:Integer}}=nothing, - polish::Bool=true, zeroslack::Bool=!polish, fname::Symbol=:cover_min) + polish::Bool=true, zeroslack::Bool=!polish, tol::Union{Real,Nothing}=nothing, + fname::Symbol=:cover_min) tt = transversal !== nothing && transversal !== false if tt linsolve in (:auto, :dense, :lsqr, :woodbury) || @@ -2101,7 +2115,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::In if eps(T) > eps(Float64) a64, b64, stats = _cover_min_abslog2(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A); κ, maxouter, maxiter, linsolve, start, boost, fillbudget, flopbudget, - transversal, polish, zeroslack, fname) + transversal, polish, zeroslack, tol, fname) # Narrowing rounds to nearest and so can round a product below its entry. a, b = T.(a64), T.(b64) boost && _certify_cover!(a, b, A, fname) @@ -2311,7 +2325,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::In s0 end x, stats = _abslog2_auglag(sys, x0; κ, maxouter, maxiter, linsolve, - boost=boost && !tt, fillbudget, flopbudget, zeroslack) + boost=boost && !tt, fillbudget, flopbudget, zeroslack, tol) if polish && !stats.converged && maxouter > 0 x, certified, nsteps, kkt = _polish_cover(x, supp, m; multipliers=stats.multipliers, budgets=(; fillbudget, flopbudget)) diff --git a/src/soft_covers.jl b/src/soft_covers.jl index d5793a5..3d332a3 100644 --- a/src/soft_covers.jl +++ b/src/soft_covers.jl @@ -88,11 +88,12 @@ julia> round.(soft_symcover([0 1; 1 0]); digits=4) """ soft_symcover(A::AbstractMatrix; kwargs...) = soft_symcover(PowerMean{2}(), A; kwargs...) -function soft_symcover(::AbsLog{2}, A::AbstractMatrix; kwargs...) +function soft_symcover(::AbsLog{2}, A::AbstractMatrix; linsolve::Symbol=:auto, + fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET) ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("soft_symcover requires a square matrix")) - a, _ = _soft_symcover_abslog2(A; kwargs...) - return a + a, nt = _soft_symcover_abslog2(A; linsolve, fillbudget, flopbudget) + return a, LeastSquaresStats(nt) end # Sole owner of the starts/σ/rng defaults for the symmetric `AbsLinear` multistart. @@ -129,8 +130,9 @@ soft_symcover!(a::AbstractVector, A::AbstractMatrix; kwargs...) = function soft_symcover!(::AbsLog{2}, a::AbstractVector, A::AbstractMatrix; kwargs...) _prepare_soft_symcover_start!(a, A) - a .= soft_symcover(AbsLog{2}(), A; kwargs...) # convex: the start is not read - return a + anew, stats = soft_symcover(AbsLog{2}(), A; kwargs...) # convex: the start is not read + a .= anew + return a, stats end """ @@ -196,9 +198,10 @@ julia> a * b' """ soft_cover(A::AbstractMatrix; kwargs...) = soft_cover(PowerMean{2}(), A; kwargs...) -function soft_cover(::AbsLog{2}, A::AbstractMatrix; kwargs...) - a, b, _ = _soft_cover_abslog2(A; kwargs...) - return a, b +function soft_cover(::AbsLog{2}, A::AbstractMatrix; linsolve::Symbol=:auto, + fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET) + a, b, nt = _soft_cover_abslog2(A; linsolve, fillbudget, flopbudget) + return a, b, LeastSquaresStats(nt) end # Sole owner of the starts/σ/rng defaults for the asymmetric `AbsLinear` multistart. @@ -231,10 +234,10 @@ soft_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) function soft_cover!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) _prepare_soft_cover_start!(a, b, A) - anew, bnew = soft_cover(AbsLog{2}(), A; kwargs...) # convex: the start is not read + anew, bnew, stats = soft_cover(AbsLog{2}(), A; kwargs...) # convex: the start is not read a .= anew b .= bnew - return a, b + return a, b, stats end # Validate a symmetric soft-cover start and clear unsupported scales. diff --git a/src/sparse_support.jl b/src/sparse_support.jl index f415269..7b4e097 100644 --- a/src/sparse_support.jl +++ b/src/sparse_support.jl @@ -147,59 +147,28 @@ end # ============================================================ # Sparse `AbsLog{2}` solvers default to LSQR; `:auto` uses support density. -function symcover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) - a, _ = _symcover_min_abslog2(A; linsolve, kwargs...) - return a -end +# Otherwise they are the generic methods, which check the other keywords. +const _SparseSym = Union{SparseMatrixCSC,Symmetric{<:Any,<:SparseMatrixCSC},Hermitian{<:Any,<:SparseMatrixCSC}} -function cover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) - a, b, _ = _cover_min_abslog2(A; linsolve, kwargs...) - return a, b -end +symcover_min(ϕ::AbsLog{2}, S::_SparseSym; linsolve::Symbol=:lsqr, kwargs...) = + invoke(symcover_min, Tuple{AbsLog{2},AbstractMatrix}, ϕ, S; linsolve, kwargs...) -function symcover_min(ϕ::AbsLog{2}, S::Symmetric{<:Any, <:SparseMatrixCSC}; linsolve::Symbol=:lsqr, kwargs...) - a, _ = _symcover_min_abslog2(S; linsolve, kwargs...) - return a -end - -function symcover_min(ϕ::AbsLog{2}, H::Hermitian{<:Any, <:SparseMatrixCSC}; linsolve::Symbol=:lsqr, kwargs...) - a, _ = _symcover_min_abslog2(H; linsolve, kwargs...) - return a -end +cover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) = + invoke(cover_min, Tuple{AbsLog{2},AbstractMatrix}, ϕ, A; linsolve, kwargs...) -function cover_transversal(A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) - a, b, _ = _cover_min_abslog2(A; transversal=true, fname=:cover_transversal, linsolve, kwargs...) - return a, b -end +cover_transversal(A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) = + invoke(cover_transversal, Tuple{AbstractMatrix}, A; linsolve, kwargs...) -# The refiners take the same sparse `linsolve` default as the solvers above. -function symcover_min!(ϕ::AbsLog{2}, a::AbstractVector, - S::Union{SparseMatrixCSC,Symmetric{<:Any,<:SparseMatrixCSC},Hermitian{<:Any,<:SparseMatrixCSC}}; - linsolve::Symbol=:lsqr, kwargs...) - _prepare_symcover_start!(a, S) - anew, _ = _symcover_min_abslog2(S; start=a, linsolve, kwargs...) - a .= anew - return a -end +symcover_min!(ϕ::AbsLog{2}, a::AbstractVector, S::_SparseSym; linsolve::Symbol=:lsqr, kwargs...) = + invoke(symcover_min!, Tuple{AbsLog{2},AbstractVector,AbstractMatrix}, ϕ, a, S; linsolve, kwargs...) -function cover_min!(ϕ::AbsLog{2}, a::AbstractVector, b::AbstractVector, - A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) - _prepare_cover_start!(a, b, A) - anew, bnew, _ = _cover_min_abslog2(A; start=(a, b), linsolve, kwargs...) - a .= anew - b .= bnew - return a, b -end +cover_min!(ϕ::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::SparseMatrixCSC; + linsolve::Symbol=:lsqr, kwargs...) = + invoke(cover_min!, Tuple{AbsLog{2},AbstractVector,AbstractVector,AbstractMatrix}, ϕ, a, b, A; + linsolve, kwargs...) -# Soft covers use the same sparse default. -function soft_symcover(ϕ::AbsLog{2}, - S::Union{SparseMatrixCSC,Symmetric{<:Any,<:SparseMatrixCSC},Hermitian{<:Any,<:SparseMatrixCSC}}; - linsolve::Symbol=:lsqr, kwargs...) - a, _ = _soft_symcover_abslog2(S; linsolve, kwargs...) - return a -end +soft_symcover(ϕ::AbsLog{2}, S::_SparseSym; linsolve::Symbol=:lsqr, kwargs...) = + invoke(soft_symcover, Tuple{AbsLog{2},AbstractMatrix}, ϕ, S; linsolve, kwargs...) -function soft_cover(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) - a, b, _ = _soft_cover_abslog2(A; linsolve, kwargs...) - return a, b -end +soft_cover(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) = + invoke(soft_cover, Tuple{AbsLog{2},AbstractMatrix}, ϕ, A; linsolve, kwargs...) diff --git a/src/stats.jl b/src/stats.jl new file mode 100644 index 0000000..a86136a --- /dev/null +++ b/src/stats.jl @@ -0,0 +1,191 @@ +# Statistics returned as the last output of the optimizing cover solvers. + +""" + AbstractCoverStats + +Supertype of the statistics that the optimizing cover solvers return as their +last output, e.g. `a, stats = symcover_min(A)` or `a, b, stats = cover_min(A)`. +Each solver family has its own subtype, documented with its fields. + +A subtype implements those of [`converged`](@ref MatrixCovers.converged), +[`iterations`](@ref MatrixCovers.iterations), +[`residual`](@ref MatrixCovers.residual), and +[`tolerance`](@ref MatrixCovers.tolerance) that apply to its solver. A solver +without a stopping criterion implements none of them, so calling one throws a +`MethodError`. +""" +abstract type AbstractCoverStats end + +""" + converged(stats::AbstractCoverStats) -> Bool + +Report whether the returned cover meets the solver's stopping criterion. +""" +function converged end + +""" + iterations(stats::AbstractCoverStats) -> Int + +Return the number of iterations of the solver's main loop; each stats type +states which loop that is. +""" +function iterations end + +""" + residual(stats::AbstractCoverStats) + +Return the residual of the returned cover that the stopping criterion compares +with [`tolerance(stats)`](@ref MatrixCovers.tolerance). +""" +function residual end + +""" + tolerance(stats::AbstractCoverStats) + +Return the tolerance that [`residual(stats)`](@ref MatrixCovers.residual) is +compared with. +""" +function tolerance end + +""" + ActiveSetStats{T} + +Result of the active-set method that finishes an `AbsLog{2}` minimal-cover +solve whose multiplier iteration stopped before converging. + +# Fields + +- `certified::Bool`: the KKT conditions hold to the tolerances below, so the + result is the minimizer to rounding error. When `false`, the solver returns + the multiplier iterate instead. +- `nsteps::Int`: active-set steps, each a graph-Laplacian solve. +- `primal::T`, `ptol::T`: the largest constraint violation of the result, in + the log domain, and its tolerance. +- `dual::T`, `dtol::T`: the magnitude of the most negative multiplier of an + active constraint, and its tolerance. + +`primal` and `dual` are `NaN` when `certified` is `false`. +""" +struct ActiveSetStats{T} + certified::Bool + nsteps::Int + primal::T + ptol::T + dual::T + dtol::T +end + +""" + AugmentedLagrangianStats{T} <: AbstractCoverStats + +Statistics of the native `AbsLog{2}` minimal-cover solver, returned by +[`symcover_min`](@ref), [`cover_min`](@ref), [`cover_transversal`](@ref), and +their mutating forms. + +The multiplier iteration stops when its KKT residual reaches `tol`. If it stops +earlier, the active-set method described by `polish` finishes the solve. The +accessors describe the returned cover: `converged` is `true` when either stage +met its criterion, and `residual` and `tolerance` are those of the stage that +produced the result (`polish.primal` and `polish.ptol` after a certified finish, +otherwise `last(kkt)` and `tol`). `iterations` counts multiplier updates. + +# Fields + +- `converged::Bool`, `tol::T`: as above. +- `nouter::Int`: multiplier updates. +- `kkt::Vector{T}`, `κs::Vector{T}`, `exits::Vector{Symbol}`, `drops::Vector{T}`: + per update, the KKT residual after the update, the penalty weight, how the + inner minimization ended, and its last relative decrease of the objective. +- `polish::Union{Nothing,ActiveSetStats{T}}`: the active-set finish, or + `nothing` if it did not run. +- `linsolve::Symbol`: the inner solver (`:dense`, `:woodbury`, or `:lsqr`). +- The remaining fields count the inner linear-algebra work and describe the + preconditioner; their names and meanings may change between releases: + `nsolves`, `lsqriters`, `lsqrtrace`, `cgiters`, `cholsolves`, `precond`, + `nrefactor`, `nforest`, `ndiagonal`, `fill_entries`, `factor_flops`, + `nzeroed`. +""" +struct AugmentedLagrangianStats{T} <: AbstractCoverStats + converged::Bool + tol::T + nouter::Int + kkt::Vector{T} + κs::Vector{T} + exits::Vector{Symbol} + drops::Vector{T} + polish::Union{Nothing,ActiveSetStats{T}} + linsolve::Symbol + nsolves::Int + lsqriters::Int + lsqrtrace::Vector{Int} + cgiters::Int + cholsolves::Int + precond::Symbol + nrefactor::Int + nforest::Int + ndiagonal::Int + fill_entries::Int + factor_flops::Float64 + nzeroed::Int +end + +function _certified(s::AugmentedLagrangianStats) + p = s.polish + return p !== nothing && p.certified +end + +converged(s::AugmentedLagrangianStats) = s.converged +iterations(s::AugmentedLagrangianStats) = s.nouter +residual(s::AugmentedLagrangianStats{T}) where {T} = + _certified(s) ? s.polish.primal : isempty(s.kkt) ? T(NaN) : last(s.kkt) +tolerance(s::AugmentedLagrangianStats) = _certified(s) ? s.polish.ptol : s.tol + +""" + LeastSquaresStats <: AbstractCoverStats + +Statistics of the soft `AbsLog{2}` covers ([`soft_symcover`](@ref) and +[`soft_cover`](@ref) with `AbsLog{2}()`), which are computed by one linear +least-squares solve. There is no stopping criterion, so none of the accessors +of [`AbstractCoverStats`](@ref MatrixCovers.AbstractCoverStats) apply. + +# Fields + +- `linsolve::Symbol`: the solver (`:dense`, `:woodbury`, or `:lsqr`). +- The remaining fields count the linear-algebra work and describe the + preconditioner, as for + [`AugmentedLagrangianStats`](@ref MatrixCovers.AugmentedLagrangianStats); their + names and meanings may change between releases: `lsqriters`, `cgiters`, + `cholsolves`, `precond`, `fill_entries`, `factor_flops`. +""" +struct LeastSquaresStats <: AbstractCoverStats + linsolve::Symbol + lsqriters::Int + cgiters::Int + cholsolves::Int + precond::Symbol + fill_entries::Int + factor_flops::Float64 +end + +# Conversions from the statistics tuple of `_abslog2_auglag`, as amended by the +# minimal-cover workers. The workers compute in at least `Float64`. +_auglag_stats_type(A::AbstractMatrix) = promote_type(float(real(eltype(A))), Float64) + +AugmentedLagrangianStats(A::AbstractMatrix, nt::NamedTuple) = + AugmentedLagrangianStats{_auglag_stats_type(A)}(nt) + +function AugmentedLagrangianStats{T}(nt::NamedTuple) where {T} + p = nt.polish + polish = p === nothing ? nothing : + ActiveSetStats{T}(p.certified, p.nsteps, p.primal, p.ptol, p.dual, p.dtol) + return AugmentedLagrangianStats{T}(nt.converged, nt.tol, nt.nouter, collect(T, nt.kkt), + collect(T, nt.κs), collect(Symbol, nt.exits), + collect(T, nt.drops), polish, nt.linsolve, nt.nsolves, + nt.lsqriters, collect(Int, nt.lsqrtrace), nt.cgiters, + nt.cholsolves, nt.precond, nt.nrefactor, nt.nforest, + nt.ndiagonal, nt.fill_entries, nt.factor_flops, nt.nzeroed) +end + +LeastSquaresStats(nt::NamedTuple) = + LeastSquaresStats(nt.linsolve, nt.lsqriters, nt.cgiters, nt.cholsolves, nt.precond, + nt.fill_entries, nt.factor_flops) diff --git a/src/transversal.jl b/src/transversal.jl index a8a830b..e6a058d 100644 --- a/src/transversal.jl +++ b/src/transversal.jl @@ -45,9 +45,13 @@ finishes a solve whose multiplier iteration stops early. See also: [`cover_transversal!`](@ref), [`cover_min`](@ref). """ -function cover_transversal(A::AbstractMatrix; kwargs...) - a, b, _ = _cover_min_abslog2(A; transversal=true, fname=:cover_transversal, kwargs...) - return a, b +function cover_transversal(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::Int=AL_MAXOUTER, + maxiter::Int=AL_MAXITER, linsolve::Symbol=:auto, + fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET, + tol::Union{Real,Nothing}=nothing) + a, b, nt = _cover_min_abslog2(A; transversal=true, fname=:cover_transversal, + κ, maxouter, maxiter, linsolve, fillbudget, flopbudget, tol) + return a, b, AugmentedLagrangianStats(A, nt) end """ @@ -61,10 +65,10 @@ and `b` and returns them. The initial contents of `a` and `b` are ignored. function cover_transversal!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) axes(A, 1) == eachindex(a) || throw(DimensionMismatch("indices of `a` must match row-indexing of `A`, got eachindex(a)=$(string(eachindex(a))), axes(A, 1)=$(string(axes(A, 1)))")) axes(A, 2) == eachindex(b) || throw(DimensionMismatch("indices of `b` must match column-indexing of `A`, got eachindex(b)=$(string(eachindex(b))), axes(A, 2)=$(string(axes(A, 2)))")) - anew, bnew = cover_transversal(A; kwargs...) + anew, bnew, stats = cover_transversal(A; kwargs...) a .= anew b .= bnew - return a, b + return a, b, stats end # Binary min-heap of `(key, index)` pairs. Entries are never decreased in place; diff --git a/test/element_types.jl b/test/element_types.jl index 6b53c6a..6928a2f 100644 --- a/test/element_types.jl +++ b/test/element_types.jl @@ -5,13 +5,13 @@ @testset "Float32 flows through the family" begin A = Float32[4 1.5; 1.5 1] for a in (symcover(A), soft_symcover(A), soft_symcover(AbsLog{1}(), A), - soft_symcover(AbsLinear{1}(), A), symcover_min(AbsLog{2}(), A), - soft_symcover(AbsLog{2}(), A)) + soft_symcover(AbsLinear{1}(), A), first(symcover_min(AbsLog{2}(), A)), + first(soft_symcover(AbsLog{2}(), A))) @test eltype(a) === Float32 @test all(isfinite, a) end @test iscover(symcover(A), A; rtol=8eps(Float32)) - @test iscover(symcover_min(AbsLog{2}(), A), A; rtol=8eps(Float32)) + @test iscover(first(symcover_min(AbsLog{2}(), A)), A; rtol=8eps(Float32)) B = Float32[4 1.5 0.5; 1.5 1 2] a, b = cover(B) @@ -25,8 +25,8 @@ X = exp.(randn(rng, 60, 60)) Asym = (X .+ X') ./ 2 Agen = exp.(randn(rng, 60, 45)) - aref = symcover_min(AbsLog{2}(), Asym) - sref = soft_symcover(AbsLog{2}(), Asym) + aref, _ = symcover_min(AbsLog{2}(), Asym) + sref, _ = soft_symcover(AbsLog{2}(), Asym) gref, href = cover_min(AbsLog{2}(), Agen) A32 = Float32.(Asym) @@ -36,26 +36,26 @@ "SparseMatrixCSC" => sparse(A32), "Symmetric{SparseMatrixCSC}" => Symmetric(sparse(triu(A32))), "Hermitian{ComplexF32}" => Hermitian(ComplexF32.(A32))) - a = symcover_min(AbsLog{2}(), M) + a, _ = symcover_min(AbsLog{2}(), M) @test a isa Vector{Float32} @test a ≈ aref rtol=1e-5 end # Offset axes survive the promotion and the conversion back. Ao = OffsetArray(A32, -1, -1) - ao = symcover_min(AbsLog{2}(), Ao) + ao, _ = symcover_min(AbsLog{2}(), Ao) @test ao isa OffsetVector{Float32} @test axes(ao, 1) == axes(Ao, 1) @test collect(ao) ≈ aref rtol=1e-5 # `Float16` pins the returned cover only to its own precision, which is what # the looser tolerance measures; the solve behind it is the same `Float64` one. - a16 = symcover_min(AbsLog{2}(), Float16.(Asym)) + a16, _ = symcover_min(AbsLog{2}(), Float16.(Asym)) @test a16 isa Vector{Float16} @test a16 ≈ aref rtol=2e-3 # The soft cover is the same worker with no continuation and no boost. - s32 = soft_symcover(AbsLog{2}(), A32) + s32, _ = soft_symcover(AbsLog{2}(), A32) @test s32 isa Vector{Float32} @test s32 ≈ sref rtol=1e-5 @@ -69,16 +69,16 @@ # A type at least as wide as `Float64` is solved in itself. Abig = BigFloat.(Asym[1:8, 1:8]) - abig = symcover_min(AbsLog{2}(), Abig) + abig, _ = symcover_min(AbsLog{2}(), Abig) @test abig isa Vector{BigFloat} - @test Float64.(abig) ≈ symcover_min(AbsLog{2}(), Float64.(Abig)) rtol=1e-6 + @test Float64.(abig) ≈ first(symcover_min(AbsLog{2}(), Float64.(Abig))) rtol=1e-6 @test abig != BigFloat.(Float32.(abig)) end @testset "BigFloat flows through the family" begin A = BigFloat[4 1.5; 1.5 1] - for a in (symcover(A), soft_symcover(A), symcover_min(AbsLog{2}(), A), - soft_symcover(AbsLog{2}(), A)) + for a in (symcover(A), soft_symcover(A), first(symcover_min(AbsLog{2}(), A)), + first(soft_symcover(AbsLog{2}(), A))) @test eltype(a) === BigFloat @test all(isfinite, a) end @@ -96,7 +96,7 @@ @test iscover(a, b, A; rtol=1e-6) for ϕ in (AbsLog{1}(), AbsLinear{1}(), AbsLinear{2}()) - x, y = soft_cover!(ϕ, copy(a0), copy(b0), A) + x, y = scales(soft_cover!(ϕ, copy(a0), copy(b0), A)) @test eltype(x) === Float64 @test eltype(y) === Float32 @test all(isfinite, x) && all(isfinite, y) diff --git a/test/extensions.jl b/test/extensions.jl index 24527b5..5b9aa96 100644 --- a/test/extensions.jl +++ b/test/extensions.jl @@ -4,7 +4,7 @@ for A in ([2.0 1.0; 1.0 3.0], [100.0 1.0; 1.0 0.01], [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0]) a_fast = symcover(A) a_lmin = symcover_min(AbsLog{1}(), A) - a_qmin = symcover_min(AbsLog{2}(), A) + a_qmin, _ = symcover_min(AbsLog{2}(), A) # qmin is a valid cover @test iscover(a_qmin, A; atol=1e-10) # qmin achieves lower or equal AbsLog{2} objective than symcover and lmin @@ -18,7 +18,7 @@ @test iscover(a, A; atol=1e-10) @test a ≈ [1, 2, 1] @test abs(cover_objective(AbsLog{1}(), a, A)) < 1e-10 - a = symcover_min(AbsLog{2}(), A) + a, _ = symcover_min(AbsLog{2}(), A) @test iscover(a, A; atol=1e-10) @test a ≈ [1, 2, 1] @test abs(cover_objective(AbsLog{2}(), a, A)) < 1e-10 @@ -42,14 +42,14 @@ # soft_symcover(AbsLog{p}): unconstrained, lower objective than constrained for A in ([2.0 1.0; 1.0 3.0], [100.0 1.0; 1.0 0.01], [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0]) for ϕ in (AbsLog{1}(), AbsLog{2}()) - a_soft = soft_symcover(ϕ, A) - a_hard = symcover_min(ϕ, A) + a_soft = scales(soft_symcover(ϕ, A)) + a_hard = scales(symcover_min(ϕ, A)) @test cover_objective(ϕ, a_soft, A) <= cover_objective(ϕ, a_hard, A) + 1e-8 end end # Rank-1 matrix: both achieve zero objective A_rank1 = [2.0 1.0 4.0; 1.0 0.5 2.0; 4.0 2.0 8.0] - @test cover_objective(AbsLog{2}(), soft_symcover(AbsLog{2}(), A_rank1), A_rank1) < 1e-8 + @test cover_objective(AbsLog{2}(), first(soft_symcover(AbsLog{2}(), A_rank1)), A_rank1) < 1e-8 @test cover_objective(AbsLog{1}(), soft_symcover(AbsLog{1}(), A_rank1), A_rank1) < 1e-6 # Non-optimal solver statuses are errors. @@ -74,12 +74,12 @@ end a_fast = symcover(AbsLinear{2}(), A) for ϕ in (AbsLinear{1}(), AbsLinear{2}()) # symcover_min: valid hard cover, at most as costly as heuristic - a_min = symcover_min(ϕ, A) + a_min = scales(symcover_min(ϕ, A)) @test iscover(a_min, A; atol=1e-6) @test cover_objective(ϕ, a_min, A) <= cover_objective(ϕ, a_fast, A) + 1e-8 # soft_symcover: lower or equal objective than constrained version - a_soft = soft_symcover(ϕ, A) + a_soft = scales(soft_symcover(ϕ, A)) @test cover_objective(ϕ, a_soft, A) <= cover_objective(ϕ, a_min, A) + 1e-8 end # AbsLinear{2} ≤ AbsLinear{1} soft objectives (p=2 is a stricter lower bound) @@ -117,8 +117,8 @@ end es = sort(map(1:12) do _ a = exp.(3 .* randn(rng, 5)) b = exp.(3 .* randn(rng, 5)) - sym ? cover_objective(ϕ, soft_symcover!(ϕ, a, A), A) : - cover_objective(ϕ, soft_cover!(ϕ, a, b, A)..., A) + sym ? cover_objective(ϕ, scales(soft_symcover!(ϕ, a, A)), A) : + cover_objective(ϕ, scales(soft_cover!(ϕ, a, b, A))..., A) end) gaps = diff(es) ./ es[begin:end-1] @test all(g -> g < 1e-8 || g > 1e-3, gaps) @@ -148,18 +148,18 @@ end for ϕ in PENALTIES # Refinement does not worsen the start beyond solver tolerance. a0 = initialize_symcover(A) - a = symcover_min!(ϕ, copy(a0), A) + a = scales(symcover_min!(ϕ, copy(a0), A)) @test iscover(a, A; rtol=1e-6) @test cover_objective(ϕ, a, A) <= cover_objective(ϕ, a0, A) * (1 + 1e-6) + 1e-8 ab0, bb0 = initialize_cover(Aasym) - ab, bb = cover_min!(ϕ, copy(ab0), copy(bb0), Aasym) + ab, bb = scales(cover_min!(ϕ, copy(ab0), copy(bb0), Aasym)) @test iscover(ab, bb, Aasym; rtol=1e-6) @test cover_objective(ϕ, ab, bb, Aasym) <= cover_objective(ϕ, ab0, bb0, Aasym) * (1 + 1e-6) + 1e-8 # The asymmetric result is pinned to the balance convention, so the start is # read only up to the gauge a -> c*a, b -> b/c that leaves a[i]*b[j] fixed. - ga, gb = cover_min!(ϕ, 4 .* ab0, bb0 ./ 4, Aasym) + ga, gb = scales(cover_min!(ϕ, 4 .* ab0, bb0 ./ 4, Aasym)) @test ga ≈ ab && gb ≈ bb end @@ -180,17 +180,17 @@ end 29.4700945830419 38.0236584552296 8.96405068596511 26.5775238859338; 0.0181293142917846 0.0279078887878805 26.5775238859338 42.6650094474717] ϕ = AbsLinear{2}() - a_hard = symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin) - a_geo = symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:geomean), Abasin) + a_hard = scales(symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin)) + a_geo = scales(symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:geomean), Abasin)) @test iscover(a_hard, Abasin; rtol=1e-6) @test iscover(a_geo, Abasin; rtol=1e-6) @test cover_objective(ϕ, a_geo, Abasin) < cover_objective(ϕ, a_hard, Abasin) - 0.5 # Offset axes survive the Ipopt position mapping. Ao = OffsetArray(A, -1, -1) - ao = symcover_min!(ϕ, initialize_symcover(Ao), Ao) + ao = scales(symcover_min!(ϕ, initialize_symcover(Ao), Ao)) @test axes(ao, 1) == axes(Ao, 1) - @test collect(ao) ≈ symcover_min!(ϕ, initialize_symcover(A), A) + @test collect(ao) ≈ scales(symcover_min!(ϕ, initialize_symcover(A), A)) end @testset "soft_symcover multistart and refiner (Ipopt)" begin @@ -202,14 +202,14 @@ end @test soft_symcover!(AbsLinear{2}(), copy(a0), A) isa AbstractVector for ϕ in (AbsLinear{1}(), AbsLinear{2}(), AbsLog{1}()) - @test soft_symcover(ϕ, A) == soft_symcover(ϕ, A) # deterministic + @test scales(soft_symcover(ϕ, A)) == scales(soft_symcover(ϕ, A)) # deterministic end # No start on the menu beats the driver. for ϕ in (AbsLinear{1}(), AbsLinear{2}()) - a = soft_symcover(ϕ, A) + a = scales(soft_symcover(ϕ, A)) for strategy in (:hardcover, :geomean, :leaveout) - single = soft_symcover!(ϕ, initialize_symcover(A; strategy, feasible=:none), A) + single = scales(soft_symcover!(ϕ, initialize_symcover(A; strategy, feasible=:none), A)) @test cover_objective(ϕ, a, A) <= cover_objective(ϕ, single, A) * (1 + 1e-6) + 1e-8 end end @@ -222,7 +222,7 @@ end 3.0721768720180109 333.53567816710364 0.92571970793600067 7.3601096491681046 1.0844635712556003] d = [20.451338935482074, 0.69212569803171398, 35.401627522529395, 15.904906661932396, 0.19696509774727827] for ϕ in (AbsLinear{1}(), AbsLinear{2}()) - @test covaries(A -> soft_symcover(ϕ, A), Abasins, d; rtol=1e-5) + @test covaries(A -> scales(soft_symcover(ϕ, A)), Abasins, d; rtol=1e-5) end @test_throws "positive scale on every supported row" soft_symcover!(AbsLinear{2}(), [1.0, -1.0, 2.0], A) @@ -241,30 +241,30 @@ end @test_throws "requires a start that covers `A`" cover_min!(AbsLinear{2}(), copy(a0), copy(b0), A) for ϕ in (AbsLinear{1}(), AbsLinear{2}(), AbsLog{1}(), AbsLog{2}()) - a, b = soft_cover(ϕ, A) - @test soft_cover(ϕ, A) == (a, b) # deterministic + a, b = scales(soft_cover(ϕ, A)) + @test scales(soft_cover(ϕ, A)) == (a, b) # deterministic # The objective depends on `a`, `b` only through their products, so the split is # fixed by the balance convention rather than left to the solver. @test balance(a, b) ≈ 0 atol=1e-7 # Gauge-invariant in the start: (a, b) and (2a, b/2) name the same point. - ag, bg = soft_cover!(ϕ, 2 .* copy(a0), copy(b0) ./ 2, A) - ah, bh = soft_cover!(ϕ, copy(a0), copy(b0), A) + ag, bg = scales(soft_cover!(ϕ, 2 .* copy(a0), copy(b0) ./ 2, A)) + ah, bh = scales(soft_cover!(ϕ, copy(a0), copy(b0), A)) @test ag ≈ ah && bg ≈ bh end # No start on the menu beats the driver. for ϕ in (AbsLinear{1}(), AbsLinear{2}()) - a, b = soft_cover(ϕ, A) + a, b = scales(soft_cover(ϕ, A)) for strategy in (:hardcover, :covariant) - sa, sb = soft_cover!(ϕ, initialize_cover(A; strategy, feasible=:none)..., A) + sa, sb = scales(soft_cover!(ϕ, initialize_cover(A; strategy, feasible=:none)..., A)) @test cover_objective(ϕ, a, b, A) <= cover_objective(ϕ, sa, sb, A) * (1 + 1e-6) + 1e-8 end # The asymmetric soft cover relaxes the symmetric one on a symmetric matrix. As = [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0] - aa, ab = soft_cover(ϕ, As) + aa, ab = scales(soft_cover(ϕ, As)) @test cover_objective(ϕ, aa, ab, As) <= - cover_objective(ϕ, soft_symcover(ϕ, As), As) * (1 + 1e-6) + 1e-8 + cover_objective(ϕ, scales(soft_symcover(ϕ, As)), As) * (1 + 1e-6) + 1e-8 end # Covariant starts make the selected cover covariant on irregular support too. @@ -290,36 +290,36 @@ end B = exp.(3 .* randn(rng, n, n)) .* (rand(rng, n, n) .> 0.25) A = (B + B') / 2 all(iszero, A) && continue - a = symcover_min(ϕ1, A) + a = scales(symcover_min(ϕ1, A)) @test iscover(a, A; rtol=1e-7) # No cover beats it on AbsLog{1}: the canonical selection does not cost optimality. # Any feasible cover is a witness; use the AbsLog{2}-minimal one, and the heuristic. - for other in (symcover_min(ϕ2, A), symcover(A)) + for other in (scales(symcover_min(ϕ2, A)), symcover(A)) @test cover_objective(ϕ1, a, A) <= cover_objective(ϕ1, other, A) * (1 + 1e-6) + 1e-8 end # Start-independent in the *result*, not merely in the objective value. for strategy in (:hardcover, :geomean, :diagfeasible) - @test symcover_min!(ϕ1, initialize_symcover(A; strategy), A) ≈ a + @test scales(symcover_min!(ϕ1, initialize_symcover(A; strategy), A)) ≈ a end # Scale-covariant: both objectives see `A` only through the residuals, which a # rescaling leaves invariant, so the selected member co-varies with the frame. D = Diagonal(exp.(randn(rng, n))) - @test symcover_min(ϕ1, D * A * D) ≈ D * a + @test scales(symcover_min(ϕ1, D * A * D)) ≈ D * a end # Gauge balancing and `AbsLog{1}` face selection are independent. Aasym = [3.0 1.0 7.0; 2.0 5.0 1.0; 8.0 1.0 4.0] - a, b = cover_min(ϕ1, Aasym) + a, b = scales(cover_min(ϕ1, Aasym)) nza = vec(count(!iszero, Aasym, dims=2)) nzb = vec(count(!iszero, Aasym, dims=1)) @test sum(nza .* log.(a)) ≈ sum(nzb .* log.(b)) atol=1e-8 @test iscover(a, b, Aasym; rtol=1e-7) # Gauge-invariant in the start, and start-independent in the result. a0, b0 = initialize_cover(Aasym; strategy=:geomean) - ag, bg = cover_min!(ϕ1, 4 .* a0, b0 ./ 4, Aasym) + ag, bg = scales(cover_min!(ϕ1, 4 .* a0, b0 ./ 4, Aasym)) @test ag ≈ a && bg ≈ b end @@ -332,43 +332,43 @@ end Aasym = [1.0 2.0 3.0; 40.0 5.0 0.6] for ϕ in (AbsLinear{1}(), AbsLinear{2}()) - a = symcover_min(ϕ, Abasin) + a = scales(symcover_min(ϕ, Abasin)) @test iscover(a, Abasin; rtol=1e-6) # No start on the menu beats the driver. for strategy in (:hardcover, :geomean, :leaveout) - single = symcover_min!(ϕ, initialize_symcover(Abasin; strategy), Abasin) + single = scales(symcover_min!(ϕ, initialize_symcover(Abasin; strategy), Abasin)) @test cover_objective(ϕ, a, Abasin) <= cover_objective(ϕ, single, Abasin) * (1 + 1e-6) + 1e-8 end - @test symcover_min(ϕ, Abasin) == a # deterministic + @test scales(symcover_min(ϕ, Abasin)) == a # deterministic - ab, bb = cover_min(ϕ, Aasym) + ab, bb = scales(cover_min(ϕ, Aasym)) @test iscover(ab, bb, Aasym; rtol=1e-6) for strategy in (:hardcover, :covariant) - sa, sb = cover_min!(ϕ, initialize_cover(Aasym; strategy)..., Aasym) + sa, sb = scales(cover_min!(ϕ, initialize_cover(Aasym; strategy)..., Aasym)) @test cover_objective(ϕ, ab, bb, Aasym) <= cover_objective(ϕ, sa, sb, Aasym) * (1 + 1e-6) + 1e-8 end - @test cover_min(ϕ, Aasym) == (ab, bb) # deterministic + @test scales(cover_min(ϕ, Aasym)) == (ab, bb) # deterministic # The asymmetric cover relaxes the symmetric one: independent row and column # scales can only do better on the same matrix. - asym_a, asym_b = cover_min(ϕ, Abasin) + asym_a, asym_b = scales(cover_min(ϕ, Abasin)) @test cover_objective(ϕ, asym_a, asym_b, Abasin) <= - cover_objective(ϕ, symcover_min(ϕ, Abasin), Abasin) * (1 + 1e-6) + 1e-8 + cover_objective(ϕ, scales(symcover_min(ϕ, Abasin)), Abasin) * (1 + 1e-6) + 1e-8 end # Restricting `strategies` changes the selected basin. ϕ = AbsLinear{2}() - @test symcover_min(ϕ, Abasin; strategies=(:hardcover,)) ≈ - symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin) - @test cover_objective(ϕ, symcover_min(ϕ, Abasin), Abasin) < - cover_objective(ϕ, symcover_min(ϕ, Abasin; strategies=(:hardcover,)), Abasin) - 0.5 + @test scales(symcover_min(ϕ, Abasin; strategies=(:hardcover,))) ≈ + scales(symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin)) + @test cover_objective(ϕ, scales(symcover_min(ϕ, Abasin)), Abasin) < + cover_objective(ϕ, scales(symcover_min(ϕ, Abasin; strategies=(:hardcover,))), Abasin) - 0.5 # A matrix whose every row carries a single support entry admits no :leaveout start; # that strategy forfeits its slot rather than failing the solve. Adiag = [4.0 0.0; 0.0 9.0] - @test symcover_min(ϕ, Adiag) ≈ [2.0, 3.0] rtol=1e-4 + @test scales(symcover_min(ϕ, Adiag)) ≈ [2.0, 3.0] rtol=1e-4 @test_throws "unknown strategy :banana" symcover_min(ϕ, Abasin; strategies=(:banana,)) @test_throws "no strategy in (:leaveout,)" symcover_min(ϕ, Adiag; strategies=(:leaveout,)) @@ -384,7 +384,7 @@ end # Do not hint when loading an extension cannot fix the argument type. A = [4.0 1.0; 1.0 4.0] e = try - symcover_min(AbsLog{2}(), "not a matrix") + first(symcover_min(AbsLog{2}(), "not a matrix")) nothing catch err err @@ -457,7 +457,7 @@ end # The widened gate still refuses to fire when no package load could help. e5 = try - symcover_min!(AbsLog{2}(), [1.0, 1.0], "not a matrix") + first(symcover_min!(AbsLog{2}(), [1.0, 1.0], "not a matrix")) nothing catch err err @@ -475,7 +475,7 @@ end dense[i, j] = dense[j, i] = v end - @test symcover_min(AbsLog{2}(), M) ≈ symcover_min(AbsLog{2}(), dense) rtol=1e-8 + @test first(symcover_min(AbsLog{2}(), M)) ≈ first(symcover_min(AbsLog{2}(), dense)) rtol=1e-8 @test MatrixCovers.symcover_min_jump(AbsLog{2}(), M) ≈ MatrixCovers.symcover_min_jump(AbsLog{2}(), dense) rtol=1e-6 @test symcover_min(AbsLog{1}(), M) ≈ symcover_min(AbsLog{1}(), dense) rtol=1e-6 @@ -487,7 +487,7 @@ end @test fh[1] ≈ fd[1] rtol=1e-6 @test fh[2] ≈ fd[2] rtol=1e-6 end - for (fh, fd) in ((cover_transversal(M), cover_transversal(dense)), + for (fh, fd) in ((scales(cover_transversal(M)), scales(cover_transversal(dense))), (MatrixCovers.cover_transversal_jump(M)[1], MatrixCovers.cover_transversal_jump(dense)[1])) @test fh[1] ≈ fd[1] rtol=1e-6 @@ -497,20 +497,20 @@ end # Ipopt: the AbsLinear kernels, whose starts the caller supplies. for p in (1, 2) ϕ = AbsLinear{p}() - sh = symcover_min!(ϕ, symcover(ϕ, M), M) - sd = symcover_min!(ϕ, symcover(ϕ, dense), dense) + sh = scales(symcover_min!(ϕ, symcover(ϕ, M), M)) + sd = scales(symcover_min!(ϕ, symcover(ϕ, dense), dense)) @test sh ≈ sd rtol=1e-5 - ah, bh = cover_min!(ϕ, cover(ϕ, M)..., M) - ad, bd = cover_min!(ϕ, cover(ϕ, dense)..., dense) + ah, bh = scales(cover_min!(ϕ, cover(ϕ, M)..., M)) + ad, bd = scales(cover_min!(ϕ, cover(ϕ, dense)..., dense)) @test ah ≈ ad rtol=1e-5 @test bh ≈ bd rtol=1e-5 - @test soft_symcover!(ϕ, soft_symcover(ϕ, M), M) ≈ - soft_symcover!(ϕ, soft_symcover(ϕ, dense), dense) rtol=1e-5 + @test scales(soft_symcover!(ϕ, scales(soft_symcover(ϕ, M)), M)) ≈ + scales(soft_symcover!(ϕ, scales(soft_symcover(ϕ, dense)), dense)) rtol=1e-5 - sah, sbh = soft_cover!(ϕ, soft_cover(ϕ, M)..., M) - sad, sbd = soft_cover!(ϕ, soft_cover(ϕ, dense)..., dense) + sah, sbh = scales(soft_cover!(ϕ, scales(soft_cover(ϕ, M))..., M)) + sad, sbd = scales(soft_cover!(ϕ, scales(soft_cover(ϕ, dense))..., dense)) @test sah ≈ sad rtol=1e-5 @test sbh ≈ sbd rtol=1e-5 end diff --git a/test/helpers.jl b/test/helpers.jl index a7c3e50..f854484 100644 --- a/test/helpers.jl +++ b/test/helpers.jl @@ -6,6 +6,13 @@ const PENALTIES = (AbsLog{1}(), AbsLog{2}(), AbsLinear{1}(), AbsLinear{2}()) # Penalties with native soft-cover solvers, for properties of the soft-cover family. const SOFT_PENALTIES = (PENALTIES..., PowerMean{1}(), PowerMean{2}(), PowerMean{3.5}()) +# The scale vectors of a cover function's result: optimizing solvers append their +# statistics, heuristics do not. Any other shape is a `MethodError`. +scales(a::AbstractVector) = a +scales(r::Tuple{AbstractVector,MatrixCovers.AbstractCoverStats}) = r[1] +scales(r::Tuple{AbstractVector,AbstractVector}) = r +scales(r::Tuple{AbstractVector,AbstractVector,MatrixCovers.AbstractCoverStats}) = (r[1], r[2]) + """ isbalanced(a, b, A) diff --git a/test/heuristic_covers.jl b/test/heuristic_covers.jl index 1959edf..41a5164 100644 --- a/test/heuristic_covers.jl +++ b/test/heuristic_covers.jl @@ -199,14 +199,14 @@ end n = rand(qrng, 5:15) B = randn(qrng, n, n) .* exp.(rand(qrng) * 3 * randn(qrng, n, n)) A = (B + B') / 2 - Emin = cover_objective(AbsLog{2}(), symcover_min(AbsLog{2}(), A), A) + Emin = cover_objective(AbsLog{2}(), first(symcover_min(AbsLog{2}(), A)), A) E3 = cover_objective(AbsLog{2}(), symcover(AbsLog{2}(), A; maxiter=3), A) iszero(Emin) || push!(gaps, log(E3 / Emin)) end for _ in 1:10 n = rand(qrng, 5:15) S = sprand(qrng, n, n, 0.3); A = Matrix(S + S') - Emin = cover_objective(AbsLog{2}(), symcover_min(AbsLog{2}(), A), A) + Emin = cover_objective(AbsLog{2}(), first(symcover_min(AbsLog{2}(), A)), A) E3 = cover_objective(AbsLog{2}(), symcover(AbsLog{2}(), A; maxiter=3), A) iszero(Emin) || push!(gaps, log(E3 / Emin)) end diff --git a/test/invariants.jl b/test/invariants.jl index 8e15bf5..ae99263 100644 --- a/test/invariants.jl +++ b/test/invariants.jl @@ -8,11 +8,11 @@ const SYM_NOTIONS = ( (name = "symcover", f = A -> symcover(A), hard = true, rtol = 1e-9), - (name = "symcover_min(AbsLog{2})", f = A -> symcover_min(AbsLog{2}(), A), hard = true, rtol = 1e-5), + (name = "symcover_min(AbsLog{2})", f = A -> first(symcover_min(AbsLog{2}(), A)), hard = true, rtol = 1e-5), (name = "symcover_min(AbsLog{1})", f = A -> symcover_min(AbsLog{1}(), A), hard = true, rtol = 1e-5), (name = "symcover_min(AbsLinear{1})", f = A -> symcover_min(AbsLinear{1}(), A), hard = true, rtol = 1e-5), (name = "symcover_min(AbsLinear{2})", f = A -> symcover_min(AbsLinear{2}(), A), hard = true, rtol = 1e-5), - (name = "soft_symcover(AbsLog{2})", f = A -> soft_symcover(AbsLog{2}(), A), hard = false, rtol = 1e-9), + (name = "soft_symcover(AbsLog{2})", f = A -> first(soft_symcover(AbsLog{2}(), A)), hard = false, rtol = 1e-9), (name = "soft_symcover(AbsLog{1})", f = A -> soft_symcover(AbsLog{1}(), A), hard = false, rtol = 1e-5), (name = "soft_symcover(AbsLinear{2})", f = A -> soft_symcover(AbsLinear{2}(), A), hard = false, rtol = 1e-5), (name = "soft_symcover(AbsLinear{1})", f = A -> soft_symcover(AbsLinear{1}(), A), hard = false, rtol = 1e-5), @@ -22,11 +22,11 @@ const SYM_NOTIONS = ( const GEN_NOTIONS = ( (name = "cover", f = A -> cover(A), hard = true, rtol = 1e-9), - (name = "cover_min(AbsLog{2})", f = A -> cover_min(AbsLog{2}(), A), hard = true, rtol = 1e-5), + (name = "cover_min(AbsLog{2})", f = A -> scales(cover_min(AbsLog{2}(), A)), hard = true, rtol = 1e-5), (name = "cover_min(AbsLog{1})", f = A -> cover_min(AbsLog{1}(), A), hard = true, rtol = 1e-5), (name = "cover_min(AbsLinear{1})", f = A -> cover_min(AbsLinear{1}(), A), hard = true, rtol = 1e-5), (name = "cover_min(AbsLinear{2})", f = A -> cover_min(AbsLinear{2}(), A), hard = true, rtol = 1e-5), - (name = "soft_cover(AbsLog{2})", f = A -> soft_cover(AbsLog{2}(), A), hard = false, rtol = 1e-9), + (name = "soft_cover(AbsLog{2})", f = A -> scales(soft_cover(AbsLog{2}(), A)), hard = false, rtol = 1e-9), (name = "soft_cover(AbsLog{1})", f = A -> soft_cover(AbsLog{1}(), A), hard = false, rtol = 1e-5), (name = "soft_cover(AbsLinear{2})", f = A -> soft_cover(AbsLinear{2}(), A), hard = false, rtol = 1e-5), (name = "soft_cover(AbsLinear{1})", f = A -> soft_cover(AbsLinear{1}(), A), hard = false, rtol = 1e-5), @@ -164,7 +164,7 @@ const GEN_NOTIONS = ( for M in mats # AbsLog{2}: native solver and the HiGHS reference model. o2 = ref_min(2, M) - @test cover_objective(AbsLog{2}(), symcover_min(AbsLog{2}(), M), M) ≈ o2 rtol=1e-5 + @test cover_objective(AbsLog{2}(), first(symcover_min(AbsLog{2}(), M)), M) ≈ o2 rtol=1e-5 @test cover_objective(AbsLog{2}(), MatrixCovers.symcover_min_jump(AbsLog{2}(), M), M) ≈ o2 rtol=1e-5 # AbsLog{1}: the HiGHS LP, whose objective is assembled from nonzero counts # rather than from the residuals directly. @@ -178,12 +178,12 @@ end certified_sym = ( A -> symcover(A), A -> symcover(A; maxiter=0), - A -> symcover_min(AbsLog{2}(), A), + A -> first(symcover_min(AbsLog{2}(), A)), ) certified_gen = ( A -> cover(A), A -> cover(A; maxiter=0), - A -> cover_min(AbsLog{2}(), A), + A -> scales(cover_min(AbsLog{2}(), A)), ) # Exercise both flattened-support and dense-grid kernels. for T in (Float64, Float32), σ in (1, 5), n in (7, 70), seed in 1:3 diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index aa9f817..c19b54c 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -12,7 +12,7 @@ end # Match HiGHS across the symmetric corpus. for (_, A) in symmetric_matrices Af = Float64.(A) - a = symcover_min(AbsLog{2}(), Af) + a, _ = symcover_min(AbsLog{2}(), Af) aj = MatrixCovers.symcover_min_jump(AbsLog{2}(), Af) @test iscover(a, Af; atol=1e-8) oj = cover_objective(AbsLog{2}(), aj, Af) @@ -27,22 +27,22 @@ end for (_, A) in symmetric_matrices[1:30] Af = Float64.(A); n = size(Af, 1) d = exp.(2 .* randn(rng, n)) - @test covaries(A -> symcover_min(AbsLog{2}(), A), Af, d; rtol=1e-6) + @test covaries(A -> first(symcover_min(AbsLog{2}(), A)), Af, d; rtol=1e-6) end # Edge cases. - @test symcover_min(AbsLog{2}(), reshape([4.0], 1, 1)) ≈ [2.0] # n = 1 + @test first(symcover_min(AbsLog{2}(), reshape([4.0], 1, 1))) ≈ [2.0] # n = 1 # [0 1; 1 0]: a₁a₂ = 1 is the (gauge-invariant) optimum, objective 0. - a = symcover_min(AbsLog{2}(), [0.0 1.0; 1.0 0.0]) + a, _ = symcover_min(AbsLog{2}(), [0.0 1.0; 1.0 0.0]) @test a[1] * a[2] ≈ 1.0 @test cover_objective(AbsLog{2}(), a, [0.0 1.0; 1.0 0.0]) < 1e-12 # Scattered zeros with an exact rank-1 cover. A = [0 0 1; 0 0 2; 1 2 1] - a = symcover_min(AbsLog{2}(), A) + a, _ = symcover_min(AbsLog{2}(), A) @test a ≈ [1, 2, 1] @test abs(cover_objective(AbsLog{2}(), a, A)) < 1e-8 # κ and maxouter keywords are accepted. - @test symcover_min(AbsLog{2}(), [2.0 1.0; 1.0 3.0]; κ=1e3, maxouter=16) isa Vector + @test first(symcover_min(AbsLog{2}(), [2.0 1.0; 1.0 3.0]; κ=1e3, maxouter=16)) isa Vector end @testset "symcover_min active-set finish on dense supports" begin @@ -75,27 +75,40 @@ end # Offset axes carry through the finish. Ao = OffsetArray(A, -5, -5) - ao = symcover_min(AbsLog{2}(), Ao; maxouter=3) + ao, _ = symcover_min(AbsLog{2}(), Ao; maxouter=3) @test axes(ao, 1) == axes(Ao, 1) - @test collect(ao) ≈ symcover_min(AbsLog{2}(), A; maxouter=3) rtol=1e-12 + @test collect(ao) ≈ first(symcover_min(AbsLog{2}(), A; maxouter=3)) rtol=1e-12 # Without the finish, the stalled iteration warns; an uncertified finish - # says so. + # says so. Any residual above the tolerance warns. @test_logs (:warn, r"Increase `maxouter` or `κ`") MatrixCovers._symcover_min_abslog2(A; maxouter=3, polish=false) - stats = (; converged=false, kkt=(1.0,), vtol=1e-13, vwarn=1e-10, nouter=3, - polish=(; certified=false, nsteps=5)) + stats = (; converged=false, kkt=(1e-12,), tol=1e-13, nouter=32, polish=nothing) + @test_logs (:warn, r"Increase `maxouter` or `κ`") MatrixCovers._warn_unconverged(:symcover_min, stats, 32) + stats = merge(stats, (; polish=(; certified=false, nsteps=5))) @test_logs (:warn, r"active-set finish .* did not certify") MatrixCovers._warn_unconverged(:symcover_min, stats, 32) - # An uncertified finish warns even when the residual is below `vwarn`. - stats = merge(stats, (; kkt=(1e-12,))) - @test_logs (:warn, r"active-set finish .* did not certify") MatrixCovers._warn_unconverged(:symcover_min, stats, 32) - # A certified finish reports a KKT check that meets its tolerances. - _, s = MatrixCovers._symcover_min_abslog2(A; maxouter=3) - @test s.converged && s.polish.certified - @test s.polish.primal <= s.polish.ptol && s.polish.dual <= s.polish.dtol + # A certified finish reports a KKT check that meets its tolerances, and the + # accessors describe it. + _, s = @test_logs symcover_min(AbsLog{2}(), A; maxouter=3) + @test s isa MatrixCovers.AugmentedLagrangianStats{Float64} + @test MatrixCovers.converged(s) && s.polish.certified + @test MatrixCovers.residual(s) == s.polish.primal <= MatrixCovers.tolerance(s) == s.polish.ptol + @test s.polish.dual <= s.polish.dtol + @test MatrixCovers.iterations(s) == s.nouter <= 3 + # A looser tolerance ends the multiplier iteration sooner, without the finish. + _, s = symcover_min(AbsLog{2}(), A; tol=1e-4) + @test MatrixCovers.converged(s) && s.polish === nothing + @test MatrixCovers.tolerance(s) == 1e-4 && MatrixCovers.residual(s) == last(s.kkt) <= 1e-4 + @test_throws "tol must be finite and nonnegative" symcover_min(AbsLog{2}(), A; tol=-1.0) # With no outer iteration nothing is checked, so the result is not # reported as converged, and there is no warning. - _, s0 = @test_logs MatrixCovers._symcover_min_abslog2(A; maxouter=0, boost=false, fname=:soft_symcover) - @test !s0.converged && s0.nouter == 0 && !haskey(s0, :polish) + _, s0 = @test_logs symcover_min(AbsLog{2}(), A; maxouter=0) + @test !MatrixCovers.converged(s0) && MatrixCovers.iterations(s0) == 0 + @test isnan(MatrixCovers.residual(s0)) && s0.polish === nothing + # The soft cover has no stopping criterion, so the accessors do not apply. + _, sl = @test_logs soft_symcover(AbsLog{2}(), A) + @test sl isa MatrixCovers.LeastSquaresStats + @test_throws MethodError MatrixCovers.converged(sl) + @test_throws MethodError MatrixCovers.residual(sl) end @testset "cover_min native AbsLog{2}" begin @@ -122,7 +135,7 @@ end for (_, A) in general_matrices[1:30] Af = Float64.(A); m, n = size(Af) dr = exp.(2 .* randn(rng, m)); dc = exp.(2 .* randn(rng, n)) - @test covaries(A -> cover_min(AbsLog{2}(), A), Af, dr, dc; rtol=1e-6) + @test covaries(A -> scales(cover_min(AbsLog{2}(), A)), Af, dr, dc; rtol=1e-6) end # Non-square matrices, both orientations (transpose swaps the roles of a, b). @@ -148,7 +161,7 @@ end @test iscover(a, b, A; atol=1e-8) @test cover_objective(AbsLog{2}(), a, b, A) ≈ 2 * log(sqrt(2))^2 # κ and maxouter keywords are accepted. - @test cover_min(AbsLog{2}(), [1.0 2.0; 3.0 4.0]; κ=1e3, maxouter=16) isa Tuple + @test scales(cover_min(AbsLog{2}(), [1.0 2.0; 3.0 4.0]; κ=1e3, maxouter=16)) isa Tuple end @testset "MMC native AbsLog{2} matrix-free LSQR path" begin @@ -160,7 +173,7 @@ end # across the committed symmetric library, and returns a feasible cover. for (_, A) in symmetric_matrices Af = Float64.(A) - a = symcover_min(AbsLog{2}(), Af; linsolve=:lsqr) + a, _ = symcover_min(AbsLog{2}(), Af; linsolve=:lsqr) aj = MatrixCovers.symcover_min_jump(AbsLog{2}(), Af) @test iscover(a, Af; atol=1e-8) @test cover_objective(AbsLog{2}(), a, Af) <= @@ -180,11 +193,11 @@ end # Gauge/edge cases the dense path handles via a ridge or v0*v0ᵀ must also work # matrix-free: bipartite support, a scalar, and a zero row/column. - a = symcover_min(AbsLog{2}(), [0.0 1.0; 1.0 0.0]; linsolve=:lsqr) + a, _ = symcover_min(AbsLog{2}(), [0.0 1.0; 1.0 0.0]; linsolve=:lsqr) @test a[1] * a[2] ≈ 1.0 - @test symcover_min(AbsLog{2}(), reshape([4.0], 1, 1); linsolve=:lsqr) ≈ [2.0] + @test first(symcover_min(AbsLog{2}(), reshape([4.0], 1, 1); linsolve=:lsqr)) ≈ [2.0] Az = [1.0 0.0 2.0; 0.0 0.0 0.0; 2.0 0.0 3.0] - a = symcover_min(AbsLog{2}(), Az; linsolve=:lsqr) + a, _ = symcover_min(AbsLog{2}(), Az; linsolve=:lsqr) @test a[2] == 0.0 @test iscover(a, Az; atol=1e-8) a, b = cover_min(AbsLog{2}(), [0.0 1.0; 1.0 0.0]; linsolve=:lsqr) @@ -272,14 +285,14 @@ end # Offset axes and views index the support through `axes(A)`, not `1:n`, on this # path as on the others. A = symlognormal(12) - aref = symcover_min(AbsLog{2}(), A; linsolve=:woodbury) + aref, _ = symcover_min(AbsLog{2}(), A; linsolve=:woodbury) Ao = OffsetArray(A, -3, -3) - ao = symcover_min(AbsLog{2}(), Ao; linsolve=:woodbury) + ao, _ = symcover_min(AbsLog{2}(), Ao; linsolve=:woodbury) @test axes(ao, 1) == axes(Ao, 1) @test collect(ao) ≈ aref rtol=1e-10 Av = view(symlognormal(16), 3:14, 3:14) - @test symcover_min(AbsLog{2}(), Matrix(Av); linsolve=:woodbury) ≈ - symcover_min(AbsLog{2}(), Av; linsolve=:woodbury) rtol=1e-10 + @test first(symcover_min(AbsLog{2}(), Matrix(Av); linsolve=:woodbury)) ≈ + first(symcover_min(AbsLog{2}(), Av; linsolve=:woodbury)) rtol=1e-10 Ag = lognormal(14, 10) Agv = view(Ag, 2:13, 2:9) av, bv = cover_min(AbsLog{2}(), Agv; linsolve=:woodbury) @@ -342,7 +355,7 @@ end A32 = Float32.(symlognormal(8)) a32, s32 = MatrixCovers._symcover_min_abslog2(A32) @test s32.linsolve === :woodbury && eltype(a32) === Float32 - @test symcover_min(AbsLog{2}(), A32; linsolve=:woodbury) ≈ symcover_min(AbsLog{2}(), A32; linsolve=:dense) rtol=1e-6 + @test first(symcover_min(AbsLog{2}(), A32; linsolve=:woodbury)) ≈ first(symcover_min(AbsLog{2}(), A32; linsolve=:dense)) rtol=1e-6 Abig = BigFloat.(symlognormal(8)) @test_throws "requires Float64 arithmetic" symcover_min(AbsLog{2}(), Abig; linsolve=:woodbury) @test MatrixCovers._symcover_min_abslog2(Abig)[2].linsolve === :dense @@ -447,7 +460,7 @@ end # The budget is a keyword of the public solvers, and the paths that never # precondition say so. - @test symcover_min(AbsLog{2}(), A; linsolve=:lsqr, fillbudget=0) ≈ ad rtol=1e-6 + @test first(symcover_min(AbsLog{2}(), A; linsolve=:lsqr, fillbudget=0)) ≈ ad rtol=1e-6 @test MatrixCovers._symcover_min_abslog2(A; linsolve=:dense)[2].precond === :none Abig = BigFloat.([4.0 1.0 0.5; 1.0 3.0 1.0; 0.5 1.0 2.5]) @test MatrixCovers._symcover_min_abslog2(Abig; linsolve=:lsqr)[2].precond === :none @@ -612,7 +625,7 @@ end @test (tl.nsolves, tl.lsqriters) == (th.nsolves, th.lsqriters) # With no outer passes the unweighted fit is the answer, not a start. - @test soft_symcover(AbsLog{2}(), A; linsolve=:lsqr) != symcover(A) + @test first(soft_symcover(AbsLog{2}(), A; linsolve=:lsqr)) != symcover(A) end # Bound iterations for Cholesky-preconditioned `Float64` LSQR. @@ -643,9 +656,9 @@ end # A working type CHOLMOD cannot factor keeps the plain matrix-free iteration. A32 = Float32.([4.0 1.0 0.5; 1.0 3.0 1.0; 0.5 1.0 2.5]) - a32 = symcover_min(AbsLog{2}(), A32; linsolve=:lsqr) + a32, _ = symcover_min(AbsLog{2}(), A32; linsolve=:lsqr) @test a32 isa Vector{Float32} - @test a32 ≈ symcover_min(AbsLog{2}(), A32; linsolve=:dense) rtol=1e-5 + @test a32 ≈ first(symcover_min(AbsLog{2}(), A32; linsolve=:dense)) rtol=1e-5 @test iscover(a32, A32; rtol=1e-5) end @@ -694,23 +707,23 @@ end rng = StableRNG(42) n = 6 M = randn(rng, ComplexF64, n, n) - Href = symcover_min(AbsLog{2}(), abs.(Matrix(Hermitian(M + M')))) + Href, _ = symcover_min(AbsLog{2}(), abs.(Matrix(Hermitian(M + M')))) Hdense = Hermitian(M + M') - @test symcover_min(AbsLog{2}(), Hdense) ≈ Href rtol = 1e-10 + @test first(symcover_min(AbsLog{2}(), Hdense)) ≈ Href rtol = 1e-10 Msp = sprandn(rng, ComplexF64, n, n, 0.5) Hsp = Hermitian(sparse(Msp + Msp')) - @test symcover_min(AbsLog{2}(), Hsp) ≈ symcover_min(AbsLog{2}(), abs.(Matrix(Hsp))) rtol = 1e-7 + @test first(symcover_min(AbsLog{2}(), Hsp)) ≈ first(symcover_min(AbsLog{2}(), abs.(Matrix(Hsp)))) rtol = 1e-7 # Symmetric{<:Complex} carries no conjugate-symmetry guarantee, but the cover # problem still only depends on abs.(A), so the same identity holds. Ssp = Symmetric(sparse(Msp + transpose(Msp))) - @test symcover_min(AbsLog{2}(), Ssp) ≈ symcover_min(AbsLog{2}(), abs.(Matrix(Ssp))) rtol = 1e-7 + @test first(symcover_min(AbsLog{2}(), Ssp)) ≈ first(symcover_min(AbsLog{2}(), abs.(Matrix(Ssp)))) rtol = 1e-7 # Real input is unaffected by deriving T from real(eltype(A)). Ar = [4.0 1.0; 1.0 4.0] - @test symcover_min(AbsLog{2}(), Ar) ≈ [2.0, 2.0] + @test first(symcover_min(AbsLog{2}(), Ar)) ≈ [2.0, 2.0] end @testset "soft_cover native AbsLog{2}" begin @@ -742,9 +755,9 @@ end # symcover_min(A) and cover_min(A) default to AbsLog{2}, matching # symcover(A)/cover(A). Solved natively, so no extension is needed. A = [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0] - @test symcover_min(A) == symcover_min(AbsLog{2}(), A) + @test first(symcover_min(A)) == first(symcover_min(AbsLog{2}(), A)) Aasym = [1.0 2.0 3.0; 4.0 5.0 6.0] - @test cover_min(Aasym) == cover_min(AbsLog{2}(), Aasym) + @test scales(cover_min(Aasym)) == scales(cover_min(AbsLog{2}(), Aasym)) end @testset "symcover_min!/cover_min! native AbsLog{2}" begin @@ -754,23 +767,23 @@ end # AbsLog{2} has a unique minimizer, so refining any valid start reproduces the # cold solve exactly, and the no-ϕ form selects the same penalty. a = initialize_symcover(A) - @test symcover_min!(AbsLog{2}(), copy(a), A) ≈ symcover_min(AbsLog{2}(), A) - @test symcover_min!(copy(a), A) == symcover_min!(AbsLog{2}(), copy(a), A) + @test first(symcover_min!(AbsLog{2}(), copy(a), A)) ≈ first(symcover_min(AbsLog{2}(), A)) + @test first(symcover_min!(copy(a), A)) == first(symcover_min!(AbsLog{2}(), copy(a), A)) ab, bb = initialize_cover(Aasym) ra, rb = cover_min!(AbsLog{2}(), copy(ab), copy(bb), Aasym) ca, cb = cover_min(AbsLog{2}(), Aasym) @test ra ≈ ca && rb ≈ cb - @test cover_min!(copy(ab), copy(bb), Aasym) == (ra, rb) + @test scales(cover_min!(copy(ab), copy(bb), Aasym)) == (ra, rb) # The refiners write through the buffers they are handed. a = initialize_symcover(A) - @test symcover_min!(AbsLog{2}(), a, A) === a + @test first(symcover_min!(AbsLog{2}(), a, A)) === a @test iscover(a, A; atol=1e-9) # The heuristics land on the coverage boundary only to within the roundoff of # their log-domain updates, so their output must be accepted as a start. - @test symcover_min!(AbsLog{2}(), symcover(A), A) ≈ symcover_min(AbsLog{2}(), A) - @test cover_min!(AbsLog{2}(), cover(Aasym)..., Aasym)[1] ≈ ca + @test first(symcover_min!(AbsLog{2}(), symcover(A), A)) ≈ first(symcover_min(AbsLog{2}(), A)) + @test scales(cover_min!(AbsLog{2}(), cover(Aasym)..., Aasym))[1] ≈ ca # Every start is read only up to the gauge a -> c*a, b -> b/c, which leaves # every product a[i]*b[j] fixed. @@ -781,7 +794,7 @@ end Anosupp = [4.0 1.0 0.0; 1.0 9.0 0.0; 0.0 0.0 0.0] a = initialize_symcover(Anosupp) a[3] = 7.0 - @test symcover_min!(AbsLog{2}(), a, Anosupp)[3] == 0.0 + @test first(symcover_min!(AbsLog{2}(), a, Anosupp))[3] == 0.0 # A start must cover A, and be positive wherever A carries support. @test_throws "requires a start that covers `A`" symcover_min!(AbsLog{2}(), fill(0.1, 3), A) @@ -795,9 +808,9 @@ end # Offset axes propagate through the start and the result. Ao = OffsetArray(A, -1, -1) ao = initialize_symcover(Ao) - @test symcover_min!(AbsLog{2}(), ao, Ao) === ao + @test first(symcover_min!(AbsLog{2}(), ao, Ao)) === ao @test axes(ao, 1) == axes(Ao, 1) - @test collect(ao) ≈ symcover_min(AbsLog{2}(), A) + @test collect(ao) ≈ first(symcover_min(AbsLog{2}(), A)) end @testset "quality vs optimal (testmatrices)" begin @@ -810,7 +823,7 @@ end a0 = symcover(AbsLog{2}(), Af / 100; maxiter=0) @test iscover(a0, Af / 100; atol=1e-12) # Covers are nearly quadratically optimal - qopt = cover_objective(AbsLog{2}(), symcover_min(AbsLog{2}(), Af), Af) + qopt = cover_objective(AbsLog{2}(), first(symcover_min(AbsLog{2}(), Af)), Af) qfast = cover_objective(AbsLog{2}(), symcover(AbsLog{2}(), Af; maxiter=10), Af) iszero(qopt) || push!(sym_ratios, qfast / qopt) end @@ -823,7 +836,7 @@ end @test iscover(a0, b0, Af; atol=1e-12) a0, b0 = cover(AbsLog{2}(), Af / 100; maxiter=0) @test iscover(a0, b0, Af / 100; atol=1e-12) - qopt = cover_objective(AbsLog{2}(), cover_min(AbsLog{2}(), Af)..., Af) + qopt = cover_objective(AbsLog{2}(), scales(cover_min(AbsLog{2}(), Af))..., Af) qfast = cover_objective(AbsLog{2}(), cover(AbsLog{2}(), Af; maxiter=10)..., Af) iszero(qopt) || push!(gen_ratios, qfast / qopt) end @@ -918,13 +931,19 @@ end # With the penalty held at the cap they are zeroed directly; left to drain by # `2(κ-1)z` per update they hold the KKT residual at 2.6e-6. # Without the active-set finish (`polish=false`), slack multipliers are zeroed. + # With the diagonal preconditioner the residual then stalls above the + # tolerance, which warns; the default solve's finish certifies the result. A = banded_sparse_sym(StableRNG(5), 10_000, 3, 1.0) for fillbudget in (MatrixCovers.LSQR_FILL_BUDGET, 0) - a, s = @test_nowarn MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, fillbudget, polish=false) + solve = () -> MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, fillbudget, polish=false) + a, s = fillbudget == 0 ? (@test_logs (:warn, r"try a larger `κ`") solve()) : (@test_nowarn solve()) + @test s.converged == (fillbudget != 0) @test s.precond === (fillbudget == 0 ? :diagonal : :factor) @test s.nzeroed > 0 @test maximum(s.κs) == 1e5 @test s.kkt[end] < 1e-9 + _, sp = @test_logs symcover_min(AbsLog{2}(), A; fillbudget) + @test MatrixCovers.converged(sp) @test cover_objective(AbsLog{2}(), a, A) ≈ 161350.449205 rtol=1e-9 @test iscover(a, A) end @@ -968,9 +987,11 @@ end @test cover_objective(AbsLog{2}(), a1, b1, A) ≈ cover_objective(AbsLog{2}(), a0, b0, A) rtol=1e-10 @test a1 .* b1' ≈ a0 .* b0' rtol=1e-7 @test iscover(a1, b1, A) - # Without the polish the iteration keeps zeroing. - _, _, s2 = MatrixCovers._cover_min_abslog2(A; polish=false) + # Without the polish the iteration keeps zeroing, and stalls above the + # tolerance. + _, _, s2 = @test_logs (:warn, r"try a larger `κ`") MatrixCovers._cover_min_abslog2(A; polish=false) @test s2.nzeroed > 0 + @test !s2.converged end @testset "MMC multiplier update is exact on an analytic problem" begin @@ -995,7 +1016,7 @@ end à = A ./ (d * d') @test à .* (d * d') == A @test symcover(Ã) .* d ≈ symcover(A) rtol=1e-12 - @test symcover_min(Ã) .* d ≈ symcover_min(A) rtol=1e-6 + @test first(symcover_min(Ã)) .* d ≈ first(symcover_min(A)) rtol=1e-6 end for seed in 1:3 rng = StableRNG(seed) @@ -1006,9 +1027,9 @@ end e = 2.0 .^ round.(Int, 4 * log2(10) .* randn(rng, n)) à = A ./ (d * d') @test symcover(Ã) .* d ≈ symcover(A) rtol=1e-12 - @test symcover_min(Ã) .* d ≈ symcover_min(A) rtol=1e-6 + @test first(symcover_min(Ã)) .* d ≈ first(symcover_min(A)) rtol=1e-6 a, b = cover_min(A) - ã, b̃ = cover_min(A ./ (d * e')) + ã, b̃ = scales(cover_min(A ./ (d * e'))) @test (ã * b̃') .* (d * e') ≈ a * b' rtol=1e-6 end end @@ -1020,7 +1041,7 @@ end a, b = cover_min(A) @test iscover(a, b, A) @test all(i -> isapprox(a[i+1] / a[i], sqrt(1 / 2); rtol=1e-6), 50:150) - @test maximum(abs, diff(log.(cover_min(toep(200, 1.0, 7.0, 1.0))[1]))) < 1e-6 + @test maximum(abs, diff(log.(scales(cover_min(toep(200, 1.0, 7.0, 1.0)))[1]))) < 1e-6 # Here the factors span about 690 decades, beyond Float64's range. B = toep(70, 1.0, 1.0, 1e20) diff --git a/test/powermean.jl b/test/powermean.jl index 6b0d006..3a043f9 100644 --- a/test/powermean.jl +++ b/test/powermean.jl @@ -63,23 +63,23 @@ end sp.nzval .= exp.(3 .* randn(rng, nnz(sp))) sp[7, :] .= 0; sp[:, 9] .= 0; dropzeros!(sp) for A in (dense, rect, wide, sp) - a, b = @test_logs soft_cover(ϕ, A) + a, b = @test_logs scales(soft_cover(ϕ, A)) @test pm_imbalance(p, a, b, A) < 1e-11 @test isbalanced(a, b, A) @test all(iszero(a[i]) == all(iszero, A[i, :]) for i in axes(A, 1)) @test all(iszero(b[j]) == all(iszero, A[:, j]) for j in axes(A, 2)) end - @test soft_cover(ϕ, sp)[1] isa Vector{Float64} + @test scales(soft_cover(ϕ, sp))[1] isa Vector{Float64} S = pm_testmatrix(rng, 35, 35; σ=3.0, zrows=(4,), zcols=(4,)) S = S + S' Ssp = sparse(S) for M in (S, Ssp, Symmetric(Ssp)) - a = soft_symcover(ϕ, M) + a = scales(soft_symcover(ϕ, M)) @test pm_imbalance(p, a, a, M) < 1e-11 @test a[4] == 0 # The unique minimizer is symmetric, so the asymmetric solver agrees. - x, y = soft_cover(ϕ, M) + x, y = scales(soft_cover(ϕ, M)) supp = findall(!iszero, S) @test (x .* y')[supp] ≈ (a .* a')[supp] rtol=1e-10 end diff --git a/test/soft_covers.jl b/test/soft_covers.jl index 5c30138..790d4fc 100644 --- a/test/soft_covers.jl +++ b/test/soft_covers.jl @@ -10,7 +10,7 @@ end # At the minimum, gradients should be near zero for A in ([2.0 1.0; 1.0 3.0], [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0]) # At the minimum, ∑_j (log a[k] + log a[j] - log|A[k,j]|) ≈ 0 for each k - a = soft_symcover(AbsLog{2}(), A) + a, _ = soft_symcover(AbsLog{2}(), A) for k in axes(A, 1) residual = sum(abs(A[k,j]) != 0 ? log(a[k]) + log(a[j]) - log(abs(A[k,j])) : 0.0 for j in axes(A, 2)) @@ -33,7 +33,7 @@ end # Rank-1 is exact: A = [4 2; 2 1] = [2;1]*[2 1], so exact cover has all ratios = 1 A = [4.0 2.0; 2.0 1.0] for ϕ in SOFT_PENALTIES - a = soft_symcover(ϕ, A) + a = scales(soft_symcover(ϕ, A)) @test cover_objective(ϕ, a, A) ≈ 0.0 atol=1e-8 end @@ -41,13 +41,13 @@ end A = [2.0 1.0; 1.0 3.0] d = [2.0, 0.5] for ϕ in SOFT_PENALTIES - @test covaries(a -> soft_symcover(ϕ, a), A, d; rtol=1e-5) + @test covaries(a -> scales(soft_symcover(ϕ, a)), A, d; rtol=1e-5) end # Lower objective than naive uniform scaling perturbations A = [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0] for ϕ in SOFT_PENALTIES - a = soft_symcover(ϕ, A) + a = scales(soft_symcover(ϕ, A)) E = cover_objective(ϕ, a, A) for scale in [0.5, 0.9, 1.1, 2.0] @test E <= cover_objective(ϕ, scale .* a, A) + 1e-8 @@ -59,8 +59,8 @@ end A_zero = [γ 1.0; 1.0 0.0] A_small = [γ 1.0; 1.0 1e-10] for ϕ in (AbsLinear{1}(), AbsLinear{2}()) - a_zero = soft_symcover(ϕ, A_zero) - a_small = soft_symcover(ϕ, A_small) + a_zero = scales(soft_symcover(ϕ, A_zero)) + a_small = scales(soft_symcover(ϕ, A_small)) @test a_small ≈ a_zero atol=1e-5 end @@ -68,8 +68,8 @@ end for ϕ in (AbsLinear{1}(), AbsLinear{2}()) for (B, d) in ((A_small, [50.0, 0.02]), ([3.0 7.6e-10; 7.6e-10 80.0], [35.0, 3400.0])) - a0 = soft_symcover(ϕ, B) - as = soft_symcover(ϕ, B .* d .* d') + a0 = scales(soft_symcover(ϕ, B)) + as = scales(soft_symcover(ϕ, B .* d .* d')) @test as .* as' ≈ (d .* d') .* (a0 .* a0') rtol=1e-6 end end @@ -177,7 +177,7 @@ end S = [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] a = soft_symcover(AbsLog{1}(), S) E = cover_objective(AbsLog{1}(), a, S) - @test E <= cover_objective(AbsLog{1}(), soft_symcover(AbsLog{2}(), S), S) + 1e-12 + @test E <= cover_objective(AbsLog{1}(), first(soft_symcover(AbsLog{2}(), S)), S) + 1e-12 for _ in 1:20 @test E <= cover_objective(AbsLog{1}(), a .* exp.(0.1 .* randn(rng, 4)), S) + 1e-8 end @@ -406,7 +406,7 @@ end asym_dense = rand(rng, 4, 6) for A in (sym_zeros, sym_dense, float.(last(symmetric_matrices[1]))) - a = soft_symcover(AbsLog{2}(), A) + a, _ = soft_symcover(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) @@ -423,34 +423,34 @@ end # Sparse support distinguishes the exact minimum from the geometric mean. @test cover_objective(AbsLog{2}(), initialize_symcover(sym_dense; strategy=:geomean, feasible=:none), sym_dense) ≈ - cover_objective(AbsLog{2}(), soft_symcover(AbsLog{2}(), sym_dense), sym_dense) rtol=1e-8 + cover_objective(AbsLog{2}(), first(soft_symcover(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(AbsLog{2}(), sym_zeros), sym_zeros) * (1 + 1e-6) + cover_objective(AbsLog{2}(), first(soft_symcover(AbsLog{2}(), sym_zeros)), sym_zeros) * (1 + 1e-6) # Convex with one minimizer, so a refiner cannot be steered by its start. for strategy in (:geomean, :leaveout, :diagfeasible, :hardcover) a0 = initialize_symcover(sym_zeros; strategy, feasible=:none) - @test soft_symcover!(AbsLog{2}(), a0, sym_zeros) ≈ soft_symcover(AbsLog{2}(), sym_zeros) + @test first(soft_symcover!(AbsLog{2}(), a0, sym_zeros)) ≈ first(soft_symcover(AbsLog{2}(), sym_zeros)) end a0, b0 = initialize_cover(asym_zeros; strategy=:geomean, feasible=:none) ar, br = soft_cover!(AbsLog{2}(), a0, b0, asym_zeros) - @test (ar, br) == soft_cover(AbsLog{2}(), asym_zeros) + @test (ar, br) == scales(soft_cover(AbsLog{2}(), asym_zeros)) # A bipartite support graph makes the signless Laplacian singular; the balanced # representative is the one reported. - @test soft_symcover(AbsLog{2}(), [0 1; 1 0]) ≈ [1.0, 1.0] + @test first(soft_symcover(AbsLog{2}(), [0 1; 1 0])) ≈ [1.0, 1.0] # The `:lsqr` and dense paths solve the same problem. - @test soft_symcover(AbsLog{2}(), sym_zeros; linsolve=:lsqr) ≈ - soft_symcover(AbsLog{2}(), sym_zeros; linsolve=:dense) rtol=1e-6 + @test first(soft_symcover(AbsLog{2}(), sym_zeros; linsolve=:lsqr)) ≈ + first(soft_symcover(AbsLog{2}(), sym_zeros; linsolve=:dense)) rtol=1e-6 # The soft cover is the hard worker's cold solve, so it reaches the Woodbury path # too: one sparse Cholesky in place of the dense factorization, same answer. rng = StableRNG(77) X = exp.(randn(rng, 40, 40)) Adense = (X .+ X') ./ 2 - @test soft_symcover(AbsLog{2}(), Adense; linsolve=:woodbury) ≈ - soft_symcover(AbsLog{2}(), Adense; linsolve=:dense) rtol=1e-8 + @test first(soft_symcover(AbsLog{2}(), Adense; linsolve=:woodbury)) ≈ + first(soft_symcover(AbsLog{2}(), Adense; linsolve=:dense)) rtol=1e-8 @test MatrixCovers._soft_symcover_abslog2(Adense)[2].linsolve === :woodbury Y = exp.(randn(rng, 40, 30)) aw, bw = soft_cover(AbsLog{2}(), Y; linsolve=:woodbury) @@ -463,20 +463,20 @@ end rng = StableRNG(19) S0 = sprand(rng, 60, 60, 0.1) Ssp = S0 + S0' + I # supported diagonal keeps every component non-bipartite - a = soft_symcover(AbsLog{2}(), Ssp) + a, _ = soft_symcover(AbsLog{2}(), Ssp) @test a isa Vector{Float64} - @test a == soft_symcover(AbsLog{2}(), Ssp; linsolve=:lsqr) - @test a ≈ soft_symcover(AbsLog{2}(), Matrix(Ssp); linsolve=:dense) rtol=1e-6 + @test a == first(soft_symcover(AbsLog{2}(), Ssp; linsolve=:lsqr)) + @test a ≈ first(soft_symcover(AbsLog{2}(), Matrix(Ssp); linsolve=:dense)) rtol=1e-6 U = Symmetric(sparse(triu(Ssp))) - @test soft_symcover(AbsLog{2}(), U) == soft_symcover(AbsLog{2}(), U; linsolve=:lsqr) + @test first(soft_symcover(AbsLog{2}(), U)) == first(soft_symcover(AbsLog{2}(), U; linsolve=:lsqr)) # Asymmetric form; the cover products are gauge-free. G = sprand(rng, 50, 40, 0.12) ga, gb = soft_cover(AbsLog{2}(), G) - @test (ga, gb) == soft_cover(AbsLog{2}(), G; linsolve=:lsqr) + @test (ga, gb) == scales(soft_cover(AbsLog{2}(), G; linsolve=:lsqr)) gad, gbd = soft_cover(AbsLog{2}(), Matrix(G); linsolve=:dense) @test ga .* gb' ≈ gad .* gbd' rtol=1e-6 a0 = copy(a) - @test soft_symcover!(AbsLog{2}(), a0, Ssp) == a + @test first(soft_symcover!(AbsLog{2}(), a0, Ssp)) == a end @testset "soft_symcover!/soft_cover! refiners" begin @@ -485,8 +485,8 @@ end @testset "refining the multistart's own start reproduces it: $ϕ" for ϕ in SOFT_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 + scales(soft_symcover!(ϕ, a, Asym)) + @test cover_objective(ϕ, a, Asym) ≈ cover_objective(ϕ, scales(soft_symcover(ϕ, Asym)), Asym) rtol=1e-6 end @testset "no-ϕ form defaults to PowerMean{2}" begin @@ -505,8 +505,8 @@ end 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)) + @test o(AbsLog{2}(), first(soft_symcover!(AbsLog{2}(), s1, Abasins))) ≈ + o(AbsLog{2}(), first(soft_symcover!(AbsLog{2}(), s2, Abasins))) s1, s2 = starts() @test !isapprox(o(AbsLinear{2}(), soft_symcover!(AbsLinear{2}(), s1, Abasins)), @@ -517,10 +517,10 @@ end @testset "start need not cover A: $ϕ" for ϕ in SOFT_PENALTIES a = fill(0.01, 3) @test !iscover(a, Asym) - @test soft_symcover!(ϕ, a, Asym) === a + @test scales(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) + @test scales(soft_cover!(ϕ, b, c, Agen)) === (b, c) end @testset "invalid starts throw, naming the function called" begin @@ -539,24 +539,24 @@ end @testset "unsupported rows are zeroed: $ϕ" for ϕ in SOFT_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 + @test scales(soft_symcover!(ϕ, a, Az))[2] == 0 end @testset "asymmetric refiners pin the balance convention: $ϕ" for ϕ in SOFT_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) + scales(soft_cover!(ϕ, a, b, Agen)) @test isbalanced(a, b, Agen) end @testset "offset axes propagate: $ϕ" for ϕ in SOFT_PENALTIES Ao = OffsetArray(Asym, -1:1, -1:1) ao = initialize_symcover(Ao; strategy=:geomean, feasible=:none) - soft_symcover!(ϕ, ao, Ao) + scales(soft_symcover!(ϕ, ao, Ao)) @test axes(ao, 1) == axes(Ao, 1) aref = initialize_symcover(Asym; strategy=:geomean, feasible=:none) - soft_symcover!(ϕ, aref, Asym) + scales(soft_symcover!(ϕ, aref, Asym)) @test collect(ao) ≈ aref rtol=1e-6 end end @@ -607,7 +607,7 @@ MatrixCovers.require_abs_symmetric(::HookOnlyMatrix, fname) = nothing @test soft_symcover(PowerMean{2}(), M) ≈ soft_symcover(PowerMean{2}(), dense) rtol=1e-12 for ϕ in (AbsLinear{1}(), AbsLinear{2}(), AbsLog{1}(), AbsLog{2}()) @test symcover(ϕ, M) ≈ symcover(ϕ, dense) rtol=1e-10 - @test soft_symcover(ϕ, M) ≈ soft_symcover(ϕ, dense) rtol=1e-6 + @test scales(soft_symcover(ϕ, M)) ≈ scales(soft_symcover(ϕ, dense)) rtol=1e-6 @test iscover(symcover(ϕ, M), dense; rtol=8eps()) @test cover_objective(ϕ, symcover(ϕ, M), M) ≈ cover_objective(ϕ, symcover(ϕ, dense), dense) rtol=1e-10 diff --git a/test/storage_types.jl b/test/storage_types.jl index d309e04..ff08a26 100644 --- a/test/storage_types.jl +++ b/test/storage_types.jl @@ -17,7 +17,7 @@ MatrixCovers.foreach_support_sym((i, j, v) -> push!(ref, (i, j, v)), Matrix(H)) @test sort(got) == sort(ref) - a = symcover_min(AbsLog{2}(), H) + a, _ = symcover_min(AbsLog{2}(), H) @test iscover(a, Matrix(H); rtol=1e-8) @test symcover(AbsLog{2}(), H) ≈ symcover(AbsLog{2}(), Matrix(H)) rtol = 1e-12 end @@ -193,11 +193,11 @@ end symdenses = [[2.0 1.0 0.0; 1.0 3.0 2.0; 0.0 2.0 5.0], [4.0 0.0 1.0; 0.0 0.0 0.0; 1.0 0.0 2.0]] # second has a zero row/column for M in symdenses - ad = symcover_min(AbsLog{2}(), M) + ad, _ = symcover_min(AbsLog{2}(), M) asd = soft_symcover(AbsLinear{2}(), M) for A in (sparse(M), Symmetric(sparse(triu(M)), :U), Symmetric(sparse(tril(M)), :L), Hermitian(sparse(triu(M)), :U)) - a = symcover_min(AbsLog{2}(), A) + a, _ = symcover_min(AbsLog{2}(), A) @test iscover(a, M; atol=1e-7) @test cover_objective(AbsLog{2}(), a, M) ≈ cover_objective(AbsLog{2}(), ad, M) rtol = 1e-7 atol = 1e-10 as = soft_symcover(AbsLinear{2}(), A) @@ -206,12 +206,12 @@ end @test as isa Vector{Float64} end end - @test symcover_min(AbsLog{2}(), sparse(symdenses[1])) isa Vector{Float64} + @test first(symcover_min(AbsLog{2}(), sparse(symdenses[1]))) isa Vector{Float64} # Every soft_symcover penalty returns a dense Vector on sparse-backed input. let Ssp = sparse(symdenses[1]) for ϕ in SOFT_PENALTIES - @test soft_symcover(ϕ, Ssp) isa Vector{Float64} - @test soft_symcover(ϕ, Symmetric(Ssp)) isa Vector{Float64} + @test scales(soft_symcover(ϕ, Ssp)) isa Vector{Float64} + @test scales(soft_symcover(ϕ, Symmetric(Ssp))) isa Vector{Float64} end end @@ -228,16 +228,16 @@ end @test cover_objective(AbsLinear{2}(), as, bs, M) ≈ cover_objective(AbsLinear{2}(), asd, bsd, M) rtol = 1e-7 atol = 1e-10 @test as isa Vector{Float64} && bs isa Vector{Float64} end - let (a, b) = cover_min(AbsLog{2}(), sparse(gendenses[1])) + let (a, b) = scales(cover_min(AbsLog{2}(), sparse(gendenses[1]))) @test a isa Vector{Float64} && b isa Vector{Float64} end for D in (Diagonal([4.0, 9.0, 1.0]), Diagonal([4.0, 0.0, 1.0])) M = Matrix(D) - a = symcover_min(AbsLog{2}(), D) + a, _ = symcover_min(AbsLog{2}(), D) @test iscover(a, M; atol=1e-7) @test cover_objective(AbsLog{2}(), a, M) ≈ - cover_objective(AbsLog{2}(), symcover_min(AbsLog{2}(), M), M) rtol = 1e-7 + cover_objective(AbsLog{2}(), first(symcover_min(AbsLog{2}(), M)), M) rtol = 1e-7 a2, b2 = cover_min(AbsLog{2}(), D) # disconnected support: exercises the gauge ridge @test iscover(a2, b2, M; atol=1e-7) @test cover_objective(AbsLog{2}(), a2, b2, M) ≈ 0.0 atol = 1e-8 # diagonal is exactly coverable @@ -248,7 +248,7 @@ end a2, b2 = cover_min(AbsLog{2}(), A) @test iscover(a2, b2, M; atol=1e-7) @test cover_objective(AbsLog{2}(), a2, b2, M) ≈ - cover_objective(AbsLog{2}(), cover_min(AbsLog{2}(), M)..., M) rtol = 1e-7 + cover_objective(AbsLog{2}(), scales(cover_min(AbsLog{2}(), M))..., M) rtol = 1e-7 end end diff --git a/test/support.jl b/test/support.jl index 8048c08..53311b9 100644 --- a/test/support.jl +++ b/test/support.jl @@ -121,7 +121,7 @@ end # Hermitian — where A[i,j] == conj(A[j,i]) — is a legitimate input. @test iscover(symcover([4.0 1.0; -1.0 4.0]), [4.0 1.0; -1.0 4.0]; rtol=8eps()) P = sparse([1, 2, 1], [1, 2, 2], ComplexF64[4.0, 2.0, 0.5+0.5im], 2, 2) - @test symcover_min(AbsLog{2}(), Hermitian(P, :U)) isa AbstractVector + @test first(symcover_min(AbsLog{2}(), Hermitian(P, :U))) isa AbstractVector # Exact symmetry is not required: a symmetric matrix that has been through # floating-point arithmetic lands a ULP or so off, and must still be accepted. diff --git a/test/transversal.jl b/test/transversal.jl index 4be0682..1b86c69 100644 --- a/test/transversal.jl +++ b/test/transversal.jl @@ -72,7 +72,7 @@ end (aj, bj), lπj = MatrixCovers.cover_transversal_jump(Af) @test lπj ≈ lπ rtol=1e-8 @test logproduct_gap(a, b, aj, bj, Af) <= 1e-5 - @test cover_objective(AbsLog{2}(), a, b, Af) >= cover_objective(AbsLog{2}(), cover_min(Af)..., Af) * (1 - 1e-6) - 1e-10 + @test cover_objective(AbsLog{2}(), a, b, Af) >= cover_objective(AbsLog{2}(), scales(cover_min(Af))..., Af) * (1 - 1e-6) - 1e-10 end # Independence from the transversal used: solve with each maximum-product @@ -99,7 +99,7 @@ end # Coincides with the minimal cover when that is already tight on the transversal. D = [4.0 1 1; 1 4 1; 1 1 4] - @test logproduct_gap(cover_transversal(D)..., cover_min(D)..., D) <= 1e-7 + @test logproduct_gap(scales(cover_transversal(D))..., scales(cover_min(D))..., D) <= 1e-7 # Scale covariance, including power-of-two rescalings, which change the # floating-point costs seen by the matching. @@ -116,7 +116,7 @@ end d = 2.0 .^ rand(rng, -30:30, n) e = 2.0 .^ rand(rng, -30:30, n) a, b = cover_transversal(A) - ã, b̃ = cover_transversal(Diagonal(d) * A * Diagonal(e)) + ã, b̃ = scales(cover_transversal(Diagonal(d) * A * Diagonal(e))) @test logproduct_gap(ã ./ d, b̃ ./ e, a, b, A) <= 1e-6 # The sparse default (`:lsqr`) and the dense solve agree. ad, bd = cover_transversal(Matrix(A); linsolve=:dense) @@ -130,7 +130,7 @@ end T3 = sparse([4.0 1.0 0.0; 1.0 4.0 1.0; 0.0 1.0 4.0]) a, b, st = _cover_min_abslog2(T3; transversal=true, linsolve=:lsqr, fname=:cover_transversal) @test st.precond === :factor - @test logproduct_gap(a, b, cover_transversal(Matrix(T3); linsolve=:dense)..., T3) <= 1e-7 + @test logproduct_gap(a, b, scales(cover_transversal(Matrix(T3); linsolve=:dense))..., T3) <= 1e-7 # Several support components, each balanced separately. Ablk = [1.0 2.0 0 0; 0.25 3.0 0 0; 0 0 1.5 2.5; 0 0 0.75 0] @@ -161,7 +161,7 @@ end # Mutating form writes through its buffers. a2, b2 = zeros(3), zeros(3) - @test cover_transversal!(a2, b2, A) === (a2, b2) + @test scales(cover_transversal!(a2, b2, A)) === (a2, b2) @test (a2, b2) == (a, b) @test_throws DimensionMismatch cover_transversal!(zeros(2), zeros(3), A) diff --git a/test/unitful.jl b/test/unitful.jl index 8a43cd4..3209cba 100644 --- a/test/unitful.jl +++ b/test/unitful.jl @@ -93,19 +93,19 @@ @test all(unit.(v) == fill(u"m^-1", 3) for v in cover(S)) @test unit.(symcover(Symmetric(S))) == fill(u"m^-1", 3) - @test unit.(symcover_min(AbsLog{2}(), S)) == fill(u"m^-1", 3) - @test unit.(symcover_min(AbsLog{2}(), Symmetric(S))) == fill(u"m^-1", 3) - @test all(unit.(v) == fill(u"m^-1", 3) for v in cover_min(AbsLog{2}(), S)) - @test unit.(soft_symcover(AbsLog{2}(), S)) == fill(u"m^-1", 3) - @test unit.(soft_symcover(AbsLog{2}(), Symmetric(S))) == fill(u"m^-1", 3) - @test all(unit.(v) == fill(u"m^-1", 3) for v in soft_cover(AbsLog{2}(), S)) + @test unit.(first(symcover_min(AbsLog{2}(), S))) == fill(u"m^-1", 3) + @test unit.(first(symcover_min(AbsLog{2}(), Symmetric(S)))) == fill(u"m^-1", 3) + @test all(unit.(v) == fill(u"m^-1", 3) for v in scales(cover_min(AbsLog{2}(), S))) + @test unit.(first(soft_symcover(AbsLog{2}(), S))) == fill(u"m^-1", 3) + @test unit.(first(soft_symcover(AbsLog{2}(), Symmetric(S)))) == fill(u"m^-1", 3) + @test all(unit.(v) == fill(u"m^-1", 3) for v in scales(soft_cover(AbsLog{2}(), S))) a = initialize_symcover(S) - @test unit.(symcover_min!(AbsLog{2}(), a, S)) == fill(u"m^-1", 3) - @test unit.(symcover_min!(AbsLog{2}(), a, Symmetric(S))) == fill(u"m^-1", 3) + @test unit.(first(symcover_min!(AbsLog{2}(), a, S))) == fill(u"m^-1", 3) + @test unit.(first(symcover_min!(AbsLog{2}(), a, Symmetric(S)))) == fill(u"m^-1", 3) a, b = initialize_cover(S) - @test all(unit.(v) == fill(u"m^-1", 3) for v in cover_min!(AbsLog{2}(), a, b, S)) - @test all(unit.(v) == fill(u"m^-1", 3) for v in cover_transversal(S)) - @test all(unit.(v) == fill(u"m^-1", 3) for v in cover_transversal!(a, b, S)) + @test all(unit.(v) == fill(u"m^-1", 3) for v in scales(cover_min!(AbsLog{2}(), a, b, S))) + @test all(unit.(v) == fill(u"m^-1", 3) for v in scales(cover_transversal(S))) + @test all(unit.(v) == fill(u"m^-1", 3) for v in scales(cover_transversal!(a, b, S))) end @testset "balance convention holds in the caller's units" begin @@ -164,7 +164,7 @@ @testset "the whole family" begin for f in (symcover, initialize_symcover, soft_symcover, symcover_min) - a = f(A) + a = scales(f(A)) @test unit.(a) == UA @test eltype(ustrip.(a)) <: AbstractFloat end @@ -174,15 +174,15 @@ end for ϕ in PENALTIES @test unit.(symcover(ϕ, A)) == UA - @test unit.(soft_symcover(ϕ, A)) == UA - @test unit.(symcover_min(ϕ, A)) == UA + @test unit.(scales(soft_symcover(ϕ, A))) == UA + @test unit.(scales(symcover_min(ϕ, A))) == UA @test all(unit.(v) == UA for v in cover(ϕ, A)) - @test all(unit.(v) == UA for v in soft_cover(ϕ, A)) - @test all(unit.(v) == UA for v in cover_min(ϕ, A)) + @test all(unit.(v) == UA for v in scales(soft_cover(ϕ, A))) + @test all(unit.(v) == UA for v in scales(cover_min(ϕ, A))) end for ϕ in (PowerMean{1}(), PowerMean{2}()) - @test unit.(soft_symcover(ϕ, A)) == UA - @test all(unit.(v) == UA for v in soft_cover(ϕ, A)) + @test unit.(scales(soft_symcover(ϕ, A))) == UA + @test all(unit.(v) == UA for v in scales(soft_cover(ϕ, A))) end end @@ -219,23 +219,23 @@ # Initializer output is valid refiner input with matching units. a = initialize_symcover(A) - @test symcover_min!(AbsLog{2}(), a, A) === a + @test first(symcover_min!(AbsLog{2}(), a, A)) === a @test unit.(a) == UA @test iscover(a, A; rtol=1e-8) a, b = initialize_cover(A) - @test cover_min!(AbsLog{2}(), a, b, A) == (a, b) + @test scales(cover_min!(AbsLog{2}(), a, b, A)) == (a, b) @test unit.(a) == unit.(b) == UA # The no-ϕ forms of the `*_min` family default their penalty, and carry # the units the same way. a = initialize_symcover(A) - @test symcover_min!(a, A) === a + @test first(symcover_min!(a, A)) === a @test unit.(a) == UA @test iscover(a, A; rtol=1e-8) a, b = initialize_cover(A) - @test cover_min!(a, b, A) == (a, b) + @test scales(cover_min!(a, b, A)) == (a, b) @test unit.(a) == unit.(b) == UA a = initialize_symcover(A) @@ -251,7 +251,7 @@ @test unit.(a) == UA a, b = initialize_cover(A) - @test soft_cover!(AbsLog{2}(), a, b, A) == (a, b) + @test scales(soft_cover!(AbsLog{2}(), a, b, A)) == (a, b) @test unit.(a) == unit.(b) == UA a = initialize_symcover(A) @@ -265,7 +265,7 @@ # A start in a dimensionally equivalent spelling is converted, not rejected. a = initialize_symcover(A) a2 = [uconvert(u, x) for (u, x) in zip((u"mm^-1", u"hr/km", u"kN^-1"), a)] - @test ustrip.(symcover_min!(AbsLog{2}(), a2, A)) ≈ ustrip.(symcover_min!(AbsLog{2}(), a, A)) + @test ustrip.(first(symcover_min!(AbsLog{2}(), a2, A))) ≈ ustrip.(first(symcover_min!(AbsLog{2}(), a, A))) # A start whose dimensions are wrong is not. a3 = Quantity{Float64}[1.0u"m", 1.0u"s", 1.0u"N"] From d3f7d4806c71be740247031630e25c3950bdc6a1 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 9 Oct 2026 05:19:52 -0500 Subject: [PATCH 3/7] Return solver statistics from PowerMean covers `soft_symcover`, `soft_cover`, and their mutating forms with a `PowerMean` penalty, including the no-penalty forms that default to `PowerMean{2}()`, now return a `PowerMeanStats` object as their last output: `a, stats = soft_symcover(A)`, `a, b, stats = soft_cover(A)`. `PowerMeanStats{T}` is a public (unexported) subtype of `AbstractCoverStats`. `converged`, `residual`, and `tolerance` report whether the largest row (and column) imbalance of the returned cover meets `tol`; `iterations` counts power-mean updates plus Newton steps. The element type is the working type of the solve, so the result is inferrable. The warning on an unconverged solve is unchanged. The public methods declare their keywords (`tol`, `maxiter`, `newton`, `maxnewton`, `linsolve`) instead of forwarding them, with the same defaults. The private solvers and the statistics tuples they return are unchanged. Assisted-by: Claude Opus 5.5 --- src/MatrixCovers.jl | 2 +- src/powermean.jl | 47 ++++++++++++------- src/stats.jl | 53 ++++++++++++++++++++++ test/element_types.jl | 4 +- test/invariants.jl | 8 ++-- test/powermean.jl | 103 ++++++++++++++++++++++++++++-------------- test/soft_covers.jl | 10 ++-- test/unitful.jl | 8 ++-- 8 files changed, 170 insertions(+), 65 deletions(-) diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index a07dec4..946e073 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -18,7 +18,7 @@ export cover_transversal, cover_transversal! # `public` is parsed as a keyword only from Julia 1.11; this package supports 1.10. @static if VERSION >= v"1.11" - eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym, scalar_type, SupportComponents, support_components, ncomponents, rowcomponent, colcomponent, SolverFailure, AbstractCoverStats, AugmentedLagrangianStats, ActiveSetStats, LeastSquaresStats, converged, iterations, residual, tolerance")) + eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym, scalar_type, SupportComponents, support_components, ncomponents, rowcomponent, colcomponent, SolverFailure, AbstractCoverStats, AugmentedLagrangianStats, ActiveSetStats, LeastSquaresStats, PowerMeanStats, converged, iterations, residual, tolerance")) end include("penalties.jl") diff --git a/src/powermean.jl b/src/powermean.jl index d9c0efc..78fc82d 100644 --- a/src/powermean.jl +++ b/src/powermean.jl @@ -43,15 +43,29 @@ _powermean_type(Ts::Type...) = float(promote_type(Ts..., Float64)) # Public methods # ============================================================ -soft_cover(ϕ::PowerMean, A::AbstractMatrix; kwargs...) = - _soft_cover_powermean(ϕ, A, :soft_cover; kwargs...) -soft_cover!(ϕ::PowerMean, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) = - _soft_cover_powermean!(ϕ, a, b, A, :soft_cover!; kwargs...) +function soft_cover(ϕ::PowerMean, A::AbstractMatrix; tol::Union{Real,Nothing}=nothing, + maxiter::Integer=POWERMEAN_MAXITER, newton::Bool=true, + maxnewton::Integer=POWERMEAN_MAXNEWTON, linsolve::Symbol=:auto) + return _soft_cover_powermean(ϕ, A, :soft_cover; tol, maxiter, newton, maxnewton, linsolve) +end + +function soft_cover!(ϕ::PowerMean, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; + tol::Union{Real,Nothing}=nothing, maxiter::Integer=POWERMEAN_MAXITER, + newton::Bool=true, maxnewton::Integer=POWERMEAN_MAXNEWTON, linsolve::Symbol=:auto) + return _soft_cover_powermean!(ϕ, a, b, A, :soft_cover!; tol, maxiter, newton, maxnewton, linsolve) +end -soft_symcover(ϕ::PowerMean, A::AbstractMatrix; kwargs...) = - _soft_symcover_powermean(ϕ, A, :soft_symcover; kwargs...) -soft_symcover!(ϕ::PowerMean, a::AbstractVector, A::AbstractMatrix; kwargs...) = - _soft_symcover_powermean!(ϕ, a, A, :soft_symcover!; kwargs...) +function soft_symcover(ϕ::PowerMean, A::AbstractMatrix; tol::Union{Real,Nothing}=nothing, + maxiter::Integer=POWERMEAN_MAXITER, newton::Bool=true, + maxnewton::Integer=POWERMEAN_MAXNEWTON, linsolve::Symbol=:auto) + return _soft_symcover_powermean(ϕ, A, :soft_symcover; tol, maxiter, newton, maxnewton, linsolve) +end + +function soft_symcover!(ϕ::PowerMean, a::AbstractVector, A::AbstractMatrix; + tol::Union{Real,Nothing}=nothing, maxiter::Integer=POWERMEAN_MAXITER, + newton::Bool=true, maxnewton::Integer=POWERMEAN_MAXNEWTON, linsolve::Symbol=:auto) + return _soft_symcover_powermean!(ϕ, a, A, :soft_symcover!; tol, maxiter, newton, maxnewton, linsolve) +end # ============================================================ # Drivers @@ -64,8 +78,8 @@ soft_symcover!(ϕ::PowerMean, a::AbstractVector, A::AbstractMatrix; kwargs...) = # The `:covariant` start makes a truncated iteration scale-covariant as well. function _soft_cover_powermean(ϕ::PowerMean, A::AbstractMatrix, fname::Symbol; kwargs...) if _abs_symmetric_exact(A) - a = _soft_symcover_powermean(ϕ, A, fname; kwargs...) - return a, copy(a) + a, stats = _soft_symcover_powermean(ϕ, A, fname; kwargs...) + return a, copy(a), stats end T = float(real(eltype(A))) a = similar(Array{T}, axes(A, 1)) @@ -77,14 +91,14 @@ end function _soft_cover_powermean!(::PowerMean{p}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix, fname::Symbol; - tol::Union{Real,Nothing}=nothing, kwargs...) where p + tol, maxiter, newton, maxnewton, linsolve) where p _prepare_soft_cover_start!(a, b, A, fname) T = _powermean_type(eltype(a), eltype(b), real(eltype(A))) rtol = _powermean_tol(T, tol) R, C = _log_support(_row_support(A, T)), _log_support(_col_support(A, T)) α = map(x -> x > 0 ? log(T(x)) : zero(T), a) β = map(x -> x > 0 ? log(T(x)) : zero(T), b) - stats = _powermean_solve!(α, β, R, C, T(p), rtol; kwargs...) + stats = _powermean_solve!(α, β, R, C, T(p), rtol; maxiter, newton, maxnewton, linsolve) _warn_powermean_unconverged(fname, stats, :sweeps) _check_representable(α, i -> !isempty(_slots(R, i)), β, j -> !isempty(_slots(C, j)), T, fname; gauge=true) for i in eachindex(a, α) @@ -93,7 +107,8 @@ function _soft_cover_powermean!(::PowerMean{p}, a::AbstractVector, b::AbstractVe for j in eachindex(b, β) b[j] = isempty(_slots(C, j)) ? zero(eltype(b)) : exp(β[j]) end - return _balance_cover!(a, b, A) + _balance_cover!(a, b, A) + return a, b, PowerMeanStats{T}(stats) end function _soft_symcover_powermean(ϕ::PowerMean, A::AbstractMatrix, fname::Symbol; kwargs...) @@ -106,20 +121,20 @@ function _soft_symcover_powermean(ϕ::PowerMean, A::AbstractMatrix, fname::Symbo end function _soft_symcover_powermean!(::PowerMean{p}, a::AbstractVector, A::AbstractMatrix, fname::Symbol; - tol::Union{Real,Nothing}=nothing, kwargs...) where p + tol, maxiter, newton, maxnewton, linsolve) where p _prepare_soft_symcover_start!(a, A, fname) T = _powermean_type(eltype(a), real(eltype(A))) rtol = _powermean_tol(T, tol) S = _log_support(_sym_support(A, T)) α = map(x -> x > 0 ? log(T(x)) : zero(T), a) - stats = _powermean_symsolve!(α, S, T(p), rtol; kwargs...) + stats = _powermean_symsolve!(α, S, T(p), rtol; maxiter, newton, maxnewton, linsolve) _warn_powermean_unconverged(fname, stats, :updates) _balance_bipartite_sym!(α, S) _check_representable(α, i -> !isempty(_slots(S, i)), nothing, nothing, T, fname; gauge=false) for i in eachindex(a, α) a[i] = isempty(_slots(S, i)) ? zero(eltype(a)) : exp(α[i]) end - return a + return a, PowerMeanStats{T}(stats) end # Whether `axes(A, 1) == axes(A, 2)` and `abs(A[i,j]) == abs(A[j,i])` for all `i, j`. diff --git a/src/stats.jl b/src/stats.jl index a86136a..d43242b 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -189,3 +189,56 @@ end LeastSquaresStats(nt::NamedTuple) = LeastSquaresStats(nt.linsolve, nt.lsqriters, nt.cgiters, nt.cholsolves, nt.precond, nt.fill_entries, nt.factor_flops) + +""" + PowerMeanStats{T} <: AbstractCoverStats + +Statistics of the `PowerMean` soft covers ([`soft_symcover`](@ref), +[`soft_cover`](@ref), and their mutating forms with a `PowerMean` penalty or no +penalty). + +The imbalance of a nonzero row (or column) is `|mean(r^p) - 1|`, the mean taken +over the ratios `r` of its nonzero entries. The solve stops when the largest +imbalance is at most `tol`. Power-mean updates run first; when they are slow, +Newton steps finish the solve. `converged`, `residual`, and `tolerance` return +the fields of the same names (`residual` returns `imbalance`); `iterations` +returns `nsweeps + nnewton`. + +# Fields + +- `converged::Bool`: `imbalance <= tol`. +- `tol::T`: the stopping tolerance. +- `imbalance::T`: the largest imbalance of the returned cover. +- `nsweeps::Int`: power-mean updates. For the asymmetric solver each counts an + update of all rows and then all columns (a sweep); for the symmetric solver, + one simultaneous update of all scales. +- `newton::Bool`: whether Newton steps ran. +- `nnewton::Int`: Newton steps. +- The remaining fields count the work of the Newton stage; their names and + meanings may change between releases: `npasses` (passes over the support), + `nfactor` (factorizations), and `linsolve` (the solver for the Newton + equations: `:dense`, `:cholesky`, or `:cg`, and `:none` when Newton did not + run). +""" +struct PowerMeanStats{T} <: AbstractCoverStats + converged::Bool + tol::T + imbalance::T + nsweeps::Int + newton::Bool + nnewton::Int + npasses::Int + nfactor::Int + linsolve::Symbol +end + +converged(s::PowerMeanStats) = s.converged +iterations(s::PowerMeanStats) = s.nsweeps + s.nnewton +residual(s::PowerMeanStats) = s.imbalance +tolerance(s::PowerMeanStats) = s.tol + +# Conversion from the statistics tuple of `_powermean_solve!` and +# `_powermean_symsolve!`; `T` is their working type. +PowerMeanStats{T}(nt::NamedTuple) where {T} = + PowerMeanStats{T}(nt.converged, nt.tol, nt.imbalance, nt.nsweeps, nt.newton, nt.nnewton, + nt.npasses, nt.nfactor, nt.linsolve) diff --git a/test/element_types.jl b/test/element_types.jl index 6928a2f..cf77277 100644 --- a/test/element_types.jl +++ b/test/element_types.jl @@ -4,7 +4,7 @@ @testset "Float32 flows through the family" begin A = Float32[4 1.5; 1.5 1] - for a in (symcover(A), soft_symcover(A), soft_symcover(AbsLog{1}(), A), + for a in (symcover(A), first(soft_symcover(A)), soft_symcover(AbsLog{1}(), A), soft_symcover(AbsLinear{1}(), A), first(symcover_min(AbsLog{2}(), A)), first(soft_symcover(AbsLog{2}(), A))) @test eltype(a) === Float32 @@ -77,7 +77,7 @@ @testset "BigFloat flows through the family" begin A = BigFloat[4 1.5; 1.5 1] - for a in (symcover(A), soft_symcover(A), first(symcover_min(AbsLog{2}(), A)), + for a in (symcover(A), first(soft_symcover(A)), first(symcover_min(AbsLog{2}(), A)), first(soft_symcover(AbsLog{2}(), A))) @test eltype(a) === BigFloat @test all(isfinite, a) diff --git a/test/invariants.jl b/test/invariants.jl index ae99263..5ac3a10 100644 --- a/test/invariants.jl +++ b/test/invariants.jl @@ -16,8 +16,8 @@ const SYM_NOTIONS = ( (name = "soft_symcover(AbsLog{1})", f = A -> soft_symcover(AbsLog{1}(), A), hard = false, rtol = 1e-5), (name = "soft_symcover(AbsLinear{2})", f = A -> soft_symcover(AbsLinear{2}(), A), hard = false, rtol = 1e-5), (name = "soft_symcover(AbsLinear{1})", f = A -> soft_symcover(AbsLinear{1}(), A), hard = false, rtol = 1e-5), - (name = "soft_symcover(PowerMean{2})", f = A -> soft_symcover(PowerMean{2}(), A), hard = false, rtol = 1e-12), - (name = "soft_symcover(PowerMean{1})", f = A -> soft_symcover(PowerMean{1}(), A), hard = false, rtol = 1e-12), + (name = "soft_symcover(PowerMean{2})", f = A -> first(soft_symcover(PowerMean{2}(), A)), hard = false, rtol = 1e-12), + (name = "soft_symcover(PowerMean{1})", f = A -> first(soft_symcover(PowerMean{1}(), A)), hard = false, rtol = 1e-12), ) const GEN_NOTIONS = ( @@ -30,8 +30,8 @@ const GEN_NOTIONS = ( (name = "soft_cover(AbsLog{1})", f = A -> soft_cover(AbsLog{1}(), A), hard = false, rtol = 1e-5), (name = "soft_cover(AbsLinear{2})", f = A -> soft_cover(AbsLinear{2}(), A), hard = false, rtol = 1e-5), (name = "soft_cover(AbsLinear{1})", f = A -> soft_cover(AbsLinear{1}(), A), hard = false, rtol = 1e-5), - (name = "soft_cover(PowerMean{2})", f = A -> soft_cover(PowerMean{2}(), A), hard = false, rtol = 1e-12), - (name = "soft_cover(PowerMean{3})", f = A -> soft_cover(PowerMean{3}(), A), hard = false, rtol = 1e-12), + (name = "soft_cover(PowerMean{2})", f = A -> scales(soft_cover(PowerMean{2}(), A)), hard = false, rtol = 1e-12), + (name = "soft_cover(PowerMean{3})", f = A -> scales(soft_cover(PowerMean{3}(), A)), hard = false, rtol = 1e-12), ) @testset "cross-notion invariants" begin diff --git a/test/powermean.jl b/test/powermean.jl index 3a043f9..d4cf58f 100644 --- a/test/powermean.jl +++ b/test/powermean.jl @@ -34,7 +34,7 @@ end R11 = (2 / (1 + κ^(p/2)))^(1/p) R12 = (2 * κ^(p/2) / (1 + κ^(p/2)))^(1/p) E = (4 / p) * log(cosh((p / 4) * log(κ))) - a = soft_symcover(PowerMean{p}(), A) + a, _ = soft_symcover(PowerMean{p}(), A) R = abs.(A) ./ (a .* a') @test R[1, 1] ≈ R11 rtol=1e-13 @test R[2, 2] ≈ R11 rtol=1e-13 @@ -91,7 +91,7 @@ end end S = pm_testmatrix(rng, 20, 20; σ=2.0, z=0.0) S = S + S' - s = @test_logs soft_symcover(PowerMean{3}(), S) + s, _ = @test_logs soft_symcover(PowerMean{3}(), S) @test pm_imbalance(3, s, s, S) < 1e-11 end @@ -102,11 +102,11 @@ end a0 = initialize_symcover(A; strategy=:geomean, feasible=:none) E = Float64[] for k in 0:12 - ak = @test_logs (:warn, r"soft_symcover!: the power-mean iteration ended after") match_mode=:any soft_symcover!(PowerMean{2}(), copy(a0), A; maxiter=k, newton=false) + ak, _ = @test_logs (:warn, r"soft_symcover!: the power-mean iteration ended after") match_mode=:any soft_symcover!(PowerMean{2}(), copy(a0), A; maxiter=k, newton=false) push!(E, cover_objective(PowerMean{2}(), ak, A)) end @test all(<(0), diff(E)) - @test cover_objective(PowerMean{2}(), soft_symcover(A), A) < E[end] + @test cover_objective(PowerMean{2}(), first(soft_symcover(A)), A) < E[end] end @testset "non-convergence warns" begin @@ -138,8 +138,8 @@ end S = pm_testmatrix(rng, 25, 25; σ=2.0, z=0.3) S = S + S' for p in (1, 2, 3) - a = soft_symcover(PowerMean{p}(), S) - @test soft_symcover(PowerMean{p}(), d1 .* S .* d1') ≈ d1 .* a rtol=1e-11 + a, _ = soft_symcover(PowerMean{p}(), S) + @test first(soft_symcover(PowerMean{p}(), d1 .* S .* d1')) ≈ d1 .* a rtol=1e-11 end end @@ -183,7 +183,7 @@ end @test pm_imbalance(2, a, b, Ssp) < 1e-11 # The symmetric iteration handles diagonally dominant banded support directly. T = SymTridiagonal(fill(10.0, 100), ones(99)) - a = @test_logs soft_symcover(T) + a, _ = @test_logs soft_symcover(T) @test pm_imbalance(2, a, a, T) < 1e-11 end @@ -195,20 +195,20 @@ end @test isbalanced(a, b, A) # Loopless bipartite symmetric support: products fixed, sides balanced. S = [0 3.0 0; 3.0 0 5.0; 0 5.0 0] - s = soft_symcover(S) + s, _ = soft_symcover(S) @test s[1] * s[2] ≈ 3 && s[2] * s[3] ≈ 5 @test 2 * log(s[2]) ≈ log(s[1]) + log(s[3]) - @test soft_symcover!(PowerMean{2}(), [7.0, 0.01, 2.0], S) ≈ s - @test soft_symcover([0 1; 1 0]) ≈ [1.0, 1.0] + @test first(soft_symcover!(PowerMean{2}(), [7.0, 0.01, 2.0], S)) ≈ s + @test first(soft_symcover([0 1; 1 0])) ≈ [1.0, 1.0] end @testset "element types" begin A = [4.0 1.5 0.5; 1.5 1.0 2.0; 0.5 2.0 3.0] - a64 = soft_symcover(A) - a32 = soft_symcover(Float32.(A)) + a64, _ = soft_symcover(A) + a32, _ = soft_symcover(Float32.(A)) @test a32 isa Vector{Float32} @test a32 ≈ a64 rtol=2eps(Float32) - abig = soft_symcover(big.(A)) + abig, _ = soft_symcover(big.(A)) @test abig isa Vector{BigFloat} Rbig = abs.(big.(A)) ./ (abig .* abig') @test maximum(abs, sum(Rbig .^ 2; dims=2) ./ 3 .- 1) < 1e-70 @@ -225,17 +225,17 @@ end @test axes(ao, 1) == axes(Ao, 1) && axes(bo, 1) == axes(Ao, 2) @test collect(ao) == a && collect(bo) == b Av = view(hcat(A, ones(6)), :, 1:4) - @test soft_cover(Av) == (a, b) + @test scales(soft_cover(Av)) == (a, b) ao2, bo2 = soft_cover!(PowerMean{2}(), OffsetArray(ones(6), -2:3), OffsetArray(ones(4), 5:8), Ao) @test collect(ao2) .* collect(bo2)' ≈ a .* b' rtol=1e-12 S = A[1:4, :] + A[1:4, :]' - s = soft_symcover(S) + s, _ = soft_symcover(S) So = OffsetArray(S, 0:3, 0:3) - so = soft_symcover(So) + so, _ = soft_symcover(So) @test axes(so, 1) == axes(So, 1) @test collect(so) == s - @test soft_symcover(view(S, 1:4, 1:4)) == s + @test first(soft_symcover(view(S, 1:4, 1:4))) == s @test_throws DimensionMismatch soft_symcover!(PowerMean{2}(), ones(3), S) @test_throws DimensionMismatch soft_cover!(PowerMean{2}(), ones(5), ones(4), A) @test_throws "soft_symcover requires a square matrix" soft_symcover(A) @@ -264,7 +264,7 @@ end @test nslow > 0 for S in [hard_sym; sp_sym] n = size(S, 1) - a = @test_logs soft_symcover(S) + a, _ = @test_logs soft_symcover(S) @test pm_imbalance(2, a, a, S) < 1e-11 # The alternating solver with Newton, on symmetric input. x, y = @test_logs soft_cover!(PowerMean{2}(), ones(n), ones(n), S) @@ -278,8 +278,8 @@ end # With `tol=0` the solve continues to roundoff, where some steps change the # objective by less than it can resolve and are judged by the imbalance. for S in hard_sym - a = @test_logs (:warn, r"tolerance 0\.0") soft_symcover(S; tol=0.0) - @test a ≈ soft_symcover(S) rtol=1e-11 + a, _ = @test_logs (:warn, r"tolerance 0\.0") soft_symcover(S; tol=0.0) + @test a ≈ first(soft_symcover(S)) rtol=1e-11 end for A in [hard_gen; sp_gen] a, b = @test_logs soft_cover(A) @@ -302,14 +302,14 @@ end x, y = @test_logs soft_cover(A; maxiter=2, linsolve) @test x .* y' ≈ a .* b' rtol=1e-11 @test isbalanced(x, y, A) - @test soft_symcover(S; maxiter=2, linsolve) ≈ soft_symcover(S; newton=false) rtol=1e-11 + @test first(soft_symcover(S; maxiter=2, linsolve)) ≈ first(soft_symcover(S; newton=false)) rtol=1e-11 end Asp = sparse(A) @test soft_cover(Asp; maxiter=2)[1] .* soft_cover(Asp; maxiter=2)[2]' ≈ a .* b' rtol=1e-11 # Bipartite symmetric support: Newton fixes the gauge and the balance holds. Sb = [0 3.0 0; 3.0 0 5.0; 0 5.0 0] - s = @test_logs soft_symcover(Sb; maxiter=0) - @test s ≈ soft_symcover(Sb; newton=false) rtol=1e-12 + s, _ = @test_logs soft_symcover(Sb; maxiter=0) + @test s ≈ first(soft_symcover(Sb; newton=false)) rtol=1e-12 # Generic indexing and element types through the Newton stage. Ao = OffsetArray(A, -3:26, 4:23) ao, bo = soft_cover(Ao; maxiter=2) @@ -317,15 +317,15 @@ end @test collect(ao) .* collect(bo)' ≈ a .* b' rtol=1e-11 @test soft_cover(view(A, :, :); maxiter=2)[1] ≈ soft_cover(A; maxiter=2)[1] rtol=1e-12 So = OffsetArray(S, 0:19, 0:19) - so = soft_symcover(So; maxiter=2) + so, _ = soft_symcover(So; maxiter=2) @test axes(so, 1) == 0:19 - @test collect(so) ≈ soft_symcover(S) rtol=1e-11 + @test collect(so) ≈ first(soft_symcover(S)) rtol=1e-11 x32, y32 = soft_cover(Float32.(A); maxiter=2) @test x32 isa Vector{Float32} && y32 isa Vector{Float32} @test x32 .* y32' ≈ a .* b' rtol=10eps(Float32) - sbig = soft_symcover(big.(S); maxiter=2) + sbig, _ = soft_symcover(big.(S); maxiter=2) @test sbig isa Vector{BigFloat} - @test sbig ≈ soft_symcover(S) rtol=1e-11 + @test sbig ≈ first(soft_symcover(S)) rtol=1e-11 # Solver statistics expose the combined cost. R = MatrixCovers._log_support(MatrixCovers._row_support(sp_gen[1], Float64)) @@ -345,10 +345,10 @@ end S0 = S - Diagonal(diag(S)) # zero diagonal for M in (S, S0, sparse(S), sparse(S0), Symmetric(sparse(triu(S))), -S) a, b = soft_cover(M) - @test a == b == soft_symcover(M) + @test a == b == first(soft_symcover(M)) @test isbalanced(a, b, M) # The alternating iteration reaches the same minimizer. - x, y = soft_cover!(PowerMean{2}(), ones(20), ones(20), M) + x, y, _ = soft_cover!(PowerMean{2}(), ones(20), ones(20), M) supp = findall(!iszero, Matrix(M)) @test (a .* b')[supp] ≈ (x .* y')[supp] rtol=1e-11 end @@ -373,7 +373,7 @@ end So = OffsetArray(S, -1:18, -1:18) ao, bo = soft_cover(So) @test axes(ao, 1) == axes(bo, 1) == -1:18 - @test collect(ao) == soft_symcover(S) + @test collect(ao) == first(soft_symcover(S)) a32, b32 = soft_cover(Float32.(S)) @test a32 isa Vector{Float32} && b32 isa Vector{Float32} && a32 == b32 # Rectangular and offset-mismatched axes are never symmetric. @@ -394,11 +394,48 @@ end @testset "default dispatch" begin A = [1.0 2.0 3.0; 6.0 5.0 4.0] S = [4.0 1.0; 1.0 2.0] - @test soft_cover(A) == soft_cover(PowerMean{2}(), A) - @test soft_symcover(S) == soft_symcover(PowerMean{2}(), S) + @test scales(soft_cover(A)) == scales(soft_cover(PowerMean{2}(), A)) + @test first(soft_symcover(S)) == first(soft_symcover(PowerMean{2}(), S)) a, b = soft_cover(A) x, y = soft_cover!(ones(2), ones(3), A) @test x .* y' ≈ a .* b' rtol=1e-12 - @test soft_symcover!(ones(2), S) ≈ soft_symcover(S) rtol=1e-12 + @test first(soft_symcover!(ones(2), S)) ≈ first(soft_symcover(S)) rtol=1e-12 + end + + @testset "statistics" begin + A = [1.0 2.0 3.0; 6.0 5.0 4.0] + S = [4.0 1.0; 1.0 2.0] + a, b, st = soft_cover(A) + @test st isa MatrixCovers.PowerMeanStats{Float64} + @test MatrixCovers.converged(st) + @test MatrixCovers.residual(st) <= MatrixCovers.tolerance(st) == 4096 * eps() + @test MatrixCovers.residual(st) ≈ pm_imbalance(2, a, b, A) atol=1e-14 + @test MatrixCovers.iterations(st) == st.nsweeps + st.nnewton > 0 + a, st = soft_symcover(Float32.(S); tol=1e-6) + @test st isa MatrixCovers.PowerMeanStats{Float64} + @test MatrixCovers.converged(st) && MatrixCovers.tolerance(st) == 1e-6 + # The Newton stage is reported. + rng = StableRNG(3) + W = exp.(3 .* randn(rng, 30, 30)) + _, _, st = soft_cover(W; maxiter=3) + @test MatrixCovers.converged(st) && st.newton && st.nsweeps == 3 + @test st.nnewton > 0 && st.npasses > 0 && st.linsolve in (:dense, :cholesky, :cg) + @test MatrixCovers.iterations(st) == 3 + st.nnewton + # A truncated solve reports failure as well as warning. + _, _, st = @test_logs (:warn, r"^soft_cover: the power-mean iteration ended") soft_cover(W; maxiter=3, newton=false) + @test !MatrixCovers.converged(st) && !st.newton && st.linsolve === :none + @test MatrixCovers.residual(st) > MatrixCovers.tolerance(st) + @test MatrixCovers.iterations(st) == 3 + _, st = @test_logs (:warn, r"^soft_symcover!: the power-mean iteration ended") soft_symcover!(ones(30), W + W'; maxiter=2, newton=false) + @test !MatrixCovers.converged(st) && MatrixCovers.iterations(st) == 2 + # Symmetric input to `soft_cover` reports the symmetric solve. + _, _, st = soft_cover(S) + @test MatrixCovers.converged(st) + @test (@inferred soft_cover(A)) isa Tuple{Vector{Float64},Vector{Float64},MatrixCovers.PowerMeanStats{Float64}} + @test (@inferred soft_symcover(PowerMean{3}(), S))[2] isa MatrixCovers.PowerMeanStats{Float64} + @inferred soft_cover(PowerMean{1}(), sparse(A)) + @inferred soft_cover!(PowerMean{2}(), ones(2), Float32[1, 1, 1], A) + @inferred soft_symcover!(ones(2), S) + @inferred soft_symcover(big.(S)) end end diff --git a/test/soft_covers.jl b/test/soft_covers.jl index 790d4fc..a1d187f 100644 --- a/test/soft_covers.jl +++ b/test/soft_covers.jl @@ -76,7 +76,7 @@ end # Default dispatch uses PowerMean{2} A = [2.0 1.0; 1.0 3.0] - @test soft_symcover(A) == soft_symcover(PowerMean{2}(), A) + @test first(soft_symcover(A)) == first(soft_symcover(PowerMean{2}(), A)) end @testset "soft_cover" begin @@ -96,7 +96,7 @@ end # Default dispatch uses PowerMean{2} A = [1.0 2.0 3.0; 6.0 5.0 4.0] - @test soft_cover(A) == soft_cover(PowerMean{2}(), A) + @test scales(soft_cover(A)) == scales(soft_cover(PowerMean{2}(), A)) # The returned objective never exceeds that of the geometric-mean start. for A in ([1.0 2.0 3.0; 6.0 5.0 4.0], @@ -491,10 +491,10 @@ end @testset "no-ϕ form defaults to PowerMean{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!(PowerMean{2}(), a2, Asym) + @test first(soft_symcover!(a1, Asym)) == first(soft_symcover!(PowerMean{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!(PowerMean{2}(), b2, c2, Agen) + @test scales(soft_cover!(b1, c1, Agen)) == scales(soft_cover!(PowerMean{2}(), b2, c2, Agen)) end # Convex refiners ignore the basin; nonconvex refiners do not. @@ -604,7 +604,7 @@ MatrixCovers.require_abs_symmetric(::HookOnlyMatrix, fname) = nothing end end - @test soft_symcover(PowerMean{2}(), M) ≈ soft_symcover(PowerMean{2}(), dense) rtol=1e-12 + @test first(soft_symcover(PowerMean{2}(), M)) ≈ first(soft_symcover(PowerMean{2}(), dense)) rtol=1e-12 for ϕ in (AbsLinear{1}(), AbsLinear{2}(), AbsLog{1}(), AbsLog{2}()) @test symcover(ϕ, M) ≈ symcover(ϕ, dense) rtol=1e-10 @test scales(soft_symcover(ϕ, M)) ≈ scales(soft_symcover(ϕ, dense)) rtol=1e-6 diff --git a/test/unitful.jl b/test/unitful.jl index 3209cba..3fef701 100644 --- a/test/unitful.jl +++ b/test/unitful.jl @@ -239,11 +239,11 @@ @test unit.(a) == unit.(b) == UA a = initialize_symcover(A) - @test soft_symcover!(a, A) === a + @test first(soft_symcover!(a, A)) === a @test unit.(a) == UA a, b = initialize_cover(A) - @test soft_cover!(a, b, A) == (a, b) + @test scales(soft_cover!(a, b, A)) == (a, b) @test unit.(a) == unit.(b) == UA a = initialize_symcover(A) @@ -255,11 +255,11 @@ @test unit.(a) == unit.(b) == UA a = initialize_symcover(A) - @test soft_symcover!(PowerMean{1}(), a, A) === a + @test first(soft_symcover!(PowerMean{1}(), a, A)) === a @test unit.(a) == UA a, b = initialize_cover(A) - @test soft_cover!(PowerMean{1}(), a, b, A) == (a, b) + @test scales(soft_cover!(PowerMean{1}(), a, b, A)) == (a, b) @test unit.(a) == unit.(b) == UA # A start in a dimensionally equivalent spelling is converted, not rejected. From 3469c7f11ad2d414075a12616f16317e228bc1bb Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 9 Oct 2026 05:34:16 -0500 Subject: [PATCH 4/7] Return solver statistics from JuMP-backed covers The covers computed through JuMP now return a statistics object as their last output, like the native solvers: `symcover_min`, `cover_min`, `soft_symcover`, `soft_cover`, and their mutating forms with `AbsLog{1}` (HiGHS) or `AbsLinear` (Ipopt). `ExternalSolverStats` records one entry per optimization stage: the solver, the termination status (as a `Symbol`), the objective value of the model, the iteration count the solver reports, and its run time. The `AbsLog{1}` covers run two stages (the LP, then the QP over its optimal face); Ipopt runs one. Failures still throw `SolverFailure`, so `converged` is always `true`; `iterations`, `residual`, and `tolerance` are not implemented. The `AbsLinear` multistarts return `MultistartStats`, which holds the label, statistics, and objective of each start whose refinement succeeded, the index of the selected one, and the labels of the random starts that failed. Its accessors delegate to the selected start. The test-only `labels`/`objs` keywords of the soft multistart are removed, since the statistics carry the same information. Assisted-by: Claude Opus 5.5 --- ext/MatrixCoversIpoptExt.jl | 46 +++++--- ext/MatrixCoversJuMPExt.jl | 68 +++++++---- src/MatrixCovers.jl | 2 +- src/minimal_covers.jl | 22 ++-- src/soft_covers.jl | 70 ++++++----- src/stats.jl | 99 ++++++++++++++++ test/element_types.jl | 6 +- test/extensions.jl | 223 +++++++++++++++++++++++++----------- test/invariants.jl | 26 ++--- test/iscover.jl | 2 +- test/soft_covers.jl | 121 +++++++++---------- test/storage_types.jl | 8 +- test/unitful.jl | 10 +- 13 files changed, 470 insertions(+), 233 deletions(-) diff --git a/ext/MatrixCoversIpoptExt.jl b/ext/MatrixCoversIpoptExt.jl index 15f5f5a..2cb0ee8 100644 --- a/ext/MatrixCoversIpoptExt.jl +++ b/ext/MatrixCoversIpoptExt.jl @@ -4,7 +4,7 @@ using JuMP: JuMP, @variable, @objective, @constraint using Ipopt: Ipopt using MatrixCovers using MatrixCovers: AbsLinear -using MatrixCovers: _edge_list, _sym_edge_list, _degrees, SolverFailure +using MatrixCovers: _edge_list, _sym_edge_list, _degrees, SolverFailure, ExternalSolverStats # Models use 1-based positions and scatter results back to `A`'s axes. Support is # gathered as an O(nnz) edge list; unsupported scales are zero. @@ -14,6 +14,14 @@ using MatrixCovers: _edge_list, _sym_edge_list, _degrees, SolverFailure check_solved(model, fname) = MatrixCovers.check_solved(JuMP.termination_status(model), "Ipopt", fname) +# Check that the solve of `model` succeeded and return its statistics. +function solved_stats(model, fname) + check_solved(model, fname) + return ExternalSolverStats(:Ipopt, [Symbol(JuMP.termination_status(model))], + [JuMP.objective_value(model)], [JuMP.barrier_iterations(model)], + [JuMP.solve_time(model)]) +end + # Suppress both solver output and Ipopt's startup banner. # # `unscaled=true` is for the AbsLinear{2} models. Their (1-r)^2 objective is nearly flat far @@ -102,13 +110,13 @@ function MatrixCovers.symcover_min!(::AbsLinear{2}, a::AbstractVector, A) @constraint(model, α[ti[k]] + α[tj[k]] >= tlog[k]) end JuMP.optimize!(model) - check_solved(model, "symcover_min!") + stats = solved_stats(model, "symcover_min!") lα = JuMP.value.(α) check_not_plateau("symcover_min!", fi, fj, flog, lα, lα, pr, pr) for (i, k) in pairs(pr) a[k] = supported[i] ? exp(lα[i]) : zero(T) end - return a + return a, stats end function MatrixCovers.symcover_min!(::AbsLinear{1}, a::AbstractVector, A) @@ -135,13 +143,13 @@ function MatrixCovers.symcover_min!(::AbsLinear{1}, a::AbstractVector, A) @constraint(model, α[ti[k]] + α[tj[k]] >= tlog[k]) end JuMP.optimize!(model) - check_solved(model, "symcover_min!") + stats = solved_stats(model, "symcover_min!") lα = JuMP.value.(α) check_not_plateau("symcover_min!", fi, fj, flog, lα, lα, pr, pr) for (i, k) in pairs(pr) a[k] = supported[i] ? exp(lα[i]) : zero(T) end - return a + return a, stats end # ============================================================ @@ -172,7 +180,7 @@ function MatrixCovers.cover_min!(::AbsLinear{2}, a::AbstractVector, b::AbstractV end _pin_gauge!(model, α, β, A, nza, nzb) JuMP.optimize!(model) - check_solved(model, "cover_min!") + stats = solved_stats(model, "cover_min!") lα, lβ = JuMP.value.(α), JuMP.value.(β) check_not_plateau("cover_min!", ei, ej, elog, lα, lβ, pr, pc) for (i, k) in pairs(pr) @@ -182,7 +190,8 @@ function MatrixCovers.cover_min!(::AbsLinear{2}, a::AbstractVector, b::AbstractV b[k] = nzb[j] > 0 ? exp(lβ[j]) : zero(T) end MatrixCovers._balance_cover!(a, b, A) - return MatrixCovers.inflate_feasible!(a, b, A) + MatrixCovers.inflate_feasible!(a, b, A) + return a, b, stats end function MatrixCovers.cover_min!(::AbsLinear{1}, a::AbstractVector, b::AbstractVector, A) @@ -211,7 +220,7 @@ function MatrixCovers.cover_min!(::AbsLinear{1}, a::AbstractVector, b::AbstractV end _pin_gauge!(model, α, β, A, nza, nzb) JuMP.optimize!(model) - check_solved(model, "cover_min!") + stats = solved_stats(model, "cover_min!") lα, lβ = JuMP.value.(α), JuMP.value.(β) check_not_plateau("cover_min!", ei, ej, elog, lα, lβ, pr, pc) for (i, k) in pairs(pr) @@ -221,7 +230,8 @@ function MatrixCovers.cover_min!(::AbsLinear{1}, a::AbstractVector, b::AbstractV b[k] = nzb[j] > 0 ? exp(lβ[j]) : zero(T) end MatrixCovers._balance_cover!(a, b, A) - return MatrixCovers.inflate_feasible!(a, b, A) + MatrixCovers.inflate_feasible!(a, b, A) + return a, b, stats end # ============================================================ @@ -246,13 +256,13 @@ function MatrixCovers.soft_symcover!(::AbsLinear{2}, a::AbstractVector, A) @objective(model, Min, sum(tw[k] * (1 - exp(tlog[k] - α[ti[k]] - α[tj[k]]))^2 for k in eachindex(ti)) + n_zeros) JuMP.optimize!(model) - check_solved(model, "soft_symcover!") + stats = solved_stats(model, "soft_symcover!") lα = JuMP.value.(α) check_not_plateau("soft_symcover!", fi, fj, flog, lα, lα, pr, pr) for (i, k) in pairs(pr) a[k] = supported[i] ? exp(lα[i]) : zero(T) end - return a + return a, stats end function MatrixCovers.soft_symcover!(::AbsLinear{1}, a::AbstractVector, A) @@ -276,13 +286,13 @@ function MatrixCovers.soft_symcover!(::AbsLinear{1}, a::AbstractVector, A) end @objective(model, Min, sum(tw[k] * t[k] for k in eachindex(ti)) + n_zeros) JuMP.optimize!(model) - check_solved(model, "soft_symcover!") + stats = solved_stats(model, "soft_symcover!") lα = JuMP.value.(α) check_not_plateau("soft_symcover!", fi, fj, flog, lα, lα, pr, pr) for (i, k) in pairs(pr) a[k] = supported[i] ? exp(lα[i]) : zero(T) end - return a + return a, stats end # ============================================================ @@ -311,7 +321,7 @@ function MatrixCovers.soft_cover!(::AbsLinear{2}, a::AbstractVector, b::Abstract sum((1 - exp(elog[e] - α[ei[e]] - β[ej[e]]))^2 for e in eachindex(ei)) + n_zeros) _pin_gauge!(model, α, β, A, nza, nzb) JuMP.optimize!(model) - check_solved(model, "soft_cover!") + stats = solved_stats(model, "soft_cover!") lα, lβ = JuMP.value.(α), JuMP.value.(β) check_not_plateau("soft_cover!", ei, ej, elog, lα, lβ, pr, pc) for (i, k) in pairs(pr) @@ -320,7 +330,8 @@ function MatrixCovers.soft_cover!(::AbsLinear{2}, a::AbstractVector, b::Abstract for (j, k) in pairs(pc) b[k] = nzb[j] > 0 ? exp(lβ[j]) : zero(T) end - return MatrixCovers._balance_cover!(a, b, A) + MatrixCovers._balance_cover!(a, b, A) + return a, b, stats end function MatrixCovers.soft_cover!(::AbsLinear{1}, a::AbstractVector, b::AbstractVector, A) @@ -347,7 +358,7 @@ function MatrixCovers.soft_cover!(::AbsLinear{1}, a::AbstractVector, b::Abstract @objective(model, Min, sum(t) + n_zeros) _pin_gauge!(model, α, β, A, nza, nzb) JuMP.optimize!(model) - check_solved(model, "soft_cover!") + stats = solved_stats(model, "soft_cover!") lα, lβ = JuMP.value.(α), JuMP.value.(β) check_not_plateau("soft_cover!", ei, ej, elog, lα, lβ, pr, pc) for (i, k) in pairs(pr) @@ -356,7 +367,8 @@ function MatrixCovers.soft_cover!(::AbsLinear{1}, a::AbstractVector, b::Abstract for (j, k) in pairs(pc) b[k] = nzb[j] > 0 ? exp(lβ[j]) : zero(T) end - return MatrixCovers._balance_cover!(a, b, A) + MatrixCovers._balance_cover!(a, b, A) + return a, b, stats end end # module MatrixCoversIpoptExt diff --git a/ext/MatrixCoversJuMPExt.jl b/ext/MatrixCoversJuMPExt.jl index 1289ee4..66135c9 100644 --- a/ext/MatrixCoversJuMPExt.jl +++ b/ext/MatrixCoversJuMPExt.jl @@ -3,13 +3,40 @@ module MatrixCoversJuMPExt using JuMP: JuMP, @variable, @objective, @constraint using HiGHS: HiGHS using MatrixCovers -using MatrixCovers: AbsLog +using MatrixCovers: AbsLog, ExternalSolverStats using MatrixCovers: _edge_list, _sym_edge_list, _degrees using LinearAlgebra: dot check_solved(model, fname) = MatrixCovers.check_solved(JuMP.termination_status(model), "HiGHS", fname) +_highs_stats() = ExternalSolverStats(:HiGHS, Symbol[], Float64[], Int[], Float64[]) + +# HiGHS reports its simplex, interior-point (with crossover), and QP iterations +# separately, and MOI exposes only the first two. +function _highs_iterations(model) + opt = JuMP.unsafe_backend(model) + n = 0 + for key in ("simplex_iteration_count", "ipm_iteration_count", "crossover_iteration_count", + "qp_iteration_count") + r = Ref{Cint}(0) + HiGHS.Highs_getIntInfoValue(opt, key, r) == HiGHS.kHighsStatusOk || + error("HiGHS did not report $key") + n += r[] + end + return n +end + +# Check that the last solve of `model` succeeded and append it to `stats` as a stage. +function record_solved!(stats::ExternalSolverStats, model, fname) + check_solved(model, fname) + push!(stats.status, Symbol(JuMP.termination_status(model))) + push!(stats.objective, JuMP.objective_value(model)) + push!(stats.niters, _highs_iterations(model)) + push!(stats.solvetime, JuMP.solve_time(model)) + return stats +end + # Models use 1-based positions and scatter results back to `A`'s axes. Support is # gathered as an O(nnz) edge list; unsupported scales are zero. @@ -43,23 +70,23 @@ MatrixCovers.symcover_min(::AbsLog{1}, A) = _symcover_min_abslog1(A, nothing) function MatrixCovers.symcover_min!(::AbsLog{1}, a::AbstractVector, A) MatrixCovers._prepare_symcover_start!(a, A) - a .= _symcover_min_abslog1(A, a) - return a + anew, stats = _symcover_min_abslog1(A, a) + a .= anew + return a, stats end # Slack for re-evaluating the `AbsLog{1}` optimum during tie-breaking. const LEX_L1_SLACK = 1e-9 # Break `AbsLog{1}` ties with `AbsLog{2}` using full-grid support weights. -function _minimize_l2_over_l1_face!(model, lin, residuals, fname) - isempty(residuals) && return nothing +function _minimize_l2_over_l1_face!(stats, model, lin, residuals, fname) + isempty(residuals) && return stats linopt = JuMP.value(lin) l1 = sum(JuMP.value, residuals) # the AbsLog{1} objective attained @constraint(model, lin <= linopt + LEX_L1_SLACK * max(one(l1), l1)) @objective(model, Min, sum(r^2 for r in residuals)) JuMP.optimize!(model) - check_solved(model, fname) - return nothing + return record_solved!(stats, model, fname) end # The soft `AbsLog{1}` objective is convex (an L1 fit in the log scales), so it @@ -68,18 +95,19 @@ MatrixCovers.soft_symcover(::AbsLog{1}, A::AbstractMatrix) = _symcover_min_abslo function MatrixCovers.soft_symcover!(::AbsLog{1}, a::AbstractVector, A::AbstractMatrix) MatrixCovers._prepare_soft_symcover_start!(a, A) - a .= _symcover_min_abslog1(A, a; soft=true) - return a + anew, stats = _symcover_min_abslog1(A, a; soft=true) + a .= anew + return a, stats end MatrixCovers.soft_cover(::AbsLog{1}, A::AbstractMatrix) = _cover_min_abslog1(A, nothing; soft=true) function MatrixCovers.soft_cover!(::AbsLog{1}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix) MatrixCovers._prepare_soft_cover_start!(a, b, A) - anew, bnew = _cover_min_abslog1(A, (a, b); soft=true) + anew, bnew, stats = _cover_min_abslog1(A, (a, b); soft=true) a .= anew b .= bnew - return a, b + return a, b, stats end # Symmetric `AbsLog{1}` LP. A second stage selects the canonical point on the @@ -124,9 +152,9 @@ function _symcover_min_abslog1(A, start; soft::Bool=false) end @objective(model, Min, lin) JuMP.optimize!(model) - check_solved(model, fname) + stats = record_solved!(_highs_stats(), model, fname) residuals = [α[ei[e]] + α[ej[e]] - elog[e] for e in eachindex(ei)] - _minimize_l2_over_l1_face!(model, lin, residuals, fname) + _minimize_l2_over_l1_face!(stats, model, lin, residuals, fname) αv = [JuMP.value(α[i]) for i in 1:n] if soft # Without the coverage constraints, bipartite components keep their gauge. @@ -141,7 +169,7 @@ function _symcover_min_abslog1(A, start; soft::Bool=false) for (i, k) in pairs(pr) a[k] = supported[i] ? exp(αv[i]) : zero(T) end - return a + return a, stats end function MatrixCovers.cover_min_jump(::AbsLog{2}, A) @@ -221,10 +249,10 @@ MatrixCovers.cover_min(::AbsLog{1}, A) = _cover_min_abslog1(A, nothing) function MatrixCovers.cover_min!(::AbsLog{1}, a::AbstractVector, b::AbstractVector, A) MatrixCovers._prepare_cover_start!(a, b, A) - anew, bnew = _cover_min_abslog1(A, (a, b)) + anew, bnew, stats = _cover_min_abslog1(A, (a, b)) a .= anew b .= bnew - return a, b + return a, b, stats end # Asymmetric `AbsLog{1}` LP. Balance globally in the model and per component @@ -273,9 +301,9 @@ function _cover_min_abslog1(A, start; soft::Bool=false) # Pin the global row/column gauge; post-processing handles components. @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n)) JuMP.optimize!(model) - check_solved(model, fname) + stats = record_solved!(_highs_stats(), model, fname) residuals = [α[ei[e]] + β[ej[e]] - elog[e] for e in eachindex(ei)] - _minimize_l2_over_l1_face!(model, lin, residuals, fname) + _minimize_l2_over_l1_face!(stats, model, lin, residuals, fname) a = similar(Array{T}, axr) b = similar(Array{T}, axc) for (i, k) in pairs(pr) @@ -285,8 +313,8 @@ function _cover_min_abslog1(A, start; soft::Bool=false) b[k] = nzb[j] > 0 ? exp(JuMP.value(β[j])) : zero(T) end MatrixCovers._balance_cover!(a, b, A) - soft && return a, b - return MatrixCovers.inflate_feasible!(a, b, A) + soft || MatrixCovers.inflate_feasible!(a, b, A) + return a, b, stats end diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index 946e073..ca96764 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -18,7 +18,7 @@ export cover_transversal, cover_transversal! # `public` is parsed as a keyword only from Julia 1.11; this package supports 1.10. @static if VERSION >= v"1.11" - eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym, scalar_type, SupportComponents, support_components, ncomponents, rowcomponent, colcomponent, SolverFailure, AbstractCoverStats, AugmentedLagrangianStats, ActiveSetStats, LeastSquaresStats, PowerMeanStats, converged, iterations, residual, tolerance")) + eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym, scalar_type, SupportComponents, support_components, ncomponents, rowcomponent, colcomponent, SolverFailure, AbstractCoverStats, AugmentedLagrangianStats, ActiveSetStats, LeastSquaresStats, PowerMeanStats, ExternalSolverStats, MultistartStats, converged, iterations, residual, tolerance")) end include("penalties.jl") diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 39faffe..aa75bb8 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -223,20 +223,26 @@ function symcover_min(ϕ::AbsLinear, A::AbstractMatrix; strategies=SYMCOVER_MIN_ starts = [similar(Array{T}, ax) for _ in strategies] # Skip strategies that cannot produce a start for `A`. built = [_initialize_symcover!(a, A, strategy, :inflate) for (a, strategy) in zip(starts, strategies)] - covers = [symcover_min!(ϕ, a, A) for (a, ok) in zip(starts, built) if ok] - isempty(covers) && + any(built) || throw(ArgumentError("symcover_min: no strategy in $(string(strategies)) yields a starting cover of `A`")) - return covers[_multistart_select([cover_objective(ϕ, a, A) for a in covers])] + labels = Symbol[strategy for (strategy, ok) in zip(strategies, built) if ok] + covers = [a for (a, ok) in zip(starts, built) if ok] + stats = [last(symcover_min!(ϕ, a, A))::ExternalSolverStats for a in covers] + objs = [cover_objective(ϕ, a, A) for a in covers] + k = _multistart_select(objs) + return covers[k], MultistartStats{ExternalSolverStats,eltype(objs)}(labels, stats, objs, k, Symbol[]) end function cover_min(ϕ::AbsLinear, A::AbstractMatrix; strategies=COVER_MIN_STRATEGIES) isempty(strategies) && throw(ArgumentError("cover_min: `strategies` must name at least one starting cover")) + labels = Symbol[strategy for strategy in strategies] covers = [initialize_cover(A; strategy) for strategy in strategies] - for (a, b) in covers - cover_min!(ϕ, a, b, A) - end - return covers[_multistart_select([cover_objective(ϕ, a, b, A) for (a, b) in covers])] + stats = [last(cover_min!(ϕ, a, b, A))::ExternalSolverStats for (a, b) in covers] + objs = [cover_objective(ϕ, a, b, A) for (a, b) in covers] + k = _multistart_select(objs) + a, b = covers[k] + return a, b, MultistartStats{ExternalSolverStats,eltype(objs)}(labels, stats, objs, k, Symbol[]) end # ============================================================ @@ -2447,7 +2453,7 @@ Base.showerror(io::IO, e::SolverFailure) = print(io, "SolverFailure: ", e.msg) # Reject all non-solved statuses, including `ALMOST_*`. Taking the status keeps # this helper independent of JuMP. function check_solved(status, solver, fname) - Symbol(status) in (:OPTIMAL, :LOCALLY_SOLVED) || + Symbol(status) in SOLVED_STATUSES || throw(SolverFailure("$fname: $solver terminated with status $status")) return nothing end diff --git a/src/soft_covers.jl b/src/soft_covers.jl index 3d332a3..168c9ae 100644 --- a/src/soft_covers.jl +++ b/src/soft_covers.jl @@ -324,27 +324,26 @@ end # Shared multistart driver: refine every start, then select. `inits_builder` returns the # labels, the starts, and the number `nfixed` of leading deterministic starts; the rest -# are random. A random start whose refinement throws `SolverFailure` gets objective -# `Inf`, while a failing deterministic start propagates the error. Optional `labels` and -# `objs` collect candidate data for tests. -function _multistart_run(inits_builder::F, refine!::G, objective::H, A::AbstractMatrix, - starts::Int, σ::Real, rng; labels=nothing, objs=nothing) where {F,G,H} +# are random. `refine!` returns statistics of type `S`. A random start whose refinement +# throws `SolverFailure` is left out of the selection, while a failing deterministic +# start propagates the error. Returns the selected start and a `MultistartStats`. +function _multistart_run(::Type{S}, inits_builder::F, refine!::G, objective::H, A::AbstractMatrix, + starts::Int, σ::Real, rng) where {S,F,G,H} labs, inits, nfixed = inits_builder(A, starts, σ, rng) - failed = falses(length(inits)) + ok = trues(length(inits)) + stats = S[] for (k, x) in enumerate(inits) - k <= nfixed && (refine!(x, A); continue) try - refine!(x, A) + push!(stats, refine!(x, A)::S) catch err - err isa SolverFailure || rethrow() - failed[k] = true + (k <= nfixed || !(err isa SolverFailure)) && rethrow() + ok[k] = false end end - E = [objective(x, A) for x in inits] - E[failed] .= Inf - labels === nothing || append!(labels, labs) - objs === nothing || append!(objs, E) - return inits[_multistart_select(E)] + kept = inits[ok] + E = [objective(x, A) for x in kept] + sel = _multistart_select(E) + return kept[sel], MultistartStats{S,eltype(E)}(labs[ok], stats, E, sel, labs[.!ok]) end # Symmetric `AbsLinear` starts in selection order, followed by log-normal @@ -353,22 +352,22 @@ function _soft_symcover_abslinear_inits(A::AbstractMatrix, starts::Int, σ::Real ax = axes(A, 1) T = float(real(eltype(A))) # The soft cover imposes no coverage constraint, so the starts are taken raw - # (`feasible=:none`); "inflate" is the one candidate that is deliberately a cover. + # (`feasible=:none`); `:inflate` is the one candidate that is deliberately a cover. ag = initialize_symcover(A; strategy=:geomean, feasible=:none) # also the perturbation base - labels = ["geomean"] + labels = [:geomean] inits = [copy(ag)] length(inits) < starts && - (push!(labels, "hardcover"); push!(inits, initialize_symcover(A; strategy=:hardcover, feasible=:none))) + (push!(labels, :hardcover); push!(inits, initialize_symcover(A; strategy=:hardcover, feasible=:none))) length(inits) < starts && - (push!(labels, "inflate"); push!(inits, initialize_symcover(A; strategy=:geomean, feasible=:inflate))) + (push!(labels, :inflate); push!(inits, initialize_symcover(A; strategy=:geomean, feasible=:inflate))) if length(inits) < starts # `:leaveout` is unavailable when no entry can be dropped; a multistart forfeits that # slot rather than fail, so it takes the start through the gated builder. lo = similar(ag) - _initialize_symcover!(lo, A, :leaveout, :none) && (push!(labels, "leaveout"); push!(inits, lo)) + _initialize_symcover!(lo, A, :leaveout, :none) && (push!(labels, :leaveout); push!(inits, lo)) end if length(inits) < starts && _nsupport(A) < length(A) - push!(labels, "feasible"); push!(inits, initialize_symcover(A; strategy=:diagfeasible, feasible=:none)) + push!(labels, :feasible); push!(inits, initialize_symcover(A; strategy=:diagfeasible, feasible=:none)) end nfixed, k = length(inits), 0 while length(inits) < starts @@ -377,18 +376,17 @@ function _soft_symcover_abslinear_inits(A::AbstractMatrix, starts::Int, σ::Real ξ = randn(rng) p[i] = ag[i] > 0 ? ag[i] * exp(T(σ) * T(ξ)) : zero(T) end - k += 1; push!(labels, "rand$k"); push!(inits, p) + k += 1; push!(labels, Symbol(:rand, k)); push!(inits, p) end return labels, inits, nfixed end # Scale-covariant symmetric `AbsLinear` multistart. -function _soft_symcover_abslinear(ϕ::AbsLinear, A::AbstractMatrix, starts::Int, σ::Real, rng; - labels=nothing, objs=nothing) - return _multistart_run(_soft_symcover_abslinear_inits, - (a, A) -> soft_symcover!(ϕ, a, A), +function _soft_symcover_abslinear(ϕ::AbsLinear, A::AbstractMatrix, starts::Int, σ::Real, rng) + return _multistart_run(ExternalSolverStats, _soft_symcover_abslinear_inits, + (a, A) -> last(soft_symcover!(ϕ, a, A)), (a, A) -> cover_objective(ϕ, a, A), - A, starts, σ, rng; labels, objs) + A, starts, σ, rng) end # Asymmetric `AbsLinear` starts: boosted covariant start, tightened cover, and @@ -396,10 +394,10 @@ end function _soft_cover_abslinear_inits(A::AbstractMatrix, starts::Int, σ::Real, rng) T = float(real(eltype(A))) ag, bg = initialize_cover(A; strategy=:covariant, feasible=:boost) - labels = ["boost"] + labels = [:boost] inits = [(copy(ag), copy(bg))] # Reuse the covariant start and boost pass when constructing the hard start. - length(inits) < starts && (push!(labels, "hardcover"); push!(inits, tighten_cover!(copy(ag), copy(bg), A))) + length(inits) < starts && (push!(labels, :hardcover); push!(inits, tighten_cover!(copy(ag), copy(bg), A))) nfixed, k = length(inits), 0 while length(inits) < starts a = similar(ag); b = similar(bg) @@ -411,16 +409,16 @@ function _soft_cover_abslinear_inits(A::AbstractMatrix, starts::Int, σ::Real, r η = randn(rng) b[j] = bg[j] > 0 ? bg[j] * exp(T(σ) * T(η)) : zero(T) end - k += 1; push!(labels, "rand$k"); push!(inits, (a, b)) + k += 1; push!(labels, Symbol(:rand, k)); push!(inits, (a, b)) end return labels, inits, nfixed end # Scale-covariant asymmetric `AbsLinear` multistart. -function _soft_cover_abslinear(ϕ::AbsLinear, A::AbstractMatrix, starts::Int, σ::Real, rng; - labels=nothing, objs=nothing) - return _multistart_run(_soft_cover_abslinear_inits, - (ab, A) -> soft_cover!(ϕ, ab[1], ab[2], A), - (ab, A) -> cover_objective(ϕ, ab[1], ab[2], A), - A, starts, σ, rng; labels, objs) +function _soft_cover_abslinear(ϕ::AbsLinear, A::AbstractMatrix, starts::Int, σ::Real, rng) + (a, b), stats = _multistart_run(ExternalSolverStats, _soft_cover_abslinear_inits, + (ab, A) -> last(soft_cover!(ϕ, ab[1], ab[2], A)), + (ab, A) -> cover_objective(ϕ, ab[1], ab[2], A), + A, starts, σ, rng) + return a, b, stats end diff --git a/src/stats.jl b/src/stats.jl index d43242b..9aaca06 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -242,3 +242,102 @@ tolerance(s::PowerMeanStats) = s.tol PowerMeanStats{T}(nt::NamedTuple) where {T} = PowerMeanStats{T}(nt.converged, nt.tol, nt.imbalance, nt.nsweeps, nt.newton, nt.nnewton, nt.npasses, nt.nfactor, nt.linsolve) + +""" + ExternalSolverStats <: AbstractCoverStats + +Statistics of the covers computed by an external optimizer through JuMP: HiGHS +for `AbsLog{1}` ([`symcover_min`](@ref), [`cover_min`](@ref), +[`soft_symcover`](@ref), [`soft_cover`](@ref), and their mutating forms), and +Ipopt for one local refinement under `AbsLinear` ([`symcover_min!`](@ref), +[`cover_min!`](@ref), [`soft_symcover!`](@ref), and [`soft_cover!`](@ref)). + +A solve may consist of several stages, each an optimization of one model. Ipopt +runs one. HiGHS runs two unless `A` has no nonzero entries: a linear program +for the `AbsLog{1}` objective, then a quadratic program that selects the point +of least `AbsLog{2}` objective among its minimizers. The fields hold one entry +per stage, in order. + +A stage that does not end with a solved status throws a +[`SolverFailure`](@ref), so every status recorded here is a solved one and +`converged` returns `true`. The accessors `iterations`, `residual`, and +`tolerance` are not implemented. + +# Fields + +- `solver::Symbol`: `:HiGHS` or `:Ipopt`. +- `status::Vector{Symbol}`: the termination status of each stage, as the name of + the `MathOptInterface.TerminationStatusCode` (`:OPTIMAL` or `:LOCALLY_SOLVED`). +- `objective::Vector{Float64}`: the optimal value of each stage's model. The + models work in log scales and may include slack variables or constant terms, + so this need not equal [`cover_objective`](@ref) of the returned cover. +- `niters::Vector{Int}`: the iterations of each stage as the solver reports them + (for HiGHS, its simplex, interior-point, and QP iterations; for Ipopt, its + interior-point iterations). +- `solvetime::Vector{Float64}`: the solver's run time of each stage, in seconds. + +`niters` and `solvetime` may change meaning between releases. +""" +struct ExternalSolverStats <: AbstractCoverStats + solver::Symbol + status::Vector{Symbol} + objective::Vector{Float64} + niters::Vector{Int} + solvetime::Vector{Float64} +end + +# Termination statuses accepted as solved, by name, so that the core package +# does not depend on MathOptInterface. +const SOLVED_STATUSES = (:OPTIMAL, :LOCALLY_SOLVED) + +converged(s::ExternalSolverStats) = all(in(SOLVED_STATUSES), s.status) + +""" + MultistartStats{S,T} <: AbstractCoverStats + +Statistics of a cover selected from several local refinements: the `AbsLinear` +forms of [`symcover_min`](@ref) and [`cover_min`](@ref), which refine one start +per entry of `strategies`, and of [`soft_symcover`](@ref) and +[`soft_cover`](@ref), which refine deterministic and random starts. The cover +with the least objective is returned. + +A failed refinement of a random start (one that throws +[`SolverFailure`](@ref)) is left out of the selection and recorded in `failed`; +any other failure throws. The accessors [`converged`](@ref MatrixCovers.converged), +[`iterations`](@ref MatrixCovers.iterations), +[`residual`](@ref MatrixCovers.residual), and +[`tolerance`](@ref MatrixCovers.tolerance) return those of the selected start, +`stats[selected]`, where implemented for `S`. + +# Fields + +- `labels::Vector{Symbol}`: the starts whose refinement succeeded, in the order + they ran: the strategy names, or for the soft covers names such as `:geomean` + and `:rand1`. +- `stats::Vector{S}`: the statistics of each of these refinements. +- `objectives::Vector{T}`: the [`cover_objective`](@ref) of each refined cover. +- `selected::Int`: the index (into the vectors above) of the returned cover. +- `failed::Vector{Symbol}`: the labels of the random starts whose refinement + failed. +""" +struct MultistartStats{S<:AbstractCoverStats,T} <: AbstractCoverStats + labels::Vector{Symbol} + stats::Vector{S} + objectives::Vector{T} + selected::Int + failed::Vector{Symbol} + + function MultistartStats{S,T}(labels, stats, objectives, selected, failed) where {S,T} + length(labels) == length(stats) == length(objectives) || + throw(DimensionMismatch("labels, stats, and objectives must have equal lengths; got $(length(labels)), $(length(stats)), and $(length(objectives))")) + selected in eachindex(labels) || + throw(ArgumentError("selected=$selected is not an index of the $(length(labels)) refined starts")) + return new{S,T}(labels, stats, objectives, selected, failed) + end +end + +_selected(s::MultistartStats) = s.stats[s.selected] +converged(s::MultistartStats) = converged(_selected(s)) +iterations(s::MultistartStats) = iterations(_selected(s)) +residual(s::MultistartStats) = residual(_selected(s)) +tolerance(s::MultistartStats) = tolerance(_selected(s)) diff --git a/test/element_types.jl b/test/element_types.jl index cf77277..50a256a 100644 --- a/test/element_types.jl +++ b/test/element_types.jl @@ -4,8 +4,8 @@ @testset "Float32 flows through the family" begin A = Float32[4 1.5; 1.5 1] - for a in (symcover(A), first(soft_symcover(A)), soft_symcover(AbsLog{1}(), A), - soft_symcover(AbsLinear{1}(), A), first(symcover_min(AbsLog{2}(), A)), + for a in (symcover(A), first(soft_symcover(A)), first(soft_symcover(AbsLog{1}(), A)), + first(soft_symcover(AbsLinear{1}(), A)), first(symcover_min(AbsLog{2}(), A)), first(soft_symcover(AbsLog{2}(), A))) @test eltype(a) === Float32 @test all(isfinite, a) @@ -96,7 +96,7 @@ @test iscover(a, b, A; rtol=1e-6) for ϕ in (AbsLog{1}(), AbsLinear{1}(), AbsLinear{2}()) - x, y = scales(soft_cover!(ϕ, copy(a0), copy(b0), A)) + x, y, _ = soft_cover!(ϕ, copy(a0), copy(b0), A) @test eltype(x) === Float64 @test eltype(y) === Float32 @test all(isfinite, x) && all(isfinite, y) diff --git a/test/extensions.jl b/test/extensions.jl index 5b9aa96..4420c71 100644 --- a/test/extensions.jl +++ b/test/extensions.jl @@ -3,7 +3,7 @@ @testset "symcover_min and cover_min (JuMP/HiGHS)" begin for A in ([2.0 1.0; 1.0 3.0], [100.0 1.0; 1.0 0.01], [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0]) a_fast = symcover(A) - a_lmin = symcover_min(AbsLog{1}(), A) + a_lmin, _ = symcover_min(AbsLog{1}(), A) a_qmin, _ = symcover_min(AbsLog{2}(), A) # qmin is a valid cover @test iscover(a_qmin, A; atol=1e-10) @@ -14,7 +14,7 @@ end # Exact case with zeros A = [0 0 1; 0 0 2; 1 2 1] - a = symcover_min(AbsLog{1}(), A) + a, _ = symcover_min(AbsLog{1}(), A) @test iscover(a, A; atol=1e-10) @test a ≈ [1, 2, 1] @test abs(cover_objective(AbsLog{1}(), a, A)) < 1e-10 @@ -42,15 +42,15 @@ # soft_symcover(AbsLog{p}): unconstrained, lower objective than constrained for A in ([2.0 1.0; 1.0 3.0], [100.0 1.0; 1.0 0.01], [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0]) for ϕ in (AbsLog{1}(), AbsLog{2}()) - a_soft = scales(soft_symcover(ϕ, A)) - a_hard = scales(symcover_min(ϕ, A)) + a_soft, _ = soft_symcover(ϕ, A) + a_hard, _ = symcover_min(ϕ, A) @test cover_objective(ϕ, a_soft, A) <= cover_objective(ϕ, a_hard, A) + 1e-8 end end # Rank-1 matrix: both achieve zero objective A_rank1 = [2.0 1.0 4.0; 1.0 0.5 2.0; 4.0 2.0 8.0] @test cover_objective(AbsLog{2}(), first(soft_symcover(AbsLog{2}(), A_rank1)), A_rank1) < 1e-8 - @test cover_objective(AbsLog{1}(), soft_symcover(AbsLog{1}(), A_rank1), A_rank1) < 1e-6 + @test cover_objective(AbsLog{1}(), first(soft_symcover(AbsLog{1}(), A_rank1)), A_rank1) < 1e-6 # Non-optimal solver statuses are errors. @test MatrixCovers.check_solved(JuMP.OPTIMAL, "HiGHS", "symcover_min") === nothing @@ -74,32 +74,32 @@ end a_fast = symcover(AbsLinear{2}(), A) for ϕ in (AbsLinear{1}(), AbsLinear{2}()) # symcover_min: valid hard cover, at most as costly as heuristic - a_min = scales(symcover_min(ϕ, A)) + a_min, _ = symcover_min(ϕ, A) @test iscover(a_min, A; atol=1e-6) @test cover_objective(ϕ, a_min, A) <= cover_objective(ϕ, a_fast, A) + 1e-8 # soft_symcover: lower or equal objective than constrained version - a_soft = scales(soft_symcover(ϕ, A)) + a_soft, _ = soft_symcover(ϕ, A) @test cover_objective(ϕ, a_soft, A) <= cover_objective(ϕ, a_min, A) + 1e-8 end # AbsLinear{2} ≤ AbsLinear{1} soft objectives (p=2 is a stricter lower bound) - a2 = soft_symcover(AbsLinear{2}(), A) - a1 = soft_symcover(AbsLinear{1}(), A) + a2, _ = soft_symcover(AbsLinear{2}(), A) + a1, _ = soft_symcover(AbsLinear{1}(), A) @test cover_objective(AbsLinear{1}(), a1, A) <= cover_objective(AbsLinear{1}(), a2, A) + 1e-8 end # Rank-1 matrix: soft cover achieves near-zero objective for all nonzero entries A_rank1 = [4.0 2.0; 2.0 1.0] - a2 = soft_symcover(AbsLinear{2}(), A_rank1) + a2, _ = soft_symcover(AbsLinear{2}(), A_rank1) @test cover_objective(AbsLinear{2}(), a2, A_rank1) < 1e-8 - a1 = soft_symcover(AbsLinear{1}(), A_rank1) + a1, _ = soft_symcover(AbsLinear{1}(), A_rank1) @test cover_objective(AbsLinear{1}(), a1, A_rank1) < 1e-8 # Matrix with zeros: zero entries contribute 1 each regardless of cover A = [0.0 1.0; 1.0 0.0] # only off-diagonal nonzero; min possible = 0 (off-diag) + 2 (diag zeros) - a2 = soft_symcover(AbsLinear{2}(), A) + a2, _ = soft_symcover(AbsLinear{2}(), A) @test cover_objective(AbsLinear{2}(), a2, A) ≈ 2.0 atol=1e-8 - a1 = soft_symcover(AbsLinear{1}(), A) + a1, _ = soft_symcover(AbsLinear{1}(), A) @test cover_objective(AbsLinear{1}(), a1, A) ≈ 2.0 atol=1e-8 end @@ -117,8 +117,8 @@ end es = sort(map(1:12) do _ a = exp.(3 .* randn(rng, 5)) b = exp.(3 .* randn(rng, 5)) - sym ? cover_objective(ϕ, scales(soft_symcover!(ϕ, a, A)), A) : - cover_objective(ϕ, scales(soft_cover!(ϕ, a, b, A))..., A) + sym ? cover_objective(ϕ, first(soft_symcover!(ϕ, a, A)), A) : + cover_objective(ϕ, soft_cover!(ϕ, a, b, A)[1:2]..., A) end) gaps = diff(es) ./ es[begin:end-1] @test all(g -> g < 1e-8 || g > 1e-3, gaps) @@ -148,25 +148,25 @@ end for ϕ in PENALTIES # Refinement does not worsen the start beyond solver tolerance. a0 = initialize_symcover(A) - a = scales(symcover_min!(ϕ, copy(a0), A)) + a, _ = symcover_min!(ϕ, copy(a0), A) @test iscover(a, A; rtol=1e-6) @test cover_objective(ϕ, a, A) <= cover_objective(ϕ, a0, A) * (1 + 1e-6) + 1e-8 ab0, bb0 = initialize_cover(Aasym) - ab, bb = scales(cover_min!(ϕ, copy(ab0), copy(bb0), Aasym)) + ab, bb, _ = cover_min!(ϕ, copy(ab0), copy(bb0), Aasym) @test iscover(ab, bb, Aasym; rtol=1e-6) @test cover_objective(ϕ, ab, bb, Aasym) <= cover_objective(ϕ, ab0, bb0, Aasym) * (1 + 1e-6) + 1e-8 # The asymmetric result is pinned to the balance convention, so the start is # read only up to the gauge a -> c*a, b -> b/c that leaves a[i]*b[j] fixed. - ga, gb = scales(cover_min!(ϕ, 4 .* ab0, bb0 ./ 4, Aasym)) + ga, gb, _ = cover_min!(ϕ, 4 .* ab0, bb0 ./ 4, Aasym) @test ga ≈ ab && gb ≈ bb end # `AbsLog{2}` tie-breaking makes the `AbsLog{1}` result start-independent. - a_cold = symcover_min(AbsLog{1}(), A) + a_cold, _ = symcover_min(AbsLog{1}(), A) for strategy in (:hardcover, :geomean, :diagfeasible) - @test symcover_min!(AbsLog{1}(), initialize_symcover(A; strategy), A) ≈ a_cold + @test first(symcover_min!(AbsLog{1}(), initialize_symcover(A; strategy), A)) ≈ a_cold end ab_cold, bb_cold = cover_min(AbsLog{1}(), Aasym) for strategy in (:hardcover, :geomean) @@ -180,17 +180,17 @@ end 29.4700945830419 38.0236584552296 8.96405068596511 26.5775238859338; 0.0181293142917846 0.0279078887878805 26.5775238859338 42.6650094474717] ϕ = AbsLinear{2}() - a_hard = scales(symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin)) - a_geo = scales(symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:geomean), Abasin)) + a_hard, _ = symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin) + a_geo, _ = symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:geomean), Abasin) @test iscover(a_hard, Abasin; rtol=1e-6) @test iscover(a_geo, Abasin; rtol=1e-6) @test cover_objective(ϕ, a_geo, Abasin) < cover_objective(ϕ, a_hard, Abasin) - 0.5 # Offset axes survive the Ipopt position mapping. Ao = OffsetArray(A, -1, -1) - ao = scales(symcover_min!(ϕ, initialize_symcover(Ao), Ao)) + ao, _ = symcover_min!(ϕ, initialize_symcover(Ao), Ao) @test axes(ao, 1) == axes(Ao, 1) - @test collect(ao) ≈ scales(symcover_min!(ϕ, initialize_symcover(A), A)) + @test collect(ao) ≈ first(symcover_min!(ϕ, initialize_symcover(A), A)) end @testset "soft_symcover multistart and refiner (Ipopt)" begin @@ -199,17 +199,17 @@ end a0 = initialize_symcover(A; strategy=:geomean, feasible=:none) @test !iscover(a0, A) @test_throws "requires a start that covers `A`" symcover_min!(AbsLinear{2}(), copy(a0), A) - @test soft_symcover!(AbsLinear{2}(), copy(a0), A) isa AbstractVector + @test first(soft_symcover!(AbsLinear{2}(), copy(a0), A)) isa AbstractVector for ϕ in (AbsLinear{1}(), AbsLinear{2}(), AbsLog{1}()) - @test scales(soft_symcover(ϕ, A)) == scales(soft_symcover(ϕ, A)) # deterministic + @test first(soft_symcover(ϕ, A)) == first(soft_symcover(ϕ, A)) # deterministic end # No start on the menu beats the driver. for ϕ in (AbsLinear{1}(), AbsLinear{2}()) - a = scales(soft_symcover(ϕ, A)) + a, _ = soft_symcover(ϕ, A) for strategy in (:hardcover, :geomean, :leaveout) - single = scales(soft_symcover!(ϕ, initialize_symcover(A; strategy, feasible=:none), A)) + single, _ = soft_symcover!(ϕ, initialize_symcover(A; strategy, feasible=:none), A) @test cover_objective(ϕ, a, A) <= cover_objective(ϕ, single, A) * (1 + 1e-6) + 1e-8 end end @@ -222,7 +222,7 @@ end 3.0721768720180109 333.53567816710364 0.92571970793600067 7.3601096491681046 1.0844635712556003] d = [20.451338935482074, 0.69212569803171398, 35.401627522529395, 15.904906661932396, 0.19696509774727827] for ϕ in (AbsLinear{1}(), AbsLinear{2}()) - @test covaries(A -> scales(soft_symcover(ϕ, A)), Abasins, d; rtol=1e-5) + @test covaries(A -> first(soft_symcover(ϕ, A)), Abasins, d; rtol=1e-5) end @test_throws "positive scale on every supported row" soft_symcover!(AbsLinear{2}(), [1.0, -1.0, 2.0], A) @@ -241,30 +241,30 @@ end @test_throws "requires a start that covers `A`" cover_min!(AbsLinear{2}(), copy(a0), copy(b0), A) for ϕ in (AbsLinear{1}(), AbsLinear{2}(), AbsLog{1}(), AbsLog{2}()) - a, b = scales(soft_cover(ϕ, A)) - @test scales(soft_cover(ϕ, A)) == (a, b) # deterministic + a, b = soft_cover(ϕ, A) + @test soft_cover(ϕ, A)[1:2] == (a, b) # deterministic # The objective depends on `a`, `b` only through their products, so the split is # fixed by the balance convention rather than left to the solver. @test balance(a, b) ≈ 0 atol=1e-7 # Gauge-invariant in the start: (a, b) and (2a, b/2) name the same point. - ag, bg = scales(soft_cover!(ϕ, 2 .* copy(a0), copy(b0) ./ 2, A)) - ah, bh = scales(soft_cover!(ϕ, copy(a0), copy(b0), A)) + ag, bg, _ = soft_cover!(ϕ, 2 .* copy(a0), copy(b0) ./ 2, A) + ah, bh, _ = soft_cover!(ϕ, copy(a0), copy(b0), A) @test ag ≈ ah && bg ≈ bh end # No start on the menu beats the driver. for ϕ in (AbsLinear{1}(), AbsLinear{2}()) - a, b = scales(soft_cover(ϕ, A)) + a, b, _ = soft_cover(ϕ, A) for strategy in (:hardcover, :covariant) - sa, sb = scales(soft_cover!(ϕ, initialize_cover(A; strategy, feasible=:none)..., A)) + sa, sb, _ = soft_cover!(ϕ, initialize_cover(A; strategy, feasible=:none)..., A) @test cover_objective(ϕ, a, b, A) <= cover_objective(ϕ, sa, sb, A) * (1 + 1e-6) + 1e-8 end # The asymmetric soft cover relaxes the symmetric one on a symmetric matrix. As = [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0] - aa, ab = scales(soft_cover(ϕ, As)) + aa, ab, _ = soft_cover(ϕ, As) @test cover_objective(ϕ, aa, ab, As) <= - cover_objective(ϕ, scales(soft_symcover(ϕ, As)), As) * (1 + 1e-6) + 1e-8 + cover_objective(ϕ, first(soft_symcover(ϕ, As)), As) * (1 + 1e-6) + 1e-8 end # Covariant starts make the selected cover covariant on irregular support too. @@ -290,36 +290,36 @@ end B = exp.(3 .* randn(rng, n, n)) .* (rand(rng, n, n) .> 0.25) A = (B + B') / 2 all(iszero, A) && continue - a = scales(symcover_min(ϕ1, A)) + a, _ = symcover_min(ϕ1, A) @test iscover(a, A; rtol=1e-7) # No cover beats it on AbsLog{1}: the canonical selection does not cost optimality. # Any feasible cover is a witness; use the AbsLog{2}-minimal one, and the heuristic. - for other in (scales(symcover_min(ϕ2, A)), symcover(A)) + for other in (first(symcover_min(ϕ2, A)), symcover(A)) @test cover_objective(ϕ1, a, A) <= cover_objective(ϕ1, other, A) * (1 + 1e-6) + 1e-8 end # Start-independent in the *result*, not merely in the objective value. for strategy in (:hardcover, :geomean, :diagfeasible) - @test scales(symcover_min!(ϕ1, initialize_symcover(A; strategy), A)) ≈ a + @test first(symcover_min!(ϕ1, initialize_symcover(A; strategy), A)) ≈ a end # Scale-covariant: both objectives see `A` only through the residuals, which a # rescaling leaves invariant, so the selected member co-varies with the frame. D = Diagonal(exp.(randn(rng, n))) - @test scales(symcover_min(ϕ1, D * A * D)) ≈ D * a + @test first(symcover_min(ϕ1, D * A * D)) ≈ D * a end # Gauge balancing and `AbsLog{1}` face selection are independent. Aasym = [3.0 1.0 7.0; 2.0 5.0 1.0; 8.0 1.0 4.0] - a, b = scales(cover_min(ϕ1, Aasym)) + a, b, _ = cover_min(ϕ1, Aasym) nza = vec(count(!iszero, Aasym, dims=2)) nzb = vec(count(!iszero, Aasym, dims=1)) @test sum(nza .* log.(a)) ≈ sum(nzb .* log.(b)) atol=1e-8 @test iscover(a, b, Aasym; rtol=1e-7) # Gauge-invariant in the start, and start-independent in the result. a0, b0 = initialize_cover(Aasym; strategy=:geomean) - ag, bg = scales(cover_min!(ϕ1, 4 .* a0, b0 ./ 4, Aasym)) + ag, bg, _ = cover_min!(ϕ1, 4 .* a0, b0 ./ 4, Aasym) @test ag ≈ a && bg ≈ b end @@ -332,43 +332,43 @@ end Aasym = [1.0 2.0 3.0; 40.0 5.0 0.6] for ϕ in (AbsLinear{1}(), AbsLinear{2}()) - a = scales(symcover_min(ϕ, Abasin)) + a, _ = symcover_min(ϕ, Abasin) @test iscover(a, Abasin; rtol=1e-6) # No start on the menu beats the driver. for strategy in (:hardcover, :geomean, :leaveout) - single = scales(symcover_min!(ϕ, initialize_symcover(Abasin; strategy), Abasin)) + single, _ = symcover_min!(ϕ, initialize_symcover(Abasin; strategy), Abasin) @test cover_objective(ϕ, a, Abasin) <= cover_objective(ϕ, single, Abasin) * (1 + 1e-6) + 1e-8 end - @test scales(symcover_min(ϕ, Abasin)) == a # deterministic + @test first(symcover_min(ϕ, Abasin)) == a # deterministic - ab, bb = scales(cover_min(ϕ, Aasym)) + ab, bb, _ = cover_min(ϕ, Aasym) @test iscover(ab, bb, Aasym; rtol=1e-6) for strategy in (:hardcover, :covariant) - sa, sb = scales(cover_min!(ϕ, initialize_cover(Aasym; strategy)..., Aasym)) + sa, sb, _ = cover_min!(ϕ, initialize_cover(Aasym; strategy)..., Aasym) @test cover_objective(ϕ, ab, bb, Aasym) <= cover_objective(ϕ, sa, sb, Aasym) * (1 + 1e-6) + 1e-8 end - @test scales(cover_min(ϕ, Aasym)) == (ab, bb) # deterministic + @test cover_min(ϕ, Aasym)[1:2] == (ab, bb) # deterministic # The asymmetric cover relaxes the symmetric one: independent row and column # scales can only do better on the same matrix. - asym_a, asym_b = scales(cover_min(ϕ, Abasin)) + asym_a, asym_b, _ = cover_min(ϕ, Abasin) @test cover_objective(ϕ, asym_a, asym_b, Abasin) <= - cover_objective(ϕ, scales(symcover_min(ϕ, Abasin)), Abasin) * (1 + 1e-6) + 1e-8 + cover_objective(ϕ, first(symcover_min(ϕ, Abasin)), Abasin) * (1 + 1e-6) + 1e-8 end # Restricting `strategies` changes the selected basin. ϕ = AbsLinear{2}() - @test scales(symcover_min(ϕ, Abasin; strategies=(:hardcover,))) ≈ - scales(symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin)) - @test cover_objective(ϕ, scales(symcover_min(ϕ, Abasin)), Abasin) < - cover_objective(ϕ, scales(symcover_min(ϕ, Abasin; strategies=(:hardcover,))), Abasin) - 0.5 + @test first(symcover_min(ϕ, Abasin; strategies=(:hardcover,))) ≈ + first(symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin)) + @test cover_objective(ϕ, first(symcover_min(ϕ, Abasin)), Abasin) < + cover_objective(ϕ, first(symcover_min(ϕ, Abasin; strategies=(:hardcover,))), Abasin) - 0.5 # A matrix whose every row carries a single support entry admits no :leaveout start; # that strategy forfeits its slot rather than failing the solve. Adiag = [4.0 0.0; 0.0 9.0] - @test scales(symcover_min(ϕ, Adiag)) ≈ [2.0, 3.0] rtol=1e-4 + @test first(symcover_min(ϕ, Adiag)) ≈ [2.0, 3.0] rtol=1e-4 @test_throws "unknown strategy :banana" symcover_min(ϕ, Abasin; strategies=(:banana,)) @test_throws "no strategy in (:leaveout,)" symcover_min(ϕ, Adiag; strategies=(:leaveout,)) @@ -478,8 +478,8 @@ end @test first(symcover_min(AbsLog{2}(), M)) ≈ first(symcover_min(AbsLog{2}(), dense)) rtol=1e-8 @test MatrixCovers.symcover_min_jump(AbsLog{2}(), M) ≈ MatrixCovers.symcover_min_jump(AbsLog{2}(), dense) rtol=1e-6 - @test symcover_min(AbsLog{1}(), M) ≈ symcover_min(AbsLog{1}(), dense) rtol=1e-6 - @test iscover(symcover_min(AbsLog{1}(), M), dense; rtol=1e-8) + @test first(symcover_min(AbsLog{1}(), M)) ≈ first(symcover_min(AbsLog{1}(), dense)) rtol=1e-6 + @test iscover(first(symcover_min(AbsLog{1}(), M)), dense; rtol=1e-8) for (fh, fd) in ((cover_min(AbsLog{1}(), M), cover_min(AbsLog{1}(), dense)), (MatrixCovers.cover_min_jump(AbsLog{2}(), M), @@ -497,20 +497,20 @@ end # Ipopt: the AbsLinear kernels, whose starts the caller supplies. for p in (1, 2) ϕ = AbsLinear{p}() - sh = scales(symcover_min!(ϕ, symcover(ϕ, M), M)) - sd = scales(symcover_min!(ϕ, symcover(ϕ, dense), dense)) + sh, _ = symcover_min!(ϕ, symcover(ϕ, M), M) + sd, _ = symcover_min!(ϕ, symcover(ϕ, dense), dense) @test sh ≈ sd rtol=1e-5 - ah, bh = scales(cover_min!(ϕ, cover(ϕ, M)..., M)) - ad, bd = scales(cover_min!(ϕ, cover(ϕ, dense)..., dense)) + ah, bh, _ = cover_min!(ϕ, cover(ϕ, M)..., M) + ad, bd, _ = cover_min!(ϕ, cover(ϕ, dense)..., dense) @test ah ≈ ad rtol=1e-5 @test bh ≈ bd rtol=1e-5 - @test scales(soft_symcover!(ϕ, scales(soft_symcover(ϕ, M)), M)) ≈ - scales(soft_symcover!(ϕ, scales(soft_symcover(ϕ, dense)), dense)) rtol=1e-5 + @test first(soft_symcover!(ϕ, first(soft_symcover(ϕ, M)), M)) ≈ + first(soft_symcover!(ϕ, first(soft_symcover(ϕ, dense)), dense)) rtol=1e-5 - sah, sbh = scales(soft_cover!(ϕ, scales(soft_cover(ϕ, M))..., M)) - sad, sbd = scales(soft_cover!(ϕ, scales(soft_cover(ϕ, dense))..., dense)) + sah, sbh, _ = soft_cover!(ϕ, soft_cover(ϕ, M)[1:2]..., M) + sad, sbd, _ = soft_cover!(ϕ, soft_cover(ϕ, dense)[1:2]..., dense) @test sah ≈ sad rtol=1e-5 @test sbh ≈ sbd rtol=1e-5 end @@ -541,9 +541,98 @@ end @testset "the Ipopt sym objective uses the full-grid weighting" begin A = Float64[4 1 0; 1 1 5; 0 5 2] a0 = fill(sqrt(maximum(A)), 3) - a = symcover_min!(AbsLinear{2}(), copy(a0), A) + a, _ = symcover_min!(AbsLinear{2}(), copy(a0), A) @test a ≈ [2.0, 1.0, 5.0] rtol=1e-6 @test iscover(a, A; rtol=1e-6) # The objective reported for the result is the one that was minimized. @test cover_objective(AbsLinear{2}(), a, A) ≈ cover_objective(AbsLinear{2}(), [2.0, 1.0, 5.0], A) rtol=1e-6 end + +@testset "solver statistics (JuMP/HiGHS/Ipopt)" begin + ESS = MatrixCovers.ExternalSolverStats + A = [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0] + Aasym = [1.0 2.0 3.0; 4.0 5.0 6.0] + + # HiGHS: the AbsLog{1} LP, then the AbsLog{2} QP over its optimal face. + a, stats = @inferred symcover_min(AbsLog{1}(), A) + @test stats isa ESS + @test stats.solver === :HiGHS + @test stats.status == [:OPTIMAL, :OPTIMAL] + @test length(stats.objective) == length(stats.niters) == length(stats.solvetime) == 2 + @test all(>=(0), stats.niters) && all(>=(0), stats.solvetime) + @test stats.objective[2] ≈ cover_objective(AbsLog{2}(), a, A) rtol=1e-6 + @test MatrixCovers.converged(stats) + @test_throws MethodError MatrixCovers.iterations(stats) + @test_throws MethodError MatrixCovers.residual(stats) + @test_throws MethodError MatrixCovers.tolerance(stats) + for r in (@inferred(cover_min(AbsLog{1}(), Aasym)), @inferred(soft_cover(AbsLog{1}(), Aasym)), + @inferred(soft_symcover(AbsLog{1}(), A)), + @inferred(symcover_min!(AbsLog{1}(), symcover(A), A)), + @inferred(soft_symcover!(AbsLog{1}(), symcover(A), A)), + @inferred(cover_min!(AbsLog{1}(), cover(Aasym)..., Aasym)), + @inferred(soft_cover!(AbsLog{1}(), cover(Aasym)..., Aasym))) + s = last(r) + @test s isa ESS && s.solver === :HiGHS && s.status == [:OPTIMAL, :OPTIMAL] + end + # Without support there is no optimal face to search, so only the LP runs. + _, _, s0 = soft_cover(AbsLog{1}(), zeros(2, 3)) + @test s0.status == [:OPTIMAL] + + # Ipopt: one local refinement. + a, stats = @inferred symcover_min!(AbsLinear{2}(), symcover(A), A) + @test stats isa ESS + @test stats.solver === :Ipopt + @test stats.status == [:LOCALLY_SOLVED] + @test only(stats.niters) > 0 + @test only(stats.objective) ≈ cover_objective(AbsLinear{2}(), a, A) rtol=1e-8 + @test MatrixCovers.converged(stats) + @test_throws MethodError MatrixCovers.iterations(stats) + for r in (@inferred(cover_min!(AbsLinear{1}(), cover(Aasym)..., Aasym)), + @inferred(soft_symcover!(AbsLinear{1}(), symcover(A), A)), + @inferred(soft_cover!(AbsLinear{2}(), cover(Aasym)..., Aasym))) + s = last(r) + @test s isa ESS && s.solver === :Ipopt && s.status == [:LOCALLY_SOLVED] + end + + # Strategy multistart: one refinement per strategy that yields a start. + Abasin = [81.892035218799 1.06622031288736 29.4700945830419 0.0181293142917846; + 1.06622031288736 0.243512973596586 38.0236584552296 0.0279078887878805; + 29.4700945830419 38.0236584552296 8.96405068596511 26.5775238859338; + 0.0181293142917846 0.0279078887878805 26.5775238859338 42.6650094474717] + ϕ = AbsLinear{2}() + a, ms = @inferred symcover_min(ϕ, Abasin) + @test ms isa MatrixCovers.MultistartStats{ESS,Float64} + @test ms.labels == [:hardcover, :geomean, :leaveout] + @test ms.labels[ms.selected] === :geomean + @test isempty(ms.failed) + @test all(s -> s.solver === :Ipopt, ms.stats) + @test ms.objectives[ms.selected] == cover_objective(ϕ, a, Abasin) + @test ms.objectives[ms.selected] <= minimum(ms.objectives) * (1 + 1e-6) + # Each record describes the refinement of its own start. + _, s_hard = symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin) + @test ms.stats[1].objective ≈ s_hard.objective + @test MatrixCovers.converged(ms) === MatrixCovers.converged(ms.stats[ms.selected]) + @test_throws MethodError MatrixCovers.iterations(ms) + # A strategy that yields no start for `A` is not recorded. + _, ms = symcover_min(ϕ, [4.0 0.0; 0.0 9.0]) + @test ms.labels == [:hardcover, :geomean] + a, b, ms = @inferred cover_min(AbsLinear{1}(), Aasym) + @test ms.labels == [:hardcover, :covariant] + @test ms.objectives[ms.selected] == cover_objective(AbsLinear{1}(), a, b, Aasym) + + # Random multistart: failed random starts are recorded and left out. + a, b, ms = @inferred soft_cover(ϕ, ones(3, 3); starts=8, σ=10.0, rng=StableRNG(2)) + @test !isempty(ms.failed) + @test all(l -> startswith(string(l), "rand"), ms.failed) + @test isempty(intersect(ms.labels, ms.failed)) + @test length(ms.labels) + length(ms.failed) == 8 + @test length(ms.stats) == length(ms.objectives) == length(ms.labels) + @test ms.objectives[ms.selected] == cover_objective(ϕ, a, b, ones(3, 3)) + a, ms = @inferred soft_symcover(ϕ, Float32.(A)) + @test ms isa MatrixCovers.MultistartStats{ESS,Float32} + @test ms.labels[1:2] == [:geomean, :hardcover] + + # The constructor checks that its records line up. + @test_throws DimensionMismatch MatrixCovers.MultistartStats{ESS,Float64}([:a, :b], [stats], [1.0], 1, Symbol[]) + @test_throws "is not an index" MatrixCovers.MultistartStats{ESS,Float64}([:a], [stats], [1.0], 2, Symbol[]) +end diff --git a/test/invariants.jl b/test/invariants.jl index 5ac3a10..45aed66 100644 --- a/test/invariants.jl +++ b/test/invariants.jl @@ -9,13 +9,13 @@ const SYM_NOTIONS = ( (name = "symcover", f = A -> symcover(A), hard = true, rtol = 1e-9), (name = "symcover_min(AbsLog{2})", f = A -> first(symcover_min(AbsLog{2}(), A)), hard = true, rtol = 1e-5), - (name = "symcover_min(AbsLog{1})", f = A -> symcover_min(AbsLog{1}(), A), hard = true, rtol = 1e-5), - (name = "symcover_min(AbsLinear{1})", f = A -> symcover_min(AbsLinear{1}(), A), hard = true, rtol = 1e-5), - (name = "symcover_min(AbsLinear{2})", f = A -> symcover_min(AbsLinear{2}(), A), hard = true, rtol = 1e-5), + (name = "symcover_min(AbsLog{1})", f = A -> first(symcover_min(AbsLog{1}(), A)), hard = true, rtol = 1e-5), + (name = "symcover_min(AbsLinear{1})", f = A -> first(symcover_min(AbsLinear{1}(), A)), hard = true, rtol = 1e-5), + (name = "symcover_min(AbsLinear{2})", f = A -> first(symcover_min(AbsLinear{2}(), A)), hard = true, rtol = 1e-5), (name = "soft_symcover(AbsLog{2})", f = A -> first(soft_symcover(AbsLog{2}(), A)), hard = false, rtol = 1e-9), - (name = "soft_symcover(AbsLog{1})", f = A -> soft_symcover(AbsLog{1}(), A), hard = false, rtol = 1e-5), - (name = "soft_symcover(AbsLinear{2})", f = A -> soft_symcover(AbsLinear{2}(), A), hard = false, rtol = 1e-5), - (name = "soft_symcover(AbsLinear{1})", f = A -> soft_symcover(AbsLinear{1}(), A), hard = false, rtol = 1e-5), + (name = "soft_symcover(AbsLog{1})", f = A -> first(soft_symcover(AbsLog{1}(), A)), hard = false, rtol = 1e-5), + (name = "soft_symcover(AbsLinear{2})", f = A -> first(soft_symcover(AbsLinear{2}(), A)), hard = false, rtol = 1e-5), + (name = "soft_symcover(AbsLinear{1})", f = A -> first(soft_symcover(AbsLinear{1}(), A)), hard = false, rtol = 1e-5), (name = "soft_symcover(PowerMean{2})", f = A -> first(soft_symcover(PowerMean{2}(), A)), hard = false, rtol = 1e-12), (name = "soft_symcover(PowerMean{1})", f = A -> first(soft_symcover(PowerMean{1}(), A)), hard = false, rtol = 1e-12), ) @@ -23,13 +23,13 @@ const SYM_NOTIONS = ( const GEN_NOTIONS = ( (name = "cover", f = A -> cover(A), hard = true, rtol = 1e-9), (name = "cover_min(AbsLog{2})", f = A -> scales(cover_min(AbsLog{2}(), A)), hard = true, rtol = 1e-5), - (name = "cover_min(AbsLog{1})", f = A -> cover_min(AbsLog{1}(), A), hard = true, rtol = 1e-5), - (name = "cover_min(AbsLinear{1})", f = A -> cover_min(AbsLinear{1}(), A), hard = true, rtol = 1e-5), - (name = "cover_min(AbsLinear{2})", f = A -> cover_min(AbsLinear{2}(), A), hard = true, rtol = 1e-5), + (name = "cover_min(AbsLog{1})", f = A -> scales(cover_min(AbsLog{1}(), A)), hard = true, rtol = 1e-5), + (name = "cover_min(AbsLinear{1})", f = A -> scales(cover_min(AbsLinear{1}(), A)), hard = true, rtol = 1e-5), + (name = "cover_min(AbsLinear{2})", f = A -> scales(cover_min(AbsLinear{2}(), A)), hard = true, rtol = 1e-5), (name = "soft_cover(AbsLog{2})", f = A -> scales(soft_cover(AbsLog{2}(), A)), hard = false, rtol = 1e-9), - (name = "soft_cover(AbsLog{1})", f = A -> soft_cover(AbsLog{1}(), A), hard = false, rtol = 1e-5), - (name = "soft_cover(AbsLinear{2})", f = A -> soft_cover(AbsLinear{2}(), A), hard = false, rtol = 1e-5), - (name = "soft_cover(AbsLinear{1})", f = A -> soft_cover(AbsLinear{1}(), A), hard = false, rtol = 1e-5), + (name = "soft_cover(AbsLog{1})", f = A -> scales(soft_cover(AbsLog{1}(), A)), hard = false, rtol = 1e-5), + (name = "soft_cover(AbsLinear{2})", f = A -> scales(soft_cover(AbsLinear{2}(), A)), hard = false, rtol = 1e-5), + (name = "soft_cover(AbsLinear{1})", f = A -> scales(soft_cover(AbsLinear{1}(), A)), hard = false, rtol = 1e-5), (name = "soft_cover(PowerMean{2})", f = A -> scales(soft_cover(PowerMean{2}(), A)), hard = false, rtol = 1e-12), (name = "soft_cover(PowerMean{3})", f = A -> scales(soft_cover(PowerMean{3}(), A)), hard = false, rtol = 1e-12), ) @@ -168,7 +168,7 @@ const GEN_NOTIONS = ( @test cover_objective(AbsLog{2}(), MatrixCovers.symcover_min_jump(AbsLog{2}(), M), M) ≈ o2 rtol=1e-5 # AbsLog{1}: the HiGHS LP, whose objective is assembled from nonzero counts # rather than from the residuals directly. - @test cover_objective(AbsLog{1}(), symcover_min(AbsLog{1}(), M), M) ≈ ref_min(1, M) rtol=1e-5 + @test cover_objective(AbsLog{1}(), first(symcover_min(AbsLog{1}(), M)), M) ≈ ref_min(1, M) rtol=1e-5 end end end diff --git a/test/iscover.jl b/test/iscover.jl index a35e447..aa4c808 100644 --- a/test/iscover.jl +++ b/test/iscover.jl @@ -87,6 +87,6 @@ # rather than forbidding it, so whether one covers is a real question. Asoft = [1.0 4.0; 4.0 1.0] @test iscover(symcover(Asoft), Asoft; rtol=8eps()) - @test !iscover(soft_symcover(AbsLinear{2}(), Asoft), Asoft) + @test !iscover(first(soft_symcover(AbsLinear{2}(), Asoft)), Asoft) end end diff --git a/test/soft_covers.jl b/test/soft_covers.jl index a1d187f..1c5170b 100644 --- a/test/soft_covers.jl +++ b/test/soft_covers.jl @@ -17,7 +17,7 @@ end @test abs(residual) < 1e-10 end # AbsLinear{2}: at the minimum, ∑_j (1 - r_{kj}) r_{kj} ≈ 0 for each k - a = soft_symcover(AbsLinear{2}(), A) + a, _ = soft_symcover(AbsLinear{2}(), A) for k in axes(A, 1) residual = sum((rj = abs(A[k,j])/(a[k]*a[j]); (1 - rj) * rj) for j in axes(A, 2)) @test abs(residual) < 1e-10 @@ -33,7 +33,7 @@ end # Rank-1 is exact: A = [4 2; 2 1] = [2;1]*[2 1], so exact cover has all ratios = 1 A = [4.0 2.0; 2.0 1.0] for ϕ in SOFT_PENALTIES - a = scales(soft_symcover(ϕ, A)) + a, _ = soft_symcover(ϕ, A) @test cover_objective(ϕ, a, A) ≈ 0.0 atol=1e-8 end @@ -41,13 +41,13 @@ end A = [2.0 1.0; 1.0 3.0] d = [2.0, 0.5] for ϕ in SOFT_PENALTIES - @test covaries(a -> scales(soft_symcover(ϕ, a)), A, d; rtol=1e-5) + @test covaries(a -> first(soft_symcover(ϕ, a)), A, d; rtol=1e-5) end # Lower objective than naive uniform scaling perturbations A = [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0] for ϕ in SOFT_PENALTIES - a = scales(soft_symcover(ϕ, A)) + a, _ = soft_symcover(ϕ, A) E = cover_objective(ϕ, a, A) for scale in [0.5, 0.9, 1.1, 2.0] @test E <= cover_objective(ϕ, scale .* a, A) + 1e-8 @@ -59,8 +59,8 @@ end A_zero = [γ 1.0; 1.0 0.0] A_small = [γ 1.0; 1.0 1e-10] for ϕ in (AbsLinear{1}(), AbsLinear{2}()) - a_zero = scales(soft_symcover(ϕ, A_zero)) - a_small = scales(soft_symcover(ϕ, A_small)) + a_zero, _ = soft_symcover(ϕ, A_zero) + a_small, _ = soft_symcover(ϕ, A_small) @test a_small ≈ a_zero atol=1e-5 end @@ -68,8 +68,8 @@ end for ϕ in (AbsLinear{1}(), AbsLinear{2}()) for (B, d) in ((A_small, [50.0, 0.02]), ([3.0 7.6e-10; 7.6e-10 80.0], [35.0, 3400.0])) - a0 = scales(soft_symcover(ϕ, B)) - as = scales(soft_symcover(ϕ, B .* d .* d')) + a0, _ = soft_symcover(ϕ, B) + as, _ = soft_symcover(ϕ, B .* d .* d') @test as .* as' ≈ (d .* d') .* (a0 .* a0') rtol=1e-6 end end @@ -130,7 +130,7 @@ end @test a isa Vector{Float64} && b isa Vector{Float64} # Determinism (the default rng is evaluated fresh per call). - @test soft_cover(AbsLinear{1}(), A1) == soft_cover(AbsLinear{1}(), A1) + @test soft_cover(AbsLinear{1}(), A1)[1:2] == soft_cover(AbsLinear{1}(), A1)[1:2] # The returned objective never exceeds that of the geometric-mean start. for M in ([1.0 2.0 3.0; 6.0 5.0 4.0], @@ -175,7 +175,7 @@ end @test a isa Vector{Float64} && b isa Vector{Float64} end S = [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] - a = soft_symcover(AbsLog{1}(), S) + a, _ = soft_symcover(AbsLog{1}(), S) E = cover_objective(AbsLog{1}(), a, S) @test E <= cover_objective(AbsLog{1}(), first(soft_symcover(AbsLog{2}(), S)), S) + 1e-12 for _ in 1:20 @@ -184,7 +184,7 @@ end # The L1 minimizer of a path is a whole face; the AbsLog{2} tie-break picks one # point, and bipartite support is balanced like the other soft covers. - @test soft_symcover(AbsLog{1}(), [0 4.0; 4.0 0]) ≈ [2.0, 2.0] + @test first(soft_symcover(AbsLog{1}(), [0 4.0; 4.0 0])) ≈ [2.0, 2.0] # A rank-1 matrix is exactly coverable, so the L1 objective reaches 0 to within # the LP solver's feasibility tolerance. @@ -193,7 +193,7 @@ end @test cover_objective(AbsLog{1}(), a, b, A1) ≈ 0 atol=1e-6 # Deterministic, and scale-covariant in the product under row/column rescaling. - @test soft_cover(AbsLog{1}(), A1) == soft_cover(AbsLog{1}(), A1) + @test soft_cover(AbsLog{1}(), A1)[1:2] == soft_cover(AbsLog{1}(), A1)[1:2] 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) @@ -208,7 +208,8 @@ end 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 + s1, _ = soft_symcover(AbsLog{1}(), S1) + @test a .* b' ≈ s1 .* s1' rtol=1e-8 # The refiners do not depend on the start. a0, b0 = initialize_cover(Ac; strategy=:geomean, feasible=:none) @@ -216,17 +217,17 @@ end ae, be = soft_cover(AbsLog{1}(), Ac) @test ar .* br' ≈ ae .* be' rtol=1e-8 s0 = initialize_symcover(S; strategy=:geomean, feasible=:none) - @test soft_symcover!(AbsLog{1}(), 2 .* s0, S) ≈ soft_symcover(AbsLog{1}(), S) rtol=1e-8 + @test first(soft_symcover!(AbsLog{1}(), 2 .* s0, S)) ≈ first(soft_symcover(AbsLog{1}(), S)) rtol=1e-8 end @testset "AbsLinear soft-cover multistart" begin A = [1.0 2.0 3.0; 6.0 5.0 4.0] As = [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0] sc(A; kwargs...) = soft_cover(AbsLinear{2}(), A; kwargs...) - ssc(A; kwargs...) = soft_symcover(AbsLinear{2}(), A; kwargs...) + ssc(A; kwargs...) = first(soft_symcover(AbsLinear{2}(), A; kwargs...)) # A caller-supplied `rng` is threaded through: identical RNG state ⇒ identical result. - @test sc(A; rng=StableRNG(7)) == sc(A; rng=StableRNG(7)) + @test sc(A; rng=StableRNG(7))[1:2] == sc(A; rng=StableRNG(7))[1:2] @test ssc(As; rng=StableRNG(7)) == ssc(As; rng=StableRNG(7)) # `starts` and `σ` are accepted; `starts=1` is a plain single start. @@ -235,9 +236,9 @@ end # `sigma` is an ASCII alias for `σ`; passing both is fine when they agree and an # error when they don't. - @test sc(A; starts=1, sigma=1.5) == sc(A; starts=1, σ=1.5) + @test sc(A; starts=1, sigma=1.5)[1:2] == sc(A; starts=1, σ=1.5)[1:2] @test ssc(As; starts=1, sigma=1.5) == ssc(As; starts=1, σ=1.5) - @test sc(A; starts=1, σ=1.5, sigma=1.5) == sc(A; starts=1, σ=1.5) + @test sc(A; starts=1, σ=1.5, sigma=1.5)[1:2] == sc(A; starts=1, σ=1.5)[1:2] @test_throws "specify only one" sc(A; σ=1.5, sigma=2.0) @test_throws "specify only one" ssc(As; σ=1.5, sigma=2.0) @@ -288,65 +289,69 @@ end end @testset "feasible start and provenance" begin - # Instrumented multistart exposes candidate labels and objectives. + # The multistart statistics expose candidate labels and objectives. function provenance(A; rng=StableRNG(0)) - labels = String[]; objs = Float64[] - a = MatrixCovers._soft_symcover_abslinear(AbsLinear{2}(), A, 5, 2.0, rng; labels, objs) - return a, labels[MatrixCovers._multistart_select(objs)], labels, objs + a, stats = soft_symcover(AbsLinear{2}(), A; rng) + return a, stats.labels[stats.selected], stats.labels, stats.objectives end # The feasible start reaches the better basin for this matrix with zeros. Afe = Float64[0 11 18 0 12; 11 0 1 0 20; 18 1 0 3 18; 0 0 3 18 0; 12 20 18 0 17] a, winner, labels, objs = provenance(Afe) - @test winner == "feasible" - # The instrumented call returns exactly what the public entry point selects. - @test a == soft_symcover(AbsLinear{2}(), Afe; rng=StableRNG(0)) + @test winner === :feasible # The winning basin is separated beyond descent tolerance. - fi = findfirst(==("feasible"), labels) + fi = findfirst(==(:feasible), labels) @test objs[fi] <= minimum(objs) * (1 + 1e-6) other_basin = minimum(o for o in objs if o > objs[fi] * (1 + 1e-6)) @test objs[fi] < other_basin * (1 - 1e-3) # Selecting the `feasible` start preserves scale-covariance of the returned cover. d = [1.5, 0.3, 4.0, 0.7, 2.2] - a_scaled = soft_symcover(AbsLinear{2}(), Afe .* d .* d'; rng=StableRNG(0)) + a_scaled, _ = soft_symcover(AbsLinear{2}(), Afe .* d .* d'; rng=StableRNG(0)) @test a_scaled * a_scaled' ≈ (d .* d') .* (a * a') rtol=1e-7 # Gate off: a fully dense matrix (no zeros) never offers the `feasible` start. _, _, dense_labels, _ = provenance(Float64[4 1 2; 1 5 3; 2 3 6]) - @test "feasible" ∉ dense_labels + @test :feasible ∉ dense_labels # The asymmetric multistart exposes provenance the same way (no `feasible` start there). - # As above, the positional arguments mirror `soft_cover`'s defaults. Ag = [1.0 2.0 3.0; 6.0 5.0 4.0] - albs = String[]; aobjs = Float64[] - ab = MatrixCovers._soft_cover_abslinear(AbsLinear{2}(), Ag, 4, 2.0, StableRNG(0); - labels=albs, objs=aobjs) - @test ab == soft_cover(AbsLinear{2}(), Ag; rng=StableRNG(0)) - @test albs[MatrixCovers._multistart_select(aobjs)] in albs - @test "feasible" ∉ albs + _, _, astats = soft_cover(AbsLinear{2}(), Ag; rng=StableRNG(0)) + @test astats.labels == [:boost, :hardcover, :rand1, :rand2] + @test astats.selected == MatrixCovers._multistart_select(astats.objectives) + @test :feasible ∉ astats.labels end @testset "multistart skips failed random starts" begin # Two deterministic starts, then two random ones; `refine!` fails on the labeled start. - builder(A, starts, σ, rng) = (["s1", "s2", "rand1", "rand2"], [[1.0], [2.0], [3.0], [0.5]], 2) + ESS = MatrixCovers.ExternalSolverStats + builder(A, starts, σ, rng) = ([:s1, :s2, :rand1, :rand2], [[1.0], [2.0], [3.0], [0.5]], 2) objective(x, A) = x[1] - function run(bad) - refine!(x, A) = x[1] == bad && throw(MatrixCovers.SolverFailure("refinement failed")) - labels = String[]; objs = Float64[] - x = MatrixCovers._multistart_run(builder, refine!, objective, ones(1, 1), 4, 2.0, nothing; - labels, objs) - return x, objs + function refine!(x, A, bad) + x[1] == bad && throw(MatrixCovers.SolverFailure("refinement failed")) + return ESS(:Ipopt, [:LOCALLY_SOLVED], [x[1]], [1], [0.0]) end + run(bad) = MatrixCovers._multistart_run(ESS, builder, (x, A) -> refine!(x, A, bad), objective, + ones(1, 1), 4, 2.0, nothing) # A failed random start is excluded from the selection, even when its start is best. - x, objs = run(0.5) + x, stats = @inferred run(0.5) @test x == [1.0] - @test objs == [1.0, 2.0, 3.0, Inf] + @test stats isa MatrixCovers.MultistartStats{ESS,Float64} + @test stats.labels == [:s1, :s2, :rand1] + @test stats.objectives == [1.0, 2.0, 3.0] + @test [only(s.objective) for s in stats.stats] == [1.0, 2.0, 3.0] + @test stats.selected == 1 + @test stats.failed == [:rand2] + @test MatrixCovers.converged(stats) + # Without failures, `failed` is empty. + _, stats = run(NaN) + @test stats.selected == 4 && isempty(stats.failed) # A failed deterministic start throws. @test_throws "refinement failed" run(2.0) # Errors other than `SolverFailure` propagate from random starts too. - @test_throws DomainError MatrixCovers._multistart_run(builder, - (x, A) -> x[1] == 0.5 && throw(DomainError(x)), objective, ones(1, 1), 4, 2.0, nothing) + @test_throws DomainError MatrixCovers._multistart_run(ESS, builder, + (x, A) -> x[1] == 0.5 ? throw(DomainError(x)) : refine!(x, A, NaN), objective, + ones(1, 1), 4, 2.0, nothing) end @testset "converged cover covariance" begin @@ -365,7 +370,7 @@ end @test covaries_objective(AbsLinear{2}(), sc, A, dr, dc; rtol=1e-12) a, b = sc(A) @test cover_objective(AbsLinear{2}(), a, b, A) < - cover_objective(AbsLinear{2}(), soft_symcover(AbsLinear{2}(), A), A) * (1 - 1e-6) + cover_objective(AbsLinear{2}(), first(soft_symcover(AbsLinear{2}(), A)), A) * (1 - 1e-6) end @testset "soft AbsLog{2} is the exact unconstrained minimum" begin @@ -485,8 +490,8 @@ end @testset "refining the multistart's own start reproduces it: $ϕ" for ϕ in SOFT_PENALTIES a = initialize_symcover(Asym; strategy=:geomean, feasible=:none) - scales(soft_symcover!(ϕ, a, Asym)) - @test cover_objective(ϕ, a, Asym) ≈ cover_objective(ϕ, scales(soft_symcover(ϕ, Asym)), Asym) rtol=1e-6 + soft_symcover!(ϕ, a, Asym) + @test cover_objective(ϕ, a, Asym) ≈ cover_objective(ϕ, first(soft_symcover(ϕ, Asym)), Asym) rtol=1e-6 end @testset "no-ϕ form defaults to PowerMean{2}" begin @@ -509,18 +514,18 @@ end o(AbsLog{2}(), first(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) + @test !isapprox(o(AbsLinear{2}(), first(soft_symcover!(AbsLinear{2}(), s1, Abasins))), + o(AbsLinear{2}(), first(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 SOFT_PENALTIES a = fill(0.01, 3) @test !iscover(a, Asym) - @test scales(soft_symcover!(ϕ, a, Asym)) === a + @test first(soft_symcover!(ϕ, a, Asym)) === a b, c = fill(0.01, 2), fill(0.01, 3) @test !iscover(b, c, Agen) - @test scales(soft_cover!(ϕ, b, c, Agen)) === (b, c) + @test soft_cover!(ϕ, b, c, Agen)[1:2] === (b, c) end @testset "invalid starts throw, naming the function called" begin @@ -539,24 +544,24 @@ end @testset "unsupported rows are zeroed: $ϕ" for ϕ in SOFT_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 scales(soft_symcover!(ϕ, a, Az))[2] == 0 + @test first(soft_symcover!(ϕ, a, Az))[2] == 0 end @testset "asymmetric refiners pin the balance convention: $ϕ" for ϕ in SOFT_PENALTIES a, b = initialize_cover(Agen; strategy=:geomean, feasible=:none) a .*= 7 # move the gauge; the objective cannot see it b ./= 7 - scales(soft_cover!(ϕ, a, b, Agen)) + soft_cover!(ϕ, a, b, Agen) @test isbalanced(a, b, Agen) end @testset "offset axes propagate: $ϕ" for ϕ in SOFT_PENALTIES Ao = OffsetArray(Asym, -1:1, -1:1) ao = initialize_symcover(Ao; strategy=:geomean, feasible=:none) - scales(soft_symcover!(ϕ, ao, Ao)) + soft_symcover!(ϕ, ao, Ao) @test axes(ao, 1) == axes(Ao, 1) aref = initialize_symcover(Asym; strategy=:geomean, feasible=:none) - scales(soft_symcover!(ϕ, aref, Asym)) + soft_symcover!(ϕ, aref, Asym) @test collect(ao) ≈ aref rtol=1e-6 end end @@ -607,7 +612,7 @@ MatrixCovers.require_abs_symmetric(::HookOnlyMatrix, fname) = nothing @test first(soft_symcover(PowerMean{2}(), M)) ≈ first(soft_symcover(PowerMean{2}(), dense)) rtol=1e-12 for ϕ in (AbsLinear{1}(), AbsLinear{2}(), AbsLog{1}(), AbsLog{2}()) @test symcover(ϕ, M) ≈ symcover(ϕ, dense) rtol=1e-10 - @test scales(soft_symcover(ϕ, M)) ≈ scales(soft_symcover(ϕ, dense)) rtol=1e-6 + @test first(soft_symcover(ϕ, M)) ≈ first(soft_symcover(ϕ, dense)) rtol=1e-6 @test iscover(symcover(ϕ, M), dense; rtol=8eps()) @test cover_objective(ϕ, symcover(ϕ, M), M) ≈ cover_objective(ϕ, symcover(ϕ, dense), dense) rtol=1e-10 diff --git a/test/storage_types.jl b/test/storage_types.jl index ff08a26..a5dc7e7 100644 --- a/test/storage_types.jl +++ b/test/storage_types.jl @@ -194,13 +194,13 @@ end [4.0 0.0 1.0; 0.0 0.0 0.0; 1.0 0.0 2.0]] # second has a zero row/column for M in symdenses ad, _ = symcover_min(AbsLog{2}(), M) - asd = soft_symcover(AbsLinear{2}(), M) + asd, _ = soft_symcover(AbsLinear{2}(), M) for A in (sparse(M), Symmetric(sparse(triu(M)), :U), Symmetric(sparse(tril(M)), :L), Hermitian(sparse(triu(M)), :U)) a, _ = symcover_min(AbsLog{2}(), A) @test iscover(a, M; atol=1e-7) @test cover_objective(AbsLog{2}(), a, M) ≈ cover_objective(AbsLog{2}(), ad, M) rtol = 1e-7 atol = 1e-10 - as = soft_symcover(AbsLinear{2}(), A) + as, _ = soft_symcover(AbsLinear{2}(), A) @test cover_objective(AbsLinear{2}(), as, M) ≈ cover_objective(AbsLinear{2}(), asd, M) rtol = 1e-7 atol = 1e-10 # Scale vectors are dense: return a plain Vector, matching cover/symcover. @test as isa Vector{Float64} @@ -210,8 +210,8 @@ end # Every soft_symcover penalty returns a dense Vector on sparse-backed input. let Ssp = sparse(symdenses[1]) for ϕ in SOFT_PENALTIES - @test scales(soft_symcover(ϕ, Ssp)) isa Vector{Float64} - @test scales(soft_symcover(ϕ, Symmetric(Ssp))) isa Vector{Float64} + @test first(soft_symcover(ϕ, Ssp)) isa Vector{Float64} + @test first(soft_symcover(ϕ, Symmetric(Ssp))) isa Vector{Float64} end end diff --git a/test/unitful.jl b/test/unitful.jl index 3fef701..c9f216a 100644 --- a/test/unitful.jl +++ b/test/unitful.jl @@ -174,11 +174,11 @@ end for ϕ in PENALTIES @test unit.(symcover(ϕ, A)) == UA - @test unit.(scales(soft_symcover(ϕ, A))) == UA - @test unit.(scales(symcover_min(ϕ, A))) == UA + @test unit.(first(soft_symcover(ϕ, A))) == UA + @test unit.(first(symcover_min(ϕ, A))) == UA @test all(unit.(v) == UA for v in cover(ϕ, A)) - @test all(unit.(v) == UA for v in scales(soft_cover(ϕ, A))) - @test all(unit.(v) == UA for v in scales(cover_min(ϕ, A))) + @test all(unit.(v) == UA for v in soft_cover(ϕ, A)[1:2]) + @test all(unit.(v) == UA for v in cover_min(ϕ, A)[1:2]) end for ϕ in (PowerMean{1}(), PowerMean{2}()) @test unit.(scales(soft_symcover(ϕ, A))) == UA @@ -247,7 +247,7 @@ @test unit.(a) == unit.(b) == UA a = initialize_symcover(A) - @test soft_symcover!(AbsLinear{2}(), a, A) === a + @test first(soft_symcover!(AbsLinear{2}(), a, A)) === a @test unit.(a) == UA a, b = initialize_cover(A) From 8944d0351a90d23f6a3018260e9cc88ec87e6bbe Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 9 Oct 2026 05:39:30 -0500 Subject: [PATCH 5/7] Count HiGHS iterations through MathOptInterface Reading HiGHS's QP iteration counter required its C interface, whose info keys are not guaranteed across the HiGHS versions allowed by `[compat]`. The counts are now the simplex and interior-point iterations that MathOptInterface reports, which do not include the QP stage. Assisted-by: Claude Opus 5.5 --- ext/MatrixCoversJuMPExt.jl | 16 ++-------------- src/stats.jl | 5 +++-- 2 files changed, 5 insertions(+), 16 deletions(-) diff --git a/ext/MatrixCoversJuMPExt.jl b/ext/MatrixCoversJuMPExt.jl index 66135c9..0aacf34 100644 --- a/ext/MatrixCoversJuMPExt.jl +++ b/ext/MatrixCoversJuMPExt.jl @@ -12,20 +12,8 @@ check_solved(model, fname) = _highs_stats() = ExternalSolverStats(:HiGHS, Symbol[], Float64[], Int[], Float64[]) -# HiGHS reports its simplex, interior-point (with crossover), and QP iterations -# separately, and MOI exposes only the first two. -function _highs_iterations(model) - opt = JuMP.unsafe_backend(model) - n = 0 - for key in ("simplex_iteration_count", "ipm_iteration_count", "crossover_iteration_count", - "qp_iteration_count") - r = Ref{Cint}(0) - HiGHS.Highs_getIntInfoValue(opt, key, r) == HiGHS.kHighsStatusOk || - error("HiGHS did not report $key") - n += r[] - end - return n -end +# Iterations reported through MathOptInterface; HiGHS's QP iterations are not among them. +_highs_iterations(model) = JuMP.simplex_iterations(model) + JuMP.barrier_iterations(model) # Check that the last solve of `model` succeeded and append it to `stats` as a stage. function record_solved!(stats::ExternalSolverStats, model, fname) diff --git a/src/stats.jl b/src/stats.jl index 9aaca06..0f7f420 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -272,8 +272,9 @@ A stage that does not end with a solved status throws a models work in log scales and may include slack variables or constant terms, so this need not equal [`cover_objective`](@ref) of the returned cover. - `niters::Vector{Int}`: the iterations of each stage as the solver reports them - (for HiGHS, its simplex, interior-point, and QP iterations; for Ipopt, its - interior-point iterations). + through MathOptInterface: for HiGHS, its simplex and interior-point iterations + (its quadratic-program iterations are not reported, so the second stage + counts 0); for Ipopt, its interior-point iterations. - `solvetime::Vector{Float64}`: the solver's run time of each stage, in seconds. `niters` and `solvetime` may change meaning between releases. From d84c43d54fc9eb43f85e54b98f6003ac990c39df Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 9 Oct 2026 06:15:32 -0500 Subject: [PATCH 6/7] Document the solver statistics outputs Update the docstrings of the optimizing solvers for the trailing stats output, the `tol` keyword of the AbsLog{2} hard covers, and the policy of warning (native solvers) versus throwing SolverFailure (JuMP-backed solvers). Port the doctests to the new return values and add a short manual section on convergence and solver statistics. Assisted-by: Claude Opus 5.5 --- docs/src/index.md | 48 +++++++++++++++++++++++++++++-- src/minimal_covers.jl | 67 ++++++++++++++++++++++++++++--------------- src/soft_covers.jl | 64 +++++++++++++++++++++++------------------ src/transversal.jl | 16 +++++++---- 4 files changed, 135 insertions(+), 60 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 2216163..040da3f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -247,7 +247,9 @@ julia> using MatrixCovers, JuMP, HiGHS julia> S = [4 1 0; 1 1 5; 0 5 2]; -julia> round.(symcover_min(AbsLog{1}(), S); digits=6) +julia> a, stats = symcover_min(AbsLog{1}(), S); + +julia> round.(a; digits=6) 3-element Vector{Float64}: 2.0 1.0 @@ -311,13 +313,17 @@ julia> using MatrixCovers, JuMP, Ipopt julia> S = [4 1 0; 1 1 5; 0 5 2]; -julia> round.(symcover_min(AbsLinear{2}(), S); digits=6) +julia> a, stats = symcover_min(AbsLinear{2}(), S); + +julia> round.(a; digits=6) 3-element Vector{Float64}: 2.0 1.0 5.0 -julia> round.(symcover_min(AbsLinear{2}(), S; strategies=(:geomean,)); digits=6) +julia> a, stats = symcover_min(AbsLinear{2}(), S; strategies=(:geomean,)); + +julia> round.(a; digits=6) 3-element Vector{Float64}: 2.0 1.0 @@ -337,6 +343,42 @@ julia> round.(a0; digits=6) For the convex `AbsLog` and `PowerMean` penalties, the start does not change the result. +## Convergence and solver statistics + +The optimizing solvers return statistics as their last output, reporting +whether the solver met its stopping criterion (the heuristics, which run a fixed +number of steps, return none): + +```jldoctest stats +julia> using MatrixCovers + +julia> A = [1 2 3; 6 5 4]; + +julia> a, b, stats = cover_min(A); + +julia> MatrixCovers.converged(stats) +true + +julia> MatrixCovers.residual(stats) <= MatrixCovers.tolerance(stats) +true +``` + +`a, b = cover_min(A)` discards the statistics, but `a = symcover_min(A)` binds +the pair; write `a, stats = symcover_min(A)`. + +Each solver family returns its own subtype of +[`AbstractCoverStats`](@ref MatrixCovers.AbstractCoverStats). The accessors +[`converged`](@ref MatrixCovers.converged), +[`iterations`](@ref MatrixCovers.iterations), +[`residual`](@ref MatrixCovers.residual), and +[`tolerance`](@ref MatrixCovers.tolerance) are public but not exported, and each +subtype implements those that apply to its solver. + +A native solver that stops before meeting its criterion warns and returns its +result; a hard cover returned this way still covers `A`. Test +`MatrixCovers.converged(stats)` to reject such results. The JuMP-backed solvers +throw a [`SolverFailure`](@ref MatrixCovers.SolverFailure) instead. + ## Consuming one factor alone: gauges and Gram covers The balanced factors of an asymmetric cover are deterministic but not diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index aa75bb8..2e6dc78 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -6,8 +6,8 @@ # ============================================================ """ - a = symcover_min(ϕ, A; kwargs...) - a = symcover_min(A; kwargs...) + a, stats = symcover_min(ϕ, A; kwargs...) + a, stats = symcover_min(A; kwargs...) Return the ϕ-minimal symmetric hard cover of `A`: the vector `a` minimizing `∑ ϕ(|A[i,j]|/(a[i]*a[j]))` subject to `a[i]*a[j] >= |A[i,j]|`. The @@ -22,11 +22,19 @@ Supported ϕ values: `AbsLog` is convex in the log scales. If `AbsLog{1}` has multiple minima, the method selects the one with the smallest `AbsLog{2}` objective. +`stats` is an [`AbstractCoverStats`](@ref MatrixCovers.AbstractCoverStats); +[`MatrixCovers.converged(stats)`](@ref MatrixCovers.converged) reports whether +`a` meets the solver's stopping criterion. The native solver warns when it does +not, and returns a cover of `A` that may not minimize the objective. JuMP-backed +solves throw a [`SolverFailure`](@ref MatrixCovers.SolverFailure) instead. + # Extended help The native solver accepts `κ` (initial augmented-Lagrangian penalty, default `1e2`), `maxouter` (multiplier updates, default `32`; `0` returns the -unconstrained fit), `maxiter` (Newton steps per update), `fillbudget` and +unconstrained fit, shifted to cover `A`, with `converged(stats) == false` and no +warning), `maxiter` (Newton steps per update), `tol` (the tolerance +on the KKT residual, default `1000*eps` of the working type), `fillbudget` and `flopbudget` (see below), and `linsolve`: - `:dense` factorizes dense normal equations at O(n³) per Newton step. @@ -56,9 +64,12 @@ exceed them are solved by conjugate gradients instead. If the solver warns that the result may not minimize the objective, follow the advice in the warning: increase `maxouter` when the update limit was reached, -and `κ` when the residual stopped contracting before the limit. +and `κ` when the residual stopped contracting before the limit. A caller who +accepts a larger residual can pass a larger `tol`, or compare +[`MatrixCovers.residual(stats)`](@ref MatrixCovers.residual) with its own +threshold. -Without a warning, the result typically matches the exact minimizer to a +When `converged(stats)` is `true`, the result typically matches the exact minimizer to a relative accuracy of about `1e-8`, and occasionally only `1e-7`: the stopping test measures feasibility and complementary slackness but not stationarity. Scale covariance holds to the same accuracy. For a diagonal rescaling @@ -74,8 +85,8 @@ function symcover_min end symcover_min(A::AbstractMatrix; kwargs...) = symcover_min(AbsLog{2}(), A; kwargs...) """ - a, b = cover_min(ϕ, A) - a, b = cover_min(A) + a, b, stats = cover_min(ϕ, A; kwargs...) + a, b, stats = cover_min(A; kwargs...) Return the ϕ-minimal asymmetric hard cover of `A`: vectors `a`, `b` minimizing `∑ ϕ(|A[i,j]|/(a[i]*b[j]))` subject to `a[i]*b[j] >= |A[i,j]|`. The @@ -94,7 +105,8 @@ distant rows or columns, such as `norm(a)`, inherit the spread. the floating-point range throw an error; converting `A` to `BigFloat` extends the range. -Supported ϕ values: +Supported ϕ values, `stats`, and the failure behavior are as for +[`symcover_min`](@ref): - `AbsLog{2}()`: native. - `AbsLog{1}()`: requires JuMP and HiGHS. - `AbsLinear{1}()` and `AbsLinear{2}()`: require JuMP and Ipopt and return the @@ -105,8 +117,8 @@ method selects the one with the smallest `AbsLog{2}` objective. # Extended help -The native solver accepts the same `κ`, `maxouter`, `maxiter`, `fillbudget`, -`flopbudget`, and `linsolve` keywords as [`symcover_min`](@ref). For +The native solver accepts the same `κ`, `maxouter`, `maxiter`, `tol`, +`fillbudget`, `flopbudget`, and `linsolve` keywords as [`symcover_min`](@ref). For `:woodbury`, an `m × n` matrix may omit at most `min(m,n) ÷ 4` entries per row or column and `4 * max(m,n)` entries in total. `:dense` costs O((m+n)³) per Newton step; sparse matrices default to `:lsqr`. @@ -120,7 +132,8 @@ continues from its result. The entries that are tight at the solution form a graph, and each step solves the objective with a spanning forest of them held tight, a sparse graph Laplacian with one unknown per tree. The method ends when the KKT conditions hold to rounding error, which makes the result the exact -minimizer; only when it fails does the solver warn. +minimizer; `stats.polish` records it. When it fails, the solver returns the +multiplier iterate and warns. See also: [`symcover_min`](@ref), [`cover`](@ref), [`cover_min!`](@ref). """ @@ -128,11 +141,12 @@ function cover_min end cover_min(A::AbstractMatrix; kwargs...) = cover_min(AbsLog{2}(), A; kwargs...) """ - a = symcover_min!(ϕ, a, A; kwargs...) - a = symcover_min!(a, A; kwargs...) + a, stats = symcover_min!(ϕ, a, A; kwargs...) + a, stats = symcover_min!(a, A; kwargs...) Refine the symmetric hard cover `a` in place. The no-ϕ form uses `AbsLog{2}()`; -supported penalties and keywords match [`symcover_min`](@ref). +supported penalties, keywords, `stats`, and failure behavior match +[`symcover_min`](@ref). `a` must cover `A` and be positive on supported rows. Unsupported scales are ignored on input and set to zero. Use [`initialize_symcover`](@ref) or @@ -149,13 +163,13 @@ symcover_min!(a::AbstractVector, A::AbstractMatrix; kwargs...) = symcover_min!(AbsLog{2}(), a, A; kwargs...) """ - a, b = cover_min!(ϕ, a, b, A; kwargs...) - a, b = cover_min!(a, b, A; kwargs...) + a, b, stats = cover_min!(ϕ, a, b, A; kwargs...) + a, b, stats = cover_min!(a, b, A; kwargs...) Refine the asymmetric hard cover `(a, b)` in place. The start must cover `A` and be positive on supported rows and columns; unsupported scales are zeroed. The -no-ϕ form uses `AbsLog{2}()`; supported penalties and keywords match -[`cover_min`](@ref). +no-ϕ form uses `AbsLog{2}()`; supported penalties, keywords, `stats`, and +failure behavior match [`cover_min`](@ref). Equivalent starts `(c*a, b/c)` give the same balanced result. @@ -2439,11 +2453,18 @@ function cover_transversal_jump end """ SolverFailure(msg) -Thrown when an external solver (Ipopt or HiGHS) stops without reaching a usable -local minimum: it reports a status other than solved, or it returns scales that -diverge on a plateau of the objective. A different start may succeed. The -`AbsLinear` multistarts of [`soft_cover`](@ref) and [`soft_symcover`](@ref) -leave out perturbed starts that throw this. +Thrown by the `AbsLog{1}` and `AbsLinear` forms of [`symcover_min`](@ref), +[`cover_min`](@ref), [`soft_symcover`](@ref), [`soft_cover`](@ref), and their +mutating forms when the external solver (HiGHS or Ipopt) stops without reaching +a usable local minimum: it reports a status other than solved, or it returns +scales that diverge on a plateau of the objective. A different start may +succeed. The `AbsLinear` multistarts of `soft_cover` and `soft_symcover` leave +out perturbed starts that throw this and record them in +[`MultistartStats`](@ref MatrixCovers.MultistartStats). + +The native solvers do not throw this: one that stops before meeting its +stopping criterion warns and returns a result whose statistics report +`converged(stats) == false`. """ struct SolverFailure <: Exception msg::String diff --git a/src/soft_covers.jl b/src/soft_covers.jl index 168c9ae..69fb3b2 100644 --- a/src/soft_covers.jl +++ b/src/soft_covers.jl @@ -5,37 +5,41 @@ # ============================================================ """ - a = soft_symcover(ϕ, A; kwargs...) - a = soft_symcover(A; kwargs...) + a, stats = soft_symcover(ϕ, A; kwargs...) + a, stats = soft_symcover(A; kwargs...) Minimize `∑ ϕ(|A[i,j]|/(a[i]*a[j]))` for symmetric `A`, without a hard coverage constraint. The default penalty is `PowerMean{2}()`. Supported penalties and their keywords are: -- `PowerMean{p}()` (default): the convex minimum, computed by damped - simultaneous power-mean updates `a[k] ← sqrt(a[k] * M_p(|A[k,j]|/a[j]))`, where - `M_p` is the `p`-power mean over the nonzeros of row `k`, followed when they - are slow by Newton steps with backtracking on the objective. Every update +- `PowerMean{p}()` (default): the convex minimum, computed by damped simultaneous power-mean updates + `a[k] ← sqrt(a[k] * M_p(|A[k,j]|/a[j]))`, where `M_p` is the `p`-power mean + over the nonzeros of row `k`, followed when they are slow by Newton steps + with backtracking on the objective. Every update decreases the objective. The iteration stops when, in every nonzero row, the mean of `r^p` over the ratios `r = |A[k,j]|/(a[k]*a[j])` of the nonzero entries is within `tol` of one. It computes in at least `Float64`, and `tol` defaults to `4096*eps` of that type, and a warning reports an iteration that ends without reaching it. The keywords `maxiter`, `newton`, `maxnewton`, and `linsolve` are described in the extended help. -- `AbsLog{2}()`: the convex minimum, computed by one linear solve; `linsolve` has - the same meaning as in [`symcover_min`](@ref). +- `AbsLog{2}()`: the convex minimum, computed by one linear solve; `linsolve` + has the same meaning as in [`symcover_min`](@ref). - `AbsLog{1}()`: the convex minimum, computed as a linear program; requires JuMP and HiGHS. The minimizer need not be unique, and among the minimizers the one with the least `AbsLog{2}` objective is returned. -- `AbsLinear{1}()`, `AbsLinear{2}()`: the objective is not convex. Several starts - are each refined to a local minimum, and the one with the least objective is - returned; requires JuMP and Ipopt. `starts` (default 5) is the number of - starts: a few deterministic starts, then log-normal perturbations of the - first deterministic start with spread `σ` (default 2.0; `sigma` is an - alias), drawn from `rng` (default `MersenneTwister(0)`). A perturbed start - whose refinement throws a [`SolverFailure`](@ref) is left out of the - selection; a failing deterministic start throws. +- `AbsLinear{1}()`, `AbsLinear{2}()`: the objective is not convex. Several starts are each refined to a local minimum, and the one with + the least objective is returned; requires JuMP and Ipopt. `starts` + (default 5) is the number of starts: a few deterministic starts, then + log-normal perturbations of the first deterministic start with spread `σ` + (default 2.0; `sigma` is an alias), drawn from `rng` (default + `MersenneTwister(0)`). A perturbed start whose refinement throws a + [`SolverFailure`](@ref) is left out of the selection; a failing + deterministic start throws. + +`stats` is an [`AbstractCoverStats`](@ref MatrixCovers.AbstractCoverStats). +The JuMP-backed solves (`AbsLog{1}`, `AbsLinear`) throw a `SolverFailure` when +the external solver fails. Unsupported rows receive zero scale. When a connected component of the support is bipartite with no diagonal entry, the products on the support do not @@ -49,7 +53,10 @@ See also: [`symcover`](@ref), [`cover_objective`](@ref), [`soft_symcover!`](@ref ```jldoctest julia> A = [4 -1; -1 1]; -julia> a = soft_symcover(A); +julia> a, stats = soft_symcover(A); + +julia> MatrixCovers.converged(stats) +true julia> R = abs.(A) ./ (a .* a'); @@ -58,7 +65,7 @@ julia> round.(sum(R .^ 2; dims=2); digits=8) # each row's mean square ratio is 2.0 2.0 -julia> round.(soft_symcover([0 1; 1 0]); digits=4) +julia> round.(first(soft_symcover([0 1; 1 0])); digits=4) 2-element Vector{Float64}: 1.0 1.0 @@ -107,8 +114,8 @@ function soft_symcover(ϕ::AbsLinear, A::AbstractMatrix; starts::Int=5, end """ - a = soft_symcover!(ϕ, a, A; kwargs...) - a = soft_symcover!(a, A; kwargs...) + a, stats = soft_symcover!(ϕ, a, A; kwargs...) + a, stats = soft_symcover!(a, A; kwargs...) Refine one symmetric soft-cover start in place. The no-ϕ form uses `PowerMean{2}()`. Build a start with [`initialize_symcover`](@ref) and @@ -136,12 +143,13 @@ function soft_symcover!(::AbsLog{2}, a::AbstractVector, A::AbstractMatrix; kwarg end """ - a, b = soft_cover(ϕ, A; kwargs...) - a, b = soft_cover(A; kwargs...) + a, b, stats = soft_cover(ϕ, A; kwargs...) + a, b, stats = soft_cover(A; kwargs...) Minimize `∑ ϕ(|A[i,j]|/(a[i]*b[j]))` without a hard coverage constraint. This is the asymmetric form of [`soft_symcover`](@ref). The default penalty is -`PowerMean{2}()`. +`PowerMean{2}()`. `stats` and the failure behavior are as for +`soft_symcover`. Supported penalties and their keywords are: @@ -161,8 +169,8 @@ Supported penalties and their keywords are: `abs.(A)` is exactly symmetric, the minimizer has `b == a` under the balance convention, and `soft_cover` instead computes `a` by the algorithm of [`soft_symcover`](@ref) (with the same keywords; `maxiter` then counts its - updates) and returns `(a, copy(a))`. The in-place [`soft_cover!`](@ref) always - uses the alternating iteration. + updates) and returns `(a, copy(a), stats)`. The in-place + [`soft_cover!`](@ref) always uses the alternating iteration. - `AbsLog{2}()`: the convex minimum, computed by one linear solve; `linsolve` has the same meaning as in [`cover_min`](@ref). - `AbsLog{1}()`: the convex minimum, computed as a linear program; requires JuMP @@ -188,7 +196,7 @@ See also: [`cover`](@ref), [`soft_symcover`](@ref), [`soft_cover!`](@ref), [`cov ```jldoctest; filter = r"(\\d+\\.\\d{4})\\d+" => s"\\1" julia> A = [1 2 3; 6 5 4]; -julia> a, b = soft_cover(A); +julia> a, b, stats = soft_cover(A); julia> a * b' 2×3 Matrix{Float64}: @@ -212,8 +220,8 @@ function soft_cover(ϕ::AbsLinear, A::AbstractMatrix; starts::Int=4, end """ - a, b = soft_cover!(ϕ, a, b, A; kwargs...) - a, b = soft_cover!(a, b, A; kwargs...) + a, b, stats = soft_cover!(ϕ, a, b, A; kwargs...) + a, b, stats = soft_cover!(a, b, A; kwargs...) Refine the starting point `(a, b)` into a soft cover of `A` in place. Scales must be finite and positive on supported rows and columns; unsupported scales are diff --git a/src/transversal.jl b/src/transversal.jl index e6a058d..ed2465b 100644 --- a/src/transversal.jl +++ b/src/transversal.jl @@ -1,7 +1,7 @@ # Transversal-tight minimal covers and the maximum-product transversal they rest on. """ - a, b = cover_transversal(A; kwargs...) + a, b, stats = cover_transversal(A; kwargs...) Return the transversal-tight minimal cover of the square matrix `A`: the hard cover minimizing `∑ log(a[i]*b[j]/|A[i,j]|)^2` over the support of `A`, among @@ -23,6 +23,10 @@ balance convention of [`cover_min`](@ref). `A` must be square and structurally nonsingular, i.e., have at least one transversal of nonzero entries; otherwise an `ArgumentError` is thrown. +`stats` is an [`AugmentedLagrangianStats`](@ref MatrixCovers.AugmentedLagrangianStats); +as for [`cover_min`](@ref), the solver warns when +[`MatrixCovers.converged(stats)`](@ref MatrixCovers.converged) is `false`. + # Extended help A maximum-product transversal `σ` is found by shortest augmenting paths on the @@ -32,7 +36,7 @@ row scales alone, with one constraint per entry off the transversal. The augmented-Lagrangian solver of [`cover_min`](@ref) solves it, starting from the dual variables of the transversal, which are a hard cover tight on `σ`. -The keywords `κ`, `maxouter`, `maxiter`, `fillbudget`, `flopbudget`, and +The keywords `κ`, `maxouter`, `maxiter`, `tol`, `fillbudget`, `flopbudget`, and `linsolve` are as for [`cover_min`](@ref). `linsolve=:woodbury` solves the row-scale system as a sparse matrix (conjugate gradients, or sparse Cholesky when poorly conditioned). For `n×n` `A`, it requires at most `n ÷ 4` zeros in @@ -55,12 +59,12 @@ function cover_transversal(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::Int end """ - a, b = cover_transversal!(a, b, A; kwargs...) + a, b, stats = cover_transversal!(a, b, A; kwargs...) Mutating counterpart of [`cover_transversal`](@ref): writes the cover into `a` -and `b` and returns them. The initial contents of `a` and `b` are ignored. -`eachindex(a)` must match `axes(A, 1)` and `eachindex(b)` must match -`axes(A, 2)`. +and `b` and returns them with the solver's statistics. The initial contents of +`a` and `b` are ignored. `eachindex(a)` must match `axes(A, 1)` and +`eachindex(b)` must match `axes(A, 2)`. """ function cover_transversal!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) axes(A, 1) == eachindex(a) || throw(DimensionMismatch("indices of `a` must match row-indexing of `A`, got eachindex(a)=$(string(eachindex(a))), axes(A, 1)=$(string(axes(A, 1)))")) From 133689e913e085df23200044bb13c73e606519f7 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 9 Oct 2026 07:08:18 -0500 Subject: [PATCH 7/7] Fix Julia 1.10 inference and rounding-dependent tests `symcover_min(::AbsLinear, A)` built its starts in a comprehension that captured the element type, which Julia 1.10 cannot infer; it now copies a typed template vector. Two tests assumed that a solve without the active-set finish stalls above the tolerance, which depends on rounding and fails on Julia 1.10. They now check that the solve warns exactly when it does not converge. The multistart test's closure no longer captures a type through a local variable. Assisted-by: Claude Opus 5.5 --- src/minimal_covers.jl | 5 +++-- test/helpers.jl | 8 ++++++++ test/minimal_covers.jl | 24 +++++++++++++++--------- test/soft_covers.jl | 3 ++- 4 files changed, 28 insertions(+), 12 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 2e6dc78..555e02c 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -234,13 +234,14 @@ function symcover_min(ϕ::AbsLinear, A::AbstractMatrix; strategies=SYMCOVER_MIN_ isempty(strategies) && throw(ArgumentError("symcover_min: `strategies` must name at least one starting cover")) T = float(real(eltype(A))) - starts = [similar(Array{T}, ax) for _ in strategies] + a0 = similar(Array{T}, ax) + starts = [similar(a0) for _ in strategies] # Skip strategies that cannot produce a start for `A`. built = [_initialize_symcover!(a, A, strategy, :inflate) for (a, strategy) in zip(starts, strategies)] any(built) || throw(ArgumentError("symcover_min: no strategy in $(string(strategies)) yields a starting cover of `A`")) labels = Symbol[strategy for (strategy, ok) in zip(strategies, built) if ok] - covers = [a for (a, ok) in zip(starts, built) if ok] + covers = starts[built] stats = [last(symcover_min!(ϕ, a, A))::ExternalSolverStats for a in covers] objs = [cover_objective(ϕ, a, A) for a in covers] k = _multistart_select(objs) diff --git a/test/helpers.jl b/test/helpers.jl index f854484..bcfe90b 100644 --- a/test/helpers.jl +++ b/test/helpers.jl @@ -13,6 +13,14 @@ scales(r::Tuple{AbstractVector,MatrixCovers.AbstractCoverStats}) = r[1] scales(r::Tuple{AbstractVector,AbstractVector}) = r scales(r::Tuple{AbstractVector,AbstractVector,MatrixCovers.AbstractCoverStats}) = (r[1], r[2]) +# `f()` and the log records it emits, for checks where whether a warning is +# expected depends on rounding. +function capture_logs(f) + logger = Test.TestLogger() + value = Base.CoreLogging.with_logger(f, logger) + return value, logger.logs +end + """ isbalanced(a, b, A) diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index c19b54c..976ebe2 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -931,13 +931,16 @@ end # With the penalty held at the cap they are zeroed directly; left to drain by # `2(κ-1)z` per update they hold the KKT residual at 2.6e-6. # Without the active-set finish (`polish=false`), slack multipliers are zeroed. - # With the diagonal preconditioner the residual then stalls above the - # tolerance, which warns; the default solve's finish certifies the result. + # Whether the residual then reaches the tolerance depends on rounding (with the + # diagonal preconditioner it can stall just above it); the solve warns exactly + # when it does not. The default solve's finish certifies the result. A = banded_sparse_sym(StableRNG(5), 10_000, 3, 1.0) for fillbudget in (MatrixCovers.LSQR_FILL_BUDGET, 0) - solve = () -> MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, fillbudget, polish=false) - a, s = fillbudget == 0 ? (@test_logs (:warn, r"try a larger `κ`") solve()) : (@test_nowarn solve()) - @test s.converged == (fillbudget != 0) + (a, s), logs = capture_logs() do + MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, fillbudget, polish=false) + end + @test isempty(logs) == s.converged + @test all(l -> l.level == Base.CoreLogging.Warn && occursin("try a larger `κ`", l.message), logs) @test s.precond === (fillbudget == 0 ? :diagonal : :factor) @test s.nzeroed > 0 @test maximum(s.κs) == 1e5 @@ -987,11 +990,14 @@ end @test cover_objective(AbsLog{2}(), a1, b1, A) ≈ cover_objective(AbsLog{2}(), a0, b0, A) rtol=1e-10 @test a1 .* b1' ≈ a0 .* b0' rtol=1e-7 @test iscover(a1, b1, A) - # Without the polish the iteration keeps zeroing, and stalls above the - # tolerance. - _, _, s2 = @test_logs (:warn, r"try a larger `κ`") MatrixCovers._cover_min_abslog2(A; polish=false) + # Without the polish the iteration keeps zeroing. Whether it then reaches the + # tolerance depends on rounding; it warns exactly when it does not. + (_, _, s2), logs = capture_logs() do + MatrixCovers._cover_min_abslog2(A; polish=false) + end @test s2.nzeroed > 0 - @test !s2.converged + @test isempty(logs) == s2.converged + @test all(l -> l.level == Base.CoreLogging.Warn && occursin("try a larger `κ`", l.message), logs) end @testset "MMC multiplier update is exact on an analytic problem" begin diff --git a/test/soft_covers.jl b/test/soft_covers.jl index 1c5170b..26a60c6 100644 --- a/test/soft_covers.jl +++ b/test/soft_covers.jl @@ -331,7 +331,8 @@ end x[1] == bad && throw(MatrixCovers.SolverFailure("refinement failed")) return ESS(:Ipopt, [:LOCALLY_SOLVED], [x[1]], [1], [0.0]) end - run(bad) = MatrixCovers._multistart_run(ESS, builder, (x, A) -> refine!(x, A, bad), objective, + run(bad) = MatrixCovers._multistart_run(MatrixCovers.ExternalSolverStats, builder, + (x, A) -> refine!(x, A, bad), objective, ones(1, 1), 4, 2.0, nothing) # A failed random start is excluded from the selection, even when its start is best. x, stats = @inferred run(0.5)