Repository navigation
Fix active-set finish stall; precondition in any element type - #85
Merged
Merged
Conversation
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 <noreply@anthropic.com>
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 <noreply@anthropic.com>
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 <noreply@anthropic.com>
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #85 +/- ##
==========================================
- Coverage 99.07% 98.79% -0.28%
==========================================
Files 22 23 +1
Lines 6132 6286 +154
==========================================
+ Hits 6075 6210 +135
- Misses 57 76 +19 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Since #83 the finish solves over-budget tree Laplacians by conjugate gradients, whose forward error (1.5e-11 on a sparse nonsymmetric lognormal matrix with n = 31623) can exceed the certificate's
ptolof 1e-12. Constraints on no cycle of the support are tight with zero multiplier at the minimizer, so that noise made a few of them look violated, and the forest shuffled among them until the step limit of 4n, about 5 h per solve.Finish (
active_set.jl)_acyclic_constraintspeels degree-one vertices of the support, and the seed activates those constraints first. Their residuals are exactly zero by construction, so they are never targeted, dropped, or blocking._laplacian_solve!on_LaplacianCGrefines the CG solution with a double-double residual until the correction is below 8 eps ‖x‖ and records the last correction as an error estimate._dual_feasible!raises the certificate's tolerances to twice that estimate, socertifiednever claims more accuracy than the solve delivers._LaplacianCGis parametric in the element type. Sparse tree Laplacians in types CHOLMOD does not factor take refined CG instead of an O(n^3) dense factorization, which remains only as the fallback when refinement stops contracting.Multiplier iteration (
minimal_covers.jl, newtree_cholesky.jl)ndiagonalstatistic.precondis now:factor,:forest, or:none.Results