LU Factorization Method

LU Factorization method is also known as Triangular method or LU Decomposition method. It uses the fact that if a square matrix, $A=[a_{ij}]$ has non zero principal minors i.e

$$a_{11} \neq 0,\ \begin{vmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{vmatrix} \neq 0,\ \begin{vmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{vmatrix} \neq 0, \text{ etc.}$$

then it can be factored into a multiple of a lower triangular matrix $L$ and an upper triangular matrix $U$ as $A=LU$. There are two main approaches for solving a linear system with LU factorization, they differ in the assumption of the diagonal elements for either $L$ or $U$ matrix. These methods are:

1. Doolittle Method
2. Crout Method

Doolittle Method

In the Doolittle method, we set the diagonal elements of the $L$ matrix to be $1$. On the other hand, the diagonal elements of the $U$ matrix may take any value. We initialize $L$ and $U$ as:

$$L=\begin{bmatrix} 1 & 0 & 0\\ l_{21} & 1 & 0\\ l_{31} & l_{32} & 1 \end{bmatrix} \text{ and } U=\begin{bmatrix} u_{11} & u_{12} & u_{13}\\ 0 & u_{22} & u_{23}\\ 0 & 0 & u_{33} \end{bmatrix}$$

Doolittle Method Working Procedure

We know a system of linear equations is represented in matrix notation as:

$$\begin{equation} AX=B \end{equation}$$

We suppose,

$$\begin{equation} A=LU \end{equation}$$

From $(1)$ and $(2)$, we can write,

$$\begin{equation} LUX=B \end{equation}$$

Again suppose,

$$\begin{equation} UX=Y \end{equation}$$ $$or, \begin{bmatrix} u_{11} & u_{12} & u_{13}\\ 0 & u_{22} & u_{23}\\ 0 & 0 & u_{33} \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix}= \begin{bmatrix} y_1\\ y_2 \\ y_3 \end{bmatrix}$$ $$\begin{equation} or, \begin{bmatrix} u_{11}x_1 + u_{12}x_2 + u_{13}x_3\\ u_{22}x_2 + u_{23}x_3\\ u_{33}x_3 \end{bmatrix}= \begin{bmatrix} y_1\\ y_2 \\ y_3 \end{bmatrix} \end{equation}$$

From $(3)$ and $(4)$, we can write,

$$LY=B$$ $$or, \begin{bmatrix} 1 & 0 & 0\\ l_{21} & 1 & 0\\ l_{31} & l_{32} & 1 \end{bmatrix} \begin{bmatrix} y_1\\ y_2 \\ y_3 \end{bmatrix}= \begin{bmatrix} b_1\\ b_2 \\ b_3 \end{bmatrix}$$ $$\begin{equation} or, \begin{bmatrix} y_1\\ l_{21}y_1+y_2\\ l_{31}y_1+l_{32}y_2+y_3 \end{bmatrix}=\begin{bmatrix} b_1\\ b_2 \\ b_3 \end{bmatrix} \end{equation}$$

Expanding $eq^n(2)$, we get,

$$\begin{bmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{bmatrix}= \begin{bmatrix} 1 & 0 & 0\\ l_{21} & 1 & 0\\ l_{31} & l_{32} & 1 \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13}\\ 0 & u_{22} & u_{23}\\ 0 & 0 & u_{33} \end{bmatrix}$$ $$or, \begin{bmatrix} u_{11} & u_{12} & u_{13}\\ l_{21}u_{11} & l_{21}u_{12}+u_{22} & l_{21}u_{13} + u_{23}\\ l_{31}u_{11} & l_{31}u_{12}+l_{32}u_{22} & l_{31}u_{13} + l_{32}u_{23} + u_{33} \end{bmatrix}=\begin{bmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{bmatrix}$$

Equating corresponding components we get,

1. $u_{11}=a_{11},\ u_{12}=a_{12},\ u_{13}=a_{13}$
2. $l_{21}u_{11}=a_{21}\rightarrow l_{21}=\frac{a_{21}}{u_{11}}$
3. $l_{31}u_{11}=a_{31}\rightarrow l_{31}=\frac{a_{31}}{u_{11}}$
4. $l_{21}u_{12}+u_{22}=a_{22}\rightarrow u_{22}=a_{22}-l_{21}u_{12}$
5. $l_{21}u_{13}+u_{23}=a_{23}\rightarrow u_{23}=a_{23}-l_{21}u_{13}$
6. $l_{31}u_{12}+l_{32}u_{22}=a_{32}\rightarrow l_{32}=\frac{a_{32}-l_{31}u_{12}}{u_{22}}$
7. $l_{31}u_{13}+l_{32}u_{23}+u_{33}=a_{33}\rightarrow u_{33}=a_{33}-l_{31}u_{13}-l_{32}u_{23}$

We can also find out the coefficients of $L$ and $U$ matrices directly with:

$$u_{ij} = a_{ij}-\sum_{k=1}^{i-1}l_{ik}u_{kj}\ \{i=1...N,j=1...N, j\geq i\}$$ $$l_{ij} = \frac{a_{ij}-\sum_{k=1}^{j-1}l_{ik}u_{kj}}{u_{jj}}\ \{i=1...N,j=1...N, j<i\}$$

Thus we know $L$ and $U$ matrices are known. Then we need to find out $Y$ by solving $LY=B$ and finally $X$ by solving $UX=Y$.

Doolittle Method Algorithm

1. Start
2. Read the size of the system
3. Read the elements of the coefficient matrix $A$
4. Read the elements of constant vector $B$
5. Calculate the matrices $L$ and $U$
6. Find $Y$ by solving $LY=B$ using forward substitution
7. Find $X$ by solving $UX=Y$ using backward substitution
8. Print $X$
9. Stop

Doolittle Method Program

C

/*
 * Doolittle LU factorization in C
 */

#include <stdio.h>

// Maximum system size
#define MS 3

void printMat(char *name, double m[][MS], int N)
{
  printf("%s:\n", name);
  for (int i = 0; i < N; i++)
  {
    for (int j = 0; j < N; j++)
      printf("%8.3f", m[i][j]);
    printf("\n");
  }
}

void printVec(char *name, double m[], int N)
{
  printf("%s:\n", name);
  for (int i = 0; i < N; i++)
    printf("%8.3f", m[i]);
  printf("\n");
}

int main()
{
  // The size of linear system
  int N = MS;

  // Preallocate enough memory for the system
  double A[MS][MS] = {0}, L[MS][MS] = {0}, U[MS][MS] = {0};
  double X[MS] = {0}, Y[MS] = {0}, B[MS] = {0};

  // Get the size of the system
  printf("Enter the size of linear system: ");
  scanf("%d", &N);

  // Limit the max system size
  if (N >= MS)
  {
    printf("Please enter size below %d", MS);
    return -1;
  }

  // Ask for the coefficients
  printf("\nEnter the coefficient matrix A:\n");
  for (int i = 0; i < N; i++)
  {
    for (int j = 0; j < N; j++)
    {
      printf("Enter A[%d][%d] = ", i + 1, j + 1);
      scanf("%lf", &A[i][j]);
    }
  }

  // Ask for the constants
  printf("\nEnter the constant vector B:\n");
  for (int i = 0; i < N; i++)
  {
    printf("Enter B[%d] = ", i + 1);
    scanf("%lf", &B[i]);
  }

  // Find out U and L matrix
  for (int i = 0; i < N; i++)
  {
    for (int j = 0; j < N; j++)
    {
      if (j >= i)
      {
        if (j == i)
          L[i][j] = 1;
        else
          L[i][j] = 0;

        U[i][j] = A[i][j];
        for (int k = 0; k < i; k++)
          U[i][j] -= L[i][k] * U[k][j];
      }
      else
      {
        L[i][j] = A[i][j];
        for (int k = 0; k < j; k++)
          L[i][j] -= L[i][k] * U[k][j];

        L[i][j] /= U[j][j];
        U[i][j] = 0;
      }
    }
  }

  printMat("L matrix is", L, N);
  printMat("U matrix is", U, N);

  // Find out Y by solving LY=B using forward substitution
  for (int i = 0; i < N; i++)
  {
    Y[i] = B[i];
    for (int j = 0; j < i; j++)
    {
      Y[i] -= L[i][j] * Y[j];
    }
  }
  printVec("Y is", Y, N);

  // Find out X by solving UX=Y using backward substitution
  for (int i = N - 1; i >= 0; i--)
  {
    X[i] = Y[i];
    for (int j = i + 1; j < N; j++)
    {
      X[i] -= X[j] * U[i][j];
    }
    X[i] /= U[i][i];
  }

  printVec("The solution X is", X, N);
  return 0;
}

Python

#
#  Doolittle LU factorization in Python
#

def print_mat(name, m):
  print(f'{name}:')
  for arr in m:
    print("\t".join([f"{el:.3f}" for el in arr]))

def print_vec(name, m):
  print(f'{name}:')
  print("\t".join([f"{el:.3f}" for el in m]))

## The size of linear system
N = int(input("Enter the size of linear system: "))
A = [[0]*N for _ in range(N)]
L = [[0]*N for _ in range(N)]
U = [[0]*N for _ in range(N)]
X = [0]*N
Y = [0]*N
B = [0]*N


# Ask for the coefficients
print("\nEnter the coefficient matrix A:");
for i in range(N):
  for j in range(N):
    A[i][j] = float(input(f"Enter A[{i+1}][{j+1}] = "))

# Ask for the constants
print("\nEnter the constant vector B:");
for i in range(N):
  B[i] = float(input(f"Enter B[{i+1}] = "))

# Find out U and L matrix
for i in range(N):
  for j in range(N):
    if j >= i:
      if j == i:
        L[i][j] = 1
      else:
        L[i][j] = 0
      U[i][j] = A[i][j]
      for k in range(i):
        U[i][j] -= L[i][k] * U[k][j]
    else:
      L[i][j] = A[i][j]
      for k in range(j):
        L[i][j] -= L[i][k] * U[k][j]
      L[i][j] /= U[j][j]
      U[i][j] = 0

print_mat("\nL", L)
print_mat("\nU", U)

# Find out Y by solving LY=B using forward substitution
for i in range(N):
  Y[i] = B[i]
  for j in range(i):
    Y[i] -= L[i][j] * Y[j]

print_vec("\nY", Y)

# Find out X by solving UX=Y using backward substitution
for i in reversed(range(N)):
  X[i] = Y[i]
  for j in range(i+1, N):
    X[i] -= X[j] * U[i][j]
  X[i] /= U[i][i]

print_vec("\nThe solution X is", X)

Crout Method

In this method we set the diagonal elements of $U$ matrix to be $1$. On the other hand, the diagonal elements of $L$ may assume any value. We initialize $L$ and $U$ as:

$$L=\begin{bmatrix} l_{11} & 0 & 0\\ l_{21} & l_{22} & 0\\ l_{31} & l_{32} & l_{33} \end{bmatrix} \text{ and } U=\begin{bmatrix} 1 & u_{12} & u_{13}\\ 0 & 1 & u_{23}\\ 0 & 0 & 1 \end{bmatrix}$$

Crout Method Working Procedure

We know a system of linear equations is represented in matrix notation as:

$$\begin{equation} AX=B \end{equation}$$

We suppose,

$$\begin{equation} A=LU \end{equation}$$

From $(1)$ and $(2)$, we can write,

$$\begin{equation} LUX=B \end{equation}$$

Again suppose,

$$\begin{equation} UX=Y \end{equation}$$ $$or, \begin{bmatrix} 1 & u_{12} & u_{13}\\ 0 & 1 & u_{23}\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix}= \begin{bmatrix} y_1\\ y_2 \\ y_3 \end{bmatrix}$$ $$\begin{equation} or, \begin{bmatrix} x_1 + u_{12}x_2 + u_{13}x_3\\ x_2 + u_{23}x_3\\ x_3 \end{bmatrix}= \begin{bmatrix} y_1\\ y_2 \\ y_3 \end{bmatrix} \end{equation}$$

From $(3)$ and $(4)$, we can write,

$$LY=B$$ $$or, \begin{bmatrix} l_{11} & 0 & 0\\ l_{21} & l_{22} & 0\\ l_{31} & l_{32} & l_{33} \end{bmatrix} \begin{bmatrix} y_1\\ y_2 \\ y_3 \end{bmatrix}= \begin{bmatrix} b_1\\ b_2 \\ b_3 \end{bmatrix}$$ $$\begin{equation} or, \begin{bmatrix} l_{11}y_1\\ l_{21}y_1+l_{22}y_2\\ l_{31}y_1+l_{32}y_2+l_{33}y_3 \end{bmatrix}=\begin{bmatrix} b_1\\ b_2 \\ b_3 \end{bmatrix} \end{equation}$$

Expanding $eq^n(2)$, we get,

$$\begin{bmatrix} l_{11} & 0 & 0\\ l_{21} & l_{22} & 0\\ l_{31} & l_{32} & l_{33} \end{bmatrix} \begin{bmatrix} 1 & u_{12} & u_{13}\\ 0 & 1 & u_{23}\\ 0 & 0 & 1 \end{bmatrix}= \begin{bmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{bmatrix}$$ $$or, \begin{bmatrix} l_{11} & l_{11}u_{12} & l_{11}u_{13}\\ l_{21} & l_{21}u_{12}+l_{22} & l_{21}u_{13} + l_{22}u_{23}\\ l_{31} & l_{31}u_{12}+l_{32} & l_{31}u_{13} + l_{32}u_{23} + l_{33} \end{bmatrix}=\begin{bmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{bmatrix}$$

Equating corresponding elements we get,

1. $l_{11}=a_{11},\ l_{21} = a_{21},\ l_{31}=a_{31}$
2. $l_{11}u_{12}=a_{12}\rightarrow u_{12}=\frac{a_{12}}{l_{11}}$
3. $l_{11}u_{13}=a_{13}\rightarrow u_{13}=\frac{a_{13}}{l_{11}}$
4. $l_{21}u_{12}+l_{22}=a_{22}\rightarrow l_{22} = a_{22}-l_{21}u_{12}$
5. $l_{31}u_{12}+l_{32}=a_{32}\rightarrow l_{32}=a_{32}-l_{31}u_{12}$
6. $l_{21}u_{13}+l_{22}u_{23}=a_{23}\rightarrow u_{23}=\frac{a_{23}-l_{21}u_{13}}{l_{22}}$
7. $l_{31}u_{13}+l_{32}u_{23}+l_{33}=a_{33}\rightarrow l_{33}=a_{33}-l_{31}u_{13}-l_{32}u_{23}$

We can directly find these coefficients for the Crout method with the following:

$$u_{ij} = \frac{a_{ij}-\sum_{k=1}^{i-1}l_{ik}u_{kj}}{l_{ii}}\ \{i=1...N,j=1...N, j> i\}$$ $$l_{ij} = a_{ij} - \sum_{k=1}^{j-1}l_{ik}u_{kj}\ \{i=1...N,j=1...N, j\leq i\}$$

Crout Method Algorithm

The algorithm is the same as that of the Doolittle method only the part for calculating $L$ and $U$ is different.

Crout Method Programs

C

/*
 * Crout LU factorization in C
 */

#include <stdio.h>

// Maximum system size
#define MS 20

void printMat(char *name, double m[][MS], int N)
{
  printf("%s:\n", name);
  for (int i = 0; i < N; i++)
  {
    for (int j = 0; j < N; j++)
      printf("%8.3f", m[i][j]);
    printf("\n");
  }
}

void printVec(char *name, double m[], int N)
{
  printf("%s:\n", name);
  for (int i = 0; i < N; i++)
    printf("%8.3f", m[i]);
  printf("\n");
}

int main()
{
  // The size of linear system
  int N = MS;

  // Preallocate enough memory for the system
  double A[MS][MS] = {0}, L[MS][MS] = {0}, U[MS][MS] = {0};
  double X[MS] = {0}, Y[MS] = {0}, B[MS] = {0};

  // Get the size of the system
  printf("Enter the size of linear system: ");
  scanf("%d", &N);

  // Limit the max system size
  if (N >= MS)
  {
    printf("Please enter size below %d", MS);
    return -1;
  }

  // Ask for the coefficients
  printf("\nEnter the coefficient matrix A:\n");
  for (int i = 0; i < N; i++)
  {
    for (int j = 0; j < N; j++)
    {
      printf("Enter A[%d][%d] = ", i + 1, j + 1);
      scanf("%lf", &A[i][j]);
    }
  }

  // Ask for the constants
  printf("\nEnter the constant vector B:\n");
  for (int i = 0; i < N; i++)
  {
    printf("Enter B[%d] = ", i + 1);
    scanf("%lf", &B[i]);
  }

  // Find out the L and U matrices
  for (int i = 0; i < N; i++)
  {
    for (int j = 0; j < N; j++)
    {
      if (j > i)
      {
        L[i][j] = 0;
        U[i][j] = A[i][j];
        for (int k = 0; k < i; k++)
          U[i][j] -= L[i][k] * U[k][j];
        U[i][j] /= L[i][i];
      }
      else
      {
        if (j == i)
          U[i][j] = 1;
        else
          U[i][j] = 0;

        L[i][j] = A[i][j];
        for (int k = 0; k < j; k++)
          L[i][j] -= L[i][k] * U[k][j];
      }
    }
  }
  printMat("L", L, N);
  printMat("U", U, N);

  // Find out Y by solving LY=B using forward substitution
  for (int i = 0; i < N; i++)
  {
    Y[i] = B[i];
    for (int j = 0; j < i; j++)
    {
      Y[i] -= L[i][j] * Y[j];
    }
    Y[i] /= L[i][i];
  }
  printVec("Y", Y, N);

  // Find out X by solving UX=Y using backward substitution
  for (int i = N - 1; i >= 0; i--)
  {
    X[i] = Y[i];
    for (int j = i + 1; j < N; j++)
    {
      X[i] -= X[j] * U[i][j];
    }
  }
  printVec("\nThe solution X is", X, N);
}

Python

#
#  Crout LU factorization in Python
#

def print_mat(name, m):
  print(f'{name}:')
  for arr in m:
    print("\t".join([f"{el:.3f}" for el in arr]))

def print_vec(name, m):
  print(f'{name}:')
  print("\t".join([f"{el:.3f}" for el in m]))

## The size of linear system
N = int(input("Enter the size of linear system: "))
A = [[0]*N for _ in range(N)]
L = [[0]*N for _ in range(N)]
U = [[0]*N for _ in range(N)]
X = [0]*N
Y = [0]*N
B = [0]*N


# Ask for the coefficients
print("\nEnter the coefficient matrix A:");
for i in range(N):
  for j in range(N):
    A[i][j] = float(input(f"Enter A[{i+1}][{j+1}] = "))

# Ask for the constants
print("\nEnter the constant vector B:");
for i in range(N):
  B[i] = float(input(f"Enter B[{i+1}] = "))

# Find out U and L matrix
for i in range(N):
  for j in range(N):
    if j > i:
      L[i][j] = 0
      U[i][j] = A[i][j]
      for k in range(i):
        U[i][j] -= L[i][k] * U[k][j]
      U[i][j] /= L[i][i]
    else:
      if j == i:
        U[i][j] = 1
      else:
        U[i][j] = 0
      
      L[i][j] = A[i][j]
      for k in range(j):
        L[i][j] -= L[i][k] * U[k][j]

print_mat("\nL", L)
print_mat("\nU", U)

# Find out Y by solving LY=B using forward substitution
for i in range(N):
  Y[i] = B[i]
  for j in range(i):
    Y[i] -= L[i][j] * Y[j]
  Y[i] /= L[i][i]

print_vec("\nY", Y)

# Find out X by solving UX=Y using backward substitution
for i in reversed(range(N)):
  X[i] = Y[i]
  for j in range(i+1, N):
    X[i] -= X[j] * U[i][j]

print_vec("\nThe solution X is", X)

Complexity

The time complexity of the LU Factorization method is $O(N^3)$. Similarly, its space complexity is $O(N^2)$. However, it uses three times the memory compared to Gauss-Jordan and Gauss-Elimination methods.

Advantages

1. Since the system is represented with extra matrices/vectors, the precision of calculation is doubled.

Disadvantages

1. Complicated to implement.
2. Uses 3 times the memory compared to Gauss-Elimination and Gauss-Jordan methods.