Definition. $A \in \mathbb{R}^{n \times n}$ is called symmetric, if $A = A^T$. In addition $A$ is called positive definite, if $x^TAx > 0$ for all $x \in \mathbb{R}^{n}$, $x \neq 0$.

Remark. If $A \in \mathbb{R}^{n \times n}$ is symmetric and positive definite then $A$ is regular and $A^{-1}$ is symmetric and positive definite, as well.

Using the Cholesky decomposition a symmetric and positive definite matrix $A \in \mathbb{R}^{n \times n}$ can be decomposed into the product $LDL^T$ of a normed lower triangular matrix $L \in \mathbb{R}^{n \times n}$ and an upper triangular matrix $R \in \mathbb{R}^{n \times n}$.

Applications.

• The problem $Ax = LDL^Tx = b$ is reduced to solving $Ly = b$ and $L^Tx = D^{-1}y$.
• The used algorithm can be used to check whether the given matrix $A$ is symmetric and positive definite.

/**
* Calculate a cholesky decomposition.
*
* @author  David Stutz
*/
class Cholesky {

const TOLERANCE = 0.00001;

/**
* @var matrix
*/
protected $_matrix; /** * Constructor: Generate cholesky deocmposition of given matrix. * * @param matrix */ public function __construct(&$matrix) {
new Assertion($matrix instanceof Matrix, 'Given matrix not of class Matrix.'); new Assertion($matrix->isSquare(), 'Given matrix is not square.');

$this->_matrix =$matrix->copy();

for ($j = 0;$j < $this->_matrix->columns();$j++) {
$d =$this->_matrix->get($j,$j);
for ($k = 0;$k < $j;$k++) {
$d -= pow($this->_matrix->get($j,$k), 2) * $this->_matrix->get($k, $k); } // Test if symmetric and positive definite can be guaranteed. new Assertion($d > Cholesky::TOLERANCE * (double)$this->_matrix->get($j, $j), 'Symmetric and positive definite can not be guaranteed: ' .$d . ' > ' . Cholesky::TOLERANCE * (double)$this->_matrix->get($j, $j));$this->_matrix->set($j,$j, $d); for ($i = $j + 1;$i < $this->_matrix->rows();$i++) {
$this->_matrix->set($i, $j,$this->_matrix->get($i,$j));
for ($k = 0;$k < $j;$k++) {
$this->_matrix->set($i, $j,$this->_matrix->get($i,$j) - $this->_matrix->get($i, $k) *$this->_matrix->get($k,$k) * $this->_matrix->get($j, $k)); }$this->_matrix->set($i,$j, $this->_matrix->get($i, $j) / ((double)$this->_matrix->get($j,$j)));
}
}
}

/**
* Get the L of the composition L^T*D*L.
*
* @return  matrix  L
*/
public function getL() {
$L =$this->_matrix->copy();

// L is the lower triangular matrix.
for ($i = 0;$i < $L->rows();$i++) {
for ($j =$i; $j <$L->columns(); $j++) { if ($j == $i) {$L->set($i,$j, 1);
}
else {
$L->set($i, $j, 0); } } } return$L;
}

/**
* Get D - the diagonal matrix.
*
* @return  matrix  D
*/
public function getD() {
$D = new Matrix($this->_matrix->rows(), $this->_matrix->columns()); for ($i = 0; $i <$D->rows(); $i++) {$D->set($i,$i, $this->_matrix->get($i, $i)); } return$D;
}
}


There are explicit formulars to compute the entries for $L$ and $D$:

• $d_{j,k} = a_{j,j} - \sum _{k = 1} ^{j - 1} l_{j,k}^2 d_{k,k}$ for the columns $j = 1, \ldots, n$.
• $l_{i,j} = \frac{a_{i,j} - \sum _{k = 1} ^{j - 1} l_{i,k} d_{k,k} l_{j,k}}{d_{j,j}}$ for the columns $j = 1, \ldots, n$.

The formulars can be seen as result of an entry-wise comparison:

$A = \left[\begin{array}{c c c} a_{1,1} & \ldots & a_{1,n} \\ \vdots & \ddots & \vdots \\ a_{n,1} & \ldots & a_{n,n} \\ \end{array} \right] = \left[\begin{array}{c c c} l_{1,1} & \ldots & l_{1,n} \\ \vdots & \ddots & \vdots \\ l_{n,1} & \ldots & l_{n,n} \\ \end{array} \right] \left[\begin{array}{c c c} d_{1,1} & & \varnothing \\ & \ddots & \\ \varnothing & & d_{n,n} \\ \end{array} \right] \left[\begin{array}{c c c} l_{1,1} & \ldots & l_{n,1} \\ \vdots & \ddots & \vdots \\ l_{1,n} & \ldots & l_{n,n} \\ \end{array} \right] = LDL^T$

Based on the symmetry of $A$ the entry-wise analysis can be restricted to the lower triangular part of the matrix.

The algorithm can easily be derived from the above formulars:

Algorithm.

• For $j = 1,2, \ldots, n$:
• $d := a _{j,j} - \sum _{k = 1} ^{j - 1} a _{j,k}^2 a_{k,k}$
• If $d > \epsilon \cdot a_{j,j}$:
• $a_{j,j} := d$
• For $i = j+1, \ldots, n$:
• $a_{i,j} := \frac{a_{i,j} - \sum _{k = 1} ^{j - 1} a_{i,k} a_{k,k} a_{j,k}}{a_{j,j}}$

For each column $j$: Computing the diagonal entry needs $j - 1$ substractions and $2(j - 1)$ multiplications. Computing the entries of $L$ needs $(n - j)$ divisions, $(n - j) (j - 1)$ additions and $2 (n - j) (j - 1)$ multiplications:

Runtime.

$T(n) = \sum _{j = 1} ^{n} 3 (j - 1) + (n - j) + 3 (n - j) (j - 1) = \ldots = 3 \frac{n(n - 1)}{2} + \frac{n(n + 1)}{2} + 3 \frac{n(n + 1)(2n + 1)}{6} \in \mathcal{O}(\frac{1}{3} n^3)$

© 2013 David Stutz - Credits - Impressum - Datenschutz