# LU Factorization

Hi. I am playing around with Julia while I revisit my old Linear Algebra textbook. One thing I found is Julia’s LU seems to give different results than python’s scipy for some matrix: e.g.:

As a newbie, it could be I am using a wrong function or function with wrong parameters or wrong package. Any thoughts about it? Thank you very much for your time.

1 Like

Please add compilable code (that people can copy-paste and test) instead of screenshots. I’d add that the error appears in Julia version 0.7 onward, older versions give the same result as Python and MATLAB (I tested).

3 Likes

executable code rather than a screenshot is definitely helpful.

I think you want “lu(A,check=false)”

I assume that you want to do LU factorization of matrix `A`:

``````julia> A = [1 3 3 2;2 6 9 7; -1 -1 3 4]
3×4 Array{Int64,2}:
1   3  3  2
2   6  9  7
-1  -1  3  4

julia> using LinearAlgebra

julia> L,U,p = lu(A)
LU{Float64,Array{Float64,2}}
L factor:
3×3 Array{Float64,2}:
1.0  0.0  0.0
-0.5  1.0  0.0
0.5  0.0  1.0
U factor:
3×4 Array{Float64,2}:
2.0  6.0   9.0   7.0
0.0  2.0   7.5   7.5
0.0  0.0  -1.5  -1.5
``````

Here, `p` is the row permutation of `A` such that `L*U = A[p,:]`:

``````julia> L*U
3×4 Array{Float64,2}:
2.0   6.0  9.0  7.0
-1.0  -1.0  3.0  4.0
1.0   3.0  3.0  2.0

julia> A[p,:]
3×4 Array{Int64,2}:
2   6  9  7
-1  -1  3  4
1   3  3  2
``````

In other words: if you seek to solve `A*x = b`, it follows that `L*U*x = b[p]`, or `U*x = L\b[p]`.
Example: `b = [1,2,3]`:

``````julia> b = [1,2,3]
3-element Array{Int64,1}:
1
2
3

julia> A\b
4-element Array{Float64,1}:
-3.3333333333333526
2.0000000000000067
-1.666666666666658
1.6666666666666567
``````

Alternatively:

``````julia> U\(L\b[p])
4-element Array{Float64,1}:
-3.333333333333332
1.9999999999999987
-1.6666666666666656
1.6666666666666674
``````
2 Likes

This should be `A = [1 3 3 2;2 6 9 7; -1 -3 3 4]` to match OP (note the third-to-last entry). With that matrix the `SingularException` is reproducible. The decomposition returned by scipy is valid (`P * L * U` is indeed `A`, `L` is lower-triangular, and `U` is upper-triangular). Note that scipy’s `U` has a zero on the diagonal.

With `check=false`, Julia returns an `LU` object that prints as

``````Failed factorization of type LU{Float64,Array{Float64,2}}
``````

but the same `L`, `U` and `P` as the ones returned by scipy’s `lu` can be obtained from the returned `LU` object.

Julia’s `lu` calls `lu!`, which calls `LAPACK.getrf!`. Note that

``````julia> A, ipiv, info = LAPACK.getrf!(Float64.(A))
([2.0 6.0 9.0 7.0; 0.5 0.0 -1.5 -1.5; -0.5 0.0 7.5 7.5], [2, 2, 3], 2)
``````

where, according to Linear Algebra · The Julia Language `info == 2` indicates that `U[2,2]` is singular (and `info < 0` indicates failure). Julia considers the LU decomposition to be successful iff `info == 0`:

Perhaps it shouldn’t?

4 Likes

Oops… problem with reading off the screenshot… do I need new glasses?

Anyway, as you indicate, the `lu` factorization still works when setting `check = false`:

``````julia> using LinearAlgebra

julia> A = [1 3 3 2;2 6 9 7;-1 -3 3 4]
3×4 Array{Int64,2}:
1   3  3  2
2   6  9  7
-1  -3  3  4

julia> rank(A)
2

julia> L,U,p = lu(A,check=false)
Failed factorization of type LU{Float64,Array{Float64,2}}

julia> U
3×4 Array{Float64,2}:
2.0  6.0   9.0   7.0
0.0  0.0  -1.5  -1.5
0.0  0.0   7.5   7.5

julia> b = [1,2,3]
3-element Array{Int64,1}:
1
2
3

julia> A\b
4-element Array{Float64,1}:
-0.0952380952380952
-0.28571428571428564
0.2190476190476191
0.3142857142857143

julia> U\(L\b[p])
4-element Array{Float64,1}:
-0.10012210012210006
-0.30036630036630013
0.20634920634920634
0.30647130647130644

``````

So matrix `A` has rank loss, as indicated by matrix `U`. Finding the unknown `x` from `A*x = b` gives slightly different result depending on how `x` is computed, but that is due to the pseudo inverse algorithm, I guess.

Yes it seems a bit misleading that this is written as a “Failed factorization”, though the `SingularException` is arguably a reasonable default (if we think people will typically use the factorization for later solves, rather than directly using the factors for something else?). The lapack documentation says

if INFO = i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular.

The function `issuccess` comes from

though it looks like there wasn’t any discussion about what to do with positive `info` there and looking in `git` history doesn’t tell whether this was carefully chosen on purpose. @andreasnoack is this just an oversight?

1 Like

We have discussed the formulation somewhere. Right now I don’t remember in which issue or if it was only on Slack. Indeed, it’s misleading to call it failed for LU so it would better to change the formulation for LU. The changes that lead to this were motivated by the Cholesky which actually is failed when a positive error message is returned.

2 Likes

Maybe you are thinking of this discussion?

2 Likes

Exactly. Thanks.