SciML · ChrisRackauckas-Claude · May 30, 2026
diff --git a/README.md b/README.md
@@ -106,6 +106,62 @@ The `refactor!` + `\` loop is allocation-free apart from the tiny `O(r)` dense p
 * **`:auto`** (default) — Woodbury, automatically falling back to the augmented system if
   `S` is singular or `C` is ill-conditioned.
 
+### QR (the numerically stable path)
+
+`qr(A; strategy=:augmented)` factors the bordered system `[S U; V -I]` with **SuiteSparseQR**
+(the stdlib sparse `qr`), returning a `SparseWithDenseRowColQRAugmented`. Like the augmented
+LU path it never forms `S⁻¹`, so it stays accurate when the sparse part `S` is ill-conditioned
+or nearly singular but the dense rows regularize `A` — the
+[FastAlmostBandedMatrices](https://github.com/SciML/FastAlmostBandedMatrices.jl) regime, where
+QR is the conventional choice. Measured forward error on `A = S + U V` with `κ(A) ≈ 6` as
+`κ(S)` sweeps from `3.6e6` to `3.6e14` (reference: BigFloat dense solve):
+
+```
+  κ(S)          qr (augmented)   factorize :woodbury, refine=2   factorize :augmented (LU)
+  3.6e6   …   3.1e-16 (flat)        8.7e-17  →  4.5e-7              1.2e-16 (flat)
+  3.6e14
+```
+
+The Woodbury path's error grows with `κ(S)·κ(C)`; here (where the *assembled* `A` stays
+well-conditioned, `κ(A) ≈ 6`) both augmented paths hold at the noise floor.
+
+**Be clear about what QR does and does not buy you.** Augmented QR is *not* more accurate than
+the augmented **LU** path. It is backward-stable, so its forward error tracks `κ(A)·eps`; when
+the *assembled* `A` is itself ill-conditioned the augmented-LU path's partial pivoting can be
+**more** accurate on these structured matrices (measured: LU near the noise floor where QR sits
+at `κ(A)·eps`, orders apart for `κ(A) ≳ 1e6`). QR's genuine, distinct properties versus
+augmented-LU are that it is **rank-revealing** (a clean numerical-rank/singularity signal) and
+that it is the QR-based interface users of almost-banded solvers expect — not higher accuracy,
+which the LU path already delivers, and not throughput (it is slower; see the trade-offs
+below). Reach for `qr` when you specifically want a QR factorization or its rank diagnostics;
+otherwise `factorize`/`lu` is the better default.
+
+There is **no** `strategy=:woodbury` for `qr`: a Woodbury-over-`qr(S)` solve shares the same
+`κ(S)·κ(C)` cancellation and is catastrophically inaccurate on ill-conditioned `S` (~`1e-1`),
+so it is deliberately not offered. `qr` throws `SingularException` when the bordered system is
+numerically rank-deficient (`κ ≳ 1/(m·eps)`); a merely ill-conditioned, accurately-solvable `A`
+is **not** rejected (SuiteSparseQR's aggressive default column deflation is disabled).
+
+Trade-offs relative to `factorize`/`lu`: QR supports **only** `Float64`/`ComplexF64`
+(SuiteSparseQR's single-precision sparse QR is version-dependent — it throws on Julia 1.11 — so
+`Float32`/`ComplexF32` route to the LU path), has **no** symbolic-reuse `refactor!`
+(SuiteSparseQR exposes no `qr!`, so it re-factors fully), and **allocates** per solve (no
+in-place `ldiv!`). It is the stability path, not the throughput path — for `Float32`/generic
+eltypes or hot Newton refactor loops use `factorize`/`lu`. The cost is real and consistent
+(banded interior + `r` dense rows, factor / solve / refactor, ms):
+
+```
+  n=2000 r=4                          n=5000 r=8
+  qr  (augmented QR) : 67 / 7.5 / 90    813 / 156 / 776     ← stability path, slowest
+  lu  (augmented LU) : 47 / 1.7 / 21    633 /  11 / 188     ← klu! symbolic reuse on refactor
+  lu  (Woodbury)     : 1.3/ 0.1 / n/a   3.6 / 0.3 / n/a     ← fast path (forms S⁻¹), ~50–60×
+```
+
+QR is slower than the augmented-LU path on every axis (≈1.4× factor, ≈4–14× solve, ≈4×
+refactor since it has no symbolic reuse) and far slower than the Woodbury fast path; forward
+error is comparable (~`1e-13`) for all three on a well-conditioned `A`. Use `qr` for the
+factorization/rank-revealing semantics, not for speed.
+
 ## Benchmarks
 
 A banded sparse interior plus `r = 8` dense rows (`n = 5000`). Factoring only `S` and
@@ -139,10 +195,12 @@ globally-scattered `S` degrades every sparse method.)
 * [x] Factorization and linear solve (`factorize`/`lu`, `\`, `ldiv!`)
 * [x] In-place refactorization with a fixed sparsity pattern (`refactor!`)
 * [x] Adjoint / transpose solves
-* [x] Generic element types (`BigFloat`, `ForwardDiff.Dual`, complex)
+* [x] Generic element types (`BigFloat`, `ForwardDiff.Dual`, complex) — LU/Woodbury path
+* [x] Rank-revealing QR (`qr`) for ill-conditioned `S` (`Float64`/`ComplexF64` only)
 
-`qr` is intentionally not provided — a general sparse `QR` causes far more fill than the
-LU-based KLU path; use `factorize`/`lu`.
+`qr(A)` builds the augmented sparse-QR factorization (above). It is the stability path, not
+the throughput path: no symbolic-reuse refactor, no zero-allocation solve, and `Float64`
+/`ComplexF64` only — use `factorize`/`lu` for `Float32`/generic eltypes or hot refactor loops.
 
 ## When is this the right tool?
 
@@ -248,15 +306,21 @@ not throw), again matching KLU. One caveat: once a singular `S` forces the augme
 the cache keeps reusing it even if later `S` values are nonsingular — `init` a fresh cache to
 get back the faster Woodbury path.
 
+For the stable QR path, pass `SparseWithDenseRowColQRFactorization(; check_pattern)` instead.
+It has no `reuse_symbolic`/`strategy`/`refine` (SuiteSparseQR has no symbolic reuse, only the
+augmented mode exists, and augmented QR needs no refinement); singular `A` likewise returns
+`ReturnCode.Infeasible`. It is opt-in — the default algorithm stays the LU one.
+
 ## Public API
 
 ```
 SparseWithDenseRowColMatrix    SelectorMatrix
 sparsepart  fillpart  lowrankfactors  exclusive_sparsepart  denserank
