Topic 18. Cholesky Factorization¶

Author: David Wood, d.wood0403@gmail.com; John Hunt, john.hunt92@gmail.com

Introduction¶

The factorization is named after Andre-Louis Cholesky, a French military officer involved in geodetics engineering. Andre-Louis Cholesky lived from 1875 to 1918, killed in World War I (WWI). He attented the renowned university École Polytechnique in France, where the first discoveries about radioactivity occurred. Cholesky survyed Crete and North Africa before WWI. Before the war Cholesky created and used the Cholesky decomposition in his survying.

The Cholesky decomposition is the decomposition of a Hermitian, positive-definite matrix into the product of a lower triangular matrix and its conjugate transpose, which is useful for efficient numerical solutions, e.g. Monte Carlo simulations and systems with multiple correlated variables (a system with independant variables would result in matrices that are not Hermitian, positive-deffinite and thus non solvable with this method).

The Cholesky decomposition is also closely related to the solution of least-squares problems, since the normal equations that characterize the least-squares solution have a symmetric positive-definite coefficient matrix. When it is applicable, the Cholesky decomposition is roughly twice as efficient as the LU decomposition for solving systems of linear equations. It can also be used to calculate and use correlation matrices in experiment statistics. In the case of a scalar (n = 1), the Cholesky factor $R$ is just the positive square root of $A$. However, $R$ should not be confused with the square roots of $A$.

Explanation of the theory¶

The Cholesky factorization can be thought of, with caution (see above), as finding the general square root of matrix $A$. Using this type of factorization, we factor $A$ into a single matrix $L$ multiplied by itself $A=L L^{T}$ or $A=R^{T} R$. In the case of a scalar (n = 1), the Cholesky factor $R$ is just the positive square root of $A$. However, $R$ should not be confused with the square roots of $A$.

Every positive definite matrix $A$ has a Cholesky factorization if there exists a lower trianglar matrix $L$ that contains positive diagonal elements such that $A=L L^{T}$. In summary, Cholesky factorization is a specific form of LU decomposition where $U=L^{T}$. This fact makes computing $L$ much easier.

Theorem : Every positive definite matrix $A$ has a Cholesky factorization which may be constructed.

Proof: Using induction on an $n-1 \space \times \space n-1$ matrix, we can prove the result for the trivial $1 \space \times \space 1$ positive definite matrix $A$ and extend it to an $n \space \times \space n$ matrix. For the $1 \space \times \space 1$ positive definite matrix $A = [a_{11}]$, $a_{11} > 0$. Therefore $L=[l_{11}]$ where $l_{11}=sqrt(a_{11})$. This result then can be extended to any positive definite matrix $A$.

Since $A$ is positive definite, it is symmetric and can be represented as

$$A = \begin{bmatrix} a_{11} & A_{21}^{T} \\ A_{21} & A_{22} \\ \end{bmatrix}$$

Since $A$ is positive definite, $a_{11} > 0$, and thus we can easily determine that $l_{11} = \sqrt{a_{11}}$ and $L_{21} = \frac{1}{l_{11}} \space A_{21}$.

We can now show that

$$K=A_{22} \space - \space L_{21} L_{21}^{T}$$

is a positive definite matrix as well. For any $n-1 \space \times \space 1$ column vector $X \neq 0$ and let

$$V = -\frac{1}{a_{11}} A_{21}^{T} X$$

Then

$$X^{T}KX \space = \space [V \space X^{T}] A \begin{bmatrix} V \\ X \\ \end{bmatrix} \space > \space 0$$

because $A$ is positive definite. Since

$$A_{22}-L_{21}L_{21}^{T}$$

is a positive definite $n-1 \space \times \space n-1$ matrix, there is a lower triangular matrix that contains all positive elements such that

$$A_{22}-L_{21}L_{21}^{T}=L_{22}L_{22}^{T}$$

Therefore

$$A = \begin{bmatrix} a_{11} & A_{21}^{T} \\ A_{21} & A_{22} \\ \end{bmatrix} = \begin{bmatrix} l_{11} & 0 \\ L_{21} & L_{22} \\ \end{bmatrix} \begin{bmatrix} l_{11} & L_{21}^{T} \\ 0 & L_{22}^{T} \\ \end{bmatrix}$$

