Gaussian elimination refactoring numpy and optimization.

[quant12345]

Describe your change:

Instead of a nested loop, vector operations (refactoring numpy and optimization) were used.
To replace a nested loop, access to array values is used through slices. With a 'coefficients' array size of up to 1900, it works faster than with two loops.

code performance tests

""" Gaussian elimination method for solving a system of linear equations. Gaussian elimination - https://en.wikipedia.org/wiki/Gaussian_elimination """ import numpy as np from numpy import float64 from numpy.typing import NDArray import datetime def retroactive_resolution( coefficients: NDArray[float64], vector: NDArray[float64] ) -> NDArray[float64]: """ This function performs a retroactive linear system resolution for triangular matrix Examples: 2x1 + 2x2 - 1x3 = 5 2x1 + 2x2 = -1 0x1 - 2x2 - 1x3 = -7 0x1 - 2x2 = -1 0x1 + 0x2 + 5x3 = 15 >>> gaussian_elimination([[2, 2, -1], [0, -2, -1], [0, 0, 5]], [[5], [-7], [15]]) array([[2.], [2.], [3.]]) >>> gaussian_elimination([[2, 2], [0, -2]], [[-1], [-1]]) array([[-1. ], [ 0.5]]) """ rows, columns = np.shape(coefficients) x: NDArray[float64] = np.zeros((rows, 1), dtype=float) for row in reversed(range(rows)): total = np.dot(coefficients[row, row + 1:], x[row + 1:]) x[row, 0] = (vector[row] - total) / coefficients[row, row] return x def gaussian_elimination_new( coefficients: NDArray[float64], vector: NDArray[float64] ) -> NDArray[float64]: """ This function performs Gaussian elimination method Examples: 1x1 - 4x2 - 2x3 = -2 1x1 + 2x2 = 5 5x1 + 2x2 - 2x3 = -3 5x1 + 2x2 = 5 1x1 - 1x2 + 0x3 = 4 >>> gaussian_elimination([[1, -4, -2], [5, 2, -2], [1, -1, 0]], [[-2], [-3], [4]]) array([[ 2.3 ], [-1.7 ], [ 5.55]]) >>> gaussian_elimination([[1, 2], [5, 2]], [[5], [5]]) array([[0. ], [2.5]]) """ # coefficients must to be a square matrix so we need to check first rows, columns = np.shape(coefficients) if rows != columns: return np.array((), dtype=float) # augmented matrix augmented_mat: NDArray[float64] = np.concatenate((coefficients, vector), axis=1) augmented_mat = augmented_mat.astype("float64") # scale the matrix leaving it triangular # factor - taken elements from 'augmented_mat' by row slice: row + 1: columns, # with 'row' column number, divided by the selected 'diagonal_element'. # [:, np.newaxis] - array is converted to two dimensions to use vector multiplication. # product of factor and the row selected from augmented_mat[row, :] is subtracted from # the value of the selected 'rows: row + 1: columns' and all columns: # augmented_mat[row + 1: columns, :] for row in range(rows - 1): diagonal_element = augmented_mat[row, row] factor = (augmented_mat[row + 1: columns, row] / diagonal_element)[:, np.newaxis] augmented_mat[row + 1: columns, :] -= factor * augmented_mat[row, :] x = retroactive_resolution( augmented_mat[:, 0:columns], augmented_mat[:, columns: columns + 1] ) return x def gaussian_elimination( coefficients: NDArray[float64], vector: NDArray[float64] ) -> NDArray[float64]: """ This function performs Gaussian elimination method Examples: 1x1 - 4x2 - 2x3 = -2 1x1 + 2x2 = 5 5x1 + 2x2 - 2x3 = -3 5x1 + 2x2 = 5 1x1 - 1x2 + 0x3 = 4 >>> gaussian_elimination([[1, -4, -2], [5, 2, -2], [1, -1, 0]], [[-2], [-3], [4]]) array([[ 2.3 ], [-1.7 ], [ 5.55]]) >>> gaussian_elimination([[1, 2], [5, 2]], [[5], [5]]) array([[0. ], [2.5]]) """ # coefficients must to be a square matrix so we need to check first rows, columns = np.shape(coefficients) if rows != columns: return np.array((), dtype=float) # augmented matrix augmented_mat: NDArray[float64] = np.concatenate((coefficients, vector), axis=1) augmented_mat = augmented_mat.astype("float64") # scale the matrix leaving it triangular for row in range(rows - 1): pivot = augmented_mat[row, row] # factor_t = augmented_mat[aaa[0], aaa[1]] / pivot for col in range(row + 1, columns): factor = augmented_mat[col, row] / pivot augmented_mat[col, :] -= factor * augmented_mat[row, :] x = retroactive_resolution( augmented_mat[:, 0:columns], augmented_mat[:, columns: columns + 1] ) return x if __name__ == "__main__": import doctest doctest.testmod() arr = [30, 50, 100, 300, 500, 700, 1000, 1500, 1900] for i in range(len(arr)): n = arr[i] coefficient = np.random.randint(low=-5, high=5, size=(n, n)) coefficient[np.diag_indices_from(coefficient)] = 10 * n vector = np.random.randint(low=-5, high=5, size=(n, 1)) now = datetime.datetime.now() original = gaussian_elimination(coefficient, vector) time_original = datetime.datetime.now() - now now = datetime.datetime.now() new = gaussian_elimination_new(coefficient, vector) time_new = datetime.datetime.now() - now word = ('array size n x n {0} time_old {1} time_new {2}' ' difference in time {3} comparison result array {4}' .format(n, time_original.total_seconds(), time_new.total_seconds(), round(time_original / time_new, 2), np.all(np.round(original, 5) == np.round(new, 5)))) print(word) 

At each iteration, arrays are generated: coefficient(n x n), vector(n x 1). Diagonal elements are made large 10 * n so that division by zero (diagonal_element) does not occur. The time spent by each algorithm in seconds, new and old, is displayed and the results are compared for equality.

Output:

array size n x n 30 time_old 0.00167 time_new 0.000522 difference in time 3.2 comparison result array True array size n x n 50 time_old 0.004336 time_new 0.000888 difference in time 4.88 comparison result array True array size n x n 100 time_old 0.017354 time_new 0.002621 difference in time 6.62 comparison result array True array size n x n 300 time_old 0.166007 time_new 0.048619 difference in time 3.41 comparison result array True array size n x n 500 time_old 0.581854 time_new 0.201526 difference in time 2.89 comparison result array True array size n x n 700 time_old 1.157788 time_new 0.6258 difference in time 1.85 comparison result array True array size n x n 1000 time_old 2.693735 time_new 1.883084 difference in time 1.43 comparison result array True array size n x n 1500 time_old 7.471813 time_new 6.979876 difference in time 1.07 comparison result array True array size n x n 1900 time_old 13.326552 time_new 13.648966 difference in time 0.98 comparison result array True 

Just in case, I’m attaching the result with a size of 5000:
array size n x n 5000 time_old 166.356578 time_new 307.890158 difference in time 0.54 comparison result array True

Need to decide what is more important: using a larger array (what size is usually used) or refactoring and using an array up to 1900 in size (if larger, the calculation will take longer)?

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

[tianyizheng02]

Thank you for your contribution, but unfortunately we're not accepting new PRs at the moment. We're currently trying to clear our backlog of existing PRs in preparation for Hacktoberfest. If you want to contribute, please wait until after October 1 to open a new PR.