From b14b186bbd35464eeba247dc911376ba34467a54 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 25 Oct 2026 03:31:05 -0500 Subject: [PATCH 1/2] Add bulk steps to the active-set finish MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The finish of `_polish_difference_qp` added one constraint per step, so a seed far from the minimizer cost hundreds of Laplacian factorizations. Bulk steps now add, in one step, the violated constraints (most violated first) that extend the active forest. A primal-dual active-set step (activate all of them, then drop the constraints whose multipliers turn negative) is kept when it raises the dual function. Otherwise the constraints are added by moving along the segment to the tight minimizer of the enlarged forest, interpolating the multipliers linearly; this keeps stationarity and raises the dual function, and the move stops where a multiplier reaches zero. Bulk steps repeat while the number of violated constraints falls, within a `maxbulk` budget, and resume after a single-constraint step that leaves more than 8 violations. Goldfarb–Idnani steps then certify the result as before. On sparse lognormal transversal problems with poor seeds, finish steps drop from 761 to 286 (n = 3162) and from 1149 to 487 (n = 10000); on random difference problems in general position, from 970 to 23. Assisted-by: Claude Opus 5.5 Assisted-by: Claude Fable 5.1 --- src/active_set.jl | 239 +++++++++++++++++++++++++++++++++++++++++---- test/active_set.jl | 16 +++ 2 files changed, 235 insertions(+), 20 deletions(-) diff --git a/src/active_set.jl b/src/active_set.jl index 793cba0..32a98f1 100644 --- a/src/active_set.jl +++ b/src/active_set.jl @@ -10,7 +10,11 @@ # most violated constraint until none is violated; the final point satisfies the # KKT conditions to rounding error, which certifies it as the minimizer. Holding # the edges of a tree tight fixes the differences of `u` along it, so each -# linear solve is a graph Laplacian with one unknown per tree. +# linear solve is a graph Laplacian with one unknown per tree. Bulk steps add +# many constraints at once: the violated constraints, most violated first, that +# extend `W` to a larger forest. Like a single addition, each bulk step raises +# the dual function; they run while the number of violated constraints falls, +# within a fixed budget, and single additions finish the solve. # Constraint supports. Constraint `e` is a position in an `EdgeList`, or a # linear index into the matrix `C` of a difference grid; grid entries outside @@ -617,9 +621,30 @@ function _most_violated(u::Vector{T}, inW::Matrix{Bool}, s::_DiffGridQP, de::T) error("internal error: the most violated constraint was not found") end +# Constraints outside `W` with residual below `-ptol`, most violated first. +function _violated(u::Vector{T}, inW, s::EdgeList, ptol::T) where {T} + viol = [e for e in 1:_nconstraints(s) if !inW[e] && _resid(u, s, e) < -ptol] + return sort!(viol; by=e -> _resid(u, s, e)) +end + +function _violated(u::Vector{T}, inW::Matrix{Bool}, s::_DiffGridQP, ptol::T) where {T} + C = s.C + L = LinearIndices(C) + ur = u[s.rowidx] + viol = Int[] + for j in axes(C, 2) + uj = u[s.colidx[j]] + for i in eachindex(ur, view(C, :, j)) + # Entries outside the support have residual `+Inf`. + !inW[i, j] && ur[i] - uj - C[i, j] < -ptol && push!(viol, L[i, j]) + end + end + return sort!(viol; by=e -> _resid(u, s, e)) +end + """ - u, certified, nsteps, kkt = _polish_difference_qp(supp, u0; τ, maxsteps, seed, budgets) - u, certified, nsteps, kkt = _polish_difference_qp(edges, cvals, u0; τ, maxsteps, seed, budgets) + u, certified, nsteps, kkt = _polish_difference_qp(supp, u0; τ, maxsteps, maxbulk, seed, budgets) + u, certified, nsteps, kkt = _polish_difference_qp(edges, cvals, u0; τ, maxsteps, maxbulk, seed, budgets) Minimize `½ ∑ₑ (u[p] - u[q] - c[e])²` subject to every term being nonnegative, starting from the near-optimal `u0`. The constraints are those of `supp`, an @@ -640,6 +665,14 @@ tightest first. Constraints on no cycle of the support graph are activated ahead of the candidates: they are tight with zero multiplier at the minimizer, so a solve error must not let them appear violated. +From the seed, each bulk step adds the violated constraints, most violated +first, that extend the active forest without closing a cycle, and drops active +constraints whose multipliers would become negative. Bulk steps repeat while +the number of violated constraints falls, up to `maxbulk` bulk steps in all; +single-constraint steps then finish the solve. A single-constraint step that +leaves more than 8 violated constraints re-enters the bulk steps while the +budget lasts. Each Laplacian solve of a bulk step counts toward `nsteps`. + `certified` reports that the returned `u` satisfies the KKT conditions to tolerances proportional to the largest magnitude in `c` and `u0`; when it is false, `u` is `u0`. @@ -660,7 +693,7 @@ function _polish_difference_qp(edges::Vector{Tuple{Int,Int}}, cvals::Vector{T}, end function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{T}; - τ::Real=T(1e-4), maxsteps::Int=max(1000, 2 * length(u0)), + τ::Real=T(1e-4), maxsteps::Int=max(1000, 2 * length(u0)), maxbulk::Int=20, seed::Union{Nothing,AbstractVector}=nothing, budgets::NamedTuple=(;)) where {T} s = _qp_support(supp) @@ -670,6 +703,8 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ dtol = ptol * max(1, _nsupported(s)) # Seed: Kruskal's algorithm on the candidate constraints, in the order the # docstring describes. + inW = _active_flags(s) + W = Int[] # the active constraints, in increasing order if seed === nothing cand = _seed_candidates(u0, s, max(T(τ), ptol)) sort!(cand; by=e -> _resid(u0, s, e)) @@ -682,20 +717,7 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ sort!(cand; by=e -> (-seed[e], _resid(u0, s, e))) end # Constraints on no cycle come first; Kruskal then skips them if they recur. - cand = vcat(_acyclic_constraints(s, nV), cand) - uf = collect(1:nV) - find(x) = (while uf[x] != x; uf[x] = uf[uf[x]]; x = uf[x]; end; x) - inW = _active_flags(s) - W = Int[] # the active constraints, in increasing order - for e in cand - p, q = _edge(s, e) - rp, rq = find(p), find(q) - rp == rq && continue - uf[rp] = rq - inW[e] = true - push!(W, e) - end - sort!(W) + _activate!(inW, W, _forest_extension(W, vcat(_acyclic_constraints(s, nV), cand), s, nV)) # Multipliers and their rates; only the entries in `W` are meaningful. E = _nconstraints(s) λ = zeros(T, E) @@ -705,12 +727,24 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ u, nsteps, ptol_eff, dtol_eff = _dual_feasible!(inW, W, λ, s, nV, ptol, nsteps, maxsteps; budgets) certified = false nokkt = (; primal=T(NaN), ptol, dual=T(NaN), dtol) + nbulk = 0 # bulk steps taken + bulkon = maxbulk > 0 + bulkmin = 0 # bulk steps run while more than this many constraints are violated while nsteps <= maxsteps - e, de = _most_violated(u, inW, s, -ptol_eff) + if bulkon + u, nsteps, ptol_eff, dtol_eff, nbulk, e, bulkon = + _bulk_steps!(inW, W, λ, s, nV, u, ptol, ptol_eff, dtol_eff, nsteps, maxsteps, nbulk, maxbulk, bulkmin; budgets) + else + e, _ = _most_violated(u, inW, s, -ptol_eff) + end + e < 0 && break if e == 0 certified = all(k -> λ[k] >= -dtol_eff, W) break end + # Single-constraint step; bulk steps resume only when it leaves more + # than `bulkmin` constraints violated. + bulkmin = 8 p, q = _edge(s, e) λe = zero(T) added = false @@ -773,7 +807,8 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ ustar = u kkt = (; primal=max(zero(T), -minimum(e -> _resid(ustar, s, e), 1:_nconstraints(s))), ptol=ptol_eff, dual=max(zero(T), -minimum(k -> λ[k], W; init=zero(T))), dtol=dtol_eff) - uf .= 1:nV + uf = collect(1:nV) + find(x) = (while uf[x] != x; uf[x] = uf[uf[x]]; x = uf[x]; end; x) _foreach_constraint(s) do p, q rp, rq = find(p), find(q) rp == rq || (uf[max(rp, rq)] = min(rp, rq)) @@ -787,6 +822,170 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ return [ustar[v] + shift[find(v)] for v in 1:nV], true, nsteps, kkt end +# Objective `½ ∑ₑ dₑ(u)²`. +_objective(u::Vector{T}, s::EdgeList) where {T} = sum(e -> _resid(u, s, e)^2, 1:_nconstraints(s); init=zero(T)) / 2 +function _objective(u::Vector{T}, s::_DiffGridQP) where {T} + C = s.C + ur = u[s.rowidx] + f = zero(T) + for j in axes(C, 2) + uj = u[s.colidx[j]] + for i in eachindex(ur, view(C, :, j)) + c = C[i, j] + isfinite(c) && (f += (ur[i] - uj - c)^2) + end + end + return f / 2 +end + +# Bulk steps from the tight minimizer `u` of the dual-feasible forest `W`, run +# while more than `minviol` constraints are violated and their number falls. +# Each step adds `V`, the violated constraints, most violated first, that +# extend `W` to a larger forest. It first tries a primal-dual active-set step: +# activate all of `V`, then drop constraints with negative multipliers +# (`_dual_feasible!`). At the tight minimizer of a dual-feasible forest the +# objective equals the dual function, so the step is kept only if it raises +# the objective, as every single-constraint step does. Otherwise it is undone +# and `V` is added by `_bulk_ascent!`, which raises the dual function by +# construction, and the remaining steps of this call use `_bulk_ascent!` +# directly. Returns the updated `u`, `nsteps`, tolerances and `nbulk`; the most +# violated constraint (0 when none is violated, -1 when `maxsteps` interrupts +# `_bulk_ascent!`); and whether further bulk steps may run. +function _bulk_steps!(inW, W::Vector{Int}, λ::Vector{T}, s, nV::Int, u::Vector{T}, ptol, ptol_eff, dtol_eff, + nsteps::Int, maxsteps::Int, nbulk::Int, maxbulk::Int, minviol::Int; budgets) where {T} + viol = _violated(u, inW, s, ptol_eff) + f = _objective(u, s) + pdas = true + while length(viol) > minviol && nbulk < maxbulk && nsteps < maxsteps + V = _forest_extension(W, viol, s, nV) + isempty(V) && break + nbulk += 1 + if pdas + Wold, λold = copy(W), λ[W] + _activate!(inW, W, V) + nsteps += 1 + u1, nsteps, ptol1, dtol1 = _dual_feasible!(inW, W, λ, s, nV, ptol, nsteps, maxsteps; budgets) + f1 = _objective(u1, s) + if f1 > f + u, f, ptol_eff, dtol_eff = u1, f1, ptol1, dtol1 + else + for e in W + inW[e] = false + end + copyto!(resize!(W, length(Wold)), Wold) + for (e, λe) in zip(Wold, λold) + inW[e] = true + λ[e] = λe + end + pdas = false + end + end + if !pdas + u, nsteps, ptol_eff, dtol_eff, done = _bulk_ascent!(inW, W, λ, V, s, nV, u, ptol, ptol_eff, dtol_eff, + nsteps, maxsteps; budgets) + done || return u, nsteps, ptol_eff, dtol_eff, nbulk, -1, false + f = _objective(u, s) + end + nprev = length(viol) + viol = _violated(u, inW, s, ptol_eff) + length(viol) < nprev || break + end + return u, nsteps, ptol_eff, dtol_eff, nbulk, isempty(viol) ? 0 : first(viol), nbulk < maxbulk +end + +# Starting from the tight minimizer `u` of the dual-feasible forest `W`, +# activate the constraints `V`, which extend `W`, keeping the multipliers +# nonnegative. The segment from `u` to the tight minimizer of `W ∪ V`, with the multipliers +# interpolated linearly, satisfies stationarity throughout and raises the dual +# function. Each pass solves for the far end and moves along the segment until +# a multiplier reaches zero; that constraint leaves `W` or `V`. Constraints +# that leave `V` remain violated and are left for later steps. Returns `u`, +# `nsteps`, the tolerances (as `_dual_feasible!` does), and whether the +# remaining `V` became active; when `maxsteps` is reached first, `u` is a point +# on a segment, not the tight minimizer of `W`. +function _bulk_ascent!(inW, W::Vector{Int}, λ::Vector{T}, V::Vector{Int}, s, nV::Int, u::Vector{T}, ptol, + ptol_eff, dtol_eff, nsteps::Int, maxsteps::Int; budgets) where {T} + E = max(1, _nsupported(s)) + λ1 = zeros(T, length(λ)) + inV = Set(V) + for e in V + λ[e] = zero(T) + end + while nsteps < maxsteps + nsteps += 1 + A = sort!(vcat(W, V)) + F = _Forest(nV, s, A) + K = _TreeLaplacian(F, s; budgets) + u1, δ = _tight_minimizer(F, K, s) + _route!(λ1, _difference_gradient(u1, s), F, s) + t, kb = one(T), 0 + for k in A + λk = max(λ[k], zero(T)) + if λ1[k] < λk + tk = λk / (λk - λ1[k]) + if tk < t + t, kb = tk, k + end + end + end + if kb == 0 + _activate!(inW, W, V) + for k in W + λ[k] = λ1[k] + end + ptol_eff = max(ptol, T(2δ)) + return u1, nsteps, ptol_eff, ptol_eff * E, true + end + if t > 0 + u = u .+ t .* (u1 .- u) + for k in A + λ[k] += t * (λ1[k] - λ[k]) + end + blocked = [kb] + else + # Every constraint at zero multiplier whose multiplier would + # decrease blocks the step; none of them moves `u`. + blocked = [k for k in A if max(λ[k], zero(T)) == 0 && λ1[k] < 0] + end + for k in blocked + λ[k] = zero(T) + k in inV ? delete!(inV, k) : (inW[k] = false) + end + filter!(in(inV), V) + filter!(k -> inW[k], W) + end + return u, nsteps, ptol_eff, dtol_eff, false +end + +# Add the constraints `V` to the active set `W`, keeping `W` sorted. +function _activate!(inW, W::Vector{Int}, V) + for e in V + inW[e] = true + end + append!(W, V) + return sort!(W) +end + +# The constraints of `cand`, in order, that extend the forest `W` without +# closing a cycle with `W` or with those before them. +function _forest_extension(W::Vector{Int}, cand, s, nV::Int) + uf = collect(1:nV) + find(x) = (while uf[x] != x; uf[x] = uf[uf[x]]; x = uf[x]; end; x) + for e in W + p, q = _edge(s, e) + uf[find(p)] = find(q) + end + ext = Int[] + for e in cand + p, q = _edge(s, e) + rp, rq = find(p), find(q) + rp == rq && continue + uf[rp] = rq + push!(ext, e) + end + return ext +end + function _deactivate!(inW, W::Vector{Int}, λ, w::Int) λ[w] = zero(eltype(λ)) inW[w] = false diff --git a/test/active_set.jl b/test/active_set.jl index 339f673..43ed5de 100644 --- a/test/active_set.jl +++ b/test/active_set.jl @@ -82,6 +82,22 @@ end @test d1 ≈ d2 atol=1e-10 @test d1 ≈ d3 atol=1e-10 end + # A poor start on a larger problem, with costs in general position so that + # the minimizer is not degenerate: the bulk steps bring the step count far + # below the number of nodes, and the minimizer matches the one reached from + # a feasible start. + nV = 2000 + edges = unique!([(p, q) for (p, q) in zip(rand(rng, 1:nV, 6nV), rand(rng, 1:nV, 6nV)) if p != q]) + φ = randn(rng, nV) + cvals = [φ[p] - φ[q] - abs(randn(rng)) for (p, q) in edges] + u1, c1, _ = _polish_difference_qp(edges, cvals, φ) + u2, c2, nsteps2, _ = _polish_difference_qp(edges, cvals, randn(rng, nV); τ=0.0) + @test c1 && c2 + @test nsteps2 < 100 + d1 = [u1[p] - u1[q] - cvals[e] for (e, (p, q)) in enumerate(edges)] + d2 = [u2[p] - u2[q] - cvals[e] for (e, (p, q)) in enumerate(edges)] + @test all(>=(-1e-12), d1) + @test d1 ≈ d2 atol=1e-10 @test_throws "seed must have axes 1:3" _polish_difference_qp([(1, 2), (2, 3), (1, 3)], [0.0, 0.0, 0.0], zeros(3); seed=ones(2)) # Infeasible constraints (a cycle with positive total cost) are not certified. edges = [(1, 2), (2, 3), (3, 1)] From 8e7b6baaabddeee455468bc1ba3bc3a9bc4a7235 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 25 Oct 2026 03:31:05 -0500 Subject: [PATCH 2/2] Solve AL stages inexactly and stop at cap stalls Each stage of `_abslog2_auglag` runs LSQR with tolerance `LSQR_FORCING` (1e-3) times the KKT residual of the previous update, capped at `LSQR_ATOL_MAX` (1e-6) and floored at the full-accuracy tolerance `_lsqr_atol`; the unweighted start solve uses the full-accuracy tolerance. The Newton steps of the multiplier iteration need no more accuracy than this: on sparse lognormal matrices (full diagonal plus five lognormal entries per row) at n = 31623, LSQR iterations per solve fall from about 35 to 8 while the number of solves rises from about 245 to 290. When the active-set finish follows (`zeroslack=false`), two consecutive updates at the penalty cap that contract the KKT residual by less than a factor of `AL_CAP_CONTRACTION` (3) end the iteration; a stage that ran out of inner iterations is skipped. The finish resolves the remaining constraints in fewer solves than further updates would. `cover_transversal` on three such matrices at n = 31623 drops from 12.5/11.3/11.2 s to 5.6/5.6/6.1 s, and at n = 10000 from 1.5 to 0.95 s; `cover_min` at n = 31623 drops from 2.7 to 2.1 s. All runs converge and certify. Assisted-by: Claude Fable 5.1 --- src/minimal_covers.jl | 36 ++++++++++++++++++++++++++++++++---- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 7b8c9dd..f95e6d7 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -423,11 +423,14 @@ function _pcg!(Bmul!, x, dg, f, r, z, d, Ad, maxiter::Int, tol) return iters, fresh && nrm <= tol end +# Full-accuracy stopping tolerance of `_lsqr`. +_lsqr_atol(::Type{T}) where {T} = 5000 * eps(T) + # Matrix-free LSQR for `min norm(M*x-b)`, warm-started from `x0`. The stopping # test uses the normal-equations residual estimated from the bidiagonalization. # Returns `(x, iters)`. function _lsqr(Amul!, Atmul!, b::AbstractVector{T}, x0::AbstractVector{T}; - atol=5000 * eps(T), maxiter::Int=2 * (length(b) + length(x0)) + 100) where {T} + atol=_lsqr_atol(T), maxiter::Int=2 * (length(b) + length(x0)) + 100) where {T} x = copy(x0) u = similar(b) Amul!(u, x) @@ -1217,6 +1220,10 @@ end const AL_PENALTY = 1e2 const AL_MAXOUTER = 32 const AL_MAXITER = 40 +# At the penalty cap the multiplier iteration contracts at a fixed rate. When +# the active-set finish follows, two consecutive updates that contract the KKT +# residual by less than this factor end the iteration. +const AL_CAP_CONTRACTION = 3 # Warn when the returned scales may not minimize the objective. With # `maxouter = 0` there is no iteration and no warning. @@ -1293,6 +1300,13 @@ const LSQR_FILL_BUDGET = 1 << 30 # LSQR iterations of spanning-forest preconditioning. const LSQR_FLOP_BUDGET = 8.0e3 +# The Newton steps of the multiplier iteration need only inexact solves: each +# stage runs LSQR with tolerance `LSQR_FORCING` times the KKT residual of the previous +# update, capped at `LSQR_ATOL_MAX` and floored at the full-accuracy tolerance +# `_lsqr_atol`. +const LSQR_ATOL_MAX = 1e-6 +const LSQR_FORCING = 1e-3 + # Return the symbolic factorization of `M`, its predicted number of values, # and its predicted numeric-factorization flop count. function _precond_analysis(M::SparseMatrixCSC{Float64,Int}) @@ -1539,6 +1553,8 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; DF = use_cg_diff ? SparseCholesky() : nothing nsolves = Ref(0) nlsqr = Ref(0) + # LSQR tolerance: full accuracy for the unweighted start, then per stage. + lsqratol = Ref(_lsqr_atol(T)) lsqrtrace = Int[] # LSQR iterations of each solve ncg = Ref(0) nchol = Ref(0) @@ -1693,7 +1709,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; end mul!(px, Msp, x) solve_ptl!(px, MF, px) - soly, it = _lsqr(Pmul!, Ptmul!, cv, px) + soly, it = _lsqr(Pmul!, Ptmul!, cv, px; atol=lsqratol[]) nlsqr[] += it push!(lsqrtrace, it) return solve_up!(soly, MF, soly) @@ -1752,7 +1768,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; return z end lt_mul!(px, HFc, x) - soly, it = _lsqr(Fmul!, Ftmul!, cv, px) + soly, it = _lsqr(Fmul!, Ftmul!, cv, px; atol=lsqratol[]) nlsqr[] += it push!(lsqrtrace, it) return solve_up!(soly, HFc, soly) @@ -1774,7 +1790,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; @. z += v0 * y[g] return z end - sol, it = _lsqr(Amul!, Atmul!, cv, x) + sol, it = _lsqr(Amul!, Atmul!, cv, x; atol=lsqratol[]) nlsqr[] += it push!(lsqrtrace, it) return sol @@ -1835,7 +1851,9 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; vbest = T(Inf) # residual that the next tenfold contraction is measured from nzeroing = 0 # direct zeroings of slack multipliers since `vbest` was set nzeroed = 0 + nslow = 0 # consecutive slow updates at the penalty cap for _ in 1:maxouter + lsqratol[] = clamp(T(LSQR_FORCING) * vprev, _lsqr_atol(T), T(LSQR_ATOL_MAX)) sscale = inv(2 * κcur) bscale = inv(2 * (κcur - oneunit(T))) fcur = _fal(x, κcur, λ, sscale, bscale, supp, symmetric) @@ -1880,6 +1898,16 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; converged = true break end + # At the penalty cap the contraction rate is fixed, so when the + # active-set finish follows, two consecutive slow updates end the + # iteration: the finish resolves the remaining constraints in fewer + # solves than further updates would. A stage that ran out of inner + # iterations is skipped: its update does not reflect the contraction + # of a minimized stage. + if !zeroslack && κcur >= κcap && exit != :maxiter + nslow = v * AL_CAP_CONTRACTION > vprev ? nslow + 1 : 0 + nslow >= 2 && break + end κnext = v > vprev / 10 ? 10 * κcur : κcur if 10 * v <= vbest vbest = v