Skip to content

Return Inf from sparse cond of a singular matrix instead of throwing - #1126

Merged
ViralBShah merged 1 commit into
mainfrom
vs/cond-singular-inf
Oct 26, 2026
Merged

ViralBShah merged 1 commit into
mainfrom
vs/cond-singular-inf

Conversation

@ViralBShah

Copy link
Copy Markdown
Member

cond(A, 1) and cond(A, Inf) of a singular sparse matrix throw SingularException. Dense cond returns Inf for these.

using SparseArrays, LinearAlgebra
for A in ([1.0 2; 3 6], [1.0 0; 2 0]), p in (1, Inf)
    println("cond(sparse($A), $p) = ", try cond(sparse(A), p) catch e; sprint(showerror, e) end)
end

Before

cond(sparse([1.0 2.0; 3.0 6.0]), 1) = SingularException: matrix is singular; factorization failed. Zero pivot found at index 0
cond(sparse([1.0 2.0; 3.0 6.0]), Inf) = SingularException: matrix is singular; factorization failed. Zero pivot found at index 0
cond(sparse([1.0 0.0; 2.0 0.0]), 1) = SingularException: matrix is singular; factorization failed. Zero pivot found at index 2
cond(sparse([1.0 0.0; 2.0 0.0]), Inf) = SingularException: matrix is singular; factorization failed. Zero pivot found at index 2

After

cond(sparse([1.0 2.0; 3.0 6.0]), 1) = Inf
cond(sparse([1.0 2.0; 3.0 6.0]), Inf) = Inf
cond(sparse([1.0 0.0; 2.0 0.0]), 1) = Inf
cond(sparse([1.0 0.0; 2.0 0.0]), Inf) = Inf

cond estimates the norm of the inverse with opnormestinv, which solves with factorize(A). A singular matrix makes lu throw. A triangular or diagonal matrix gets no factorization object, so its solve throws. issuccess on lu(A; check = false) would only cover the lu branch, so cond catches SingularException instead, as dense cond does. It then returns convert(float(real(eltype(A))), Inf), the same value and type as dense. A Float32 matrix gives Inf32. The inferred return types are unchanged.

A singular Hermitian matrix still throws ZeroPivotException from the LDLt in factorize. That exception is not caught, because on main a nonsingular indefinite matrix such as [0.0 1; 1 0] throws it too. #1116 makes that path fall back to lu, which throws SingularException. With both PRs a singular Hermitian matrix also gives Inf; I checked this on #1116's branch with this change applied. This PR does not depend on #1116.

🤖 Generated with Claude Code

@ViralBShah ViralBShah added backport 1.13 Change should be backported to release-1.13 backport 1.14 labels Oct 26, 2026
@codecov

codecov Bot commented Oct 26, 2026 •

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 94.70%. Comparing base (66a7b5c) to head (0b47749).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files
@@            Coverage Diff             @@
##             main    #1126      +/-   ##
==========================================
- Coverage   94.70%   94.70%   -0.01%     
==========================================
  Files          16       16              
  Lines       10387    10385       -2     
==========================================
- Hits         9837     9835       -2     
  Misses        550      550              

☔ 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
ViralBShah force-pushed the vs/cond-singular-inf branch from 2f901cc to cc8551a Compare October 26, 2026 14:58
`cond(A, 1)` and `cond(A, Inf)` caught no `SingularException` from the
factorization or the triangular solve in `opnormestinv`; dense `cond`
returns `Inf` for a singular matrix, and now so does sparse.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@ViralBShah
ViralBShah force-pushed the vs/cond-singular-inf branch from cc8551a to 0b47749 Compare October 26, 2026 21:58
@ViralBShah
ViralBShah enabled auto-merge (squash) October 26, 2026 21:58
@ViralBShah
ViralBShah merged commit 4c87861 into main Oct 26, 2026
10 checks passed
@ViralBShah
ViralBShah deleted the vs/cond-singular-inf branch October 26, 2026 22:08
ViralBShah added a commit that referenced this pull request Oct 27, 2026
#1126 made sparse `cond` return `Inf` for a singular matrix, but its
docstring does not say so, nor that rounding can hide the singularity.
Raised in the post-merge review of #1126.

The material is the behaviour on `main`:

```julia
julia> cond(sparse([1.0 2; 3 6]), 1)
Inf

julia> A = [1.0 2 3; 4 5 6; 7 8 9];   # singular

julia> cond(sparse(A), 1), cond(A, 1)
(1.0808639105689187e17, Inf)
```

The UMFPACK LU of the second matrix ends with a pivot that rounding
leaves nonzero, while LAPACK's partial pivoting meets an exact zero, so
the docstring says that a finite but very large result can still mean a
singular matrix, and that sparse and dense can disagree on which case
they report.

Written by Claude Code.

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

Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
ViralBShah added a commit that referenced this pull request Oct 27, 2026
…1132)

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 <noreply@anthropic.com>
ViralBShah added a commit that referenced this pull request Oct 27, 2026
…gular matrices (#1131)

Sparse `cond` of a `Float32` or `ComplexF32` matrix returned `Float64`
when the matrix went through the LU and was nonsingular, but `Float32`
when it was singular (since #1126) or triangular, diagonal or positive
definite. Dense `cond` returns `Float32` in every case. Raised in the
post-merge review of #1126.

**Before**
```julia
julia> cond(sparse(Float32[1 2; 3 4]), 1)
20.99999999999999

julia> cond(sparse(Float32[1 2; 3 6]), 1)
Inf32

julia> Base.return_types(cond, (SparseMatrixCSC{Float32,Int}, Int))
1-element Vector{Any}:
 Union{Float32, Float64}
```

**After**
```julia
julia> cond(sparse(Float32[1 2; 3 4]), 1)
21.0f0

julia> cond(sparse(Float32[1 2; 3 6]), 1)
Inf32

julia> Base.return_types(cond, (SparseMatrixCSC{Float32,Int}, Int))
1-element Vector{Any}:
 Float32
```

UMFPACK works in double precision, so the estimate of the inverse's norm
is a `Float64` for any eltype that takes the LU (a nonsymmetric matrix,
or an indefinite Hermitian one where `factorize` falls back to `lu`),
and `opnorm(A, p) * normAinv` promoted to it. `cond` now computes `T =
float(real(eltype(A)))` once and returns `zero(T)`, `convert(T, Inf)` or
`convert(T, opnorm(A, p) * normAinv)`. A sweep against dense `cond` over
`Int`, `Float32`, `Float64`, `ComplexF32`, `ComplexF64`,
`Rational{Int}`, `BigFloat` and `Complex{BigFloat}`, on general,
singular, triangular, diagonal and Hermitian matrices, gives the dense
result type in every case that runs; only the `Float32` and `ComplexF32`
LU cases changed.

Potential breakage: `cond` of a nonsingular `Float32` or `ComplexF32`
sparse matrix that is not triangular, diagonal or positive definite now
returns `Float32` instead of `Float64`, as dense `cond` does.

Written by Claude Code.

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

---------

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

Labels

backport 1.13 Change should be backported to release-1.13 backport 1.14

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant