From 90064ed9363f27b755b07341bd7da5e0af0b267a Mon Sep 17 00:00:00 2001 From: MajaByr Date: Tue, 31 Mar 2026 21:36:32 +0200 Subject: [PATCH] :sparkles: Add LU decomposition --- CHANGELOG.md | 6 +++ matrix_calculus/lu_decomposition.py | 71 +++++++++++++++++++++++++++++ pyproject.toml | 2 +- test.sh | 2 + tests/lu_decomposition_test.py | 54 ++++++++++++++++++++++ 5 files changed, 134 insertions(+), 1 deletion(-) create mode 100644 matrix_calculus/lu_decomposition.py create mode 100644 tests/lu_decomposition_test.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 09de071..93b934d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,9 @@ +## [1.4.0] - 2026-03-31 + +### Added + +- LU decomposition + ## [1.3.0] - 2026-03-20 ### Added diff --git a/matrix_calculus/lu_decomposition.py b/matrix_calculus/lu_decomposition.py new file mode 100644 index 0000000..43f0ef9 --- /dev/null +++ b/matrix_calculus/lu_decomposition.py @@ -0,0 +1,71 @@ +from __future__ import annotations + +import numpy as np + + +class LUDecomposition: + @staticmethod + def decompose( + matrix: list[list[float]], + pivoting=False, + ) -> tuple: + if pivoting is True: + return LUDecomposition.decompose_with_pivot(matrix=matrix) + else: + return LUDecomposition.decompose_no_pivot(matrix=matrix) + + @staticmethod + def decompose_no_pivot( + matrix: list[list[float]], + ) -> tuple[list[list[float]], list[list[float]]]: + n = len(matrix) + L = [[0.0] * n for _ in range(n)] + U = [[0.0] * n for _ in range(n)] + + for i in range(n): + for k in range(i, n): + sum_u = sum(L[i][j] * U[j][k] for j in range(i)) + U[i][k] = matrix[i][k] - sum_u + L[i][i] = 1.0 + for k in range(i + 1, n): + sum_l = sum(L[k][j] * U[j][i] for j in range(i)) + if U[i][i] == 0: + raise ZeroDivisionError( + "U diagonal element equals zero. Cannot decompose without pivoting.", + ) + L[k][i] = (matrix[k][i] - sum_l) / U[i][i] + + return L, U + + @staticmethod + def decompose_with_pivot( + matrix: list[list[float]], + ) -> tuple[list[list[float]], list[list[float]], list[list[float]]]: + A = np.array(matrix).copy().astype(float) + n = A.shape[0] + + L = np.zeros((n, n)) + U = A.copy() + P = np.eye(n) + + for k in range(n): + pivot = np.argmax(np.abs(U[k:n, k])) + k + if U[pivot, k] < 1e-5: + raise ZeroDivisionError( + "Matrix is singular. Cannot use LU decomposition with pivoting.", + ) + + # Switch rows + U[[k, pivot], :] = U[[pivot, k], :] + P[[k, pivot], :] = P[[pivot, k], :] + if k > 0: + L[[k, pivot], :k] = L[[pivot, k], :k] + + for j in range(k + 1, n): + L[j, k] = U[j, k] / U[k, k] + U[j, :] = U[j, :] - L[j, k] * U[k, :] + + for i in range(n): + L[i, i] = 1.0 + + return P, L, U diff --git a/pyproject.toml b/pyproject.toml index a302f88..de8156f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "matrix_calculus" -version = "1.3.0" +version = "1.4.0" description = "Matrix calculus toolkit" requires-python = ">=3.10" dependencies = ["numpy", "seaborn", "matplotlib", "weasyprint", "markdown"] diff --git a/test.sh b/test.sh index 62b8159..8d09de1 100755 --- a/test.sh +++ b/test.sh @@ -5,6 +5,8 @@ echo "=== Activating virtual environment ===" source .venv/bin/activate echo "=== Running tests ===" + +# python tests TEST_DIR="tests" pytest -v "$TEST_DIR" --maxfail=1 --disable-warnings --tb=short --cov=. --cov-report=term-missing diff --git a/tests/lu_decomposition_test.py b/tests/lu_decomposition_test.py new file mode 100644 index 0000000..6b0c294 --- /dev/null +++ b/tests/lu_decomposition_test.py @@ -0,0 +1,54 @@ +from __future__ import annotations + +import numpy as np +import pytest # noqa: F401 + +from matrix_calculus.lu_decomposition import LUDecomposition +from matrix_calculus.matrix_multiplier import MatrixMultiplier + + +def test_lu_no_pivot(): + matrix = [[1, 2, 3], [5, 6, 7], [9, 10, 11]] + L, U = LUDecomposition.decompose(matrix=matrix, pivoting=False) + assert L == [[1, 0, 0], [5, 1, 0], [9, 2, 1]] + assert U == [[1, 2, 3], [0, -4, -8], [0, 0, 0]] + + +def test_lu_with_pivot_singular(): + # Test singular matrix (determinant equal to 0) + with pytest.raises(ZeroDivisionError): + matrix = [[1, 2, 3], [5, 6, 7], [9, 10, 11]] + LUDecomposition.decompose(matrix=matrix, pivoting=True) + + +def test_lu_with_pivot_correct(): + # Test correct decomposition + matrix = [[1, 2, 3], [4, 5, 6], [8, 10, 10]] + P, L, U = LUDecomposition.decompose(matrix=matrix, pivoting=True) + L == [[1, 0, 0], [4, 1, 0], [8, 2, 1]] + assert np.allclose( + MatrixMultiplier.traditional_multiplication(P, matrix), + MatrixMultiplier.traditional_multiplication(L, U), + ) + + +def test_lu_example(): + date = 8 + 24 + np.random.seed(date) + matrix = np.random.rand(date, date) * 10 // 1 + # Make sure that generated matrix is not singular + matrix += np.eye(date) * 100 + + # LU without pivoting + L, U = LUDecomposition.decompose(matrix, pivoting=False) + assert np.allclose( + MatrixMultiplier.traditional_multiplication(L, U), + matrix, + ) + + # LU with pivoting + P, L, U = LUDecomposition.decompose(matrix, pivoting=True) + assert np.allclose( + MatrixMultiplier.traditional_multiplication(L, U), + MatrixMultiplier.traditional_multiplication(P, matrix), + )