From 09be1ea7cf437616551c4511a807e47d2085336a Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Tue, 8 Sep 2026 10:47:30 +0000 Subject: [PATCH 1/3] Expose CHOLMOD's reciprocal condition number estimate as `CHOLMOD.rcond` `cholmod_rcond` was already wrapped in `wrappers.jl` but unreachable from Julia. Add a thin `rcond(F::Factor)` in the index-type loop and document what the estimate actually is: min/max of the factor diagonal, squared for `LL'`, so that `cholesky` and `ldlt` of the same matrix agree. It is exact for diagonal matrices and an upper bound on `1 / cond(A, 2)` otherwise -- verified over 300 random SPD matrices, where it never dipped below the true value and was up to 91x optimistic. Documented as a way to detect a badly conditioned or singular factorization rather than to measure conditioning. Fixes #118. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01L8HenXTVrASo1sBZVGhpDQ --- docs/src/solvers.md | 1 + src/solvers/cholmod.jl | 39 +++++++++++++++++++++++++++++++++++++++ test/cholmod_ops.jl | 17 ++++++++++++++++- 3 files changed, 56 insertions(+), 1 deletion(-) diff --git a/docs/src/solvers.md b/docs/src/solvers.md index c2624b81..2d6830bc 100644 --- a/docs/src/solvers.md +++ b/docs/src/solvers.md @@ -29,6 +29,7 @@ SparseArrays.CHOLMOD.lowrankdowndate SparseArrays.CHOLMOD.lowrankdowndate! SparseArrays.CHOLMOD.lowrankupdowndate! SparseArrays.CHOLMOD.ldlt +SparseArrays.CHOLMOD.rcond SparseArrays.SPQR.qr SparseArrays.UMFPACK.lu ``` diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 7f18b509..7f51bba4 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -29,6 +29,8 @@ export Factor, Sparse +public rcond + import SparseArrays: AbstractSparseMatrix, SparseMatrixCSC, indtype, sparse, spzeros, nnz, sparsevec @@ -547,6 +549,9 @@ for TI ∈ IndexTypes function check_factor(F::Factor{Tv, $TI}) where Tv<:VTypes $(cholname(:check_factor, TI))(F, getcommon($TI)) != 0 end + function rcond(F::Factor{Tv, $TI}) where Tv<:VTypes + $(cholname(:rcond, TI))(F, getcommon($TI)) + end nnz(A::Sparse{<:VTypes, $TI}) = $(cholname(:nnz, TI))(A, getcommon($TI)) function speye(m::Integer, n::Integer, ::Type{Tv}, ::Type{$TI}) where Tv<:VTypes @@ -2037,6 +2042,40 @@ end det(L::Factor) = exp(logdet(L)) +""" + rcond(F::CHOLMOD.Factor) -> Float64 + +Return CHOLMOD's rough estimate of the reciprocal condition number of the +factorized matrix, computed from the diagonal of the factor alone: the smallest +entry of `abs.(diag(F))` divided by the largest, squared when `F` is an `LL'` +factorization so that the result estimates the reciprocal condition number of +the factorized matrix rather than of its factor. + +This is much cheaper than a norm-based estimate such as `cond(A, 1)`, but also +much cruder. For positive definite `A` it is exact when `A` is diagonal, and +otherwise an upper bound on `1 / cond(A, 2)`, so it can report a matrix as far +better conditioned than it is. Use it to detect a badly conditioned or singular +factorization, not to measure conditioning accurately. + +Returns `0` if the matrix is singular or the factor has a zero or `NaN` on its +diagonal, and `1` if the matrix is 1-by-1. `NaN` is never returned. + +# Examples +```jldoctest +julia> A = sparse(Diagonal([1.0, 2.0, 4.0])); + +julia> SparseArrays.CHOLMOD.rcond(cholesky(A)) +0.25 + +julia> SparseArrays.CHOLMOD.rcond(ldlt(A)) +0.25 + +julia> SparseArrays.CHOLMOD.rcond(cholesky(sparse(Diagonal([1.0, 0.0])); check=false)) +0.0 +``` +""" +rcond + function issuccess(F::Factor) s = unsafe_load(pointer(F)) return s.minor == size(F, 1) diff --git a/test/cholmod_ops.jl b/test/cholmod_ops.jl index acad5e7d..ca252c8d 100644 --- a/test/cholmod_ops.jl +++ b/test/cholmod_ops.jl @@ -14,7 +14,7 @@ using SparseArrays.CHOLMOD: getcommon using Random using Serialization using LinearAlgebra: - I, cholesky, cholesky!, det, diag, eigmax, ishermitian, isposdef, issuccess, + I, cholesky, cholesky!, cond, det, diag, eigmax, ishermitian, isposdef, issuccess, issymmetric, ldiv!, ldlt, ldlt!, logdet, norm, opnorm, Diagonal, Hermitian, Symmetric, PosDefException, ZeroPivotException, RowMaximum using SparseArrays @@ -710,6 +710,21 @@ end @test Array(U) == Tv[20 0 0; 0 30 0; 10 0 0] end +@testset "rcond (#118)" begin + D = SparseMatrixCSC{Tv,Ti}(sparse(Diagonal(Tv[1, 2, 4]))) + # exact for a diagonal matrix, and the same estimate from LL' and LDL' + @test CHOLMOD.rcond(cholesky(D)) === 0.25 + @test CHOLMOD.rcond(ldlt(D)) === 0.25 + # 1-by-1 and singular special cases + @test CHOLMOD.rcond(cholesky(SparseMatrixCSC{Tv,Ti}(sparse(Diagonal(Tv[3]))))) === 1.0 + S = SparseMatrixCSC{Tv,Ti}(sparse(Diagonal(Tv[1, 0]))) + @test CHOLMOD.rcond(cholesky(S; check=false)) === 0.0 + # the estimate never reports a matrix as worse conditioned than it is + B = sprandn(Tv, 20, 20, 0.4) + C = SparseMatrixCSC{Tv,Ti}(B*B' + 20I) + @test CHOLMOD.rcond(cholesky(C)) >= 1/cond(Matrix{Float64}(C), 2) - sqrt(eps(Tv)) +end + end # for Tv ∈ (Float32, Float64) end # Base.USE_GPL_LIBS From 62d4b9769208ece1c15ec664dfe33d13e23f0684 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Wed, 9 Sep 2026 05:14:50 -0400 Subject: [PATCH 2/3] Expose UMFPACK's reciprocal condition number estimate as `UMFPACK.rcond` UMFPACK writes min(abs(diag(U))) / max(abs(diag(U))) into `Info[UMFPACK_RCOND]` during the numeric factorization, but nothing read it back. Add `rcond(F::UmfpackLU)` as the LU counterpart of `CHOLMOD.rcond`, with a docstring, `public` declaration, docs entry and tests. Two things worth noting in the docstring: UMFPACK row-scales the matrix before factorizing by default, so the estimate describes the scaled matrix and every diagonal matrix reports 1; and unlike the Cholesky version it is not a bound on 1 / cond(A, 2) in either direction. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_011BA2YkbzRZcTbTrc3S8odn --- docs/src/solvers.md | 1 + src/solvers/cholmod.jl | 3 ++- src/solvers/umfpack.jl | 41 +++++++++++++++++++++++++++++++++++++++++ test/umfpack.jl | 20 +++++++++++++++++++- 4 files changed, 63 insertions(+), 2 deletions(-) diff --git a/docs/src/solvers.md b/docs/src/solvers.md index 2d6830bc..dce8f4a1 100644 --- a/docs/src/solvers.md +++ b/docs/src/solvers.md @@ -32,6 +32,7 @@ SparseArrays.CHOLMOD.ldlt SparseArrays.CHOLMOD.rcond SparseArrays.SPQR.qr SparseArrays.UMFPACK.lu +SparseArrays.UMFPACK.rcond ``` ```@meta diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 7f51bba4..1f1bd462 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -2055,7 +2055,8 @@ This is much cheaper than a norm-based estimate such as `cond(A, 1)`, but also much cruder. For positive definite `A` it is exact when `A` is diagonal, and otherwise an upper bound on `1 / cond(A, 2)`, so it can report a matrix as far better conditioned than it is. Use it to detect a badly conditioned or singular -factorization, not to measure conditioning accurately. +factorization, not to measure conditioning accurately. The LU counterpart is +[`UMFPACK.rcond`](@ref SparseArrays.UMFPACK.rcond). Returns `0` if the matrix is singular or the factor has a zero or `NaN` on its diagonal, and `1` if the matrix is 1-by-1. `NaN` is never returned. diff --git a/src/solvers/umfpack.jl b/src/solvers/umfpack.jl index 9e1836ef..5abbc53d 100644 --- a/src/solvers/umfpack.jl +++ b/src/solvers/umfpack.jl @@ -4,6 +4,8 @@ module UMFPACK export UmfpackLU +public rcond + import Base: (\), getproperty, show, size using LinearAlgebra using LinearAlgebra: AdjOrTrans @@ -43,6 +45,8 @@ import ..LibSuiteSparse: ## Sizes of Control and Info arrays for returning information from solver UMFPACK_INFO, UMFPACK_CONTROL, + # index of the info array in ZERO BASED indexing + UMFPACK_RCOND, # index of the control arrays in ZERO BASED indexing UMFPACK_PRL, UMFPACK_DENSE_ROW, @@ -95,6 +99,7 @@ const JL_UMFPACK_SCALE = UMFPACK_SCALE + 1 const JL_UMFPACK_FRONT_ALLOC_INIT = UMFPACK_FRONT_ALLOC_INIT + 1 const JL_UMFPACK_DROPTOL = UMFPACK_DROPTOL + 1 const JL_UMFPACK_IRSTEP = UMFPACK_IRSTEP + 1 +const JL_UMFPACK_RCOND = UMFPACK_RCOND + 1 struct MatrixIllConditionedException <: Exception msg::String @@ -928,6 +933,42 @@ end LinearAlgebra.issuccess(lu::UmfpackLU) = lu.status == UMFPACK_OK +""" + rcond(F::UmfpackLU) -> Float64 + +Return UMFPACK's rough estimate of the reciprocal condition number of the +factorized matrix, computed from the diagonal of the factor alone: the smallest +entry of `abs.(diag(F.U))` divided by the largest. + +This is much cheaper than a norm-based estimate such as `cond(A, 1)`, but also +much cruder, and it describes the matrix UMFPACK actually factorized rather +than `A` itself. UMFPACK scales the rows of `A` before factorizing by default +(see `F.Rs`), so for instance every diagonal matrix reports `1`. Unlike the +Cholesky-based [`CHOLMOD.rcond`](@ref SparseArrays.CHOLMOD.rcond), the value +is neither an upper nor a lower bound on `1 / cond(A, 2)`. Use it to detect a +singular or badly pivoted factorization, not to measure conditioning. + +Returns `0` if the matrix is singular, and `1` if the matrix is 1-by-1. + +# Examples +```jldoctest +julia> F = lu(sparse([1.0 3.0; 0.0 1.0])); + +julia> SparseArrays.UMFPACK.rcond(F) +0.25 + +julia> minimum(abs, diag(F.U)) / maximum(abs, diag(F.U)) +0.25 + +julia> SparseArrays.UMFPACK.rcond(lu(sparse([1.0 2.0; 0.0 0.0]); check=false)) +0.0 +``` +""" +function rcond(F::UmfpackLU) + umfpack_numeric!(F) # ensure the numeric decomposition exists + return F.info[JL_UMFPACK_RCOND] +end + ### Solve with Factorization ldiv!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = diff --git a/test/umfpack.jl b/test/umfpack.jl index 0a4525e2..b79aede6 100644 --- a/test/umfpack.jl +++ b/test/umfpack.jl @@ -11,7 +11,7 @@ using Random using SparseArrays using Serialization using LinearAlgebra: - LinearAlgebra, I, det, issuccess, ldiv!, lu, lu!, Transpose, SingularException, Diagonal, logabsdet + LinearAlgebra, I, det, diag, issuccess, ldiv!, lu, lu!, Transpose, SingularException, Diagonal, logabsdet using SparseArrays: nnz, sparse, sprand, sprandn, SparseMatrixCSC, UMFPACK, increment! function umfpack_report(l::UMFPACK.UmfpackLU) @@ -369,6 +369,24 @@ end end end + @testset "rcond (#118) for $Tv, $Ti" for Tv in (Float64, ComplexF64), Ti in (Int32, Int64) + # the number is min/max of |diag(U)| of the row-scaled matrix UMFPACK factorized + F = lu(SparseMatrixCSC{Tv,Ti}(sparse(Tv[1 3; 0 1]))) + @test UMFPACK.rcond(F) === 0.25 + @test UMFPACK.rcond(F) === minimum(abs, diag(F.U)) / maximum(abs, diag(F.U)) + # row scaling is on by default, so a diagonal matrix is perfectly conditioned + @test UMFPACK.rcond(lu(SparseMatrixCSC{Tv,Ti}(sparse(Diagonal(Tv[1, 2, 4]))))) === 1.0 + # 1-by-1 and singular special cases + @test UMFPACK.rcond(lu(SparseMatrixCSC{Tv,Ti}(sparse(Diagonal(Tv[3]))))) === 1.0 + @test UMFPACK.rcond(lu(SparseMatrixCSC{Tv,Ti}(sparse(Tv[1 2; 0 0])); check=false)) === 0.0 + # a factor without a numeric decomposition gets one on demand + G = UMFPACK.UmfpackLU(SparseMatrixCSC{Tv,Ti}(sparse(Tv[1 3; 0 1]))) + @test UMFPACK.rcond(G) === 0.25 + # lu! refreshes the estimate + lu!(F, SparseMatrixCSC{Tv,Ti}(sparse(Tv[1 1; 0 1]))) + @test UMFPACK.rcond(F) === 0.5 + end + @testset "deserialization" begin A = 10*I + sprandn(10, 10, 0.4) F1 = lu(A) From ea16779bb4ebc8dca17f5e2efda825322c6410c0 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Wed, 9 Sep 2026 05:28:52 -0400 Subject: [PATCH 3/3] Loop the UMFPACK rcond tests over the index types UMFPACK supports On 32-bit platforms UMFPACK only supports Int32 indices, so the testset now iterates Base.uniontypes(UMFPACK.UMFITypes) like the rest of the file. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_011BA2YkbzRZcTbTrc3S8odn --- test/umfpack.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/umfpack.jl b/test/umfpack.jl index b79aede6..b1c40b78 100644 --- a/test/umfpack.jl +++ b/test/umfpack.jl @@ -369,7 +369,7 @@ end end end - @testset "rcond (#118) for $Tv, $Ti" for Tv in (Float64, ComplexF64), Ti in (Int32, Int64) + @testset "rcond (#118) for $Tv, $Ti" for Tv in (Float64, ComplexF64), Ti in Base.uniontypes(UMFPACK.UMFITypes) # the number is min/max of |diag(U)| of the row-scaled matrix UMFPACK factorized F = lu(SparseMatrixCSC{Tv,Ti}(sparse(Tv[1 3; 0 1]))) @test UMFPACK.rcond(F) === 0.25