Add a PowerMean off-transversal, target t<1 minimizer

Project.toml

@@ -1,6 +1,6 @@
 name = "MatrixCovers"
 uuid = "727e6139-ff52-4636-a344-ed1d23e73ffc"
-version = "1.1.0"
+version = "1.2.0"
 authors = ["Tim Holy <tim.holy@gmail.com> and contributors"]
 
 [deps]

docs/src/index.md

@@ -276,7 +276,10 @@ minimize a chosen penalty, at greater cost (see
 for conditioning: on the matrix of
 [Conditioning and the error in `A \ y`](@ref), [`cover_min`](@ref) leaves
 `cond(Ahat)` near `3e8`, because it is not tight on the transversal.
-[`cover_transversal`](@ref) minimizes the penalty among the covers that are.
+[`cover_transversal`](@ref) minimizes the penalty among the covers that are;
+`cover_transversal(PowerMean{2}(), A; target)` instead pushes the ratios of the
+entries off the transversal toward `target`, which separates the blocks of a
+partly decomposable `A`.
 
 For applications that standardize data rather than precondition a solver,
 consider [`soft_cover`](@ref) and [`soft_symcover`](@ref). A soft cover need not

docs/src/objectives.md

@@ -58,15 +58,22 @@ To override a solver's default penalty, pass a penalty as its first argument.
 | [`cover`](@ref) | no | hard (`r ≤ 1`) | heuristic | — |
 | [`symcover_min`](@ref) | yes | hard (`r ≤ 1`) | `AbsLog{2}` (or `AbsLog{1}`, `AbsLinear`) | native for `AbsLog{2}`; else JuMP |
 | [`cover_min`](@ref) | no | hard (`r ≤ 1`) | `AbsLog{2}` (or `AbsLog{1}`, `AbsLinear`) | native for `AbsLog{2}`; else JuMP |
-| [`cover_transversal`](@ref) | no | hard, tight on a maximum-product transversal | `AbsLog{2}` | — |
+| [`cover_transversal`](@ref) | no | hard, tight on a maximum-product transversal | `AbsLog{2}` (or `PowerMean{p}` with a target ratio, optionally max-balanced within blocks) | — |
 | [`soft_symcover`](@ref) | yes | soft (penalized) | `PowerMean{2}` (or `PowerMean{p}`, `AbsLog`, `AbsLinear`) | native for `PowerMean`, `AbsLog{2}`; else JuMP |
 | [`soft_cover`](@ref) | no | soft (penalized) | `PowerMean{2}` (or `PowerMean{p}`, `AbsLog`, `AbsLinear`) | native for `PowerMean`, `AbsLog{2}`; else JuMP |
 
 For hard covers, [`symcover`](@ref) and [`cover`](@ref) are fast heuristics;
 their `_min` counterparts minimize the selected objective.
 [`cover_transversal`](@ref) minimizes the `AbsLog{2}` objective among the hard
 covers that are tight on a maximum-product transversal (see
-[Transversals and conditioning](@ref)). A soft cover has no constraint to
+[Transversals and conditioning](@ref)); with a `PowerMean{p}` penalty and a
+`target` ratio `t`, it instead pushes the ratios of the entries off the
+transversal toward `t`, minimizing `∑ t^p φ(r/t)` over them with `φ` the
+`PowerMean{p}` penalty. Adding `maxbalance=true` replaces that objective within
+each fully indecomposable block by its `p → ∞` limit, max-balancing: across
+every partition of a block's rows, the largest ratio in one direction equals
+the largest in the other, which also minimizes the block's largest ratio. The
+`PowerMean{p}` objective then fixes only the scales between blocks. A soft cover has no constraint to
 satisfy, so there is no separate heuristic: [`soft_symcover`](@ref) and
 [`soft_cover`](@ref) always minimize the objective.
 
@@ -78,6 +85,11 @@ The method used to minimize the objective depends on the penalty:
   converge slowly.
 - `AbsLog{2}` (native) is a single linear least-squares solve for soft covers;
   for hard covers, see [Objective-minimal covers](@ref).
+- `PowerMean{p}` with a target ratio, for [`cover_transversal`](@ref) (native),
+  uses an active-set Newton method that keeps every iterate a hard cover
+  tight on the transversal. With `maxbalance=true`, max-plus policy iteration
+  max-balances each block first, and the Newton method acts on one scale per
+  block.
 - `AbsLog{1}` is a linear program solved with JuMP and HiGHS.
 - `AbsLinear` is not convex; several starts are each refined to a local minimum
   with JuMP and Ipopt, and the best is returned.

docs/src/solvers.md

@@ -26,7 +26,8 @@ and selection among several starts:
 
 In general, plain forms choose their own starts, and `!` forms refine the
 supplied start. The exceptions are [`symcover!`](@ref), [`cover!`](@ref), and
-[`cover_transversal!`](@ref), which write their result in place.
+[`cover_transversal!`](@ref) (with or without a penalty), which write their
+result in place.
 
 The three calls below (the default multistart solver, the solver restricted to
 the `:geomean` start, and an explicit initialization followed by refinement)

ext/MatrixCoversIpoptExt.jl

@@ -152,6 +152,58 @@ function MatrixCovers.symcover_min!(::AbsLinear{1}, a::AbstractVector, A)
     return a, stats
 end
 
