Report solver stats

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

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

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