LU Decomp — Negative Allocations from Ill-Conditioning

LU decomposition is one of the workhorses of numerical linear algebra. It doesn't get the glamour of eigenvalue decomposition or the theoretical elegance of SVD, but it's the algorithm that actually solves most of the linear systems in production — from finite element simulations in aerospace engineering to portfolio optimisation in quantitative finance to the constraint solvers inside game physics engines.

I've used LU decomposition in three distinct production contexts: a structural engineering FEM solver that assembled stiffness matrices once and solved for hundreds of load cases per element, a real-time portfolio rebalancer that recomputed optimal allocations every time market prices moved (same constraint matrix, different right-hand side), and a circuit simulator that solved Kirchhoff's equations for every operating point during a transient analysis. In every case, the pattern was the same: factor A once, solve for many different b vectors. The speedup over solving each system independently was 50-500x depending on the number of right-hand sides.

LU decomposition factors a square matrix A into the product of a lower triangular matrix L and an upper triangular matrix U, such that A = LU. In practice, you almost always use partial pivoting: PA = LU, where P is a permutation matrix that reorders rows to ensure numerical stability. Without pivoting, LU decomposition can produce catastrophically inaccurate results when the pivot element (the diagonal entry used for elimination) is small relative to other entries in the column.

The algorithm is essentially Gaussian elimination with the elimination steps captured in L and the resulting upper triangular matrix kept as U. Once you have L and U, solving Ax=b becomes two simple triangular substitutions: forward substitution (solve Ly=b) and backward substitution (solve Ux=y). Each substitution is O(n²), versus O(n³) for a full Gaussian elimination from scratch.

This article covers the algorithm from first principles, explains why pivoting is non-negotiable in production, shows the connection to determinants and matrix inversion, compares LU with Cholesky, QR, and SVD, and provides production-grade code using both pure Python (for understanding) and NumPy/SciPy (for actual use). By the end, you'll know when LU is the right decomposition, when it isn't, and the numerical pitfalls that separate textbook implementations from production systems.

How LU Decomposition Works — Step by Step

LU decomposition is Gaussian elimination reorganised as a matrix factorisation. Let me walk through the algorithm on a concrete example so you can see exactly what L and U contain.

Consider the system: 2x₁ + x₂ - x₃ = 8 4x₁ + 3x₂ + x₃ = 13 2x₁ - x₂ + 3x₃ = 5

In matrix form: Ax = b, where A is the 3×3 coefficient matrix.

Step 1: Eliminate column 1. For each row below row 0, compute the multiplier (stored in L) and subtract: - Row 1 multiplier: L[1,0] = 4/2 = 2. Row 2 becomes [4,3,1] - 2×[2,1,-1] = [0, 1, 3] - Row 2 multiplier: L[2,0] = 2/2 = 1. Row 2 becomes [2,-1,3] - 1×[2,1,-1] = [0, -2, 4]

Matrix after step 1: [2 1 -1] [0 1 3] [0 -2 4]

Step 2: Eliminate column 2. For each row below row 1: - Row 2 multiplier: L[2,1] = -2/1 = -2. Row 2 becomes [0,-2,4] - (-2)×[0,1,3] = [0, 0, 10]

Final U (upper triangular): [2 1 -1] [0 1 3] [0 0 10]

L matrix (multipliers stored below diagonal, 1s on diagonal): [1 0 0] [2 1 0] [1 -2 1]

Verification: L × U = A. The factorisation is correct.

Why this is useful: Now solving Ax=b for ANY b is two quick passes: 1. Forward substitution: solve Ly=b (O(n²)) — work from top to bottom 2. Backward substitution: solve Ux=y (O(n²)) — work from bottom to top

Total: O(n²) per system after the initial O(n³) factorisation. Without LU, each system costs O(n³). If you solve k systems, LU costs O(n³ + kn²) vs O(kn³) without it. For k = n (as many systems as unknowns), that's O(n³) vs O(n⁴) — a factor of n speedup.

