LU Decomposition

LU Decomposition is a fundamental technique in linear algebra that factors a matrix into the product of a lower triangular matrix and an upper triangular matrix. This decomposition is particularly useful in numerical analysis due to its numerical stability and efficiency, especially when combined with partial pivoting. LU Decomposition is widely employed in solving linear systems, inverting matrices, and computing determinants. In this article, we will explore the mathematical foundation of LU Decomposition, how to compute it, and its applications, providing practical examples to illustrate its use.

Understanding LU Decomposition

What is LU Decomposition?

LU Decomposition (or LU factorization) is the process of decomposing a given square matrix $\mathbf{A}$ into the product of two matrices:

$$\mathbf{A} = \mathbf{L} \mathbf{U}$$

Where:

• $\mathbf{L}$ is a lower triangular matrix (i.e., all elements above the main diagonal are zero).
• $\mathbf{U}$ is an upper triangular matrix (i.e., all elements below the main diagonal are zero).

Properties:

• Lower Triangular Matrix ($\mathbf{L}$):

  $$\mathbf{L} = \begin{pmatrix} l_{11} & 0 & 0 & \dots & 0 \\ l_{21} & l_{22} & 0 & \dots & 0 \\ l_{31} & l_{32} & l_{33} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & 0 \\ l_{n1} & l_{n2} & l_{n3} & \dots & l_{nn} \end{pmatrix}$$

• Upper Triangular Matrix ($\mathbf{U}$):

  $$\mathbf{U} = \begin{pmatrix} 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 & u_{nn} \\ 0 & 0 & 0 & \dots & u_{nn} \end{pmatrix}$$

Why Use LU Decomposition?

LU Decomposition is particularly useful in numerical analysis for several reasons:

1. Solving Linear Systems: LU Decomposition provides an efficient and numerically stable method for solving systems of linear equations, especially when dealing with multiple right-hand sides.

2. Matrix Inversion: LU Decomposition can be used to compute the inverse of a matrix by solving multiple linear systems.

3. Determinant Calculation: With $\mathbf{P}\mathbf{A}=\mathbf{L}\mathbf{U}$ and unit diagonal $\mathbf{L}$, $\det(\mathbf{A})=\det(\mathbf{P})\prod_i u_{ii}$; the permutation sign must not be omitted.

4. Numerical Stability: LU Decomposition, especially when combined with partial pivoting, provides a numerically stable method for solving linear systems, reducing the risk of computational errors.

5. Reusability: Once a matrix is decomposed into $\mathbf{L}$ and $\mathbf{U}$, it can be reused to solve multiple linear systems with different right-hand sides efficiently.

Mathematical Foundation of LU Decomposition

Lower and Upper Triangular Matrices

To understand LU Decomposition, it's important to understand the properties of lower and upper triangular matrices:

• Lower Triangular Matrix ($\mathbf{L}$): A matrix is lower triangular if all elements above the main diagonal are zero:

  $$\mathbf{L} = \begin{pmatrix} l_{11} & 0 & 0 & \dots & 0 \\ l_{21} & l_{22} & 0 & \dots & 0 \\ l_{31} & l_{32} & l_{33} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & 0 \\ l_{n1} & l_{n2} & l_{n3} & \dots & l_{nn} \end{pmatrix}$$

  Properties:

  • Simplifies matrix operations like solving linear systems via forward substitution.
  • The determinant of a lower triangular matrix is the product of its diagonal elements.

• Upper Triangular Matrix ($\mathbf{U}$): A matrix is upper triangular if all elements below the main diagonal are zero:

  $$\mathbf{U} = \begin{pmatrix} 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 & u_{nn} \\ 0 & 0 & 0 & \dots & u_{nn} \end{pmatrix}$$

  Properties:

  • Simplifies matrix operations like solving linear systems via back substitution.
  • The determinant of an upper triangular matrix is the product of its diagonal elements.

Existence of LU Decomposition

An LU factorization without row exchanges exists with nonzero elimination pivots when the required leading principal minors are nonzero. Nonsingularity alone is not sufficient to avoid row exchanges. When row exchanges are needed, use LU Decomposition with Partial Pivoting ($\mathbf{P}\mathbf{A}=\mathbf{L}\mathbf{U}$), where $\mathbf{P}$ records the swaps.

LU Decomposition with Partial Pivoting:

$$\mathbf{PA} = \mathbf{LU}$$

Where:

• $\mathbf{P}$ is a permutation matrix.
• $\mathbf{L}$ is a lower triangular matrix.
• $\mathbf{U}$ is an upper triangular matrix.

