check kwarg for factorizations #27336

fredrikekre · 2018-05-31T00:11:31Z

Still some cleanup to do.

andreasnoack

Looks good to me. Only thing I'm wondering is if throw=true is better than check=false.

fredrikekre · 2018-05-31T10:54:13Z

Stefan pointed out (and I agree) that throw = true looks weird if it does not actually throw (i.e. if the matrix is PD).

Also, should we do the same for CHOLMOD.cholesky[!]? I think that is a trivial extension and would be nice with a consistent interface, but then we should probably do the same for UMFPACK.lu[!] but that seemed a little bit more involved; we basically check issuccess there when we use the factorization. Is it possible to know if the factorization was successful before calling downstream functions in UMFPACK?

andreasnoack · 2018-05-31T11:08:42Z

I don't have strong opinions about the name so let's just stick with check.

It would be great to handle this consistently for the sparse versions as well. For LU, I guess it is a matter of storing the status value in the struct ~~instead of throwing right away~~ in

julia/stdlib/SuiteSparse/src/umfpack.jl

Lines 230 to 237 in 22590d5

    
           status = ccall(($num_r, :libumfpack), $itype, 
        
                          (Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Cvoid}, Ptr{Cvoid}, 
        
                           Ptr{Float64}, Ptr{Float64}), 
        
                          U.colptr, U.rowval, U.nzval, U.symbolic, tmp, 
        
                          umf_ctrl, umf_info) 
        
           if status != UMFPACK_WARNING_singular_matrix 
        
               umferror(status) 
        
           end

.

mbauman · 2018-05-31T21:36:39Z

So to use UMFPACK, we need to first precompute a symbolic factorization and then do the numeric thing. Is it possible for the umfpack_*_numeric to fail when the symbolic one didn't? The docs aren't quite clear on this front — it can return a singular matrix warning, but I'm not sure if that will happen in any cases where it didn't already do so in the symbolic pre-computation.

andreasnoack · 2018-05-31T22:14:42Z

I think the singular warning is only returned in the numerical factorization. The symbolic factorization doesn't fail on structural singularity. Instead of storing the status value, as I suggested in the last comment, we might just make the status != UMFPACK_WARNING_singular_matrix check conditional on !check.

fredrikekre · 2018-06-01T06:28:23Z

we might just make the status != UMFPACK_WARNING_singular_matrix check conditional on !check.

That means we have to propagate the check kwarg all the way to the function where the ccall happens though, probably better to store status in the struct and do the check "higher up"?

andreasnoack · 2018-06-01T06:51:43Z

That means we have to propagate the check kwarg all the way to the function where the ccall happens though

I think you'd only have to add a (positional?) argument to umfpack_numeric! and then, of course, the keyword argument to lu so it's a small change. I think it is simpler than extending the struct.

Sacha0 · 2018-06-02T16:48:08Z

(Apologies for review latency; I hope to eke out review time this weekend! :) )

Sacha0 · 2018-06-03T01:57:30Z

stdlib/LinearAlgebra/src/bunchkaufman.jl

 rook pivoting is not used.

+When `check = true` the matrix is checked for singularity, and makes sure
+the factorization is successful, or else throw an error. When `check = false`


Perhaps slightly more compact "When check = true, an error is thrown if the decomposition fails. When check = false, responsibility for checking the decomposition's validity (via issuccess) lies with the user." ?

Sacha0 · 2018-06-03T02:04:13Z

stdlib/LinearAlgebra/src/cholesky.jl

    end
+    check && checkpositivedefinite(info)
+    return C
 end


Perhaps slightly simpler?

function cholesky!(A::RealHermSymComplexHerm, ::Val{false}=Val(false); check::Bool = true) C, info = _chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular) check && checkpositivedefinite(info) return Cholesky(C.data, A.uplo, info) end

Sacha0 · 2018-06-03T02:08:11Z

stdlib/LinearAlgebra/src/cholesky.jl

+When `check = true` the matrix will be checked for positive definiteness, and make sure
+the factorization was successful, or else throw an error. When `check = false`
+the user is responsible to make sure the result is valid, with
+[`issuccess`](@ref), before passing the result to downstream functions.


