From 5ce63b5966b6ad9c7ce36c449fb31112c3e1d084 Mon Sep 17 00:00:00 2001 From: Tianyi Zheng Date: Sat, 1 Apr 2023 01:11:24 -0400 Subject: [PATCH] Fix `mypy` errors in `lu_decomposition.py` (attempt 2) (#8100) * updating DIRECTORY.md * Fix mypy errors in lu_decomposition.py * Replace for-loops with comprehensions * Add explanation of LU decomposition and extra doctests Add an explanation of LU decomposition with conditions for when an LU decomposition exists Add extra doctests to handle each of the possible conditions for when a decomposition exists/doesn't exist * updating DIRECTORY.md * updating DIRECTORY.md --------- Co-authored-by: github-actions <${GITHUB_ACTOR}@users.noreply.github.com> --- arithmetic_analysis/lu_decomposition.py | 91 ++++++++++++++++++------- 1 file changed, 65 insertions(+), 26 deletions(-) diff --git a/arithmetic_analysis/lu_decomposition.py b/arithmetic_analysis/lu_decomposition.py index 217719cf4..941c1dadf 100644 --- a/arithmetic_analysis/lu_decomposition.py +++ b/arithmetic_analysis/lu_decomposition.py @@ -1,62 +1,101 @@ -"""Lower-Upper (LU) Decomposition. +""" +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. -Reference: -- https://en.wikipedia.org/wiki/LU_decomposition +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 -from numpy import float64 -from numpy.typing import ArrayLike -def lower_upper_decomposition( - table: ArrayLike[float64], -) -> tuple[ArrayLike[float64], ArrayLike[float64]]: - """Lower-Upper (LU) Decomposition - - Example: - +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]]) - >>> outcome = lower_upper_decomposition(matrix) - >>> outcome[0] + >>> lower_mat, upper_mat = lower_upper_decomposition(matrix) + >>> lower_mat array([[1. , 0. , 0. ], [0. , 1. , 0. ], [2.5, 8. , 1. ]]) - >>> outcome[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_upper_decomposition(matrix) + >>> 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 """ - # Table that contains our data - # Table has to be a square array so we need to check first + # Ensure that table is a square array 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}" + f"'table' has to be of square shaped array but got a " + f"{rows}x{columns} 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] + total = sum(lower[i][k] * upper[k][j] for k in range(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 = 0 - for k in range(i): - total += lower[i][k] * upper[k][j] + total = sum(lower[i][k] * upper[k][j] for k in range(j)) upper[i][j] = table[i][j] - total return lower, upper