Complete Algebra Using a 3x3 matrix, we will show the Cholesky factorization in its entirety which may be coded. As an example, we have to solve the following system of equations:

$$A = \begin{bmatrix} a_{11} & a_{21} & a_{31} \\ a_{21} & a_{22} & a_{32} \\ a_{31} & a_{32} & a_{33} \\ \end{bmatrix} = \begin{bmatrix} l_{11} & 0 & 0\\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \\ \end{bmatrix} \begin{bmatrix} l_{11} & l_{21} & l_{31} \\ 0 & l_{22} & l_{32} \\ 0 & 0 & l_{33} \\ \end{bmatrix} = LL^{T}$$

$$= \begin{bmatrix} l_{11}^{2} & l_{11} l_{21} & l_{11} a_{31} \\ l_{11} l_{21} & l_{21}^{2} + l_{22}^{2} & l_{31} l_{21} + l_{32} l_{22} \\ l_{11} l_{31} & l_{31} l_{21} + l_{32} l_{22} & l_{31}^{2} + l_{32}^{2} + l_{33}^{2} \\ \end{bmatrix}$$

We can easily see that the diagonal elements have a pattern. In general:

$$l_{ii} = \sqrt{a_{ii} - \sum_{j=1}^{i-1} l_{ij}^{2}}$$

For off-diagonal elements ($l_{ij}$) there is also a pattern:

$$l_{21}=\frac{1}{l_{11}} a_{21}$$ $$l_{31}=\frac{1}{l_{11}} a_{31}$$ $$l_{32}=\frac{1}{l_{22}} (a_{32} - l_{21}l_{31})$$

which in general terms is:

$$l_{ij} = \frac{1}{l_{jj}}(a_{ij} - \sum_{k=1}^{j-1} l_{ik}l_{jk}) \qquad i > j$$

Use in linear algebra To solve linear systems of equations, Cholesky factorization can be used when computing $A^{-1}$ is difficult or computationally expensive - as is the case for large systems. To solve equations in $Ax=b$ form, we only have to perform one intermediate step before solving for $\textbf{x}$. This comes from the fact that once $A$ is factored into $L L^{T}$ we can redefine the equation with a temporary vector $\textbf{y}$

$$Ax \space = \space L(L^{T}x) \space = \space Ly \space = \space b$$

Once $A$ is factored into $L L^{T}$ we must use forward substitution to solve for the temporary vector $\textbf{y}$ using

$$Ly \space = \space \begin{bmatrix} l_{11} & 0 \\ L_{21} & L_{22} \\ \end{bmatrix} \begin{bmatrix} y_{1} \\ y_{2} \\ \end{bmatrix} \space = \space \begin{bmatrix} b_{1} \\ b_{2} \\ \end{bmatrix} \space = \space b$$

Then, having found $\textbf{y}$, the next step is use backward substitution to solve for $\textbf{x}$ in

$$L^{T}x \space = \space \begin{bmatrix} l_{11} & L_{21} \\ 0 & L_{22} \\ \end{bmatrix} \begin{bmatrix} x_{1} \\ x_{2} \\ \end{bmatrix} \space = \space \begin{bmatrix} y_{1} \\ y_{2} \\ \end{bmatrix} \space = \space y$$

Both forward substitution and backward substitution can be performed by directly solving each equation in the system sequentially.

Additionally, the Cholesky factorization could also be written as $A=U^{H} U$. Where $U$ = $L^{H}$. Here $U$ is an upper triangle matrix instead of a lower triangle matrix. This could be useful if a solver for an upper triangle has already been implemented then the order of the matrix can change my taking the Hermitian of the L matrix.

Example:

$$A = \begin{bmatrix} 1 & 2 & 4 & 1 \\ 2 & 13 & 17 & 8 \\ 4 & 17 & 29 & 16 \\ 1 & 8 & 16 & 30 \\ \end{bmatrix} = LL^{T}$$ $$LL^{H} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & 3 & 0 & 0 \\ 4 & 3 & 2 & 0 \\ 1 & 2 & 3 & 4 \\ \end{bmatrix} \begin{bmatrix} 1 & 2 & 4 & 1 \\ 0 & 3 & 3 & 2 \\ 0 & 0 & 2 & 3 \\ 0 & 0 & 0 & 4 \\ \end{bmatrix}$$