Partial pivoting prevents zero pivots for a nonsingular square matrix and is backward stable for most practical matrices, but it is not an unconditional stability guarantee: rare matrices exhibit large element growth. Singular matrices do not have a unique solve even if elimination returns factors.

Computing LU Decomposition

There are several methods for computing LU Decomposition, with Doolittle’s Method and LU Decomposition with Partial Pivoting being among the most commonly used.

1. Doolittle’s Method (Without Pivoting)

Doolittle's method is a popular algorithm for computing LU Decomposition, where the diagonal elements of $\mathbf{L}$ are set to 1, and $\mathbf{U}$ is computed directly.

Step-by-Step Process

Given a matrix $\mathbf{A}$, the steps to compute its LU Decomposition using Doolittle’s method are as follows:

1. Initialize $\mathbf{L}$ and $\mathbf{U}$:

   • Set $\mathbf{L}$ as an identity matrix (diagonal elements are 1, and all other elements are 0).
   • Set $\mathbf{U}$ as a zero matrix.

2. Compute Elements of $\mathbf{U}$:

   • For each element $u_{ij}$ in $\mathbf{U}$, compute it using the formula:

     $$u_{ij} = a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}, \qquad j\ge i$$

3. Compute Elements of $\mathbf{L}$:

   • For each element $l_{ij}$ in $\mathbf{L}$ (where $i > j$), compute it using the formula:

     $$l_{ij} = \frac{a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj}}{u_{jj}}$$

4. Repeat Until All Elements Are Computed:

   • Continue this process for all rows and columns until all elements of $\mathbf{L}$ and $\mathbf{U}$ are filled.

Example: LU Decomposition Using Doolittle’s Method

Problem Setup:

Consider the matrix $\mathbf{A}$:

$$\mathbf{A} = \begin{pmatrix} 2 & 3 & 1 \\ 4 & 7 & -1 \\ -2 & 4 & 5 \end{pmatrix}$$

We want to find its LU Decomposition using Doolittle’s method.

Step 1: Initialize $\mathbf{L}$ and $\mathbf{U}$

Start with:

$$\mathbf{L} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad \mathbf{U} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}$$

Step 2: Compute $\mathbf{U}$ and $\mathbf{L}$

1. Compute $u_{11}$, $u_{12}$, and $u_{13}$:

   $$u_{11} = 2$$ $$u_{12} = 3$$ $$u_{13} = 1$$

2. Compute $u_{22}$ and $u_{23}$:

   $$u_{22} = 7 - l_{21}u_{12} = 7 - 2 \times 3 = 7 - 6 = 1$$ $$u_{23} = -1 - l_{21}u_{13} = -1 - 2 \times 1 = -1 - 2 = -3$$

3. Compute $u_{33}$:

   $$u_{33} = 5 - l_{31}u_{13} - l_{32}u_{23}$$

   First, compute $l_{31}$ and $l_{32}$:

   $$l_{21} = 2 \quad (\text{from Step 2})$$ $$l_{31} = \frac{-2}{u_{11}} = \frac{-2}{2} = -1$$ $$l_{32} = \frac{4 - l_{31}u_{12}}{u_{22}} = \frac{4 - (-1) \times 3}{1} = \frac{4 + 3}{1} = 7$$

   Now,

   $$u_{33} = 5 - (-1) \times 1 - 7 \times (-3) = 5 + 1 + 21 = 27$$

4. Update $\mathbf{U}$:

   $$\mathbf{U} = \begin{pmatrix} 2 & 3 & 1 \\ 0 & 1 & -3 \\ 0 & 0 & 27 \end{pmatrix}$$

5. Update $\mathbf{L}$:

   $$\mathbf{L} = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & 7 & 1 \end{pmatrix}$$

Verification:

Compute $\mathbf{L} \mathbf{U}$ to verify that it equals $\mathbf{A}$:

$$\mathbf{L} \mathbf{U} = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & 7 & 1 \end{pmatrix} \begin{pmatrix} 2 & 3 & 1 \\ 0 & 1 & -3 \\ 0 & 0 & 27 \end{pmatrix} = \begin{pmatrix} 2 & 3 & 1 \\ 4 & 7 & -1 \\ -2 + 0 + 0 & -3 + 7 \times 1 + 0 & -1 \times 1 + 7 \times (-3) + 1 \times 27 \end{pmatrix} = \begin{pmatrix} 2 & 3 & 1 \\ 4 & 7 & -1 \\ -2 & 4 & 5 \end{pmatrix}$$

Thus, $\mathbf{A} = \mathbf{L} \mathbf{U}$ is verified.

2. LU Decomposition with Partial Pivoting (PA = LU)