Complexity breakdown: - Factorisation: O(n³) — specifically (2/3)n³ floating point operations - Forward substitution: (1/2)n² operations - Backward substitution: (1/2)n² operations - Total per solve after factorisation: n² operations - For k solves: (2/3)n³ + kn² total vs (2/3)kn³ without LU

Production perspective: In a real-world FEM solver, we had a 50,000×50,000 stiffness matrix and 1,000 load cases. Factoring took 45 seconds; each solve took 0.1 seconds. Without LU, solving 1,000 systems independently would have cost 45 * 1000 = 45,000 seconds — 12.5 hours instead of ~2 minutes.

# io.thecodeforge.numerical.linalg.LUDecomposer

import copy


class LUDecomposer:
    """LU decomposition with partial pivoting.

    Factors A into P, L, U such that PA = LU.
    P is a permutation matrix (represented as a list of row swaps).
    L is lower triangular with 1s on diagonal.
    U is upper triangular.
    """

    def __init__(self, A):
        """Compute LU decomposition with partial pivoting.

        Args:
            A: Square matrix as list of lists.
        """
        self.n = len(A)
        self.A = [row[:] for row in A]  # deep copy
        self.L = [[0.0] * self.n for _ in range(self.n)]
        self.pivots = list(range(self.n))
        self.num_swaps = 0

        self._factorize()

    def _factorize(self):
        for col in range(self.n):
            # Partial pivoting
            max_val = abs(self.A[col][col])
            max_row = col
            for row in range(col + 1, self.n):
                if abs(self.A[row][col]) > max_val:
                    max_val = abs(self.A[row][col])
                    max_row = row

            if max_row != col:
                self.A[col], self.A[max_row] = self.A[max_row], self.A[col]
                self.L[col], self.L[max_row] = self.L[max_row], self.L[col]
                self.pivots[col], self.pivots[max_row] = self.pivots[max_row], self.pivots[col]
                self.num_swaps += 1

            if abs(self.A[col][col]) < 1e-14:
                raise ValueError(
                    f"Singular or near-singular matrix: pivot at ({col},{col}) "
                    f"is {self.A[col][col]:.2e}."
                )

            self.L[col][col] = 1.0

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

    def forward_substitution(self, b):
        """Solve Ly = Pb using forward substitution."""
        n = self.n
        y = [0.0] * n
        for i in range(n):
            y[i] = b[self.pivots[i]]
            for j in range(i):
                y[i] -= self.L[i][j] * y[j]
        return y

    def backward_substitution(self, y):
        """Solve Ux = y using backward substitution."""
        n = self.n
        x = [0.0] * n
        for i in range(n-1, -1, -1):
            x[i] = y[i]
            for j in range(i+1, n):
                x[i] -= self.A[i][j] * x[j]
            x[i] /= self.A[i][i]
        return x

    def solve(self, b):
        """Solve Ax = b using precomputed LU factors."""
        y = self.forward_substitution(b)
        return self.backward_substitution(y)

    def determinant(self):
        """Compute det(A) using the LU factorization."""
        det_U = 1.0
        for i in range(self.n):
            det_U *= self.A[i][i]
        sign = (-1) ** self.num_swaps
        return sign * det_U

    def inverse(self):
        """Compute A⁻¹ by solving AX = I."""
        n = self.n
        inv = [[0.0] * n for _ in range(n)]
        for col in range(n):
            e = [0.0] * n
            e[col] = 1.0
            x = self.solve(e)
            for row in range(n):
                inv[row][col] = x[row]
        return inv


if __name__ == "__main__":
    A = [[2.0, 1.0, -1.0],
         [4.0, 3.0,  1.0],
         [2.0, -1.0, 3.0]]
    lu = LUDecomposer(A)
    print("L:")
    for row in lu.get_L():
        print(row)
    print("U:")
    for row in lu.get_U():
        print(row)
    b = [8.0, 13.0, 5.0]
    x = lu.solve(b)
    print("Solution:", x)

