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/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..0aacf34 100644
--- a/ext/MatrixCoversJuMPExt.jl
+++ b/ext/MatrixCoversJuMPExt.jl
@@ -3,13 +3,28 @@ 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[])
+
+# 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)
+    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 +58,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 +83,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 +140,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 +157,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 +237,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 +289,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 +301,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/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..ca96764 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, PowerMeanStats, ExternalSolverStats, MultistartStats, 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/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..555e02c 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.
 
@@ -172,33 +186,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.
@@ -208,23 +234,30 @@ 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)]
-    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 = 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)
+    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
 
 # ============================================================
@@ -1170,18 +1203,19 @@ end
 # Augmented-Lagrangian defaults.
 const AL_PENALTY = 1e2
 const AL_MAXOUTER = 32
+const AL_MAXITER = 40
 
-# 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. With
+# `maxouter = 0` there is no iteration and no warning.
 function _warn_unconverged(fname::Symbol, stats, maxouter::Int)
-    stats.converged && 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 ?
+    (stats.converged || maxouter == 0) && 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
 
@@ -1343,11 +1377,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] <= 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`
+# (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
@@ -1359,8 +1394,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
@@ -1822,8 +1860,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))
@@ -1835,7 +1872,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
@@ -1927,7 +1964,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)
@@ -1938,12 +1975,17 @@ 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,
+                               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.
@@ -1956,7 +1998,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)
@@ -2043,11 +2085,11 @@ 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 = _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,13 +2110,15 @@ 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,
+                            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) ||
@@ -2092,7 +2136,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)
@@ -2302,11 +2346,11 @@ 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 = _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
@@ -2410,11 +2454,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
@@ -2424,7 +2475,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/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/soft_covers.jl b/src/soft_covers.jl
index d5793a5..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
@@ -88,11 +95,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.
@@ -106,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
@@ -129,17 +137,19 @@ 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
 
 """
-    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:
 
@@ -159,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
@@ -186,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}:
@@ -196,9 +206,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.
@@ -209,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
@@ -231,10 +242,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.
@@ -321,27 +332,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
@@ -350,22 +360,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
@@ -374,18 +384,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
@@ -393,10 +402,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)
@@ -408,16 +417,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/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..0f7f420
--- /dev/null
+++ b/src/stats.jl
@@ -0,0 +1,344 @@
+# 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)
+
+"""
+    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)
+
+"""
+    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
+  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.
+"""
+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/src/transversal.jl b/src/transversal.jl
index a8a830b..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
@@ -45,26 +49,30 @@ 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
 
 """
-    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)))"))
     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/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/element_types.jl b/test/element_types.jl
index 6b53c6a..50a256a 100644
--- a/test/element_types.jl
+++ b/test/element_types.jl
@@ -4,14 +4,14 @@
 
     @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))
+        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)
         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), first(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, _ = 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..4420c71 100644
--- a/test/extensions.jl
+++ b/test/extensions.jl
@@ -3,8 +3,8 @@
 @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_qmin  = symcover_min(AbsLog{2}(), 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)
         # qmin achieves lower or equal AbsLog{2} objective than symcover and lmin
@@ -14,11 +14,11 @@
     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
-    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,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 = soft_symcover(ϕ, A)
-            a_hard = 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}(), 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{2}(), first(soft_symcover(AbsLog{2}(), A_rank1)), A_rank1) < 1e-8
+    @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 = 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 = 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(ϕ, soft_symcover!(ϕ, a, A), A) :
-                  cover_objective(ϕ, 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 = 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 = 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 = 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 = symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin)
-    a_geo = 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 = symcover_min!(ϕ, initialize_symcover(Ao), Ao)
+    ao, _ = 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) ≈ 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 soft_symcover(ϕ, A) == 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 = soft_symcover(ϕ, A)
+        a, _ = soft_symcover(ϕ, A)
         for strategy in (:hardcover, :geomean, :leaveout)
-            single = 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 -> 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)
@@ -242,29 +242,29 @@ end
 
     for ϕ in (AbsLinear{1}(), AbsLinear{2}(), AbsLog{1}(), AbsLog{2}())
         a, b = soft_cover(ϕ, A)
-        @test soft_cover(ϕ, A) == (a, b)                          # deterministic
+        @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 = soft_cover!(ϕ, 2 .* copy(a0), copy(b0) ./ 2, A)
-        ah, bh = 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 = soft_cover(ϕ, A)
+        a, b, _ = soft_cover(ϕ, A)
         for strategy in (:hardcover, :covariant)
-            sa, sb = 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 = soft_cover(ϕ, As)
+        aa, ab, _ = soft_cover(ϕ, As)
         @test cover_objective(ϕ, aa, ab, As) <=
-              cover_objective(ϕ, 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 = 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 (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 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 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 = 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 = 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 = 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 = 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 symcover_min(ϕ, Abasin) == a   # deterministic
+        @test first(symcover_min(ϕ, Abasin)) == a   # deterministic
 
-        ab, bb = cover_min(ϕ, Aasym)
+        ab, bb, _ = 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, _ = 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 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 = cover_min(ϕ, Abasin)
+        asym_a, asym_b, _ = cover_min(ϕ, Abasin)
         @test cover_objective(ϕ, asym_a, asym_b, Abasin) <=
-              cover_objective(ϕ, 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 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 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 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,))
@@ -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,11 +475,11 @@ 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
-    @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),
@@ -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, _ = symcover_min!(ϕ, symcover(ϕ, M), M)
+        sd, _ = 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, _ = 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 soft_symcover!(ϕ, soft_symcover(ϕ, M), M) ≈
-              soft_symcover!(ϕ, 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 = soft_cover!(ϕ, soft_cover(ϕ, M)..., M)
-        sad, sbd = soft_cover!(ϕ, 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/helpers.jl b/test/helpers.jl
index a7c3e50..bcfe90b 100644
--- a/test/helpers.jl
+++ b/test/helpers.jl
@@ -6,6 +6,21 @@ 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])
+
+# `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/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..45aed66 100644
--- a/test/invariants.jl
+++ b/test/invariants.jl
@@ -8,30 +8,30 @@
 
 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{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{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 = "symcover_min(AbsLog{2})",         f = A -> first(symcover_min(AbsLog{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 -> 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),
 )
 
 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{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{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 = "cover_min(AbsLog{2})",           f = A -> scales(cover_min(AbsLog{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 -> 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),
 )
 
 @testset "cross-notion invariants" begin
@@ -164,11 +164,11 @@ 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.
-            @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
@@ -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/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/minimal_covers.jl b/test/minimal_covers.jl
index b8b6c4e..976ebe2 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,16 +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)
+    # 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 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
@@ -111,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).
@@ -137,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
@@ -149,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) <=
@@ -169,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)
@@ -261,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)
@@ -331,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
@@ -436,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
@@ -601,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.
@@ -632,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
 
@@ -683,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
@@ -731,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
@@ -743,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.
@@ -770,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)
@@ -784,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
@@ -799,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
@@ -812,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
@@ -907,13 +931,22 @@ 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.
+    # 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)
-        a, s = @test_nowarn MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, fillbudget, polish=false)
+        (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
         @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
@@ -957,9 +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.
-    _, _, s2 = 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 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
@@ -984,7 +1022,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)
@@ -995,9 +1033,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
@@ -1009,7 +1047,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..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
@@ -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
@@ -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 5c30138..26a60c6 100644
--- a/test/soft_covers.jl
+++ b/test/soft_covers.jl
@@ -10,14 +10,14 @@ 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))
             @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 = 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 -> 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 = 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  = soft_symcover(ϕ, A_zero)
