SciPy - linalg.lu() Function

The scipy.linalg.lu() is a function in scipy which performs LU decomposition of a given matrix ( A ) by breaking it into three components namely, a permutation matrix ( P ), a lower triangular matrix ( L ) and an upper triangular matrix ( U ) such that ( PA = LU ).

This function is used to analyze and solve systems of linear equations, invert matrices and compute determinants. With the permute_l=True option it returns a combined permuted lower triangular matrix ( L_{text{permuted}} ) and ( U ) by skipping ( P ). It is a powerful tool for numerical stability and efficient computation in linear algebra.

Syntax

Following is the syntax of the function scipy.linalg.lu() to perform LU Decomposition −

scipy.linalg.lu(a, permute_l=False, overwrite_a=False, check_finite=True, p_indices=False)

Parameters

Below are the parameters of the scipy.linalg.lu() function −

• a (array-like, shape = (M, M)): This is the input square matrix to decompose. It should be a 2D array.
• permute_l (bool, default: False): If True then the function will return the matrix PA=LU where P is the permutation matrix. The lower triangular matrix L will be returned with the row permutations already applied. If False then the function returns P, L, U.
• overwrite_a (bool, default: False): If True then the function will modify the input matrix a in place. If False then a copy of a will be used to perform the decomposition by leaving a unchanged.
• check_finite (bool, default: True): If True then checks whether the input matrix contains only finite numbers i.e., no NaN or inf values and if False then disables this check. This can be useful for performance if we are certain the input is finite.
• p_indices (bool, default: False):: If True then the function will return the permutation indices instead of the permutation matrix P and if False then the function returns the permutation matrix P directly.

Return Value

The scipy.linalg.lu() function returns the following values depends upon the permute_l parameter −

• If permute_l=False (default) then returns P, L, U, if permute_l=True then it returns L, U. If p_indices=Trueand permute_l=False then P as an array of pivot indices instead of the permutation matrix. If p_indices=True andpermute_l=True then returns the permuted lower triangular matrix L and the upper triangular matrix U but the permutation is applied through pivot indices.

Standard LU Decomposition (permute_l=False)

Following is the example which uses scipy.linalg.lu() function to perform basic LU Decomposition of a square matrix into P, L and U−

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)

print("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))

Here is the output of the scipy.linalg.lu() function −

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

LU Decomposition Without Permutation Matrix

To perform LU decomposition without returning a separate permutation matrix we can use the permute_l=True parameter in scipy.linalg.lu() function. This function returns a combined permuted lower triangular matrix (PL) and the upper triangular matrix (U) −

import numpy as np
from scipy.linalg import lu

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

# Perform LU decomposition with permute_l=True
PL, U = lu(A, permute_l=True)

# Display the results
print("Original Matrix A:")
print(A)

print("\nPermuted Lower Triangular Matrix (PL):")
print(PL)

print("\nUpper Triangular Matrix (U):")
print(U)

# Verify PL @ U = A
reconstructed_A = PL @ U
print("\nReconstructed Matrix (PL @ U):")
print(reconstructed_A)

print("\nCheck if PL @ U equals A:")
print(np.allclose(reconstructed_A, A))  # Should return True

Here is the output of the scipy.linalg.lu() function which is used to perform LU decomposition without permutation matrix −

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

Permuted Lower Triangular Matrix (PL):
[[0.66666667 1.        ]
 [1.         0.        ]]

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

Reconstructed Matrix (PL @ U):
[[4. 3.]
 [6. 3.]]

Check if PL @ U equals A:
True

LU Decomposition With Pivot Indices

In LU decomposition with pivoting the matrix A is factorized into three matrices namely, a permutation matrix P, a lower triangular matrix L and an upper triangular matrix U. Additionally the pivot indices are used to track the row swaps applied during the decomposition to improve numerical stability. Below is the example which retrieves pivot indices instead of the permutation matrix −

import numpy as np
from scipy.linalg import lu

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

# Perform LU decomposition with pivot indices
piv_indices, L, U = lu(A, p_indices=True)

print("Matrix A:")
print(A)
print("\nPivot Indices:")
print(piv_indices)
print("\nLower Triangular Matrix L:")
print(L)
print("\nUpper Triangular Matrix U:")
print(U)

Here is the output of the scipy.linalg.lu() function used to perform Lu Decomposition with pivot indices −

Matrix A:
[[3 1]
 [6 2]]

Pivot Indices:
[1 0]

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

Upper Triangular Matrix U:
[[6. 2.]
 [0. 0.]]