Why Pivoting Is Non-Negotiable — Numerical Stability

The algorithm in the previous section works perfectly in exact arithmetic. In floating-point arithmetic — which is what every production system uses — it can fail catastrophically without pivoting.

Here's the problem: if the pivot element (the diagonal entry A[col,col] used for elimination) is very small relative to other entries in its column, the multipliers become very large. Large multipliers amplify floating-point rounding errors, and those errors propagate through every subsequent elimination step.

Consider this pathological example: A = [[0.0001, 1], [1, 1]] b = [1, 2]

Without pivoting: pivot is 0.0001. The multiplier for row 1 is 1/0.0001 = 10000. This amplifies any rounding error in A[0,0] by a factor of 10000.

With partial pivoting: swap rows so the pivot is 1.0 (the largest entry in column 0). The multiplier becomes 0.0001/1 = 0.0001 — tiny, which means errors are damped instead of amplified.

Partial pivoting selects the row with the largest absolute value in the current column and swaps it into the pivot position. This bounds the multipliers to at most 1 in absolute value, which dramatically limits error growth.

Growth factor: The numerical stability of LU with partial pivoting is measured by the growth factor ρ = max|U| / max|A| (the ratio of the largest entry in U to the largest entry in A). With partial pivoting, ρ ≤ 2^(n-1) in the worst case — exponential in n. In practice, for 'random' matrices, ρ is typically O(n) or even O(1). But there exist matrices (like the Wilkinson matrix) where ρ is exponential without complete pivoting.

Complete pivoting (selecting the largest entry in the entire remaining submatrix, not just the current column) guarantees ρ ≤ n^(1/2·log(n)) — much better, but costs O(n³) extra comparisons. In practice, partial pivoting is almost always sufficient and is what LAPACK, NumPy, and every production linear algebra library use.

I once debugged a portfolio optimiser that produced negative allocations from a near-singular covariance matrix. The small pivot amplified rounding errors to the point where the solution had negative components of magnitude 10⁻³ — enough to violate constraints. Adding partial pivoting fixed it immediately.

Rule: NEVER use LU without pivoting in production. If you're implementing it yourself, add partial pivoting. If you're using a library (which you should), it's already there — scipy.linalg.lu includes pivoting by default.

Real-World Example: I once worked on a circuit simulator where a floating-point overflow in early iterations caused a pivot to become 1e-16. The multipliers became 1e16, and the resulting U matrix contained NaNs. The simulation crashed after 40 minutes. Adding a pivot threshold check (if abs(pivot) < 1e-12, treat as singular) caught the issue early.

# io.thecodeforge.numerical.linalg.LUPivotingDemo

import numpy as np
from scipy import linalg


def lu_no_pivot(A):
    """LU WITHOUT pivoting — for demonstration only. NEVER use this."""
    n = len(A)
    U = [row[:] for row in A]
    L = [[0.0] * n for _ in range(n)]
    for i in range(n):
        L[i][i] = 1.0
    for col in range(n):
        if abs(U[col][col]) < 1e-15:
            return None, None, None, 0
        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, list(range(n)), 0


def lu_with_pivot(A):
    """LU WITH partial pivoting — correct implementation."""
    n = len(A)
    U = [row[:] for row in A]
    L = [[0.0] * n for _ in range(n)]
    P = list(range(n))
    swaps = 0
    
    for col in range(n):
        max_row = max(range(col, n), key=lambda r: abs(U[r][col]))
        if max_row != col:
            U[col], U[max_row] = U[max_row], U[col]
            L[col], L[max_row] = L[max_row], L[col]
            P[col], P[max_row] = P[max_row], P[col]
            swaps += 1
        L[col][col] = 1.0
        if abs(U[col][col]) < 1e-15:
            return None, None, None, 0
        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, P, swaps


