|
1 |
| -"""Lower-Upper (LU) Decomposition. |
| 1 | +""" |
| 2 | +Lower–upper (LU) decomposition factors a matrix as a product of a lower |
| 3 | +triangular matrix and an upper triangular matrix. A square matrix has an LU |
| 4 | +decomposition under the following conditions: |
| 5 | + - If the matrix is invertible, then it has an LU decomposition if and only |
| 6 | + if all of its leading principal minors are non-zero (see |
| 7 | + https://en.wikipedia.org/wiki/Minor_(linear_algebra) for an explanation of |
| 8 | + leading principal minors of a matrix). |
| 9 | + - If the matrix is singular (i.e., not invertible) and it has a rank of k |
| 10 | + (i.e., it has k linearly independent columns), then it has an LU |
| 11 | + decomposition if its first k leading principal minors are non-zero. |
| 12 | +
|
| 13 | +This algorithm will simply attempt to perform LU decomposition on any square |
| 14 | +matrix and raise an error if no such decomposition exists. |
2 | 15 |
|
3 |
| -Reference: |
4 |
| -- https://en.wikipedia.org/wiki/LU_decomposition |
| 16 | +Reference: https://en.wikipedia.org/wiki/LU_decomposition |
5 | 17 | """
|
6 | 18 | from __future__ import annotations
|
7 | 19 |
|
8 | 20 | import numpy as np
|
9 |
| -from numpy import float64 |
10 |
| -from numpy.typing import ArrayLike |
11 |
| - |
12 | 21 |
|
13 |
| -def lower_upper_decomposition( |
14 |
| - table: ArrayLike[float64], |
15 |
| -) -> tuple[ArrayLike[float64], ArrayLike[float64]]: |
16 |
| - """Lower-Upper (LU) Decomposition |
17 |
| -
|
18 |
| - Example: |
19 | 22 |
|
| 23 | +def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray]: |
| 24 | + """ |
| 25 | + Perform LU decomposition on a given matrix and raises an error if the matrix |
| 26 | + isn't square or if no such decomposition exists |
20 | 27 | >>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
|
21 |
| - >>> outcome = lower_upper_decomposition(matrix) |
22 |
| - >>> outcome[0] |
| 28 | + >>> lower_mat, upper_mat = lower_upper_decomposition(matrix) |
| 29 | + >>> lower_mat |
23 | 30 | array([[1. , 0. , 0. ],
|
24 | 31 | [0. , 1. , 0. ],
|
25 | 32 | [2.5, 8. , 1. ]])
|
26 |
| - >>> outcome[1] |
| 33 | + >>> upper_mat |
27 | 34 | array([[ 2. , -2. , 1. ],
|
28 | 35 | [ 0. , 1. , 2. ],
|
29 | 36 | [ 0. , 0. , -17.5]])
|
30 | 37 |
|
| 38 | + >>> matrix = np.array([[4, 3], [6, 3]]) |
| 39 | + >>> lower_mat, upper_mat = lower_upper_decomposition(matrix) |
| 40 | + >>> lower_mat |
| 41 | + array([[1. , 0. ], |
| 42 | + [1.5, 1. ]]) |
| 43 | + >>> upper_mat |
| 44 | + array([[ 4. , 3. ], |
| 45 | + [ 0. , -1.5]]) |
| 46 | +
|
| 47 | + # Matrix is not square |
31 | 48 | >>> matrix = np.array([[2, -2, 1], [0, 1, 2]])
|
32 |
| - >>> lower_upper_decomposition(matrix) |
| 49 | + >>> lower_mat, upper_mat = lower_upper_decomposition(matrix) |
33 | 50 | Traceback (most recent call last):
|
34 | 51 | ...
|
35 | 52 | ValueError: 'table' has to be of square shaped array but got a 2x3 array:
|
36 | 53 | [[ 2 -2 1]
|
37 | 54 | [ 0 1 2]]
|
| 55 | +
|
| 56 | + # Matrix is invertible, but its first leading principal minor is 0 |
| 57 | + >>> matrix = np.array([[0, 1], [1, 0]]) |
| 58 | + >>> lower_mat, upper_mat = lower_upper_decomposition(matrix) |
| 59 | + Traceback (most recent call last): |
| 60 | + ... |
| 61 | + ArithmeticError: No LU decomposition exists |
| 62 | +
|
| 63 | + # Matrix is singular, but its first leading principal minor is 1 |
| 64 | + >>> matrix = np.array([[1, 0], [1, 0]]) |
| 65 | + >>> lower_mat, upper_mat = lower_upper_decomposition(matrix) |
| 66 | + >>> lower_mat |
| 67 | + array([[1., 0.], |
| 68 | + [1., 1.]]) |
| 69 | + >>> upper_mat |
| 70 | + array([[1., 0.], |
| 71 | + [0., 0.]]) |
| 72 | +
|
| 73 | + # Matrix is singular, but its first leading principal minor is 0 |
| 74 | + >>> matrix = np.array([[0, 1], [0, 1]]) |
| 75 | + >>> lower_mat, upper_mat = lower_upper_decomposition(matrix) |
| 76 | + Traceback (most recent call last): |
| 77 | + ... |
| 78 | + ArithmeticError: No LU decomposition exists |
38 | 79 | """
|
39 |
| - # Table that contains our data |
40 |
| - # Table has to be a square array so we need to check first |
| 80 | + # Ensure that table is a square array |
41 | 81 | rows, columns = np.shape(table)
|
42 | 82 | if rows != columns:
|
43 | 83 | raise ValueError(
|
44 |
| - f"'table' has to be of square shaped array but got a {rows}x{columns} " |
45 |
| - + f"array:\n{table}" |
| 84 | + f"'table' has to be of square shaped array but got a " |
| 85 | + f"{rows}x{columns} array:\n{table}" |
46 | 86 | )
|
| 87 | + |
47 | 88 | lower = np.zeros((rows, columns))
|
48 | 89 | upper = np.zeros((rows, columns))
|
49 | 90 | for i in range(columns):
|
50 | 91 | for j in range(i):
|
51 |
| - total = 0 |
52 |
| - for k in range(j): |
53 |
| - total += lower[i][k] * upper[k][j] |
| 92 | + total = sum(lower[i][k] * upper[k][j] for k in range(j)) |
| 93 | + if upper[j][j] == 0: |
| 94 | + raise ArithmeticError("No LU decomposition exists") |
54 | 95 | lower[i][j] = (table[i][j] - total) / upper[j][j]
|
55 | 96 | lower[i][i] = 1
|
56 | 97 | for j in range(i, columns):
|
57 |
| - total = 0 |
58 |
| - for k in range(i): |
59 |
| - total += lower[i][k] * upper[k][j] |
| 98 | + total = sum(lower[i][k] * upper[k][j] for k in range(j)) |
60 | 99 | upper[i][j] = table[i][j] - total
|
61 | 100 | return lower, upper
|
62 | 101 |
|
|
0 commit comments