Table of contents
1.
Introduction
2.
Doolittle Algorithm
2.1.
Example
2.2.
Code Example
3.
FAQs
4.
Key Takeaways
Last Updated: Mar 27, 2024

Doolittle Algorithm

Author Rajat Agrawal
0 upvote
Career growth poll
Do you think IIT Guwahati certified course can help you in your career?

Introduction

The LU decomposition is the factorization of a square matrix into two triangular matrices, one upper triangular matrix, and one lower triangular matrix, with the product of these two matrices equaling the original matrix.

LU decomposition has a variety of applications, including solving a system of equations, which is an integral part of many other applications, such as finding current in a circuit, finding the inverse of a matrix, and finding the determinant of a matrix.

Doolittle’s Algorithm provides a way to factor a matrix A into an LU decomposition.

Doolittle Algorithm

A square matrix can always be factored into a lower triangular matrix and an upper triangular matrix. In other words, [A] = [L][U].

Doolittle's method allows us to factor A into a LU decomposition without going through the complexity of Gaussian Elimination.

We assume that a LU decomposition exists for a general n×n matrix A and write the form of L and U directly. We then systematically use the equations from the multiplications required for A = LU to solve the entries in L and U.

A = LU can be expanded as follow:

<math xmlns="http://www.w3.org/1998/Math/MathML"><mfenced open="[" close="]"><mtable><mtr><mtd><msub><mi>A</mi><mn>00</mn></msub></mtd><mtd><msub><mi>A</mi><mn>01</mn></msub></mtd><mtd><msub><mi>A</mi><mn>02</mn></msub></mtd></mtr><mtr><mtd><msub><mi>A</mi><mn>10</mn></msub></mtd><mtd><msub><mi>A</mi><mn>11</mn></msub></mtd><mtd><msub><mi>A</mi><mn>12</mn></msub></mtd></mtr><mtr><mtd><msub><mi>A</mi><mn>20</mn></msub></mtd><mtd><msub><mi>A</mi><mn>21</mn></msub></mtd><mtd><msub><mi>A</mi><mn>22</mn></msub></mtd></mtr></mtable></mfenced><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><msub><mi>L</mi><mn>10</mn></msub></mtd><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><msub><mi>L</mi><mn>20</mn></msub></mtd><mtd><msub><mi>L</mi><mn>21</mn></msub></mtd><mtd><mn>1</mn></mtd></mtr></mtable></mfenced><mo>&#xA0;</mo><mfenced open="[" close="]"><mtable><mtr><mtd><msub><mi>U</mi><mn>00</mn></msub></mtd><mtd><msub><mi>U</mi><mn>01</mn></msub></mtd><mtd><msub><mi>U</mi><mn>02</mn></msub></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><msub><mi>U</mi><mn>11</mn></msub></mtd><mtd><msub><mi>U</mi><mn>12</mn></msub></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><msub><mi>U</mi><mn>22</mn></msub></mtd></mtr></mtable></mfenced><mo>&#xA0;</mo><mspace linebreak="newline"/></math>

The terms of the U matrix is given by:

<math xmlns="http://www.w3.org/1998/Math/MathML"><mo>&#xA0;</mo><mo>&#xA0;</mo><mo>&#xA0;</mo><mo>&#xA0;</mo><mo>&#xA0;</mo><mo>&#xA0;</mo><mo>&#x2200;</mo><mi>j</mi><mo>&#xA0;</mo><mspace linebreak="newline"/><mi>i</mi><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mn>0</mn><mo>&#xA0;</mo><mo>&#x2192;</mo><msub><mi>U</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><msub><mi>A</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>&#xA0;</mo><mspace linebreak="newline"/><mi>i</mi><mo>&#xA0;</mo><mo>&gt;</mo><mo>&#xA0;</mo><mn>0</mn><mo>&#xA0;</mo><mo>&#x2192;</mo><mo>&#xA0;</mo><msub><mi>U</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><msub><mi>A</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>&#xA0;</mo><mo>-</mo><mo>&#xA0;</mo><mstyle displaystyle="false"><munderover><mo>&#x2211;</mo><mrow><mi>k</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>i</mi><mo>-</mo><mn>1</mn></mrow></munderover></mstyle><mo>&#xA0;</mo><msub><mi>L</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub><msub><mi>U</mi><mrow><mi>k</mi><mi>j</mi></mrow></msub></math>

