diff --git a/docs/src/objectives.md b/docs/src/objectives.md index c945fcf..9d30aba 100644 --- a/docs/src/objectives.md +++ b/docs/src/objectives.md @@ -87,9 +87,9 @@ The method used to minimize the objective depends on the penalty: for hard covers, see [Objective-minimal covers](@ref). - `PowerMean{p}` with a target ratio, for [`cover_transversal`](@ref) (native), uses an active-set Newton method that keeps every iterate a hard cover - tight on the transversal. With `maxbalance=true`, max-plus policy iteration - max-balances each block first, and the Newton method acts on one scale per - block. + tight on the transversal. With `maxbalance=true`, a parametric shortest-path + algorithm max-balances each block first, and the Newton method acts on one + scale per block. - `AbsLog{1}` is a linear program solved with JuMP and HiGHS. - `AbsLinear` is not convex; several starts are each refined to a local minimum with JuMP and Ipopt, and the best is returned. diff --git a/src/active_set.jl b/src/active_set.jl index 183bd73..793cba0 100644 --- a/src/active_set.jl +++ b/src/active_set.jl @@ -386,23 +386,50 @@ end function _laplacian_factor(L::SparseMatrixCSC{T,Int}; fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET, cgmaxiter::Int=10 * size(L, 1) + 100) where {T} size(L, 1) == 0 && return _laplacian_factor(Matrix(L)) + return _laplacian_cg(L; cgmaxiter) +end +# Conjugate-gradient solver for `L`, built without a symbolic analysis, for a +# pattern presumed to exceed the budgets. The fallback factorization is analyzed +# only if CG fails (see `_cg_fallback`). +function _laplacian_cg(L::SparseMatrixCSC{T,Int}; cgmaxiter::Int=10 * size(L, 1) + 100) where {T} S = L + triu(L, 1)' return _LaplacianCG{T}(L, S, Vector(diag(S)), nothing, Ref{Any}(nothing), cgmaxiter, Ref(zero(T))) end +# The factorization that replaces conjugate gradients once they fail. In +# `Float64` it reuses the symbolic analysis in `C.F` when there is one. +function _cg_fallback(C::_LaplacianCG{Float64}) + F = C.F + if F === nothing + F = SparseCholesky() + analyze!(F, C.L) + end + return factorize!(F, C.L) +end +_cg_fallback(C::_LaplacianCG) = _laplacian_factor(Matrix(C.L)) _laplacian_factor(L::Matrix) = LinearAlgebra.cholesky!(Symmetric(L, :U)) -_laplacian_solve!(x::Vector{Float64}, F::SparseCholesky, b::Vector{Float64}) = +# Solves `L x = b`. The keywords apply only to a `_LaplacianCG`: `rtol` is the +# relative residual at which conjugate gradients stop, and `refine` enables the +# iterative refinement that follows, which reduces the forward error to the +# resolution of the arithmetic. `_solve_error` then estimates the remaining +# error; it is `NaN` when refinement is skipped. +_laplacian_solve!(x::Vector{Float64}, F::SparseCholesky, b::Vector{Float64}; rtol=nothing, refine=true) = isempty(b) ? x : solve!(x, F, CHOLMOD_A, b) -_laplacian_solve!(x, F::Cholesky, b) = ldiv!(x, F, b) +_laplacian_solve!(x, F::Cholesky, b; rtol=nothing, refine=true) = ldiv!(x, F, b) -function _laplacian_solve!(x::Vector{T}, C::_LaplacianCG{T}, b::Vector{T}) where {T} +function _laplacian_solve!(x::Vector{T}, C::_LaplacianCG{T}, b::Vector{T}; rtol::Real=sqrt(eps(T)), + refine::Bool=true) where {T} if C.fallback[] === nothing n = length(b) fill!(x, zero(T)) rdd, r, z, d, Ad, e = ntuple(_ -> zeros(T, n), 6) Smul! = (y, v) -> mul!(y, C.S, v) - η = sqrt(eps(T)) + η = T(rtol) _, ok = _pcg!(Smul!, x, C.dg, b, r, z, d, Ad, C.maxiter, η * norm(b)) + if ok && !refine + C.err[] = T(NaN) + return x + end # Iterative refinement with the residual in compensated arithmetic. # CG leaves a forward error of its normwise backward error times the # condition number, and the multipliers, routed sums of the solution @@ -426,7 +453,7 @@ function _laplacian_solve!(x::Vector{T}, C::_LaplacianCG{T}, b::Vector{T}) where δ <= δprev / 100 || break δprev = δ end - C.fallback[] = C.F === nothing ? _laplacian_factor(Matrix(C.L)) : factorize!(C.F, C.L) + C.fallback[] = _cg_fallback(C) end C.err[] = zero(T) return _laplacian_solve!(x, C.fallback[], b) diff --git a/src/balanced.jl b/src/balanced.jl index 0820693..489087c 100644 --- a/src/balanced.jl +++ b/src/balanced.jl @@ -15,6 +15,11 @@ 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 +# Relative residual of a conjugate-gradient solve for a Newton direction. Every +# iterate of conjugate gradients from zero is a descent direction, and the +# reduced gradient is recomputed at each step, so the direction need not be +# accurate. +const BALANCED_CG_RTOL = 1e-4 # 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 @@ -93,12 +98,17 @@ 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. +With `maxbalance=true`, each block is max-balanced as in Algorithm 3.10 of Hook +et al.: a cycle of largest geometric-mean ratio in the block is found, 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 cycles are found by the parametric +shortest-path algorithm of Young, Tarjan & Orlin (Networks 21:205–221, 1991). +For a block with `n` rows and `m` entries, its cost is at most +`O(nm + n² log n)`. The cost is roughly linear in `m` when the block's +shortest-path tree stays shallow, as it does for most sparse matrices, and +about `O(nm)` for a banded block, whose tree is a path. 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 @@ -110,7 +120,8 @@ nonnegative but within the roundoff of the groups' derivatives. `maxiter` 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)` +gradients to a loose tolerance, since a descent direction need not be +accurate. 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. @@ -663,7 +674,7 @@ function _balanced_direction!(δv, cache, d, gid, uidx, nu, gG, P, Q, p::T, tp:: 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) + _laplacian_solve!(x, F, rhs; rtol=BALANCED_CG_RTOL, refine=false) kind = F isa SparseCholesky ? :cholesky : :cg break end @@ -682,17 +693,27 @@ end # 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. +# Factor the pinned Laplacian `L` (upper triangle), or return a +# conjugate-gradient solver for it when its factor would exceed the budgets. +# Returns `nothing` when `L` is not numerically positive definite. `cache` holds +# one of: +# - the analyzed factorization and its pattern, reused when the pattern of `L` +# is unchanged; +# - the order and number of nonzeros of the last pattern whose factor exceeded +# the budgets. A pattern more than half that size in both is presumed to +# exceed them too and gets conjugate gradients without a symbolic analysis, +# which would cost more than the solve. 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) + if c isa Tuple{SparseCholesky,Vector{Int},Vector{Int}} + c[2] == SparseArrays.getcolptr(L) && c[3] == rowvals(L) && return factorize!(c[1], L) + elseif c isa Tuple{Int,Int} + 2 * size(L, 1) > c[1] && 2 * nnz(L) > c[2] && return _laplacian_cg(L) end F = _laplacian_factor(L; fillbudget, flopbudget) - cache[] = F isa SparseCholesky ? (F, copy(SparseArrays.getcolptr(L)), copy(rowvals(L))) : nothing + cache[] = F isa SparseCholesky ? (F, copy(SparseArrays.getcolptr(L)), copy(rowvals(L))) : + (size(L, 1), nnz(L)) return F catch err err isa LinearAlgebra.PosDefException || rethrow() diff --git a/src/maxbalance.jl b/src/maxbalance.jl index 8aee864..c14fad8 100644 --- a/src/maxbalance.jl +++ b/src/maxbalance.jl @@ -1,6 +1,6 @@ # 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). +# connected components: the max-balanced Hungarian scaling of Hook, Pestana, +# Tisseur & Hogg (SIAM J. Matrix Anal. Appl. 40:320–346, 2019). # # 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 @@ -30,172 +30,408 @@ function _out_arcs(nv::Int, P::AbstractVector{Int}, keep) 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. +# Indexed binary max-heap over the integers `1:n`: `key[e]` is the key of `e`, +# `heap` holds the members, and `pos[e]` is the position of `e` in `heap`, or 0. +struct IndexedHeap{T} + key::Vector{T} + heap::Vector{Int} + pos::Vector{Int} +end +IndexedHeap{T}(n::Int) where {T} = IndexedHeap{T}(Vector{T}(undef, n), Int[], zeros(Int, n)) + +function _heap_swap!(h::IndexedHeap, i::Int, j::Int) + h.heap[i], h.heap[j] = h.heap[j], h.heap[i] + h.pos[h.heap[i]] = i + h.pos[h.heap[j]] = j + return nothing +end +function _heap_up!(h::IndexedHeap, i::Int) + while i > 1 + p = i >> 1 + h.key[h.heap[p]] >= h.key[h.heap[i]] && break + _heap_swap!(h, p, i) + i = p + end + return nothing +end +function _heap_down!(h::IndexedHeap, i::Int) + n = length(h.heap) + while true + c = 2i + c <= n || break + c < n && h.key[h.heap[c+1]] > h.key[h.heap[c]] && (c += 1) + h.key[h.heap[c]] > h.key[h.heap[i]] || break + _heap_swap!(h, c, i) + i = c + end + return nothing +end +# Inserts `e` with key `k`, or changes its key. +function _heap_set!(h::IndexedHeap, e::Int, k) + h.key[e] = k + i = h.pos[e] + if i == 0 + push!(h.heap, e) + h.pos[e] = length(h.heap) + _heap_up!(h, length(h.heap)) + else + _heap_up!(h, i) + _heap_down!(h, h.pos[e]) + end + return nothing +end +# Removes `e` if present. +function _heap_remove!(h::IndexedHeap, e::Int) + i = h.pos[e] + i == 0 && return nothing + last = pop!(h.heap) + h.pos[e] = 0 + if i <= length(h.heap) + h.heap[i] = last + h.pos[last] = i + _heap_up!(h, i) + _heap_down!(h, h.pos[last]) + end + return nothing +end +_heap_pop!(h::IndexedHeap) = (e = h.heap[1]; _heap_remove!(h, e); e) + +# The parametric shortest-path tree of a digraph whose arcs have length `λ - w` +# for a decreasing parameter `λ`, with each cycle contracted when its length +# reaches zero. The nodes are the vertices `1:nv`; a contracted cycle is +# represented by its node nearest the root, and the other nodes of the cycle are +# no longer live. # -# 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 +# - Arc `e` joins the live nodes `src[e] → dst[e]` with weight `w[e]`; +# `alive[e]` is false once both ends lie in one node. +# - Live node `v` has depth `a[v]`, and the arcs of its tree path have total +# weight `b[v]`, so its distance from the root is `a[v] * λ - b[v]`. `tarc[v]` +# is the tree arc entering `v`, or 0 at the root. +# - The children of `v` form a doubly linked list with first member +# `fchild[v]`, last member `lchild[v]`, and links `nsib`, `psib`. +# - The arcs incident to `v` form a doubly linked list with first member +# `lhead[v]`, last member `ltail[v]`, and links `lnext`, `lprev`. Its members +# are slots: slot `e` is arc `e` as it leaves `src[e]`, and slot +# `e + length(w)` is arc `e` as it enters `dst[e]`. +# - `heap` holds every live non-tree arc `p → q` with `a[p] + 1 > a[q]`, keyed +# by the `λ` at which it becomes tight. +struct ParametricTree{T} + src::Vector{Int} + dst::Vector{Int} + w::Vector{T} + alive::BitVector + a::Vector{Int} + b::Vector{T} + tarc::Vector{Int} + fchild::Vector{Int} + lchild::Vector{Int} + nsib::Vector{Int} + psib::Vector{Int} + lhead::Vector{Int} + ltail::Vector{Int} + lnext::Vector{Int} + lprev::Vector{Int} + heap::IndexedHeap{T} +end +function ParametricTree(nv::Int, src::Vector{Int}, dst::Vector{Int}, w::Vector{T}) where {T} + ne = length(w) + z() = zeros(Int, nv) + return ParametricTree{T}(copy(src), copy(dst), copy(w), BitVector(src[e] != dst[e] for e in 1:ne), + z(), zeros(T, nv), z(), z(), z(), z(), z(), z(), z(), + zeros(Int, 2ne), zeros(Int, 2ne), IndexedHeap{T}(ne)) +end + +function _child_append!(S::ParametricTree, p::Int, c::Int) + S.psib[c] = S.lchild[p] + S.nsib[c] = 0 + S.lchild[p] == 0 ? (S.fchild[p] = c) : (S.nsib[S.lchild[p]] = c) + S.lchild[p] = c + return nothing +end +function _child_remove!(S::ParametricTree, p::Int, c::Int) + pr, nx = S.psib[c], S.nsib[c] + pr == 0 ? (S.fchild[p] = nx) : (S.nsib[pr] = nx) + nx == 0 ? (S.lchild[p] = pr) : (S.psib[nx] = pr) + return nothing +end +# Moves the children of `z` to the end of the children of `p`. +function _child_concat!(S::ParametricTree, p::Int, z::Int) + S.fchild[z] == 0 && return nothing + if S.lchild[p] == 0 + S.fchild[p] = S.fchild[z] + else + S.nsib[S.lchild[p]] = S.fchild[z] + S.psib[S.fchild[z]] = S.lchild[p] 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 + S.lchild[p] = S.lchild[z] + S.fchild[z] = S.lchild[z] = 0 + return nothing +end +function _arcs_append!(S::ParametricTree, v::Int, s::Int) + S.lprev[s] = S.ltail[v] + S.lnext[s] = 0 + S.ltail[v] == 0 ? (S.lhead[v] = s) : (S.lnext[S.ltail[v]] = s) + S.ltail[v] = s + return nothing +end +function _arcs_concat!(S::ParametricTree, v::Int, z::Int) + S.lhead[z] == 0 && return nothing + if S.lhead[v] == 0 + S.lhead[v] = S.lhead[z] + else + S.lnext[S.ltail[v]] = S.lhead[z] + S.lprev[S.lhead[z]] = S.ltail[v] + end + S.ltail[v] = S.ltail[z] + S.lhead[z] = S.ltail[z] = 0 + return nothing +end +function _arcs_unlink!(S::ParametricTree, v::Int, s::Int) + pr, nx = S.lprev[s], S.lnext[s] + pr == 0 ? (S.lhead[v] = nx) : (S.lnext[pr] = nx) + nx == 0 ? (S.ltail[v] = pr) : (S.lprev[nx] = pr) + return nothing +end +# The arc of slot `s` and the node at the other end of it. +function _slot_arc(S::ParametricTree, s::Int) + ne = length(S.w) + return s > ne ? (s - ne, S.src[s-ne]) : (s, S.dst[s]) +end + +# Updates the heap entry of the live arc `e` after its ends have moved. The arc +# `p → q` has reduced length `(a[p] + 1 - a[q]) * λ - (b[p] + w[e] - b[q])`, +# which falls to zero as `λ` decreases only when `a[p] + 1 > a[q]`. Such a +# non-tree arc is keyed by the `λ` at which it does; any other arc is removed. +function _key!(S::ParametricTree, e::Int) + p, q = S.src[e], S.dst[e] + c = S.a[p] + 1 - S.a[q] + if S.tarc[q] == e || c <= 0 + _heap_remove!(S.heap, e) + else + _heap_set!(S.heap, e, (S.b[p] + S.w[e] - S.b[q]) / c) + end + return nothing +end + +# Adds `Δa` to the depth and `Δb` to the path weight of each node `y` of the +# subtree listed in `sub`, marks it with `stamp[y] = st`, and rekeys the arcs +# between the subtree and the rest of the digraph. Arcs within the subtree keep +# their keys, since both ends shift alike. +function _shift_subtree!(S::ParametricTree, sub::Vector{Int}, Δa::Int, Δb, stamp::Vector{Int}, st::Int) + for y in sub + S.a[y] += Δa + S.b[y] += Δb + stamp[y] = st + end + for y in sub + s = S.lhead[y] + while s != 0 + f, o = _slot_arc(S, s) + stamp[o] == st || _key!(S, f) + s = S.lnext[s] 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 + return nothing +end + +# Appends the descendants of `v` to `sub`, with `stack` as workspace. +function _subtree!(sub::Vector{Int}, S::ParametricTree, v::Int, stack::Vector{Int}) + empty!(stack) + c = S.fchild[v] + while c != 0 + push!(stack, c) + c = S.nsib[c] + end + while !isempty(stack) + y = pop!(stack) + push!(sub, y) + c = S.fchild[y] + while c != 0 + push!(stack, c) + c = S.nsib[c] end end - error("internal error: policy iteration did not converge") + return sub 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`. +# `src[e] → dst[e]` of weight `w[e]`, by the minimum-balance algorithm of Young, +# Tarjan & Orlin (Networks 21:205–221, 1991) with the signs reversed; the shift +# is unique up to a constant (Schneider & Schneider, Math. Oper. Res. +# 16:208–222, 1991). Returns `x`, with `x[1] == 0`. +# +# Give arc `e` the length `λ - w[e]` and let `d` be the shortest-path distances +# from vertex 1. Then `w[e] - d[src[e]] + d[dst[e]] <= λ` on every arc, with +# equality on the tree arcs. As `λ` decreases from `+∞`, the tree changes each +# time a non-tree arc becomes tight. If the tail of the arc is not a descendant +# of its head, the arc replaces the tree arc into its head. Otherwise the arc +# and the tree path from its head to its tail form a cycle of zero length, +# whose mean weight `λ` is the largest cycle mean of the digraph left by the +# earlier contractions. The distances within the cycle are then fixed relative +# to one another, and the cycle is contracted into its node nearest the root. +# When node `z` is absorbed, `uf[z]` records the node that absorbed it and +# `uΔ[z]` its distance relative to that node; the shift of `z` is the sum of +# `uΔ` along the chain `z`, `uf[z]`, `uf[uf[z]]`, … up to vertex 1. +# +# The tree starts as the shortest-path tree for `λ → ∞`: fewest arcs, then +# largest weight. Each time an arc becomes tight, the work is proportional to +# the size of the subtrees whose depths change plus the number of arcs incident +# to them: `O(nm + n² log n)` in all, and far less when the tree stays shallow. 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)) + ne = length(w) + S = ParametricTree(nv, src, dst, w) + a, b, tarc, heap = S.a, S.b, S.tarc, S.heap + csrc, cdst, cw = S.src, S.dst, S.w + fill!(a, -1) + ptr, adj = _out_arcs(nv, src, e -> S.alive[e]) + queue = Int[1] + a[1] = 0 + i = 1 + while i <= length(queue) + v = queue[i] + i += 1 + for s in ptr[v]:ptr[v+1]-1 + e = adj[s] + q = dst[e] + if a[q] == -1 + a[q] = a[v] + 1 + push!(queue, q) + end + if a[q] == a[v] + 1 && (tarc[q] == 0 || b[v] + w[e] > b[q]) + b[q] = b[v] + w[e] + tarc[q] = e + end 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 + length(queue) == nv || error("internal error: the digraph is not strongly connected") + for q in 2:nv + _child_append!(S, src[tarc[q]], q) + end + for e in 1:ne + S.alive[e] || continue + _arcs_append!(S, src[e], e) + _arcs_append!(S, dst[e], e + ne) + _key!(S, e) + end + uf = collect(1:nv) + uΔ = zeros(T, nv) + stamp = zeros(Int, nv) + st = 0 + sub = Int[] + stack = Int[] + path = Int[] + touched = Int[] + nlive = nv + while nlive > 1 + isempty(heap.heap) && error("internal error: no arc can become tight on a strongly connected digraph") + e = _heap_pop!(heap) + λ = heap.key[e] + u, v = csrc[e], cdst[e] + # Decide whether `u` is a descendant of `v` by alternating a step up the + # tree from `u` with a step of listing the subtree of `v`, so that a + # cycle costs the smaller of its length and that subtree. If `u` is not + # a descendant, `sub` ends up listing the descendants of `v`. + cycle = false + z = u + empty!(sub) + empty!(stack) + push!(stack, v) + while true + if z != 0 + z == v && (cycle = true; break) + t = tarc[z] + z = t == 0 ? 0 : csrc[t] + end + isempty(stack) && break + y = pop!(stack) + y == v || push!(sub, y) + y == u && (cycle = true; break) + c = S.fchild[y] + while c != 0 + push!(stack, c) + c = S.nsib[c] + end 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)]] + if !cycle + # `e` replaces the tree arc into `v`: the subtree of `v` moves under + # `u` and becomes deeper. + Δa = a[u] + 1 - a[v] + Δb = b[u] + cw[e] - b[v] + _child_remove!(S, csrc[tarc[v]], v) + _child_append!(S, u, v) + tarc[v] = e + push!(sub, v) + _shift_subtree!(S, sub, Δa, Δb, stamp, st += 1) + else + # Contract the cycle `v → ⋯ → u → v` into `v`. + empty!(path) + z = u + while z != v + push!(path, z) + z = csrc[tarc[z]] + end + for z in path + uf[z] = v + uΔ[z] = (a[z] - a[v]) * λ - (b[z] - b[v]) + _child_remove!(S, csrc[tarc[z]], z) + end + empty!(touched) + for z in path + δ = uΔ[z] + s = S.lhead[z] + while s != 0 + f, _ = _slot_arc(S, s) + if s > ne + cdst[f] = v + cw[f] += δ + else + csrc[f] = v + cw[f] -= δ + end + push!(touched, f) + s = S.lnext[s] + end + _arcs_concat!(S, v, z) + # The children of `z` move under `v`, which is shallower, and + # their subtrees rise by `a[z] - a[v]`. + Δa = a[v] - a[z] + _shift_subtree!(S, _subtree!(empty!(sub), S, z, stack), Δa, Δa * λ, stamp, st += 1) + _child_concat!(S, v, z) + end + for f in touched + S.alive[f] || continue + if csrc[f] == cdst[f] + S.alive[f] = false + _arcs_unlink!(S, v, f) + _arcs_unlink!(S, v, f + ne) + _heap_remove!(heap, f) + else + _key!(S, f) + end + end + nlive -= length(path) end - for (k, r) in enumerate(reps) - pol[r] = cpol[k] + end + # Sum `uΔ` along each vertex's chain of contractions, compressing the + # chains as they are summed. + for z in 1:nv + empty!(stack) + y = z + while uf[y] != y + push!(stack, y) + y = uf[y] end - root = reps[cyc[1]] - for k in cyc - parent[reps[k]] = root + r = y + acc = zero(T) + for i in length(stack):-1:1 + y = stack[i] + acc += uΔ[y] + uΔ[y] = acc + uf[y] = r end - pol[root] = 0 - nnodes -= length(cyc) - 1 + x[z] = acc end x .-= x[1] return x diff --git a/test/maxbalance.jl b/test/maxbalance.jl index 7a79d47..53c4569 100644 --- a/test/maxbalance.jl +++ b/test/maxbalance.jl @@ -3,7 +3,7 @@ # `minimax_ratio`, `offtransversal_ratios`, and `check_balanced_cover` from # balanced.jl. -using MatrixCovers: _cover_transversal_balanced, _transversal_system, _howard!, _out_arcs, TIGHT_ULPS +using MatrixCovers: _cover_transversal_balanced, _transversal_system, _maxbalance # Whether the square `A`, whose diagonal is nonzero, is fully indecomposable: # whether every row reaches every other through the nonzero entries. @@ -53,18 +53,50 @@ function random_indecomposable(rng, n) 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)) +# Whether the shift `x` max-balances the strongly connected digraph with arcs +# `src[e] → dst[e]` of weight `w[e]`: for every partition of its vertices into +# two nonempty sets, the largest shifted weight leaving one set equals the +# largest entering it. +function maxbalanced(nv, src, dst, w, x; tol=1e-10) + wx = [w[e] - x[src[e]] + x[dst[e]] for e in eachindex(w)] + for mask in 1:(2^nv-2) + inset(v) = isodd(mask >> (v - 1)) + out = maximum((wx[e] for e in eachindex(w) if inset(src[e]) && !inset(dst[e])); init=-Inf) + inn = maximum((wx[e] for e in eachindex(w) if !inset(src[e]) && inset(dst[e])); init=-Inf) + abs(out - inn) <= tol || return false + end + return true +end + +@testset "_maxbalance" begin + # Random strongly connected digraphs with parallel arcs, self loops, and + # (on even trials) many tied weights. + rng = StableRNG(3) + for trial in 1:200 + nv = rand(rng, 2:9) + perm = randperm(rng, nv) + extra = rand(rng, 0:3nv) + src = vcat(perm, rand(rng, 1:nv, extra)) + dst = vcat(perm[mod1.(2:nv+1, nv)], rand(rng, 1:nv, extra)) + w = isodd(trial) ? randn(rng, length(src)) : Float64.(rand(rng, -2:2, length(src))) + x = _maxbalance(nv, src, dst, w) + @test x[1] == 0 + @test maxbalanced(nv, src, dst, w, x) + end + # A banded digraph, whose shortest-path tree is deep and moves as a whole at + # every contraction. It has too many vertices to check every partition, so + # only the cuts around single vertices are checked. + nv = 300 + src = vcat(1:nv-1, 2:nv, 1:nv-2) + dst = vcat(2:nv, 1:nv-1, 3:nv) + w = randn(rng, length(src)) + x = _maxbalance(nv, src, dst, w) + wx = [w[e] - x[src[e]] + x[dst[e]] for e in eachindex(w)] + for v in 1:nv + out = maximum(wx[e] for e in eachindex(w) if src[e] == v) + inn = maximum(wx[e] for e in eachindex(w) if dst[e] == v) + @test out ≈ inn atol = 1e-9 + end end @testset "cover_transversal(::PowerMean; maxbalance=true)" begin