LU Decomposition of a Square Matrix

#333 · Linear Algebra · Medium

Problem

Perform LU decomposition of a square matrix. Decompose matrix A into a lower triangular matrix L and an upper triangular matrix U such that A = L * U (or PA = LU with partial pivoting).

Solution

from typing import List, Tuple

def lu_decomposition(
    matrix: List[List[float]]
) -> Tuple[List[List[float]], List[List[float]]]:
    n = len(matrix)
    L = [[0.0] * n for _ in range(n)]
    U = [row[:] for row in matrix]

    for i in range(n):
        L[i][i] = 1.0

    for col in range(n):
        # Partial pivoting
        max_row = col
        for row in range(col + 1, n):
            if abs(U[row][col]) > abs(U[max_row][col]):
                max_row = row

        if max_row != col:
            U[col], U[max_row] = U[max_row], U[col]
            # Swap L entries for columns already processed
            for k in range(col):
                L[col][k], L[max_row][k] = L[max_row][k], L[col][k]

        if abs(U[col][col]) < 1e-12:
            continue

        for row in range(col + 1, n):
            factor = U[row][col] / U[col][col]
            L[row][col] = factor
            for j in range(col, n):
                U[row][j] -= factor * U[col][j]

    return L, U

Explanation

Initialize L as the identity matrix and U as a copy of the input matrix.
For each column, find the row with the largest absolute value for partial pivoting to improve numerical stability.
Swap rows in both U and the already-computed portion of L.
Compute the multiplier for each row below the pivot and subtract the scaled pivot row from it. Store multipliers in L.
After processing all columns, A = L * U (with row permutations applied).

Complexity

Time: O(n^3)
Space: O(n^2) for L and U matrices

← #332 #334 →