✨ Add LU decomposition

CHANGELOG.md

@@ -1,3 +1,9 @@
+## [1.4.0] - 2026-03-31
+
+### Added
+
+- LU decomposition
+
 ## [1.3.0] - 2026-03-20
 
 ### Added

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

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

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
 

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