How to check if LU factorization failed?

Specific Domains Numerics
,
21

This also does not seem like the best solution, even if you could use issuccess()

22

Thank you @fredrikekre, what oneliner should I use for now?

Yes, I am not sure what I could do instead, I need a check inside of a loop, at least issuccess() would be more performant?

23

For LU-factorization: issuccess(F::LU) = F.info == 0

I don’t know your use case so its hard to tell. Why are the matrices possibly singular in the first place? What do you do if the factorization fails and what do you do if it succeeds?

24

@fredrikekre, I think I misinterpreted the interface then. I thought that LU.info returns 0 when the LU factorization succeeds and that LinAlg.issuccess stores the status after trying to use the already factorized object to solve the system LU \ b?

My use case is still under development, but you can see the try/catch block here: https://github.com/juliohm/GeoStats.jl/blob/master/src/solvers/sgsim_solver.jl#L143-L160

Please ignore that line of code that says “fix possible numerical issues”, I will remove it soon I fix the estimator interface. Basically, I need to add some preconditioning code in the estimator before I can do anything else.

25

yes

What do you mean? LinAlg.issuccess is just a function with the exact definition from my last post.

26

I thought that LinAlg.SingularException was replaced by LinAlg.issuccess, and that any form of exception handling was removed from LinAlg in Julia v0.7. So basically, I if I want to check for failure, I need to use exceptions, right?

27

Perhaps an example:

julia> A = zeros(2,2); b = rand(2);

julia> F = lufact(A); # this doesn't fail even though the matrix is singular

julia> LinAlg.issuccess(F) # false, since the factorization F is failed
false

julia> F\b # in the backsolve we throw if F is failed
ERROR: Base.LinAlg.SingularException(1)
[...]

So we still throw an SingularException, but issuccess can be used to know if we will throw before doing the backsolve.

1 Like
28

What about the success check for Cholesky? I don’t see a field info for it so that I could write my own issuccess in Julia <= v0.6.

29

That field was implemented in https://github.com/JuliaLang/julia/pull/21976 so it does not exist on 0.6. You can probably do some local backporting of that PR if you really need it for 0.6.

But if you sometimes will use matrices that are not positive definite, why not just use LU? Are you switching to LU in case of a PosDefException?

Edit: On 0.6 we actually have issuccess(F::Cholesky) = true since if you manage to obtain a Cholesky object it is always a success.

30

Ok, so I can just write issuccess(F::Cholesky) = true for now on my package and wait for Julia v0.7 where I can remove these definitions and use the ones in LinAlg.