From 82a0880ffa6558501584a489d26c9ad17c3525c7 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 23:03:44 +0000 Subject: [PATCH 1/2] Return the same float type from sparse `cond` for singular and nonsingular matrices Co-Authored-By: Claude Opus 5.5 --- src/linalg.jl | 8 +++++--- test/solvers/issues.jl | 2 ++ 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 446f071d..caf696c3 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -988,8 +988,10 @@ default, throws an `ArgumentError`, so `cond(A)` does too: use `cond(Array(A))` small enough to store dense. """ function cond(A::AbstractSparseMatrixCSC, p::Real=2) + # the type of dense `cond`, whichever precision the factorization of `A` works in + T = float(real(eltype(A))) # as for a dense matrix; the estimator of the inverse's norm needs at least one column - (p == 1 || p == Inf) && checksquare(A) == 0 && return zero(real(float(eltype(A)))) + (p == 1 || p == Inf) && checksquare(A) == 0 && return zero(T) if p == 1 || p == Inf # a singular matrix has an infinite condition number, as for a dense matrix; a # triangular or diagonal `factorize(A)` reports it only in the solve @@ -997,9 +999,9 @@ function cond(A::AbstractSparseMatrixCSC, p::Real=2) p == 1 ? opnormestinv(A) : opnormestinv(copy(A')) catch e e isa SingularException || rethrow() - return convert(float(real(eltype(A))), Inf) + return convert(T, Inf) end - return opnorm(A, p) * normAinv + return convert(T, opnorm(A, p) * normAinv) elseif p == 2 throw(ArgumentError("only 1- and Inf-norm condition numbers are implemented for sparse matrices, for 2-norm try cond(Array(A), 2) instead")) else diff --git a/test/solvers/issues.jl b/test/solvers/issues.jl index 60e094f9..0d50f795 100644 --- a/test/solvers/issues.jl +++ b/test/solvers/issues.jl @@ -258,6 +258,8 @@ end @test cond(sparse([1.0 2; 3 6]), 1) === Inf # the triangular solve, not the factorization, finds this one singular @test cond(sparse([1.0 0; 2 0]), Inf) === Inf + # the LU works in double precision, but the result is of the matrix's precision + @test cond(sparse(Float32[1 2; 3 4]), 1) isa Float32 end end # module From ae3041dd293065b8dad045ca7b14b43009ae6671 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 20:13:04 -0400 Subject: [PATCH 2/2] Catch only the factorization's `SingularException` in sparse `cond` (#1132) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Stacked on #1131. Since #1126 sparse `cond` wraps the whole `opnormestinv` estimator in a `try`, so a `SingularException` from anywhere inside it becomes `Inf`, and a singular triangular or diagonal matrix is only found once the estimator has allocated its blocks and run its first solve. Raised in the post-merge review of #1126. **Before** (#1131), `cond(A, 1)` of a singular `A` with `n = 10^6` and one zero on the diagonal | `A` | time | allocated | |---|---|---| | diagonal | 64.6 ms | 125.3 MiB | | lower bidiagonal | 16.1 ms | 68.7 MiB | **After** | `A` | time | allocated | |---|---|---| | diagonal | 9.4 ms | 56.6 MiB | | lower bidiagonal | 2.6 ms | 0 | Both return `Inf`, before and after. `factorize` throws a `SingularException` for a singular LU, and for a singular Hermitian matrix once Cholesky and LDLt have failed and it falls back to the LU (#1116); only the triangular and diagonal factors it returns leave the zero pivot to the solve. So `cond` now catches around `factorize` alone and checks the diagonal of those factors with `_issingularfactor`, which reuses the triangular solve's `_trisingular`, before it calls the estimator; `opnormestinv(A, t)` keeps its behaviour and calls the new `_opnormestinv(A, F, t)` with the factorization. A sweep of `cond` against dense over `Int`, `Float32`, `Float64`, `ComplexF32`, `ComplexF64`, `Rational{Int}`, `BigFloat` and `Complex{BigFloat}`, on general, singular, triangular, diagonal and Hermitian matrices, gives the same results as #1131 with a seeded RNG. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 --- src/linalg.jl | 30 +++++++++++++++++++++--------- test/solvers/issues.jl | 5 ++++- 2 files changed, 25 insertions(+), 10 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index caf696c3..9caf66d0 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -993,15 +993,16 @@ function cond(A::AbstractSparseMatrixCSC, p::Real=2) # as for a dense matrix; the estimator of the inverse's norm needs at least one column (p == 1 || p == Inf) && checksquare(A) == 0 && return zero(T) if p == 1 || p == Inf - # a singular matrix has an infinite condition number, as for a dense matrix; a - # triangular or diagonal `factorize(A)` reports it only in the solve - normAinv = try - p == 1 ? opnormestinv(A) : opnormestinv(copy(A')) + B = p == 1 ? A : copy(A') + # a singular matrix has an infinite condition number, as for a dense matrix + F = try + factorize(B) catch e e isa SingularException || rethrow() return convert(T, Inf) end - return convert(T, opnorm(A, p) * normAinv) + _issingularfactor(F) && return convert(T, Inf) + return convert(T, opnorm(A, p) * _opnormestinv(B, F, min(2, size(B, 1)))) elseif p == 2 throw(ArgumentError("only 1- and Inf-norm condition numbers are implemented for sparse matrices, for 2-norm try cond(Array(A), 2) instead")) else @@ -1009,11 +1010,22 @@ function cond(A::AbstractSparseMatrixCSC, p::Real=2) end end -function opnormestinv(A::AbstractSparseMatrixCSC{T}, t::Integer = min(2,maximum(size(A)))) where T +# `factorize` throws for a singular LU or Hermitian factorization, but leaves a triangular +# or diagonal one to the solve +_issingularfactor(F) = false +_issingularfactor(D::Diagonal) = any(iszero, D.diag) +_issingularfactor(L::LowerTriangular{<:Any,<:AbstractSparseMatrixCSC}) = !iszero(_trisingular(false, parent(L))) +_issingularfactor(U::UpperTriangular{<:Any,<:AbstractSparseMatrixCSC}) = !iszero(_trisingular(true, parent(U))) + +function opnormestinv(A::AbstractSparseMatrixCSC, t::Integer = min(2,maximum(size(A)))) + checksquare(A) + return _opnormestinv(A, factorize(A), t) +end + +# the estimate for the factorization `F` of `A` +function _opnormestinv(A::AbstractSparseMatrixCSC{T}, F, t::Integer) where T maxiter = 5 - # Check the input - n = checksquare(A) - F = factorize(A) + n = size(A, 1) if t <= 0 throw(ArgumentError("number of blocks must be a positive integer")) end diff --git a/test/solvers/issues.jl b/test/solvers/issues.jl index 0d50f795..2eff6dfa 100644 --- a/test/solvers/issues.jl +++ b/test/solvers/issues.jl @@ -256,8 +256,11 @@ end @testset "cond of a singular matrix is Inf" begin @test cond(sparse([1.0 2; 3 6]), 1) === Inf - # the triangular solve, not the factorization, finds this one singular + # `factorize` leaves the zero on the diagonal of a triangular or diagonal factor @test cond(sparse([1.0 0; 2 0]), Inf) === Inf + @test cond(sparse([1.0 0; 0 0]), 1) === Inf + # Hermitian, so `factorize` tries Cholesky and LDLt before the LU finds it singular + @test cond(sparse([1.0 1; 1 1]), 1) === Inf # the LU works in double precision, but the result is of the matrix's precision @test cond(sparse(Float32[1 2; 3 4]), 1) isa Float32 end