Sparse Cholesky decomposition (:mod:`sksparse.cholmod`)
=======================================================
.. module:: sksparse.cholmod
:synopsis: Sparse Cholesky decomposition using CHOLMOD
.. versionadded:: 0.1
Overview
--------
This module provides efficient implementations of all the basic linear
algebra operations for sparse, symmetric, positive-definite matrices
(as, for instance, commonly arise in least squares problems).
Specifically, it exposes most of the capabilities of the `CHOLMOD
`_ package,
including:
* Computation of the `Cholesky decomposition
`_ :math:`LL' =
A` or :math:`LDL' = A` (with fill-reducing permutation) for both
real and complex sparse matrices :math:`A`, in any format supported
by :mod:`scipy.sparse`. (However, CSC matrices will be most
efficient.)
* A convenient and efficient interface for using this decomposition to
solve problems of the form :math:`Ax = b`.
* The ability to perform the costly fill-reduction analysis once, and
then re-use it to efficiently decompose many matrices with the same
pattern of non-zero entries.
* In-place 'update' and 'downdate' operations, for computing the
Cholesky decomposition of a rank-k update of :math:`A` and of
product :math:`AA'`. So, the result is the Cholesky decomposition of
:math:`A + CC'` (or :math:`AA' + CC'`). The last case is useful when the
columns of `A` become available incrementally (e.g., due to memory
constraints), or when many matrices with similar but non-identical
columns must be factored.
* Convenience functions for computing the (log) determinant of the
matrix that has been factored.
* A convenience function for explicitly computing the inverse of the
matrix that has been factored (though this is rarely useful).
Quickstart
----------
If :math:`A` is a sparse, symmetric, positive-definite matrix, and
:math:`b` is a matrix or vector (either sparse or dense), then the
following code solves the equation :math:`Ax = b`::
from sksparse.cholmod import cholesky
factor = cholesky(A)
x = factor(b)
If we just want to compute its determinant::
factor = cholesky(A)
ld = factor.logdet()
(This returns the log of the determinant, rather than the determinant
itself, to avoid issues with underflow/overflow. See :meth:`logdet`,
:meth:`log`.)
If you have a least-squares problem to solve, minimizing :math:`||Mx -
b||^2`, and :math:`M` is a sparse matrix, the `solution
`_
is :math:`x = (M'M)^{-1} M'b`, which can be efficiently calculated
as::
from sksparse.cholmod import cholesky_AAt
# Notice that CHOLMOD computes AA' and we want M'M, so we must set A = M'!
factor = cholesky_AAt(M.T)
x = factor(M.T * b)
However, you should be aware that for least squares problems, the
Cholesky method is usually faster but somewhat less numerically stable
than QR- or SVD-based techniques.
Top-level functions
-------------------
All usage of this module starts by calling one of four functions, all
of which return a :class:`Factor` object, documented below.
Most users will want one of the ``cholesky`` functions, which perform
a fill-reduction analysis and decomposition together:
.. autofunction:: cholesky(A, beta=0, mode="auto", ordering_method="default", use_long=None)
.. autofunction:: cholesky_AAt(A, beta=0, mode="auto", ordering_method="default", use_long=None)
However, some users may want to break the fill-reduction analysis and
actual decomposition into separate steps, and instead begin with one
of the ``analyze`` functions, which perform only fill-reduction:
.. autofunction:: analyze(A, mode="auto", ordering_method="default", use_long=None)
.. autofunction:: analyze_AAt(A, mode="auto", ordering_method="default", use_long=None)
.. note:: Even if you used :func:`cholesky` or :func:`cholesky_AAt`,
you can still call :meth:`cholesky_inplace()
` or :meth:`cholesky_AAt_inplace()
` on the resulting :class:`Factor` to
quickly factor another matrix with the same non-zero pattern as your
original matrix.
:class:`Factor` objects
-----------------------
.. class:: Factor
A :class:`Factor` object represents the Cholesky decomposition of some
matrix :math:`A` (or :math:`AA'`). Each :class:`Factor` fixes:
* A specific fill-reducing permutation
* A choice of which Cholesky algorithm to use (see :func:`analyze`)
* Whether we are currently working with real numbers or complex
Given a :class:`Factor` object, you can:
* Compute new Cholesky decompositions of matrices that have the same
pattern of non-zeros
* Perform 'updates' or 'downdates'
* Access the various Cholesky factors
* Solve equations involving those factors
Factoring new matrices
++++++++++++++++++++++
.. automethod:: Factor.cholesky_inplace(A, beta=0)
.. automethod:: Factor.cholesky_AAt_inplace(A, beta=0)
.. automethod:: Factor.cholesky(A, beta=0)
.. automethod:: Factor.cholesky_AAt(A, beta=0)
Updating/Downdating
+++++++++++++++++++
.. automethod:: Factor.update_inplace(C, subtract=False)
Accessing Cholesky factors explicitly
+++++++++++++++++++++++++++++++++++++
.. note:: When possible, it is generally more efficient to use the
``solve_...`` functions documented below rather than extracting the
Cholesky factors explicitly.
.. automethod:: Factor.P
.. automethod:: Factor.D
.. automethod:: Factor.L
.. automethod:: Factor.LD
.. automethod:: Factor.L_D
Solving equations
+++++++++++++++++
All methods in this section accept both sparse and dense matrices (or
vectors) ``b``, and return either a sparse or dense ``x``
accordingly.
All methods in this section act on :math:`LDL'` factorizations by default.
Thus `L` refers by default to the matrix returned by :meth:`L_D`, not that
returned by :meth:`L` (though conversion is not performed unless necessary).
.. automethod:: Factor.solve_A(b)
.. automethod:: Factor.__call__(b)
.. automethod:: Factor.solve_LDLt(b)
.. automethod:: Factor.solve_LD(b)
.. automethod:: Factor.solve_DLt(b)
.. automethod:: Factor.solve_L(b)
.. automethod:: Factor.solve_Lt(b)
.. automethod:: Factor.solve_D(b)
.. automethod:: Factor.apply_P(b)
.. automethod:: Factor.apply_Pt(b)
Convenience methods
-------------------
.. automethod:: Factor.logdet()
.. automethod:: Factor.det()
.. automethod:: Factor.slogdet()
.. automethod:: Factor.inv()
.. automethod:: Factor.copy()
Error handling
--------------
.. class:: CholmodError
.. class:: CholmodNotPositiveDefiniteError
.. class:: CholmodNotInstalledError
.. class:: CholmodOutOfMemoryError
.. class:: CholmodTooLargeError
.. class:: CholmodNotPositiveDefiniteError
.. class:: CholmodInvalidError
.. class:: CholmodGpuProblemError
Errors detected by CHOLMOD or by our wrapper code are converted into
exceptions of type :class:`CholmodError` or an appropriated subclass.
.. class:: CholmodWarning
Warnings issued by CHOLMOD are converted into Python warnings of
type :class:`CholmodWarning`.
.. class:: CholmodTypeConversionWarning
CHOLMOD itself supports matrices in CSC form with 32-bit integer
indices and 'double' precision floats (64-bits, or 128-bits total
for complex numbers). If you pass some other sort of matrix, then
the wrapper code will convert it for you before passing it to
CHOLMOD, and issue a warning of type
:class:`CholmodTypeConversionWarning` to let you know that your
efficiency is not as high as it might be.
.. warning:: Not all conversions currently produce warnings. This is
a bug.
Child of :class:`CholmodWarning`.