mirror of
https://github.com/TheAlgorithms/Python.git
synced 2024-12-21 10:40:14 +00:00
4700297b3e
* Enable ruff RUF002 rule * Fix --------- Co-authored-by: Christian Clauss <cclauss@me.com>
113 lines
3.8 KiB
Python
113 lines
3.8 KiB
Python
"""
|
|
Lower-upper (LU) decomposition factors a matrix as a product of a lower
|
|
triangular matrix and an upper triangular matrix. A square matrix has an LU
|
|
decomposition under the following conditions:
|
|
- If the matrix is invertible, then it has an LU decomposition if and only
|
|
if all of its leading principal minors are non-zero (see
|
|
https://en.wikipedia.org/wiki/Minor_(linear_algebra) for an explanation of
|
|
leading principal minors of a matrix).
|
|
- If the matrix is singular (i.e., not invertible) and it has a rank of k
|
|
(i.e., it has k linearly independent columns), then it has an LU
|
|
decomposition if its first k leading principal minors are non-zero.
|
|
|
|
This algorithm will simply attempt to perform LU decomposition on any square
|
|
matrix and raise an error if no such decomposition exists.
|
|
|
|
Reference: https://en.wikipedia.org/wiki/LU_decomposition
|
|
"""
|
|
|
|
from __future__ import annotations
|
|
|
|
import numpy as np
|
|
|
|
|
|
def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
|
|
"""
|
|
Perform LU decomposition on a given matrix and raises an error if the matrix
|
|
isn't square or if no such decomposition exists
|
|
>>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
|
|
>>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
|
|
>>> lower_mat
|
|
array([[1. , 0. , 0. ],
|
|
[0. , 1. , 0. ],
|
|
[2.5, 8. , 1. ]])
|
|
>>> upper_mat
|
|
array([[ 2. , -2. , 1. ],
|
|
[ 0. , 1. , 2. ],
|
|
[ 0. , 0. , -17.5]])
|
|
|
|
>>> matrix = np.array([[4, 3], [6, 3]])
|
|
>>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
|
|
>>> lower_mat
|
|
array([[1. , 0. ],
|
|
[1.5, 1. ]])
|
|
>>> upper_mat
|
|
array([[ 4. , 3. ],
|
|
[ 0. , -1.5]])
|
|
|
|
# Matrix is not square
|
|
>>> matrix = np.array([[2, -2, 1], [0, 1, 2]])
|
|
>>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
|
|
Traceback (most recent call last):
|
|
...
|
|
ValueError: 'table' has to be of square shaped array but got a 2x3 array:
|
|
[[ 2 -2 1]
|
|
[ 0 1 2]]
|
|
|
|
# Matrix is invertible, but its first leading principal minor is 0
|
|
>>> matrix = np.array([[0, 1], [1, 0]])
|
|
>>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
|
|
Traceback (most recent call last):
|
|
...
|
|
ArithmeticError: No LU decomposition exists
|
|
|
|
# Matrix is singular, but its first leading principal minor is 1
|
|
>>> matrix = np.array([[1, 0], [1, 0]])
|
|
>>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
|
|
>>> lower_mat
|
|
array([[1., 0.],
|
|
[1., 1.]])
|
|
>>> upper_mat
|
|
array([[1., 0.],
|
|
[0., 0.]])
|
|
|
|
# Matrix is singular, but its first leading principal minor is 0
|
|
>>> matrix = np.array([[0, 1], [0, 1]])
|
|
>>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
|
|
Traceback (most recent call last):
|
|
...
|
|
ArithmeticError: No LU decomposition exists
|
|
"""
|
|
# Ensure that table is a square array
|
|
rows, columns = np.shape(table)
|
|
if rows != columns:
|
|
msg = (
|
|
"'table' has to be of square shaped array but got a "
|
|
f"{rows}x{columns} array:\n{table}"
|
|
)
|
|
raise ValueError(msg)
|
|
|
|
lower = np.zeros((rows, columns))
|
|
upper = np.zeros((rows, columns))
|
|
|
|
# in 'total', the necessary data is extracted through slices
|
|
# and the sum of the products is obtained.
|
|
|
|
for i in range(columns):
|
|
for j in range(i):
|
|
total = np.sum(lower[i, :i] * upper[:i, j])
|
|
if upper[j][j] == 0:
|
|
raise ArithmeticError("No LU decomposition exists")
|
|
lower[i][j] = (table[i][j] - total) / upper[j][j]
|
|
lower[i][i] = 1
|
|
for j in range(i, columns):
|
|
total = np.sum(lower[i, :i] * upper[:i, j])
|
|
upper[i][j] = table[i][j] - total
|
|
return lower, upper
|
|
|
|
|
|
if __name__ == "__main__":
|
|
import doctest
|
|
|
|
doctest.testmod()
|