$$L^{H} = U = \begin{bmatrix} 1 & 2 & 4 & 1 \\ 0 & 3 & 3 & 2 \\ 0 & 0 & 2 & 3 \\ 0 & 0 & 0 & 4 \\ \end{bmatrix}$$

$$A = \begin{bmatrix} 1 & 2 & 4 & 1 \\ 2 & 13 & 17 & 8 \\ 4 & 17 & 29 & 16 \\ 1 & 8 & 16 & 30 \\ \end{bmatrix} = U^{H}U = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & 3 & 0 & 0 \\ 4 & 3 & 2 & 0 \\ 1 & 2 & 3 & 4 \\ \end{bmatrix} \begin{bmatrix} 1 & 2 & 4 & 1 \\ 0 & 3 & 3 & 2 \\ 0 & 0 & 2 & 3 \\ 0 & 0 & 0 & 4 \\ \end{bmatrix}$$

By inspection it can be seen that $U^{H}U$ = $LL^{H}$

Simple Numerical Examples¶

Provide some simple python code and examples that emphasize the basic concepts.

python code of Cholesky decomposition¶

In its entirely, explicitly coding Cholesky factorization can be an arduous task and greatly taxes a computer's computational capacity if $A$ becomes too large. Thankfully, the SciPy library has a built in Cholesky decomposition function - scipy.linalg.cholesky - that does the calculations while requiring only a fraction of the computation power. Below is the same example using SciPy's built in functions.

As you can see, both implementations of Cholesky factorization gave the same answer. However, in the example were we created the Cholesky function we only have the lower triangle matrix. The transposed lower triangle matrix may also be calculated with additional coding.

Matlab code of Cholesky decomposition (copy and paste in Matlab for application)¶

%% Cholesky Factorization

% -- User Information (Matrix and olution vector)-- %

B = [1 2 4 1; 2 13 17 8; 4 17 29 16; 1 8 16 30];

b = [2;2;2;2];

% B = [6 15 55; 15 55 225; 55 225 979];

[row, col] = size(B);

[V, D] = eig(B);

d = [];

for i=1:length(D)

d = [d;D(i,i)];

end