-factorize   lu   \\   ldiv!   refactor!   lu!   update_lowrank!
-SparseWithDenseRowColWoodbury  SparseWithDenseRowColAugmented
+factorize   lu   qr   \\   ldiv!   refactor!   lu!   qr!   update_lowrank!
+SparseWithDenseRowColWoodbury  SparseWithDenseRowColAugmented  SparseWithDenseRowColQRAugmented
 recommend_lowrank_peel   PeelRecommendation
-SparseWithDenseRowColFactorization   # LinearSolve.jl algorithm (extension)
+SparseWithDenseRowColFactorization     # LinearSolve.jl algorithm (LU/Woodbury, default)
+SparseWithDenseRowColQRFactorization   # LinearSolve.jl algorithm (stable QR, opt-in)
 ```
 
 ## Installation

diff --git a/ext/SparseWithDenseRowColMatricesLinearSolveExt.jl b/ext/SparseWithDenseRowColMatricesLinearSolveExt.jl
@@ -2,10 +2,11 @@ module SparseWithDenseRowColMatricesLinearSolveExt
 
 using SparseWithDenseRowColMatrices
 using SparseWithDenseRowColMatrices: SparseWithDenseRowColMatrix, SparseWithDenseRowColWoodbury,
-    SparseWithDenseRowColAugmented, denserank, refactor!, SparseWithDenseRowColFactorization
+    SparseWithDenseRowColAugmented, SparseWithDenseRowColQRAugmented, denserank, refactor!,
+    SparseWithDenseRowColFactorization, SparseWithDenseRowColQRFactorization
 using LinearSolve
 using LinearSolve: LinearCache, OperatorAssumptions, LinearVerbosity, AbstractSparseFactorization
-using LinearAlgebra: LinearAlgebra, ldiv!, issuccess, Factorization, SingularException, factorize
+using LinearAlgebra: LinearAlgebra, ldiv!, issuccess, Factorization, SingularException, factorize, qr
 
 # LinearSolve `@reexport using SciMLBase`, so SciMLBase is reachable through it.
 const SciMLBase = LinearSolve.SciMLBase
@@ -105,4 +106,67 @@ function SciMLBase.solve!(cache::LinearCache, alg::SWDRCFactorizationAlg; kwargs
     end
 end
 
+# ------------------
+# QR algorithm (numerically stable, opt-in)
+# ------------------
+# Mirrors SWDRCFactorizationAlg but factors via the augmented sparse QR. There is no
+# `reuse_symbolic`/`strategy`/`refine`: SuiteSparseQR has no symbolic-reuse refactor (a value
+# update always re-`qr`s), only the augmented mode exists, and augmented QR needs no refinement.
+struct SWDRCQRFactorizationAlg <: AbstractSparseFactorization
+    check_pattern::Bool
+end
+
+function SparseWithDenseRowColMatrices.SparseWithDenseRowColQRFactorization(; check_pattern::Bool = true)
+    return SWDRCQRFactorizationAlg(check_pattern)
+end
+
+function LinearSolve.init_cacheval(
+        ::SWDRCQRFactorizationAlg, A::SparseWithDenseRowColMatrix{T}, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    ) where {T}
+    return SWDRCCacheval{T}(nothing)
+end
+
+# Refresh the cached QR factorization. `refactor!` re-`qr`s in place (reusing the buffers and
+# the bordered-matrix shape check); it throws on a singular `A` or a changed pattern, in which
+# case we rebuild with `qr`. A genuinely singular `A` makes `qr` throw `SingularException`,
+# which we map to a `nothing` cache → `Infeasible` (matching the LU algorithm's contract).
+function _refresh_qr!(wrap::SWDRCCacheval, A::SparseWithDenseRowColMatrix, alg::SWDRCQRFactorizationAlg)
+    f = wrap.fact
+    if f isa SparseWithDenseRowColQRAugmented &&
+            size(f) == size(A) && denserank(f) == denserank(A)
+        try
+            refactor!(f, A; check = alg.check_pattern)
+            return wrap
+        catch e
+            (e isa SingularException || e isa ArgumentError || e isa DimensionMismatch) || rethrow(e)
+        end
+    end
+    try
+        wrap.fact = qr(A)
+    catch e
+        e isa SingularException || rethrow(e)
+        wrap.fact = nothing
+    end
+    return wrap
+end
+
+function SciMLBase.solve!(cache::LinearCache, alg::SWDRCQRFactorizationAlg; kwargs...)
+    if cache.isfresh
+        wrap = LinearSolve.@get_cacheval(cache, :SWDRCQRFactorizationAlg)
+        _refresh_qr!(wrap, cache.A, alg)
+        cache.isfresh = false
+    end
+    F = LinearSolve.@get_cacheval(cache, :SWDRCQRFactorizationAlg).fact
+    return if F !== nothing && issuccess(F)
+        y = ldiv!(cache.u, F, cache.b)
+        SciMLBase.build_linear_solution(alg, y, nothing, cache; retcode = SciMLBase.ReturnCode.Success)
+    else
+        SciMLBase.build_linear_solution(
+            alg, cache.u, nothing, cache; retcode = SciMLBase.ReturnCode.Infeasible
+        )
+    end
+end
+
 end # module
diff --git a/src/SparseWithDenseRowColMatrices.jl b/src/SparseWithDenseRowColMatrices.jl
@@ -1,8 +1,8 @@
 module SparseWithDenseRowColMatrices
 
 using LinearAlgebra
-using LinearAlgebra: Factorization, SingularException, ldiv!, lu, lu!, mul!, opnorm,
-    issuccess, factorize
+using LinearAlgebra: Factorization, SingularException, ldiv!, lu, lu!, qr, qr!, mul!, opnorm,
+    issuccess, factorize, diag
 using SparseArrays
 using SparseArrays: AbstractSparseMatrixCSC, SparseMatrixCSC, nonzeros, rowvals, getcolptr,
     nnz, sparse
@@ -15,13 +15,16 @@ include("matrix.jl")
 include("matvec.jl")
 include("woodbury.jl")
 include("augmented.jl")
+include("qr.jl")
 include("factorize.jl")
 include("detect.jl")
 
 export SparseWithDenseRowColMatrix, SelectorMatrix,
     sparsepart, fillpart, lowrankfactors, exclusive_sparsepart, denserank,
-    SparseWithDenseRowColWoodbury, SparseWithDenseRowColAugmented, refactor!, update_lowrank!,
-    recommend_lowrank_peel, PeelRecommendation, SparseWithDenseRowColFactorization
+    SparseWithDenseRowColWoodbury, SparseWithDenseRowColAugmented,
+    SparseWithDenseRowColQRAugmented, refactor!, update_lowrank!,
+    recommend_lowrank_peel, PeelRecommendation,
+    SparseWithDenseRowColFactorization, SparseWithDenseRowColQRFactorization
 
 # ---------------------
 # Precompilation

diff --git a/src/factorize.jl b/src/factorize.jl
@@ -70,10 +70,50 @@ Alias for [`refactor!`](@ref) — refactor `F` in place with the new values of `
 LinearAlgebra.lu!(F::SparseWithDenseRowColFactorization, A::SparseWithDenseRowColMatrix; kwargs...) =
     refactor!(F, A; kwargs...)
 
