Skip to content

Give sparse \ and qr(A) \ b the eltype of dense \ - #1124

Open
ViralBShah wants to merge 2 commits into
mainfrom
vs/backslash-eltype-shape
Open

ViralBShah wants to merge 2 commits into
mainfrom
vs/backslash-eltype-shape

Conversation

@ViralBShah

@ViralBShah ViralBShah commented Oct 26, 2026 •

Copy link
Copy Markdown
Member

The eltype of sparse A \ B depended on the shape of A: a tall or wide Float32 matrix with a Float64 right-hand side gave Float32, where a square one and dense \ give Float64, and an Int tall matrix with a Float32 right-hand side gave Float64 where dense gives Float32. So \ was not inferable, and the least squares solve rounded b to Float32. The same rounding happened at the factorization level in qr(A) \ b and qr(A)' \ b, and transpose(A) \ B of a complex Hermitian ComplexF32 matrix with a matrix B gave ComplexF64.

using SparseArrays, LinearAlgebra
for A in (Float32[1 0; 0 1; 1 1], Float32[2 0; 1 1], Float32[1 0 1; 0 1 1])
    b = collect(1.0:size(A, 1))
    println(size(A), ": sparse ", eltype(sparse(A) \ b), ", dense ", eltype(A \ b))
end
println(Base.infer_return_type(\, Tuple{SparseMatrixCSC{Float32,Int}, Vector{Float64}}))
A = sparse(Float32[1 0; 0 1; 1 1]); b = [1.0, 2.0, 3.0 + 1e-7]
x = A \ b
println(norm(Matrix{Float64}(A)' * (Matrix{Float64}(A) * x - b)))   # normal-equations residual

Before

(3, 2): sparse Float32, dense Float64
(2, 2): sparse Float64, dense Float64
(2, 3): sparse Float32, dense Float64
Union{Vector{Float32}, Vector{Float64}}
3.6563210698207683e-7

After

(3, 2): sparse Float64, dense Float64
(2, 2): sparse Float64, dense Float64
(2, 3): sparse Float64, dense Float64
Vector{Float64}
1.1102230246251565e-15
using SparseArrays, LinearAlgebra
A = sparse(ComplexF32[2 1im; -1im 2]); B = Float32[1 3; 2 4]
println("transpose(A) \\ B: sparse ", eltype(transpose(A) \ B), ", dense ", eltype(transpose(Matrix(A)) \ B))
println(Base.infer_return_type(\, Tuple{typeof(transpose(A)), Matrix{Float32}}))
T = sparse(Float32[1 0; 0 1; 1 1]); b = [1.0, 2.0, 3.0 + 1e-7]
println("qr(A) \\ b: sparse ", eltype(qr(T) \ b), ", dense ", eltype(qr(Matrix(T)) \ b))
println("qr(A)' \\ b: sparse ", eltype(qr(T)' \ b[1:2]), ", dense ", eltype(qr(Matrix(T))' \ b[1:2]))

Before

transpose(A) \ B: sparse ComplexF64, dense ComplexF32
Union{Matrix{ComplexF32}, Matrix{ComplexF64}}
qr(A) \ b: sparse Float32, dense Float64
qr(A)' \ b: sparse Float32, dense Float64

After

transpose(A) \ B: sparse ComplexF32, dense ComplexF32
Matrix{ComplexF32}
qr(A) \ b: sparse Float64, dense Float64
qr(A)' \ b: sparse Float64, dense Float64

Only the LU branch of _sparse_factorize converted its solution to typeof(one(eltype(A)) \ one(eltype(B))). The QR and LQ branches, including the qr(A') branch of the adjoint and transpose methods, returned the eltype of the factor, which SPQR gives the eltype of the matrix. They now convert the same way. Converting the result alone would keep the Float32 solve, since qr of a Float32 matrix factorizes in double precision and rounds the factor back, and its solve converts B to the factor's eltype. So before the QR or LQ factorization, a matrix of a narrower eltype than the solution is converted to the wider one. That costs nothing extra, because qr makes that copy anyway. A matrix and right-hand side of the same eltype factorize as before. The conversion is skipped for an abstract eltype, which has no one, so tall \ Any[...] still works.

At the factorization level, qr(A) \ b and qr(A)' \ b (and transpose(qr(A)) \ b, which goes through the adjoint) converted b to the eltype of the factor. When dense \ gives a wider solution, they now convert the factorization to that precision with the existing QRSparse{T}(F) constructor, as dense \ does with Factorization{T}(F). This copies the factor's arrays exactly and does not refactorize, so the accuracy is that of the Float32 factorization, as for dense. A real factorization stays real for a complex right-hand side and keeps the reinterpreting solve. The factor is not converted when the target is not a BLAS float, such as BigFloat, or when b has an abstract eltype; those cases behave as before. A \ b does not take this path, because its factor already has the solution's precision.

The transposed complex Hermitian \ method was defined for a Vector only, so a matrix right-hand side went to LinearAlgebra's generic \. It now takes VecOrMat, and the matrix goes through the sparse Hermitian solve as the vector does. This touches dispatch.

Potential breakage: code that works today gets a different result eltype. A tall or wide Float32 matrix with a Float64 right-hand side returned Vector{Float32} from A \ b, qr(A) \ b and qr(A)' \ b, and now returns Vector{Float64}, with the solve done in the promoted precision. transpose(A) \ B of a complex Hermitian ComplexF32 matrix with a matrix B returned Matrix{ComplexF64} and now returns Matrix{ComplexF32}.

Left out: at the factorization level, lu and cholesky still differ from dense. lu of a Float32 or ComplexF32 matrix is an UmfpackLU{Float64} or UmfpackLU{ComplexF64}, so lu(A) \ b returns double precision where dense keeps single. The factor does not record the original eltype, so this cannot be fixed by converting the solution. cholesky(A) \ b of a Float32 matrix keeps a single-precision CHOLMOD factor and returns Float32 for a Float64 b, where dense returns Float64. Converting a CHOLMOD factor is not the same small change. A \ b is not affected by either. #999 changes the sparse right-hand side \ beside this code. Its rectangular branch calls the dense right-hand side \, so it gets this fix, and the diffs touch different lines.

🤖 Generated with Claude Code

@codecov

codecov Bot commented Oct 26, 2026 •

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 91.30435% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.74%. Comparing base (e4a42b4) to head (7407317).
⚠️ Report is 4 commits behind head on main.

Files with missing lines Patch % Lines
src/solvers/spqr.jl 75.00% 2 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main    #1124      +/-   ##
==========================================
+ Coverage   94.73%   94.74%   +0.01%     
==========================================
  Files          16       16              
  Lines       10325    10334       +9     
==========================================
+ Hits         9781     9791      +10     
+ Misses        544      543       -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 October 26, 2026 14:55
The QR and LQ branches of `A \ B` returned the eltype of the factor, so the
result type depended on the shape of `A` and was not inferable. Convert their
solution as the LU branch does, and factorize a matrix of a narrower eltype
than the solution in the wider one so that the solve keeps the precision of `B`.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ense `\`

A QR factorization of a narrower eltype than the solution is converted to the
wider one before the solve, as dense `\` converts its factorization, rather
than rounding the right-hand side. The transposed complex Hermitian `\` takes
a matrix right-hand side as well as a vector.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@ViralBShah ViralBShah changed the title Give sparse \ the eltype of dense \ for tall and wide matrices Give sparse \ and qr(A) \ b the eltype of dense \ Oct 26, 2026
@ViralBShah
ViralBShah force-pushed the vs/backslash-eltype-shape branch from 97ad8e7 to 7407317 Compare October 26, 2026 15:04

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant