The Python – Scipy linalg LU decomposition provides different results for my textbook

Scipy linalg LU decomposition provides different results for my textbook… here is a solution to the problem.

Scipy linalg LU decomposition provides different results for my textbook

I

am studying the 5th edition of linear programming by Saul I. Gass.

He gave the following text and examples:
“Given an n x n non-singular matrix A, then A can be expressed as the product A = LU If we make the diagonal elements of U (or L) all equal to 1, the LU decomposition will be unique

A=LU

I

managed to get a reversible up-down decomposition with this code I found from this SO problem: Is there a built-in/easy LDU decomposition method in Numpy?

But I still don’t know what’s going on and why L and U are so different from my textbooks. Can anyone explain it to me?

So this code:

import numpy as np
import scipy.linalg as la
a = np.array([[1, 1, -1],
              [-2, 1, 1],
              [1, 1, 1]])
(P, L, U) = la.lu(a)

print(P)
print(L)
print(U)

D = np.diag(np.diag(U))   # D is just the diagonal of U
U /= np.diag(U)[:, None]  # Normalize rows of U
print(P.dot(L.dot(D.dot(U))))    # Check

Give this output:

[[ 0.  1.  0.]
 [ 1.  0.  0.]
 [ 0.  0.  1.]]
[[ 1.   0.   0. ]
 [-0.5  1.   0. ]
 [-0.5  1.   1. ]]
[[-2.   1.   1. ]
 [ 0.   1.5 -0.5]
 [ 0.   0.   2. ]]
[[ 1.  1. -1.]
 [-2.  1.  1.]
 [ 1.  1.  1.]]

Solution

You can choose which matrix (L or U) should have one diagonally. The textbook example chooses U, but the implementation of scipy chooses L. This explains the difference.

To illustrate this, we can turn things around:

(P, L, U) = la.lu(a.T)

print(P.T)
# [[ 1.  0.  0.]
#  [ 0.  1.  0.]
#  [ 0.  0.  1.]]
print(L.T)
# [[ 1.          1.         -1.        ]
#  [ 0.          1.         -0.33333333]
#  [ 0.          0.          1.        ]]
print(U.T)
# [[ 1.  0.  0.]
#  [-2.  3.  0.]
#  [ 1.  0.  2.]]

By transposing the matrix, we basically swap you and L so that the other matrix gets 1 diagonally. And, lo and behold, the result is the same as in the textbooks.

(Note that if the permutation matrix P is not an identity matrix, the result will look a little different.) )

Related Problems and Solutions