When a matrix requires row exchanges to ensure numerical stability or when certain elements of $\mathbf{U}$ become zero or very small, LU Decomposition with partial pivoting is used. In this method, a permutation matrix $\mathbf{P}$ is introduced such that:

$$\mathbf{PA} = \mathbf{LU}$$

Where:

• $\mathbf{P}$ is a permutation matrix representing the row swaps.
• $\mathbf{L}$ is a lower triangular matrix.
• $\mathbf{U}$ is an upper triangular matrix.

Advantages of Partial Pivoting:

• Numerical Stability: Prevents division by very small numbers, reducing computational errors.
• Applicability: Avoids a zero elimination pivot for a nonsingular matrix. A singular system still needs rank detection and may have no unique solution.

Example: LU Decomposition with Partial Pivoting

Problem Setup:

Consider the matrix $\mathbf{A}$:

$$\mathbf{A} = \begin{pmatrix} 0 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix}$$

Attempting LU Decomposition without pivoting fails because $a_{11} = 0$, leading to division by zero. Using partial pivoting:

1. Identify the Pivot:

   • In the first column, the largest absolute value is 7 (in the third row). Swap the first and third rows.

2. Permutation Matrix $\mathbf{P}$: Partial pivoting also swaps the remaining second and third rows at the next step, giving

   $$\mathbf{P} = \begin{pmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix}$$

3. Compute $\mathbf{LU}$:

   • Perform LU Decomposition on the permuted matrix $\mathbf{PA}$.

4. Result:

   $$\mathbf{PA} = \mathbf{LU} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \frac{4}{7} & \frac{3}{14} & 1 \end{pmatrix} \begin{pmatrix} 7 & 8 & 9 \\ 0 & 2 & 3 \\ 0 & 0 & \frac{3}{14} \end{pmatrix}$$

   Multiplying the displayed factors reproduces $\mathbf{P}\mathbf{A}$ exactly.

Applications of LU Decomposition

1. Solving Linear Systems

Given a linear system of equations $\mathbf{A}\mathbf{x} = \mathbf{b}$, LU Decomposition can be used to solve for $\mathbf{x}$ efficiently. By decomposing $\mathbf{A}$ into $\mathbf{L}\mathbf{U}$, the system becomes:

$$\mathbf{L}\mathbf{U}\mathbf{x} = \mathbf{b}$$

This can be solved in two steps:

1. Solve $\mathbf{L}\mathbf{y} = \mathbf{b}$ for $\mathbf{y}$ using forward substitution.
2. Solve $\mathbf{U}\mathbf{x} = \mathbf{y}$ for $\mathbf{x}$ using back substitution.

Advantages:

• Efficiency: Especially beneficial when solving multiple systems with the same $\mathbf{A}$ but different $\mathbf{b}$ vectors.
• Numerical Stability: Provides a stable method for solving linear systems compared to direct inversion.

Example: Solving a Linear System Using LU Decomposition

Problem Setup:

Solve the system:

$$\begin{pmatrix} 2 & 3 & 1 \\ 4 & 7 & -1 \\ -2 & 4 & 5 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 3 \\ 5 \\ 1 \end{pmatrix}$$

Solution:

1. Perform LU Decomposition of $\mathbf{A}$:

   Using the corrected example:

   $$\mathbf{A} = \mathbf{L} \mathbf{U} = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & 7 & 1 \end{pmatrix} \begin{pmatrix} 2 & 3 & 1 \\ 0 & 1 & -3 \\ 0 & 0 & 27 \end{pmatrix}$$

2. Solve $\mathbf{L}\mathbf{y} = \mathbf{b}$:

   $$\begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & 7 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} 3 \\ 5 \\ 1 \end{pmatrix}$$

   Solving:

   • $y_1 = 3$
   • $2y_1 + y_2 = 5 \Rightarrow y_2 = 5 - 2(3) = -1$
   • $-y_1 + 7y_2 + y_3 = 1 \Rightarrow -3 + 7(-1) + y_3 = 1 \Rightarrow y_3 = 1 + 3 + 7 = 11$

   Thus,

   $$\mathbf{y} = \begin{pmatrix} 3 \\ -1 \\ 11 \end{pmatrix}$$

3. Solve $\mathbf{U}\mathbf{x} = \mathbf{y}$:

   $$\begin{pmatrix} 2 & 3 & 1 \\ 0 & 1 & -3 \\ 0 & 0 & 27 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 3 \\ -1 \\ 11 \end{pmatrix}$$

   Solving:

   • $27x_3 = 11 \Rightarrow x_3 = \frac{11}{27} \approx 0.4074$
   • $x_2 - 3x_3 = -1 \Rightarrow x_2 = -1 + 3(0.4074) \approx -1 + 1.2222 = 0.2222$
   • $2x_1 + 3x_2 + x_3 = 3 \Rightarrow 2x_1 + 3(0.2222) + 0.4074 = 3 \Rightarrow 2x_1 + 0.6666 + 0.4074 = 3 \Rightarrow 2x_1 = 3 - 1.074 = 1.926 \Rightarrow x_1 \approx 0.963$

   Solution Vector:

   $$\mathbf{x} \approx \begin{pmatrix} 0.963 \\ 0.222 \\ 0.407 \end{pmatrix}$$

