SciPy - Matrix LU Decomposition

---

---

What is LU Decomposition?

LU decomposition is also called as LU factorization. This is a mathematical technique used to decompose a matrix into the product of two simpler matrices namely a lower triangular matrix L and an upper triangular matrix U. This is particularly useful in solving systems of linear equations, inverting matrices and calculating determinants.

We can explain the LU decomposition as, when given a square matrix A then Lu decomposition finds two matrices as follows −

• L: A lower triangular matrix with ones on the diagonal l_ii = 1.
• U: An upper triangular matrix.

Such that −

A = L . U

If partial pivoting is used to improve numerical stability then the decomposition includes a permutation matrix P which results in −

P. A = L . U

Here P is a permutation matrix representing row swaps.

Steps in Lu Decomposition

Following are the steps in Lu Decomposition in scipy −

• Start with the matrix A.
• Decompose A into L and U by systematically eliminating elements below the diagonal i.e., Gaussian elimination.
• Optionally include row permutations to form P by ensuring the pivot elements i.e., diagonal entries of U which are nonzero and numerically stable.

Implementation of Matrix LU Decomposition in scipy

In Python's SciPy library LU decomposition is performed using scipy.linalg.lu which provides −

• P, L and U matrices directly.
• Functionality for solving equations using LU decomposition via scipy.linalg.lu_solve in combination with scipy.linalg.lu_factor.

Decomposing a Matrix into P,L,U

Following is the example of performing LU decomposition of a square matrix into permutation, lower and upper triangular matrices by using the function scipy.linalg.lu() −

import numpy as np
from scipy.linalg import lu

# Define a square matrix A
A = np.array([[4, 3], [6, 3]])

# Perform LU decomposition
P, L, U = lu(A)

# Display the results
print("Original Matrix A:")
print(A)
print("\nPermutation Matrix P:")
print(P)
print("\nLower Triangular Matrix L:")
print(L)
print("\nUpper Triangular Matrix U:")
print(U)

# Verify PA = LU
PA = P @ A
LU = L @ U
print("\nCheck if PA equals LU:")
print(np.allclose(PA, LU))

Following is the output of decomposing a matrix into P, L, U in scipy −

Original Matrix A:
[[4 3]
 [6 3]]

Permutation Matrix P:
[[0. 1.]
 [1. 0.]]

Lower Triangular Matrix L:
[[1.         0.        ]
 [0.66666667 1.        ]]

Upper Triangular Matrix U:
[[6. 3.]
 [0. 1.]]

Check if PA equals LU:
True

Solving a Linear System Using LU Decomposition

We can solve the linear system using the matrix LU Decomposition in scipy with the help of the functions scipy.linalg.lu() and scipy.linalg.lu_solve(). In this example we are solving AX = B for a given A and B −

import numpy as np
from scipy.linalg import lu_factor, lu_solve

# Define matrix A and right-hand side vector B
A = np.array([[3, 1], [1, 2]])
B = np.array([9, 8])

# Factorize A into LU components
lu, piv = lu_factor(A)

# Solve for X in AX = B
X = lu_solve((lu, piv), B)

# Display the results
print("Solution X:")
print(X)

# Verify the solution
print("\nCheck if A @ X equals B:")
print(np.allclose(A @ X, B))

Here is the output of solving the linear system AX = B −

Solution X:
[2. 3.]

Check if A @ X equals B:
True

LU Decomposition Without Explicit Permutation Matrix

Here is the example of finding LU Decomposition without Explicit Permutation Matrix and return PA directly without separate P −

import numpy as np
from scipy.linalg import lu

# Define a matrix A
A = np.array([[2, 3, 1], [4, 7, 1], [6, 18, 1]])

# Perform LU decomposition with permuted L
L_permuted, U = lu(A, permute_l=True)

# Display results
print("Matrix A:")
print(A)
print("\nPermuted Lower Triangular Matrix (L_permuted):")
print(L_permuted)
print("\nUpper Triangular Matrix U:")
print(U)

# Verify PA = LU
PA = L_permuted @ U
print("\nReconstructed Matrix (PA):")
print(PA)
print("\nCheck if PA equals A:")
print(np.allclose(PA, A))

Here is the output of the LU Decompositton without explicit Permutation matrix −

Matrix A:
[[ 2  3  1]
 [ 4  7  1]
 [ 6 18  1]]

Permuted Lower Triangular Matrix (L_permuted):
[[0.33333333 0.6        1.        ]
 [0.66666667 1.         0.        ]
 [1.         0.         0.        ]]

Upper Triangular Matrix U:
[[ 6.         18.          1.        ]
 [ 0.         -5.          0.33333333]
 [ 0.          0.          0.46666667]]

Reconstructed Matrix (PA):
[[ 2.  3.  1.]
 [ 4.  7.  1.]
 [ 6. 18.  1.]]

Check if PA equals A:
True