Fork me on GitHub

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
 * @license http://www.gnu.org/licenses/gpl-3.0
 */
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