mathematics//linear algebra//matrix factorization//LU decomposition

The LU decomposition is the factorization of a square matrix into a lower triangular matrix \(L\) and an upper triangular matrix \(U\), usually after reordering its rows, and it is the engine under `np.linalg.solve` and every general-purpose linear solver: Gaussian elimination from school, recorded so that it can be reused. With the row permutation \(P\) written in,

The LU decomposition is the factorization of a square matrix into a lower triangular matrix LLL and an upper triangular matrix UUU, usually after reordering its rows, and it is the engine under np.linalg.solve and every general-purpose linear solver: Gaussian elimination from school, recorded so that it can be reused. With the row permutation PPP written in,

PA=LU.PA=LU .PA=LU.

Solving Ax=bAx=bAx=b then takes two substitutions. First Ly=PbLy=PbLy=Pb, top to bottom, then Ux=yUx=yUx=y, bottom to top, each one unknown at a time. The factorization costs about 23n3\tfrac23n^332​n3 operations; each substitution pair afterwards costs about n2n^2n2. For a 1,000 × 1,000 system that is roughly 670 million operations once and two million per new right-hand side, which is why a simulator or a filter that solves with the same matrix many times factors it once and keeps LLL and UUU.

This is why solve beats inv.

Forming A−1A^{-1}A−1 costs about three times as much as the factorization, and multiplying by it is less accurate than substituting, because the inverse carries rounding of its own into every product (linear system of equations).

The row reordering is what keeps it stable. Without pivoting, a tiny number in a pivot position would be divided into everything below it and wreck the result; partial pivoting swaps the largest available entry into place at each step, and in practice this is stable for almost every matrix met in engineering.

It assumes nothing about the matrix beyond being square and invertible, which is its strength and its waste. A symmetric positive definite matrix gets the same answer at half the cost through the Cholesky decomposition; a rectangular least-squares problem needs the QR decomposition instead.

The determinant falls out for free as the product of the diagonal of UUU (with a sign from the permutation), which is how libraries compute it. Large sparse systems, the meshes of a structural model or a power-flow Jacobian, use sparse LU variants that reorder the unknowns to keep the factors sparse, and that ordering matters more than the arithmetic.