Sparse Cholesky decomposition (sksparse.cholmod
)¶
New in version 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 \(LL' = A\) or \(LDL' = A\) (with fill-reducing permutation) for both real and complex sparse matrices \(A\), in any format supported by
scipy.sparse
. (However, CSC matrices will be most efficient.)A convenient and efficient interface for using this decomposition to solve problems of the form \(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 \(A\) and of product \(AA'\). So, the result is the Cholesky decomposition of \(A + CC'\) (or \(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 \(A\) is a sparse, symmetric, positive-definite matrix, and \(b\) is a matrix or vector (either sparse or dense), then the following code solves the equation \(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 logdet()
,
log()
.)
If you have a least-squares problem to solve, minimizing \(||Mx - b||^2\), and \(M\) is a sparse matrix, the solution is \(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 Factor
object, documented below.
Most users will want one of the cholesky
functions, which perform
a fill-reduction analysis and decomposition together:
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:
Note
Even if you used cholesky()
or cholesky_AAt()
,
you can still call cholesky_inplace()
or cholesky_AAt_inplace()
on the resulting Factor
to
quickly factor another matrix with the same non-zero pattern as your
original matrix.
Factor
objects¶
-
class
sksparse.cholmod.
Factor
¶ A
Factor
object represents the Cholesky decomposition of some matrix \(A\) (or \(AA'\)). EachFactor
fixes:A specific fill-reducing permutation
A choice of which Cholesky algorithm to use (see
analyze()
)Whether we are currently working with real numbers or complex
Given a
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¶
Updating/Downdating¶
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.
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 \(LDL'\) factorizations by default.
Thus L refers by default to the matrix returned by L_D()
, not that
returned by L()
(though conversion is not performed unless necessary).
Convenience methods¶
Error handling¶
-
class
sksparse.cholmod.
CholmodError
¶
-
class
sksparse.cholmod.
CholmodNotPositiveDefiniteError
¶
-
class
sksparse.cholmod.
CholmodNotInstalledError
¶
-
class
sksparse.cholmod.
CholmodOutOfMemoryError
¶
-
class
sksparse.cholmod.
CholmodTooLargeError
¶
-
class
sksparse.cholmod.
CholmodNotPositiveDefiniteError
-
class
sksparse.cholmod.
CholmodInvalidError
¶
-
class
sksparse.cholmod.
CholmodGpuProblemError
¶ Errors detected by CHOLMOD or by our wrapper code are converted into exceptions of type
CholmodError
or an appropriated subclass.
-
class
sksparse.cholmod.
CholmodWarning
¶ Warnings issued by CHOLMOD are converted into Python warnings of type
CholmodWarning
.
-
class
sksparse.cholmod.
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
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
CholmodWarning
.