Ref. comment on similar snippet above :).

Sacha0 · 2018-06-03T02:09:44Z

stdlib/LinearAlgebra/src/cholesky.jl

 inv(C::Cholesky{<:BlasFloat,<:StridedMatrix}) = inv!(copy(C))

 function inv(C::CholeskyPivoted)
-    chkfullrank(C)


I feel like we can afford the vowel in check 😄. Perhaps some day...

Sacha0 · 2018-06-03T02:10:44Z

stdlib/LinearAlgebra/src/factorization.jl

-   :($(esc(info)) == 0 ? $(esc(A)) : throw(SingularException($(esc(info)))))
-end
+checkpositivedefinite(info) = info == 0 || throw(PosDefException(info))
+checknonsingular(info) = info == 0 || throw(SingularException(info))


Sacha0 · 2018-06-03T02:12:12Z

stdlib/LinearAlgebra/src/generic.jl

        return convert(S, det(UpperTriangular(A)))
    end
-    return det(lu(A))
+    return det(lu(A; check = false))


Does this check = false belong here?

Yes, otherwise we throw, and det deals with failed factorizations:

julia/stdlib/LinearAlgebra/src/lu.jl

Line 333 in c8ce43a

issuccess(F) || return zero(T)

Sacha0 · 2018-06-03T02:15:30Z

stdlib/LinearAlgebra/src/lu.jl

    AA = similar(A, S)
    copyto!(AA, A)
-    F = lu!(AA, Val(false))
+    F = lu!(AA, Val(false); check = false) # ?


Re. the # ?, looks sound in combination with the downstream logic.

Sacha0 · 2018-06-03T02:18:50Z

stdlib/LinearAlgebra/test/bunchkaufman.jl

+    end
+    F = bunchkaufman(R; check = false)
+    @test sprint(show, "text/plain", F) == "Failed factorization of type $(typeof(F))"
+end


Do these tests preserve all that can/should be preserved from the block removed above? If so, carry on :). (Edit: E.g. eltype combinations.)

Yes, should hit all relevant code paths.

Sacha0 · 2018-06-03T02:24:23Z

stdlib/SuiteSparse/src/cholmod.jl

        Hermitian{T,SparseMatrixCSC{T,SuiteSparse_long}}};
-    shift = 0.0) where {T<:Real} =
-    cholesky!(F, Sparse(A); shift = shift)
+    shift = 0.0, check::Bool = true) where {T<:Real} =


Though the indentation was off before, perhaps indent this part of the signature to match the preceding part?

Sacha0 · 2018-06-03T02:31:59Z

stdlib/SuiteSparse/src/umfpack.jl

+function LinearAlgebra.issuccess(lu::UmfpackLU)
+    U, status = umfpack_numeric!(lu, false)
+    return status == UMFPACK_OK
+end


I haven't worked through the umfpack_numeric! logic, but presume that calling this issuccess method on a completed UmfpackLU will not result in an additional numerical decomposition phase?

I was wondering the same when I wrote this, my tests suggested it was an quick return, but I tried again now (with larger matrices) and that does not seem to be the case... Might have to store status in the struct after all!

Browsed through the source and indeed it doesn't look like it is returning early if the factorization is already computed so storing the status is probably the way to go. UMFPACK stores an estimate of the reciprocal condition number and this number is used to determine singularity but I don't think we have a way to extract that estimate.

Done; I think the diff is smaller for this implementation as well.

Sacha0

Looks great! Wonderful to cleanup to the use-sites :). Much thanks @fredrikekre!

Caveat: I did not examine the logic for the sparse decompositions carefully.

StefanKarpinski · 2018-06-07T19:45:16Z

FreeBSD failure is OOM; Windows failures are timeouts. Do we have reason to be suspicious of this PR on non-Linux platforms? It does pass on macOS at least, so it's fine on at least two different OSes.

fredrikekre · 2018-06-07T19:53:48Z

It passed on FreeBSD before the rebase as well.

StefanKarpinski · 2018-06-07T20:08:01Z

