LU Decomposition — Factoring Matrices for Efficient Solving

def lu_decompose(A):
    """LU decomposition without pivoting. A = L @ U."""
    n = len(A)
    L = [[1.0 if i==j else 0.0 for j in range(n)] for i in range(n)]
    U = [row[:] for row in A]
    for col in range(n):
        for row in range(col+1, n):
            if abs(U[col][col]) < 1e-12:
                raise ValueError('Singular matrix — use pivoting')
            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

def forward_sub(L, b):
    n = len(b)
    y = [0.0] * n
    for i in range(n):
        y[i] = b[i] - sum(L[i][j]*y[j] for j in range(i))
    return y

def back_sub(U, y):
    n = len(y)
    x = [0.0] * n
    for i in range(n-1, -1, -1):
        x[i] = (y[i] - sum(U[i][j]*x[j] for j in range(i+1,n))) / U[i][i]
    return x

A = [[2,1,-1],[4,3,1],[2,-1,3]]
L, U = lu_decompose(A)
b = [8,13,5]
y = forward_sub(L, b)   # Solve Ly = b
x = back_sub(U, y)       # Solve Ux = y
print(f'Solution: {[round(v,4) for v in x]}')

from functools import reduce
from operator import mul

def determinant_lu(A):
    """det(A) = product of U's diagonal elements."""
    _, U = lu_decompose(A)
    return reduce(mul, (U[i][i] for i in range(len(U))), 1.0)

import numpy as np
A = [[2,1,-1],[4,3,1],[2,-1,3]]
print(f'det via LU: {determinant_lu(A):.4f}')
print(f'numpy det:  {np.linalg.det(A):.4f}')