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/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 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)]