|
1 |
| -"""Lower-Upper (LU) Decomposition.""" |
| 1 | +"""Lower-Upper (LU) Decomposition. |
2 | 2 |
|
3 |
| -# lower–upper (LU) decomposition - https://en.wikipedia.org/wiki/LU_decomposition |
4 |
| -import numpy |
| 3 | +Reference: |
| 4 | +- https://en.wikipedia.org/wiki/LU_decomposition |
| 5 | +""" |
| 6 | +from typing import Tuple |
5 | 7 |
|
| 8 | +import numpy as np |
| 9 | +from numpy import ndarray |
6 | 10 |
|
7 |
| -def LUDecompose(table): |
| 11 | + |
| 12 | +def lower_upper_decomposition(table: ndarray) -> Tuple[ndarray, ndarray]: |
| 13 | + """Lower-Upper (LU) Decomposition |
| 14 | +
|
| 15 | + Example: |
| 16 | +
|
| 17 | + >>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]]) |
| 18 | + >>> outcome = lower_upper_decomposition(matrix) |
| 19 | + >>> outcome[0] |
| 20 | + array([[1. , 0. , 0. ], |
| 21 | + [0. , 1. , 0. ], |
| 22 | + [2.5, 8. , 1. ]]) |
| 23 | + >>> outcome[1] |
| 24 | + array([[ 2. , -2. , 1. ], |
| 25 | + [ 0. , 1. , 2. ], |
| 26 | + [ 0. , 0. , -17.5]]) |
| 27 | +
|
| 28 | + >>> matrix = np.array([[2, -2, 1], [0, 1, 2]]) |
| 29 | + >>> lower_upper_decomposition(matrix) |
| 30 | + Traceback (most recent call last): |
| 31 | + ... |
| 32 | + ValueError: 'table' has to be of square shaped array but got a 2x3 array: |
| 33 | + [[ 2 -2 1] |
| 34 | + [ 0 1 2]] |
| 35 | + """ |
8 | 36 | # Table that contains our data
|
9 | 37 | # Table has to be a square array so we need to check first
|
10 |
| - rows, columns = numpy.shape(table) |
11 |
| - L = numpy.zeros((rows, columns)) |
12 |
| - U = numpy.zeros((rows, columns)) |
| 38 | + rows, columns = np.shape(table) |
13 | 39 | if rows != columns:
|
14 |
| - return [] |
| 40 | + raise ValueError( |
| 41 | + f"'table' has to be of square shaped array but got a {rows}x{columns} " |
| 42 | + + f"array:\n{table}" |
| 43 | + ) |
| 44 | + lower = np.zeros((rows, columns)) |
| 45 | + upper = np.zeros((rows, columns)) |
15 | 46 | for i in range(columns):
|
16 | 47 | for j in range(i):
|
17 |
| - sum = 0 |
| 48 | + total = 0 |
18 | 49 | for k in range(j):
|
19 |
| - sum += L[i][k] * U[k][j] |
20 |
| - L[i][j] = (table[i][j] - sum) / U[j][j] |
21 |
| - L[i][i] = 1 |
| 50 | + total += lower[i][k] * upper[k][j] |
| 51 | + lower[i][j] = (table[i][j] - total) / upper[j][j] |
| 52 | + lower[i][i] = 1 |
22 | 53 | for j in range(i, columns):
|
23 |
| - sum1 = 0 |
| 54 | + total = 0 |
24 | 55 | for k in range(i):
|
25 |
| - sum1 += L[i][k] * U[k][j] |
26 |
| - U[i][j] = table[i][j] - sum1 |
27 |
| - return L, U |
| 56 | + total += lower[i][k] * upper[k][j] |
| 57 | + upper[i][j] = table[i][j] - total |
| 58 | + return lower, upper |
28 | 59 |
|
29 | 60 |
|
30 | 61 | if __name__ == "__main__":
|
31 |
| - matrix = numpy.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]]) |
32 |
| - L, U = LUDecompose(matrix) |
33 |
| - print(L) |
34 |
| - print(U) |
| 62 | + import doctest |
| 63 | + |
| 64 | + doctest.testmod() |
0 commit comments