The terms of the L matrix is given by:

<math xmlns="http://www.w3.org/1998/Math/MathML"><mo>&#xA0;</mo><mo>&#xA0;</mo><mo>&#xA0;</mo><mo>&#xA0;</mo><mo>&#xA0;</mo><mo>&#xA0;</mo><mo>&#x2200;</mo><mi>i</mi><mspace linebreak="newline"/><mi>j</mi><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mn>0</mn><mo>&#xA0;</mo><mo>&#x2192;</mo><mo>&#xA0;</mo><msub><mi>L</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mfrac><msub><mi>A</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mi>U</mi><mrow><mi>j</mi><mi>j</mi></mrow></msub></mfrac><mspace linebreak="newline"/><mi>j</mi><mo>&#xA0;</mo><mo>&gt;</mo><mo>&#xA0;</mo><mn>0</mn><mo>&#xA0;</mo><mo>&#x2192;</mo><mo>&#xA0;</mo><msub><mi>L</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mfrac><msub><mi>A</mi><mrow><mi>i</mi><mi>j</mi><mo>&#xA0;</mo><mo>-</mo><mo>&#xA0;</mo><munderover><mo>&#x2211;</mo><mrow><mi>k</mi><mo>=</mo><mn>0</mn><mo>&#xA0;</mo></mrow><mrow><mi>j</mi><mo>-</mo><mn>1</mn></mrow></munderover><mo>&#xA0;</mo><msub><mi>L</mi><mrow><mi>i</mi><mi>k</mi><mo>&#xA0;</mo></mrow></msub><msub><mi>U</mi><mrow><mi>k</mi><mi>j</mi></mrow></msub></mrow></msub><msub><mi>U</mi><mrow><mi>j</mi><mi>j</mi></mrow></msub></mfrac></math>

Example

Given a matrix A, decompose it to Lower and Upper triangular matrix using the Doolittle algorithm.

<math xmlns="http://www.w3.org/1998/Math/MathML"><mi>A</mi><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>2</mn></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mn>2</mn></mtd><mtd><mn>2</mn></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mn>2</mn></mtd><mtd><mn>3</mn></mtd></mtr></mtable></mfenced><mo>&#xA0;</mo></math>

Solution: 

Using Doolittle’s algorithm, we can write A = LU, where L is the lower triangular matrix and U is the upper triangular matrix.

<math xmlns="http://www.w3.org/1998/Math/MathML"><mfenced open="[" close="]"><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mn>2</mn></mtd><mtd><mn>2</mn></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mn>2</mn></mtd><mtd><mn>3</mn></mtd></mtr></mtable></mfenced><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><msub><mi>l</mi><mn>10</mn></msub></mtd><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><msub><mi>l</mi><mn>20</mn></msub></mtd><mtd><msub><mi>l</mi><mn>21</mn></msub></mtd><mtd><mn>1</mn></mtd></mtr></mtable></mfenced><mo>&#xA0;</mo><mfenced open="[" close="]"><mtable><mtr><mtd><msub><mi>u</mi><mn>00</mn></msub></mtd><mtd><msub><mi>u</mi><mn>01</mn></msub></mtd><mtd><msub><mi>u</mi><mn>02</mn></msub></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><msub><mi>u</mi><mn>11</mn></msub></mtd><mtd><msub><mi>u</mi><mn>12</mn></msub></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><msub><mi>u</mi><mn>22</mn></msub></mtd></mtr></mtable></mfenced><mspace linebreak="newline"/></math>

Now, multiply the L and U matrix.

