Source code for spectrum.cholesky

.. topic:: Cholesky methods

    .. autosummary::


    .. 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



    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