+# ============================================================
+# Transversal-tight reference: cover_transversal_jump(::PowerMean{p}, A; target)
+# Row and column log scales with every entry of a maximum-product transversal
+# tight; the objective sums `f(d) = (exp(-p d) - 1)/p + target^p d` over the
+# log slacks `d = α[i] + β[j] - log|A[i,j]| >= 0` of the other entries. Returns
+# `((a, b), logπ)` with `logπ` the maximum log product of a transversal.
+# ============================================================
+
+function MatrixCovers.cover_transversal_jump(::PowerMean{p}, A; target::Real=1e-2) where {p}
+    axr, axc = axes(A, 1), axes(A, 2)
+    length(axr) == length(axc) || throw(ArgumentError("cover_transversal_jump requires a square matrix"))
+    T = float(real(eltype(A)))
+    m = n = length(axr)
+    σ, α0, β0, logπ = MatrixCovers._max_product_transversal(MatrixCovers._row_support(A, T), axc,
+                                                            :cover_transversal_jump)
+    ei, ej, elog = _edge_list(A, T)
+    nza = _degrees(ei, m)
+    nzb = _degrees(ej, n)
+    tp = T(target)^p
+    model = _ipopt_model()
+    JuMP.set_attribute(model, "tol", 1e-11)
+    @variable(model, α[i=1:m], start = α0[i])
+    @variable(model, β[j=1:n], start = β0[j])
+    off = Int[]
+    for e in eachindex(ei, ej, elog)
+        if ej[e] == σ[ei[e]]
+            @constraint(model, α[ei[e]] + β[ej[e]] == elog[e])
+        else
+            @constraint(model, α[ei[e]] + β[ej[e]] >= elog[e])
+            push!(off, e)
+        end
+    end
+    if !isempty(off)
+        @objective(model, Min,
+            sum((exp(-p * (α[ei[e]] + β[ej[e]] - elog[e])) - 1) / p + tp * (α[ei[e]] + β[ej[e]] - elog[e])
+                for e in off))
+    end
+    _pin_gauge!(model, α, β, A, nza, nzb)
+    JuMP.optimize!(model)
+    check_solved(model, "cover_transversal_jump")
+    a = similar(Array{T}, axr)
+    b = similar(Array{T}, axc)
+    for (i, k) in enumerate(axr)
+        a[k] = exp(JuMP.value(α[i]))
+    end
+    for (j, k) in enumerate(axc)
+        b[k] = exp(JuMP.value(β[j]))
+    end
+    MatrixCovers._balance_cover!(a, b, A)
+    return MatrixCovers.inflate_feasible!(a, b, A), logπ
+end
+
 # ============================================================
 # Hard cover: cover_min!(::AbsLinear{p}, a, b, A)
 # Asymmetric hard-cover model in row and column log scales. The model pins the

ext/MatrixCoversUnitfulExt.jl

@@ -323,7 +323,9 @@ for P in PENALTIES
     end
 end
 
-# `PowerMean` is soft-only and parametric in its exponent.
+# `PowerMean` is parametric in its exponent.
+MC.cover_transversal(ϕ::PowerMean, A::QMatrix; kwargs...) = asym(MC.cover_transversal, A, ϕ; kwargs...)
+MC.cover_transversal!(ϕ::PowerMean, a::QVector, b::QVector, A::QMatrix; kwargs...) = asym!(MC.cover_transversal!, a, b, A, ϕ; kwargs...)
 MC.soft_symcover(ϕ::PowerMean, A::QMatrix; kwargs...) = sym(MC.soft_symcover, A, ϕ; kwargs...)
 MC.soft_symcover!(ϕ::PowerMean, a::QVector, A::QMatrix; kwargs...) = symstart!(MC.soft_symcover!, a, A, ϕ; kwargs...)
 MC.soft_cover(ϕ::PowerMean, A::QMatrix; kwargs...) = asym(MC.soft_cover, A, ϕ; kwargs...)

src/MatrixCovers.jl

@@ -18,7 +18,7 @@ export cover_transversal, cover_transversal!
 
 # `public` is parsed as a keyword only from Julia 1.11; this package supports 1.10.
 @static if VERSION >= v"1.11"
-    eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym, scalar_type, SupportComponents, support_components, ncomponents, rowcomponent, colcomponent, SolverFailure, AbstractCoverStats, AugmentedLagrangianStats, ActiveSetStats, LeastSquaresStats, PowerMeanStats, ExternalSolverStats, MultistartStats, converged, iterations, residual, tolerance"))
+    eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym, scalar_type, SupportComponents, support_components, ncomponents, rowcomponent, colcomponent, SolverFailure, AbstractCoverStats, AugmentedLagrangianStats, ActiveSetStats, LeastSquaresStats, PowerMeanStats, BalancedCoverStats, ExternalSolverStats, MultistartStats, converged, iterations, residual, tolerance"))
 end
 
 include("penalties.jl")
@@ -37,6 +37,8 @@ include("tree_cholesky.jl")    # fill-free factorization of forest-patterned mat
 include("minimal_covers.jl")
 include("active_set.jl")       # exact solution of the AbsLog{2} problem from a near-optimal point
 include("transversal.jl")      # transversal-tight minimal covers
+include("balanced.jl")         # transversal-tight covers pushed toward a target ratio
+include("maxbalance.jl")       # max-balancing within the fully indecomposable blocks
 include("transversal_heuristic.jl") # the transversal-tight start of `cover`
 include("powermean_newton.jl") # Newton refinement of the PowerMean soft covers
 include("sparse_support.jl")  # sparse traversal and the sparse solver defaults