-function LinearAlgebra.qr(::SparseWithDenseRowColMatrix)
-    throw(ArgumentError("SparseWithDenseRowColMatrix is solved by an LU-based factorization; use `factorize`/`lu`, not `qr`."))
+"""
+    qr(A::SparseWithDenseRowColMatrix; strategy=:auto) -> SparseWithDenseRowColQRAugmented
+
+QR-factorize `A = S + U V` for numerically stable repeated solves. Builds a single sparse QR
+(SuiteSparseQR) of the bordered system `[S U; V -I]`, returning a
+[`SparseWithDenseRowColQRAugmented`](@ref). Unlike [`factorize`](@ref)/[`lu`](@ref) (the
+LU/Woodbury throughput path) this never forms `S⁻¹`, so it stays accurate even when the sparse
+part `S` is ill-conditioned or nearly singular, as long as `A` itself is nonsingular — the
+FastAlmostBandedMatrices regime.
+
+* `strategy = :auto` (default) or `:augmented`: the augmented bordered-system QR (the only mode).
+* `strategy = :woodbury`: **not supported** — a Woodbury-over-qr(S) approach shares the
+  κ(S)·κ(C) cancellation of the LU Woodbury path and is catastrophically inaccurate on
+  ill-conditioned `S`, defeating the purpose of using QR. Throws.
+
+Only `Float64`/`ComplexF64` are supported (SuiteSparseQR's single-precision sparse QR is
+unavailable on some Julia versions); use `factorize`/`lu` for `Float32`/`ComplexF32` and
+generic eltypes. The factorization has no symbolic-reuse refactor and no
+zero-allocation solve (SuiteSparseQR limitations); it is the stability path, not the
+throughput path. Solve with `F \\ b` / `ldiv!(F, b)`; update values with [`refactor!`](@ref) /
+[`qr!`](@ref) (a full re-`qr`).
+"""
+function LinearAlgebra.qr(A::SparseWithDenseRowColMatrix; strategy::Symbol = :auto)
+    strategy in (:auto, :augmented) && return _augmented_qr(A)
+    strategy === :woodbury && throw(
+        ArgumentError(
+            "qr does not support strategy = :woodbury: a Woodbury-over-qr(S) solve shares the \
+         κ(S)·κ(C) cancellation of the LU Woodbury path and is inaccurate on ill-conditioned S, \
+         defeating the purpose of QR. Use `qr(A; strategy = :augmented)` (the default), or \
+         `factorize(A; strategy = :woodbury)` for the LU Woodbury path."
+        )
+    )
+    throw(ArgumentError("strategy must be :auto or :augmented for qr; got :$strategy"))
 end
 
