# How to check if LU factorization failed?

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.

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

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.

2 Likes

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?

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.

2 Likes

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?

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.

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