-        a_small = 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,15 +68,15 @@ 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, _ = soft_symcover(ϕ, B)
+            as, _ = soft_symcover(ϕ, B .* d .* d')
             @test as .* as' ≈ (d .* d') .* (a0 .* a0') rtol=1e-6
         end
     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],
@@ -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,16 +175,16 @@ 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}(), 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
 
     # 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,70 @@ 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(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, 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 +371,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
@@ -406,7 +412,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 +429,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 +469,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
@@ -486,15 +492,15 @@ 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
+        @test cover_objective(ϕ, a, Asym) ≈ cover_objective(ϕ, first(soft_symcover(ϕ, Asym)), Asym) rtol=1e-6
     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.
@@ -505,22 +511,22 @@ 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)),
-                        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 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 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,7 +545,7 @@ 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 first(soft_symcover!(ϕ, a, Az))[2] == 0
     end
 
     @testset "asymmetric refiners pin the balance convention: $ϕ" for ϕ in SOFT_PENALTIES
@@ -604,10 +610,10 @@ 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 soft_symcover(ϕ, M) ≈ 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 d309e04..a5dc7e7 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,25 +193,25 @@ 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)
-        asd = soft_symcover(AbsLinear{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)
+            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}
         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 first(soft_symcover(ϕ, Ssp)) isa Vector{Float64}
+            @test first(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..c9f216a 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.(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 soft_cover(ϕ, A))
-            @test all(unit.(v) == UA for v in 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.(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,53 +219,53 @@
 
         # 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)
-        @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)
-        @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)
-        @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)
-        @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.
         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"]