2. Matrix Inversion

LU Decomposition can be used to compute the inverse of a matrix $\mathbf{A}$. By decomposing $\mathbf{A}$ into $\mathbf{L}\mathbf{U}$, the inverse can be found by solving the linear systems $\mathbf{A}\mathbf{X} = \mathbf{I}$, where $\mathbf{I}$ is the identity matrix. This involves solving multiple linear systems for each column of $\mathbf{X}$.

Advantages:

• Efficiency: Reusing the LU Decomposition for multiple systems simplifies the inversion process.
• Numerical practice: Reuses the factorization, but explicit inversion should be avoided when the actual goal is to solve $\mathbf{A}\mathbf{x}=\mathbf{b}$; triangular solves are cheaper and usually more accurate.

Example: Computing the Inverse Using LU Decomposition

Problem Setup:

Find the inverse of matrix $\mathbf{A}$ using its LU Decomposition.

$$\mathbf{A} = \begin{pmatrix} 2 & 3 & 1 \\ 4 & 7 & -1 \\ -2 & 4 & 5 \end{pmatrix} = \mathbf{L} \mathbf{U}$$

Solution:

1. Perform LU Decomposition:

   Using the corrected example:

   $$\mathbf{A} = \mathbf{L} \mathbf{U} = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & 7 & 1 \end{pmatrix} \begin{pmatrix} 2 & 3 & 1 \\ 0 & 1 & -3 \\ 0 & 0 & 27 \end{pmatrix}$$

2. Compute $\mathbf{A}^{-1}$:

   Solve $\mathbf{A}\mathbf{X} = \mathbf{I}$:

   • For each column $\mathbf{e}_i$ of the identity matrix $\mathbf{I}$, solve $\mathbf{A}\mathbf{x}_i = \mathbf{e}_i$ using LU Decomposition.

   • This involves solving $\mathbf{L}\mathbf{y}_i = \mathbf{e}_i$ and then $\mathbf{U}\mathbf{x}_i = \mathbf{y}_i$ for each $\mathbf{x}_i$.

3. Assemble the Inverse Matrix:

   Collect all solution vectors $\mathbf{x}_i$ to form $\mathbf{A}^{-1}$.

Result:

$$\mathbf{A}^{-1} = \begin{pmatrix} \frac{13}{18} & -\frac{11}{54} & -\frac{5}{27} \\ -\frac{1}{3} & \frac{2}{9} & \frac{1}{9} \\ \frac{5}{9} & -\frac{7}{27} & \frac{1}{27} \end{pmatrix}$$

Multiplication by $\mathbf{A}$ gives the identity; rounding only affects a decimal presentation, not the exact inverse.

3. Determinant Calculation

The determinant of a matrix $\mathbf{A}$ can be easily computed using LU Decomposition. Since $\mathbf{A} = \mathbf{L}\mathbf{U}$, and the determinant of a triangular matrix is the product of its diagonal elements, we have:

$$\text{det}(\mathbf{A}) = \text{det}(\mathbf{L}) \times \text{det}(\mathbf{U})$$

Given that the determinant of a lower triangular matrix $\mathbf{L}$ with 1s on the diagonal is 1, the determinant of $\mathbf{A}$ is simply the product of the diagonal elements of $\mathbf{U}$.

Example: Determinant Calculation Using LU Decomposition

Given:

$$\mathbf{U} = \begin{pmatrix} 2 & 3 & 1 \\ 0 & 1 & -3 \\ 0 & 0 & 27 \end{pmatrix}$$

Compute the determinant:

$$\text{det}(\mathbf{A}) = (2) \times (1) \times (27) = 54$$

Conclusion

LU Decomposition is a fundamental technique in linear algebra with numerous applications in solving linear systems, inverting matrices, and computing determinants. By understanding both the mathematical foundation and practical computation methods (such as Doolittle’s method and LU Decomposition with partial pivoting), you can leverage LU Decomposition to solve complex problems in data science, engineering, and beyond.

Whether you are dealing with small matrices or large-scale problems, mastering LU Decomposition will equip you with a powerful technique for analyzing and manipulating matrices, leading to more efficient and stable solutions in your work.