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