LU Decomposition

Introduction

A generic linear system looks like this:

whose short (matrix) form is

where is a 2D matrix with rows and columns whose values are known, are the unknowns and the r.h.s. is a -dimensional vector.

LU Decomposition

When solving systems of linear equations, one powerful method is to decompose the matrix A into a product of two matrices: a lower triangular matrix L and an upper triangular matrix U:

where:

  • L is the lower triangular matrix, meaning it has non-zero elements only on the diagonal and below.
  • U is the upper triangular matrix, with non-zero elements only on the diagonal and above.

This decomposition allows us to break the problem into simpler parts.

For example, for a 4x4 matrix A, the decomposition in Eq. (1) would look like this:

With this decomposition, we can solve the linear system:

First, we solve for y using forward substitution:

Then, we solve for x by backsubstitution:

[!NOTE]

Why Decompose?

The advantage of this approach is that triangular systems (like L and U) are easier to solve.

In a triangular system:

  • L (lower triangular) only has non-zero elements on the diagonal and below. This makes it possible to solve for the unknowns sequentially using a method called forward substitution, starting from the first row and moving downwards.
  • U (upper triangular) has non-zero elements on the diagonal and above. For U, you use a method called backsubstitution, which solves the system from the last row upwards.

These two methods —- forward substitution and back substitution -— are computationally efficient. Instead of having to deal with the full matrix at once, you solve the system in steps, reducing the overall complexity. This stepwise approach results in fewer calculations compared to methods like Gaussian elimination, making it faster and less prone to numerical errors for large systems.

Forward Substitution

The forward substitution for solving L is:

For the remaining elements of y:

Click here to see an example for a 4x4 case

Back Substitution

Once y is found, we can use backsubstitution to solve U for x:

For the remaining elements of x:

Click here to see an example for a 4x4 case

By decomposing the matrix A into L and U, we can solve the system of equations step-by-step, first with forward substitution and then with backsubstitution, making the process more efficient.

Performing the LU decomposition

In the previous section, we saw that once we have the lower and upper triangular matrices ( and , respectively) such that , we can easily solve a linear system like . Now, we see how to compute and , i.e., how to perform the LU decomposition.

Let's start writing:

and the associated system of equations is:

As you can see from this example, we have equations and unknowns. Since the number of unknowns is greater than the number of equations, we are invited to specify of the unknowns arbitrarily and then try to solve for the others. We therefore choose:

We now apply the following procedure:

Crout’s algorithm provides a very efficient way to solve the system of equations generated by the matrix decomposition . The method tackles the set of equations, which include both the elements of and , by cleverly rearranging them so that fewer unknowns are solved at each step. This reordering allows us to sequentially solve for each and , as outlined below:

  1. First, we set the diagonal elements of to unity:

  1. Then, for each , we follow these two steps:
  • Use the following equation to solve for each where (in the below equation, the summation term is zero when ):

  • Use the following equation to solve for where :

As you proceed through these iterations, you will notice that each and is determined by previously computed values, allowing for an “in-place” decomposition. This means that no extra storage is required for these values, as they overwrite the corresponding locations in the matrix.

The algorithm essentially fills in the matrix by columns from left to right and within each column, from top to bottom.

Thus, after completing the process, the resulting matrix looks like this:

This combined matrix stores both the ’s and ’s, effectively performing decomposition in a highly efficient manner.

Here an example of the LU decomposition

Pivoting

Pivoting is a technique used to improve the numerical stability of the factorization process: without it, LU decomposition can suffer from numerical instability, especially when small or zero elements are encountered on the diagonal.

Pivoting refers to reordering the rows (or sometimes columns) of a matrix during the decomposition to avoid dividing by small or zero pivot elements, which can lead to large numerical errors.

In the context of LU decomposition, partial pivoting (the most common form) means swapping rows to ensure that the largest absolute value in each column becomes the pivot element (diagonal element of the matrix). This helps reduce round-off errors and increases the numerical stability of the factorization.

Without pivoting, if a small pivot element is encountered during the factorization process, division by that small value can amplify rounding errors, making the solution inaccurate. Pivoting minimizes this risk by ensuring that the pivot element is large enough to avoid such issues.

Pivoting in the Crout Algorithm

In the Crout method, you decompose the matrix  into:

where has ones on the diagonal (in some variations), and is an upper triangular matrix.

Standard Crout Algorithm (Without Pivoting):

  1. Initialize as a lower triangular matrix with arbitrary entries and as an upper triangular matrix.
  2. For each column of the matrix , compute the elements of and by solving systems of equations based on the known entries of , , and .

However, this approach is prone to numerical issues if a small or zero pivot element appears during the decomposition.

Crout Algorithm with Partial Pivoting:

In Crout’s method with partial pivoting, the algorithm includes an additional step where the rows of the matrix are swapped to ensure that the pivot element (the diagonal element of ) is the largest element in its column. Here’s how the process works:

  1. Initialize the LU Decomposition: Start with matrix , and begin the process of decomposing it into .
  2. At each step (for each column):
    • Before computing the current column of and , search for the largest absolute value in the current column (from the diagonal element down) and swap rows so that the largest value becomes the pivot element. 3 Row Swapping:
    • When a row swap is performed, you must update both and to reflect the new row ordering. Typically, these row swaps are recorded in a permutation matrix , and the decomposition becomes:

 where is the permutation matrix that tracks row swaps. 4. Compute the LU Factors: - After pivoting, compute the entries of the and matrices as usual. Ensure that the pivot element (diagonal element of ) is now large enough to avoid division by small numbers.

Example of LU Decomposition with Pivoting:

Consider a matrix :

Without pivoting, you would start by dividing by 0, which is not possible. With partial pivoting:

  1. Find the largest element in the first column (7 in this case).
  2. Swap rows 1 and 3:

  1. Continue the decomposition as normal, ensuring that each diagonal element is large enough to avoid numerical issues.

Matrix Inverse with LU Decomposition

LU decomposition is useful for solving a series of problems with the same matrix and different matrices. This is advantageous for computing the inverse of , as only one decomposition is required.

The inverse of a matrix is defined such that:

where (I) is the identity matrix.

For a general 3x3 matrix, this looks like:

This can be set up as () -problems with different matrices, where the matrix becomes the -th column in the inverse matrix.

For example, the first problem is:

The second column of the inverse can be computed by changing to , and the third column with , and so on.

This method is efficient because only back- and forward-substitution is required after the initial LU decomposition.

[!WARNING] To compute the inverse of the matrix, you need to solve () -problems, one for each column of the inverse. If you're trying to compute the inverse just to solve an problem, you introduce more numerical error and slow down your computation time. In short, avoid computing the matrix inverse unless absolutely necessary, and instead solve the system of equations directly.

Applications

Here a list of examples involving the matrix inversion:

  1. Extraction of spectral densities from lattice correlators