# How to check if LU factorization failed?

**URL:** <https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553>\
**Category:** Numerics\
**Tags:** question, linearalgebra\
**Created:** [August 25, 2017, 3:50am UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553 "2017-08-25T03:50:26Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 25, 2017, 3:50am UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/1 "2017-08-25T03:50:26Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 25, 2017, 4:04am UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/2 "2017-08-25T04:04:43Z")

</div>

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`?

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [August 25, 2017, 4:16am UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/3 "2017-08-25T04:16:28Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 25, 2017, 4:29am UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/4 "2017-08-25T04:29:21Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [August 25, 2017, 8:26am UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/5 "2017-08-25T08:26:15Z")

</div>

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
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](https://github.com/JuliaLang/julia/pull/23315)).

`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](https://github.com/JuliaLang/julia/pull/22561) 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.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [August 25, 2017, 8:46am UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/6 "2017-08-25T08:46:56Z")

</div>

> 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](https://www.math.utah.edu/software/lapack/lapack-d/dsytrs.html)?

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [August 25, 2017, 9:16am UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/7 "2017-08-25T09:16:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [August 25, 2017, 10:07am UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/8 "2017-08-25T10:07:04Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [August 25, 2017, 1:42pm UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/9 "2017-08-25T13:42:06Z")

</div>

The basic idea is this:

```julia
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?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 25, 2017, 1:48pm UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/10 "2017-08-25T13:48:49Z")

</div>

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

```julia
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?

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [August 25, 2017, 2:04pm UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/11 "2017-08-25T14:04:02Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 25, 2017, 2:11pm UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/12 "2017-08-25T14:11:59Z")

</div>

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:

```julia
"""
    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:

```julia
"""
    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.

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [August 25, 2017, 3:03pm UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/13 "2017-08-25T15:03:53Z")

</div>

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](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](https://github.com/JuliaLang/julia/pull/22478)

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 25, 2017, 3:48pm UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/14 "2017-08-25T15:48:10Z")

</div>

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](https://discourse.julialang.org/t/symmetric-matrices/4086/9):
>
> 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](http://www.netlib.org/lapack/lapack-3.2.html), about haflway down the page). If these capabilities fit the bill, the price to be paid is in developer time.

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [August 25, 2017, 3:57pm UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/15 "2017-08-25T15:57:30Z")

</div>

> [@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?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 25, 2017, 4:04pm UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/16 "2017-08-25T16:04:34Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [August 25, 2017, 4:09pm UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/17 "2017-08-25T16:09:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [August 26, 2017, 3:26am UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/18 "2017-08-26T03:26:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [September 20, 2017, 5:07pm UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/19 "2017-09-20T17:07:17Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [September 20, 2017, 6:03pm UTC](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553/20 "2017-09-20T18:03:06Z")

</div>

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.

[Next page](https://discourse.julialang.org/t/how-to-check-if-lu-factorization-failed/5553.md?page=2)