LU Decomposition with Partial Pivoting
Factor a small matrix into permuted lower- and upper-triangular factors, then solve a system by forward and back substitution.
Problem
LU decomposition writes P A = L U, with P a row-permutation, L unit lower triangular, and U upper triangular. Partial pivoting swaps rows to put the largest available pivot on the diagonal, which keeps the factorization numerically stable. Solving Ax=b then costs only two triangular solves.
Factorization
Gaussian elimination with pivoting eliminates entries below each pivot, recording the multipliers in L and the row swaps in P. Once factored, any right-hand side is solved cheaply by reusing L and U.
import numpy as np
from scipy.linalg import lu, solve_triangular
A=np.array([[2.,1.,1.],[4.,3.,3.],[8.,7.,9.]])
b=np.array([4.,10.,26.])
P,L,U=lu(A)
y=solve_triangular(L,P.T@b,lower=True)
x=solve_triangular(U,y)
print('x',np.round(x,3))
print('check',np.allclose(A@x,b))Result
The solver factors A once, applies the permutation to b, solves L y = P'b downward, then U x = y upward. The recovered x satisfies Ax=b. Partial pivoting chose the row with the largest leading entry (8) as the first pivot, which prevents small pivots from amplifying rounding error. The same L and U solve any number of right-hand sides without re-factoring.
- Factoring is O(n^3) but each subsequent solve is only O(n^2), so LU pays off with many right-hand sides.
- Without pivoting, a zero or tiny pivot makes elimination fail or lose accuracy.
- Direct LU solvers handle the dense sub-blocks in Kronos coupled electromagnetic-structural models.