How to check if LU factorization failed?

The basic idea is this:

F = lufact(A)
(F.info == 0) || error("factorization failed")
kappa = cond(F,1)
relative_error_estimate = kappa * max(expected_error(A),expected_error(b),eps())

This uses a condition estimate based on the factorization, which is much cheaper than the SVD-based precise condition you get from cond(A). Thanks to your question, I learned that such an estimate is not readily available for bkfact results. (The relevant routine from LAPACK doesn’t yet have a friendly wrapper.) That, and the poorer performance of the back-solver pointed out by others, are the disadvantages.

Thanks @fredrikekre for correcting my recollection about when the exception would be thrown. Is the error check in my snippet above reliable?

1 Like

Thank you @fredrikekre and @Ralph_Smith, very helpful information. Can we conclude that a good way of checking for failure is catching the exception?

try
  A \ b
catch e
  if e isa SingularException
    # remediate
  end
end

Is it ok to have a try/catch block when I am solving thousands of RHS? What is the overhead?

Related question, do I increase performance by passing lufact(Symmetric(A)) instead of lufact(A) when A is value symmetric?

Internally in the linear algebra code, we have gone away from catching exceptions. Instead, we check the factorizations with LinAlg.issuccess.

Related question, do I increase performance by passing lufact(Symmetric(A)) instead of lufact(A) when A is value symmetric?

lufact works on a full matrix so there is nothing to be gained from lufact(Symmetric(A)) compared to lufact(A).

1 Like

Thank you @andreasnoack, can you confirm LinAlg.issuccess will be available in Julia v0.7?

Regarding the Symmetric(A) question, I am planning to rewrite a utility function from:

"""
    pairwise(metric, X)

Evaluate symmetric `metric` between all n² pairs of columns in
a m-by-n matrix `X` efficiently.
"""
function pairwise(metric::Function, X::AbstractMatrix)
  m, n = size(X)
  D = zeros(n, n)
  for j=1:n
    xj = view(X, :, j)
    for i=j+1:n
      xi = view(X, :, i)
      @inbounds D[i,j] = metric(xi, xj)
    end
    @inbounds D[j,j] = metric(xj, xj)
    for i=1:j-1
      @inbounds D[i,j] = D[j,i] # leveraging the symmetry
    end

  end

  D
end

to:

"""
    pairwise(metric, X)

Evaluate symmetric `metric` between all n² pairs of columns in
a m-by-n matrix `X` efficiently.
"""
function pairwise(metric::Function, X::AbstractMatrix)
  m, n = size(X)
  D = zeros(n, n)
  for j=1:n
    xj = view(X, :, j)
    for i=j+1:n
      xi = view(X, :, i)
      @inbounds D[i,j] = metric(xi, xj)
    end
    @inbounds D[j,j] = metric(xj, xj)
  end

  Symmetric(D, :L)
end

assuming that all the solvers like lufact that expect a symmetric matrix won’t touch the upper part of the internal array. Is this a safe assumption? I can save some memory access in this case.

First thing that happens when calling lufact with a matrix wrapped in Symmetric is that the upper/lower part is copied to the other one. You don’t save memory, but you call lufact(Symmetric(A)) and let julia do the copying. See https://github.com/JuliaLang/julia/blob/90a33006980bf982f223d209671e9dd87641b2f2/base/linalg/lu.jl#L22-L25

Note that lufact(::Symmetric) does not work on julia v0.6. Implemented in https://github.com/JuliaLang/julia/pull/22478

Thank you @fredrikekre, very helpful information again. From what I understood, if I stick with the idea of replacing my pairwise implementation above to return a Symmetric(D), I will get into trouble in Julia v0.6.

I won’t change it for now, but it would be great to have these fixed and announced somewhere (perhaps here).

Regarding memory savings, I remember posting a question some time ago about the fact that Symmetric is using the full array instead of just half, @Ralph_Smith did mention that there is support in LAPACK, but no one has had the time to implement it yet:

Well, you can just emulate the 0.7 behaviour from that PR (you basically only need to copytri! and call lufact on the full matrix).

Not sure what you mean, fix what?

I saw that the PR was merged, I meant to maybe back port it to Julia v0.6?

Probably not. It is easy enough to do something like lufact(full(Symmetric(A))) which should behave equivalently and will work on both versions.

Edit: You can even do an inplace lufact since full will allocate a new matrix: lufact!(full(Symmetric(A)))

Just a small comment, if you are going with lufact, then actually no effort is made to preserve symmetry during pivoting (in fact quite the opposite) so your L and U will not be symmetric (not even with some minor re-scaling), and hence a full matrix is needed to represent the L and U. This means that no memory savings can happen at the output. Here is an example and a sanity check:

julia> a = rand(5,5)
5×5 Array{Float64,2}:
 0.16469   0.0148345  0.371097    0.579042  0.694243
 0.95235   0.746349   0.00304103  0.483745  0.888242
 0.41297   0.525047   0.947458    0.175485  0.906909
 0.538816  0.0940616  0.80986     0.464461  0.780789
 0.633868  0.558475   0.4369      0.856279  0.815061

julia> for i in 1:5
           for j in i+1:5
               a[i,j] = a[j,i]
           end
       end

julia> a
5×5 Array{Float64,2}:
 0.16469   0.95235    0.41297   0.538816   0.633868
 0.95235   0.746349   0.525047  0.0940616  0.558475
 0.41297   0.525047   0.947458  0.80986    0.4369
 0.538816  0.0940616  0.80986   0.464461   0.856279
 0.633868  0.558475   0.4369    0.856279   0.815061

julia> fact = lufact(a);

julia> fact[:L]
5×5 Array{Float64,2}:
 1.0        0.0       0.0        0.0        0.0
 0.17293    1.0       0.0        0.0        0.0
 0.565775  -0.398652  1.0        0.0        0.0
 0.665583   0.074965  0.0986935  1.0        0.0
 0.433632   0.244637  0.999576   0.0316445  1.0

julia> fact[:U]
5×5 Array{Float64,2}:
 0.95235  0.746349  0.525047  0.0940616   0.558475
 0.0      0.823284  0.322174  0.52255     0.537291
 0.0      0.0       0.641237  0.619559    0.7545
 0.0      0.0       0.0       0.693353    0.328607
 0.0      0.0       0.0       0.0        -0.701292

Symmetric triangular decompositions only use symmetric permutations to preserve symmetry.

So faced by the weird evidence that LU is faster than BK (a clear contradiction to the implementation-free theory), and assuming Symmetric(D) can be made to save only one side of the matrix as it is proposed in the other post, it seems like this is one of those running time vs memory decisions you have to make. Although in theory, BK should strictly dominate LU in both memory usage and running time, but clearly it is quite the opposite where we stand right now.

@andreasnoack is there any chance that LinAlg.issuccess will be available in Julia v0.6? If not, do you have any estimate of when Julia v0.7 is coming out? I really need to go away from exception handling, I have a try/catch block that is being executed thousands of times. I am guessing this is not good for performance.

Also, any chance that LL.info for Cholesky will be available in Julia v0.6? Currently I can only access this field for LU factorization. It is very hard to provide a consistent interface currently.

Those won’t be backported. LinAlg.issuccess could be added to Compat I guess, but you might as well define your own local function, it is only one line.

1 Like

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

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?

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?

@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.

yes

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

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?

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

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.