Doolittle Algorithm | LU Decomposition

Doolittle Algorithm: The Doolittle Algorithm is a method for performing LU Decomposition, where a given matrix is decomposed into a lower triangular matrix L and an upper triangular matrix U. This decomposition is widely used in solving systems of linear equations, inverting matrices, and computing determinants.

LU Decomposition

Let A be a square matrix. An LU factorization refers to the factorization of A, with proper row and/or column orderings or permutations, into two factors, a lower triangular matrix L and an upper triangular matrix U, A=LU.

Read More: LU Decomposition

Doolittle Algorithm

It is always possible to factor a square matrix into a lower triangular matrix and an upper triangular matrix. That is, [A] = [L][U]
Doolittle's method provides an alternative way to factor A into an LU decomposition without going through the hassle of Gaussian Elimination.
For a general n×n matrix A, we assume that an LU decomposition exists, and write the form of L and U explicitly. We then systematically solve for the entries in L and U from the equations that result from the multiplications necessary for A=LU.

Terms of the U matrix are given by:

\forall j\\ i=0 \rightarrow U_{ij} = A_{ij}\\ i>0 \rightarrow U_{ij} = A_{ij} - \sum_{k=0}^{i-1}L_{ik}U_{kj}            

And the terms for L matrix:

 \forall i\\ j=0 \rightarrow L_{ij} = \frac{A_{ij}}{U_{jj}}\\ j>0 \rightarrow L_{ij} = \frac{A_{ij} - \sum_{k=0}^{j-1}L_{ik}U_{kj}}{U_{jj}}         

Example:

Input : 

Output :

Implementation:

C++

// C++ Program to decompose a matrix into
// lower and upper triangular matrix
#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 nicely
    cout << setw(6) 
         << "      Lower Triangular" 
         << setw(32)
         << "Upper Triangular" << endl;

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

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

// Driver code
int main()
{
    int mat[][MAX]
        = { { 2, -1, -2 }, { -4, 6, 3 }, { -4, -2, 8 } };

    luDecomposition(mat, 3);
    return 0;
}

Java

// Java Program to decompose a matrix into
// lower and upper triangular matrix
class GFG {

    static int MAX = 100;
    static String s = "";
    static void luDecomposition(int[][] mat, int n)
    {
        int[][] lower = new int[n][n]; 
        int[][] upper = new int[n][n];

        // 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 nicely
        System.out.println(setw(2) + "     Lower Triangular"
                           + setw(10) + "Upper Triangular");

        // Displaying the result :
        for (int i = 0; i < n; i++) 
        {
            // Lower
            for (int j = 0; j < n; j++)
                System.out.print(setw(4) + lower[i][j]
                                 + "\t");
            System.out.print("\t");

            // Upper
            for (int j = 0; j < n; j++)
                System.out.print(setw(4) + upper[i][j]
                                 + "\t");
            System.out.print("\n");
        }
    }
    static String setw(int noOfSpace)
    {
        s = "";
        for (int i = 0; i < noOfSpace; i++)
            s += " ";
        return s;
    }

    // Driver code
    public static void main(String arr[])
    {
        int mat[][] = { { 2, -1, -2 },
                        { -4, 6, 3 },
                        { -4, -2, 8 } };

        luDecomposition(mat, 3);
    }
}

/* This code contributed by PrinciRaj1992 */

Python3

# Python3 Program to decompose
# a matrix into lower and
# upper triangular matrix
MAX = 100


def luDecomposition(mat, n):

    lower = [[0 for x in range(n)]
             for y in range(n)]
    upper = [[0 for x in range(n)]
             for y in range(n)]

    # Decomposing matrix into Upper
    # and Lower triangular matrix
    for i in range(n):

        # Upper Triangular
        for k in range(i, n):

            # Summation of L(i, j) * U(j, k)
            sum = 0
            for j in range(i):
                sum += (lower[i][j] * upper[j][k])

            # Evaluating U(i, k)
            upper[i][k] = mat[i][k] - sum

        # Lower Triangular
        for k in range(i, n):
            if (i == k):
                lower[i][i] = 1  # Diagonal as 1
            else:

                # Summation of L(k, j) * U(j, i)
                sum = 0
                for j in range(i):
                    sum += (lower[k][j] * upper[j][i])

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

    # setw is for displaying nicely
    print("Lower Triangular\t\tUpper Triangular")

    # Displaying the result :
    for i in range(n):

        # Lower
        for j in range(n):
            print(lower[i][j], end="\t")
        print("", end="\t")

        # Upper
        for j in range(n):
            print(upper[i][j], end="\t")
        print("")


# Driver code
mat = [[2, -1, -2],
       [-4, 6, 3],
       [-4, -2, 8]]

luDecomposition(mat, 3)

# This code is contributed by mits

C#

// C# Program to decompose a matrix into
// lower and upper triangular matrix
using System;

class GFG {

