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

@@ -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