From b9879dcc98c2ae7a3599df61cc6fd820440e1b96 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 23 Oct 2026 11:18:29 -0500 Subject: [PATCH 1/5] Factor out the transversal-tight setup and finish `_transversal_system` builds the transversal-tight difference system (maximum-product transversal, transversal log magnitudes, the edge list of off-transversal entries, tight components, and the feasible start), and `_contract_system` reduces it to one unknown per tight component. `_finish_cover` applies the balance convention, checks that the scales are representable, exponentiates, and certifies the cover. `_cover_min_abslog2` now calls these helpers; its results are unchanged. Assisted-by: Claude Opus 5.5 --- src/minimal_covers.jl | 66 +++++++++++++++-------------------------- src/transversal.jl | 68 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 91 insertions(+), 43 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 8fd8aa8..4d471e4 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -2163,64 +2163,30 @@ function _cover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::In Grid{T}(C) else G = _row_support(A, T) - edges = Tuple{Int,Int}[] - cvals = T[] if tt - σ, α0, _, logπ = _max_product_transversal(G, axc, fname) - if transversal isa AbstractVector - σ = _check_transversal(transversal, G, axc, logπ, fname) - end - # Tightness on the transversal fixes `β[σ[i]] = ctrans[σ[i]] - α[i]`, - # so entry `(i, j)` off the transversal has residual - # `α[i] - α[k] - (log|A[i,j]| - ctrans[j])` with `k = τ[j]`. - τ = invperm(σ) - ctrans = zeros(T, n) - for (ip, i) in enumerate(axr) - for s in _slots(G, i) - G.idx[s] - first(axc) + 1 == σ[ip] && (ctrans[σ[ip]] = log(G.val[s])) - end - end - for (ip, i) in enumerate(axr) - for s in _slots(G, i) - jp = G.idx[s] - first(axc) + 1 - jp == σ[ip] && continue - push!(edges, (ip, τ[jp])) - push!(cvals, log(G.val[s]) - ctrans[jp]) - end - end - dsupp = EdgeList{T}(edges, cvals, -oneunit(T)) - # Solve for one scale per component of rows forced to move together. - comp, offs, nc = _tight_components(edges, cvals, α0) - y0 = zeros(T, nc) - for p in 1:m - y0[comp[p]] = α0[p] - offs[p] - end - x0 = y0 + ts = _transversal_system(G, axc, T, fname, transversal) + (; σ, τ, ctrans, dsupp, comp, offs, nc) = ts + x0 = ts.y0 if use_woodbury # Entry `(i, j)` constrains components `comp[i]` and `comp[τ[j]]`. Cd = fill(T(-Inf), m, n) - ke = 0 # index into `edges`, which follows the same traversal + ke = 0 # index into `dsupp.edges`, which follows the same traversal for (ip, i) in enumerate(axr) for s in _slots(G, i) jp = G.idx[s] - first(axc) + 1 jp == σ[ip] && continue ke += 1 comp[ip] == comp[τ[jp]] && continue - Cd[ip, jp] = cvals[ke] - offs[ip] + offs[τ[jp]] + Cd[ip, jp] = dsupp.cvals[ke] - offs[ip] + offs[τ[jp]] end end DiffGrid{T}(Cd, comp, comp[τ]) else - cedges = Tuple{Int,Int}[] - ccvals = T[] - for (e, (p, q)) in enumerate(edges) - comp[p] == comp[q] && continue - push!(cedges, (comp[p], comp[q])) - push!(ccvals, cvals[e] - offs[p] + offs[q]) - end - EdgeList{T}(cedges, ccvals, -oneunit(T)) + _contract_system(dsupp, comp, offs) end else + edges = Tuple{Int,Int}[] + cvals = T[] for (ip, i) in enumerate(axr) for s in _slots(G, i) push!(edges, (ip, m + G.idx[s] - first(axc) + 1)) @@ -2313,6 +2279,20 @@ function _cover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::In if tt x = _transversal_scales(x[comp] .+ offs, dsupp, σ, ctrans, boost) end + a, b = _finish_cover(A, x, nzrow, nzcol, T, fname, boost) + return a, b, stats +end + +# Cover factors of `A` from the stacked log scales `x = (α; β)` (row positions, +# then column positions), with the balance convention applied independently to +# each support component, and, with `boost`, certified to cover `A`. `nzrow` +# and `nzcol` count the support entries of each row and column position. +function _finish_cover(A::AbstractMatrix, x::Vector{T}, nzrow::Vector{Int}, nzcol::Vector{Int}, + ::Type{T}, fname::Symbol, boost::Bool) where {T} + axr, axc = axes(A, 1), axes(A, 2) + m, n = length(axr), length(axc) + hasrow = nzrow .> 0 + hascol = nzcol .> 0 # Apply the balance convention independently to each support component. rowcomp, colcomp, ncomp, _, _ = _support_components(A) Lα = zeros(T, ncomp) @@ -2347,7 +2327,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::In b[j] = hascol[jp] ? exp(lβ[jp]) : zero(T) end boost && _certify_cover!(a, b, A, fname) - return a, b, stats + return a, b end # Stacked log scales `(α; β)` from the row scales `α` of the difference system diff --git a/src/transversal.jl b/src/transversal.jl index ae3a58e..c216e8b 100644 --- a/src/transversal.jl +++ b/src/transversal.jl @@ -301,6 +301,74 @@ end @noinline _throw_singular(fname) = throw(ArgumentError("$fname requires a structurally nonsingular matrix, but no transversal of nonzero entries exists")) +# The transversal-tight difference system of the square support grouped by row +# in `G`, whose partners are the column indices `axc`. A cover tight on a +# maximum-product transversal `σ` has `β[σ[i]] = ctrans[σ[i]] - α[i]`, so entry +# `(i, j)` off the transversal, with `k = τ[j]`, has residual +# `α[i] - α[k] - (log|A[i,j]| - ctrans[j])`, which must be nonnegative. +# +# `transversal` is `true` to use the transversal found by +# `_max_product_transversal`, or a vector naming the column position matched to +# each row position, which must also be a maximum-product transversal. Returns a +# named tuple in positions: +# +# - `σ`, `τ = invperm(σ)`: the transversal and its inverse; +# - `ctrans`: the log magnitude of the transversal entry of each column; +# - `dsupp`: the residuals as an `EdgeList` with `qsign = -1`, one entry per +# support entry off the transversal, in the order of `G`; +# - `comp`, `offs`, `nc`: the components of `_tight_components`, so that +# `α[p] = y[comp[p]] + offs[p]` for component scales `y`; +# - `y0`: the component scales of the transversal duals, which are feasible and +# tight on `σ`. +function _transversal_system(G::GroupedSupport, axc::AbstractUnitRange, ::Type{T}, fname::Symbol, + transversal::Union{Bool,AbstractVector{<:Integer}}) where {T} + axr = G.ax + m = length(axr) + n = length(axc) + σ, α0, _, logπ = _max_product_transversal(G, axc, fname) + if transversal isa AbstractVector + σ = _check_transversal(transversal, G, axc, logπ, fname) + end + τ = invperm(σ) + ctrans = zeros(T, n) + for (ip, i) in enumerate(axr) + for s in _slots(G, i) + G.idx[s] - first(axc) + 1 == σ[ip] && (ctrans[σ[ip]] = log(G.val[s])) + end + end + edges = Tuple{Int,Int}[] + cvals = T[] + for (ip, i) in enumerate(axr) + for s in _slots(G, i) + jp = G.idx[s] - first(axc) + 1 + jp == σ[ip] && continue + push!(edges, (ip, τ[jp])) + push!(cvals, log(G.val[s]) - ctrans[jp]) + end + end + dsupp = EdgeList{T}(edges, cvals, -oneunit(T)) + comp, offs, nc = _tight_components(edges, cvals, α0) + y0 = zeros(T, nc) + for p in 1:m + y0[comp[p]] = α0[p] - offs[p] + end + return (; σ, τ, ctrans, dsupp, comp, offs, nc, y0) +end + +# The difference system `dsupp` on the component scales `y`, where row `p` has +# scale `y[comp[p]] + offs[p]`. Entries within one component have fixed +# residuals and are dropped. +function _contract_system(dsupp::EdgeList{T}, comp::Vector{Int}, offs::Vector{T}) where {T} + cedges = Tuple{Int,Int}[] + ccvals = T[] + for (e, (p, q)) in enumerate(dsupp.edges) + comp[p] == comp[q] && continue + push!(cedges, (comp[p], comp[q])) + push!(ccvals, dsupp.cvals[e] - offs[p] + offs[q]) + end + return EdgeList{T}(cedges, ccvals, -oneunit(T)) +end + # Rows joined by entries that lie on some maximum-product transversal. # # In the difference system on the row log scales `α`, entry `e = (p, q)` has From e175f620c1267b1e865846e7dfbaf9f7043d1c44 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 24 Oct 2026 05:48:32 -0500 Subject: [PATCH 2/5] Add cover_transversal with a PowerMean penalty MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `cover_transversal(PowerMean{p}(), A; target)` returns the hard cover tight on the maximum-product transversals that minimizes ∑ (r^p - 1 - p*log(r/target)*target^p)/p over the coverage ratios r of the other entries, pushing them toward `target`. The solver is a primal active-set Newton method on the contracted row-scale difference system: rows joined by zero-slack entries move together as groups, Newton steps solve weighted graph-Laplacian systems on the groups (dense or sparse Cholesky, or conjugate gradients over the fill budget), a ratio test and Armijo backtracking keep every iterate feasible, and at stationary points a maximum-flow test on each group's active entries either certifies the KKT conditions or yields a set of rows whose scales decrease, releasing entries. Every iterate is a hard cover tight on the transversal, so an unconverged solve warns and still returns one. Statistics are returned as `BalancedCoverStats`. The objective's terms can span many orders of magnitude, so the solver works at the resolution of each group's terms. The gradient and the line-search slope sum the `target^p` part by group with integer counts, so that the slack changes of a group with as many entries leaving as entering cancel exactly. A ridge for a Newton matrix that roundoff makes indefinite is relative to each diagonal entry, so it does not swamp the curvature of groups whose terms are all small. In each weakly connected part, the group pinned to remove the common shift is the one whose gradient has the largest terms: the pinned group's gradient carries the others' roundoff, which is then smallest relative to its own terms. A Newton step whose slope is nonnegative but within the roundoff of the groups' gradients is at the resolution of the arithmetic: the groups are taken to be stationary and the multiplier test decides whether to continue, so `converged` can hold with `gradnorm > tol`. `cover_transversal(AbsLog{2}(), A)` and its mutating form forward to the existing solver. The Ipopt extension adds the reference `cover_transversal_jump(::PowerMean, A; target)`, and the Unitful extension forwards the new methods. test/balanced.jl checks the maximum-flow multiplier test on small digraphs, the line search where slack changes cancel, where the slope is at roundoff, and where an entry at zero slack blocks the step, the covariance of a solve whose groups' terms differ greatly in size, and the solver against the Ipopt reference (committed 5×5 matrices, random sparse patterns, 0/1 matrices), against an extended-precision solve for p = 4, the KKT identity across a block-triangular cut, a group holding an inactive entry, the minimax ratio for large p, covariance, permutation equivariance, independence of the transversal used, sparse and conjugate-gradient solves, element types, generic axes, the iteration limit, and the argument errors. The manual describes the new method. Assisted-by: Claude Opus 5.5 --- docs/src/index.md | 5 +- docs/src/objectives.md | 10 +- docs/src/solvers.md | 3 +- ext/MatrixCoversIpoptExt.jl | 52 +++ ext/MatrixCoversUnitfulExt.jl | 4 +- src/MatrixCovers.jl | 3 +- src/balanced.jl | 769 ++++++++++++++++++++++++++++++++++ src/minimal_covers.jl | 3 +- src/stats.jl | 64 +++ src/transversal.jl | 7 +- test/balanced.jl | 440 +++++++++++++++++++ test/runtests.jl | 4 +- test/unitful.jl | 3 + 13 files changed, 1357 insertions(+), 10 deletions(-) create mode 100644 src/balanced.jl create mode 100644 test/balanced.jl diff --git a/docs/src/index.md b/docs/src/index.md index 6d24afe..9fd42fd 100644 --- a/docs/src/index.md +++ b/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 diff --git a/docs/src/objectives.md b/docs/src/objectives.md index d3d8c61..e1a3c14 100644 --- a/docs/src/objectives.md +++ b/docs/src/objectives.md @@ -58,7 +58,7 @@ 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) | — | | [`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 | @@ -66,7 +66,10 @@ 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. 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 +81,9 @@ 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. - `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. diff --git a/docs/src/solvers.md b/docs/src/solvers.md index 457e765..5c1b120 100644 --- a/docs/src/solvers.md +++ b/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) diff --git a/ext/MatrixCoversIpoptExt.jl b/ext/MatrixCoversIpoptExt.jl index 2cb0ee8..c5f6882 100644 --- a/ext/MatrixCoversIpoptExt.jl +++ b/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 diff --git a/ext/MatrixCoversUnitfulExt.jl b/ext/MatrixCoversUnitfulExt.jl index 6653fbf..09dff7e 100644 --- a/ext/MatrixCoversUnitfulExt.jl +++ b/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...) diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index 2fd892d..431e2b1 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -18,7 +18,7 @@ export cover_transversal, cover_transversal! # `public` is parsed as a keyword only from Julia 1.11; this package supports 1.10. @static if VERSION >= v"1.11" - eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym, scalar_type, SupportComponents, support_components, ncomponents, rowcomponent, colcomponent, SolverFailure, AbstractCoverStats, AugmentedLagrangianStats, ActiveSetStats, LeastSquaresStats, PowerMeanStats, 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,7 @@ 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("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 diff --git a/src/balanced.jl b/src/balanced.jl new file mode 100644 index 0000000..63ac4f6 --- /dev/null +++ b/src/balanced.jl @@ -0,0 +1,769 @@ +# Transversal-tight hard covers that push the coverage ratios of the entries off +# the transversal toward a target ratio. + +# Stopping tolerance (see `BalancedCoverStats`). +const BALANCED_TOL = 1e-9 +# Sufficient-decrease constant and the most halvings a line search may try. +const BALANCED_ARMIJO = 1e-4 +const BALANCED_NBACKTRACK = 100 +# A first trial step changes no slack by more than this many units of `1/p`, so +# that `exp(-p d)` stays in range. +const BALANCED_MAXSTEP = 64 +# Roundoff allowance, in units of `eps`, for a change in the objective too small +# to resolve. +const BALANCED_NOISE = 64 +# Newton systems on at most this many unknowns, or at least an eighth full, are +# factored densely. +const BALANCED_DENSE_MAX = 128 +# Ridge retries for a Newton matrix that roundoff makes indefinite; each retry +# scales the ridge by `BALANCED_RIDGE_GROWTH`. The ridge is relative to each +# diagonal entry, so that it does not swamp the curvature of a group whose +# entries all have small `f''`. +const BALANCED_NRIDGE = 4 +const BALANCED_RIDGE_GROWTH = 1e4 + +""" + a, b, stats = cover_transversal(ϕ::PowerMean{p}, A; target=1e-2, tol=1e-9, maxiter, fillbudget, flopbudget) + +Return the hard cover of the square matrix `A` that is tight on its +maximum-product transversals and pushes the coverage ratios +`r = |A[i,j]| / (a[i]*b[j])` of the other nonzero entries toward `target`: it +minimizes + + ∑ (r^p - 1 - p*log(r/target)*target^p) / p, + +the sum taken over the nonzero entries off the transversals, among the hard +covers with `prod(a) * prod(b) == π*` (see [`cover_transversal(A)`](@ref +cover_transversal)). Each term is `target^p * φ(r/target)` up to a constant, with +`φ` the [`PowerMean{p}`](@ref) penalty, so it is least at `r = target`; the hard +constraint caps every ratio at one, and the transversal entries have ratio one. +`0 < target <= 1` and `0 < p < ∞` are required; otherwise an `ArgumentError` is +thrown. `cover_transversal(AbsLog{2}(), A)` is `cover_transversal(A)`. + +The objective is strictly convex in the products `a[i]*b[j]`, so they are +unique on the support of `A`. As for `cover_transversal(A)`, they do not depend +on which maximum-product transversal is used, they are scale-covariant, and the +factors use the balance convention of [`cover_min`](@ref). + +With `target = 1`, every ratio is pushed up to one, as far as the constraints +allow. With small `target`, the ratios within each block of a block-triangular +form of `A` balance in the sense of the `p`-power mean, and the entries joining +different blocks are pushed toward `target`; at the minimizer, the `p`-power +mean of the ratios of the entries crossing such a cut equals `target` whenever +none of them has ratio one. Large `p` approaches the cover that minimizes the +largest ratio within each block. + +`stats` is a [`BalancedCoverStats`](@ref MatrixCovers.BalancedCoverStats). The +solver warns when [`MatrixCovers.converged(stats)`](@ref MatrixCovers.converged) +is `false`; the result is then still a hard cover tight on the maximum-product +transversals. + +# Extended help + +As for `cover_transversal(A)`, tightness on a maximum-product transversal `σ` +leaves the row log scales as unknowns, with one constraint `d ≥ 0` on the log +slack `d = -log(r)` of each entry off the transversal. Starting from the +transversal's dual variables, which are feasible, a primal active-set Newton +method keeps every iterate feasible. Entries at zero slack are *active*, and the +rows they join form groups whose scales move together; the entries at zero +slack at the start are active from the outset. Each Newton step solves a +weighted graph-Laplacian system on the groups, a ratio test stops a step where +another entry reaches zero slack, and a backtracking line search ensures +descent. Where the groups are stationary, a maximum-flow problem on each group's +active entries either proves the KKT conditions or yields a set of rows whose +scales can decrease, releasing the active entries that enter it. + +`tol` (default `1e-9`) is the stopping tolerance: the groups are stationary when +the derivative of the objective with respect to each group's scale is at most +`tol` times the sum of the magnitudes of its terms, or when roundoff keeps the +Newton step from reducing the objective: the step changes no log slack by more +than `tol`, or the directional derivative of the objective along it is +nonnegative but within the roundoff of the groups' derivatives. `maxiter` +(default `100 + 10n` for `n×n` `A`) bounds the Newton steps plus the changes of +the active set. +`fillbudget` and `flopbudget` are as for [`cover_min`](@ref): a Newton system +whose sparse Cholesky factor would exceed them is solved by conjugate +gradients. The Newton systems have condition numbers of order `target^(-p)` +times that of the graph, and of order `r^(-p)` where cycles of entries hold +ratios `r` far below `target`, so a very small `target^p` or a large `p` +degrades the linear solves and the convergence. + +`A` must be square and structurally nonsingular. Narrower element types than +`Float64` are computed in `Float64` and narrowed, with the coverage certified +afterwards. + +See also: [`cover_transversal!`](@ref), [`PowerMean`](@ref). +""" +function cover_transversal(ϕ::PowerMean{p}, A::AbstractMatrix; target::Real=1e-2, + tol::Real=BALANCED_TOL, maxiter::Union{Integer,Nothing}=nothing, + fillbudget::Real=LSQR_FILL_BUDGET, + flopbudget::Real=LSQR_FLOP_BUDGET) where {p} + a, b, nt = _cover_transversal_balanced(A, p, target; tol, maxiter, fillbudget, flopbudget) + return a, b, BalancedCoverStats{_auglag_stats_type(A)}(nt) +end + +cover_transversal(::AbsLog{2}, A::AbstractMatrix; kwargs...) = cover_transversal(A; kwargs...) + +""" + a, b, stats = cover_transversal!(ϕ::PowerMean, a, b, A; kwargs...) + +Mutating counterpart of [`cover_transversal(ϕ, A)`](@ref cover_transversal): +writes the cover into `a` 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!(ϕ::PowerMean, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) + _check_cover_axes(a, b, A) + anew, bnew, stats = cover_transversal(ϕ, A; kwargs...) + a .= anew + b .= bnew + return a, b, stats +end + +cover_transversal!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) = + cover_transversal!(a, b, A; kwargs...) + +# Worker for `cover_transversal(::PowerMean{p}, A; target)`, returning `(a, b, nt)` +# with the statistics as a named tuple. `transversal` is as for +# `_transversal_system`. +function _cover_transversal_balanced(A::AbstractMatrix, p::Real, target::Real; tol::Real=BALANCED_TOL, + maxiter::Union{Integer,Nothing}=nothing, + fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET, + transversal::Union{Bool,AbstractVector{<:Integer}}=true, + fname::Symbol=:cover_transversal) + isfinite(p) || throw(ArgumentError("$fname requires a finite exponent, got PowerMean{$p}")) + 0 < target <= 1 || throw(ArgumentError("$fname: target must satisfy 0 < target <= 1, got $target")) + tol > 0 || throw(ArgumentError("$fname: tol must be positive, got $tol")) + maxiter === nothing || maxiter >= 0 || + throw(ArgumentError("$fname: maxiter must be nonnegative, got $maxiter")) + T = float(real(eltype(A))) + # The tolerances require at least Float64 resolution. + if eps(T) > eps(Float64) + a64, b64, nt = _cover_transversal_balanced(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A), + p, target; tol, maxiter, fillbudget, flopbudget, + transversal, fname) + # Narrowing rounds to nearest and so can round a product below its entry. + a, b = T.(a64), T.(b64) + _certify_cover!(a, b, A, fname) + return a, b, nt + end + pT = T(p) + tp = T(target)^pT + tp > 0 || throw(ArgumentError("$fname: target^p = $target^$p underflows in $T")) + axr, axc = axes(A, 1), axes(A, 2) + or, oc = first(axr) - 1, first(axc) - 1 + nzrow = zeros(Int, length(axr)) + nzcol = zeros(Int, length(axc)) + foreach_support(A) do i, j, v + nzrow[i-or] += 1 + nzcol[j-oc] += 1 + end + G = _row_support(A, T) + ts = _transversal_system(G, axc, T, fname, transversal) + sys = _contract_system(ts.dsupp, ts.comp, ts.offs) + maxit = maxiter === nothing ? 100 + 10 * ts.nc : Int(maxiter) + y, nt = _balanced_newton!(copy(ts.y0), sys, pT, tp; tol=T(tol), maxiter=maxit, fillbudget, flopbudget) + _warn_balanced_unconverged(fname, nt, maxit) + # The uniform shift of `_transversal_scales` absorbs roundoff violations. + x = _transversal_scales(y[ts.comp] .+ ts.offs, ts.dsupp, ts.σ, ts.ctrans, true) + a, b = _finish_cover(A, x, nzrow, nzcol, T, fname, true) + return a, b, nt +end + +function _check_cover_axes(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) + 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)))")) + return nothing +end + +function _warn_balanced_unconverged(fname::Symbol, nt::NamedTuple, maxiter::Int) + nt.converged && return nothing + reason = nt.stalled ? + "the line search could not reduce the objective; a larger `target` or smaller `p` improves the conditioning of the Newton systems" : + "it reached maxiter=$maxiter; increase `maxiter`" + @warn "$fname: the active-set Newton iteration stopped after $(nt.nnewton) Newton steps with relative reduced gradient $(nt.gradnorm) (tolerance $(nt.tol)) because $reason. The result covers `A` and is tight on a maximum-product transversal, but may not minimize the objective." + return nothing +end + +# The objective sums, over the entries off the transversal, `f` of the log slack `d`: +# +# f(d) = (exp(-p d) - 1)/p + tp*d, f'(d) = tp - exp(-p d), f''(d) = p exp(-p d), +# +# with `tp = target^p`; `f` is least at `d = -log(target)`. + +# Gradient `g[k]` of the objective with respect to a common shift of the scales +# of the rows `v` with `label[v] == k`, over the entries that join different +# labels when `between` is true and over all entries otherwise; the entries +# `e` with both rows labeled `k` cancel in the former. `f'` is `tp - w` with +# `w = exp(-p d)`, so `g[k] = tp * cnt[k] - (sum of ±w)`, where `cnt[k]` counts +# the entries leaving `k` minus those entering it. Summing the two parts +# separately avoids cancellation between them; `s[k] = tp * |cnt[k]| + (sum of +# w)` bounds the magnitudes summed and so scales the roundoff in `g[k]`. +function _balanced_gradient!(g, s, cnt, d, label, P, Q, p::T, tp::T, between::Bool) where {T} + fill!(g, zero(T)) + fill!(s, zero(T)) + fill!(cnt, 0) + for e in eachindex(d) + G, H = label[P[e]], label[Q[e]] + between && G == H && continue + w = exp(-p * d[e]) + g[G] -= w + g[H] += w + s[G] += w + s[H] += w + cnt[G] += 1 + cnt[H] -= 1 + end + for k in eachindex(g, s, cnt) + g[k] += tp * cnt[k] + s[k] += tp * abs(cnt[k]) + end + return g +end + +# Roundoff allowance for a zero slack of entry `e` at the scales `y`. +_balanced_tight(y, C, P, Q, e) = + TIGHT_ULPS * eps(eltype(y)) * max(oneunit(eltype(y)), abs(y[P[e]]), abs(y[Q[e]]), abs(C[e])) + +# Label the connected components of the graph on `1:nv` formed by the edges +# `(P[e], Q[e])` with `keep(e)`, numbered by their first vertex. Returns the +# number of components. +function _edge_components!(lab::Vector{Int}, keep, P::Vector{Int}, Q::Vector{Int}) + nv = length(lab) + parent = collect(1:nv) + function find(v) + while parent[v] != v + parent[v] = parent[parent[v]] + v = parent[v] + end + return v + end + for e in eachindex(P, Q) + keep(e) || continue + rp, rq = find(P[e]), find(Q[e]) + rp == rq || (parent[max(rp, rq)] = min(rp, rq)) + end + fill!(lab, 0) + k = 0 + for v in 1:nv + r = find(v) + lab[r] == 0 && (lab[r] = (k += 1)) + lab[v] = lab[r] + end + return k +end + +# Feasibility of a nonnegative flow on the uncapacitated arcs `tails[k] → heads[k]` +# among the vertices `1:nv` whose divergence (outflow minus inflow) at `v` is +# `g[v]`. Solved as a maximum flow from a source joined to each `v` with +# `g[v] > 0` (capacity `g[v]`) to a sink joined from each `v` with `g[v] < 0` +# (capacity `-g[v]`), by shortest augmenting paths (Edmonds–Karp). +# +# Returns `(deficit, S)`: `deficit` is the positive supply the maximum flow +# leaves unrouted, zero exactly when such a flow exists, and `S` marks the +# vertices reachable from the source in the final residual graph. No arc leaves +# `S`, no flow enters it, and `deficit == sum(g[S])` up to roundoff. +function _flow_deficit(nv::Int, tails::Vector{Int}, heads::Vector{Int}, g::AbstractVector{T}) where {T} + length(g) == nv || throw(DimensionMismatch("g has length $(length(g)) but there are $nv vertices")) + src, snk = nv + 1, nv + 2 + to = Int[] + cap = T[] + from = Int[] + function addarc!(u, v, c) + push!(from, u); push!(to, v); push!(cap, c) # arc `2k-1` + push!(from, v); push!(to, u); push!(cap, zero(T)) # its reverse, `2k` + end + for k in eachindex(tails, heads) + addarc!(tails[k], heads[k], T(Inf)) + end + supply = zero(T) + for v in 1:nv + if g[v] > 0 + addarc!(src, v, g[v]) + supply += g[v] + elseif g[v] < 0 + addarc!(v, snk, -g[v]) + end + end + # Arcs leaving each node, in compressed form. + ptr = zeros(Int, nv + 3) + for u in from + ptr[u+1] += 1 + end + ptr[1] = 1 + cumsum!(ptr, ptr) + adj = Vector{Int}(undef, length(from)) + cursor = ptr[1:end-1] + for (a, u) in enumerate(from) + adj[cursor[u]] = a + cursor[u] += 1 + end + pred = zeros(Int, nv + 2) # arc reaching each node in the search; 0 = unreached + queue = Int[] + flow = zero(T) + while true + fill!(pred, 0) + empty!(queue) + push!(queue, src) + pred[src] = -1 + head = 1 + while head <= length(queue) && pred[snk] == 0 + u = queue[head] + head += 1 + for s in ptr[u]:ptr[u+1]-1 + a = adj[s] + v = to[a] + (pred[v] == 0 && cap[a] > 0) || continue + pred[v] = a + push!(queue, v) + end + end + pred[snk] == 0 && break + bottleneck = T(Inf) + v = snk + while v != src + a = pred[v] + bottleneck = min(bottleneck, cap[a]) + v = from[a] + end + v = snk + while v != src + a = pred[v] + cap[a] -= bottleneck + cap[isodd(a) ? a + 1 : a - 1] += bottleneck + v = from[a] + end + flow += bottleneck + end + S = BitVector(pred[v] != 0 for v in 1:nv) + return max(supply - flow, zero(T)), S +end + +# Primal active-set Newton method for +# +# minimize ∑ₑ f(y[P[e]] - y[Q[e]] - C[e]) subject to y[P[e]] - y[Q[e]] - C[e] >= 0, +# +# over the component scales `y` of the contracted transversal-tight system `sys` +# (edges `(P[e], Q[e])`, constants `C[e]`), from a feasible `y`, which is +# updated in place. Every iterate is feasible up to roundoff. +# +# Active entries are held at zero slack. Groups are the connected components of +# the active entries; a step moves all scales of a group together, so an entry +# within a group keeps its slack. Newton steps on the group scales use the +# Laplacian of the entries between groups, weighted by `f''`. At a point where +# each group's reduced gradient vanishes, the KKT conditions require +# multipliers `μ >= 0` on the active entries of each group with divergence equal +# to the gradient `g[v]` of the objective at each of its rows. `g[v]` sums `f'` +# over every entry at row `v`: the inactive ones, including those within the +# group, and the active ones, where `f'(0) = tp - 1`. `_flow_deficit` decides +# whether the multipliers exist. If they do not, its set +# `S` has no active entry leaving it and `sum(g[S]) > 0`, so lowering the scales +# of `S` decreases the objective; the active entries entering `S` are released. +# +# Returns `(y, nt)`, with `nt` the fields of `BalancedCoverStats`. +function _balanced_newton!(y::Vector{T}, sys::EdgeList{T}, p::T, tp::T; tol::T, maxiter::Int, + fillbudget::Real=LSQR_FILL_BUDGET, + flopbudget::Real=LSQR_FLOP_BUDGET) where {T} + nc = length(y) + ne = length(sys.edges) + P = first.(sys.edges) + Q = last.(sys.edges) + C = sys.cvals + d = zeros(T, ne) + Δd = zeros(T, ne) + # Entries at zero slack at the start are active; the multiplier test + # releases those that should not be. + active = falses(ne) + for e in 1:ne + y[P[e]] - y[Q[e]] - C[e] <= _balanced_tight(y, C, P, Q, e) && (active[e] = true) + end + gid = zeros(Int, nc) + wcomp = zeros(Int, nc) + nw = _edge_components!(wcomp, _ -> true, P, Q) + gv = zeros(T, nc) # gradient over the inactive entries, by row + sv = zeros(T, nc) # sum of the magnitudes of its terms + δv = zeros(T, nc) + gG = T[] + sG = T[] + cnt = Int[] + uidx = Int[] + ng = 0 + nu = 0 + regroup = repin = true + cache = Ref{Any}(nothing) # the sparse factorization and its pattern, for reuse + nnewton = nfactor = nreleases = iters = 0 + nactivations = count(active) + converged = stalled = false + gnorm = T(Inf) + step = T(NaN) # largest change of a slack in the last Newton step + atfloor = false # the last Newton step was below the resolution of the arithmetic + linsolve = :none + while true + if regroup + ng = _edge_components!(gid, e -> active[e], P, Q) + resize!(gG, ng) + resize!(sG, ng) + resize!(cnt, ng) + regroup = false + repin = true + end + for e in 1:ne + d[e] = y[P[e]] - y[Q[e]] - C[e] + end + # Reduced gradient: the entries between groups. + _balanced_gradient!(gG, sG, cnt, d, gid, P, Q, p, tp, true) + if repin + # Pin, in each weak component, the group whose gradient has the + # largest terms; the others are the unknowns of the Newton systems. + # Once their gradients vanish, the pinned group's is minus their + # sum, and so holds their roundoff; that is smallest relative to + # the pinned group's own terms when those are the largest. + pin = zeros(Int, nw) + for v in 1:nc + G, k = gid[v], wcomp[v] + (pin[k] == 0 || sG[G] > sG[pin[k]]) && (pin[k] = G) + end + uidx = ones(Int, ng) + for k in 1:nw + uidx[pin[k]] = 0 + end + nu = 0 + for G in 1:ng + uidx[G] == 0 || (uidx[G] = (nu += 1)) + end + repin = false + end + gnorm = zero(T) + for G in 1:ng + sG[G] > 0 && (gnorm = max(gnorm, abs(gG[G]) / sG[G])) + end + if gnorm <= tol || atfloor + atfloor = false + # Entries at zero slack may carry multipliers; make them active. + k = 0 + for e in 1:ne + if !active[e] && d[e] <= _balanced_tight(y, C, P, Q, e) + active[e] = true + k += 1 + end + end + if k > 0 + (iters += 1) > maxiter && break + nactivations += k + regroup = true + continue + end + inS = _balanced_release_set!(gv, sv, d, active, gid, ng, P, Q, p, tp, tol) + if inS === nothing + converged = true + break + end + (iters += 1) > maxiter && break + for e in 1:ne + active[e] && !inS[P[e]] && inS[Q[e]] && (active[e] = false) + end + nreleases += 1 + regroup = true + for v in 1:nc + δv[v] = inS[v] ? -one(T) : zero(T) + end + status = _balanced_step!(y, δv, d, Δd, active, 1:nc, P, Q, C, p, tp, nothing) + (status === :stalled || status === :floor) && (stalled = true; break) + status isa Int && (nactivations += status) + continue + end + (iters += 1) > maxiter && break + linsolve = _balanced_direction!(δv, cache, d, gid, uidx, nu, gG, P, Q, p, tp, + fillbudget, flopbudget) + nnewton += 1 + nfactor += 1 + step = zero(T) + for e in 1:ne + step = max(step, abs(δv[P[e]] - δv[Q[e]])) + end + status = _balanced_step!(y, δv, d, Δd, active, gid, P, Q, C, p, tp, one(T)) + if status === :floor || status === :stalled + # A Newton step whose slope is zero to roundoff, or that changes no + # slack by more than `tol` but does not reduce the objective, is at + # the resolution of the arithmetic: the reduced gradient is zero to + # roundoff. + status === :floor || step <= tol || (stalled = true; break) + atfloor = true + elseif status isa Int + nactivations += status + status > 0 && (regroup = true) + end + end + return y, (; converged, stalled, tol, gradnorm=gnorm, step, nnewton, nfactor, nactivations, + nreleases, nactive=count(active), linsolve) +end + +# The multiplier test at a point where every group's reduced gradient vanishes. +# Returns `nothing` when multipliers exist for every group, and otherwise the +# union of the sets `S` of `_flow_deficit` over the groups that lack them. +function _balanced_release_set!(gv, sv, d, active, gid, ng, P, Q, p::T, tp::T, tol::T) where {T} + nc = length(gv) + _balanced_gradient!(gv, sv, zeros(Int, nc), d, 1:nc, P, Q, p, tp, false) + # Rows and active entries of each group, in compressed form. + vptr = zeros(Int, ng + 1) + for v in 1:nc + vptr[gid[v]+1] += 1 + end + vptr[1] = 1 + cumsum!(vptr, vptr) + verts = Vector{Int}(undef, nc) + local_index = Vector{Int}(undef, nc) + cursor = vptr[1:ng] + for v in 1:nc + G = gid[v] + verts[cursor[G]] = v + local_index[v] = cursor[G] - vptr[G] + 1 + cursor[G] += 1 + end + eptr = zeros(Int, ng + 1) + for e in eachindex(active) + active[e] && (eptr[gid[P[e]]+1] += 1) + end + eptr[1] = 1 + cumsum!(eptr, eptr) + aedges = Vector{Int}(undef, eptr[end] - 1) + cursor = eptr[1:ng] + for e in eachindex(active) + active[e] || continue + G = gid[P[e]] + aedges[cursor[G]] = e + cursor[G] += 1 + end + inS = nothing + for G in 1:ng + er = eptr[G]:eptr[G+1]-1 + isempty(er) && continue + vr = vptr[G]:vptr[G+1]-1 + tails = [local_index[P[aedges[k]]] for k in er] + heads = [local_index[Q[aedges[k]]] for k in er] + g = [gv[verts[k]] for k in vr] + ctol = tol * sum(k -> sv[verts[k]], vr) + deficit, S = _flow_deficit(length(vr), tails, heads, g) + deficit <= ctol && continue + inS === nothing && (inS = falses(nc)) + for (k, v) in enumerate(view(verts, vr)) + S[k] && (inS[v] = true) + end + end + return inS +end + +# Newton direction on the group scales, scattered to the rows in `δv`. Returns +# the kind of linear solve. +function _balanced_direction!(δv, cache, d, gid, uidx, nu, gG, P, Q, p::T, tp::T, + fillbudget, flopbudget) where {T} + rhs = zeros(T, nu) + for G in eachindex(uidx) + uidx[G] > 0 && (rhs[uidx[G]] = -gG[G]) + end + nb = 0 + for e in eachindex(d) + gid[P[e]] == gid[Q[e]] || (nb += 1) + end + weight(e) = p * exp(-p * d[e]) + x = zeros(T, nu) + kind = :none + if nu == 0 + elseif nu <= BALANCED_DENSE_MAX || 8 * nb >= nu^2 + H = zeros(T, nu, nu) + for e in eachindex(d) + G, K = uidx[gid[P[e]]], uidx[gid[Q[e]]] + gid[P[e]] == gid[Q[e]] && continue + w = weight(e) + G > 0 && (H[G, G] += w) + K > 0 && (H[K, K] += w) + G > 0 && K > 0 && (H[min(G, K), max(G, K)] -= w) + end + dmax = maximum(diag(H)) + dfloor = eps(T) * (dmax > 0 ? dmax : one(T)) + ridge = zero(T) + for _ in 0:BALANCED_NRIDGE + Hr = copy(H) + for k in 1:nu + Hr[k, k] += ridge * max(H[k, k], dfloor) + end + F = LinearAlgebra.cholesky!(Symmetric(Hr, :U); check=false) + if LinearAlgebra.issuccess(F) + ldiv!(x, F, rhs) + kind = :dense + break + end + ridge = _balanced_ridge(ridge) + end + kind === :none && throw(LinearAlgebra.PosDefException(1)) + else + I = Int[] + J = Int[] + V = T[] + dg = zeros(T, nu) + for e in eachindex(d) + gid[P[e]] == gid[Q[e]] && continue + G, K = uidx[gid[P[e]]], uidx[gid[Q[e]]] + w = weight(e) + G > 0 && (dg[G] += w) + K > 0 && (dg[K] += w) + if G > 0 && K > 0 + push!(I, min(G, K)); push!(J, max(G, K)); push!(V, -w) + end + end + dmax = maximum(dg) + dfloor = eps(T) * (dmax > 0 ? dmax : one(T)) + ridge = zero(T) + for _ in 0:BALANCED_NRIDGE + L = sparse([I; 1:nu], [J; 1:nu], [V; dg .+ ridge .* max.(dg, dfloor)], nu, nu) + F = _balanced_factor!(cache, L, fillbudget, flopbudget) + if F !== nothing + _laplacian_solve!(x, F, rhs) + kind = F isa SparseCholesky ? :cholesky : :cg + break + end + ridge = _balanced_ridge(ridge) + end + kind === :none && throw(LinearAlgebra.PosDefException(1)) + end + for v in eachindex(δv) + k = uidx[gid[v]] + δv[v] = k > 0 ? x[k] : zero(T) + end + return kind +end + +# The next ridge, relative to the diagonal, for a Newton matrix that roundoff +# made indefinite. +_balanced_ridge(ridge::T) where {T} = iszero(ridge) ? eps(T) : ridge * T(BALANCED_RIDGE_GROWTH) + +# Factor the pinned Laplacian `L` (upper triangle), reusing the symbolic analysis +# in `cache` when the pattern is unchanged. Returns `nothing` when `L` is not +# numerically positive definite. +function _balanced_factor!(cache, L::SparseMatrixCSC{Float64,Int}, fillbudget, flopbudget) + c = cache[] + try + if c !== nothing && c[2] == SparseArrays.getcolptr(L) && c[3] == rowvals(L) + return factorize!(c[1], L) + end + F = _laplacian_factor(L; fillbudget, flopbudget) + cache[] = F isa SparseCholesky ? (F, copy(SparseArrays.getcolptr(L)), copy(rowvals(L))) : nothing + return F + catch err + err isa LinearAlgebra.PosDefException || rethrow() + return nothing + end +end +_balanced_factor!(cache, L::SparseMatrixCSC, fillbudget, flopbudget) = _laplacian_factor(L) + +# Move the scales `y` along `δv`, which leaves the slack of every active entry +# unchanged, by a step that keeps every slack nonnegative and decreases the +# objective. `α0` is the first trial step (`nothing` for the minimizer of the +# quadratic model along `δv`), reduced by the ratio test. Returns the number of +# entries made active (an `Int`, possibly zero), `:floor` when the slope along +# `δv` is nonnegative but within the roundoff of the gradients of the labels, or +# `:stalled` when no step decreases the objective. Inactive entries at zero +# slack that `δv` would make negative are made active without moving. `δv` is +# constant on the sets of rows with equal `label`. +function _balanced_step!(y, δv, d, Δd, active, label, P, Q, C, p::T, tp::T, α0) where {T} + # `f'` and `f(d + h) - f(d)` split into the parts in `tp` and in `exp(-p d)`, + # summed separately: the first is `tp` times the sum `sΔ` of the slack + # changes. As in `_balanced_gradient!`, `sΔ` sums `δv` on each label times + # the count of inactive entries leaving it minus those entering it, so that + # the changes of a label with as many entries leaving as entering cancel + # exactly; a sum by entry leaves a rounding error proportional to their size. + δl = zeros(T, length(δv)) + cnt = zeros(Int, length(δv)) + for v in eachindex(δv) + δl[label[v]] = δv[v] + end + for v in eachindex(δv) + δv[v] == δl[label[v]] || error("internal error: the step is not constant on a label") + end + sw = gscale = curv = zero(T) + αmax = T(Inf) + Δmax = zero(T) + nblock = 0 + for e in eachindex(d) + Δd[e] = δv[P[e]] - δv[Q[e]] + if active[e] + iszero(Δd[e]) || error("internal error: the step changes the slack of an active entry") + continue + end + Δ = Δd[e] + w = exp(-p * d[e]) + cnt[label[P[e]]] += 1 + cnt[label[Q[e]]] -= 1 + sw += w * Δ + label[P[e]] == label[Q[e]] || (gscale += (tp + w) * (abs(δv[P[e]]) + abs(δv[Q[e]]))) + curv += p * w * Δ^2 + Δmax = max(Δmax, abs(Δ)) + if Δ < 0 + if d[e] <= _balanced_tight(y, C, P, Q, e) + nblock += 1 + else + αmax = min(αmax, d[e] / -Δ) + end + end + end + if nblock > 0 + for e in eachindex(d) + if !active[e] && Δd[e] < 0 && d[e] <= _balanced_tight(y, C, P, Q, e) + active[e] = true + end + end + return nblock + end + sΔ = sΔabs = zero(T) + for k in eachindex(cnt, δl) + cnt[k] == 0 && continue + sΔ += cnt[k] * δl[k] + sΔabs += abs(cnt[k] * δl[k]) + end + slope = tp * sΔ - sw + if slope >= 0 + # The slope sums the gradient of each label times its change, and that + # gradient is known only to within the roundoff of its terms; a slope + # within the sum of that roundoff does not resolve the direction. + slope <= BALANCED_NOISE * eps(T) * gscale && return :floor + return :stalled + end + α = α0 === nothing ? -slope / curv : T(α0) + α = min(α, αmax, BALANCED_MAXSTEP / (p * Δmax)) + for _ in 1:BALANCED_NBACKTRACK + ΔF = scale = sw1 = zero(T) + for e in eachindex(d) + active[e] && continue + h = α * Δd[e] + iszero(h) && continue + w = exp(-p * d[e]) + q = w * expm1(-p * h) / p + ΔF += q + scale += abs(q) + sw1 += w * exp(-p * h) * Δd[e] + end + ΔF += tp * α * sΔ + scale += tp * α * sΔabs + slope1 = tp * sΔ - sw1 + accepted = ΔF <= BALANCED_ARMIJO * α * slope || + # The objective cannot resolve this step: accept it if the + # directional derivative shrinks. + (abs(ΔF) <= BALANCED_NOISE * eps(T) * scale && abs(slope1) < abs(slope)) + if accepted + y .+= α .* δv + α == αmax || return 0 + k = 0 + for e in eachindex(d) + active[e] && continue + if Δd[e] < 0 && (d[e] / -Δd[e] == αmax || + y[P[e]] - y[Q[e]] - C[e] <= _balanced_tight(y, C, P, Q, e)) + active[e] = true + k += 1 + end + end + return k + end + α /= 2 + end + return :stalled +end diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 4d471e4..7b8c9dd 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -2385,7 +2385,8 @@ function symcover_min_jump end # JuMP reference used to test the native asymmetric solver. function cover_min_jump end -# JuMP reference used to test `cover_transversal`; also returns `log π*`. +# JuMP references used to test `cover_transversal`, with HiGHS for `AbsLog{2}` and +# Ipopt for `PowerMean`; each also returns `log π*`. function cover_transversal_jump end """ diff --git a/src/stats.jl b/src/stats.jl index 107f20f..76c6fab 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -241,6 +241,70 @@ 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) +""" + BalancedCoverStats{T} <: AbstractCoverStats + +Statistics of the transversal-tight covers with a `PowerMean` penalty and a +target ratio ([`cover_transversal(::PowerMean, A)`](@ref cover_transversal) and +its mutating form). + +The solver is an active-set Newton method on groups of rows joined by entries +at zero slack. The *relative reduced gradient* of a group is the derivative of +the objective with respect to its common log scale, divided by the sum of the +magnitudes of the terms of that derivative; `gradnorm` is its largest value +over the groups. The groups are stationary when `gradnorm <= tol`, or when +roundoff keeps the Newton step from reducing the objective: the step changes no +log slack `-log(r)` by more than `tol`, or the directional derivative of the +objective along it is nonnegative but within the roundoff of the groups' +derivatives. This happens when roundoff limits `gradnorm`, so `converged` can +hold with `gradnorm > tol`. The solve stops at a stationary point where every +group's multiplier test succeeds. +`converged`, `residual`, and `tolerance` return `converged`, `gradnorm`, and +`tol`; `iterations` returns `nnewton`. + +# Fields + +- `converged::Bool`: the stopping criterion holds. +- `stalled::Bool`: the line search could not reduce the objective (and so + `converged` is `false`); otherwise an unconverged solve reached `maxiter`. +- `tol::T`, `gradnorm::T`: as above. +- `step::T`: the largest change of a log slack in the last Newton step (`NaN` + if none was computed). +- `nnewton::Int`: Newton steps. +- `nactivations::Int`: entries that reached zero slack and joined the active + set. +- `nreleases::Int`: multiplier tests that released active entries. +- `nactive::Int`: entries at zero slack at the solution, not counting those + that lie on some maximum-product transversal. +- The remaining fields describe the linear algebra; their names and meanings + may change between releases: `nfactor` (factorizations) and `linsolve` (the + solver of the last Newton system: `:dense`, `:cholesky`, `:cg`, or `:none` + when no system was solved). +""" +struct BalancedCoverStats{T} <: AbstractCoverStats + converged::Bool + stalled::Bool + tol::T + gradnorm::T + step::T + nnewton::Int + nfactor::Int + nactivations::Int + nreleases::Int + nactive::Int + linsolve::Symbol +end + +converged(s::BalancedCoverStats) = s.converged +iterations(s::BalancedCoverStats) = s.nnewton +residual(s::BalancedCoverStats) = s.gradnorm +tolerance(s::BalancedCoverStats) = s.tol + +# Conversion from the statistics tuple of `_balanced_newton!`. +BalancedCoverStats{T}(nt::NamedTuple) where {T} = + BalancedCoverStats{T}(nt.converged, nt.stalled, nt.tol, nt.gradnorm, nt.step, nt.nnewton, nt.nfactor, + nt.nactivations, nt.nreleases, nt.nactive, nt.linsolve) + """ ExternalSolverStats <: AbstractCoverStats diff --git a/src/transversal.jl b/src/transversal.jl index c216e8b..caf2df2 100644 --- a/src/transversal.jl +++ b/src/transversal.jl @@ -53,6 +53,10 @@ grid, `:woodbury` when its requirements hold, and `:dense` otherwise. Sparse matrices default to `:lsqr`. As for [`cover_min`](@ref), an active-set method finishes a solve whose multiplier iteration stops early. +`cover_transversal(AbsLog{2}(), A)` is the same. With a [`PowerMean`](@ref) +penalty, [`cover_transversal(ϕ, A; target)`](@ref cover_transversal) instead +pushes the coverage ratios of the entries off the transversal toward `target`. + See also: [`cover_transversal!`](@ref), [`cover_min`](@ref). """ function cover_transversal(A::AbstractMatrix; κ::Real=AL_PENALTY, maxouter::Int=AL_MAXOUTER, @@ -73,8 +77,7 @@ and `b` and returns them with the solver's statistics. The initial contents of `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)))")) + _check_cover_axes(a, b, A) anew, bnew, stats = cover_transversal(A; kwargs...) a .= anew b .= bnew diff --git a/test/balanced.jl b/test/balanced.jl new file mode 100644 index 0000000..aaa22f0 --- /dev/null +++ b/test/balanced.jl @@ -0,0 +1,440 @@ +# Transversal-tight covers with a PowerMean penalty and a target ratio. +# Uses `max_product_transversals`, `logproduct_gap`, `_max_product_transversal`, and +# `_row_support` from transversal.jl. + +using MatrixCovers: _flow_deficit, _cover_transversal_balanced, _balanced_step! + +if !isdefined(@__MODULE__, :symmetric_matrices) + include("testmatrices.jl") +end + +# Coverage ratios of the entries of `A` off the transversal `σ` (by row), as +# `(i, j) => r` pairs. +function offtransversal_ratios(a, b, A, σ) + return [(i, j) => abs(A[i, j]) / (a[i] * b[j]) for i in axes(A, 1), j in axes(A, 2) + if !iszero(A[i, j]) && j != σ[i]] +end + +# The objective of `cover_transversal(PowerMean{p}(), A; target=t)` over the +# entries off `σ`. +function balanced_objective(a, b, A, σ, p, t) + return sum(((r / t)^p - 1 - p * log(r / t)) * t^p / p for (_, r) in offtransversal_ratios(a, b, A, σ); + init=0.0) +end + +# Largest geometric-mean coverage ratio of a cycle of entries off the unique +# maximum-product transversal: the least possible largest ratio of a cover tight +# on the transversal. Entry `(i, j)` off the transversal joins row `i` to the row +# `k` matched to column `j` with weight `log|A[i,j]| - log|A[k,j]|`; Karp's +# algorithm gives the largest mean weight of a cycle. +function minimax_ratio(A) + σs, _ = max_product_transversals(A) + length(σs) == 1 || error("the maximum-product transversal is not unique") + σ = only(σs) + n = size(A, 1) + τ = invperm(σ) + E = [(i, τ[j], log(abs(A[i, j])) - log(abs(A[τ[j], j]))) for i in 1:n, j in 1:n + if !iszero(A[i, j]) && j != σ[i]] + D = fill(-Inf, n + 1, n) # D[k+1, v]: heaviest walk of k edges ending at v + D[1, :] .= 0 + for k in 1:n, (u, v, c) in E + D[k+1, v] = max(D[k+1, v], D[k, u] + c) + end + λ = -Inf + for v in 1:n + isfinite(D[n+1, v]) || continue + λ = max(λ, minimum((D[n+1, v] - D[k+1, v]) / (n - k) for k in 0:n-1)) + end + return exp(λ), σ +end + +# Checks shared by every solve: a certified hard cover whose log product is the +# largest log product of a transversal, so that it is tight on every +# maximum-product transversal. For small `A`, the transversals are enumerated +# and the tightness is checked entry by entry. +function check_balanced_cover(a, b, A; atol=1e-9) + iscover(a, b, A) || return false + if size(A, 1) <= 7 + σs, lπ = max_product_transversals(A) + all(σ -> all(i -> abs(log(a[i] * b[σ[i]] / abs(A[i, σ[i]]))) <= atol, axes(A, 1)), σs) || + return false + else + _, _, _, lπ = _max_product_transversal(_row_support(A, Float64), axes(A, 2), :test) + end + return abs(sum(log, a) + sum(log, b) - lπ) <= atol * max(1, abs(lπ)) +end + +@testset "maximum-flow multiplier test" begin + # Opposite supplies on a 2-cycle: either orientation routes. + for g in ([1.0, -1.0], [-1.0, 1.0]) + deficit, _ = _flow_deficit(2, [1, 2], [2, 1], g) + @test deficit == 0 + end + # One arc: feasible along it, infeasible against it, with `S` the stranded supply. + deficit, _ = _flow_deficit(2, [1], [2], [1.0, -1.0]) + @test deficit == 0 + deficit, S = _flow_deficit(2, [1], [2], [-1.0, 1.0]) + @test deficit == 1 + @test S == [false, true] + # Redundant arcs forming an undirected cycle: two routes from 1 to 3. + deficit, _ = _flow_deficit(3, [1, 1, 2], [2, 3, 3], [2.0, 0.0, -2.0]) + @test deficit == 0 + deficit, S = _flow_deficit(3, [1, 1, 2], [2, 3, 3], [-2.0, 0.0, 2.0]) + @test deficit == 2 + @test S == [false, false, true] + # Partly routable: the chain 1 → 2 → 3 absorbs only 0.5 of the supply. + g = [1.0, 0.5, -0.5, -1.0] + deficit, S = _flow_deficit(4, [1, 2], [2, 3], g) + @test deficit ≈ 1 + @test S == [true, true, true, false] + @test deficit ≈ sum(g[S]) + # No arcs and no supply. + deficit, _ = _flow_deficit(2, Int[], Int[], [0.0, 0.0]) + @test deficit == 0 + @test_throws DimensionMismatch _flow_deficit(3, [1], [2], [1.0, -1.0]) + # On random digraphs, no arc leaves `S` and the deficit is the supply in `S`. + rng = StableRNG(4) + for _ in 1:50 + nv = rand(rng, 2:8) + na = rand(rng, 0:12) + tails, heads = rand(rng, 1:nv, na), rand(rng, 1:nv, na) + g = randn(rng, nv) + g .-= sum(g) / nv + deficit, S = _flow_deficit(nv, tails, heads, g) + @test deficit >= 0 + @test !any(k -> S[tails[k]] && !S[heads[k]], eachindex(tails)) + @test deficit ≈ max(sum(g[S]; init=0.0), 0) atol = 1e-12 + end +end + +@testset "line search with cancelling slack changes" begin + # Row 2 has one entry entering and one leaving, so moving it changes the sum + # of the slacks by exactly zero. The descent comes from the change `ε` of + # the slack of the first entry, which is far below the resolution of the + # change `δ` of row 2: summed by entry, `ε` is lost and the slope is + # positive. + p, tp = 2.0, 1e-4 + ε, δ = -1e-17, 0.25 + P, Q = [4, 1, 2], [3, 2, 3] + d = [30.0, -log(2e-21) / p, -log(1e-21) / p] + y, C = zeros(4), -d + δv = [0, δ, 0, ε] + status = _balanced_step!(y, δv, d, zeros(3), falses(3), 1:4, P, Q, C, p, tp, 1.0) + @test status == 0 + @test y == δv + @test_throws "not constant on a label" _balanced_step!(zeros(4), δv, d, zeros(3), falses(3), [1, 1, 2, 3], + P, Q, C, p, tp, 1.0) +end + +@testset "line search at the resolution of the arithmetic" begin + # Moving row 2, which has one entry entering and one leaving, is stationary + # to first order when their slacks are equal, and uphill when the entering + # one has the smaller slack. + p, tp = 2.0, 1e-4 + P, Q = [1, 2], [2, 3] + for (d, expected) in (([5.0, 5.0], :floor), ([1.0, 2.0], :stalled)) + y = zeros(3) + @test _balanced_step!(y, [0, 0.5, 0], d, zeros(2), falses(2), 1:3, P, Q, -d, p, tp, 1.0) === expected + @test iszero(y) + end + # An inactive entry at zero slack that the step would make negative is + # made active, and the scales do not move. + y, active = zeros(3), falses(2) + @test _balanced_step!(y, [0, 0.5, 0], [0.0, 1.0], zeros(2), active, 1:3, P, Q, [0.0, -1.0], p, tp, 1.0) == 1 + @test active == [true, false] && iszero(y) +end + +@testset "groups whose terms differ greatly in size" begin + # Entries spanning about 1e±11, in units spanning a further 2^±30: some + # groups' gradients have terms many orders of magnitude smaller than + # others', and the solve must still reach the tight tolerance. + rng = StableRNG(16) + A = sprandn(rng, 12, 12, 0.15) + A = A .* exp.(6 .* sprandn(rng, 12, 12, 1.0) .* (A .!= 0)) + A = dropzeros(A + spdiagm(0 => exp.(6 .* randn(rng, 12)))) + a, b, st = cover_transversal(PowerMean{2}(), A) + @test MatrixCovers.converged(st) + @test check_balanced_cover(a, b, A) + d, e = 2.0 .^ rand(rng, -30:30, 12), 2.0 .^ rand(rng, -30:30, 12) + a12, b12 = cover_transversal(PowerMean{2}(), A; tol=1e-12) + ã, b̃ = cover_transversal(PowerMean{2}(), Diagonal(d) * A * Diagonal(e); tol=1e-12) + @test logproduct_gap(ã ./ d, b̃ ./ e, a12, b12, A) <= 1e-9 +end + +@testset "cover_transversal(::PowerMean)" begin + # Agreement with the Ipopt reference on the committed matrices, whose many + # repeated magnitudes leave degenerate active sets. + idx = round.(Int, range(1, length(general_matrices), length=40)) + nreleases = 0 + for k in idx + Af = Float64.(general_matrices[k][2]) + for (p, t) in ((1, 1e-2), (2, 1e-2), (2, 0.3), (4, 0.3), (1, 1.0), (2, 1.0), (4, 1.0)) + a, b, st = cover_transversal(PowerMean{p}(), Af; target=t) + @test MatrixCovers.converged(st) + @test check_balanced_cover(a, b, Af) + (aj, bj), lπ = MatrixCovers.cover_transversal_jump(PowerMean{p}(), Af; target=t) + @test logproduct_gap(a, b, aj, bj, Af) <= 1e-6 + nreleases += st.nreleases + end + end + @test nreleases > 0 + + # Generic magnitudes with sparse patterns. + rng = StableRNG(3) + for _ in 1:30 + n = rand(rng, 3:7) + A = exp.(2 .* randn(rng, n, n)) .* (rand(rng, n, n) .< 0.6) + for i in 1:n + A[i, i] = exp(2 * randn(rng)) + 0.1 + end + for (p, t, rtol) in ((1, 1e-2, 1e-6), (2, 1e-2, 1e-6), (2, 0.3, 1e-6), (2, 1.0, 1e-5)) + a, b, st = cover_transversal(PowerMean{p}(), A; target=t) + @test MatrixCovers.converged(st) + @test check_balanced_cover(a, b, A) + (aj, bj), _ = MatrixCovers.cover_transversal_jump(PowerMean{p}(), A; target=t) + @test logproduct_gap(a, b, aj, bj, A) <= rtol + end + # For p = 4 and a small target the objective is too flat for Ipopt to + # resolve the products; compare with an extended-precision solve. + a, b, st = cover_transversal(PowerMean{4}(), A; target=1e-2) + @test MatrixCovers.converged(st) + ab, bb, sb = _cover_transversal_balanced(BigFloat.(A), 4, BigFloat(1e-2); tol=1e-30) + @test sb.converged + @test logproduct_gap(a, b, ab, bb, A) <= 1e-8 + end + + # KKT conditions across a cut: with every entry of `A12` covered loosely, + # shifting the second block's scales is unconstrained, so the p-power mean + # of the ratios of the entries crossing the cut equals the target. + A11 = [4.0 1 0.5; 1 4 1; 0.5 1 4] + A22 = [3.0 0.5 1; 1 5 0.5; 2 1 3] + for (p, t, A12) in ((2, 1e-2, [1e-3 5.0 0; 2.0 0 1e3; 0 30.0 0.1]), + (1, 0.1, [1.0 0 0; 0 0 1e-4; 2.0 0 0]), + (3, 0.5, [7.0 1 1; 1 1 1; 1 1 0.01])) + A = [A11 A12; zeros(3, 3) A22] + a, b, st = cover_transversal(PowerMean{p}(), A; target=t) + @test MatrixCovers.converged(st) + rcut = [A[i, j] / (a[i] * b[j]) for i in 1:3, j in 4:6 if A[i, j] != 0] + @test maximum(rcut) < 1 + @test sum(r -> r^p, rcut) / length(rcut) ≈ t^p rtol = 1e-8 + end + + # A group with an inactive entry inside it: rows 4:end are tied to row 3 + # (all of the transversals through them are maximal), and the entries + # `A[1, 4:end]`, far below their targets, pull row 1 down against row 3 + # until the cycle 1 → 2 → 3 → 1 has two tight entries; its third entry keeps + # the cycle's slack. + function pulled_cycle(npull) + n = 3 + npull + A = Matrix{Float64}(I, n, n) + A[1, 2] = A[2, 3] = A[3, 1] = 0.5 + for k in 4:n + A[3, k] = A[k, 3] = 1.0 + A[1, k] = 1e-3 + end + return A + end + for (npull, p, t) in ((2, 1, 0.5), (4, 2, 0.5)) + A = pulled_cycle(npull) + a, b, st = cover_transversal(PowerMean{p}(), A; target=t) + @test MatrixCovers.converged(st) + R = A ./ (a .* b') + @test R[1, 2] ≈ 1 && R[2, 3] ≈ 1 + @test R[3, 1] ≈ 0.125 + @test st.nactive == 2 + (aj, bj), _ = MatrixCovers.cover_transversal_jump(PowerMean{p}(), A; target=t) + @test logproduct_gap(a, b, aj, bj, A) <= 1e-6 + end + # With fewer pulling entries, no entry of the cycle is tight. + A = pulled_cycle(1) + a, b, st = cover_transversal(PowerMean{1}(), A; target=0.5) + @test st.nactive == 0 + @test maximum(A[i, j] / (a[i] * b[j]) for (i, j) in ((1, 2), (2, 3), (3, 1))) < 0.99 + + # 0/1 matrices: every off-transversal ratio is tied at the start. + rng = StableRNG(5) + for _ in 1:30 + n = rand(rng, 3:7) + A = Float64.(rand(rng, n, n) .< 0.5) + for i in 1:n + A[i, i] = 1 + end + for (p, t, rtol) in ((1, 0.3, 1e-6), (2, 0.3, 1e-6), (2, 1.0, 1e-5)) + a, b, st = cover_transversal(PowerMean{p}(), A; target=t) + @test MatrixCovers.converged(st) + @test check_balanced_cover(a, b, A) + (aj, bj), _ = MatrixCovers.cover_transversal_jump(PowerMean{p}(), A; target=t) + @test logproduct_gap(a, b, aj, bj, A) <= rtol + end + a, b, st = cover_transversal(PowerMean{2}(), A; target=1e-2) + @test MatrixCovers.converged(st) + @test check_balanced_cover(a, b, A) + end + # Every transversal is maximal: one tight component and nothing to solve. + a, b, st = cover_transversal(PowerMean{2}(), ones(4, 4)) + @test a .* b' ≈ ones(4, 4) + @test MatrixCovers.converged(st) && st.nnewton == 0 + + # Large p approaches the cover with the least largest ratio off the + # transversal, which equals the largest geometric-mean ratio of a cycle. + rng = StableRNG(21) + for _ in 1:8 + A = exp.(randn(rng, 5, 5)) + θ, σ = minimax_ratio(A) + a, b, st = cover_transversal(PowerMean{16}(), A; target=1e-2) + @test MatrixCovers.converged(st) + rmax = maximum(last, offtransversal_ratios(a, b, A, σ)) + @test θ * (1 - 1e-12) <= rmax <= 1.02 * θ + end + + # A 2×2 triangular matrix: the off-diagonal ratio is the target, so the + # scaled matrix approaches the identity as the target shrinks. + A = [1.0 1e6; 0 1.0] + for t in (1e-1, 1e-2, 1e-4) + a, b, st = cover_transversal(PowerMean{2}(), A; target=t) + @test A[1, 2] / (a[1] * b[2]) ≈ t + @test cond(A ./ (a .* b')) <= 1 + 2t + end + + # The minimizer depends on `target`, and the objective is minimal among + # nearby covers tight on the transversal. + A = Float64.(general_matrices[17][2]) + σ = first(max_product_transversals(A)[1]) + M = exp.(randn(StableRNG(8), 5, 5)) .* (rand(StableRNG(9), 5, 5) .< 0.5) + 4I + a1, b1 = cover_transversal(PowerMean{2}(), M; target=0.5) + a2, b2 = cover_transversal(PowerMean{2}(), M; target=0.05) + @test logproduct_gap(a1, b1, a2, b2, M) > 1e-3 + a, b = cover_transversal(PowerMean{2}(), A; target=0.1) + F = balanced_objective(a, b, A, σ, 2, 0.1) + rng = StableRNG(9) + for _ in 1:20 + # Perturb the row scales, keep the transversal tight, and restore coverage. + x = 1e-3 .* randn(rng, 5) + ã = a .* exp.(x) + b̃ = copy(b) + for i in 1:5 + b̃[σ[i]] = b[σ[i]] * exp(-x[i]) + end + iscover(ã, b̃, A) || continue + @test balanced_objective(ã, b̃, A, σ, 2, 0.1) >= F - 1e-12 + end + + # Independence from the maximum-product transversal used. + ntie = 0 + for (_, A) in general_matrices + ntie >= 10 && break + Af = Float64.(A) + σs, _ = max_product_transversals(Af) + length(σs) > 1 || continue + ntie += 1 + a, b = cover_transversal(PowerMean{2}(), Af) + for σ in σs + aσ, bσ, _ = _cover_transversal_balanced(Af, 2, 1e-2; transversal=σ) + @test logproduct_gap(a, b, aσ, bσ, Af) <= 1e-8 + end + end + @test ntie == 10 + + # Scale covariance and permutation equivariance. + rng = StableRNG(1) + for (_, A) in general_matrices[1:20] + Af = Float64.(A) + dr = exp.(2 .* randn(rng, 5)); dc = exp.(2 .* randn(rng, 5)) + @test covaries(M -> scales(cover_transversal(PowerMean{2}(), M)), Af, dr, dc; rtol=1e-7) + pr, pc = randperm(rng, 5), randperm(rng, 5) + a, b = cover_transversal(PowerMean{2}(), Af) + ap, bp = cover_transversal(PowerMean{2}(), Af[pr, pc]) + @test (ap .* bp') ≈ (a .* b')[pr, pc] rtol = 1e-8 + end + + # Sparse matrices: power-of-two rescaling, sparse and dense storage, and the + # sparse and conjugate-gradient Newton solves. + for seed in 1:3 + rng = StableRNG(seed) + n = 200 + A = sprandn(rng, n, n, 0.02) + spdiagm(0 => exp.(3 .* randn(rng, n))) + a, b, st = cover_transversal(PowerMean{2}(), A) + @test MatrixCovers.converged(st) + @test st.linsolve === :cholesky + @test check_balanced_cover(a, b, A; atol=1e-8) + @test isbalanced(a, b, A) + d = 2.0 .^ rand(rng, -30:30, n) + e = 2.0 .^ rand(rng, -30:30, n) + ã, b̃ = scales(cover_transversal(PowerMean{2}(), Diagonal(d) * A * Diagonal(e))) + @test logproduct_gap(ã ./ d, b̃ ./ e, a, b, A) <= 1e-8 + ad, bd, sd = cover_transversal(PowerMean{2}(), Matrix(A)) + @test logproduct_gap(a, b, ad, bd, A) <= 1e-8 + ac, bc, sc = cover_transversal(PowerMean{2}(), A; flopbudget=0) + @test sc.linsolve === :cg + @test MatrixCovers.converged(sc) + @test logproduct_gap(a, b, ac, bc, A) <= 1e-8 + end + + # 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] + a, b = cover_transversal(PowerMean{2}(), Ablk) + @test check_balanced_cover(a, b, Ablk) + @test isbalanced(a, b, Ablk) + + # The AbsLog{2} penalty is the default solver. + A = [1.0 1 0; 0.25 1 1; 1 0 0] + @test scales(cover_transversal(AbsLog{2}(), A)) == scales(cover_transversal(A)) + a2, b2 = zeros(3), zeros(3) + @test scales(cover_transversal!(AbsLog{2}(), a2, b2, A)) === (a2, b2) + @test (a2, b2) == scales(cover_transversal(A)) + + # Generic indexing, element types, and the mutating form. + A = [1.0 0.5 0.2; 0.3 1.0 0.5; 0.4 0.1 1.0] + a, b, st = cover_transversal(PowerMean{2}(), A) + @test st isa MatrixCovers.BalancedCoverStats{Float64} + @test MatrixCovers.iterations(st) == st.nnewton + @test MatrixCovers.residual(st) == st.gradnorm + @test MatrixCovers.tolerance(st) == st.tol == 1e-9 + Ao = OffsetArray(A, -1:1, 5:7) + ao, bo = cover_transversal(PowerMean{2}(), Ao) + @test axes(ao, 1) == axes(Ao, 1) && axes(bo, 1) == axes(Ao, 2) + @test parent(ao) .* parent(bo)' ≈ a .* b' + B = [9.0 9 9 9; 9 1 0.5 0.2; 9 0.3 1 0.5; 9 0.4 0.1 1] + av, bv = cover_transversal(PowerMean{2}(), view(B, 2:4, 2:4)) + @test av .* bv' ≈ a .* b' + ao, bo = cover_transversal!(PowerMean{2}(), OffsetArray(zeros(3), -1:1), OffsetArray(zeros(3), 5:7), Ao) + @test parent(ao) .* parent(bo)' ≈ a .* b' + a32, b32, st32 = cover_transversal(PowerMean{2}(), Float32.(A)) + @test eltype(a32) == eltype(b32) == Float32 + @test st32 isa MatrixCovers.BalancedCoverStats{Float64} + @test iscover(a32, b32, Float32.(A)) + @test a32 .* b32' ≈ a .* b' rtol = 1e-6 + aC, bC = cover_transversal(PowerMean{2}(), A .* cis.(reshape(1:9, 3, 3))) + @test aC .* bC' ≈ a .* b' + aB, bB, stB = cover_transversal(PowerMean{2}(), BigFloat.(A)) + @test eltype(aB) == BigFloat + @test stB isa MatrixCovers.BalancedCoverStats{BigFloat} + @test aB .* bB' ≈ a .* b' rtol = 1e-8 + a2, b2 = zeros(3), zeros(3) + @test scales(cover_transversal!(PowerMean{2}(), a2, b2, A)) === (a2, b2) + @test (a2, b2) == (a, b) + @test_throws DimensionMismatch cover_transversal!(PowerMean{2}(), zeros(2), zeros(3), A) + + # An iteration limit: a warning, and still a hard cover tight on the transversal. + M = exp.(randn(StableRNG(2), 6, 6)) + (a, b, st), logs = capture_logs(() -> cover_transversal(PowerMean{2}(), M; maxiter=1)) + @test !MatrixCovers.converged(st) + @test !st.stalled + @test any(l -> l.level == Base.CoreLogging.Warn && occursin("maxiter=1", l.message), logs) + @test check_balanced_cover(a, b, M) + a, b, st = @test_logs (:warn, r"may not minimize the objective") cover_transversal(PowerMean{2}(), M; maxiter=0) + @test check_balanced_cover(a, b, M) + @test isbalanced(a, b, M) + + # Errors. + @test_throws "requires a square matrix" cover_transversal(PowerMean{2}(), [1.0 2.0 3.0; 4.0 5.0 6.0]) + @test_throws "structurally nonsingular" cover_transversal(PowerMean{2}(), [1.0 1 1; 1 0 0; 1 0 0]) + @test_throws "structurally nonsingular" cover_transversal(PowerMean{2}(), sparse([1.0 1 1; 1 0 0; 1 0 0])) + @test_throws "finite entries" cover_transversal(PowerMean{2}(), [1.0 Inf; 1 1]) + @test_throws "target must satisfy 0 < target <= 1, got 0" cover_transversal(PowerMean{2}(), A; target=0) + @test_throws "target must satisfy 0 < target <= 1, got 1.5" cover_transversal(PowerMean{2}(), A; target=1.5) + @test_throws "requires a finite exponent" cover_transversal(PowerMean{Inf}(), A) + @test_throws "underflows" cover_transversal(PowerMean{2}(), A; target=1e-200) + @test_throws "tol must be positive" cover_transversal(PowerMean{2}(), A; tol=0) + @test_throws "maxiter must be nonnegative" cover_transversal(PowerMean{2}(), A; maxiter=-1) +end diff --git a/test/runtests.jl b/test/runtests.jl index 96e6e7e..8cc9b98 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,6 +30,7 @@ include("helpers.jl") # isbalanced, covaries, PENALTIES include("tree_cholesky.jl") # the fill-free forest factorization behind the LSQR preconditioner include("minimal_covers.jl") # the *_min family (native solvers) include("transversal.jl") # cover_transversal and the maximum-product transversal + include("balanced.jl") # cover_transversal with a PowerMean penalty and a target include("active_set.jl") # exact solution of stalled AbsLog{2} minimal covers include("storage_types.jl") # sparse/structured/wrapped storage vs dense reference include("element_types.jl") # Float32/BigFloat: eltype-scaled internal tolerances @@ -49,7 +50,8 @@ include("helpers.jl") # isbalanced, covaries, PENALTIES :require_abs_symmetric, :_edge_list, :_sym_edge_list, :_degrees, :_support_components, :_balance_cover!, :inflate_feasible!, - :_balance_bipartite_sym!, :_sym_support) + :_balance_bipartite_sym!, :_sym_support, + :_max_product_transversal, :_row_support) # External non-public names with no usable public equivalent. foreign = (:FreeUnits, :Unit, :Units, :Optimizer, :Experimental, :register_error_hint, :RefValue, :getcolptr, :LibSuiteSparse, :TRUE, diff --git a/test/unitful.jl b/test/unitful.jl index 33d7ebb..a7f1f40 100644 --- a/test/unitful.jl +++ b/test/unitful.jl @@ -106,6 +106,8 @@ @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))) + @test all(unit.(v) == fill(u"m^-1", 3) for v in scales(cover_transversal(PowerMean{2}(), S))) + @test all(unit.(v) == fill(u"m^-1", 3) for v in scales(cover_transversal!(PowerMean{2}(), a, b, S))) end @testset "balance convention holds in the caller's units" begin @@ -183,6 +185,7 @@ for ϕ in (PowerMean{1}(), PowerMean{2}()) @test unit.(scales(soft_symcover(ϕ, A))) == UA @test all(unit.(v) == UA for v in scales(soft_cover(ϕ, A))) + @test all(unit.(v) == UA for v in scales(cover_transversal(ϕ, A))) end end From 53774bf6a7ad73eac82bd0b19b66ae5a9400ff90 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 23 Oct 2026 12:59:03 -0500 Subject: [PATCH 3/5] Factor Tarjan's algorithm out of _tight_components `_strong_components(ptr, heads)` labels the strongly connected components of a digraph given in compressed adjacency form, numbered in reverse topological order. `_tight_components` calls it on the zero-residual entries. Assisted-by: Claude Opus 5.5 --- src/transversal.jl | 97 ++++++++++++++++++++++++++-------------------- 1 file changed, 54 insertions(+), 43 deletions(-) diff --git a/src/transversal.jl b/src/transversal.jl index caf2df2..e337ce7 100644 --- a/src/transversal.jl +++ b/src/transversal.jl @@ -406,49 +406,8 @@ function _tight_components(edges::Vector{Tuple{Int,Int}}, cvals::Vector{T}, α0: outedge[cursor[p]] = e cursor[p] += 1 end - # Iterative Tarjan. - comp = zeros(Int, m) - index = zeros(Int, m) - low = zeros(Int, m) - onstack = falses(m) - stack = Int[] - callstack = Tuple{Int,Int}[] # (row, next slot in its adjacency list) - counter = 0 - nc = 0 - for r in 1:m - index[r] == 0 || continue - push!(callstack, (r, outptr[r])) - counter += 1 - index[r] = low[r] = counter - push!(stack, r); onstack[r] = true - while !isempty(callstack) - p, s = callstack[end] - if s < outptr[p+1] - callstack[end] = (p, s + 1) - q = edges[outedge[s]][2] - if index[q] == 0 - counter += 1 - index[q] = low[q] = counter - push!(stack, q); onstack[q] = true - push!(callstack, (q, outptr[q])) - elseif onstack[q] - low[p] = min(low[p], index[q]) - end - else - pop!(callstack) - if low[p] == index[p] - nc += 1 - while true - q = pop!(stack) - onstack[q] = false - comp[q] = nc - q == p && break - end - end - isempty(callstack) || (low[callstack[end][1]] = min(low[callstack[end][1]], low[p])) - end - end - end + heads = [edges[e][2] for e in outedge] + comp, nc = _strong_components(outptr, heads) # Offsets from a breadth-first spanning tree of each component, crossing # tight entries within the component in either direction. inptr = zeros(Int, m + 1) @@ -494,5 +453,57 @@ function _tight_components(edges::Vector{Tuple{Int,Int}}, cvals::Vector{T}, α0: return comp, o, nc end +# Strongly connected components of the digraph on `1:nv`, `nv = length(ptr) - 1`, +# with arcs `v → heads[s]` for `s in ptr[v]:ptr[v+1]-1` (iterative Tarjan). +# Returns `(comp, nc)`: vertex `v` belongs to component `comp[v]` in `1:nc`. The +# components are numbered in reverse topological order: every arc between +# components goes from a higher label to a lower one. +function _strong_components(ptr::Vector{Int}, heads::Vector{Int}) + m = length(ptr) - 1 + comp = zeros(Int, m) + index = zeros(Int, m) + low = zeros(Int, m) + onstack = falses(m) + stack = Int[] + callstack = Tuple{Int,Int}[] # (vertex, next slot in its adjacency list) + counter = 0 + nc = 0 + for r in 1:m + index[r] == 0 || continue + push!(callstack, (r, ptr[r])) + counter += 1 + index[r] = low[r] = counter + push!(stack, r); onstack[r] = true + while !isempty(callstack) + p, s = callstack[end] + if s < ptr[p+1] + callstack[end] = (p, s + 1) + q = heads[s] + if index[q] == 0 + counter += 1 + index[q] = low[q] = counter + push!(stack, q); onstack[q] = true + push!(callstack, (q, ptr[q])) + elseif onstack[q] + low[p] = min(low[p], index[q]) + end + else + pop!(callstack) + if low[p] == index[p] + nc += 1 + while true + q = pop!(stack) + onstack[q] = false + comp[q] = nc + q == p && break + end + end + isempty(callstack) || (low[callstack[end][1]] = min(low[callstack[end][1]], low[p])) + end + end + end + return comp, nc +end + # Roundoff allowance, in units of `eps`, for a zero residual at the transversal duals. const TIGHT_ULPS = 64 From de549f9eaebd2b3f94bfed377f8a4acb6ffb5645 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 24 Oct 2026 04:38:57 -0500 Subject: [PATCH 4/5] Add maxbalance keyword to cover_transversal MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `cover_transversal(ϕ::PowerMean{p}, A; target, maxbalance=true)` and its mutating form max-balance the scales within each fully indecomposable block of `A` (Hook, Pestana, Tisseur & Hogg 2019, Algorithm 3.10: max-plus policy iteration for a cycle of largest mean, a shift by its bias, contraction of the cycle) and then minimize the PowerMean objective with a target ratio over the entries joining different blocks, by the active-set Newton method on one scale per block. The products are unique and scale-covariant. A cycle mean accumulates the roundoff of up to `nv` weights, so the policy iteration treats cycle means within `nv * tol` of each other as equal; with `tol` as the threshold the comparisons are not transitive and the iteration can cycle among tied cycles. `BalancedCoverStats` records `maxbalanced` and the number of fully indecomposable blocks, `nblocks`, for every solve. test/maxbalance.jl covers Example 3.11 of Hook et al. (2019), policy iteration on nearly equal cycle means, the cut balance by brute force and the minimax ratio on random fully indecomposable matrices, agreement with an independent implementation of their Algorithm 3.10, block-triangular matrices (the KKT identity across each cut, an active entry between blocks, covariance, the products of each block alone), tied transversals, sparse and wrapped inputs, element types, and the mutating form. Assisted-by: Claude Opus 5.5 --- docs/src/objectives.md | 12 +- src/MatrixCovers.jl | 1 + src/balanced.jl | 55 ++++++- src/maxbalance.jl | 250 +++++++++++++++++++++++++++++++ src/stats.jl | 9 +- test/maxbalance.jl | 323 +++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 7 files changed, 642 insertions(+), 9 deletions(-) create mode 100644 src/maxbalance.jl create mode 100644 test/maxbalance.jl diff --git a/docs/src/objectives.md b/docs/src/objectives.md index e1a3c14..c945fcf 100644 --- a/docs/src/objectives.md +++ b/docs/src/objectives.md @@ -58,7 +58,7 @@ 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}` (or `PowerMean{p}` with a target ratio) | — | +| [`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 | @@ -69,7 +69,11 @@ covers that are tight on a maximum-product transversal (see [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. A soft cover has no constraint to +`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. @@ -83,7 +87,9 @@ The method used to minimize the objective depends on the penalty: 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. + 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. diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index 431e2b1..42840dd 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -38,6 +38,7 @@ 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 diff --git a/src/balanced.jl b/src/balanced.jl index 63ac4f6..0820693 100644 --- a/src/balanced.jl +++ b/src/balanced.jl @@ -23,7 +23,7 @@ const BALANCED_NRIDGE = 4 const BALANCED_RIDGE_GROWTH = 1e4 """ - a, b, stats = cover_transversal(ϕ::PowerMean{p}, A; target=1e-2, tol=1e-9, maxiter, fillbudget, flopbudget) + a, b, stats = cover_transversal(ϕ::PowerMean{p}, A; target=1e-2, maxbalance=false, tol=1e-9, maxiter, fillbudget, flopbudget) Return the hard cover of the square matrix `A` that is tight on its maximum-product transversals and pushes the coverage ratios @@ -53,6 +53,26 @@ mean of the ratios of the entries crossing such a cut equals `target` whenever none of them has ratio one. Large `p` approaches the cover that minimizes the largest ratio within each block. +With `maxbalance=true`, the criterion within each fully indecomposable block of +`A` is the `p → ∞` limit of the objective (Hook, Pestana, Tisseur & Hogg, SIAM +J. Matrix Anal. Appl. 40:320–346, 2019, Theorem 3.2): the scales are +*max-balanced* within each block. Let `σ` be a maximum-product transversal, +and let an entry `A[i,j]` off it lead from row `i` to the row `k` with +`σ[k] == j`. For every partition of the rows of a block into two nonempty sets +`I` and `K`, the largest ratio of an entry leading from `I` to `K` equals the +largest ratio of an entry leading from `K` to `I`. This determines the products +`a[i]*b[j]` within each block uniquely. When `σ` is the only maximum-product +transversal, the largest ratio within a block equals the block's `θ`, the +largest geometric-mean ratio of a cycle of entries off `σ`, which is the least +largest ratio that any cover tight on `σ` can have there. `ϕ` and `target` fix +the scales between blocks: they minimize the objective above summed over the +entries joining different blocks, with each block's scales moving together. At +the minimizer, for every set of blocks that the entries joining it to the other +blocks all lead out of, the `p`-power mean of the ratios of those entries +equals `target` whenever none of them has ratio one. +The result is a hard cover tight on the maximum-product transversals, and its +products are unique on the support and scale-covariant. + `stats` is a [`BalancedCoverStats`](@ref MatrixCovers.BalancedCoverStats). The solver warns when [`MatrixCovers.converged(stats)`](@ref MatrixCovers.converged) is `false`; the result is then still a hard cover tight on the maximum-product @@ -73,6 +93,13 @@ descent. Where the groups are stationary, a maximum-flow problem on each group's active entries either proves the KKT conditions or yields a set of rows whose scales can decrease, releasing the active entries that enter it. +With `maxbalance=true`, each block is max-balanced by Algorithm 3.10 of Hook et +al.: max-plus policy iteration finds a cycle of largest geometric-mean ratio in +the block, a diagonal similarity brings every ratio of the block down to that +mean, and the cycle's rows, whose relative scales are then fixed, are merged +into one; this repeats until one row remains. The active-set Newton method then +runs with one scale per block. + `tol` (default `1e-9`) is the stopping tolerance: the groups are stationary when the derivative of the objective with respect to each group's scale is at most `tol` times the sum of the magnitudes of its terms, or when roundoff keeps the @@ -95,10 +122,11 @@ afterwards. See also: [`cover_transversal!`](@ref), [`PowerMean`](@ref). """ function cover_transversal(ϕ::PowerMean{p}, A::AbstractMatrix; target::Real=1e-2, + maxbalance::Bool=false, tol::Real=BALANCED_TOL, maxiter::Union{Integer,Nothing}=nothing, fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET) where {p} - a, b, nt = _cover_transversal_balanced(A, p, target; tol, maxiter, fillbudget, flopbudget) + a, b, nt = _cover_transversal_balanced(A, p, target; maxbalance, tol, maxiter, fillbudget, flopbudget) return a, b, BalancedCoverStats{_auglag_stats_type(A)}(nt) end @@ -126,7 +154,8 @@ cover_transversal!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::Abstrac # Worker for `cover_transversal(::PowerMean{p}, A; target)`, returning `(a, b, nt)` # with the statistics as a named tuple. `transversal` is as for # `_transversal_system`. -function _cover_transversal_balanced(A::AbstractMatrix, p::Real, target::Real; tol::Real=BALANCED_TOL, +function _cover_transversal_balanced(A::AbstractMatrix, p::Real, target::Real; maxbalance::Bool=false, + tol::Real=BALANCED_TOL, maxiter::Union{Integer,Nothing}=nothing, fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET, transversal::Union{Bool,AbstractVector{<:Integer}}=true, @@ -140,7 +169,7 @@ function _cover_transversal_balanced(A::AbstractMatrix, p::Real, target::Real; t # The tolerances require at least Float64 resolution. if eps(T) > eps(Float64) a64, b64, nt = _cover_transversal_balanced(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A), - p, target; tol, maxiter, fillbudget, flopbudget, + p, target; maxbalance, tol, maxiter, fillbudget, flopbudget, transversal, fname) # Narrowing rounds to nearest and so can round a product below its entry. a, b = T.(a64), T.(b64) @@ -162,7 +191,23 @@ function _cover_transversal_balanced(A::AbstractMatrix, p::Real, target::Real; t ts = _transversal_system(G, axc, T, fname, transversal) sys = _contract_system(ts.dsupp, ts.comp, ts.offs) maxit = maxiter === nothing ? 100 + 10 * ts.nc : Int(maxiter) - y, nt = _balanced_newton!(copy(ts.y0), sys, pT, tp; tol=T(tol), maxiter=maxit, fillbudget, flopbudget) + # The strongly connected components of the entries off the transversal are + # the fully indecomposable blocks of `A`. + ptr, arcs = _out_arcs(ts.nc, first.(sys.edges), _ -> true) + scc, nblocks = _strong_components(ptr, [sys.edges[e][2] for e in arcs]) + if maxbalance + # Each block's scales are fixed relative to one another by + # max-balancing; the blocks' common scales minimize the objective over + # the entries between blocks. + yb = _maxbalance_components(sys, ts.y0, scc, nblocks) + bsys = _contract_system(sys, scc, yb) + Y, nt = _balanced_newton!(_feasible_component_scales(bsys, nblocks), bsys, pT, tp; + tol=T(tol), maxiter=maxit, fillbudget, flopbudget) + y = Y[scc] .+ yb + else + y, nt = _balanced_newton!(copy(ts.y0), sys, pT, tp; tol=T(tol), maxiter=maxit, fillbudget, flopbudget) + end + nt = (; nt..., maxbalanced=maxbalance, nblocks) _warn_balanced_unconverged(fname, nt, maxit) # The uniform shift of `_transversal_scales` absorbs roundoff violations. x = _transversal_scales(y[ts.comp] .+ ts.offs, ts.dsupp, ts.σ, ts.ctrans, true) diff --git a/src/maxbalance.jl b/src/maxbalance.jl new file mode 100644 index 0000000..8aee864 --- /dev/null +++ b/src/maxbalance.jl @@ -0,0 +1,250 @@ +# Max-balancing of the transversal-tight difference system within its strongly +# connected components (Hook, Pestana, Tisseur & Hogg, SIAM J. Matrix Anal. +# Appl. 40:320–346, 2019, Algorithm 3.10). +# +# In the difference system on the component scales `y`, entry `e` joins +# `P[e] → Q[e]` with log slack `d[e] = y[P[e]] - y[Q[e]] - C[e] >= 0`, and so log +# coverage ratio `w[e] = -d[e] <= 0`. A shift `x` of the scales changes `w[e]` to +# `w[e] - x[P[e]] + x[Q[e]]`, which leaves the sum of `w` around every cycle +# unchanged. On a strongly connected digraph, the max-balanced shift makes the +# largest `w` leaving every set of vertices equal to the largest `w` entering it. +# It is unique up to a constant, and every `w` is then at most the largest mean +# weight of a cycle. + +# Arcs `e` with `keep(e)` grouped by tail `P[e]` among the vertices `1:nv`: +# `adj[ptr[v]:ptr[v+1]-1]` are the arcs leaving `v`. +function _out_arcs(nv::Int, P::AbstractVector{Int}, keep) + ptr = zeros(Int, nv + 1) + for e in eachindex(P) + keep(e) && (ptr[P[e]+1] += 1) + end + ptr[1] = 1 + cumsum!(ptr, ptr) + adj = Vector{Int}(undef, ptr[end] - 1) + cursor = ptr[1:nv] + for e in eachindex(P) + keep(e) || continue + adj[cursor[P[e]]] = e + cursor[P[e]] += 1 + end + return ptr, adj +end + +# Max-plus policy iteration (Cochet-Terrasson, Cohen, Gaubert, McGettrick & +# Quadrat, 1998) on the strongly connected digraph on `1:nv` whose arcs leaving +# `v` are `adj[ptr[v]:ptr[v+1]-1]`, arc `e` going to `dst[e]` with weight `w[e]`. +# `pol[v]` is the arc chosen at `v`, or 0 for none; it is the starting policy and +# is updated in place. `tol` is the resolution of the weights; a cycle mean +# accumulates the roundoff of up to `nv` of them, so cycle means within +# `nv * tol` of each other are treated as equal. +# +# Returns `(λ, s, cycle)`: the largest mean weight `λ` of a cycle, a vector `s` +# with `w[e] - s[v] + s[dst[e]] <= λ` on every arc `e` leaving `v`, up to +# roundoff, and the vertices of a cycle of mean `λ`, in order. +function _howard!(pol::Vector{Int}, nv::Int, ptr::Vector{Int}, adj::Vector{Int}, + dst::Vector{Int}, w::Vector{T}, tol::T) where {T} + for v in 1:nv + ptr[v] < ptr[v+1] || error("internal error: a vertex has no arc leaving it") + pol[v] == 0 || continue + best = adj[ptr[v]] + for s in ptr[v]:ptr[v+1]-1 + w[adj[s]] > w[best] && (best = adj[s]) + end + pol[v] = best + end + ηtol = nv * tol + η = zeros(T, nv) # mean weight of the cycle each vertex's policy path reaches + s = zeros(T, nv) # bias, relative to that cycle + state = zeros(Int, nv) + path = Int[] + cyc = Int[] + best = Int[] + for _ in 1:(100 + 10 * nv) + # Value determination on the functional graph `v → dst[pol[v]]`. + fill!(state, 0) # 0 unvisited, 1 on the current path, 2 done + λ = T(-Inf) + for r in 1:nv + state[r] == 0 || continue + empty!(path) + v = r + while state[v] == 0 + state[v] = 1 + push!(path, v) + v = dst[pol[v]] + end + if state[v] == 1 + # A new cycle through `v`: its mean fixes `η`, and `s[v]` + # keeps its value. + empty!(cyc) + u = v + μ = zero(T) + while true + push!(cyc, u) + μ += w[pol[u]] + u = dst[pol[u]] + u == v && break + end + μ /= length(cyc) + for u in cyc + η[u] = μ + end + for k in length(cyc):-1:2 + u = cyc[k] + s[u] = w[pol[u]] - μ + s[dst[pol[u]]] + end + for u in cyc + state[u] = 2 + end + if μ > λ + λ = μ + copy!(best, cyc) + end + end + for u in Iterators.reverse(path) + state[u] == 2 && continue + t = dst[pol[u]] + η[u] = η[t] + s[u] = w[pol[u]] - η[t] + s[t] + state[u] = 2 + end + end + # Policy improvement: first the cycle means, then the bias. + changed = false + for v in 1:nv + ηv = η[dst[pol[v]]] + e = pol[v] + for k in ptr[v]:ptr[v+1]-1 + a = adj[k] + η[dst[a]] > ηv + ηtol && (ηv = η[dst[a]]; e = a) + end + e == pol[v] || (pol[v] = e; changed = true) + end + if !changed + for v in 1:nv + cur = w[pol[v]] + s[dst[pol[v]]] + e = pol[v] + for k in ptr[v]:ptr[v+1]-1 + a = adj[k] + abs(η[dst[a]] - η[v]) <= ηtol || continue + val = w[a] + s[dst[a]] + val > cur + tol + TIGHT_ULPS * eps(T) * abs(cur) && (cur = val; e = a) + end + e == pol[v] || (pol[v] = e; changed = true) + end + end + if !changed + ηmin, ηmax = extrema(η) + ηmax - ηmin <= ηtol || + error("internal error: the cycle means differ ($ηmin to $ηmax) on a strongly connected digraph") + return λ, s, best + end + end + error("internal error: policy iteration did not converge") +end + +# Max-balancing shift `x` of the strongly connected digraph on `1:nv` with arcs +# `src[e] → dst[e]` of weight `w[e]`: repeatedly find a cycle of largest mean in +# the contracted digraph, shift by the bias of `_howard!`, which lowers every +# arc to at most that mean, and contract the cycle, whose relative shifts are +# then fixed. Returns `x`, with `x[1] == 0`. +function _maxbalance(nv::Int, src::Vector{Int}, dst::Vector{Int}, w::Vector{T}) where {T} + x = zeros(T, nv) + nv == 1 && return x + tol = TIGHT_ULPS * eps(T) * max(one(T), maximum(abs, w)) + parent = collect(1:nv) + function find(v) + while parent[v] != v + parent[v] = parent[parent[v]] + v = parent[v] + end + return v + end + pos = zeros(Int, nv) # position of each representative in the contracted digraph + pol = zeros(Int, nv) # policy arc of each representative, by arc of `src` + ne = length(src) + live = falses(ne) + cs = zeros(Int, ne) + cd = zeros(Int, ne) + cw = zeros(T, ne) + reps = Int[] + nnodes = nv + while nnodes > 1 + empty!(reps) + for v in 1:nv + find(v) == v && (push!(reps, v); pos[v] = length(reps)) + end + for e in 1:ne + r, q = find(src[e]), find(dst[e]) + live[e] = r != q + live[e] || continue + cs[e], cd[e] = pos[r], pos[q] + cw[e] = w[e] - x[src[e]] + x[dst[e]] + end + nr = length(reps) + ptr, adj = _out_arcs(nr, cs, e -> live[e]) + cpol = [pol[r] != 0 && live[pol[r]] ? pol[r] : 0 for r in reps] + _, sh, cyc = _howard!(cpol, nr, ptr, adj, cd, cw, tol) + for v in 1:nv + x[v] += sh[pos[find(v)]] + end + for (k, r) in enumerate(reps) + pol[r] = cpol[k] + end + root = reps[cyc[1]] + for k in cyc + parent[reps[k]] = root + end + pol[root] = 0 + nnodes -= length(cyc) - 1 + end + x .-= x[1] + return x +end + +# The scales `y` of the difference system `sys` shifted by the max-balancing +# shift of each strongly connected component, `scc` and `nscc` as returned by +# `_strong_components`. The shift of each component is fixed by making its +# first vertex keep its scale. +function _maxbalance_components(sys::EdgeList{T}, y::Vector{T}, scc::Vector{Int}, nscc::Int) where {T} + P = first.(sys.edges) + Q = last.(sys.edges) + C = sys.cvals + ynew = copy(y) + vptr, verts = _out_arcs(nscc, scc, _ -> true) + eptr, edges = _out_arcs(nscc, [scc[P[e]] for e in eachindex(P)], e -> scc[P[e]] == scc[Q[e]]) + loc = zeros(Int, length(y)) + for k in 1:nscc + vr = vptr[k]:vptr[k+1]-1 + length(vr) > 1 || continue + for (t, s) in enumerate(vr) + loc[verts[s]] = t + end + er = eptr[k]:eptr[k+1]-1 + src = [loc[P[edges[s]]] for s in er] + dst = [loc[Q[edges[s]]] for s in er] + w = [C[edges[s]] - y[P[edges[s]]] + y[Q[edges[s]]] for s in er] + x = _maxbalance(length(vr), src, dst, w) + for (t, s) in enumerate(vr) + ynew[verts[s]] += x[t] + end + end + return ynew +end + +# Feasible scales `Y` of the components of the difference system `csys` between +# strongly connected components, whose edges go from higher to lower labels: +# each component, from the sources down, takes the largest scale at most zero +# that covers the entries entering it. +function _feasible_component_scales(csys::EdgeList{T}, nscc::Int) where {T} + Y = zeros(T, nscc) + iptr, inarcs = _out_arcs(nscc, last.(csys.edges), _ -> true) + for k in nscc:-1:1 + for s in iptr[k]:iptr[k+1]-1 + e = inarcs[s] + j, _ = csys.edges[e] + j > k || error("internal error: an edge between components does not follow their order") + Y[k] = min(Y[k], Y[j] - csys.cvals[e]) + end + end + return Y +end diff --git a/src/stats.jl b/src/stats.jl index 76c6fab..8042f1f 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -276,6 +276,11 @@ group's multiplier test succeeds. - `nreleases::Int`: multiplier tests that released active entries. - `nactive::Int`: entries at zero slack at the solution, not counting those that lie on some maximum-product transversal. +- `maxbalanced::Bool`: the scales within each fully indecomposable block are + max-balanced (`maxbalance=true`); the Newton steps, active-set changes, and + `nactive` then concern only the entries joining different blocks. +- `nblocks::Int`: the fully indecomposable blocks of `A` (the diagonal blocks of + its finest block-triangular form). - The remaining fields describe the linear algebra; their names and meanings may change between releases: `nfactor` (factorizations) and `linsolve` (the solver of the last Newton system: `:dense`, `:cholesky`, `:cg`, or `:none` @@ -292,6 +297,8 @@ struct BalancedCoverStats{T} <: AbstractCoverStats nactivations::Int nreleases::Int nactive::Int + maxbalanced::Bool + nblocks::Int linsolve::Symbol end @@ -303,7 +310,7 @@ tolerance(s::BalancedCoverStats) = s.tol # Conversion from the statistics tuple of `_balanced_newton!`. BalancedCoverStats{T}(nt::NamedTuple) where {T} = BalancedCoverStats{T}(nt.converged, nt.stalled, nt.tol, nt.gradnorm, nt.step, nt.nnewton, nt.nfactor, - nt.nactivations, nt.nreleases, nt.nactive, nt.linsolve) + nt.nactivations, nt.nreleases, nt.nactive, nt.maxbalanced, nt.nblocks, nt.linsolve) """ ExternalSolverStats <: AbstractCoverStats diff --git a/test/maxbalance.jl b/test/maxbalance.jl new file mode 100644 index 0000000..7a79d47 --- /dev/null +++ b/test/maxbalance.jl @@ -0,0 +1,323 @@ +# Transversal-tight covers max-balanced within the fully indecomposable blocks. +# Uses `max_product_transversals` and `logproduct_gap` from transversal.jl, and +# `minimax_ratio`, `offtransversal_ratios`, and `check_balanced_cover` from +# balanced.jl. + +using MatrixCovers: _cover_transversal_balanced, _transversal_system, _howard!, _out_arcs, TIGHT_ULPS + +# Whether the square `A`, whose diagonal is nonzero, is fully indecomposable: +# whether every row reaches every other through the nonzero entries. +function reach(A) + n = size(A, 1) + R = (A .!= 0) .| Matrix(I, n, n) + for k in 1:n, i in 1:n, j in 1:n + R[i, j] |= R[i, k] & R[k, j] + end + return R +end +fully_indecomposable(A) = all(reach(A)) +# The fully indecomposable blocks of `A` (nonzero diagonal), as lists of rows. +function indecomposable_blocks(A) + R = reach(A) + return unique([findall(k -> R[i, k] && R[k, i], axes(A, 1)) for i in axes(A, 1)]) +end + +# Largest discrepancy, over every partition of `rows` into two nonempty sets, of +# the log of the largest ratio of an entry leading from one set to the other +# and of the largest leading back. Entry `A[i,j]` off the transversal `σ` leads +# from row `i` to the row `k` with `σ[k] == j`; only entries between `rows` +# count. +function cut_imbalance(a, b, A, σ, rows) + τ = invperm(σ) + E = [(i, τ[j], log(abs(A[i, j]) / (a[i] * b[j]))) for i in rows, j in axes(A, 2) + if !iszero(A[i, j]) && j != σ[i] && τ[j] in rows] + m = length(rows) + worst = 0.0 + for mask in 1:(2^m-2) + S = Set(rows[k] for k in 1:m if isodd(mask >> (k - 1))) + out = maximum((w for (u, v, w) in E if u in S && !(v in S)); init=-Inf) + inn = maximum((w for (u, v, w) in E if !(u in S) && v in S); init=-Inf) + worst = max(worst, out == inn ? 0.0 : abs(out - inn)) + end + return worst +end + +# A random fully indecomposable `n×n` matrix with a nonzero diagonal. +function random_indecomposable(rng, n) + while true + A = exp.(2 .* randn(rng, n, n)) .* (rand(rng, n, n) .< 0.5) + for i in 1:n + A[i, i] = exp(2 * randn(rng)) + end + fully_indecomposable(A) && return A + end +end + +@testset "policy iteration with nearly equal cycle means" begin + # Row 1 leads to self-loops whose means differ by less than the resolution + # `tol` pairwise but not overall, with the bias favoring the lowest; each + # loop returns to row 1 only through a loop of nearly equal mean. + src = [1, 1, 1, 2, 3, 4, 2, 3, 4] + dst = [2, 3, 4, 2, 3, 4, 1, 1, 3] + tol = TIGHT_ULPS * eps() * 5 + w = [-2.0, -1.5, -1.2, -1 + 1.2tol, -1 + 0.6tol, -1.0, -5.0, -5.0, -5.0] + ptr, adj = _out_arcs(4, src, _ -> true) + λ, s, cyc = _howard!(zeros(Int, 4), 4, ptr, adj, dst, w, tol) + @test λ == w[4] && cyc == [2] + @test all(e -> w[e] - s[src[e]] + s[dst[e]] <= λ + 4tol, eachindex(w)) +end + +@testset "cover_transversal(::PowerMean; maxbalance=true)" begin + # Example 3.11 of Hook, Pestana, Tisseur & Hogg (SIAM J. Matrix Anal. Appl. + # 40:320–346, 2019): the max-balanced Hungarian scaling of `H`. + H = [1 1 1; exp(-1) 1 exp(-2); 0 exp(-4) 1] + M = exp.([0 -0.5 -2.25; -0.5 0 -3.75; -Inf -2.25 0]) + a, b, st = cover_transversal(PowerMean{2}(), H; maxbalance=true) + @test MatrixCovers.converged(st) + @test st.maxbalanced && st.nblocks == 1 && st.nnewton == 0 + @test H ./ (a .* b') ≈ M rtol = 1e-14 + # Covariance: the same ratios in other units. + D1, D2 = Diagonal([1e3, 2.0, 1e-5]), Diagonal([7.0, 1e-4, 3.0]) + a, b = cover_transversal(PowerMean{2}(), D1 * H * D2; maxbalance=true) + @test (D1 * H * D2) ./ (a .* b') ≈ M rtol = 1e-12 + + # Fully indecomposable matrices: the ratios balance across every cut, the + # largest ratio is the minimax ratio `θ`, and no cover pushed toward a + # target has a smaller largest ratio. + rng = StableRNG(11) + for _ in 1:100 + n = rand(rng, 2:8) + A = random_indecomposable(rng, n) + a, b, st = cover_transversal(PowerMean{2}(), A; maxbalance=true) + @test MatrixCovers.converged(st) + @test st.maxbalanced && st.nblocks == 1 + @test check_balanced_cover(a, b, A) + θ, σ = minimax_ratio(A) + @test cut_imbalance(a, b, A, σ, 1:n) <= 1e-10 + @test maximum(last, offtransversal_ratios(a, b, A, σ)) ≈ θ rtol = 1e-10 + # The products do not depend on `ϕ` or `target`. + a1, b1 = cover_transversal(PowerMean{1}(), A; target=0.5, maxbalance=true) + @test logproduct_gap(a, b, a1, b1, A) <= 1e-12 + a2, b2 = cover_transversal(PowerMean{2}(), A) + @test maximum(last, offtransversal_ratios(a2, b2, A, σ)) >= θ * (1 - 1e-12) + end + + # Agreement with the max-balanced Hungarian scaling of Hook et al. as + # computed by an independent implementation of their Algorithm 3.10 + # (`maxbalanced_hungarian` in `src/maxbalance.jl` of the paper "Numerical + # methods for scale-invariant analysis"); these are its products `a .* b'`. + refs = [([4.0 1 0 2; 0 3 5 0; 1 0 2 7; 6 1 0 1], + [4.898979485566356 1.0 1.1665046746253358 2.857539627437798; + 20.998542021016053 4.286309441156698 4.999999999999999 12.248299083566028; + 12.000833188694234 2.4496598167132055 2.8575396274377978 6.999999999999999; + 6.0 1.2247448713915892 1.4286706177016946 3.4997570035026757]), + ([1.0 8 0 0 2; 0 1 3 0 0; 0 0 2 9 0; 5 0 0 1 4; 0 3 0 1 6], + [2.714417616594907 7.999999999999998 5.009957159878615 4.7061882644660615 5.894450397824619; + 1.625413669202318 4.790460124529586 3.000000000000001 2.818100902431721 3.52964119834955; + 5.190986245452685 15.299005469805346 9.580920249059176 9.000000000000002 11.272403609726895; + 5.0 14.73612599456154 9.228419991915873 8.66887290241235 10.857670466379625; + 2.7630236239802906 8.143252849784718 5.0996684899351115 4.790460124529582 6.0]), + ([2.0 1e3 0; 1e-3 5 7; 3 1e-2 1], + [25.61158329974988 999.9999999999998 109.32553318667111; + 1.6398829978000684 64.02895824937472 6.999999999999997; + 3.000000000000001 117.13449984286204 12.805791649874939])] + for (A, P) in refs + a, b = cover_transversal(PowerMean{2}(), A; maxbalance=true) + @test a .* b' ≈ P rtol = 1e-10 + end + + # Partly decomposable: three blocks in block upper-triangular form. + A11 = [4.0 1 0.5; 1 4 1; 0.5 1 4] + A22 = [3.0 0.5; 1 5] + A33 = [2.0 0 1; 1 3 0; 0 0.25 1] + A12 = [1e-3 5.0; 2.0 0; 0 30.0] + A13 = [0 1e2 0; 0 0 0; 4.0 0 0.1] + A23 = [0.5 0 0; 0 1e-2 2.0] + A = [A11 A12 A13; zeros(2, 3) A22 A23; zeros(3, 5) A33] + blocks = (1:3, 4:5, 6:8) + for (p, t) in ((2, 1e-2), (1, 0.1), (3, 0.5)) + a, b, st = cover_transversal(PowerMean{p}(), A; target=t, maxbalance=true) + @test MatrixCovers.converged(st) + @test st.maxbalanced && st.nblocks == 3 + @test check_balanced_cover(a, b, A) + R = A ./ (a .* b') + # The KKT conditions across the cuts after the first and second blocks. + for k in (3, 5) + rcut = [R[i, j] for i in 1:k, j in k+1:8 if A[i, j] != 0] + @test maximum(rcut) < 1 + @test sum(r -> r^p, rcut) / length(rcut) ≈ t^p rtol = 1e-8 + end + # Within each block, the products are those of the block alone. + for r in blocks + ab, bb = cover_transversal(PowerMean{p}(), A[r, r]; target=t, maxbalance=true) + @test logproduct_gap(a[r], b[r], ab, bb, A[r, r]) <= 1e-12 + @test cut_imbalance(a[r], b[r], A[r, r], collect(1:length(r)), 1:length(r)) <= 1e-10 + end + # Covariance under power-of-two changes of units, to roundoff with a + # stopping tolerance near it. + rng = StableRNG(p) + d, e = 2.0 .^ rand(rng, -30:30, 8), 2.0 .^ rand(rng, -30:30, 8) + a12, b12 = cover_transversal(PowerMean{p}(), A; target=t, maxbalance=true, tol=1e-12) + ã, b̃ = cover_transversal(PowerMean{p}(), Diagonal(d) * A * Diagonal(e); target=t, maxbalance=true, + tol=1e-12) + @test logproduct_gap(ã ./ d, b̃ ./ e, a12, b12, A) <= 1e-12 + # Permutation equivariance. + pr, pc = randperm(rng, 8), randperm(rng, 8) + ap, bp = cover_transversal(PowerMean{p}(), A[pr, pc]; target=t, maxbalance=true) + @test (ap .* bp') ≈ (a .* b')[pr, pc] rtol = 1e-10 + end + + # A cut whose largest entry would exceed its cover at the unconstrained + # minimizer: it is held at ratio one, and the p-power mean across the cut + # is below the target, so that raising the ratios further would decrease + # the objective. + Bk = [2.0 1; 1 2] + A = [Bk [1.0 0; 1e-6 1e-6]; zeros(2, 2) Bk] + for (p, t) in ((2, 0.9), (1, 0.5)) + a, b, st = cover_transversal(PowerMean{p}(), A; target=t, maxbalance=true) + @test MatrixCovers.converged(st) + @test st.nblocks == 2 && st.nactive == 1 + @test check_balanced_cover(a, b, A) + rcut = [A[i, j] / (a[i] * b[j]) for i in 1:2, j in 3:4 if A[i, j] != 0] + @test maximum(rcut) ≈ 1 rtol = 1e-12 + @test sum(r -> r^p, rcut) / length(rcut) < t^p + end + + # Ties: 0/1 matrices, whose maximum-product transversals are many and + # whose rows are joined into tight components. + rng = StableRNG(5) + ntie = ntight = 0 + for _ in 1:40 + n = rand(rng, 3:7) + A = Float64.(rand(rng, n, n) .< 0.4) + for i in 1:n + A[i, i] = 1 + end + a, b, st = cover_transversal(PowerMean{2}(), A; maxbalance=true) + @test MatrixCovers.converged(st) + @test check_balanced_cover(a, b, A) + blks = indecomposable_blocks(A) + @test st.nblocks == length(blks) + for r in blks + @test cut_imbalance(a, b, A, collect(1:n), r) <= 1e-10 + end + σs, _ = max_product_transversals(A) + length(σs) > 1 && (ntie += 1) + _transversal_system(_row_support(A, Float64), axes(A, 2), Float64, :test, true).nc < n && (ntight += 1) + for σ in σs[1:min(end, 4)] + aσ, bσ, _ = _cover_transversal_balanced(A, 2, 1e-2; transversal=σ, maxbalance=true) + @test logproduct_gap(a, b, aσ, bσ, A) <= 1e-12 + end + end + @test ntie > 10 && ntight > 10 + # Generic magnitudes with tied transversals. + ntie = 0 + for (_, A) in general_matrices + ntie >= 10 && break + Af = Float64.(A) + σs, _ = max_product_transversals(Af) + length(σs) > 1 || continue + ntie += 1 + a, b, st = cover_transversal(PowerMean{2}(), Af; maxbalance=true) + @test MatrixCovers.converged(st) + @test check_balanced_cover(a, b, Af) + for σ in σs + aσ, bσ, _ = _cover_transversal_balanced(Af, 2, 1e-2; transversal=σ, maxbalance=true) + @test logproduct_gap(a, b, aσ, bσ, Af) <= 1e-10 + end + end + @test ntie == 10 + + # Small cases: a 1×1 matrix, a diagonal matrix, and a triangular one. + a, b, st = cover_transversal(PowerMean{2}(), fill(3.0, 1, 1); maxbalance=true) + @test a[1] * b[1] ≈ 3 + @test st.nblocks == 1 && st.nnewton == 0 && MatrixCovers.converged(st) + D = [2.0 0 0; 0 3 0; 0 0 5] + a, b, st = cover_transversal(PowerMean{2}(), D; maxbalance=true) + @test a .* b' .* (D .!= 0) ≈ D + @test st.nblocks == 3 && st.nactive == 0 && MatrixCovers.converged(st) + A = [1.0 1e6; 0 1.0] + a, b, st = cover_transversal(PowerMean{2}(), A; target=1e-2, maxbalance=true) + @test A[1, 2] / (a[1] * b[2]) ≈ 1e-2 + @test st.nblocks == 2 + + # Sparse matrices: sparse and dense storage, changes of units, and the + # balance of each row of the blocks between its largest ratio leaving and + # entering. + for seed in 1:2 + rng = StableRNG(seed) + n = 200 + A = sprandn(rng, n, n, 0.01) + spdiagm(0 => exp.(3 .* randn(rng, n))) + a, b, st = cover_transversal(PowerMean{2}(), A; maxbalance=true) + @test MatrixCovers.converged(st) + @test st.nblocks > 1 + @test check_balanced_cover(a, b, A; atol=1e-8) + @test isbalanced(a, b, A) + ad, bd = cover_transversal(PowerMean{2}(), Matrix(A); maxbalance=true) + @test logproduct_gap(a, b, ad, bd, A) <= 1e-10 + d, e = 2.0 .^ rand(rng, -30:30, n), 2.0 .^ rand(rng, -30:30, n) + ã, b̃ = cover_transversal(PowerMean{2}(), Diagonal(d) * A * Diagonal(e); maxbalance=true) + @test logproduct_gap(ã ./ d, b̃ ./ e, a, b, A) <= 1e-10 + σ, _ = MatrixCovers._max_product_transversal(_row_support(A, Float64), axes(A, 2), :test) + τ = invperm(σ) + R = reach(Matrix(A[:, σ]) .!= 0) + lout, lin = fill(-Inf, n), fill(-Inf, n) + for (i, j, v) in zip(findnz(A)...) + k = τ[j] + (k != i && R[i, k] && R[k, i]) || continue + w = log(abs(v) / (a[i] * b[j])) + lout[i] = max(lout[i], w) + lin[k] = max(lin[k], w) + end + @test all(k -> lout[k] == lin[k] || abs(lout[k] - lin[k]) <= 1e-10, 1:n) + @test count(isfinite, lout) > n ÷ 2 + end + + # Generic indexing, element types, and the mutating form. + A = [1.0 0.5 0 0.3; 0.2 1.0 0 0; 0 0 1.0 0.5; 0 0 0.4 1.0] + a, b, st = cover_transversal(PowerMean{2}(), A; maxbalance=true) + @test st isa MatrixCovers.BalancedCoverStats{Float64} + @test st.maxbalanced && st.nblocks == 2 + @test !cover_transversal(PowerMean{2}(), A)[3].maxbalanced + @test cover_transversal(PowerMean{2}(), A)[3].nblocks == 2 + Ao = OffsetArray(A, -1:2, 5:8) + ao, bo = cover_transversal(PowerMean{2}(), Ao; maxbalance=true) + @test axes(ao, 1) == axes(Ao, 1) && axes(bo, 1) == axes(Ao, 2) + @test parent(ao) .* parent(bo)' ≈ a .* b' rtol = 1e-12 + B = [9.0 9 9 9 9; 9 1 0.5 0 0.3; 9 0.2 1 0 0; 9 0 0 1 0.5; 9 0 0 0.4 1] + av, bv = cover_transversal(PowerMean{2}(), view(B, 2:5, 2:5); maxbalance=true) + @test av .* bv' ≈ a .* b' rtol = 1e-12 + as, bs = cover_transversal(PowerMean{2}(), sparse(A); maxbalance=true) + @test as .* bs' ≈ a .* b' rtol = 1e-12 + ao, bo = cover_transversal!(PowerMean{2}(), OffsetArray(zeros(4), -1:2), OffsetArray(zeros(4), 5:8), Ao; + maxbalance=true) + @test parent(ao) .* parent(bo)' ≈ a .* b' rtol = 1e-12 + a32, b32, st32 = cover_transversal(PowerMean{2}(), Float32.(A); maxbalance=true) + @test eltype(a32) == eltype(b32) == Float32 + @test st32.maxbalanced + @test iscover(a32, b32, Float32.(A)) + @test a32 .* b32' ≈ a .* b' rtol = 1e-6 + aB, bB, stB = cover_transversal(PowerMean{2}(), BigFloat.(A); maxbalance=true) + @test eltype(aB) == BigFloat + @test aB .* bB' ≈ a .* b' rtol = 1e-12 + a2, b2 = zeros(4), zeros(4) + _, _, st2 = cover_transversal!(PowerMean{2}(), a2, b2, A; maxbalance=true) + @test (a2, b2) == (a, b) + @test st2.maxbalanced + + # An iteration limit on the problem between blocks: a warning, and still a + # hard cover tight on the transversal. + A = [Bk [1.0 0.3; 0.2 0.1]; zeros(2, 2) Bk] + a, b, st = @test_logs (:warn, r"may not minimize the objective") cover_transversal(PowerMean{2}(), A; + maxbalance=true, maxiter=0) + @test !MatrixCovers.converged(st) + @test check_balanced_cover(a, b, A) + + # Errors. + @test_throws "requires a finite exponent" cover_transversal(PowerMean{Inf}(), A; maxbalance=true) + @test_throws "target must satisfy 0 < target <= 1, got 0" cover_transversal(PowerMean{2}(), A; target=0, + maxbalance=true) + @test_throws "structurally nonsingular" cover_transversal(PowerMean{2}(), [1.0 1 1; 1 0 0; 1 0 0]; + maxbalance=true) + @test_throws "requires a square matrix" cover_transversal(PowerMean{2}(), ones(2, 3); maxbalance=true) +end diff --git a/test/runtests.jl b/test/runtests.jl index 8cc9b98..1f2af28 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,6 +31,7 @@ include("helpers.jl") # isbalanced, covaries, PENALTIES include("minimal_covers.jl") # the *_min family (native solvers) include("transversal.jl") # cover_transversal and the maximum-product transversal include("balanced.jl") # cover_transversal with a PowerMean penalty and a target + include("maxbalance.jl") # cover_transversal(PowerMean; maxbalance=true): max-balanced blocks include("active_set.jl") # exact solution of stalled AbsLog{2} minimal covers include("storage_types.jl") # sparse/structured/wrapped storage vs dense reference include("element_types.jl") # Float32/BigFloat: eltype-scaled internal tolerances From e251706906faa514935f16e0c2925b4439a8252d Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 23 Oct 2026 14:45:28 -0500 Subject: [PATCH 5/5] Bump version to 1.2.0 cover_transversal gains the PowerMean method and its maxbalance keyword, new public functionality since 1.1.0. Assisted-by: Claude Fable 5.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 50a5b7b..e7109cb 100644 --- a/Project.toml +++ b/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 and contributors"] [deps]