+"""
+    qr!(F::SparseWithDenseRowColQRAugmented, A::SparseWithDenseRowColMatrix; kwargs...) -> F
+
+Alias for [`refactor!`](@ref) — re-`qr` `F` in place with the new values of `A`. Note this is a
+full re-factorization (SuiteSparseQR has no symbolic reuse), unlike [`lu!`](@ref).
+"""
+LinearAlgebra.qr!(F::SparseWithDenseRowColQRAugmented, A::SparseWithDenseRowColMatrix; kwargs...) =
+    refactor!(F, A; kwargs...)
+
 Base.:\(A::SparseWithDenseRowColMatrix, b::AbstractVecOrMat) = factorize(A) \ b
 
 """
@@ -97,3 +137,20 @@ It is also the default algorithm `LinearSolve` picks for a `SparseWithDenseRowCo
 Only available when `LinearSolve` is loaded (provided by a package extension).
 """
 function SparseWithDenseRowColFactorization end
+
+"""
+    SparseWithDenseRowColQRFactorization(; check_pattern=true)
+
+A [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) algorithm that solves
+`SparseWithDenseRowColMatrix` systems with the numerically stable augmented sparse QR
+([`SparseWithDenseRowColQRAugmented`](@ref)) through LinearSolve's caching interface. Opt-in
+(the default LinearSolve algorithm for the matrix type remains the LU
+[`SparseWithDenseRowColFactorization`](@ref)); choose this when `S` is ill-conditioned.
+
+* `check_pattern` — validate the pattern on each refactor (`false` skips the check).
+
+There is no `reuse_symbolic` option (SuiteSparseQR has no symbolic-reuse refactor, so a value
+update always re-`qr`s), no `strategy` (only the augmented mode exists), and no `refine` (the
+augmented QR needs no Woodbury refinement). Only available when `LinearSolve` is loaded.
+"""
+function SparseWithDenseRowColQRFactorization end