    static int MAX = 100;
    static String s = "";
    static void luDecomposition(int[, ] mat, int n)
    {
        int[, ] lower = new int[n, n];
        int[, ] upper = new int[n, n];

        // 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 nicely
        Console.WriteLine(setw(2) + "     Lower Triangular"
                          + setw(10) + "Upper Triangular");

        // Displaying the result :
        for (int i = 0; i < n; i++) 
        {
            // Lower
            for (int j = 0; j < n; j++)
                Console.Write(setw(4) + lower[i, j] + "\t");
            Console.Write("\t");

            // Upper
            for (int j = 0; j < n; j++)
                Console.Write(setw(4) + upper[i, j] + "\t");
            Console.Write("\n");
        }
    }
    static String setw(int noOfSpace)
    {
        s = "";
        for (int i = 0; i < noOfSpace; i++)
            s += " ";
        return s;
    }

    // Driver code
    public static void Main(String[] arr)
    {
        int[, ] mat = { { 2, -1, -2 },
                        { -4, 6, 3 },
                        { -4, -2, 8 } };

        luDecomposition(mat, 3);
    }
}

// This code is contributed by Princi Singh

JavaScript

<script>

// Javascript Program to decompose a matrix
// into lower and upper triangular matrix    
// function MAX = 100;
var s = "";

function luDecomposition(mat, n)
{
    var lower = Array(n).fill(0).map(
           x => Array(n).fill(0));
    var upper = Array(n).fill(0).map(
           x => Array(n).fill(0));

    // Decomposing matrix into Upper and
    // Lower triangular matrix
    for(var i = 0; i < n; i++) 
    {
        
        // Upper Triangular
        for(var k = i; k < n; k++)
        {
            
            // Summation of L(i, j) * U(j, k)
            var sum = 0;
            for(var 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(var k = i; k < n; k++)
        {
            if (i == k)
            
                // Diagonal as 1
                lower[i][i] = 1; 
            else 
            {
                
                // Summation of L(k, j) * U(j, i)
                var sum = 0;
                for(var j = 0; j < i; j++)
                    sum += (lower[k][j] * upper[j][i]);

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

    // Setw is for displaying nicely
    document.write(setw(2) + "Lower Triangular" + 
                   setw(10) + "Upper Triangular<br>");

    // Displaying the result :
    for(var i = 0; i < n; i++) 
    {
       
        // Lower
        for(var j = 0; j < n; j++)
            document.write(setw(4) + 
                           lower[i][j] + "\t");
                           
        document.write(setw(10));

        // Upper
        for(var j = 0; j < n; j++)
            document.write(setw(4) + 
                           upper[i][j] + "\t");
                           
        document.write("<br>");
    }
}

function setw(noOfSpace) 
{
    var s = "";
    for(i = 0; i < noOfSpace; i++)
        s += "&nbsp;";
        
    return s;
}

// Driver code
var mat = [ [ 2, -1, -2 ], 
            [ -4, 6, 3 ], 
            [ -4, -2, 8 ] ];

luDecomposition(mat, 3);

// This code is contributed by Amit Katiyar 

</script>

PHP

<?php
// PHP Program to decompose 
// a matrix into lower and
// upper triangular matrix
$MAX = 100;

function luDecomposition($mat, $n)
{
    $lower;
    $upper;
    for($i = 0; $i < $n; $i++)
    for($j = 0; $j < $n; $j++)
    {
        $lower[$i][$j]= 0;
        $upper[$i][$j]= 0;
    }
    // Decomposing matrix 
    // into Upper and Lower
    // triangular matrix
    for ($i = 0; $i < $n; $i++) 
    {

        // Upper Triangular
        for ($k = $i; $k < $n; $k++) 
        {

            // Summation of 
            // L(i, j) * U(j, k)
            $sum = 0;
            for ($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 ($k = $i; $k < $n; $k++) 
        {
            if ($i == $k)
                $lower[$i][$i] = 1; // Diagonal as 1
            else
            {

                // Summation of L(k, j) * U(j, i)
                $sum = 0;
                for ($j = 0; $j < $i; $j++)
                    $sum += ($lower[$k][$j] * 
                             $upper[$j][$i]);

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

    // setw is for 
    // displaying nicely
    echo "\t\tLower Triangular";
    echo "\t\t\tUpper Triangular\n";

    // Displaying the result :
    for ($i = 0; $i < $n; $i++) 
    {
        // Lower
        for ($j = 0; $j < $n; $j++)
            echo "\t" . $lower[$i][$j] . "\t"; 
        echo "\t";

        // Upper
        for ($j = 0; $j < $n; $j++)
        echo $upper[$i][$j] . "\t";
        echo "\n";
    }
}

// Driver code
$mat = array(array(2, -1, -2),
             array(-4, 6, 3),
             array(-4, -2, 8));

luDecomposition($mat, 3);

// This code is contributed by mits
?>

Output

      Lower Triangular                Upper Triangular
     1         0         0             2        -1        -2    
    -2         1         0             0         4        -1    
    -2        -1         1             0         0         3    

Applications on Doolittle Algorithm

1. Solving Linear Equations: LU Decomposition simplifies solving AX = B by solving LY = B and then UX = Y.
2. Matrix Inversion: Decompose A into LU, then invert L and U.
3. Determinant Calculation: The determinant of A is the product of the diagonal elements of U.

Conclusion - Doolittle Algorithm

The Doolittle Algorithm is an efficient method for LU Decomposition, enabling various applications in numerical analysis and linear algebra.