In linear algebra, a matrix decomposition or matrix factorization is a factorization of a matrix into a product of matrices. There are many different matrix decompositions. One of them is Cholesky Decomposition.
What is Cholesky Decomposition?
The Cholesky decomposition or Cholesky factorization is a decomposition of a Hermitian, positive-definite matrix into the product of a lower triangular matrix and its conjugate transpose. The Cholesky decomposition is roughly twice as efficient as the LU decomposition for solving systems of linear equations.
The Cholesky decomposition of a Hermitian positive-definite matrix A is a decomposition of the form A = [L][L]T, where L is a real lower triangular matrix with positive diagonal entries, and LT denotes the conjugate transpose of L.
Every Hermitian positive-definite matrix (and thus also every real-valued symmetric positive-definite matrix) has a unique Cholesky decomposition.
Every symmetric, positive definite matrix A can be decomposed into a product of a unique lower triangular matrix L and its transpose: A = L LT
The following formulas are obtained by solving above lower triangular matrix and its transpose. These are the basis of Cholesky Decomposition Algorithm :
L_{j,j}=\sqrt {A_{j,j}-\sum_{k=0}^{j-1}(L_{j,k})^{2}}
Example :
Input :
\left[\begin{array}{ccc} 4 & 12 & -16 \\ 12 & 37 & -43 \\ -16 & -43 & 98 \end{array}\right] Output :
\left[\begin{array}{ccc} 2 & 0 & 0 \\ 6 & 1 & 0 \\ -8 & 5 & 3 \end{array}\right]\left[\begin{array}{ccc} 2 & 6 & -8 \\ 0 & 1 & 5 \\ 0 & 0 & 3 \end{array}\right]
Below is the implementation of Cholesky Decomposition:
// CPP program to find decompose of a
// matrix using Cholesky Decomposition
#include <bits/stdc++.h>
using namespace std;
void choleskyDecomposition(vector<vector<int>> matrix){
int n = matrix.size();
// to store the lower triangular matrix
vector<vector<int>> lower(n, vector<int>(n, 0));
// Decomposing a matrix into Lower Triangular
for (int i = 0; i < n; i++) {
for (int j = 0; j <= i; j++) {
int sum = 0;
// summation for diagonals
if (j == i) {
for (int k = 0; k < j; k++)
sum += pow(lower[j][k], 2);
lower[j][j] = sqrt(matrix[j][j] - sum);
} else {
// Evaluating L(i, j) using L(j, j)
for (int k = 0; k < j; k++)
sum += (lower[i][k] * lower[j][k]);
lower[i][j] = (matrix[i][j] - sum) / lower[j][j];
}
}
}
// Displaying Lower Triangular Matrix
for (int i = 0; i < n; i++) {
// Lower Triangular
for (int j = 0; j < n; j++)
cout << lower[i][j] << " ";
cout << endl;
}
cout<<endl;
// Displaying Transpose of Lower Triangular Matrix
for (int i = 0; i < n; i++) {
// Lower Triangular
for (int j = 0; j < n; j++)
cout << lower[j][i] << " ";
cout << endl;
}
}
int main() {
vector<vector<int>> matrix = {
{ 4, 12, -16 },
{ 12, 37, -43 },
{ -16, -43, 98 }
};
choleskyDecomposition(matrix);
return 0;
}
// Java program to find decompose of a
// matrix using Cholesky Decomposition
import java.util.*;
class GfG {
static void choleskyDecomposition(List<List<Integer>> matrix) {
int n = matrix.size();
// to store the lower triangular matrix
int[][] lower = new int[n][n];
// Decomposing a matrix into Lower Triangular
for (int i = 0; i < n; i++) {
for (int j = 0; j <= i; j++) {
int sum = 0;
// summation for diagonals
if (j == i) {
for (int k = 0; k < j; k++)
sum += Math.pow(lower[j][k], 2);
lower[j][j] = (int) Math.sqrt(matrix.get(j).get(j) - sum);
} else {
// Evaluating L(i, j) using L(j, j)
for (int k = 0; k < j; k++)
sum += (lower[i][k] * lower[j][k]);
lower[i][j] = (matrix.get(i).get(j) - sum) / lower[j][j];
}
}
}
// Displaying Lower Triangular Matrix
for (int i = 0; i < n; i++) {
// Lower Triangular
for (int j = 0; j < n; j++)
System.out.print(lower[i][j] + " ");
System.out.println();
}
System.out.println();
// Displaying Transpose of Lower Triangular Matrix
for (int i = 0; i < n; i++) {
// Lower Triangular
for (int j = 0; j < n; j++)
System.out.print(lower[j][i] + " ");
System.out.println();
}
}
public static void main(String[] args) {
List<List<Integer>> matrix = Arrays.asList(
Arrays.asList(4, 12, -16),
Arrays.asList(12, 37, -43),
Arrays.asList(-16, -43, 98)
);
choleskyDecomposition(matrix);
}
}
# Python program to find decompose of a
# matrix using Cholesky Decomposition
def choleskyDecomposition(matrix):
n = len(matrix)
# to store the lower triangular matrix
lower = [[0] * n for _ in range(n)]
# Decomposing a matrix into Lower Triangular
for i in range(n):
for j in range(i + 1):
sum = 0
# summation for diagonals
if j == i:
for k in range(j):
sum += lower[j][k]**2
lower[j][j] = int((matrix[j][j] - sum)**0.5)
else:
# Evaluating L(i, j) using L(j, j)
for k in range(j):
sum += lower[i][k] * lower[j][k]
lower[i][j] = (matrix[i][j] - sum) // lower[j][j]
# Displaying Lower Triangular Matrix
for i in range(n):
print(" ".join(map(str, lower[i])))
print()
# Displaying Transpose of Lower Triangular Matrix
for i in range(n):
print(" ".join(map(str, [lower[j][i] for j in range(n)])))
if __name__ == "__main__":
matrix = [
[4, 12, -16],
[12, 37, -43],
[-16, -43, 98]
]
choleskyDecomposition(matrix)
// C# program to find decompose of a
// matrix using Cholesky Decomposition
using System;
using System.Collections.Generic;
class GfG {
static void choleskyDecomposition(List<List<int>> matrix) {
int n = matrix.Count;
// to store the lower triangular matrix
int[,] lower = new int[n, n];
// Decomposing a matrix into Lower Triangular
for (int i = 0; i < n; i++) {
for (int j = 0; j <= i; j++) {
int sum = 0;
// summation for diagonals
if (j == i) {
for (int k = 0; k < j; k++)
sum += (int)Math.Pow(lower[j, k], 2);
lower[j, j] = (int)Math.Sqrt(matrix[j][j] - sum);
} else {
// Evaluating L(i, j) using L(j, j)
for (int k = 0; k < j; k++)
sum += (lower[i, k] * lower[j, k]);
lower[i, j] = (matrix[i][j] - sum) / lower[j, j];
}
}
}
// Displaying Lower Triangular Matrix
for (int i = 0; i < n; i++) {
// Lower Triangular
for (int j = 0; j < n; j++)
Console.Write(lower[i, j] + " ");
Console.WriteLine();
}
Console.WriteLine();
// Displaying Transpose of Lower Triangular Matrix
for (int i = 0; i < n; i++) {
// Lower Triangular
for (int j = 0; j < n; j++)
Console.Write(lower[j, i] + " ");
Console.WriteLine();
}
}
static void Main(string[] args) {
List<List<int>> matrix = new List<List<int>> {
new List<int> { 4, 12, -16 },
new List<int> { 12, 37, -43 },
new List<int> { -16, -43, 98 }
};
choleskyDecomposition(matrix);
}
}
// JavaScript program to find decompose of a
// matrix using Cholesky Decomposition
function choleskyDecomposition(matrix) {
let n = matrix.length;
// to store the lower triangular matrix
let lower = Array.from({ length: n }, () => Array(n).fill(0));
// Decomposing a matrix into Lower Triangular
for (let i = 0; i < n; i++) {
for (let j = 0; j <= i; j++) {
let sum = 0;
// summation for diagonals
if (j === i) {
for (let k = 0; k < j; k++)
sum += Math.pow(lower[j][k], 2);
lower[j][j] = Math.sqrt(matrix[j][j] - sum) | 0;
} else {
// Evaluating L(i, j) using L(j, j)
for (let k = 0; k < j; k++)
sum += lower[i][k] * lower[j][k];
lower[i][j] = ((matrix[i][j] - sum) / lower[j][j]) | 0;
}
}
}
// Displaying Lower Triangular Matrix
for (let i = 0; i < n; i++) {
console.log(lower[i].join(" "));
}
console.log();
// Displaying Transpose of Lower Triangular Matrix
for (let i = 0; i < n; i++) {
let row = [];
for (let j = 0; j < n; j++) {
row.push(lower[j][i]);
}
console.log(row.join(" "));
}
}
let matrix = [
[4, 12, -16],
[12, 37, -43],
[-16, -43, 98]
];
choleskyDecomposition(matrix);
Output
2 0 0 6 1 0 -8 5 3 2 6 -8 0 1 5 0 0 3
Time Complexity: O(n^3) as 3 nested for loops from 0 to n-1 are being used.
Auxiliary Space: O(n^2) to store the values of the lower triangular matrix.