if sum(d>0)==length(d) && row*row == sum(sum(B - B' == 0))% -- Check if matrix is PSD & Symmetry -- %

L = zeros(row,col);

for i = 1:row

    for j = 1:col

        if j > i

            L(i,j) = 0;

        elseif j == i

            num = 0;

            for k = 1:row-1

                num = num + L(i,k)^2;

            end

            L(i,j) = sqrt(B(i,j) - num);

        else

            num = 0;

            for k = 1:row-1

                num = num + L(j,k)*L(i,k);

            end

            L(i,j) = (B(i,j)-num)/L(j,j);

        end

    end

end

else

disp('Cannot use Cholesky Factorization b/c it is not symmetric or it is not positive semi definite')

end

% -- Solve -- %

y = L(b);

x = L'\y

ANSWER

L = [1,0,0,0;2,3,0,0;4,3,2,0;1,2,3,4]

x = [5.97222222222222; 1.15972222222222; -1.68750000000000; 0.458333333333333]

In addition $Matlab$ has prebuilt functions as well, below is an example using those functions.

A = [1 2 4 1; 2 13 17 8; 4 17 29 16; 1 8 16 30]

(or you can read in data from a file and continue)

L = chol(A)

y = L(b);

x = L'\y

The answer is the same as the example above.

An Engineering Application¶

Two main uses of the Cholesky decomposition in engineering are to solve a least squares problem - because the Cholesky factorization requires less computational power to compute and solve than QR or LU - and random sampling of a multivariate Normal distribution. It can be used to decompose a covariance matrix because these matrices may only be positive semi-definite - one of the prerequisites to using the Cholesky decomposition.

Here we will solve a problem for $\textbf{x}$ that is defined by a system of 3 equations.

Example 2: A system is defined by the equations below

$$\begin{bmatrix} 4 & 2 & 14 \\ 2 & 17 & -5 \\ 14 & -5 & 83 \\ \end{bmatrix} \begin{bmatrix} x_{1} \\ x_{2} \\ x_{3} \\ \end{bmatrix} = \begin{bmatrix} 14 \\ -101 \\ 155 \\ \end{bmatrix}$$

Using Cholesky decomposition, we find

$$L = \begin{bmatrix} 2 & 0 & 0 \\ 1 & 4 & 0 \\ 7 & -3 & 5 \\ \end{bmatrix}$$ $$L^{T} = \begin{bmatrix} 2 & 1 & 7 \\ 0 & 4 & -3 \\ 0 & 0 & 5 \\ \end{bmatrix}$$

First solving $Ly=b$ with forward substitution

$$\textbf{y} = \begin{bmatrix} 7 & -27 & 5 \\ \end{bmatrix}$$

then solving $L^{T}y=b$ with backward substitution, finally gives a result of

$$\textbf{x} = \begin{bmatrix} 3 & -6 & 1 \\ \end{bmatrix}$$

Example 2: We can also use the factorization to generate a random sampling from the multivariate Normal distribution $Z \space \approx \space N(0,I)$ where the affine transformation

$$Y \space = \space QX \space + \space \mu$$

is to be approximated. The variable $Y$ are the random samples that are generated, $X$ are the samples from the known univariate Normal distribution, and Q is defined as follows:

$$Q \space = \space \lambda^{1/2} \Phi$$

where $\lambda$ is a diagonal matrix containing the eigenvalues of the target matrix $\Sigma$, and $\Phi$ contains the associated eigenvectors of $\Sigma$ arranged to match their eigenvalues in $\lambda$.

Note that $Y$ has the distribution

$$Y \space \approx \space N(\mu,QQ')$$

We use the Cholesky decompositon because we want $\Sigma = QQ'$. Therefore, to simulate (or sample) $Z$ we can use $X \space \approx \space N(\mu,\Sigma)$ which is less computationally expensive to calculate.

Future Homework¶

Use the Cholesky decomposition to transform 10 normal distributed independent (and thus uncorrelated) random variables with mean of 0 and variance of 1 - N(0,1) to correlated variables in using Matlab or python libraries with the following covariance matrix (please copy and paste this in Matlab or python)

$$K_{xx} = \begin{bmatrix} 1 & 2 & 3 & 5 & -1 & 1 & -2 & 5 & 8 & 1 \\ 2 & 13 & 18 & 37 & 4 & 8 & 5 & 37 & 22 & 5 \\ 3 & 18 & 50 & 46 & 20 & 41 & 26 & 56 & 37 & 17 \\ 5 & 37 & 46 & 111 & 18 & 27 & 23 & 121 & 75 & 16 \\ -1 & 4 & 20 & 18 & 55 & 61 & 50 & 63 & 50 & 40 \\ 1 & 8 & 41 & 27 & 61 & 91 & 88 & 90 & 96 & 60 \\ -2 & 5 & 26 & 23 & 50 & 88 & 143 & 104 & 116 & 86 \\ 5 & 37 & 56 & 121 & 63 & 90 & 104 & 223 & 199 & 110 \\ 8 & 22 & 37 & 75 & 50 & 96 & 116 & 199 & 285 & 175 \\ 1 & 5 & 17 & 16 & 40 & 60 & 86 & 110 & 175 & 230 \\ \end{bmatrix}$$

and

$$SD = \begin{bmatrix} 1\\ 3.60555127546399 \\ 7.07106781186548 \\ 10.5356537528527 \\ 7.41619848709566 \\ 9.53939201416946 \\ 11.9582607431014 \\ 14.9331845230681 \\ 16.8819430161341 \\ 15.1657508881031 \\ \end{bmatrix}$$

 (Maybe put a file on LS?)

Please turn in the following:

1. A plot of the uncorrelated variables as well as the correlated variables (using subplot in Matlab or pairplot from seaborn library in python will be helpful).

2. An explination of why your plots shows that applying the cholesky decompotision correlates the random variables.

Future Homework Solutions¶

MATLAB CODE

% Solution

mu = 0; sigma = 1; sz = 10;

x_uncor = normrnd(mu,sigma,sz);

% Uncorrelated Plots

figure(1); count = 1;

for i = 1:10

for j = 1:10

    subplot(10,10,count)

    scatter(x_uncor(:,i),x_uncor(:,j))

    count = count + 1;

end

end

L = chol(Kxx);

x_cor = L*Kxx;

% Correlated Plots

figure(2); count = 1;

for i = 1:10

for j = 1:10

    subplot(10,10,count)

    scatter(x_cor(:,i),x_cor(:,j))

    count = count + 1;

end

end