# LU Factorization 

**URL:** <https://discourse.julialang.org/t/lu-factorization/20763>\
**Category:** New to Julia\
**Created:** [February 13, 2019, 8:02pm UTC](https://discourse.julialang.org/t/lu-factorization/20763 "2019-02-13T20:02:32Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![ZQ\_Li](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zq_li/32/4954_2.png) [@ZQ\_Li](https://discourse.julialang.org/u/ZQ_Li)\
**Post date:** [February 13, 2019, 8:02pm UTC](https://discourse.julialang.org/t/lu-factorization/20763/1 "2019-02-13T20:02:32Z")

</div>

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

 ![combined](https://global.discourse-cdn.com/julialang/original/3X/a/a/aade2e152fa931fe7ccf1ccf5dbb5a8db6fb48bb.jpeg)

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.

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [February 13, 2019, 8:31pm UTC](https://discourse.julialang.org/t/lu-factorization/20763/2 "2019-02-13T20:31:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![chrisvwx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisvwx/32/45289_2.png) [@chrisvwx](https://discourse.julialang.org/u/chrisvwx)\
**Post date:** [February 13, 2019, 8:33pm UTC](https://discourse.julialang.org/t/lu-factorization/20763/3 "2019-02-13T20:33:52Z")

</div>

executable code rather than a screenshot is definitely helpful.

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

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [February 13, 2019, 10:37pm UTC](https://discourse.julialang.org/t/lu-factorization/20763/4 "2019-02-13T22:37:19Z")

</div>

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

```julia
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
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
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
julia> U\(L\b[p])
4-element Array{Float64,1}:
 -3.333333333333332
  1.9999999999999987
 -1.6666666666666656
  1.6666666666666674

```

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [February 14, 2019, 2:41am UTC](https://discourse.julialang.org/t/lu-factorization/20763/5 "2019-02-14T02:41:29Z")

</div>

> [@BLI](#):
>
> ```julia
> A = [1 3 3 2; 2 6 9 7; -1 -1 3 4]
> 
> ```

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

```julia
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
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](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.LAPACK.getrf!), `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`:

> <https://github.com/JuliaLang/julia/blob/38d624785216bcb1374f651831436dc886a6bb4b/stdlib/LinearAlgebra/src/lu.jl#L304>

Perhaps it shouldn’t?

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [February 14, 2019, 7:22am UTC](https://discourse.julialang.org/t/lu-factorization/20763/6 "2019-02-14T07:22:43Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![c42f](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/c42f/32/52842_2.png) [@c42f](https://discourse.julialang.org/u/c42f)\
**Post date:** [February 14, 2019, 7:34am UTC](https://discourse.julialang.org/t/lu-factorization/20763/7 "2019-02-14T07:34:36Z")

</div>

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](http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga0019443faea08275ca60a734d0593e60.html) 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

> <https://github.com/JuliaLang/julia/pull/22345>
>
> Introduce \`issuccess\` function to test if a \`Factorization\` succeeded.
> 
> Fixes …#22335

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?

---

<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:** [February 15, 2019, 2:38pm UTC](https://discourse.julialang.org/t/lu-factorization/20763/8 "2019-02-15T14:38:01Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [February 17, 2019, 7:14pm UTC](https://discourse.julialang.org/t/lu-factorization/20763/9 "2019-02-17T19:14:56Z")

</div>

Maybe you are thinking of this discussion?

[https://github.com/JuliaLang/julia/issues/27657](https://github.com/JuliaLang/julia/issues/27657)

---

<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:** [February 17, 2019, 7:39pm UTC](https://discourse.julialang.org/t/lu-factorization/20763/10 "2019-02-17T19:39:08Z")

</div>

Exactly. Thanks.