QR via SparseColumnPivotedQR + minimum-norm least squares (lstsq)#6
Merged
Conversation
Two new capabilities for SparseWithDenseRowColMatrix (A = S + U V), plus the
structured matvec they need.
QR — `qr(A)` factors the bordered system [S U; V -I] with SparseColumnPivotedQR
(pure-Julia, rank-revealing, column-pivoted; the QR analogue of PureKLU). It is
the numerically stable analogue of the augmented-LU path (never forms S⁻¹), and
the backend gives it an allocation-free solve, a symbolic-reuse `refactor!`/`qr!`
(usable in Newton/time-stepping loops), genuine rank-revealing singularity
detection, and generic element types (BLAS floats plus BigFloat/Dual, on every
Julia version — no SuiteSparseQR single-precision/1.11 issues). No
strategy=:woodbury (Woodbury-over-qr(S) shares the κ(S)·κ(C) cancellation and is
catastrophically inaccurate on ill-conditioned S). Verified vs BigFloat/pinv.
lstsq — `lstsq(A, b)` returns the minimum-norm least-squares solution A⁺b for
rank-deficient/inconsistent systems where \\/factorize/qr throw. Default
alg=:auto picks a structure-exploiting DIRECT solve when applicable, else the
exact dense COD:
:structured — S nonsingular: A = S(I+ZV), rank deficiency collapses into the
r×r capacitance C = I + V S⁻¹U; A⁺b from one PureKLU factorization of S plus
small dense SVD/QR, never densifying A (O(factor(S)+n·r²), >100× faster than
dense; handles consistent and inconsistent b). A selector-singular S (BVP
boundary case) is peeled to a sparse nonsingular S̃ = S + U Uᴴ.
:dense — LAPACK gelsy COD, exact; the fallback for a general singular S.
:iterative — LSQR/LSMR via the structured matvec/adjoint (IterativeSolvers
extension), for large n.
Every formula was verified against pinv(Matrix(A))*b before implementation.
Also adds a structured adjoint/transpose matvec (Aᴴu = Sᴴu + Vᴴ(Uᴴu)) to the
core, fixing the ~1000× slower generic getindex fallback for A'*u.
LinearSolve integration: SparseWithDenseRowColFactorization (default) and an
opt-in SparseWithDenseRowColQRFactorization. New deps: SparseColumnPivotedQR +
SparseMatricesCSR (QR backend), IterativeSolvers (weakdep, lstsq iterative).
620/620 tests pass on Julia 1.10, 1.11, and 1.12. Runic-formatted, typos clean.
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
This was referenced May 31, 2026
…eated solves Adds a reusable structured least-squares factorization so a Newton / time-stepping loop (many solves, same A) does the expensive work once. `SparseWithDenseRowColLeastSquares(A)` caches S's PureKLU factorization and the small dense SVD/QR decompositions; each `F \ b` / `ldiv!(F, b)` is then a cheap back-solve that reuses them. Measured at n=2000, r=6: a reused solve is 0.05 ms / 16 KB vs the one-shot lstsq(A, b) at 2.1 ms / 2.6 MB — ~44× faster and ~160× less allocation, since it skips re-factoring S. The one-shot `lstsq(A, b)` is refactored to share the same setup/apply path (its post-hoc LS-optimality guard is preserved), so its behavior is unchanged. The constructor errors when the structured method does not apply (general singular S) — use lstsq(A, b; alg=:dense). 639/639 tests pass on Julia 1.10 and 1.11 (1.12 in progress; CI covers all three). Runic, typos clean. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
This was referenced May 31, 2026
Matrix-RHS adjoint/transpose matvec falls through to the generic dense fallback — ~3000x slowdown
#8
Open
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
2 participants
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.
Stable QR (SparseColumnPivotedQR backend) + minimum-norm least squares (
lstsq)Two new capabilities for
SparseWithDenseRowColMatrix(A = S + U V), plus the structured matvec they rely on. Everything was verified against apinv/BigFloat oracle before implementation; 620/620 tests pass on Julia 1.10, 1.11, and 1.12.qr(A)— stable, rank-revealing QR (now on SparseColumnPivotedQR)Factors the bordered system
[S U; V -I]with SparseColumnPivotedQR — pure-Julia, rank-revealing, column-pivoted (the QR analogue of PureKLU). The numerically stable analogue of the augmented-LU path (never formsS⁻¹); onA = S+UVwithκ(A)≈6andκ(S)swept to3.6e14it holds at ~3e-16while the Woodbury path degrades to~4.5e-7.The backend gives the QR path the engineering niceties it lacked under SuiteSparseQR:
ldiv!refactor!tol=0workaround)rankNo
strategy=:woodbury(Woodbury-over-qr(S)shares theκ(S)·κ(C)cancellation; ~1e-1error).lstsq(A, b)— minimum-norm least squaresA⁺bFor rank-deficient / inconsistent systems where
\/factorize/qrthrow.alg=:auto(default) picks a structure-exploiting direct solve when applicable, else the exact dense COD::structured— whenSis nonsingular,A = S(I+ZV)and the entire rank deficiency collapses into ther×rcapacitanceC = I + V S⁻¹U(nullity(A)=nullity(C)).A⁺bfrom one PureKLU factorization ofS+ small dense SVD/QR —Ais never densified (O(factor(S)+n·r²), >100× faster than dense at n=2000; handles consistent and inconsistentb). A selector-singularS(BVP boundary case) is peeled to a sparse nonsingularS̃ = S + UUᴴ.:dense— LAPACKgelsyCOD, exact; the fallback for a general singularS(whose null-space correction is dense — a genuine limit, not a gap).:iterative— LSQR/LSMR via the structured matvec/adjoint (IterativeSolvers extension), for largen.Also adds a structured adjoint/transpose matvec (
Aᴴu = Sᴴu + Vᴴ(Uᴴu)) to the core, fixing a ~1000× slower genericgetindexfallback forA'*u.Notes for review
Pkg.developstep can drop).refactor!allocates nothing (the kernel does; my wrapper rebuilds the bordered matrix — fixed the docstring); and the singular-SBVP case does have a structure-exploiting direct method (the peel), contrary to an earlier "falls back to dense" statement.Runic-formatted, typos clean.
🤖 Generated with Claude Code
Co-Authored-By: Chris Rackauckas accounts@chrisrackauckas.com