Skip to content

Solve with the adjoint of a sparse QR factorization - #804

Merged
ViralBShah merged 8 commits into
mainfrom
vs/spqr-adjoint-ldiv
Sep 16, 2026
Merged

ViralBShah merged 8 commits into
mainfrom
vs/spqr-adjoint-ldiv

Conversation

@ViralBShah

@ViralBShah ViralBShah commented Sep 9, 2026

Copy link
Copy Markdown
Member

Fixes #656, and completes the remaining half of #115 (qr(A') already works).

qr(A)' \ b previously threw

ERROR: MethodError: no method matching ldiv!(::AdjointFactorization{Float64, QRSparse{Float64, Int64}}, ::Vector{Float64})

This adds ldiv!(X, ::AdjointFactorization{<:Any,<:QRSparse}, B) plus the corresponding \ methods, so a single factorization of A can be used to solve with both A and A' — the use case from #115.

What it computes

With A[prow, pcol] == Q*R we have A' == Pcol*R'*Q'*Prow, so solving A'x = b is a forward substitution with R' followed by a multiply by Q and the inverse row permutation. It reuses the permutations, workspace and lock of the existing ldiv!, so a solve allocates only the result.

For a tall A the system A'x = b is underdetermined and the minimum-norm solution is returned (the trailing components of Q'x are zeroed), matching qr(A)' \ b for dense A and pinv(A') * b:

julia> A = sprand(3, 100, 0.1);

julia> b = randn(3);