Ok, good to go then?

* check = (true|false) for LinearAlgebra.cholesky[!] * check = (true|false) for LinearAlgebra.lu[!] * check = (true|false) for LinearAlgebra.bunchkaufman[!] * check = (true|false) for CHOLMOD.cholesky[!] and CHOLMOD.ldlt[!] * check = (true|false) for UMFPACK.lu

fredrikekre added the linear algebra Linear algebra label May 31, 2018

fredrikekre requested review from Sacha0, andreasnoack and mbauman May 31, 2018 00:11

fredrikekre force-pushed the fe/check-or-not-to-check-that-is-the-question branch from db3ed25 to b8aac1f Compare May 31, 2018 09:52

andreasnoack reviewed May 31, 2018

View reviewed changes

andreasnoack closed this May 31, 2018

andreasnoack reopened this May 31, 2018

fredrikekre force-pushed the fe/check-or-not-to-check-that-is-the-question branch 2 times, most recently from 1e57e83 to e970dd6 Compare June 1, 2018 19:24

Sacha0 reviewed Jun 3, 2018

View reviewed changes

Sacha0 approved these changes Jun 3, 2018

View reviewed changes

fredrikekre force-pushed the fe/check-or-not-to-check-that-is-the-question branch 2 times, most recently from 30b7db2 to 112731f Compare June 4, 2018 06:37

fredrikekre requested review from Sacha0 and andreasnoack June 4, 2018 11:17

fredrikekre force-pushed the fe/check-or-not-to-check-that-is-the-question branch from 112731f to 6238946 Compare June 5, 2018 12:20

fredrikekre added 5 commits June 7, 2018 08:14

check = (true|false) for LinearAlgebra.cholesky[!]

6f95317

check = (true|false) for LinearAlgebra.lu[!]

b35d22a

check = (true|false) for LinearAlgebra.bunchkaufman[!]

8e7020c

check = (true|false) for CHOLMOD.cholesky[!] and CHOLMOD.ldlt[!]

937423c

check = (true|false) for UMFPACK.lu

fa6f586

fredrikekre force-pushed the fe/check-or-not-to-check-that-is-the-question branch from 6238946 to fa6f586 Compare June 7, 2018 06:16

fredrikekre closed this Jun 7, 2018

fredrikekre reopened this Jun 7, 2018

andreasnoack merged commit 37d70a5 into master Jun 7, 2018

andreasnoack deleted the fe/check-or-not-to-check-that-is-the-question branch June 7, 2018 20:20

haampie mentioned this pull request Jun 19, 2018

[LinearAlgebra] LU-factorization failure (0.7.0-alpha.0) JuliaLang/LinearAlgebra.jl#535

Closed

tpapp mentioned this pull request Apr 21, 2020

Add singularity checks for lu, method for issuccess. JuliaArrays/StaticArrays.jl#766

Merged

Uh oh!

check kwarg for factorizations #27336

check kwarg for factorizations #27336

Uh oh!

Conversation

fredrikekre commented May 31, 2018

Uh oh!

andreasnoack left a comment

Choose a reason for hiding this comment

Uh oh!

fredrikekre commented May 31, 2018

Uh oh!

andreasnoack commented May 31, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mbauman commented May 31, 2018

Uh oh!

andreasnoack commented May 31, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

fredrikekre commented Jun 1, 2018

Uh oh!

andreasnoack commented Jun 1, 2018

Uh oh!

Sacha0 commented Jun 2, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Sacha0 Jun 3, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

andreasnoack Jun 3, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Sacha0 left a comment

Choose a reason for hiding this comment

Uh oh!

StefanKarpinski commented Jun 7, 2018

Uh oh!

fredrikekre commented Jun 7, 2018

Uh oh!

StefanKarpinski commented Jun 7, 2018

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

6 participants

andreasnoack commented May 31, 2018 •

edited

Loading

andreasnoack commented May 31, 2018 •

edited

Loading

Sacha0 Jun 3, 2018 •

edited

Loading

andreasnoack Jun 3, 2018 •

edited

Loading