How to check if LU factorization failed?

[juliohm]

I am currently using lufact() in a project, but I couldn’t figure out by just reading the docs how one can check if the factorization failed. Normally, one would check if the pivots are zero, but how to do that with the LU object returned by lufact()?

Also, I was reading the Linear Algebra docs and encountered the bkfact() factorization. It is apparently the second option after Cholesky that the operator \ picks up when the matrix is not positive definite, but only symmetric. Anyone familiar with this factorization could elaborate on the advantages of using it instead of LU for example? First time reading about it.

[juliohm]

I think for the first question, I can assume that all square matrices have a LU factorization with partial pivoting (and that is the default algorithm in lufact()? My matrices are square matrices, what we can conclude in terms of the system Ax=b?

[Ralph_Smith]

If you work with BLAS-compatible types (Float64, Complex128, etc), lufact() should throw an LAPACK exception in case of failure. For other types, you may have to poke around in the specific method code to make sure the checks are done (but the LinAlg developers are usually careful about this sort of thing).

The Bunch-Kauffman scheme requires about half as many operations as LU with pivoting for a given matrix.

The experts say that problems with partial pivoting are “extremely rare”. You can always check your residual if you are nervous, and it’s a good idea to get a condition number for A for a handle on achievable accuracy.

[juliohm]

Thank you @Ralph_Smith, very helpful. Would you mind sharing some code to illustrate the workflow with the condition number and the factorization types in Julia? Also, what is the disadvantage of Bunch-Kauffman compared to LU?

[fredrikekre]

You will indeed get an error when using the LU object. Note though, that the error is not thrown when doing the factorization, but rather when the factorization is used for backsolve for instance.

julia> F = lufact(zeros(3,3)); # no error here

julia> F\rand(3); # this will throw
ERROR: Base.LinAlg.SingularException(1)
Stacktrace:
 [1] A_ldiv_B! at ./linalg/lu.jl:238 [inlined]
 [2] \(::Base.LinAlg.LU{Float64,Array{Float64,2}}, ::Array{Float64,1}) at ./linalg/factorization.jl:48

This is also the case for bkfact (and on 0.7 and above this will also be the case for Cholesky factorization, see this PR).

bkfact will only work if your matrix is symmetric, in which case you can do F = bkfact(Symmetric(A)) which will avoid checking if the matrix is indeed symmetric, this will ignore either the upper/lower part of the data. Regarding performance, see some discussion here where it was essentially concluded that if you have multiple right hand sides it is faster to use lufact (even for the symmetric case). But if you solve for one right hand side bkfact seems to be faster.

[antoine-levitt]

Regarding performance, see some discussion here1 where it was essentially concluded that if you have multiple right hand sides it is faster to use lufact (even for the symmetric case). But if you solve for one right hand side bkfact seems to be faster.

The reason given there was that BK is not able to use BLAS3 for multiple RHS solves. Any reason I’m missing why that does not just call lapack-d/dsytrs.html?

[fredrikekre]

We call that already, see https://github.com/JuliaLang/julia/blob/ceeb591003642f5c6e708df2ee4ef35d561394e2/base/linalg/bunchkaufman.jl#L243

[antoine-levitt]

Interesting. Running the benchmark in that thread, I get consistently worse results for Bunch-Kaufman than for LU, whether it’s for one or many RHS. This is probably strongly LAPACK/BLAS dependent.

[Ralph_Smith]

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?

[juliohm]

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?

[andreasnoack]

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

[juliohm]

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.

[fredrikekre]

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

[juliohm]

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:

Symmetric matrices

For some operations (esp. Cholesky and triangular solves) LAPACK has routines for rectangular full packed format (newer than the format @nalimilan and @stevengj address above), which enables the use of half the full storage while maintaining efficiency by using Level 3 BLAS/LAPACK kernels. (detaiils here, about haflway down the page). If these capabilities fit the bill, the price to be paid is in developer time.

[fredrikekre]

@juliohm

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.

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

@juliohm

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

Not sure what you mean, fix what?

[juliohm]

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

[fredrikekre]

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

[mohamed82008]

@juliohm

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

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.

[juliohm]

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

[fredrikekre]

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.