Why does scipy.linalg.lu() return the wrong decomposition matrices of square matrix B in this code?

Question

I have the following code of which the outcome is very confusing:

import numpy as np
import scipy.linalg
B = np.array([[2,-1,-2],[-4,6,3],[-4,-2,8]])
P,L,U = scipy.linalg.lu(B)
print(L)

Which returns the following:

[[ 1. 0. 0. ] [ 1. 1. 0. ] [-0.5 -0.25 1. ]]

But this is not the correct L matrix in the LU-decomposition of B. As far as I know the command scipy.linalg.lu(matrix) simply computes the LU decomposition matrices of whatever matrix you put in. However, in this case the L matrix is not correct. What is going on here? Any help is appreciated.

Why is this not the correct lower-matrix in the decomposition? What is the correct matrix in your opinion? — ShlomiF, Dec 06 '18 at 20:45
*"...the command scipy.linalg.lu(matrix) simply computes the LU decomposition..."* More precisely, it computes the factorization P@L@U, where P is a permutation matrix, L is lower trianguler and U is upper triangular. Did you take P into account when you assessed the correctness of the result? — Warren Weckesser, Dec 06 '18 at 22:33
@ShlomiF The correct matrix in my opinion, using an online LU decomposition tool like this one: https://matrixcalc.org/en/ is the following: ([[1,0,0],[-2,1,0],[-2,-1-,1]]). This is not the L-matrix my code returns, hence my confusion. — Charlie Shuffler, Dec 07 '18 at 07:53

ShlomiF · Answer 1 · 2018-12-08T21:22:14.810

I think you might be misunderstanding what the decomposition should be doing. 'cause it looks correct to me... Let's go through your example's results, with some extra details and comments:

import numpy as np
import scipy.linalg
B = np.array([[2,-1,-2],[-4,6,3],[-4,-2,8]])
P,L,U = scipy.linalg.lu(B)  

# Now let's see if P is a permutation matrix (a single 1 in each row and column, all the rest should be zero):
print(P)  
>> [[ 0.  0.  1.]
   [ 1.  0.  0.]
   [ 0.  1.  0.]]
# Yup! That's a successful permutation matrix  

# Now let's see if L is a lower triangular matrix (anything above the main diagonal should be zero):
print(L)
>> [[ 1.    0.    0.  ]
   [ 1.    1.    0.  ]
   [-0.5  -0.25  1.  ]]
# Yup! still doing good.  

# Now let's see if U is a upper triangular matrix (anything below the main diagonal should be zero):
print(U)
>> [[-4.    6.    3.  ]
   [ 0.   -8.    5.  ]
   [ 0.    0.    0.75]]

# And last but not least, let's see if this works as a decomposition of B (i.e. PLU==B):  
print(np.matmul(P, np.matmul(L, U)))
>> [[ 2. -1. -2.]
   [-4.  6.  3.]
   [-4. -2.  8.]]

# :-)

I hope that clarifies things a bit. If you're still not sure, then maybe reread about permutation matrices, triangular matrices, lu-decomposition, scipy.linalg.lu, and closely related topics.
Good luck!

It seems a clarification is in place:
LU-decomposition is not - in the general case - necessarily unique.
If you want details then aside for the relevant sub-chapter in the above wikipedia link, I recommend the first and third answers on this stack-exchange question.
So if you happen to get two different answers from to different implementations or methods, that doesn't mean that one or the other is wrong. If you've got a permutation matrix P (even if it's the trivial identity matrix), a lower-matrix L, an upper-matrix U, and they decompose your matrix, then you've got yourself a decomposition. Hope that clears things up!

But how come the L matrix it returns is not the L matrix you would see when using an online LU decomposition tool such as: https://matrixcalc.org/en/ ? The tool, and a computation by hand, return the matrix ([[1,0,0],[-2,1,0],[-2,-1,1]]). How do I obtain this L matrix using the matrices that scipy.linalg.lu returns? — Charlie Shuffler, Dec 07 '18 at 07:56
I'm not 100% sure about this but I do know that LU decomposition is not unique and that permutation matrices are used to make the algorithm run more efficiently so I am guessing the difference is that the site you used & hand-done exercise is the long way, but for efficiency the scipy function incorporates the permutation matrix and thereby extracts another possible LU decomposition. I would be happy to be corrected if I am mistaken here. — Alxmrphi, Dec 07 '18 at 15:14

Why does scipy.linalg.lu() return the wrong decomposition matrices of square matrix B in this code?

1 Answers1