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] 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..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}` | — | +| [`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 | @@ -66,7 +66,14 @@ For hard covers, [`symcover`](@ref) and [`cover`](@ref) are fast heuristics; their `_min` counterparts minimize the selected objective. [`cover_transversal`](@ref) minimizes the `AbsLog{2}` objective among the hard covers that are tight on a maximum-product transversal (see -[Transversals and conditioning](@ref)). A soft cover has no constraint to +[Transversals and conditioning](@ref)); with a `PowerMean{p}` penalty and a +`target` ratio `t`, it instead pushes the ratios of the entries off the +transversal toward `t`, minimizing `∑ t^p φ(r/t)` over them with `φ` the +`PowerMean{p}` penalty. Adding `maxbalance=true` replaces that objective within +each fully indecomposable block by its `p → ∞` limit, max-balancing: across +every partition of a block's rows, the largest ratio in one direction equals +the largest in the other, which also minimizes the block's largest ratio. The +`PowerMean{p}` objective then fixes only the scales between blocks. A soft cover has no constraint to satisfy, so there is no separate heuristic: [`soft_symcover`](@ref) and [`soft_cover`](@ref) always minimize the objective. @@ -78,6 +85,11 @@ The method used to minimize the objective depends on the penalty: converge slowly. - `AbsLog{2}` (native) is a single linear least-squares solve for soft covers; for hard covers, see [Objective-minimal covers](@ref). +- `PowerMean{p}` with a target ratio, for [`cover_transversal`](@ref) (native), + uses an active-set Newton method that keeps every iterate a hard cover + tight on the transversal. With `maxbalance=true`, max-plus policy iteration + max-balances each block first, and the Newton method acts on one scale per + block. - `AbsLog{1}` is a linear program solved with JuMP and HiGHS. - `AbsLinear` is not convex; several starts are each refined to a local minimum with JuMP and Ipopt, and the best is returned. 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..42840dd 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,8 @@ include("tree_cholesky.jl") # fill-free factorization of forest-patterned mat include("minimal_covers.jl") include("active_set.jl") # exact solution of the AbsLog{2} problem from a near-optimal point include("transversal.jl") # transversal-tight minimal covers +include("balanced.jl") # transversal-tight covers pushed toward a target ratio +include("maxbalance.jl") # max-balancing within the fully indecomposable blocks include("transversal_heuristic.jl") # the transversal-tight start of `cover` include("powermean_newton.jl") # Newton refinement of the PowerMean soft covers include("sparse_support.jl") # sparse traversal and the sparse solver defaults diff --git a/src/balanced.jl b/src/balanced.jl new file mode 100644 index 0000000..0820693 --- /dev/null +++ b/src/balanced.jl @@ -0,0 +1,814 @@ +# 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, 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 +`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. + +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 +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. + +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 +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, + 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; maxbalance, 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; 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, + 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; 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) + _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) + # 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) + 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/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/minimal_covers.jl b/src/minimal_covers.jl index 8fd8aa8..7b8c9dd 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 @@ -2405,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..8042f1f 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -241,6 +241,77 @@ 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. +- `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` + 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 + maxbalanced::Bool + nblocks::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.maxbalanced, nt.nblocks, nt.linsolve) + """ ExternalSolverStats <: AbstractCoverStats diff --git a/src/transversal.jl b/src/transversal.jl index ae3a58e..e337ce7 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 @@ -301,6 +304,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 @@ -335,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) @@ -423,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 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/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 96e6e7e..1f2af28 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,6 +30,8 @@ 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("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 @@ -49,7 +51,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