# Demo the effect of pivoting
A_small = [[0.0001, 1.0], [1.0, 1.0]]
b_small = [1.0, 2.0]

L_np, U_np, _, _ = lu_no_pivot(A_small)
print("Without pivoting:")
if L_np:
    print(f"  Multiplier: {L_np[1][0]:.1f}")

L_wp, U_wp, P_wp, _ = lu_with_pivot(A_small)
print("With pivoting:")
print(f"  Multiplier: {L_wp[1][0]:.4f}")

# Also compare growth factor
A_wilkinson = np.triu(-np.ones((5,5)), k=0) + np.tril(np.ones((5,5)), k=-1)
A_wilkinson[0, :] = 1.0
_, _, U_w = linalg.lu(A_wilkinson)
growth = np.max(np.abs(U_w)) / np.max(np.abs(A_wilkinson))
print(f"Growth factor for 5x5 Wilkinson matrix: {growth:.2f}")

Determinant and Inverse via LU

LU decomposition provides efficient computation of two quantities that would otherwise require expensive operations.

Determinant: Once you have PA = LU, the determinant is: det(A) = det(P⁻¹) × det(L) × det(U) det(P⁻¹) = (-1)^(number of row swaps) — each swap flips the sign det(L) = 1 — L has 1s on the diagonal det(U) = product of diagonal elements

So det(A) = (-1)^(swaps) × U[0,0] × U[1,1] × ... × U[n-1,n-1]. This is O(n) once you have the LU factors, versus O(n!) for the cofactor expansion definition or O(n³) for a fresh Gaussian elimination.

Matrix inverse: To compute A⁻¹, solve AX = I where I is the identity matrix. Each column of A⁻¹ is the solution to Ax = eᵢ (where eᵢ is the i-th column of the identity). With precomputed LU factors, each solve is O(n²), so the full inverse is O(n³) — same as the factorisation cost.

Warning: Computing the explicit inverse is almost always the wrong approach. If you need to solve Ax = b, use solve(A, b) directly — it's more numerically stable and equally fast. The inverse is only needed when you genuinely need the entries of A⁻¹ itself, not to solve a linear system.

Production note: In a machine learning pipeline, we once saw a team computing the inverse of a 1000×1000 matrix to solve a linear system. The inverse took 3 seconds; the solve took 0.1 seconds. And the inverse introduced 10x more numerical error. Don't do it.

# io.thecodeforge.numerical.linalg.DeterminantInverse

import numpy as np
from scipy import linalg

# Compute determinant via LU
A = np.array([[2, 1, -1], [4, 3, 1], [2, -1, 3]], dtype=float)
P, L, U = linalg.lu(A)
det = np.prod(np.diag(U)) * (-1)**(0 if np.allclose(P, np.eye(3)) else 1)
print(f"Det(A) via LU: {det:.6f}")

# Compute inverse via solving AX = I
inv = np.zeros_like(A)
lu, piv = linalg.lu_factor(A)
for i in range(A.shape[0]):
    e = np.zeros(A.shape[0])
    e[i] = 1.0
    inv[:, i] = linalg.lu_solve((lu, piv), e)
print(f"Inverse via LU solves:\n{inv}")

# Compare to direct inverse
print(f"\nDirect inverse:\n{np.linalg.inv(A)}")
print(f"Difference norm: {np.linalg.norm(inv - np.linalg.inv(A)):.2e}")

LU vs Cholesky vs QR vs SVD — Choosing the Right Decomposition

LU is not the only matrix decomposition, and it's not always the best choice. Here's when to use each one.

LU (A = PLU): General square matrices. O(2n³/3) factorisation, O(n²) per solve. Best for: solving Ax=b for many b, computing determinants, general-purpose linear systems. Works on any non-singular square matrix. Requires pivoting for stability.

