Cholesky decomposition: Python & Matlab

Cholesky decomposition is an efficient method for inversion of symmetric positive-definite matrices. Let’s demonstrate the method in Python and Matlab.

Cholesky factor

Any symmetric positive definite matrix can be factored as

where is lower triangular matrix. The lower triangular matrix is often called “Cholesky Factor of ”. The matrix can be interpreted as square root of the positive definite matrix .

Basic Algorithm to find Cholesky Factorization:

Note: In the following text, the variables represented in Greek letters represent scalar values, the variables represented in small Latin letters are column vectors and the variables represented in capital Latin letters are Matrices.

Given a positive definite matrix , it is partitioned as follows.

first element of and respectively
column vector at first column starting from second row of and respectively
remaining lower part of the matrix of and respectively of size

Having partitioned the matrix as shown above, the Cholesky factorization can be computed by the following iterative procedure.

Steps in computing the Cholesky factorization:

Step 1: Compute the scalar:
Step 2: Compute the column vector:
Step 3: Compute the matrix :
Step 4: Replace with , i.e,
Step 5: Repeat from step 1 till the matrix size at Step 4 becomes .

Matlab Program (implementing the above algorithm):

Function 1: [F]=cholesky(A,option)

function [F]=cholesky(A,option)
%Function to find the Cholesky factor of a Positive Definite matrix A
%Author Mathuranathan for https://www.gaussianwaves.com
%Licensed under Creative Commons: CC-NC-BY-SA 3.0
%A = positive definite matrix
%Option can be one of the following 'Lower','Upper'
%L = Cholesky factorizaton of A such that A=LL^T
%If option ='Lower', then it returns the Cholesky factored matrix L in
%lower triangular form
%If option ='Upper', then it returns the Cholesky factored matrix L in
%upper triangular form    

%Test for positive definiteness (symmetricity need to satisfy)
%Check if the matrix is symmetric
if ~isequal(A,A'),
    error('Input Matrix is not Symmetric');
end
    
if isPositiveDefinite(A),
    [m,n]=size(A);
    L=zeros(m,m);%Initialize to all zeros
    row=1;col=1;
    j=1;
    for i=1:m,
        a11=sqrt(A(1,1));
        L(row,col)=a11;
        if(m~=1), %Reached the last partition
            L21=A(j+1:m,1)/a11;
            L(row+1:end,col)=L21;
            A=(A(j+1:m,j+1:m)-L21*L21');
            [m,n]=size(A);
            row=row+1;
            col=col+1;
        end
    end
    switch nargin
        case 2
            if strcmpi(option,'upper'),F=L';
            else
                if strcmpi(option,'lower'),F=L;
                else error('Invalid option');
                end
            end
        case 1
            F=L;
        otherwise
            error('Not enough input arguments')
    end
else
        error('Given Matrix A is NOT Positive definite');
end
end

Function 2: x=isPositiveDefinite(A)

function x=isPositiveDefinite(A)
%Function to check whether a given matrix A is positive definite
%Author Mathuranathan for https://www.gaussianwaves.com
%Licensed under Creative Commons: CC-NC-BY-SA 3.0
%Returns x=1, if the input matrix is positive definite
%Returns x=0, if the input matrix is not positive definite
[m,~]=size(A);
    
%Test for positive definiteness
x=1; %Flag to check for positiveness
for i=1:m
    subA=A(1:i,1:i); %Extract upper left kxk submatrix
    if(det(subA)<=0); %Check if the determinent of the kxk submatrix is +ve
        x=0;
        break;
    end
end
%For debug
%if x, display('Given Matrix is Positive definite');
%else display('Given Matrix is NOT positive definite'); end    
end

Sample Run:

A is a randomly generated positive definite matrix. To generate a random positive definite matrix check the link in “external link” section below.

>> A=[3.3821 ,0.8784,0.3613,-2.0349; 0.8784, 2.0068, 0.5587, 0.1169; 0.3613, 0.5587, 3.6656, 0.7807; -2.0349, 0.1169, 0.7807, 2.5397];

>> cholesky(A,’Lower’) 

>> cholesky(A,’upper’) 

Python (numpy)

Let us verify the above results using Python’s Numpy package. The numpy package numpy.linalg contains the cholesky function for computing the Cholesky decomposition (returns in lower triangular matrix form). It can be summoned as follows

>>> import numpy as np

>>> A=[[3.3821 ,0.8784,0.3613,-2.0349],\
       [0.8784, 2.0068, 0.5587, 0.1169],\
       [0.3613, 0.5587, 3.6656, 0.7807],\
       [-2.0349, 0.1169, 0.7807, 2.5397]]
>>> np.linalg.cholesky(A)
>>>
array([[ 1.83904867,  0.        ,  0.        ,  0.        ],
       [ 0.47763826,  1.33366476,  0.        ,  0.        ],
       [ 0.19646027,  0.34856065,  1.87230041,  0.        ],
       [-1.106496  ,  0.48393333,  0.44298574,  0.94071184]])

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Related Topics

[1]	An Introduction to Estimation Theory
[2]	Bias of an Estimator
[3]	Minimum Variance Unbiased Estimators (MVUE)
[4]	Maximum Likelihood Estimation
[5]	Maximum Likelihood Decoding
[6]	Probability and Random Process
[7]	Likelihood Function and Maximum Likelihood Estimation (MLE)
[8]	Score, Fisher Information and Estimator Sensitivity
[9]	Introduction to Cramer Rao Lower Bound (CRLB)
[10]	Cramer Rao Lower Bound for Scalar Parameter Estimation
[11]	Applying Cramer Rao Lower Bound (CRLB) to find a Minimum Variance Unbiased Estimator (MVUE)
[12]	Efficient Estimators and CRLB
[13]	Cramer Rao Lower Bound for Phase Estimation
[14]	Normalized CRLB - an alternate form of CRLB and its relation to estimator sensitivity
[15]	Cramer Rao Lower Bound (CRLB) for Vector Parameter Estimation
[16]	The Mean Square Error – Why do we use it for estimation problems
[17]	How to estimate unknown parameters using Ordinary Least Squares (OLS)
[18]	Essential Preliminary Matrix Algebra for Signal Processing
[19]	Why Cholesky Decomposition ? A sample case:
[20]	Tests for Positive Definiteness of a Matrix
[21]	Solving a Triangular Matrix using Forward & Backward Substitution
[22]	Cholesky Factorization - Matlab and Python
[23]	LTI system models for random signals – AR, MA and ARMA models
[24]	Comparing AR and ARMA model - minimization of squared error
[25]	Yule Walker Estimation
[26]	AutoCorrelation (Correlogram) and persistence – Time series analysis
[27]	Linear Models - Least Squares Estimator (LSE)
[28]	Best Linear Unbiased Estimator (BLUE)

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. $14.99 – Add to Cart Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. $14.99 – Add to Cart Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. $14.99 – Add to Cart Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing