From 79cecbff10833f3c020bca0279c366effee59322 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 10 Oct 2026 02:58:37 -0500 Subject: [PATCH 1/3] Keep the active-set finish from cycling on CG noise MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Since #83 the finish solves over-budget tree Laplacians by conjugate gradients. CG leaves a forward error of its normwise backward error times the condition number: on a sparse nonsymmetric lognormal matrix with n = 31623 the per-tree values were off by 1.5e-11 while the certificate's ptol is 1e-12. Constraints on no cycle of the support (single-entry rows and columns, chains hanging off the rest) are tight with zero multiplier at the minimizer, so the noise made a few of them look violated; adding one gave a step of order 1e-12 and a zero multiplier, the ratio test then blocked at zero, and the forest shuffled among the same handful of constraints until the step limit of 4n, about 5 h per solve. The paper's timing sweep stalled on two such matrices. Two changes, both in the finish: - `_acyclic_constraints` peels degree-one vertices of the support, on the edge-list and dense-grid layouts, and the seed activates those constraints first. They lie on no cycle, their residuals are exactly zero by construction of the tree offsets, and their routed multipliers and rates are sums of exact zeros, so they are never targeted, dropped, or blocking. - `_laplacian_solve!` on `_LaplacianCG` refines the CG solution with the residual computed in double-double arithmetic until the correction is below 8 eps ‖x‖, and records the last correction as the solve's error estimate. `_dual_feasible!` raises the certificate's tolerances to twice that estimate, so `certified` never claims more accuracy than the solve delivers. On the stalling matrix the first correction equals the measured error and the second is 3e-16; the refined CG and Cholesky solutions then agree bit for bit. All three n = 31623 matrices now certify in 0 steps in about 2 s under the CG path. Finishes that took steps run about 3x faster; finishes that already took 0 steps pay 3-8% for the extra CG iterations. Assisted-by: Claude Fable 5.1 --- src/active_set.jl | 237 ++++++++++++++++++++++++++++++++++++--------- test/active_set.jl | 61 ++++++++++++ 2 files changed, 253 insertions(+), 45 deletions(-) diff --git a/src/active_set.jl b/src/active_set.jl index e829ee9..ae9ab31 100644 --- a/src/active_set.jl +++ b/src/active_set.jl @@ -346,11 +346,11 @@ function _tree_laplacian_triplets(F::_Forest{T}, supp) where {T} return idx, sparse(I, J, V, nfree, nfree) end -# Jacobi-preconditioned CG on a pinned tree Laplacian whose sparse Cholesky -# factor would exceed the fill or flop budget. Such Laplacians arise from -# expander-like supports, which fill in under elimination but are well -# conditioned, so CG converges in few iterations. A solve that CG does not -# finish falls back to the factorization. +# Jacobi-preconditioned CG, with iterative refinement, on a pinned tree +# Laplacian whose sparse Cholesky factor would exceed the fill or flop budget. +# Such Laplacians arise from expander-like supports, which fill in under +# elimination but have clustered spectra, so CG converges in few iterations. A +# solve that CG does not finish falls back to the factorization. struct _LaplacianCG L::SparseMatrixCSC{Float64,Int} # upper triangle, as analyzed S::SparseMatrixCSC{Float64,Int} # both triangles, for products @@ -358,8 +358,15 @@ struct _LaplacianCG F::SparseCholesky factored::Base.RefValue{Bool} maxiter::Int # CG iterations per pass + err::Base.RefValue{Float64} # estimated 2-norm forward error of the last solve end +# Estimated 2-norm forward error of the last `_laplacian_solve!`. A +# factorization is taken as exact at the certification floor of the finish. +_solve_error(::SparseCholesky) = 0.0 +_solve_error(::Cholesky) = 0.0 +_solve_error(C::_LaplacianCG) = C.err[] + function _laplacian_factor(L::SparseMatrixCSC{Float64,Int}; fillbudget::Real=LSQR_FILL_BUDGET, flopbudget::Real=LSQR_FLOP_BUDGET, cgmaxiter::Int=10 * size(L, 1) + 100) F = SparseCholesky() @@ -367,7 +374,7 @@ function _laplacian_factor(L::SparseMatrixCSC{Float64,Int}; fillbudget::Real=LSQ analyze!(F, L) if sizeof(Float64) * factor_entries(F) > fillbudget || factor_flops(F) > flopbudget * nnz(L) S = L + triu(L, 1)' - return _LaplacianCG(L, S, Vector(diag(S)), F, Ref(false), cgmaxiter) + return _LaplacianCG(L, S, Vector(diag(S)), F, Ref(false), cgmaxiter, Ref(0.0)) end factorize!(F, L) return F @@ -384,24 +391,69 @@ function _laplacian_solve!(x::Vector{Float64}, C::_LaplacianCG, b::Vector{Float6 if !C.factored[] n = length(b) fill!(x, 0.0) - work = ntuple(_ -> zeros(n), 4) - # The multipliers are routed sums of the solution along the trees, so - # certification needs a normwise backward error of a few units of - # roundoff, as from a factorization. `2 max(dg)` bounds `‖S‖₂` for a - # Laplacian; the first pass estimates `‖x‖`. + rdd, r, z, d, Ad, e = ntuple(_ -> zeros(n), 6) Smul! = (y, v) -> mul!(y, C.S, v) - _pcg!(Smul!, x, C.dg, b, work..., C.maxiter, sqrt(eps(Float64)) * norm(b)) - tol = 8 * eps(Float64) * (2 * maximum(C.dg) * norm(x) + norm(b)) - _, ok = _pcg!(Smul!, x, C.dg, b, work..., C.maxiter, tol) - ok && return x + η = sqrt(eps(Float64)) + _, ok = _pcg!(Smul!, x, C.dg, b, r, z, d, Ad, C.maxiter, η * norm(b)) + # Iterative refinement with the residual in double-double arithmetic. + # CG leaves a forward error of its normwise backward error times the + # condition number, and the multipliers, routed sums of the solution + # along the trees, need the solution to the certification floor. A + # correction solved to relative residual `η` shrinks the error by a + # factor of about `κ η`, so the size of the last correction estimates + # the error that remains. Refinement that stops contracting hands the + # system to the factorization. + δprev = Inf + while ok + _residual_dd!(rdd, b, C.S, x) + fill!(e, 0.0) + _, ok = _pcg!(Smul!, e, C.dg, rdd, r, z, d, Ad, C.maxiter, η * norm(rdd)) + ok || break + x .+= e + δ = norm(e) + if δ <= 8 * eps(Float64) * norm(x) + C.err[] = δ + return x + end + δ <= δprev / 100 || break + δprev = δ + end factorize!(C.F, C.L) C.factored[] = true end + C.err[] = 0.0 return solve!(x, C.F, CHOLMOD_A, b) end -# Per-tree values `y` minimizing `½‖L y - rhs‖` with pinned trees at zero; `rhs` -# is indexed by tree. +# `r = b - S x` with every product and sum carried in double-double arithmetic, +# so that `r` is accurate to a few units of roundoff in its own magnitude +# rather than in that of `b` and `S x`. +function _residual_dd!(r::Vector{Float64}, b::Vector{Float64}, S::SparseMatrixCSC{Float64,Int}, x::Vector{Float64}) + axes(r) == axes(b) == axes(x) == (axes(S, 1),) == (axes(S, 2),) || + throw(DimensionMismatch("residual, right-hand side, matrix, and solution must share axes")) + copyto!(r, b) + lo = zeros(length(r)) + rv = rowvals(S) + nz = nonzeros(S) + for j in axes(S, 2) + xj = x[j] + for k in nzrange(S, j) + i = rv[k] + p = nz[k] * xj + pe = fma(nz[k], xj, -p) # `p + pe` is the exact product + hi = r[i] + s = hi - p + t = s - hi + lo[i] += (hi - (s - t)) + (-p - t) - pe + r[i] = s + end + end + r .+= lo + return r +end + +# Per-tree values `y` minimizing `½‖L y - rhs‖` with pinned trees at zero, and +# the estimated 2-norm error of the solve; `rhs` is indexed by tree. function _tree_solve(K::_TreeLaplacian, rhs::Vector{T}) where {T} b = zeros(T, K.nfree) for (t, i) in enumerate(K.idx) @@ -413,10 +465,11 @@ function _tree_solve(K::_TreeLaplacian, rhs::Vector{T}) where {T} for (t, i) in enumerate(K.idx) i > 0 && (y[t] = x[i]) end - return y + return y, _solve_error(K.fac) end -# Minimizer of the objective with every edge of `F` tight. +# Minimizer of the objective with every edge of `F` tight, and the estimated +# error of the Laplacian solve behind it. function _tight_minimizer(F::_Forest{T}, K::_TreeLaplacian, supp::EdgeList) where {T} rhs = zeros(T, F.ntree) for (e, (p, q)) in enumerate(supp.edges) @@ -426,8 +479,8 @@ function _tight_minimizer(F::_Forest{T}, K::_TreeLaplacian, supp::EdgeList) wher rhs[tp] -= k rhs[tq] += k end - y = _tree_solve(K, rhs) - return [y[F.tree[v]] + F.off[v] for v in eachindex(F.tree)] + y, δ = _tree_solve(K, rhs) + return [y[F.tree[v]] + F.off[v] for v in eachindex(F.tree)], δ end function _tight_minimizer(F::_Forest{T}, K::_TreeLaplacian, s::_DiffGridQP) where {T} @@ -438,8 +491,8 @@ function _tight_minimizer(F::_Forest{T}, K::_TreeLaplacian, s::_DiffGridQP) wher for v in eachindex(F.tree, g) rhs[F.tree[v]] -= g[v] end - y = _tree_solve(K, rhs) - return [y[F.tree[v]] + F.off[v] for v in eachindex(F.tree)] + y, δ = _tree_solve(K, rhs) + return [y[F.tree[v]] + F.off[v] for v in eachindex(F.tree)], δ end # Add `∑ₑ dₑ gₑ` to `g`, where `dₑ = v[p] - v[q] - c[e]` when `withc` and @@ -544,16 +597,22 @@ vector of nonnegative weights indexed like the constraints (edge positions of an `EdgeList`, linear indices of a `DiffGrid`'s `C`), such as the multipliers of an augmented-Lagrangian iteration, the candidates are the constraints with positive weight or residual below `ptol`, largest weight first and then -tightest first. `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`. `kkt = (; primal, ptol, dual, dtol)` gives the -largest constraint violation of the returned `u` and the most negative active -multiplier (as nonnegative magnitudes), with the tolerances they were certified -against; `primal` and `dual` are `NaN` when `certified` is false. Each step -solves a graph-Laplacian system, and `maxsteps` bounds their number. `budgets`, a `NamedTuple` with fields -`fillbudget` and/or `flopbudget`, sets the limits above which these systems are -solved by conjugate gradients rather than a sparse Cholesky factorization. -The constant on each connected component of the graph is taken from `u0`. +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. `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`. +`kkt = (; primal, ptol, dual, dtol)` gives the largest constraint violation of +the returned `u` and the most negative active multiplier (as nonnegative +magnitudes), with the tolerances they were certified against; `primal` and +`dual` are `NaN` when `certified` is false. Each step solves a graph-Laplacian +system, and `maxsteps` bounds their number. `budgets`, a `NamedTuple` with +fields `fillbudget` and/or `flopbudget`, sets the limits above which these +systems are solved by conjugate gradients rather than a sparse Cholesky +factorization; the tolerances of the certificate then rise to twice the +estimated error of the solve, so that `certified` never claims more accuracy +than the solve delivers. The constant on each connected component of the graph +is taken from `u0`. """ function _polish_difference_qp(edges::Vector{Tuple{Int,Int}}, cvals::Vector{T}, u0::Vector{T}; kwargs...) where {T} return _polish_difference_qp(EdgeList{T}(edges, cvals, -oneunit(T)), u0; kwargs...) @@ -581,6 +640,8 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ _seed_candidates(u0, s, ptol)) 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) @@ -600,13 +661,13 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ r = zeros(T, E) nsteps = 0 u = u0 - u, nsteps = _dual_feasible!(inW, W, λ, s, nV, dtol, nsteps, maxsteps; budgets) + 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) while nsteps <= maxsteps - e, de = _most_violated(u, inW, s, -ptol) + e, de = _most_violated(u, inW, s, -ptol_eff) if e == 0 - certified = all(k -> λ[k] >= -dtol, W) + certified = all(k -> λ[k] >= -dtol_eff, W) break end p, q = _edge(s, e) @@ -635,7 +696,7 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ rhs = zeros(T, F.ntree) rhs[F.tree[p]] += one(T) rhs[F.tree[q]] -= one(T) - y = _tree_solve(K, rhs) + y, _ = _tree_solve(K, rhs) z = [y[F.tree[v]] for v in 1:nV] a = z[p] - z[q] a > 0 || return u0, false, nsteps, nokkt @@ -663,14 +724,14 @@ function _polish_difference_qp(supp::Union{EdgeList{T},DiffGrid{T}}, u0::Vector{ end added || break # Recompute from `W` alone, so rounding does not accumulate. - u, nsteps = _dual_feasible!(inW, W, λ, s, nV, dtol, nsteps, maxsteps; budgets) + u, nsteps, ptol_eff, dtol_eff = _dual_feasible!(inW, W, λ, s, nV, ptol, nsteps, maxsteps; budgets) end certified || return u0, false, nsteps, nokkt # Restore each component's constant from `u0`. (`u` is reassigned above, so # capturing it in a closure would box it; `ustar` is bound once.) ustar = u - kkt = (; primal=max(zero(T), -minimum(e -> _resid(ustar, s, e), 1:_nconstraints(s))), ptol, - dual=max(zero(T), -minimum(k -> λ[k], W; init=zero(T))), dtol) + 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 _foreach_constraint(s) do p, q rp, rq = find(p), find(q) @@ -693,20 +754,28 @@ function _deactivate!(inW, W::Vector{Int}, λ, w::Int) end # Drop the edges of `W` with negative multipliers until none remain, and return -# the tight minimizer of `W` with the step count; `λ` holds the multipliers. -function _dual_feasible!(inW, W::Vector{Int}, λ::Vector{T}, supp, nV::Int, dtol, +# the tight minimizer of `W`, the step count, and the primal and dual +# tolerances that minimizer can be certified against; `λ` holds the +# multipliers. The tolerances start from `ptol` and `ptol` times the number of +# constraints and rise with the error of the Laplacian solve: an error `δ` in +# the per-tree values moves a residual by at most `2δ`, and a routed +# multiplier, a sum of residuals over a subtree, by at most `2δ` per constraint. +function _dual_feasible!(inW, W::Vector{Int}, λ::Vector{T}, supp, nV::Int, ptol, nsteps::Int, maxsteps::Int; budgets::NamedTuple=(;)) where {T} + E = max(1, _nsupported(supp)) while true F = _Forest(nV, supp, W) K = _TreeLaplacian(F, supp; budgets) - u = _tight_minimizer(F, K, supp) + u, δ = _tight_minimizer(F, K, supp) + ptol_eff = max(ptol, T(2δ)) + dtol_eff = ptol_eff * E _route!(λ, _difference_gradient(u, supp), F, supp) nW = length(W) for e in W - λ[e] < -dtol && (inW[e] = false) + λ[e] < -dtol_eff && (inW[e] = false) end filter!(e -> inW[e], W) - (length(W) < nW && nsteps < maxsteps) || return u, nsteps + (length(W) < nW && nsteps < maxsteps) || return u, nsteps, ptol_eff, dtol_eff nsteps += 1 end end @@ -820,3 +889,81 @@ function _unfold_symcover(x::Vector{T}, (u, certified, nsteps, kkt)) where {T} xnew = (u[1:n] .- u[n+1:2n]) ./ 2 return xnew, certified, nsteps, kkt end + +# Constraints on no cycle of the support graph, found by peeling vertices of +# degree one. Each is tight with zero multiplier at the minimizer: the vertices +# beyond it enter no other constraint, so their scales make it tight at no cost +# and the gradient vanishes there. Seeding them keeps the finish from chasing +# roundoff-level violations of constraints that are tight by structure. +function _acyclic_constraints(s::EdgeList, nV::Int) + E = length(s.edges) + ptr = zeros(Int, nV + 1) + for (p, q) in s.edges + ptr[p+1] += 1 + ptr[q+1] += 1 + end + ptr[1] = 1 + cumsum!(ptr, ptr) + adj = zeros(Int, ptr[end] - 1) + cursor = ptr[1:nV] + for (e, (p, q)) in enumerate(s.edges) + adj[cursor[p]] = e + cursor[p] += 1 + adj[cursor[q]] = e + cursor[q] += 1 + end + deg = [ptr[v+1] - ptr[v] for v in 1:nV] + alive = trues(E) + stack = [v for v in 1:nV if deg[v] == 1] + out = Int[] + while !isempty(stack) + v = pop!(stack) + deg[v] == 1 || continue + for k in ptr[v]:ptr[v+1]-1 + e = adj[k] + alive[e] || continue + alive[e] = false + push!(out, e) + deg[v] -= 1 + p, q = s.edges[e] + w = p == v ? q : p + deg[w] -= 1 + deg[w] == 1 && push!(stack, w) + end + end + return out +end + +function _acyclic_constraints(s::_DiffGridQP, nV::Int) + C = s.C + L = LinearIndices(C) + m, n = size(C) + degr = copy(s.nfr) + degc = copy(s.nfc) + alive = isfinite.(C) + out = Int[] + rows = [i for i in 1:m if degr[i] == 1] + cols = [j for j in 1:n if degc[j] == 1] + while !(isempty(rows) && isempty(cols)) + if !isempty(rows) + i = pop!(rows) + degr[i] == 1 || continue + j = findfirst(view(alive, i, :)) + alive[i, j] = false + push!(out, L[i, j]) + degr[i] -= 1 + degc[j] -= 1 + degc[j] == 1 && push!(cols, j) + else + j = pop!(cols) + degc[j] == 1 || continue + i = findfirst(view(alive, :, j)) + alive[i, j] = false + push!(out, L[i, j]) + degc[j] -= 1 + degr[i] -= 1 + degr[i] == 1 && push!(rows, i) + end + end + return out +end diff --git a/test/active_set.jl b/test/active_set.jl index 708de4d..c445325 100644 --- a/test/active_set.jl +++ b/test/active_set.jl @@ -162,6 +162,7 @@ end @test !Fcg.factored[] @test xcg ≈ xch rtol=1e-13 @test norm(S * xcg - b) <= 8 * eps() * (opnorm(Matrix(S)) * norm(xcg) + norm(b)) + @test Fcg.err[] <= 8 * eps() * norm(xcg) # A solve that CG does not finish falls back to the factorization, which # later solves reuse. Fcg = MatrixCovers._laplacian_factor(L; flopbudget=0, cgmaxiter=0) @@ -190,6 +191,66 @@ end @test acg .* bcg' ≈ aref .* bref' rtol=1e-10 end +@testset "refined conjugate-gradient Laplacian solves" begin + # A path Laplacian pinned at one end has condition number of order `nV²`. + # CG alone leaves a forward error of that times roundoff; refinement with a + # double-double residual brings the solve to the roundoff level of the + # solution. The right-hand side is exact, so the error is known. + rng = StableRNG(23) + function path_laplacian(nV) + S = spdiagm(-1 => fill(-1.0, nV - 1), 0 => [fill(2.0, nV - 1); 1.0], 1 => fill(-1.0, nV - 1)) + return S, sparse(triu(S)) + end + for nV in (300, 3000) + S, L = path_laplacian(nV) + xtrue = Float64.(rand(rng, -5:5, nV)) + b = S * xtrue + Fcg = MatrixCovers._laplacian_factor(L; flopbudget=0) + @test Fcg isa MatrixCovers._LaplacianCG + x = MatrixCovers._laplacian_solve!(zeros(nV), Fcg, b) + @test !Fcg.factored[] + @test maximum(abs, x - xtrue) <= 64 * eps() * norm(xtrue) + @test Fcg.err[] <= 8 * eps() * norm(x) + # The residual helper: exact at the solution, accurate for a perturbed one. + @test all(iszero, MatrixCovers._residual_dd!(zeros(nV), b, S, xtrue)) + y = xtrue .+ 1e-9 .* randn(rng, nV) + @test MatrixCovers._residual_dd!(zeros(nV), b, S, y) ≈ Vector(S * (xtrue - y)) rtol=1e-5 + end +end + +@testset "constraints on no cycle of the support" begin + # Rows and columns with a single stored entry, and chains hanging off the + # rest of the support, give constraints that are tight with zero multiplier + # at the minimizer. The finish activates them before anything else. + rng = StableRNG(29) + m, n = 14, 12 + core = [(i, m + j) for j in 1:10 for i in 1:10] + # Single-entry rows 11 and 12, a single-entry column 11, and the chain + # row 13 – column 12 – row 14 – column 5 into the core. + tail = [(11, m + 1), (12, m + 2), (3, m + 11), (13, m + 12), (14, m + 12), (14, m + 5)] + edges = vcat(core, tail) + cvals = randn(rng, length(edges)) + nV = m + n + s = MatrixCovers.EdgeList{Float64}(edges, cvals, -1.0) + acy = MatrixCovers._acyclic_constraints(s, nV) + @test sort(acy) == collect(length(core)+1:length(edges)) + C = fill(-Inf, m, n) + for (e, (p, q)) in enumerate(edges) + C[p, q-m] = cvals[e] + end + grid = MatrixCovers._qp_support(MatrixCovers.DiffGrid{Float64}(C, collect(1:m), collect(m+1:m+n))) + @test Set(C[MatrixCovers._acyclic_constraints(grid, nV)]) == Set(cvals[length(core)+1:end]) + # The finish certifies the minimizer of a cover whose support has these + # appendages, whether the Laplacians are factored or solved by CG. + A = sparse(first.(edges), last.(edges) .- m, exp.(cvals), m, n) + aj, bj = MatrixCovers.cover_min_jump(AbsLog{2}(), Matrix(A)) + for flopbudget in (Inf, 0) + a, b, st = _cover_min_abslog2(A; maxouter=2, linsolve=:lsqr, flopbudget) + @test st.polish.certified + @test a .* b' ≈ aj .* bj' rtol=1e-6 + end +end + @testset "difference-grid active-set solver" begin # The grid must reproduce the edge list of its supported entries, with rows # and columns sharing unknowns and some entries outside the support. From 17fb457812e1585502fb4058350d3ca604cee74c Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 10 Oct 2026 03:48:58 -0500 Subject: [PATCH 2/3] Precondition and refine for every element type The LSQR preconditioner was switched off for element types other than Float64 because its factorizations are CHOLMOD's. Only those are now gated on Float64; the diagonal preconditioner, weighted degrees, is generic and serves other types throughout. The finish's `_LaplacianCG` is now parametric in the element type. A sparse tree Laplacian in a type CHOLMOD does not factor was densified and factored at O(n^3); it now takes CG with compensated-residual refinement like Float64, falling back to a dense Cholesky factorization only when refinement stops contracting. The error-free transformations in the residual hold for any type with a correctly rounded `fma`. On a 3162x3162 sparse lognormal matrix in Double64, a solve whose multiplier iteration stops early drops from 204 s to 16 s, and the default solve's LSQR iterations fall from 67699 to 43901. The results agree with Float64 to 1e-15 in the log scales. Assisted-by: Claude Fable 5.1 --- src/active_set.jl | 74 +++++++++++++++++++++++------------------- src/minimal_covers.jl | 17 ++++++---- test/active_set.jl | 26 +++++++++++++-- test/element_types.jl | 5 +++ test/minimal_covers.jl | 7 ++-- 5 files changed, 83 insertions(+), 46 deletions(-) diff --git a/src/active_set.jl b/src/active_set.jl index ae9ab31..e6fe98b 100644 --- a/src/active_set.jl +++ b/src/active_set.jl @@ -223,7 +223,7 @@ function _TreeLaplacian(F::_Forest{T}, supp; budgets::NamedTuple=(;)) where {T} else idx, L = _tree_laplacian_triplets(F, supp) end - fac = L isa SparseMatrixCSC{Float64,Int} ? _laplacian_factor(L; budgets...) : _laplacian_factor(L) + fac = L isa SparseMatrixCSC ? _laplacian_factor(L; budgets...) : _laplacian_factor(L) return _TreeLaplacian{T,typeof(fac)}(idx, fac, size(L, 1)) end @@ -347,18 +347,21 @@ function _tree_laplacian_triplets(F::_Forest{T}, supp) where {T} end # Jacobi-preconditioned CG, with iterative refinement, on a pinned tree -# Laplacian whose sparse Cholesky factor would exceed the fill or flop budget. -# Such Laplacians arise from expander-like supports, which fill in under -# elimination but have clustered spectra, so CG converges in few iterations. A -# solve that CG does not finish falls back to the factorization. -struct _LaplacianCG - L::SparseMatrixCSC{Float64,Int} # upper triangle, as analyzed - S::SparseMatrixCSC{Float64,Int} # both triangles, for products - dg::Vector{Float64} - F::SparseCholesky - factored::Base.RefValue{Bool} +# Laplacian: in `Float64`, one whose sparse Cholesky factor would exceed the +# fill or flop budget; in other element types, every sparse one, since CHOLMOD +# factors only `Float64`. Such Laplacians arise from expander-like supports, +# which fill in under elimination but have clustered spectra, so CG converges +# in few iterations. A solve that refinement does not finish falls back to a +# factorization: the analyzed CHOLMOD factor `F` in `Float64`, a dense Cholesky +# factorization otherwise. +struct _LaplacianCG{T} + L::SparseMatrixCSC{T,Int} # upper triangle + S::SparseMatrixCSC{T,Int} # both triangles, for products + dg::Vector{T} + F::Union{SparseCholesky,Nothing} # analyzed, `Float64` only + fallback::Base.RefValue{Any} # the factorization once CG gave up; `nothing` before maxiter::Int # CG iterations per pass - err::Base.RefValue{Float64} # estimated 2-norm forward error of the last solve + err::Base.RefValue{T} # estimated 2-norm forward error of the last solve end # Estimated 2-norm forward error of the last `_laplacian_solve!`. A @@ -374,28 +377,33 @@ function _laplacian_factor(L::SparseMatrixCSC{Float64,Int}; fillbudget::Real=LSQ analyze!(F, L) if sizeof(Float64) * factor_entries(F) > fillbudget || factor_flops(F) > flopbudget * nnz(L) S = L + triu(L, 1)' - return _LaplacianCG(L, S, Vector(diag(S)), F, Ref(false), cgmaxiter, Ref(0.0)) + return _LaplacianCG{Float64}(L, S, Vector(diag(S)), F, Ref{Any}(nothing), cgmaxiter, Ref(0.0)) end factorize!(F, L) return F end -# Wider types: dense Cholesky of the upper triangle. -_laplacian_factor(L::SparseMatrixCSC) = _laplacian_factor(Matrix(L)) +# Element types CHOLMOD does not factor; the budgets do not apply. +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)) + S = L + triu(L, 1)' + return _LaplacianCG{T}(L, S, Vector(diag(S)), nothing, Ref{Any}(nothing), cgmaxiter, Ref(zero(T))) +end _laplacian_factor(L::Matrix) = LinearAlgebra.cholesky!(Symmetric(L, :U)) _laplacian_solve!(x::Vector{Float64}, F::SparseCholesky, b::Vector{Float64}) = isempty(b) ? x : solve!(x, F, CHOLMOD_A, b) _laplacian_solve!(x, F::Cholesky, b) = ldiv!(x, F, b) -function _laplacian_solve!(x::Vector{Float64}, C::_LaplacianCG, b::Vector{Float64}) - if !C.factored[] +function _laplacian_solve!(x::Vector{T}, C::_LaplacianCG{T}, b::Vector{T}) where {T} + if C.fallback[] === nothing n = length(b) - fill!(x, 0.0) - rdd, r, z, d, Ad, e = ntuple(_ -> zeros(n), 6) + 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(Float64)) + η = sqrt(eps(T)) _, ok = _pcg!(Smul!, x, C.dg, b, r, z, d, Ad, C.maxiter, η * norm(b)) - # Iterative refinement with the residual in double-double arithmetic. + # 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 # along the trees, need the solution to the certification floor. A @@ -403,36 +411,36 @@ function _laplacian_solve!(x::Vector{Float64}, C::_LaplacianCG, b::Vector{Float6 # factor of about `κ η`, so the size of the last correction estimates # the error that remains. Refinement that stops contracting hands the # system to the factorization. - δprev = Inf + δprev = T(Inf) while ok _residual_dd!(rdd, b, C.S, x) - fill!(e, 0.0) + fill!(e, zero(T)) _, ok = _pcg!(Smul!, e, C.dg, rdd, r, z, d, Ad, C.maxiter, η * norm(rdd)) ok || break x .+= e δ = norm(e) - if δ <= 8 * eps(Float64) * norm(x) + if δ <= 8 * eps(T) * norm(x) C.err[] = δ return x end δ <= δprev / 100 || break δprev = δ end - factorize!(C.F, C.L) - C.factored[] = true + C.fallback[] = C.F === nothing ? _laplacian_factor(Matrix(C.L)) : factorize!(C.F, C.L) end - C.err[] = 0.0 - return solve!(x, C.F, CHOLMOD_A, b) + C.err[] = zero(T) + return _laplacian_solve!(x, C.fallback[], b) end -# `r = b - S x` with every product and sum carried in double-double arithmetic, -# so that `r` is accurate to a few units of roundoff in its own magnitude -# rather than in that of `b` and `S x`. -function _residual_dd!(r::Vector{Float64}, b::Vector{Float64}, S::SparseMatrixCSC{Float64,Int}, x::Vector{Float64}) +# `r = b - S x` with every product and sum carried in twice the working +# precision by error-free transformations, so that `r` is accurate to a few +# units of roundoff in its own magnitude rather than in that of `b` and `S x`. +# `fma` must be correctly rounded in `T`. +function _residual_dd!(r::Vector{T}, b::Vector{T}, S::SparseMatrixCSC{T,Int}, x::Vector{T}) where {T} axes(r) == axes(b) == axes(x) == (axes(S, 1),) == (axes(S, 2),) || throw(DimensionMismatch("residual, right-hand side, matrix, and solution must share axes")) copyto!(r, b) - lo = zeros(length(r)) + lo = zeros(T, length(r)) rv = rowvals(S) nz = nonzeros(S) for j in axes(S, 2) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 555e02c..ce8ba70 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -55,7 +55,8 @@ number of stored entries in `A`'s support, where a symmetric off-diagonal entry counts twice, as `nnz` counts it for a fully stored matrix. Otherwise it factors a sparser matrix that keeps the couplings only among entries whose constraints are active and the diagonal contributions of the rest; whenever -that factor also exceeds the budgets, it uses a diagonal preconditioner. +that factor also exceeds the budgets, it uses a diagonal preconditioner. Other +element types use the diagonal preconditioner throughout. `fillbudget=Inf` and `flopbudget=Inf` together always select the full Cholesky preconditioner; `fillbudget=0` or `flopbudget=0` always select the diagonal one. The same budgets apply to the active-set method that finishes an @@ -1415,8 +1416,10 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; use_cg = use_woodbury || use_cg_diff ne = supp isa EdgeList ? length(supp.edges) : 0 use_lsqr = linsolve === :lsqr - # CHOLMOD preconditioning is limited to Float64. - use_precond = use_lsqr && T === Float64 + use_precond = use_lsqr + # The factor preconditioners are CHOLMOD factorizations, so other element + # types use the diagonal preconditioner throughout. + use_cholmod = T === Float64 # Residual multiplicity of each stored entry. symmetric = sys.symmetric mult = (p, q) -> (symmetric && p != q) ? 2 : 1 @@ -1477,7 +1480,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; # stages. Use its diagonal if the predicted Cholesky factor exceeds the fill # or flop budget. Msp = _precond_pattern(T, supp, v0, use_precond ? N : 0, mult) - MF, fill_entries, flops = use_precond ? _precond_analysis(Msp) : (nothing, 0, 0.0) + MF, fill_entries, flops = use_precond && use_cholmod ? _precond_analysis(Msp) : (nothing, 0, 0.0) # Number of stored entries of the user's support, counted as `nnz` counts # the input matrix: a symmetric off-diagonal entry is stored in both # triangles and so counts twice, and every other stored entry counts once. @@ -1488,7 +1491,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; end end flopcap = flopbudget * nstored # largest affordable predicted flop count - use_factor = use_precond && sizeof(T) * fill_entries <= fillbudget && + use_factor = use_precond && use_cholmod && sizeof(T) * fill_entries <= fillbudget && (flopbudget == Inf || flops <= flopcap) # Positions of the entries each factored solve overwrites: the diagonal, and # both copies of each off-diagonal support entry. @@ -1507,7 +1510,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; # contribute only to the diagonal. Its pattern follows the violation # pattern, so each new `(κ, vpat)` gets its own symbolic analysis and budget # check; `forestok` records whether the current one passed. - HF = use_precond && !use_factor ? SparseCholesky() : nothing + HF = use_precond && use_cholmod && !use_factor ? SparseCholesky() : nothing HP = Ref(spzeros(T, 0, 0)) forestok = Ref(false) hI = Int[] @@ -1686,7 +1689,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; push!(lsqrtrace, it) return solve_up!(soly, MF, soly) elseif use_precond - if prevκ[] != κl || vpat != prevpat + if HF !== nothing && (prevκ[] != κl || vpat != prevpat) # Heavy-forest matrix: the violated part of `RᵀWR`, the # diagonal of the rest, the diagonal of `v0 v0ᵀ`, and a # ridge. `sparse` sums duplicates. diff --git a/test/active_set.jl b/test/active_set.jl index c445325..339f673 100644 --- a/test/active_set.jl +++ b/test/active_set.jl @@ -159,7 +159,7 @@ end b = randn(rng, nV - 1) xch = MatrixCovers._laplacian_solve!(zeros(nV - 1), Fch, b) xcg = MatrixCovers._laplacian_solve!(zeros(nV - 1), Fcg, b) - @test !Fcg.factored[] + @test Fcg.fallback[] === nothing @test xcg ≈ xch rtol=1e-13 @test norm(S * xcg - b) <= 8 * eps() * (opnorm(Matrix(S)) * norm(xcg) + norm(b)) @test Fcg.err[] <= 8 * eps() * norm(xcg) @@ -167,7 +167,7 @@ end # later solves reuse. Fcg = MatrixCovers._laplacian_factor(L; flopbudget=0, cgmaxiter=0) xfb = MatrixCovers._laplacian_solve!(zeros(nV - 1), Fcg, b) - @test Fcg.factored[] + @test Fcg.fallback[] !== nothing @test xfb ≈ xch rtol=1e-13 @test MatrixCovers._laplacian_solve!(zeros(nV - 1), Fcg, 2b) ≈ 2xch rtol=1e-13 end @@ -208,7 +208,7 @@ end Fcg = MatrixCovers._laplacian_factor(L; flopbudget=0) @test Fcg isa MatrixCovers._LaplacianCG x = MatrixCovers._laplacian_solve!(zeros(nV), Fcg, b) - @test !Fcg.factored[] + @test Fcg.fallback[] === nothing @test maximum(abs, x - xtrue) <= 64 * eps() * norm(xtrue) @test Fcg.err[] <= 8 * eps() * norm(x) # The residual helper: exact at the solution, accurate for a perturbed one. @@ -216,6 +216,26 @@ end y = xtrue .+ 1e-9 .* randn(rng, nV) @test MatrixCovers._residual_dd!(zeros(nV), b, S, y) ≈ Vector(S * (xtrue - y)) rtol=1e-5 end + # Element types CHOLMOD does not factor take the same path, to their own + # roundoff, and fall back to a dense Cholesky factorization. + nV = 300 + S, L = path_laplacian(nV) + Sb, Lb = BigFloat.(S), BigFloat.(L) + xtrue = BigFloat.(rand(rng, -5:5, nV)) + b = Sb * xtrue + Fcg = MatrixCovers._laplacian_factor(Lb) + @test Fcg isa MatrixCovers._LaplacianCG{BigFloat} + @test Fcg.F === nothing + x = MatrixCovers._laplacian_solve!(zeros(BigFloat, nV), Fcg, b) + @test Fcg.fallback[] === nothing + @test maximum(abs, x - xtrue) <= 64 * eps(BigFloat) * norm(xtrue) + @test Fcg.err[] <= 8 * eps(BigFloat) * norm(x) + Fcg = MatrixCovers._laplacian_factor(Lb; cgmaxiter=0) + x = MatrixCovers._laplacian_solve!(zeros(BigFloat, nV), Fcg, b) + @test Fcg.fallback[] isa Cholesky + @test Fcg.err[] == 0 + @test maximum(abs, x - xtrue) <= 1e-6 * norm(xtrue) + @test MatrixCovers._laplacian_factor(spzeros(BigFloat, 0, 0)) isa Cholesky end @testset "constraints on no cycle of the support" begin diff --git a/test/element_types.jl b/test/element_types.jl index 50a256a..ee08b8b 100644 --- a/test/element_types.jl +++ b/test/element_types.jl @@ -73,6 +73,11 @@ @test abig isa Vector{BigFloat} @test Float64.(abig) ≈ first(symcover_min(AbsLog{2}(), Float64.(Abig))) rtol=1e-6 @test abig != BigFloat.(Float32.(abig)) + # LSQR on such a type is preconditioned by the diagonal, since the + # factor preconditioners are CHOLMOD's. + abig_lsqr, sbig = symcover_min(AbsLog{2}(), sparse(Abig); linsolve=:lsqr) + @test MatrixCovers.converged(sbig) && sbig.precond === :diagonal + @test Float64.(abig_lsqr) ≈ Float64.(abig) rtol=1e-6 end @testset "BigFloat flows through the family" begin diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 976ebe2..ed0c442 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -458,12 +458,13 @@ end @test iscover(gf, hf, G) && iscover(gd, hd, G) @test cover_objective(AbsLog{2}(), gf, hf, G) ≈ cover_objective(AbsLog{2}(), gd, hd, G) rtol=1e-8 - # The budget is a keyword of the public solvers, and the paths that never - # precondition say so. + # The budget is a keyword of the public solvers, the dense path never + # preconditions, and an element type CHOLMOD does not factor uses the + # diagonal. @test first(symcover_min(AbsLog{2}(), A; linsolve=:lsqr, fillbudget=0)) ≈ ad rtol=1e-6 @test MatrixCovers._symcover_min_abslog2(A; linsolve=:dense)[2].precond === :none Abig = BigFloat.([4.0 1.0 0.5; 1.0 3.0 1.0; 0.5 1.0 2.5]) - @test MatrixCovers._symcover_min_abslog2(Abig; linsolve=:lsqr)[2].precond === :none + @test MatrixCovers._symcover_min_abslog2(Abig; linsolve=:lsqr)[2].precond === :diagonal end """ From a5d8a6f939b905020975beb97a7b9c807314c848 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 10 Oct 2026 05:05:24 -0500 Subject: [PATCH 3/3] Precondition LSQR with a spanning forest in any element type The heavy-forest preconditioner of the multiplier iteration kept every coupling among violated entries and factored the result with CHOLMOD, with a symbolic analysis and a fill-budget check at each change of penalty or violation pattern, and fell back to the diagonal when that factor was over budget or the element type was not Float64. In practice the violated couplings form a forest: over eight sparse lognormal matrices of the paper's sweep, up to n = 100000, at most 45 of about two million violated couplings per solve closed a cycle, and restricting the pattern to a maximum-weight spanning forest changed the LSQR iteration totals by at most 8 %. The preconditioner is now that forest, factored by eliminating from the leaves (`tree_cholesky.jl`): no fill, no analysis, no budget, and the arithmetic of the element type. Its diagonal is every entry's weighted degree plus the ridge, so it is strictly diagonally dominant and the pivots stay positive. The diagonal regime is retired, with its `ndiagonal` statistic; `precond` is now `:factor`, `:forest`, or `:none`. Float64 solves on the sweep matrices run 20-30 % faster, the analysis per rebuild being gone. A Double64 solve on a 3162x3162 sparse matrix drops from 84 s and 43901 LSQR iterations to 12 s and 2052. Assisted-by: Claude Fable 5.1 --- src/MatrixCovers.jl | 1 + src/minimal_covers.jl | 130 ++++++++++++----------------------------- src/stats.jl | 6 +- src/tree_cholesky.jl | 125 +++++++++++++++++++++++++++++++++++++++ test/element_types.jl | 15 ++++- test/minimal_covers.jl | 66 +++++++++++---------- test/runtests.jl | 1 + test/tree_cholesky.jl | 43 ++++++++++++++ 8 files changed, 254 insertions(+), 133 deletions(-) create mode 100644 src/tree_cholesky.jl create mode 100644 test/tree_cholesky.jl diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index ca96764..12465fa 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -33,6 +33,7 @@ include("initializers.jl") # the start menu; consumed by both solver families include("soft_covers.jl") include("powermean.jl") # the PowerMean soft covers, the soft-cover default include("sparse_cholesky.jl") # CHOLMOD factorizations for the minimal-cover solvers +include("tree_cholesky.jl") # fill-free factorization of forest-patterned matrices, for preconditioning include("minimal_covers.jl") include("active_set.jl") # exact solution of the AbsLog{2} problem from a near-optimal point include("transversal.jl") # transversal-tight minimal covers diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index ce8ba70..61008ff 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -52,11 +52,11 @@ For `Float64`, `:lsqr` uses a Cholesky preconditioner when its predicted storage does not exceed `fillbudget` bytes (default `2^30`) and its predicted numeric- factorization flop count does not exceed `flopbudget` (default `8e3`) times the number of stored entries in `A`'s support, where a symmetric off-diagonal entry -counts twice, as `nnz` counts it for a fully stored matrix. Otherwise it -factors a sparser matrix that keeps the couplings only among entries whose -constraints are active and the diagonal contributions of the rest; whenever -that factor also exceeds the budgets, it uses a diagonal preconditioner. Other -element types use the diagonal preconditioner throughout. +counts twice, as `nnz` counts it for a fully stored matrix. Otherwise, and for +every other element type, it preconditions with a maximum-weight spanning +forest of the couplings among entries whose constraints are active, together +with the diagonal contributions of all entries; that factorization has no fill +and needs no budget. `fillbudget=Inf` and `flopbudget=Inf` together always select the full Cholesky preconditioner; `fillbudget=0` or `flopbudget=0` always select the diagonal one. The same budgets apply to the active-set method that finishes an @@ -1383,12 +1383,10 @@ _dense_factor_type(::Type{T}) where {T} = LinearAlgebra.LU{T,Matrix{T},Vector{In # the `linsolve` and `precond` choices, the full Cholesky preconditioner's # predicted entry count `fill_entries` and flop count `factor_flops` (`0` and # `0.0` without a factor), the number of numeric factorizations `nrefactor` -# (full or heavy-forest), the numbers of LSQR solves -# preconditioned by the heavy-forest factor (`nforest`) and by the diagonal -# (`ndiagonal`) when the full factor is over budget, and the number `nzeroed` -# of positive multipliers zeroed directly on entries with slack while the -# penalty sat at its cap. `precond` is `:factor` (full factor), `:forest` (the -# heavy-forest factor served at least one solve), `:diagonal`, or `:none`. +# (full or forest), the number of LSQR solves preconditioned by the forest +# factor (`nforest`), and the number `nzeroed` of positive multipliers zeroed +# directly on entries with slack while the penalty sat at its cap. `precond` is +# `:factor` (full factor), `:forest`, or `:none`. # With `zeroslack=false` the iteration stops at a stall at the penalty cap # instead of zeroing such multipliers, for callers that finish the solve by the # active-set method of `_polish_cover` or `_polish_symcover`. @@ -1417,8 +1415,8 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; ne = supp isa EdgeList ? length(supp.edges) : 0 use_lsqr = linsolve === :lsqr use_precond = use_lsqr - # The factor preconditioners are CHOLMOD factorizations, so other element - # types use the diagonal preconditioner throughout. + # The full-factor preconditioner is a CHOLMOD factorization, so other + # element types use the forest preconditioner throughout. use_cholmod = T === Float64 # Residual multiplicity of each stored entry. symmetric = sys.symmetric @@ -1439,7 +1437,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; cv = zeros(T, ne + 1) # √weight · log|A_ij|, with a trailing 0 gauge target # Violated entries under the current frozen weights. vpat = _violation_pattern(supp) - # The dense and factor-preconditioned solves (full or heavy-forest factor) + # The dense and factor-preconditioned solves (full or forest factor) # refactor only when the weights change, and the weights depend only on # (κ, active set): cache that key. prevκ = Ref(zero(T)) @@ -1504,18 +1502,18 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; epos[2 * e] = _nzindex(Msp, q, p) end end - psqrt = zeros(T, use_factor ? 0 : (use_precond ? N : 0)) # `K` of the diagonal preconditioner - # Heavy-forest preconditioner, used when the full factor is over budget: - # violated entries keep their off-diagonal couplings and the others - # contribute only to the diagonal. Its pattern follows the violation - # pattern, so each new `(κ, vpat)` gets its own symbolic analysis and budget - # check; `forestok` records whether the current one passed. - HF = use_precond && use_cholmod && !use_factor ? SparseCholesky() : nothing - HP = Ref(spzeros(T, 0, 0)) - forestok = Ref(false) + # Forest preconditioner, used when the full factor is over budget or the + # element type is not CHOLMOD's: a maximum-weight spanning forest of the + # couplings among violated entries, every entry's contribution to the + # diagonal, the diagonal of `v0 v0ᵀ`, and a ridge. The couplings among + # violated entries form a forest up to a few edges in practice, so little is + # lost, and the factorization has no fill, no symbolic analysis, and no + # budget. It is rebuilt for each new `(κ, vpat)`. + HF = Ref{Union{Nothing,TreeCholesky{T}}}(nothing) hI = Int[] hJ = Int[] hV = T[] + hW = T[] rhs = zeros(T, use_woodbury ? N : 0, size(U, 2) + 1) wsol = similar(rhs) WF = use_woodbury ? SparseCholesky() : nothing @@ -1534,7 +1532,6 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; nchol = Ref(0) nrefactor = Ref(0) nforest = Ref(0) - ndiagonal = Ref(0) solve_weighted = function (x, κ) nsolves[] += 1 if supp isa Grid @@ -1689,11 +1686,12 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; push!(lsqrtrace, it) return solve_up!(soly, MF, soly) elseif use_precond - if HF !== nothing && (prevκ[] != κl || vpat != prevpat) - # Heavy-forest matrix: the violated part of `RᵀWR`, the - # diagonal of the rest, the diagonal of `v0 v0ᵀ`, and a - # ridge. `sparse` sums duplicates. - empty!(hI); empty!(hJ); empty!(hV) + if prevκ[] != κl || vpat != prevpat + # Forest matrix: a maximum-weight spanning forest of the + # violated couplings of `RᵀWR`, the diagonal of all of + # `RᵀWR` and of `v0 v0ᵀ`, and a ridge. The diagonal + # dominates the kept couplings, so the pivots stay positive. + empty!(hI); empty!(hJ); empty!(hV); empty!(hW) fill!(mdiag, zero(T)) for (e, (p, q)) in enumerate(edges) w = ws[e]^2 @@ -1703,7 +1701,7 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; mdiag[p] += w mdiag[q] += w if vpat[e] - push!(hI, p, q); push!(hJ, q, p); push!(hV, qs * w, qs * w) + push!(hI, p); push!(hJ, q); push!(hV, qs * w); push!(hW, w) end end end @@ -1713,27 +1711,15 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; dmax = max(dmax, mdiag[p]) end ρ = _precond_ridge(dmax) - for p in 1:N - push!(hI, p); push!(hJ, p); push!(hV, mdiag[p] + ρ) - end - HP[] = sparse(hI, hJ, hV, N, N) - HFa = HF::SparseCholesky - analyze!(HFa, HP[]) - forestok[] = sizeof(T) * factor_entries(HFa) <= fillbudget && - (flopbudget == Inf || factor_flops(HFa) <= flopcap) - if forestok[] - factorize!(HFa, HP[]) - nrefactor[] += 1 - end + keep = max_weight_forest(N, hI, hJ, hW) + HF[] = tree_cholesky(N, hI[keep], hJ[keep], hV[keep], mdiag .+ ρ) + nrefactor[] += 1 prevκ[] = κl copyto!(prevpat, vpat) end - end - if use_precond && !use_factor && forestok[] nforest[] += 1 - HFc = HF::SparseCholesky - HPc = HP[] - # `K = P'L` from the heavy-forest factor, as in the full-factor branch. + HFc = HF[]::TreeCholesky{T} + # `K = L` from the forest factor, as `K = P'L` in the full-factor branch. Fmul! = function (y, yv) solve_up!(pxv, HFc, yv) for (e, (p, q)) in enumerate(edges) @@ -1753,55 +1739,11 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; solve_ptl!(z, HFc, pg) return z end - mul!(px, HPc, x) - solve_ptl!(px, HFc, px) + lt_mul!(px, HFc, x) soly, it = _lsqr(Fmul!, Ftmul!, cv, px) nlsqr[] += it push!(lsqrtrace, it) return solve_up!(soly, HFc, soly) - elseif use_precond && !use_factor - # The heavy-forest factor is over budget for this pattern. - ndiagonal[] += 1 - # Weighted degrees, the diagonal of `RᵀWR`. - fill!(mdiag, zero(T)) - for (e, (p, q)) in enumerate(edges) - w = ws[e]^2 - if p == q - mdiag[p] += 4 * w - else - mdiag[p] += w - mdiag[q] += w - end - end - # An unknown outside the support gets an identity row. - for p in 1:N - d = mdiag[p] + v0[p]^2 - psqrt[p] = sqrt(d > 0 ? d : oneunit(T)) - end - # Diagonal `K` needs only elementwise scaling. - Dmul! = function (y, yv) - @. px = yv / psqrt - for (e, (p, q)) in enumerate(edges) - y[e] = ws[e] * (px[p] + qs * px[q]) - end - y[g] = dot(v0, px) - return y - end - Dtmul! = function (z, y) - fill!(pg, zero(T)) - for (e, (p, q)) in enumerate(edges) - t = ws[e] * y[e] - pg[p] += t - pg[q] += qs * t - end - @. pg += v0 * y[g] - @. z = pg / psqrt - return z - end - soly, it = _lsqr(Dmul!, Dtmul!, cv, psqrt .* x) - nlsqr[] += it - push!(lsqrtrace, it) - return soly ./ psqrt end Amul! = function (y, xx) for (e, (p, q)) in enumerate(edges) @@ -1969,8 +1911,8 @@ function _abslog2_auglag(sys::SupportSystem{T}, x0; drops=Tuple(drops), kkt=Tuple(viols), κs=Tuple(κtrace), converged, tol=vtol, polish=nothing, linsolve=(use_lsqr ? :lsqr : use_cg ? :woodbury : :dense), - precond=(!use_precond ? :none : use_factor ? :factor : nforest[] > 0 ? :forest : :diagonal), - nrefactor=nrefactor[], nforest=nforest[], ndiagonal=ndiagonal[], fill_entries, factor_flops=flops, nzeroed) + precond=(!use_precond ? :none : use_factor ? :factor : :forest), + nrefactor=nrefactor[], nforest=nforest[], fill_entries, factor_flops=flops, nzeroed) end # Worker for `symcover_min(::AbsLog{2})`, returning `(a, stats)`. A supplied diff --git a/src/stats.jl b/src/stats.jl index 0f7f420..107f20f 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -102,8 +102,7 @@ otherwise `last(kkt)` and `tol`). `iterations` counts multiplier updates. - The remaining fields count the inner linear-algebra work and describe the preconditioner; their names and meanings may change between releases: `nsolves`, `lsqriters`, `lsqrtrace`, `cgiters`, `cholsolves`, `precond`, - `nrefactor`, `nforest`, `ndiagonal`, `fill_entries`, `factor_flops`, - `nzeroed`. + `nrefactor`, `nforest`, `fill_entries`, `factor_flops`, `nzeroed`. """ struct AugmentedLagrangianStats{T} <: AbstractCoverStats converged::Bool @@ -123,7 +122,6 @@ struct AugmentedLagrangianStats{T} <: AbstractCoverStats precond::Symbol nrefactor::Int nforest::Int - ndiagonal::Int fill_entries::Int factor_flops::Float64 nzeroed::Int @@ -183,7 +181,7 @@ function AugmentedLagrangianStats{T}(nt::NamedTuple) where {T} collect(T, nt.drops), polish, nt.linsolve, nt.nsolves, nt.lsqriters, collect(Int, nt.lsqrtrace), nt.cgiters, nt.cholsolves, nt.precond, nt.nrefactor, nt.nforest, - nt.ndiagonal, nt.fill_entries, nt.factor_flops, nt.nzeroed) + nt.fill_entries, nt.factor_flops, nt.nzeroed) end LeastSquaresStats(nt::NamedTuple) = diff --git a/src/tree_cholesky.jl b/src/tree_cholesky.jl new file mode 100644 index 0000000..234fbbf --- /dev/null +++ b/src/tree_cholesky.jl @@ -0,0 +1,125 @@ +# Cholesky factorization of a symmetric positive definite matrix whose +# off-diagonal pattern is a forest. Eliminating each tree from its leaves +# produces no fill, so the factor `L` of `M = L Lᵀ` has one diagonal entry per +# vertex and one off-diagonal entry per non-root vertex, and factorization, +# solves, and products are linear in the number of vertices. The arithmetic is +# that of the element type, so any element type serves. +struct TreeCholesky{T} + parent::Vector{Int} # parent vertex; 0 for a root + order::Vector{Int} # the vertices, each after all of its children + d::Vector{T} # `L[v, v]` + l::Vector{T} # `L[parent[v], v]`; 0 for a root +end + +# Factor the `N×N` matrix with diagonal `diag` and off-diagonal entries `V[k]` +# at `(I[k], J[k])` and `(J[k], I[k])`. The edges must form a forest, and the +# matrix must be positive definite; either failure throws. +function tree_cholesky(N::Int, I::AbstractVector{Int}, J::AbstractVector{Int}, + V::AbstractVector{T}, diag::AbstractVector{T}) where {T} + axes(I) == axes(J) == axes(V) || + throw(DimensionMismatch("edge endpoints and values must share axes, got $(axes(I)), $(axes(J)), $(axes(V))")) + length(diag) == N || throw(DimensionMismatch("the diagonal has $(length(diag)) entries; expected $N")) + ptr = zeros(Int, N + 1) + for k in eachindex(I) + ptr[I[k]+1] += 1 + ptr[J[k]+1] += 1 + end + ptr[1] = 1 + cumsum!(ptr, ptr) + adj = zeros(Int, ptr[end] - 1) + cursor = ptr[1:N] + for k in eachindex(I) + adj[cursor[I[k]]] = k + cursor[I[k]] += 1 + adj[cursor[J[k]]] = k + cursor[J[k]] += 1 + end + # Breadth-first search from each root lists parents before children. + parent = zeros(Int, N) + pedge = zeros(Int, N) + seen = falses(N) + bfs = Int[] + sizehint!(bfs, N) + for r in 1:N + seen[r] && continue + seen[r] = true + head = length(bfs) + 1 + push!(bfs, r) + while head <= length(bfs) + v = bfs[head] + head += 1 + for s in ptr[v]:ptr[v+1]-1 + k = adj[s] + k == pedge[v] && continue + w = I[k] == v ? J[k] : I[k] + seen[w] && throw(ArgumentError("the edges contain a cycle through vertex $w")) + seen[w] = true + parent[w] = v + pedge[w] = k + push!(bfs, w) + end + end + end + order = reverse!(bfs) + d = zeros(T, N) + l = zeros(T, N) + s = collect(T, diag) # each pivot, less the squares of its children's factor entries + for v in order + s[v] > 0 || throw(ArgumentError("the matrix is not positive definite: pivot $(s[v]) at vertex $v")) + d[v] = sqrt(s[v]) + p = parent[v] + if p != 0 + l[v] = V[pedge[v]] / d[v] + s[p] -= l[v]^2 + end + end + return TreeCholesky{T}(parent, order, d, l) +end + +# `x = L \ b`; `x` may be `b`. +function solve_ptl!(x::AbstractVector, F::TreeCholesky, b::AbstractVector) + x === b || copyto!(x, b) + for v in F.order + x[v] /= F.d[v] + p = F.parent[v] + p == 0 || (x[p] -= F.l[v] * x[v]) + end + return x +end + +# `x = Lᵀ \ b`; `x` may be `b`. +function solve_up!(x::AbstractVector, F::TreeCholesky, b::AbstractVector) + x === b || copyto!(x, b) + for v in Iterators.reverse(F.order) + p = F.parent[v] + p == 0 || (x[v] -= F.l[v] * x[p]) + x[v] /= F.d[v] + end + return x +end + +# `y = Lᵀ x`; `y` may be `x`. +function lt_mul!(y::AbstractVector, F::TreeCholesky, x::AbstractVector) + for v in F.order + p = F.parent[v] + y[v] = F.d[v] * x[v] + (p == 0 ? zero(eltype(y)) : F.l[v] * x[p]) + end + return y +end + +# Indices of the edges of a maximum-weight spanning forest of the graph on +# `1:N` with edges `(I[k], J[k])` and weights `W[k]`, by Kruskal's algorithm. +function max_weight_forest(N::Int, I::AbstractVector{Int}, J::AbstractVector{Int}, W::AbstractVector) + axes(I) == axes(J) == axes(W) || + throw(DimensionMismatch("edge endpoints and weights must share axes, got $(axes(I)), $(axes(J)), $(axes(W))")) + uf = collect(1:N) + find(x) = (while uf[x] != x; uf[x] = uf[uf[x]]; x = uf[x]; end; x) + keep = Int[] + for k in sortperm(W; rev=true) + rp, rq = find(I[k]), find(J[k]) + rp == rq && continue + uf[rp] = rq + push!(keep, k) + end + return keep +end diff --git a/test/element_types.jl b/test/element_types.jl index ee08b8b..6cad9aa 100644 --- a/test/element_types.jl +++ b/test/element_types.jl @@ -73,11 +73,20 @@ @test abig isa Vector{BigFloat} @test Float64.(abig) ≈ first(symcover_min(AbsLog{2}(), Float64.(Abig))) rtol=1e-6 @test abig != BigFloat.(Float32.(abig)) - # LSQR on such a type is preconditioned by the diagonal, since the - # factor preconditioners are CHOLMOD's. + # LSQR on such a type is preconditioned by the forest factor, since the + # full factor is CHOLMOD's. abig_lsqr, sbig = symcover_min(AbsLog{2}(), sparse(Abig); linsolve=:lsqr) - @test MatrixCovers.converged(sbig) && sbig.precond === :diagonal + @test MatrixCovers.converged(sbig) && sbig.precond === :forest @test Float64.(abig_lsqr) ≈ Float64.(abig) rtol=1e-6 + # Double64 goes through the same path at its own precision, and the + # forest keeps the iteration count near Float64's. + Adbl = Double64.(sparse(Asym)) + adbl, sdbl = symcover_min(AbsLog{2}(), Adbl; linsolve=:lsqr) + a64, s64 = symcover_min(AbsLog{2}(), Float64.(sparse(Asym)); linsolve=:lsqr, fillbudget=0) + @test adbl isa Vector{Double64} + @test MatrixCovers.converged(sdbl) && sdbl.precond === :forest + @test Float64.(adbl) ≈ a64 rtol=1e-8 + @test sdbl.lsqriters <= 4 * s64.lsqriters end @testset "BigFloat flows through the family" begin diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index ed0c442..596e7ce 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -423,8 +423,8 @@ end end # The LSQR preconditioner runs in two regimes: a Cholesky factor of the weighted -# normal matrix while the fill budget allows it, and that matrix's diagonal -# beyond it. Both must reach the same cover. +# normal matrix while the fill budget allows it, and a spanning forest of the +# active couplings beyond it. Both must reach the same cover. @testset "MMC :lsqr preconditioner regimes" begin rng = StableRNG(23) n = 200 @@ -437,34 +437,36 @@ end af, sf = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) ad, sd = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr, fillbudget=0) @test sf.precond === :factor - @test sd.precond === :diagonal - # The factored regime (re)factorizes the preconditioner at least once and - # reports its predicted fill; the diagonal regime never factorizes. - @test sf.nrefactor >= 1 + @test sd.precond === :forest + # Both regimes (re)factorize their preconditioner at least once; only the + # full factor has predicted fill, and the forest regime still analyzes the + # full pattern for the budget check. + @test sf.nrefactor >= 1 && sd.nrefactor >= 1 @test sf.fill_entries > 0 @test sf.factor_flops > 0 - @test sd.nrefactor == 0 - @test sd.factor_flops > 0 # analyzed for the fill check even though the diagonal path is taken + @test sd.factor_flops > 0 + @test sf.nforest == 0 && sd.nforest >= 1 @test af ≈ ad rtol=1e-6 @test iscover(af, af, A) && iscover(ad, ad, A) @test cover_objective(AbsLog{2}(), af, af, A) ≈ cover_objective(AbsLog{2}(), ad, ad, A) rtol=1e-8 - # The factored regime is the one that converges in a few iterations per solve. + # The full factor is the one that converges in a few iterations per solve. @test sf.lsqriters < sd.lsqriters gf, hf, tf = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr) gd, hd, td = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr, fillbudget=0) - @test (tf.precond, td.precond) === (:factor, :diagonal) + @test (tf.precond, td.precond) === (:factor, :forest) @test gf .* hf' ≈ gd .* hd' rtol=1e-6 @test iscover(gf, hf, G) && iscover(gd, hd, G) @test cover_objective(AbsLog{2}(), gf, hf, G) ≈ cover_objective(AbsLog{2}(), gd, hd, G) rtol=1e-8 # The budget is a keyword of the public solvers, the dense path never # preconditions, and an element type CHOLMOD does not factor uses the - # diagonal. + # forest throughout. @test first(symcover_min(AbsLog{2}(), A; linsolve=:lsqr, fillbudget=0)) ≈ ad rtol=1e-6 @test MatrixCovers._symcover_min_abslog2(A; linsolve=:dense)[2].precond === :none Abig = BigFloat.([4.0 1.0 0.5; 1.0 3.0 1.0; 0.5 1.0 2.5]) - @test MatrixCovers._symcover_min_abslog2(Abig; linsolve=:lsqr)[2].precond === :diagonal + sbig = MatrixCovers._symcover_min_abslog2(Abig; linsolve=:lsqr)[2] + @test sbig.precond === :forest && sbig.fill_entries == 0 end """ @@ -502,7 +504,7 @@ end # refactorization cost, not just storage: a banded matrix has cheap fill-in and # stays on the factor path even at large `n`, while a matrix with the same # density but random sparsity fills in enough to cost more to refactorize than -# the heavy-forest preconditioner it falls back to. +# the forest preconditioner it falls back to. @testset "MMC :lsqr preconditioner cost rule" begin rngb = StableRNG(31) Aband = banded_sparse_sym(rngb, 2000, 3, 1.0) @@ -522,19 +524,18 @@ end # The in-solver ratio of predicted flops to stored support entries matches # `factor_flops / nnz(A)` of the input for a symmetric matrix with both # triangles stored: a flop budget just above that ratio selects the factor - # path, and just below it selects the heavy-forest path. + # path, and just below it selects the forest path. ratio = sf.factor_flops / nnz(Arand) @test MatrixCovers._symcover_min_abslog2(Arand; linsolve=:lsqr, flopbudget=1.01 * ratio)[2].precond === :factor @test MatrixCovers._symcover_min_abslog2(Arand; linsolve=:lsqr, flopbudget=0.99 * ratio)[2].precond === :forest - @test MatrixCovers._symcover_min_abslog2(Arand; linsolve=:lsqr, flopbudget=0)[2].precond === :diagonal + @test MatrixCovers._symcover_min_abslog2(Arand; linsolve=:lsqr, flopbudget=0)[2].precond === :forest end -# Random sparse supports have expensive full factors. The heavy-forest -# preconditioner (violated entries with their couplings, the other entries on -# the diagonal) must reach the same cover as the full factor, and fall back to -# the diagonal for patterns whose forest factor is over budget. -@testset "MMC :lsqr heavy-forest preconditioner" begin +# Random sparse supports have expensive full factors. The forest +# preconditioner (a spanning forest of the couplings among violated entries, +# every entry on the diagonal) must reach the same cover as the full factor. +@testset "MMC :lsqr forest preconditioner" begin rng = MersenneTwister(5) n = 400 S = sprand(rng, n, n, 5 / n) + sparse(randperm(rng, n), 1:n, ones(n), n, n) @@ -549,7 +550,7 @@ end @test tf.precond === :factor @test ts.precond === :forest @test ts.nforest > 0 - @test ts.nforest + ts.ndiagonal <= ts.nsolves + @test ts.nforest <= ts.nsolves @test ts.nrefactor >= 1 @test gs[Gi] .* hs[Gj] ≈ gf[Gi] .* hf[Gj] rtol=1e-8 @test iscover(gs, hs, G) @@ -567,8 +568,8 @@ end @test sf.precond === :factor @test ss.precond === :forest @test ss.nforest > 0 - @test ss.nforest + ss.ndiagonal <= ss.nsolves - @test as[Ai] .* as[Aj] ≈ af[Ai] .* af[Aj] rtol=1e-12 + @test ss.nforest <= ss.nsolves + @test as[Ai] .* as[Aj] ≈ af[Ai] .* af[Aj] rtol=1e-10 @test iscover(as, as, A) @test cover_objective(AbsLog{2}(), as, as, A) ≈ cover_objective(AbsLog{2}(), af, af, A) rtol=1e-12 cfn = M -> MatrixCovers._symcover_min_abslog2(M; flopbudget=fb)[1] @@ -596,14 +597,15 @@ end end @test af[Ai] .* af[Aj] ≈ ga[Ai] .* gb[Aj] rtol=1e-12 - # A flop budget just at what the forest factors need sends some patterns - # to the diagonal preconditioner, and the cover is unchanged. + # The forest factor needs no budget: any budget below the full factor's + # gives the forest regime for every solve, and the cover is unchanged. gf, hf, _ = MatrixCovers._cover_min_abslog2(G; flopbudget=Inf) - gm, hm, tm = MatrixCovers._cover_min_abslog2(G; flopbudget=1.0) - @test tm.nforest > 0 && tm.ndiagonal > 0 - @test tm.nforest + tm.ndiagonal <= tm.nsolves - @test gm[Gi] .* hm[Gj] ≈ gf[Gi] .* hf[Gj] rtol=1e-8 - @test MatrixCovers._cover_min_abslog2(G; flopbudget=0)[3].nforest == 0 + for budget in (1.0, 0) + gm, hm, tm = MatrixCovers._cover_min_abslog2(G; flopbudget=budget) + @test tm.precond === :forest + @test 0 < tm.nforest <= tm.nsolves + @test gm[Gi] .* hm[Gj] ≈ gf[Gi] .* hf[Gj] rtol=1e-8 + end end # LSQR continuation starts from the heuristic cover. @@ -933,7 +935,7 @@ end # `2(κ-1)z` per update they hold the KKT residual at 2.6e-6. # Without the active-set finish (`polish=false`), slack multipliers are zeroed. # Whether the residual then reaches the tolerance depends on rounding (with the - # diagonal preconditioner it can stall just above it); the solve warns exactly + # forest preconditioner it can stall just above it); the solve warns exactly # when it does not. The default solve's finish certifies the result. A = banded_sparse_sym(StableRNG(5), 10_000, 3, 1.0) for fillbudget in (MatrixCovers.LSQR_FILL_BUDGET, 0) @@ -942,7 +944,7 @@ end end @test isempty(logs) == s.converged @test all(l -> l.level == Base.CoreLogging.Warn && occursin("try a larger `κ`", l.message), logs) - @test s.precond === (fillbudget == 0 ? :diagonal : :factor) + @test s.precond === (fillbudget == 0 ? :forest : :factor) @test s.nzeroed > 0 @test maximum(s.κs) == 1e5 @test s.kkt[end] < 1e-9 diff --git a/test/runtests.jl b/test/runtests.jl index 68e6a87..96e6e7e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -27,6 +27,7 @@ include("helpers.jl") # isbalanced, covaries, PENALTIES include("powermean.jl") # the PowerMean soft covers (the soft default) include("initializers.jl") # initialize_symcover/initialize_cover strategies include("sparse_cholesky.jl") # the CHOLMOD driver and the multi-array eachindex helper + include("tree_cholesky.jl") # the fill-free forest factorization behind the LSQR preconditioner include("minimal_covers.jl") # the *_min family (native solvers) include("transversal.jl") # cover_transversal and the maximum-product transversal include("active_set.jl") # exact solution of stalled AbsLog{2} minimal covers diff --git a/test/tree_cholesky.jl b/test/tree_cholesky.jl new file mode 100644 index 0000000..f6bbef6 --- /dev/null +++ b/test/tree_cholesky.jl @@ -0,0 +1,43 @@ +using MatrixCovers: TreeCholesky, tree_cholesky, max_weight_forest, solve_ptl!, solve_up!, lt_mul! + +@testset "tree Cholesky" begin + rng = StableRNG(31) + for T in (Float64, BigFloat), trial in 1:5 + N = 40 + # A random forest: most vertices attach to an earlier one. + I = Int[] + J = Int[] + for v in 2:N + rand(rng) < 0.8 && (push!(I, v); push!(J, rand(rng, 1:v-1))) + end + V = T.(randn(rng, length(I))) + # Diagonal dominance makes the matrix positive definite. + diag = T.(abs.(randn(rng, N)) .+ 1) + for k in eachindex(I) + diag[I[k]] += abs(V[k]) + diag[J[k]] += abs(V[k]) + end + M = Matrix(sparse([I; J; 1:N], [J; I; 1:N], [V; V; diag], N, N)) + F = tree_cholesky(N, I, J, V, diag) + @test F isa TreeCholesky{T} + b = T.(randn(rng, N)) + x = solve_up!(zeros(T, N), F, solve_ptl!(zeros(T, N), F, b)) + @test M * x ≈ b + # `L \ (M x) == Lᵀ x` shows `M == L Lᵀ`. + @test solve_ptl!(zeros(T, N), F, M * x) ≈ lt_mul!(zeros(T, N), F, x) + # Every operation works in place. + xa = copy(b) + solve_ptl!(xa, F, xa) + solve_up!(xa, F, xa) + @test xa ≈ x + ya = copy(x) + lt_mul!(ya, F, ya) + @test ya ≈ lt_mul!(zeros(T, N), F, x) + end + @test_throws "contain a cycle" tree_cholesky(3, [1, 2, 3], [2, 3, 1], ones(3), fill(10.0, 3)) + @test_throws "not positive definite" tree_cholesky(2, [1], [2], [5.0], [1.0, 1.0]) + @test_throws DimensionMismatch tree_cholesky(2, [1], [2], [0.5], [1.0, 1.0, 1.0]) + # Kruskal keeps the two heaviest edges of the triangle and the pendant edge. + @test sort(max_weight_forest(4, [1, 2, 3, 1], [2, 3, 1, 4], [1.0, 3.0, 2.0, 0.5])) == [2, 3, 4] + @test max_weight_forest(3, Int[], Int[], Float64[]) == Int[] +end