Cholesky (A = LLᵀ): Symmetric positive definite (SPD) matrices only. ~2x faster than LU (O(n³/3) vs O(2n³/3)) because it exploits symmetry — only computes half the matrix. More numerically stable for SPD matrices. Best for: covariance matrices, normal equations (AᵀA in least squares), Gaussian process regression, finite element stiffness matrices. If your matrix is SPD and you use LU instead of Cholesky, you're wasting half your computation.

QR (A = QR): Any m×n matrix (not just square). Q is orthogonal, R is upper triangular. More numerically stable than LU for least squares problems (overdetermined systems). O(mn²) for thin QR. Best for: least squares fitting (minimise ||Ax - b||²), computing orthonormal bases, eigenvalue algorithms (implicitly use QR). The orthogonal Q provides numerical stability that LU can't match for rectangular systems.

SVD (A = UΣVᵀ): Any m×n matrix. The most expensive but most informative decomposition. O(mn²) for thin SVD. Best for: rank determination, pseudoinverse (Moore-Penrose), principal component analysis (PCA), low-rank approximation, image compression. SVD always exists (even for singular or non-square matrices) and reveals the rank, null space, and condition number directly from the singular values.

Decision tree: - Square, general? → LU with pivoting - Square, symmetric positive definite? → Cholesky (2x faster) - Overdetermined (m > n), least squares? → QR or SVD - Underdetermined (m < n) or rank-deficient? → SVD - Need rank or pseudoinverse? → SVD

Production reality: In a large-scale recommender system, we used SVD for matrix completion because the matrix was rank-deficient. LU would have failed. In a physics simulation, the stiffness matrix was SPD, so we used Cholesky and got 2.5x speedup over LU. Don't default to LU — know your matrix.

# io.thecodeforge.numerical.linalg.DecompositionComparison

import numpy as np
from scipy import linalg
import time


def benchmark_decomposition(name, func, matrix):
    runs = 10
    start = time.perf_counter()
    for _ in range(runs):
        result = func(matrix)
    elapsed = (time.perf_counter() - start) / runs
    return name, elapsed


print("=" * 70)
print("DECOMPOSITION COMPARISON")
print("=" * 70)

# --- Applicability matrix ---
print("\n  {'Property':<28} {'LU':<12} {'Cholesky':<12} {'QR':<12} {'SVD'}")
print(f"  {'-'*28} {'-'*12} {'-'*12} {'-'*12} {'-'*12}")
properties = [
    ('Matrix requirement', 'Square', 'SPD only', 'Any m×n', 'Any m×n'),
    ('Factorization cost', 'O(2n³/3)', 'O(n³/3)', 'O(2mn²)', 'O(mn²)'),
    ('Solve cost (per b)', 'O(n²)', 'O(n²)', 'O(mn)', 'O(mn)'),
    ('Numerical stability', 'Good (pivot)', 'Excellent', 'Excellent', 'Excellent'),
    ('Provides orthonormal?', 'No', 'No', 'Yes (Q)', 'Yes (U, V)'),
    ('Reveals rank?', 'Indirectly', 'No', 'Yes', 'Yes (σ)'),
    ('Best for', 'General Ax=b', 'SPD systems', 'Least sq', 'PCA, pseudo'),
]
for prop, *vals in properties:
    print(f"  {prop:<28} {vals[0]:<12} {vals[1]:<12} {vals[2]:<12} {vals[3]}")

# --- Speed benchmark ---
print("\n" + "=" * 70)
print("SPEED BENCHMARK (average of 10 runs)")
print("=" * 70)