julia> qr(sparse(A'))' \ b  Matrix(A) \ b
true

When A is rank deficient, zeroing the free variables also drops the equations that the leading rank block of R cannot represent — the counterpart of the basic solution that qr(A) \ b already returns.

Wide A is rejected with the same DimensionMismatch("overdetermined systems are not supported") that dense QR throws, since that solve needs a least-squares factorization of A' rather than of A.

Also fixed

\(::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, B) and the Transpose variant already routed non-square matrices through adjoint(qr(A)) \ B, so A' \ b and transpose(A) \ b threw a MethodError for any non-square sparse A. They now work for every shape: a tall A uses the new method above and returns the minimum-norm solution, while a wide A factorizes the sparse transpose and returns the least squares solution, like the dense path.

Relation to #301

#301 wants A \ b to return the minimum-norm solution for wide A, which needs exactly this method (qr(sparse(A'))' \ b). That is a behavior change to \ and is left for a separate PR; this one supplies the missing piece.

Tests

Added to the existing spqr.jl testsets: minimum-norm agreement with the dense Array(A)' \ B over the element-type matrix (real/complex factorization × Int/real/complex rhs, vector and matrix), the Transpose and A' \ b paths, coverage across all SPQR.ORDERINGS (which exercises the empty-cpiv ORDERING_FIXED case), and the dimension-mismatch errors.

🤖 Generated with Claude Code

https://claude.ai/code/session_01CSC6hQqcWDvAuvdkXJstBM

`qr(A)' \ b` and `ldiv!(x, qr(A)', b)` now solve the underdetermined
system `A'x = b` for a tall `A`, returning the minimum-norm solution.
With `A[prow, pcol] == Q*R` we have `A' == Pcol*R'*Q'*Prow`, so the
solve is a forward substitution with `R'` followed by a multiplication
by `Q`, reusing the same permutations, workspace and lock as the
existing `ldiv!`. Free variables are zeroed to select the minimum-norm
solution, which also drops the equations that the leading rank block of
`R` cannot represent when `A` is rank deficient.

Wide `A` is rejected with the same `DimensionMismatch` that dense QR
throws, since solving the resulting overdetermined system needs a
factorization of `A'` rather than of `A`.

This also fixes `A' \ b` and `transpose(A) \ b` for non-square sparse
`A`, which routed through the missing method and threw a `MethodError`.

Fixes #656.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CSC6hQqcWDvAuvdkXJstBM
@codecov

codecov Bot commented Sep 9, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 92.42%. Comparing base (f4b0d50) to head (c62f30f).

Additional details and impacted files
@@            Coverage Diff             @@
##             main     #804      +/-   ##
==========================================
+ Coverage   92.36%   92.42%   +0.05%     
==========================================
  Files          12       12              
  Lines        8632     8682      +50     
==========================================
+ Hits         7973     8024      +51     
+ Misses        659      658       -1     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

ViralBShah and others added 2 commits September 11, 2026 10:12
`\(::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, B)` and the `Transpose`
variant routed every non-square `A` through `adjoint(qr(A)) \ B`, which
rejects a wide `A` since solving `A'x = b` is then an overdetermined
problem needing a factorization of `A'`. Factorize the (cheap) sparse
transpose in that case, so `A' \ b` now works for every shape: tall `A`
keeps the minimum-norm solution via the factorization of `A`, wide `A`
gets the least squares solution like the dense path.

Also drop the dead `rpivinv` branch and the redundant view in the
adjoint `ldiv!`, matching the forward solve, and test the wide case.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Adhxb2mGC8BXYRA5CwgAwQ
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Adhxb2mGC8BXYRA5CwgAwQ
@ViralBShah

Copy link
Copy Markdown
Member Author

This change was accidentally included in the squash merge of #810. #817 reverts that portion from main; once it lands, this PR merges cleanly against main again and can be reviewed as intended.

ViralBShah added a commit that referenced this pull request Sep 11, 2026
#810 was meant to only bump the version and julia compat to 1.14, but it
was accidentally based on the branch of #804 and squash-merged with the
"Solve with the adjoint of a sparse QR factorization" change included
(see
#810 (comment)).
This reverts the `src/solvers/spqr.jl` and `test/spqr.jl` portion of
a248d20 so that #804 can be reviewed and merged on its own. The
`Project.toml` bump is kept.

Both files are restored byte-for-byte to their state before a248d20, and
`test/spqr.jl` passes locally on the 1.14-DEV build.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

https://claude.ai/code/session_01RFzXKGsJuSogARQpDGipPr

Co-authored-by: Claude Fable 5.1 <noreply@anthropic.com>
ViralBShah and others added 5 commits September 11, 2026 11:55
The copy dates from the first version of #676, which had no lock and copied
`F` in `\` to keep the non-mutating solve thread safe. The lock was added
before that PR merged but the copy stayed, so every `\` allocated a new lock
and workspace for nothing. `ldiv!` takes `F._lock` around all workspace
access, which is the same model UMFPACK and CHOLMOD use, so both `\` methods
now solve with the caller's factorization directly.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SgWK99c6dYwt3gxH44MxCs
The gather, scatter and zeroing loops are linear in the problem size while
the solve is dominated by the Q multiplication and the triangular solve, so
eliding the bounds checks is not measurable: within noise from 2000x500 up,
and about 40ns on a 100x10 problem. Every index is already guaranteed in
range by the dimension checks and the SPQR permutations, so the annotations
only served to turn a BoundsError into memory corruption for a malformed
QRSparse.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SgWK99c6dYwt3gxH44MxCs
@ViralBShah
ViralBShah merged commit 1b358c2 into main Sep 16, 2026
10 checks passed
@ViralBShah
ViralBShah deleted the vs/spqr-adjoint-ldiv branch September 16, 2026 21:01
ViralBShah added a commit that referenced this pull request Sep 17, 2026
…e adjoint (#833)

Fixes #114.

```julia
julia> A = sparse([1.0 0 1 0; 0 1 0 1]);

julia> lq(A)
ERROR: MethodError: no method matching lq!(::SparseMatrixCSC{Float64, Int64})
```

`lq(A)` now returns `qr(A')'`, the `AdjointQRSparse` from #804, computed
by SPQR on the sparse transpose, with `tol` and `ordering` passed
through. It gets `L`, `Q`, `prow` and `pcol` properties so that
`A[F.prow, F.pcol] == F.L * F.Q`, with `F.L` a lower triangular
`SparseMatrixCSC` and `F.Q` the adjoint of the sparse `Q`; `F \ b` is
the minimum-norm solution for a wide `A` and throws `DimensionMismatch`
for a tall one, as dense `lq` does; `F'` is `qr(A')`, `lq(A')` reuses
`qr(A)` without a copy, and `rank` and `show` work. `qr` also accepts
adjoint and transpose wrappers of a `SparseMatrixCSC` directly, and the
`\` for a tall `A'` uses that instead of spelling the copy. The identity
`F.L * F.Q` exposed a bug in the padded products with a sparse `Q`: they
took the thin width from the number of stored Householder vectors, which
SPQR keeps below the column count when the matrix is sparse enough, so
`qr(A).Q * qr(A).R` threw `DimensionMismatch` for such an `A`. The width
now comes from the column count that `QRSparseQ` records.

Tests in `test/spqr.jl` cover the factorization identity, properties,
solves, tolerance, the tall error, the adjoint and transpose inputs and
the thin product regression, for real and complex eltypes and both index
types; `spqr`, `linalg`, Aqua, whitespace and the doctests pass on
1.14-DEV. `docs/src/solvers.md` lists `lq`. #831 will use `lq` for the
wide `\`.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-authored-by: Claude Fable 5.1 <noreply@anthropic.com>
ViralBShah added a commit that referenced this pull request Sep 17, 2026
#831)

Fixes #301. Stacked on #833, which adds the sparse `lq`.

`A \ b` for a wide sparse `A` returned a basic solution, while dense `\`
returns the minimum-norm one:

```julia
julia> A = sparse([1.0 0 1 0; 0 1 0 1]); b = [1.0, 2.0];

julia> A \ b
4-element Vector{Float64}:
 1.0
 2.0
 0.0
 0.0

julia> Matrix(A) \ b
4-element Vector{Float64}:
 0.5
 1.0
 0.5
 1.0
```

The wide branch of `\` called `qr(A) \ B`, which SPQR solves by
back-substituting with the leading `rank` columns of `R`. It now calls
`lq(A) \ B`, the minimum-norm solve through the adjoint of the QR
factorization of `A'` from #804, as the `Adjoint`/`Transpose` methods
already do for a tall parent. `factorize(A)` returns `lq(A)` for a wide
`A`, so `factorize(A) \ b == A \ b`. `qr(A) \ b` still returns the basic
solution and its docstring says so.

Tests added next to the existing underdetermined check in the SPQR `\`
testset; `spqr`, `linalg`, `linalg_solvers`, whitespace and the doctests
pass on 1.14-DEV.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-authored-by: Claude Fable 5.1 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

missing ldiv for transpose of sparse QR

1 participant