Sparse Cholesky decomposition (scikits.sparse.cholmod)

New in version 0.1.

Overview

This module exposes most of the capabilities of the CHOLMOD package. Specifically, it provides:

  • 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 product AA' 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.

The most common use is probably for solving least squares problems.

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 scikits.sparse.cholmod import cholesky
factor = cholesky(A)
x = factor(b)

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 scikits.sparse.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)

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 scikits.sparse.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; L always refers to the matrix returned by L_D(), not that returned by L() (though conversion is not performed unless necessary).

Note

If you need an efficient implementation of solve_L() or solve_Lt() that works with the LL' factorization, then drop us a line, it’d be easy to add.

Error handling

class scikits.sparse.cholmod.CholmodError

Errors detected by CHOLMOD or by our wrapper code are converted into exceptions of type CholmodError.

class scikits.sparse.cholmod.CholmodWarning

Warnings issued by CHOLMOD are converted into Python warnings of type CholmodWarning.

class scikits.sparse.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.

Citing CHOLMOD

Tim Davies, the author of CHOLMOD, asks that if you use CHOLMOD in the production of published scientific research, you cite one or more of the following papers:

  • Dynamic supernodes in sparse Cholesky update/downdate and triangular solves, T. A. Davis and W. W. Hager, ACM Trans. Math. Software, Vol 35, No. 4, 2009.
  • Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate , Y. Chen, T. A. Davis, W. W. Hager, and S. Rajamanickam, ACM Trans. Math. Software, Vol 35, No. 3, 2009.
  • Row modifications of a sparse Cholesky factorization, T. A. Davis and W. W. Hager, SIAM Journal on Matrix Analysis and Applications, vol 26, no 3, pp. 621-639, 2005.
  • Multiple-rank modifications of a sparse Cholesky factorization, T. A. Davis and W. W. Hager, SIAM Journal on Matrix Analysis and Applications, vol. 22, no. 4, pp. 997-1013, 2001.
  • Modifying a sparse Cholesky factorization, T. A. Davis and W. W. Hager, SIAM Journal on Matrix Analysis and Applications, vol. 20, no. 3, pp. 606-627, 1999.