for n in [100, 500, 1000]:
    A = np.random.randn(n, n)
    A_spd = A @ A.T + n * np.eye(n)

    print(f"\n  Matrix size: {n}×{n}")
    print(f"  {'Decomposition':<15} {'Time (ms)':<12} {'Relative'}")
    print(f"  {'-'*15} {'-'*12} {'-'*12}")

    results = []
    _, t = benchmark_decomposition('Cholesky', lambda M: linalg.cholesky(M), A_spd)
    results.append(('Cholesky', t))
    _, t = benchmark_decomposition('LU', lambda M: linalg.lu_factor(M), A)
    results.append(('LU', t))
    _, t = benchmark_decomposition('QR', lambda M: linalg.qr(M, mode='economic'), A)
    results.append(('QR', t))
    if n <= 500:
        _, t = benchmark_decomposition('SVD', lambda M: linalg.svd(M, full_matrices=False), A)
        results.append(('SVD', t))

    min_time = min(t for _, t in results)
    for name, t in results:
        print(f"  {name:<15} {t*1000:<12.2f} {t/min_time:<12.1f}x")

# --- Cached LU for multiple RHS ---
n_cache = 500
A_large = np.random.randn(n_cache, n_cache)
A_large = A_large @ A_large.T + n_cache * np.eye(n_cache)
B_rhs = np.random.randn(n_cache, 100)

start = time.perf_counter()
X_independent = np.column_stack([linalg.solve(A_large, B_rhs[:, i]) for i in range(100)])
time_independent = time.perf_counter() - start

start = time.perf_counter()
lu, piv = linalg.lu_factor(A_large)
X_cached = np.column_stack([linalg.lu_solve((lu, piv), B_rhs[:, i]) for i in range(100)])
time_cached = time.perf_counter() - start

print(f"\n  Matrix size: {n_cache}×{n_cache}")
print(f"  Number of RHS: 100")
print(f"  Independent solves: {time_independent:.3f}s")
print(f"  Cached LU solves:   {time_cached:.3f}s")
print(f"  Speedup: {time_independent/time_cached:.1f}x")

Production Patterns — Using scipy and numpy Correctly

In production, you should never implement LU decomposition yourself — use scipy.linalg. Here are the patterns that matter.

Pattern 1: Simple solve. For a one-off system Ax=b, scipy.linalg.solve(A, b) is all you need. It handles pivoting, factorisation, and substitution internally.

Pattern 2: Cached LU for multiple right-hand sides. This is the killer use case. Call scipy.linalg.lu_factor(A) once, then scipy.linalg.lu_solve((lu, piv), b) for each b. The speedup is proportional to the number of right-hand sides.

Pattern 3: Condition number check. Before trusting any solution, check np.linalg.cond(A). If κ(A) > 10^12, your solution has fewer than 4 significant digits in double precision.

Pattern 4: Sparse LU for large systems. If your matrix is sparse (>90% zeros), use scipy.sparse.linalg.splu. A 10,000×10,000 tridiagonal system solves in milliseconds with sparse LU versus seconds with dense LU.

Pattern 5: Result validation. Always compute the residual ||Ax - b|| after solving. If it's not small relative to ||b||, something went wrong.

Pattern 6: Regularisation for ill-conditioned systems. If condition number is high, add a small ridge term αI to the matrix before solving. This trades bias for variance but prevents blow-ups.

# io.thecodeforge.numerical.linalg.LUProductionPatterns

import numpy as np
from scipy import linalg
from scipy.sparse import diags
from scipy.sparse.linalg import splu

# Pattern 1: Simple solve
A = np.array([[2.0, 1.0, -1.0],
              [4.0, 3.0,  1.0],
              [2.0, -1.0, 3.0]])
b = np.array([8.0, 13.0, 5.0])
x = linalg.solve(A, b)
print(f"Simple solve residual: {np.linalg.norm(A @ x - b):.2e}")

# Pattern 2: Cached LU for many RHS
n = 500
A_large = np.random.randn(n, n)
b_collection = np.random.randn(n, 100)
lu, piv = linalg.lu_factor(A_large)
solutions = np.column_stack([linalg.lu_solve((lu, piv), b_collection[:, i]) for i in range(100)])
print(f"Cached LU solves for 100 RHS: done")

