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.

Author: Codebuff Contributor

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