<math xmlns="http://www.w3.org/1998/Math/MathML"><mfenced open="[" close="]"><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><msub><mi>l</mi><mn>10</mn></msub></mtd><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><msub><mi>l</mi><mn>20</mn></msub></mtd><mtd><msub><mi>l</mi><mn>21</mn></msub></mtd><mtd><mn>1</mn></mtd></mtr></mtable></mfenced><mo>&#xA0;</mo><mfenced open="[" close="]"><mtable><mtr><mtd><msub><mi>u</mi><mn>00</mn></msub></mtd><mtd><msub><mi>u</mi><mn>01</mn></msub></mtd><mtd><msub><mi>u</mi><mn>02</mn></msub></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><msub><mi>u</mi><mn>11</mn></msub></mtd><mtd><msub><mi>u</mi><mn>12</mn></msub></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><msub><mi>u</mi><mn>22</mn></msub></mtd></mtr></mtable></mfenced><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mfenced open="[" close="]"><mtable><mtr><mtd><msub><mi>u</mi><mn>00</mn></msub></mtd><mtd><msub><mi>u</mi><mn>01</mn></msub></mtd><mtd><msub><mi>u</mi><mn>02</mn></msub></mtd></mtr><mtr><mtd><msub><mi>l</mi><mn>10</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mrow><mn>00</mn><mo>&#xA0;</mo></mrow></msub></mtd><mtd><msub><mi>l</mi><mn>10</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>01</mn></msub><mo>+</mo><msub><mi>u</mi><mn>11</mn></msub></mtd><mtd><mo>&#xA0;</mo><msub><mi>l</mi><mn>10</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>02</mn></msub><mo>+</mo><msub><mi>u</mi><mn>12</mn></msub></mtd></mtr><mtr><mtd><msub><mi>l</mi><mn>20</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>00</mn></msub></mtd><mtd><mo>&#xA0;</mo><mo>&#xA0;</mo><msub><mi>l</mi><mn>20</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>01</mn></msub><mo>+</mo><msub><mi>l</mi><mn>21</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>11</mn></msub><mo>&#xA0;</mo></mtd><mtd><mo>&#xA0;</mo><msub><mi>l</mi><mn>20</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>02</mn></msub><mo>+</mo><msub><mi>l</mi><mn>21</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>12</mn></msub><mo>+</mo><msub><mi>u</mi><mn>22</mn></msub></mtd></mtr></mtable></mfenced><mspace linebreak="newline"/></math>

Now, compare the multiplied LU matrix with the A matrix row by row.

<math xmlns="http://www.w3.org/1998/Math/MathML"><mfenced open="[" close="]"><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mn>2</mn></mtd><mtd><mn>2</mn></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mn>2</mn></mtd><mtd><mn>3</mn></mtd></mtr></mtable></mfenced><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mfenced open="[" close="]"><mtable><mtr><mtd><msub><mi>u</mi><mn>00</mn></msub></mtd><mtd><msub><mi>u</mi><mn>01</mn></msub></mtd><mtd><msub><mi>u</mi><mn>02</mn></msub></mtd></mtr><mtr><mtd><msub><mi>l</mi><mn>10</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>00</mn></msub></mtd><mtd><msub><mi>l</mi><mn>10</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>01</mn></msub><mo>+</mo><msub><mi>u</mi><mn>11</mn></msub></mtd><mtd><msub><mi>l</mi><mn>10</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>02</mn></msub><mo>+</mo><msub><mi>u</mi><mn>12</mn></msub></mtd></mtr><mtr><mtd><msub><mi>l</mi><mn>20</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>00</mn></msub></mtd><mtd><mo>&#xA0;</mo><mo>&#xA0;</mo><msub><mi>l</mi><mn>20</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>01</mn></msub><mo>+</mo><msub><mi>l</mi><mn>21</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>11</mn></msub><mo>&#xA0;</mo></mtd><mtd><msub><mi>l</mi><mn>20</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>02</mn></msub><mo>+</mo><msub><mi>l</mi><mn>21</mn></msub><mo>&#xB7;</mo><msub><mi>u</mi><mn>12</mn></msub><mo>+</mo><msub><mi>u</mi><mn>22</mn></msub></mtd></mtr></mtable></mfenced><mspace linebreak="newline"/></math>

After comparing each row element in both the above matrices, we will get the values of unknowns as follows:

<math xmlns="http://www.w3.org/1998/Math/MathML"><mspace linebreak="newline"/><mo>&#xA0;</mo><msub><mi>u</mi><mn>00</mn></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mn>1</mn><mo>,</mo><mo>&#xA0;</mo><msub><mi>u</mi><mn>01</mn></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mn>1</mn><mo>,</mo><mo>&#xA0;</mo><msub><mi>u</mi><mn>02</mn></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mn>1</mn><mo>,</mo><mo>&#xA0;</mo><msub><mi>u</mi><mn>11</mn></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mn>1</mn><mo>,</mo><mo>&#xA0;</mo><msub><mi>u</mi><mn>12</mn></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mn>1</mn><mo>,</mo><mo>&#xA0;</mo><msub><mi>u</mi><mn>22</mn></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mn>1</mn><mspace linebreak="newline"/><mo>&#xA0;</mo><msub><mi>l</mi><mn>10</mn></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mn>1</mn><mo>,</mo><mo>&#xA0;</mo><msub><mi>l</mi><mn>20</mn></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mn>1</mn><mo>,</mo><mo>&#xA0;</mo><msub><mi>l</mi><mn>21</mn></msub><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mn>1</mn></math>

