Expose the CHOLMOD and UMFPACK reciprocal condition number estimates as rcond

docs/src/solvers.md

@@ -29,8 +29,10 @@ SparseArrays.CHOLMOD.lowrankdowndate
 SparseArrays.CHOLMOD.lowrankdowndate!
 SparseArrays.CHOLMOD.lowrankupdowndate!
 SparseArrays.CHOLMOD.ldlt
+SparseArrays.CHOLMOD.rcond
 SparseArrays.SPQR.qr
 SparseArrays.UMFPACK.lu
+SparseArrays.UMFPACK.rcond
 ```
 
 ```@meta

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,41 @@ 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. 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.
+
+# 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)

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} =

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

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 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
+        @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)