.. _cholesky:

Cholesky Decomposition for Positive Definite Systems
====================================================

Systems of linear equations of the form
:math:`\mathbf{A}\,\bm{x} = \bm{b}` are often symmetric and *positive
definite*, especially equations derived from physical applications such
as electrical, mechanical, and structural systems. We will see in
:ref:`projection` that the solutions to over-determined systems, such
as found with least squares regression problems (:ref:`leastSquare`),
may make use of positive definite matrices of the form
:math:`\mathbf{A}^T\mathbf{A}`. The diagonal values of a positive
definite matrix are always positive. The special symmetry of positive
definite matrices allows for a simple and numerically accurate factoring
called the *Cholesky decomposition* that returns an upper triangular
matrix, :math:`\bf{R}`. The transpose of :math:`\bf{R}` gives us a lower
triangular matrix. The product of the lower and upper triangular
matrices is the decomposition of :math:`\bf{A}` that we will use as we
did with LU decomposition to find the solution using a triangular
solver.

.. math::
   :label: eq-cholesky

   \mathbf{A} = \mathbf{R}^T\mathbf{R}.

The Cholesky decomposition can only be used for positive definite
matrices. A *positive definite* matrix :math:`\bf{A}` is real,
symmetric, and satisfies the property

.. math::
   :label: eq-pos-definite

   \bm{x}^T\mathbf{A}\,\bm{x} > 0, \; \text{for all nonzero } 
   \bm{x} \in \mathbb{R}^n.

The left side of equation :eq:`eq-pos-definite` reduces to a dot
product that yields a real, scalar number. Unfortunately, applying the
test from equation :eq:`eq-pos-definite` is difficult because it
allows for an infinite set of :math:`\bm{x}` vectors. An alternate
criterion for evaluating whether a matrix is positive definite consists
of three tests.

#. Is the matrix symmetric?

#. Are all elements along the matrix’s diagonal positive?

#. Are all of the eigenvalues of the matrix positive?

The procedure for finding the eigenvalues of a matrix is addressed in
:ref:`eig`, but we will not use the third test. A faster final test to
see if the Cholesky decomposition can be used is to attempt to use it.
The decomposition fails if the algorithm needs to take the square root
of a number that is zero or negative. So we can attempt to use the
Cholesky decomposition algorithm if the matrix is symmetric and has only
positive values along the diagonal. 

The development of the Cholesky decomposition algorithm comes from the
algebra of how :math:`\bf{A}` is computed as the product
:math:`\mathbf{R}^T\, \bf{R}`.

.. math::

   \mathbf{A} =  \begin{bmatrix}
                   r_{00}  & {} & {} & {} \\ 
                   r_{01}  & r_{11} & {} & {} \\ 
                   r_{02}  & r_{12} & r_{22} & {} \\ 
                   r_{03}  & r_{13} & r_{23} & r_{33}
                \end{bmatrix} \,
            \begin{bmatrix}
              r_{00}  & r_{01} & r_{02} & r_{03} \\ 
              {}  & r_{11} & r_{12} & r_{13} \\ 
              {}  & {} & r_{22} & r_{23} \\ 
              {}  & {} & {} & r_{33}
           \end{bmatrix}

The first row is simple because each element of :math:`\bf{A}` comes
from only one product.

.. math::
   :label: eq-cholesky1

   \begin{array}{lcl}
           a_{00} = r_{00}\,r_{00} &\quad   &r_{00} = \sqrt{a_{00}} \\
           a_{01} = r_{00}\,r_{01} &\quad   &r_{01} = a_{01} / r_{00} \\
           a_{02} = r_{00}\,r_{02} &\quad   &r_{02} = a_{02} / r_{00} \\
           a_{03} = r_{00}\,r_{03} &\quad   &r_{03} = a_{03} / r_{00}
       \end{array}

In the subsequent rows, we find the value of the diagonal element using
previously found values in the column of the diagonal element.

.. math::

   \begin{array}{lcl}
           a_{11} = r_{01}^2 + r_{11}^2 &\quad   
               &r_{11} = \sqrt{a_{11} - r_{01}^2} \\
           a_{22} = r_{02}^2 + r_{12}^2 + r_{22}^2   &\quad       
               &r_{22} = \sqrt{a_{22} - r_{02}^2 - r_{12}^2}
       \end{array}

.. math::
   :label: eq-cholesky-diag

   r_{ii} = \sqrt{a_{ii} - \sum_{k = 0}^{i-1}r_{ki}^2},
       \quad \text{for } i = 1 \ldots n-1

The algorithm uses equation :eq:`eq-cholesky-diag` to check if the
matrix might still be positive definite. If the argument to the square
root is zero or a negative number, the algorithm is finished because the
matrix is not positive definite. The remaining values of each row are
found after finding the diagonal element.

.. math::

   \begin{array}{lcl}
           a_{12} = r_{01}\,r_{02} + r_{11}\,r_{12}    
              &\quad     &r_{12} = \left(a_{12} - r_{01}\,r_{02}\right) / r_{11} \\
           a_{23} = r_{02}\,r_{03} + r_{12}\,r_{13} + r_{22}\,r_{23} 
          &\quad     &r_{23} = \left(a_{23} - \left(r_{02}\,r_{03} +
         r_{12}\,r_{13}\right)\right) / r_{22}
       \end{array}

