From 96b65ffb8e4c18f2fa144c1c4d249d7cd646d77f Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 13:21:13 +0000 Subject: [PATCH 1/2] Give sparse `\` the eltype of dense `\` for tall and wide matrices 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 --- src/linalg.jl | 38 ++++++++++++++++++++++++-------------- test/solvers/solvers.jl | 6 +++++- 2 files changed, 29 insertions(+), 15 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index f4feb7c2..29bbc685 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1436,10 +1436,11 @@ end _hermitian_solve(A::AbstractSparseMatrixCSC) = !(eltype(A) <: Union{Integer, Complex{<:Integer}}) && ishermitian(A) -# `\` and `factorize` choose the same method from the structure of `A`; `f` (or `flu`, for -# the LU factorization) is applied in each branch so that the result stays inferrable. A -# Hermitian matrix is only wrapped, because `\` and `factorize` finish it differently. -function _sparse_factorize(f, A::AbstractSparseMatrixCSC, flu = f) +# `\` and `factorize` choose the same method from the structure of `A`; `f` (or `fsolver`, +# for the LU, QR and LQ factorizations) is applied in each branch so that the result stays +# inferrable, and the QR and LQ factorizations are of `fqr(A)`. A Hermitian matrix is only +# wrapped, because `\` and `factorize` finish it differently. +function _sparse_factorize(f, A::AbstractSparseMatrixCSC, fsolver = f, fqr = identity) m, n = size(A) if m == n if istril(A) @@ -1449,12 +1450,12 @@ function _sparse_factorize(f, A::AbstractSparseMatrixCSC, flu = f) elseif _hermitian_solve(A) return f(Hermitian(A)) end - return flu(lu(A)) + return fsolver(lu(A)) elseif m > n - return f(qr(A)) + return fsolver(qr(fqr(A))) else # A is wide, so the LQ factorization gives the minimum-norm solution - return f(lq(A)) + return fsolver(lq(fqr(A))) end end @@ -1462,15 +1463,24 @@ end _solve_factor(F) = F _solve_factor(D::Diagonal) = Diagonal(Vector(diag(D))) _sparse_solve(op, F, B) = \(op(_solve_factor(F)), B) -# the sparse LU works in double precision, so its solution is converted to the eltype that -# dense `\` gives -_sparse_lusolve(op, F, A, B) = - convert(AbstractArray{typeof(one(eltype(A)) \ one(eltype(B)))}, _sparse_solve(op, F, B)) +# The eltype that dense `\` gives. It is `nothing` for an abstract eltype, which has no `one` +# and whose solution is returned as the solver gives it. +_ldiv_eltype(A, B) = isconcretetype(eltype(A)) && isconcretetype(eltype(B)) ? + typeof(one(eltype(A)) \ one(eltype(B))) : nothing +# The sparse LU works in double precision and the QR and LQ in the precision of the matrix, +# so their solution is converted to the eltype that dense `\` gives, and for the QR and LQ a +# matrix of a narrower eltype than that is converted before it is factorized. +_sparse_solversolve(op, F, A, B) = _convert_solution(_ldiv_eltype(A, B), _sparse_solve(op, F, B)) +_convert_solution(T, X) = convert(AbstractArray{T}, X) +_convert_solution(::Nothing, X) = X +_qr_matrix(A, B) = _qr_matrix(A, _ldiv_eltype(A, B)) +_qr_matrix(A, T::Type) = convert(AbstractMatrix{promote_type(float(eltype(A)), real(T))}, A) +_qr_matrix(A, ::Nothing) = A function \(A::AbstractSparseMatrixCSC, B::AbstractVecOrMat) require_one_based_indexing(A, B) return _sparse_factorize(F -> _sparse_solve(identity, F, B), A, - F -> _sparse_lusolve(identity, F, A, B)) + F -> _sparse_solversolve(identity, F, A, B), A -> _qr_matrix(A, B)) end for (xformtype, xformop) in ((:Adjoint, :adjoint), (:Transpose, :transpose)) @eval function \(xformA::($xformtype){<:Any,<:AbstractSparseMatrixCSC}, B::AbstractVecOrMat) @@ -1479,12 +1489,12 @@ for (xformtype, xformop) in ((:Adjoint, :adjoint), (:Transpose, :transpose)) m, n = size(A) if m < n # A' is tall, so the least squares solve needs a factorization of A' itself - return _sparse_solve(identity, qr($xformop(A)), B) + return _sparse_solversolve(identity, qr($xformop(_qr_matrix(A, B))), A, B) end # Otherwise the transformed choice for A serves. For a wide A' that is the # factorization of A, which gives the minimum-norm solution. return _sparse_factorize(F -> _sparse_solve($xformop, F, B), A, - F -> _sparse_lusolve($xformop, F, A, B)) + F -> _sparse_solversolve($xformop, F, A, B), A -> _qr_matrix(A, B)) end end diff --git a/test/solvers/solvers.jl b/test/solvers/solvers.jl index 7f763e38..57391d57 100644 --- a/test/solvers/solvers.jl +++ b/test/solvers/solvers.jl @@ -129,7 +129,11 @@ end # the adjoint solve converts its result to the eltype the plain solve gives x = square' \ Float32.(b) @test x isa Vector{Float32} && x ≈ Matrix(square') \ b - # only the LU solve converts its result, so a least squares solve takes any right-hand side + # the least squares solve gives the eltype of dense `\`, and keeps the precision of `b` + x = @inferred Float32.(tall) \ b + @test x isa Vector{Float64} && x ≈ Matrix(tall) \ b + # the solution for an abstract eltype is not converted, so a least squares solve takes + # any right-hand side @test tall \ Any[1.0, 2, 3] ≈ Matrix(tall) \ b end # when the Cholesky factorization of a Hermitian matrix fails, `factorize` returns From 74073172293d8146ec7f6e8d7c6c248ef0888b42 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 15:03:53 +0000 Subject: [PATCH 2/2] Give `qr(A) \ b` and transposed complex Hermitian `\` the eltype of dense `\` 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 --- src/linalg.jl | 2 +- src/solvers/spqr.jl | 14 ++++++++++++-- test/solvers/solvers.jl | 3 +++ test/solvers/spqr.jl | 2 ++ 4 files changed, 18 insertions(+), 3 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 29bbc685..03807a5b 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -554,7 +554,7 @@ function _dot(x::AbstractSparseVector, A::SparseMatrixCSCOrColumnSubset, y::Abst end ## end of symmetric/Hermitian -\(A::Transpose{<:Complex,<:Hermitian{<:Complex,<:AbstractSparseMatrixCSC}}, B::Vector) = copy(A) \ B +\(A::Transpose{<:Complex,<:Hermitian{<:Complex,<:AbstractSparseMatrixCSC}}, B::VecOrMat) = copy(A) \ B ## diagonal scaling diff --git a/src/solvers/spqr.jl b/src/solvers/spqr.jl index cbfd69ef..6c565685 100644 --- a/src/solvers/spqr.jl +++ b/src/solvers/spqr.jl @@ -654,6 +654,15 @@ LinearAlgebra.rank(F::AdjointQRSparse) = rank(parent(F)) _rhs_eltype(::Type{T}, ::Type{<:Complex}) where {T<:LinearAlgebra.BlasReal} = Complex{T} _rhs_eltype(::Type{T}, ::Type) where {T} = T _rhs_array(F, B::AbstractVecOrMat) = convert(Array{_rhs_eltype(eltype(F), eltype(B))}, B) +# Rather than round the right-hand side to the precision of the factorization, a solve whose +# solution dense `\` gives a wider eltype converts the factorization to that precision, as +# dense `\` does. A real factorization stays real for a complex right-hand side. +_solve_eltype(F, B) = isconcretetype(eltype(B)) ? + promote_type(eltype(F), real(typeof(zero(eltype(B)) / oneunit(eltype(F))))) : eltype(F) +_promoted_factor(F::QRSparse, B) = _promoted_factor(F, _solve_eltype(F, B)) +_promoted_factor(F::QRSparse{T}, ::Type{T}) where {T} = F +_promoted_factor(F::QRSparse, ::Type{T}) where {T<:LinearAlgebra.BlasFloat} = QRSparse{T}(F) +_promoted_factor(F::QRSparse, ::Type) = F ## Two helper methods _ret_size(F::Union{QRSparse,AdjointQRSparse}, b::AbstractVector) = (size(F, 2),) @@ -733,7 +742,7 @@ julia> qr(A)\\fill(1.0, 4) 0.0 ``` """ -(\)(F::QRSparse, B::AbstractVecOrMat) = F\_rhs_array(F, B) +(\)(F::QRSparse, B::AbstractVecOrMat) = (G = _promoted_factor(F, B); G\_rhs_array(G, B)) function LinearAlgebra.ldiv!(X::StridedVecOrMat{T}, F::QRSparse{T}, B::StridedVecOrMat{T}; workspace::Union{Nothing,SpqrWS{T}} = nothing) where {T} @@ -842,7 +851,8 @@ julia> A'x 0.0 ``` """ -(\)(Fadj::AdjointQRSparse, B::AbstractVecOrMat) = Fadj\_rhs_array(Fadj, B) +(\)(Fadj::AdjointQRSparse, B::AbstractVecOrMat) = + (G = adjoint(_promoted_factor(parent(Fadj), B)); G\_rhs_array(G, B)) function LinearAlgebra.ldiv!(X::StridedVecOrMat{T}, Fadj::AdjointQRSparse{T}, B::StridedVecOrMat{T}; workspace::Union{Nothing,SpqrWS{T}} = nothing) where {T} diff --git a/test/solvers/solvers.jl b/test/solvers/solvers.jl index 57391d57..21cbb8d1 100644 --- a/test/solvers/solvers.jl +++ b/test/solvers/solvers.jl @@ -132,6 +132,9 @@ end # the least squares solve gives the eltype of dense `\`, and keeps the precision of `b` x = @inferred Float32.(tall) \ b @test x isa Vector{Float64} && x ≈ Matrix(tall) \ b + # so does the solve of a transposed complex Hermitian matrix with a matrix right-hand side + Ht = transpose(sparse(ComplexF32[4 im; -im 3])) + @test @inferred(Ht \ Float32[1 2; 3 4]) isa Matrix{ComplexF32} # the solution for an abstract eltype is not converted, so a least squares solve takes # any right-hand side @test tall \ Any[1.0, 2, 3] ≈ Matrix(tall) \ b diff --git a/test/solvers/spqr.jl b/test/solvers/spqr.jl index 2417d4a4..8adbd6ca 100644 --- a/test/solvers/spqr.jl +++ b/test/solvers/spqr.jl @@ -254,6 +254,8 @@ end b, B = collect(1.0:m), reshape(collect(1.0:3m), 3, m) eltyA <: Real && @test F.Q * b ≈ F.Q * eltyA.(b) eltyA <: Real && @test B * F.Q' ≈ eltyA.(B) * F.Q' + # a solve converts F to the precision of the solution, as dense `\` does + @test @inferred(F \ b) isa Vector{widen(eltyA)} end end