Refactoring and optimization of the lu_decomposition algorithm

[quant12345]

Describe your change:

Replacing the generator on numpy vector operations from lu_decomposition.

before: total = sum(lower[i][k] * upper[k][j] for k in range(j))

after: total = np.sum(lower[i, :i] * upper[:i, j])

In 'total', the necessary data is extracted through slices and the sum of the products is obtained.
Moreover, as the array size increases, the ratio of the calculation time of the original algorithm
growing towards the new. As a result, with an array size of n x n 1000, the calculation time for the
original algorithm is almost 7 minutes, and the new one is 12 seconds(time in tests is seconds).

The tests generate n x n arrays. To prevent it from triggering:

raise ArithmeticError("No LU decomposition exists")

diagonal elements are enlarged. Two arrays are extracted from the resulting tuple
rounded to the fifth digit and compared for identity.

code performance tests

"""
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.

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 datetime

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

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


def lower_upper_decomposition_new(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()

    arr = [18, 30, 50, 100, 200, 300, 500, 700, 1000]

    for i in range(len(arr)):
        n = arr[i]

        matrix = np.random.randint(low=-5, high=5, size=(n, n))
        matrix[np.diag_indices_from(matrix)] = 10 * n

        now = datetime.datetime.now()
        original = lower_upper_decomposition(matrix)
        time_original = datetime.datetime.now() - now

        now = datetime.datetime.now()
        new = lower_upper_decomposition_new(matrix)
        time_new = datetime.datetime.now() - now

        word = ('array size n x n {0} time_old {1} time_new {2}'
                ' difference in time {3} {4} array 1 {5} array 2 {6}'
                .format(n, time_original.total_seconds(),
                        time_new.total_seconds(), round(time_original / time_new, 2),
                        'comparison result',
                        np.all(np.round(original[0], 5) == np.round(new[0], 5)),
                        np.all(np.round(original[1], 5) == np.round(new[1], 5))))


        print(word)

Output:

array size n x n 18 time_old 0.006474 time_new 0.006151 difference in time 1.05 comparison result array 1 True array 2 True
array size n x n 30 time_old 0.012381 time_new 0.008019 difference in time 1.54 comparison result array 1 True array 2 True
array size n x n 50 time_old 0.055401 time_new 0.022404 difference in time 2.47 comparison result array 1 True array 2 True
array size n x n 100 time_old 0.418046 time_new 0.085158 difference in time 4.91 comparison result array 1 True array 2 True
array size n x n 200 time_old 3.175042 time_new 0.352767 difference in time 9.0 comparison result array 1 True array 2 True
array size n x n 300 time_old 10.780385 time_new 0.834303 difference in time 12.92 comparison result array 1 True array 2 True
array size n x n 500 time_old 50.109723 time_new 2.55773 difference in time 19.59 comparison result array 1 True array 2 True
array size n x n 700 time_old 138.859871 time_new 6.515611 difference in time 21.31 comparison result array 1 True array 2 True
array size n x n 1000 time_old 408.079446 time_new 11.721369 difference in time 34.81 comparison result array 1 True array 2 True 

• Add an algorithm?
• Fix a bug or typo in an existing algorithm?
• Documentation change?

Checklist:

• I have read CONTRIBUTING.md.
• This pull request is all my own work -- I have not plagiarized.
• I know that pull requests will not be merged if they fail the automated tests.
• This PR only changes one algorithm file. To ease review, please open separate PRs for separate algorithms.
• All new Python files are placed inside an existing directory.
• All filenames are in all lowercase characters with no spaces or dashes.
• All functions and variable names follow Python naming conventions.
• All function parameters and return values are annotated with Python type hints.
• All functions have doctests that pass the automated testing.
• All new algorithms include at least one URL that points to Wikipedia or another similar explanation.
• If this pull request resolves one or more open issues then the description above includes the issue number(s) with a closing keyword: "Fixes #ISSUE-NUMBER".