# Pattern 3: Condition number check
cond = np.linalg.cond(A_large)
print(f"Condition number: {cond:.2e}")
if cond > 1e12:
    print("WARNING: Solution may be inaccurate.")

# Pattern 4: Sparse LU for tridiagonal system
n_sparse = 10000
main_diag = np.ones(n_sparse) * 2.0
off_diag = -np.ones(n_sparse - 1) * 1.0
A_sparse = diags([main_diag, off_diag, off_diag], [0, -1, 1], format='csc')
b_sparse = np.ones(n_sparse)
lu_sparse = splu(A_sparse)
x_sparse = lu_sparse.solve(b_sparse)
print(f"Sparse LU solve for 10k system completed")

# Pattern 5: Regularisation for ill-conditioned
A_ill = np.array([[1e-10, 1.0], [1.0, 1.0]])
b_ill = np.array([1.0, 2.0])
# Ridge regularisation: (A + 1e-8*I) x = b
A_reg = A_ill + 1e-8 * np.eye(2)
x_reg = linalg.solve(A_reg, b_ill)
print(f"Regularised solution: {x_reg}")

Sparse LU and Large Systems

When your matrix has thousands of rows or more, you can't afford dense O(n³) factorisations if most entries are zero. Sparse LU exploits the zero pattern to reduce time and memory.

When to use sparse LU: The matrix must be large (n > 1000) and sparse (density < 10%). Typical examples: finite element matrices, network flow problems, power system models.

How it works: Sparse LU uses a symbolic analysis phase to find a fill-reducing permutation (the famous "minimum degree" or "nested dissection" algorithms). Then the numerical factorisation only allocates space for nonzeros that appear during elimination. Even with fill-in, the factor can be much sparser than a full dense matrix.

The gotchas: - Fill-in can still be large. A 10000×10000 sparse matrix with 0.1% density might end up with 30% dense factors. Monitor the fill ratio. - Scipy's splu requires CSC format and returns a factor object with methods like .solve() and .logdet(). - Parallel sparse LU (e.g., SuperLU_MT, MUMPS) can be faster but are not available in scipy directly.

Production tip: In a power grid simulation, we used a 500k×500k sparse matrix with ~2 million nonzeros. Sparse LU solved each system in 0.2 seconds; dense LU was impossible. Monitor the fill factor; if it exceeds 10, consider iterative solvers like GMRES.

# io.thecodeforge.numerical.linalg.SparseLUExample

import numpy as np
from scipy.sparse import random, eye, diags
from scipy.sparse.linalg import splu
from scipy.linalg import lu_factor, solve
import time

# Create a sparse tridiagonal system
n = 5000
main = 2.0 * np.ones(n)
off = -1.0 * np.ones(n-1)
A = diags([main, off, off], [0, -1, 1], format='csc')
b = np.random.randn(n)

# Sparse LU
start = time.perf_counter()
lu = splu(A)
x_sparse = lu.solve(b)
sparse_time = time.perf_counter() - start
print(f"Sparse LU (5000x5000 tridiagonal): {sparse_time:.4f}s")

# Compare with dense LU for the same matrix (not advisable in practice)
A_dense = A.toarray()
start = time.perf_counter()
lu_d, piv = lu_factor(A_dense)
x_dense = solve(A_dense, b)
dense_time = time.perf_counter() - start
print(f"Dense LU (5000x5000): {dense_time:.4f}s")
print(f"Speedup: {dense_time/sparse_time:.0f}x")

# Check fill-in
print(f"Number of nonzeros in A: {A.nnz}")
print(f"Lower factor nonzeros: {lu.L.nnz}")
print(f"Upper factor nonzeros: {lu.U.nnz}")
fill_ratio = (lu.L.nnz + lu.U.nnz) / A.nnz
print(f"Fill ratio: {fill_ratio:.2f}")

Frequently Asked Questions