From c4269e067c28773b970d7b21183a4c3d50eb24a9 Mon Sep 17 00:00:00 2001 From: Codebuff Contributor Date: Mon, 18 May 2026 21:36:21 +0600 Subject: [PATCH 1/4] feat(maths): add LU decomposition algorithm for matrix factorization This commit introduces a complete implementation of the LU decomposition algorithm in the maths module, providing a fundamental numerical linear algebra tool for matrix factorization and solving systems of linear equations efficiently. LU Decomposition Overview: LU decomposition factorizes a square matrix A into the product of two triangular matrices: a lower triangular matrix L and an upper triangular matrix U, such that A = L * U. This decomposition is named after the mathematician Tadeusz Banachiewicz who formalized it in 1938, though the underlying algorithm traces back to the work of Alan Turing in 1948. The Doolittle algorithm implementation used in this contribution produces an L matrix with ones on the main diagonal (unit lower triangular) and a general upper triangular U matrix. This is one of the most commonly used variants of LU decomposition and is particularly well-suited for solving systems of linear equations. Mathematical Foundation: Given a square matrix A of size n x n, the decomposition computes: A = L * U where: - L is a lower triangular matrix with L[i][i] = 1 for all i - U is an upper triangular matrix with U[i][j] = 0 for all i > j The algorithm computes elements of L and U iteratively: For each column k from 0 to n-1: 1. Compute the k-th row of U: U[k][j] = A[k][j] - sum(L[k][s] * U[s][j] for s < k) 2. Compute the k-th column of L: L[i][k] = (A[i][k] - sum(L[i][s] * U[s][k] for s < k)) / U[k][k] The algorithm requires that all pivot elements U[k][k] be non-zero. If a zero pivot is encountered, the matrix may be singular or may require partial pivoting (row exchanges) for numerical stability. Implementation Details: The module provides two main functions: 1. lu_decomposition(matrix): Takes a square matrix as input and returns the (L, U) tuple. Includes comprehensive input validation: - Checks that the matrix is square - Detects zero pivots and raises informative errors - Converts all input values to floats for consistent arithmetic 2. solve_with_lu(lower, upper, b): Solves the system Ax = b given the LU decomposition of A. Uses a two-step approach: - Forward substitution to solve L * y = b - Back substitution to solve U * x = y Both functions include detailed docstrings following PEP 257 conventions, with comprehensive doctests that verify correctness across multiple test cases including: - Standard 2x2 and 3x3 matrices - Identity matrix (trivial case) - Singular matrix detection (zero pivot) - Non-square matrix validation - System of equations solving Algorithmic Complexity: Time Complexity: O(n^3) for the decomposition, O(n^2) for solving a system given the decomposition. This is the same asymptotic complexity as Gaussian elimination, but LU decomposition has the advantage that once the factorization is computed, solving multiple systems with the same coefficient matrix only requires O(n^2) work per system instead of O(n^3). Space Complexity: O(n^2) for storing the L and U matrices. Practical Applications: LU decomposition is used extensively in: - Engineering simulations (finite element analysis, CFD) - Computer graphics (transformation matrices) - Machine learning (linear regression, covariance computation) - Economics (input-output models, Leontief models) - Circuit analysis (nodal and mesh analysis) - Structural engineering (stiffness matrix solutions) This implementation provides an educational reference for understanding the algorithm while being functional enough for small to medium-sized matrices in practical applications. --- maths/lu_decomposition.py | 169 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 169 insertions(+) create mode 100644 maths/lu_decomposition.py diff --git a/maths/lu_decomposition.py b/maths/lu_decomposition.py new file mode 100644 index 000000000000..2009009d634b --- /dev/null +++ b/maths/lu_decomposition.py @@ -0,0 +1,169 @@ +""" +LU Decomposition + +Decomposes a square matrix into a lower triangular matrix (L) and an +upper triangular matrix (U) such that A = L * U. + +This decomposition is useful for: +- Solving systems of linear equations efficiently +- Computing matrix determinants +- Finding matrix inverses +- Repeated solving with the same coefficient matrix + +Reference: https://en.wikipedia.org/wiki/LU_decomposition +""" + + +def lu_decomposition(matrix: list[list[float]]) -> tuple[list[list[float]], list[list[float]]]: + """Perform LU decomposition on a square matrix. + + Decomposes the input matrix A into L (lower triangular) and U (upper + triangular) matrices such that A = L * U. The diagonal of L contains + all ones (Doolittle algorithm). + + The algorithm proceeds by iterating through each column and computing + the elements of U and L using the following formulas: + + For U (upper triangular): + U[k][j] = A[k][j] - sum(L[k][s] * U[s][j] for s in range(k)) + + For L (lower triangular): + L[i][k] = (A[i][k] - sum(L[i][s] * U[s][k] for s in range(k))) / U[k][k] + + Args: + matrix: A square matrix represented as a list of lists of floats. + + Returns: + A tuple (L, U) where L is a lower triangular matrix with ones on + the diagonal, and U is an upper triangular matrix. + + Raises: + ValueError: If the matrix is not square or if a zero pivot is + encountered (matrix is singular or requires pivoting). + + Examples: + >>> lu_decomposition([[2, 1, 1], [4, 3, 3], [8, 7, 9]]) + ([[1.0, 0.0, 0.0], [2.0, 1.0, 0.0], [4.0, 3.0, 1.0]], [[2.0, 1.0, 1.0], [0.0, 1.0, 1.0], [0.0, 0.0, 2.0]]) + + >>> lu_decomposition([[1, 2], [3, 4]]) + ([[1.0, 0.0], [3.0, 1.0]], [[1.0, 2.0], [0.0, -2.0]]) + + >>> lu_decomposition([[4, 3], [6, 3]]) + ([[1.0, 0.0], [1.5, 1.0]], [[4.0, 3.0], [0.0, -1.5]]) + + >>> lu_decomposition([[1, 0], [0, 1]]) + ([[1.0, 0.0], [0.0, 1.0]], [[1.0, 0.0], [0.0, 1.0]]) + + >>> lu_decomposition([[0, 1], [1, 0]]) + Traceback (most recent call last): + ... + ValueError: Zero pivot encountered. Matrix may be singular or require partial pivoting. + + >>> lu_decomposition([[1, 2, 3], [4, 5, 6]]) + Traceback (most recent call last): + ... + ValueError: Matrix must be square. + """ + n = len(matrix) + + # Validate that the matrix is square + if any(len(row) != n for row in matrix): + raise ValueError("Matrix must be square.") + + # Convert matrix to float values + a: list[list[float]] = [[float(val) for val in row] for row in matrix] + + # Initialize L with zeros and set diagonal to 1 (Doolittle algorithm) + lower: list[list[float]] = [[0.0] * n for _ in range(n)] + for i in range(n): + lower[i][i] = 1.0 + + # Initialize U with zeros + upper: list[list[float]] = [[0.0] * n for _ in range(n)] + + for k in range(n): + # Compute the k-th row of U + for j in range(k, n): + sum_val = sum(lower[k][s] * upper[s][j] for s in range(k)) + upper[k][j] = a[k][j] - sum_val + + # Check for zero pivot + if upper[k][k] == 0: + raise ValueError( + "Zero pivot encountered. Matrix may be singular or require partial pivoting." + ) + + # Compute the k-th column of L + for i in range(k + 1, n): + sum_val = sum(lower[i][s] * upper[s][k] for s in range(k)) + lower[i][k] = (matrix[i][k] - sum_val) / upper[k][k] + + return lower, upper + + +def solve_with_lu( + lower: list[list[float]], upper: list[list[float]], b: list[float] +) -> list[float]: + """Solve a system of linear equations Ax = b using LU decomposition. + + Given the LU decomposition of A (where A = L * U), this function + solves the system in two steps: + + 1. Forward substitution: Solve L * y = b for y + 2. Back substitution: Solve U * x = y for x + + Args: + lower: The lower triangular matrix L from LU decomposition. + upper: The upper triangular matrix U from LU decomposition. + b: The right-hand side vector of the system. + + Returns: + The solution vector x. + + Examples: + >>> l, u = lu_decomposition([[2, 1], [1, 3]]) + >>> solve_with_lu(l, u, [5, 7]) + [1.6, 1.8] + + >>> l, u = lu_decomposition([[1, 1, 1], [2, 3, 1], [1, 1, 2]]) + >>> solve_with_lu(l, u, [6, 11, 7]) + [5.0, 0.0, 1.0] + """ + n = len(b) + + # Forward substitution: solve L * y = b + y: list[float] = [0.0] * n + for i in range(n): + y[i] = b[i] - sum(lower[i][j] * y[j] for j in range(i)) + + # Back substitution: solve U * x = y + x: list[float] = [0.0] * n + for i in range(n - 1, -1, -1): + x[i] = (y[i] - sum(upper[i][j] * x[j] for j in range(i + 1, n))) / upper[i][i] + + return x + + +if __name__ == "__main__": + import doctest + + doctest.testmod() + + # Demonstration: solve a system of equations + # 2x + y + z = 8 + # 4x + 3y + 3z = 20 + # 8x + 7y + 9z = 46 + A = [[2, 1, 1], [4, 3, 3], [8, 7, 9]] + b = [8, 20, 46] + + L, U = lu_decomposition(A) + print("L matrix:") + for row in L: + print([f"{val:.2f}" for val in row]) + + print("\nU matrix:") + for row in U: + print([f"{val:.2f}" for val in row]) + + solution = solve_with_lu(L, U, b) + print(f"\nSolution: x = {[f'{val:.2f}' for val in solution]}") From 0135a70da18cd6309c9804ed335af1fc2e3c5e68 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 18 May 2026 15:43:23 +0000 Subject: [PATCH 2/4] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- maths/lu_decomposition.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/maths/lu_decomposition.py b/maths/lu_decomposition.py index 2009009d634b..3e4de571a9c0 100644 --- a/maths/lu_decomposition.py +++ b/maths/lu_decomposition.py @@ -14,7 +14,9 @@ """ -def lu_decomposition(matrix: list[list[float]]) -> tuple[list[list[float]], list[list[float]]]: +def lu_decomposition( + matrix: list[list[float]], +) -> tuple[list[list[float]], list[list[float]]]: """Perform LU decomposition on a square matrix. Decomposes the input matrix A into L (lower triangular) and U (upper From b8a2aa587e87e8318973d943ec935f21b7095f4b Mon Sep 17 00:00:00 2001 From: Md Mushfiqur Rahim <20mahin2020@gmail.com> Date: Mon, 18 May 2026 22:09:46 +0600 Subject: [PATCH 3/4] feat(maths): add LU decomposition algorithm for matrix factorization --- maths/lu_decomposition.py | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/maths/lu_decomposition.py b/maths/lu_decomposition.py index 3e4de571a9c0..65724da7f439 100644 --- a/maths/lu_decomposition.py +++ b/maths/lu_decomposition.py @@ -14,9 +14,12 @@ """ +Matrix = list[list[float]] + + def lu_decomposition( - matrix: list[list[float]], -) -> tuple[list[list[float]], list[list[float]]]: + matrix: Matrix, +) -> tuple[Matrix, Matrix]: """Perform LU decomposition on a square matrix. Decomposes the input matrix A into L (lower triangular) and U (upper @@ -44,8 +47,13 @@ def lu_decomposition( encountered (matrix is singular or requires pivoting). Examples: - >>> lu_decomposition([[2, 1, 1], [4, 3, 3], [8, 7, 9]]) - ([[1.0, 0.0, 0.0], [2.0, 1.0, 0.0], [4.0, 3.0, 1.0]], [[2.0, 1.0, 1.0], [0.0, 1.0, 1.0], [0.0, 0.0, 2.0]]) + >>> l, u = lu_decomposition( + ... [[2, 1, 1], [4, 3, 3], [8, 7, 9]] + ... ) # doctest: +NORMALIZE_WHITESPACE + >>> l + [[1.0, 0.0, 0.0], [2.0, 1.0, 0.0], [4.0, 3.0, 1.0]] + >>> u + [[2.0, 1.0, 1.0], [0.0, 1.0, 1.0], [0.0, 0.0, 2.0]] >>> lu_decomposition([[1, 2], [3, 4]]) ([[1.0, 0.0], [3.0, 1.0]], [[1.0, 2.0], [0.0, -2.0]]) @@ -56,10 +64,10 @@ def lu_decomposition( >>> lu_decomposition([[1, 0], [0, 1]]) ([[1.0, 0.0], [0.0, 1.0]], [[1.0, 0.0], [0.0, 1.0]]) - >>> lu_decomposition([[0, 1], [1, 0]]) + >>> lu_decomposition([[0, 1], [1, 0]]) # doctest: +ELLIPSIS Traceback (most recent call last): ... - ValueError: Zero pivot encountered. Matrix may be singular or require partial pivoting. + ValueError: Zero pivot encountered... >>> lu_decomposition([[1, 2, 3], [4, 5, 6]]) Traceback (most recent call last): @@ -92,13 +100,14 @@ def lu_decomposition( # Check for zero pivot if upper[k][k] == 0: raise ValueError( - "Zero pivot encountered. Matrix may be singular or require partial pivoting." + "Zero pivot encountered. " + "Matrix may be singular or require pivoting." ) # Compute the k-th column of L for i in range(k + 1, n): sum_val = sum(lower[i][s] * upper[s][k] for s in range(k)) - lower[i][k] = (matrix[i][k] - sum_val) / upper[k][k] + lower[i][k] = (a[i][k] - sum_val) / upper[k][k] return lower, upper From f21408763258b95bb1247c2bdd5b933e3e379e1f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 19 May 2026 01:15:43 +0000 Subject: [PATCH 4/4] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- maths/lu_decomposition.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/maths/lu_decomposition.py b/maths/lu_decomposition.py index 65724da7f439..8c9b2c5ff227 100644 --- a/maths/lu_decomposition.py +++ b/maths/lu_decomposition.py @@ -13,7 +13,6 @@ Reference: https://en.wikipedia.org/wiki/LU_decomposition """ - Matrix = list[list[float]] @@ -100,8 +99,7 @@ def lu_decomposition( # Check for zero pivot if upper[k][k] == 0: raise ValueError( - "Zero pivot encountered. " - "Matrix may be singular or require pivoting." + "Zero pivot encountered. Matrix may be singular or require pivoting." ) # Compute the k-th column of L