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.
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).
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
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?
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?
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.
Maybe you are thinking of this discussion?
Exactly. Thanks.