Lowerupper 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 nonzero see https:en.wikipedia.orgwikiMinorlinearalgebra 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 nonzero. 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.orgwikiLUdecomposition 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.array2, 2, 1, 0, 1, 2, 5, 3, 1  lowermat, uppermat  lowerupperdecompositionmatrix  lowermat array1. , 0. , 0. , 0. , 1. , 0. , 2.5, 8. , 1.   uppermat array  2. ,  2. ,   1. ,   0. ,   1. ,   2. ,   0. ,   0. , 17.5  matrix  np.array4, 3, 6, 3  lowermat, uppermat  lowerupperdecompositionmatrix  lowermat array1. , 0. , 1.5, 1.   uppermat array 4. ,  3. ,  0. , 1.5 Matrix is not square  matrix  np.array2, 2, 1, 0, 1, 2  lowermat, uppermat  lowerupperdecompositionmatrix 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.array0, 1, 1, 0  lowermat, uppermat  lowerupperdecompositionmatrix Traceback most recent call last: ... ArithmeticError: No LU decomposition exists Matrix is singular, but its first leading principal minor is 1  matrix  np.array1, 0, 1, 0  lowermat, uppermat  lowerupperdecompositionmatrix  lowermat array1., 0., 1., 1.  uppermat array1., 0., 0., 0. Matrix is singular, but its first leading principal minor is 0  matrix  np.array0, 1, 0, 1  lowermat, uppermat  lowerupperdecompositionmatrix Traceback most recent call last: ... ArithmeticError: No LU decomposition exists Ensure that table is a square array in 'total', the necessary data is extracted through slices and the sum of the products is obtained.
from __future__ import annotations

import numpy as np


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]])
    >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
    >>> lower_mat
    array([[1. , 0. , 0. ],
           [0. , 1. , 0. ],
           [2.5, 8. , 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_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
    """
    # Ensure that table is a square array
    rows, columns = np.shape(table)
    if rows != columns:
        msg = (
            "'table' has to be of square shaped array but got a "
            f"{rows}x{columns} array:\n{table}"
        )
        raise ValueError(msg)

    lower = np.zeros((rows, columns))
    upper = np.zeros((rows, columns))

    # in 'total', the necessary data is extracted through slices
    # and the sum of the products is obtained.

    for i in range(columns):
        for j in range(i):
            total = np.sum(lower[i, :i] * upper[:i, 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 = np.sum(lower[i, :i] * upper[:i, j])
            upper[i][j] = table[i][j] - total
    return lower, upper


if __name__ == "__main__":
    import doctest

    doctest.testmod()