So finally, the L and U matrix becomes:

<math xmlns="http://www.w3.org/1998/Math/MathML"><mi>L</mi><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd></mtr></mtable></mfenced><mspace linebreak="newline"/><mspace linebreak="newline"/><mi>U</mi><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd></mtr></mtable></mfenced></math>

Code Example

The C++ program to calculate the Upper and Lower Triangular Matrix is given below:

#include <bits/stdc++.h>
using namespace std;

const int MAX = 100;

void LUDecomposition(int mat[][MAX], int n)
{
    int lower[n][n], upper[n][n];
    memset(lower, 0, sizeof(lower));
    memset(upper, 0, sizeof(upper));

    // Decomposing matrix into Upper and Lower triangular matrix
    for (int i = 0; i < n; i++)
    {
        // Upper Triangular
        for (int k = i; k < n; k++)
        {
            // Summation of L(i, j) * U(j, k)
            int sum = 0;
            for (int j = 0; j < i; j++)
                sum += (lower[i][j] * upper[j][k]);

            // Evaluating U(i, k)
            upper[i][k] = mat[i][k] - sum;
        }

        // Lower Triangular
        for (int k = i; k < n; k++)
        {
            if (i == k)
                lower[i][i] = 1; // Diagonal as 1
            else
            {
                // Summation of L(k, j) * U(j, i)
                int sum = 0;
                for (int j = 0; j < i; j++)
                    sum += (lower[k][j] * upper[j][i]);

                // Evaluating L(k, i)
                lower[k][i] = (mat[k][i] - sum) / upper[i][i];
            }
        }
    }

    // setw is for displaying in correct manner
    cout << setw(6)
        << "      Lower Triangular"
        << setw(32)
        << "   Upper Triangular" << endl;

    // Printing the result
    for (int i = 0; i < n; i++)
    {
        // Lower Triangular Matrix
        for (int j = 0; j < n; j++)
            cout << setw(6) << lower[i][j] << "\t";
        cout << "\t";

        // Upper Triangular Matrix
        for (int j = 0; j < n; j++)
            cout << setw(6) << upper[i][j] << "\t";
        cout << endl;
    }
}


int main()
{
    int mat[][MAX] = {{1, 1, 1}, {1, 2, 2}, {1, 2, 3}};

    LUDecomposition(mat, 3);
    return 0;
}
You can also try this code with Online C++ Compiler
Run Code

Output: 

The matrix we decomposed is:

<math xmlns="http://www.w3.org/1998/Math/MathML"><mi>A</mi><mo>&#xA0;</mo><mo>=</mo><mo>&#xA0;</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>2</mn></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mn>2</mn></mtd><mtd><mn>2</mn></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mn>2</mn></mtd><mtd><mn>3</mn></mtd></mtr></mtable></mfenced><mo>&#xA0;</mo></math>

 

FAQs

1. What is LU Decomposition?

The LU decomposition is the factorization of a square matrix into two triangular matrices, one upper triangular matrix, and one lower triangular matrix, with the product of these two matrices equaling the original matrix.

2. State the Doolittle Algorithm.

According to Doolittle Algorithm, a square matrix can always be factored into a lower triangular matrix and an upper triangular matrix. In other words, [A] = [L][U].

3. Is LU decomposition can be applied to any matrix?

No, the LU decomposition can only be applied on Square Matrices.

Key Takeaways

In this article, we have extensively discussed the LU Decomposition, Doolittle Algorithm, and how to decompose a square matrix using the Doolittle algorithm.

We hope that this blog has helped you enhance your knowledge regarding Linear Algebra. If you want to learn more, check out our articles on the Matrix and Types of Matrix.

Do upvote our blog to help other ninjas grow.

Happy Learning!

Live masterclass