From ca6b248e2f09fb53d8c835aa2c2f3531e312ddd5 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 24 Oct 2026 15:54:59 -0500 Subject: [PATCH 1/2] Max-balance blocks by parametric shortest paths MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the policy iteration in _maxbalance, which found one cycle per pass, with the minimum-balance algorithm of Young, Tarjan & Orlin (1991), adapted to maximum cycle means by reversing signs. The algorithm maintains a shortest-path tree for arc lengths λ - w as λ decreases, keeps the non-tree arcs in an indexed heap keyed by the λ at which each becomes tight, and contracts each zero-length cycle in place. Each time an arc becomes tight, the work is proportional to the size of the subtrees whose depths change, so the total is near-linear for most sparse blocks and at most O(nm + n² log n). _howard! and its tie tolerance are removed. The shifts agree with the previous implementation to roundoff. Assisted-by: Claude Fable 5.1 Assisted-by: Claude Opus 5.5 --- docs/src/objectives.md | 6 +- src/balanced.jl | 17 +- src/maxbalance.jl | 540 +++++++++++++++++++++++++++++------------ test/maxbalance.jl | 58 ++++- 4 files changed, 447 insertions(+), 174 deletions(-) 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/balanced.jl b/src/balanced.jl index 0820693..e36e512 100644 --- a/src/balanced.jl +++ b/src/balanced.jl @@ -93,12 +93,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 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 From 0be88eae25e5e85c0bf28617d70b598bf2f5a861 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 24 Oct 2026 16:20:55 -0500 Subject: [PATCH 2/2] Use loose CG solves for balanced Newton steps MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit _laplacian_solve! gains rtol and refine keywords. _balanced_direction! now solves its conjugate-gradient systems to a relative residual of 1e-4 and skips iterative refinement. The Newton direction does not need to be accurate: every CG iterate started from zero is a descent direction, and the reduced gradient is recomputed at each step. The defaults keep the refinement, which the AbsLog{2} minimal-cover solver needs to certify its KKT conditions. _balanced_factor! also records the order and number of nonzeros of the last pattern whose Cholesky factor exceeded the fill or flop budget. While later patterns stay more than half that size in both, it skips the CHOLMOD symbolic analysis and goes straight to CG; the fallback factorization is analyzed only if CG fails. For cover_transversal(PowerMean{2}(), A; target=1e-2) on a random sparse 16000×16000 matrix with five entries per row, the run time drops from 90 s to 11 s. Assisted-by: Claude Fable 5.1 Assisted-by: Claude Opus 5.5 --- src/active_set.jl | 37 ++++++++++++++++++++++++++++++++----- src/balanced.jl | 32 ++++++++++++++++++++++++-------- 2 files changed, 56 insertions(+), 13 deletions(-) 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 e36e512..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 @@ -115,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. @@ -668,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 @@ -687,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()