.. math::
   :label: eq-cholesky-other

   r_{ij} = \left( a_{ij} - \sum_{k = 0}^{i-1}r_{ki}\,r_{kj}\right) /
       r_{ii},
       \quad \text{for } i = 0 \ldots n-1, \text{ and } j = i+1 \ldots n-1

The ``chol_decomp`` function uses equation :eq:`eq-cholesky1`,
equation :eq:`eq-cholesky-diag`, and equation :eq:`eq-cholesky-other` 
to find the Cholesky decomposition. 
It raises an exception if the matrix is not positive definite because it
is not square, symmetric, or has non-positive values on the diagonal.
It returns a boolean flag indicating if the algorithm could complete the
decomposition, which also tells us if the matrix is positive definite or
not.
A matrix may pass the obvious tests for positive definite pass, but may still
not be positive definite.

.. index:: chol_decomp

::

    # -*- coding: utf-8 -*-
    """
     File: chol_decomp.py
    """

    import numpy as np
    from scipy.linalg import issymmetric
        
    def chol_decomp(A: np.array) -> tuple[np.array, np.array]: 
        # Cholesky decomposition. If A is posive definite, A = R'*R.
        # isPositiveDefinite is a boolean variable. Use R only if
        # isPositiveDefinite is true.
        
        m,n = A.shape
        if m != n:
            raise Exception('Matrix is not square')
        else:
            if not issymmetric(A):
                raise Exception('Matrix is not symmetric')
        
        if np.any(np.diag(A <= 0)):
            msg = 'Matrix is not positive definite\n'
            msg = msg + 'Diagonal not all positive'
            raise Exception(msg)
        
        B = A.copy() # A is mutable, don't change it.
        isPositiveDefinite = True
        for i in np.arange(0,n):
            for k in np.arange(0,i): # skip i = 0
                B[i,i] = B[i,i] - B[k,i]**2
            if A[i,i] <= 0:
                isPositiveDefinite = False
                R = np.array([])
                return R,isPositiveDefinite
            B[i,i] = np.sqrt(B[i,i])
            for j in np.arange(i+1,n):
                for k in np.arange(0,i):
                    B[i,j] = B[i,j] - B[k,i] * B[k,j]
                B[i,j] = B[i,j] / B[i,i]
        
        R = np.triu(B)
        return R,isPositiveDefinite

    if __name__ == "__main__":
        A = np.random.rand(3,3)
        try:
            R,_ = chol_decomp(A)
        except Exception as e:
            print(f"Expected error occurred: {e}")
        try:
            B = A.T @ A
            R,isPD = chol_decomp(B)
            print(isPD)
            print(np.allclose(B, (R.T @ R)))
        except Exception as e:
            print(f"An error occurred: {e}")
        try:
            B[1,1] = -1
            R,isPD = chol_decomp(B)
            print(isPD)
        except Exception as e:
            print(f"Expected error occurred: {e}")
            

The Python ``scipy.linalg`` module has a few functions that make
use of the Cholesky decomposition. The ``cholesky`` function computes
the decomposition. It raises an exception if the decomposition fails. 
A triangular solver can be used when the decomposition is successful.


An alternative for when the factor is not needed is ``cho_factor``
and ``cho_solve``. The ``cho_factor`` returns the Cholesky factor and
Boolen variable indicating if the factor is lower or upper triangular. 

::

    # -*- coding: utf-8 -*-
    """
    File: cholesky_test.py
    Script to test our Cholesky decomposition and similar functions
    in the scipy.linalg module.
    """
    import numpy as np
    from scipy.linalg import cholesky, cho_factor, cho_solve, solve_triangular
    from chol_decomp import chol_decomp

    A = np.random.rand(4,4)
    
    # not positive definite case
    try:
        R,_ = chol_decomp(A)
    except Exception as e:
        print(f"Expected error occurred: {e}")
    try:
        R = cholesky(A)
    except Exception as e:
        print(f"Expected error occurred: {e}") 

    # positive definite case
    b = np.random.rand(4,1)
    B = A.T @ A
    try:
        R = cholesky(B)
        R,isPD = chol_decomp(B)
        print(isPD)
        print(np.allclose(B, (R.T @ R)))
        if isPD:
            x = solve_triangular(R, solve_triangular(R.T, b, lower=True))
            # verify correct
            print(np.allclose((B@x - b), np.zeros((4, 1))))
    except Exception as e:
        print(f"An error occurred: {e}")
        
    # Use cho_factor and cho_solve
    try:
        R,lower = cho_factor(B)
        x = cho_solve((R,lower), b)
        # verify correct
        print(np.allclose((B@x - b), np.zeros((4, 1))))
    except Exception as e:
        print(f"An error occurred: {e}")

.. index:: chol
.. index:: positive definite matrix
.. index:: Cholesky decomposition