# Source code for spectrum.cholesky

```"""
.. topic:: Cholesky methods

.. autosummary::

CHOLESKY

.. codeauthor:: Thomas Cokelaer, 2011
"""
import numpy

__all__ = ["CHOLESKY"]

def _numpy_cholesky(A, B):
"""Solve Ax=B using numpy cholesky solver

A = LU

in the case where A is square and Hermitian, A = L.L* where L* is
transpoed and conjugate matrix

Ly = b

where

Ux=y

so x = U^{-1} y
where U = L*
and y = L^{-1} B
"""
L = numpy.linalg.cholesky(A)
# A=L*numpy.transpose(L).conjugate()
# Ly = b
y = numpy.linalg.solve(L,B)
# Ux = y
x = numpy.linalg.solve(L.transpose().conjugate(),y)
return x, L

def _numpy_solver(A, B):
"""This function solve Ax=B directly without taking care of the input
matrix properties.
"""
x = numpy.linalg.solve(A, B)
return x

[docs]def CHOLESKY(A, B, method='scipy'):
"""Solve linear system `AX=B` using CHOLESKY method.

:param A: an input Hermitian matrix
:param B: an array
:param str method: a choice of method in [numpy, scipy, numpy_solver]

* `numpy_solver` relies entirely on numpy.solver (no cholesky decomposition)
* `numpy` relies on the numpy.linalg.cholesky for the decomposition and
numpy.linalg.solve for the inversion.
* `scipy` uses scipy.linalg.cholesky for the decomposition and
scipy.linalg.cho_solve for the inversion.

.. rubric:: Description

When a matrix is square and Hermitian (symmetric with lower part being
the complex conjugate of the upper one), then the usual triangular
factorization takes on the special form:

.. math:: A = R R^H

where :math:`R` is a lower triangular matrix with nonzero real principal
diagonal element. The input matrix can be made of complex data. Then, the
inversion to find :math:`x` is made as follows:

.. math::  Ry = B

and

.. math::   Rx = y

.. doctest::

>>> import numpy
>>> from spectrum import CHOLESKY
>>> A = numpy.array([[ 2.0+0.j ,  0.5-0.5j, -0.2+0.1j],
...    [ 0.5+0.5j,  1.0+0.j ,  0.3-0.2j],
...    [-0.2-0.1j,  0.3+0.2j,  0.5+0.j ]])
>>> B = numpy.array([ 1.0+3.j ,  2.0-1.j ,  0.5+0.8j])
>>> CHOLESKY(A, B)
array([ 0.95945946+5.25675676j,  4.41891892-7.04054054j,
-5.13513514+6.35135135j])

"""
if method == 'numpy_solver':
X = _numpy_solver(A,B)
return X
elif method == 'numpy':
X, _L = _numpy_cholesky(A, B)
return X
elif method == 'scipy':
import scipy.linalg
L = scipy.linalg.cholesky(A)
X = scipy.linalg.cho_solve((L, False), B)
else:
raise ValueError('method must be numpy_solver, numpy_cholesky or cholesky_inplace')
return X
```