While Singular Value Decomposition (SVD) offers a powerful way to analyze any matrix, LU decomposition provides a different, efficient factorization specifically tailored for square matrices and solving linear systems. It's a workhorse algorithm often employed behind the scenes in numerical software.

LU decomposition factors a square matrix $A$ into the product of two matrices: a lower triangular matrix $L$ and an upper triangular matrix $U$ .

$A = LU$

Here:

$L$ is a lower triangular matrix. This means all entries above the main diagonal are zero. Often, $L$ is constructed to have 1s on its main diagonal (this specific form is called Doolittle decomposition). $L = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 \\ l_{21} & 1 & 0 & \dots & 0 \\ l_{31} & l_{32} & 1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & l_{n3} & \dots & 1 \end{bmatrix}$
$U$ is an upper triangular matrix. This means all entries below the main diagonal are zero. $U = \begin{bmatrix} u_{11} & u_{12} & u_{13} & \dots & u_{1n} \\ 0 & u_{22} & u_{23} & \dots & u_{2n} \\ 0 & 0 & u_{33} & \dots & u_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & u_{nn} \end{bmatrix}$

Why Decompose A into L and U?

The primary advantage of LU decomposition lies in solving systems of linear equations of the form $Ax = b$ . If we can write $A = LU$ , then the original equation becomes:

$LUx = b$

Solving this equation might seem like adding complexity, but it actually simplifies the process by breaking it into two easier steps:

Solve $Ly = b$ for $y$ : Define an intermediate vector $y = Ux$ . Since $L$ is lower triangular, we can solve for the elements of $y$ sequentially starting from $y_1$ down to $y_n$ using a method called forward substitution. This is computationally inexpensive.
Solve $Ux = y$ for $x$ : Now that we have $y$ , we solve for the original unknown vector $x$ . Since $U$ is upper triangular, we can solve for the elements of $x$ sequentially starting from $x_n$ up to $x_1$ using backward substitution. This is also computationally inexpensive.

The LU decomposition transforms the problem of solving $Ax=b$ into two sequential solves involving triangular matrices, which are much faster using forward and backward substitution.

This two-step process is particularly efficient if you need to solve $Ax = b$ for the same matrix $A$ but multiple different vectors $b$ . You perform the computationally heavier LU decomposition of $A$ just once, and then repeat the fast forward and backward substitution steps for each new $b$ .

Existence and Variations

LU decomposition doesn't exist for every square matrix. For example, if a zero appears in a pivot position during the factorization process (which is related to Gaussian elimination), the basic algorithm fails.

In practice, LU decomposition is often performed with pivoting. This involves swapping rows of the matrix $A$ during the decomposition to ensure numerical stability and avoid division by small or zero pivots. Row swapping can be represented by a permutation matrix $P$ . The resulting decomposition is often written as:

$PA = LU$

where $P$ is a matrix that reorders the rows of $A$ . Solving $Ax=b$ then becomes solving $PAx = Pb$ , which proceeds similarly: solve $Ly = Pb$ using forward substitution, then $Ux = y$ using backward substitution. Most numerical libraries implement LU decomposition with pivoting.

While less directly visible in high-level machine learning model design compared to SVD or eigen-decomposition, LU decomposition is a fundamental tool in numerical linear algebra. It's frequently used under the hood in libraries (like NumPy and SciPy) when solving the linear systems that arise in various optimization problems and intermediate calculations within ML algorithms, such as solving the normal equations in Linear Regression. Understanding its principle helps appreciate the efficiency of these underlying solvers.