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'\)). Each Factor 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.