Python/arithmetic_analysis/lu_decomposition.py
arfy slowy 6732fa0131
[fixed] module 'numpy' is imported with both 'import' and 'import from' (#4544)
* [fixed] module 'numy' is imported with both 'import' and 'import from'

* remove commented
2021-07-26 14:52:52 +02:00

64 lines
1.8 KiB
Python

"""Lower-Upper (LU) Decomposition.
Reference:
- https://en.wikipedia.org/wiki/LU_decomposition
"""
from typing import Tuple
import numpy as np
def lower_upper_decomposition(table: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
"""Lower-Upper (LU) Decomposition
Example:
>>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
>>> outcome = lower_upper_decomposition(matrix)
>>> outcome[0]
array([[1. , 0. , 0. ],
[0. , 1. , 0. ],
[2.5, 8. , 1. ]])
>>> outcome[1]
array([[ 2. , -2. , 1. ],
[ 0. , 1. , 2. ],
[ 0. , 0. , -17.5]])
>>> matrix = np.array([[2, -2, 1], [0, 1, 2]])
>>> 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]]
"""
# Table that contains our data
# Table has to be a square array so we need to check first
rows, columns = np.shape(table)
if rows != columns:
raise ValueError(
f"'table' has to be of square shaped array but got a {rows}x{columns} "
+ f"array:\n{table}"
)
lower = np.zeros((rows, columns))
upper = np.zeros((rows, columns))
for i in range(columns):
for j in range(i):
total = 0
for k in range(j):
total += lower[i][k] * upper[k][j]
lower[i][j] = (table[i][j] - total) / upper[j][j]
lower[i][i] = 1
for j in range(i, columns):
total = 0
for k in range(i):
total += lower[i][k] * upper[k][j]
upper[i][j] = table[i][j] - total
return lower, upper
if __name__ == "__main__":
import doctest
doctest.testmod()