6.3 Positive Definite Linear Equations (cvxopt.cholmod)

cvxopt.cholmod is an interface to the Cholesky factorization routines of the CHOLMOD package. It includes functions for Cholesky factorization of sparse positive definite matrices, and for solving sparse sets of linear equations with positive definite matrices. The routines can also be used for computing LDLT (or LDLH) factorizations of symmetric indefinite matrices (with L unit lower-triangular and D diagonal and nonsingular) if such a factorization exists.

See also: CHOLMOD code, documentation, copyright and license.

linsolve(A, B[, p=None[, uplo=’L’]])

Solves

AX = B

with A sparse and real symmetric or complex Hermitian. B is a dense matrix of the same type as A. On exit it is overwritten with the solution. The argument p is an integer matrix with length equal to the order of A, and specifies an optional reordering of A. If p is not specified, CHOLMOD used a reordering from the AMD library. Raises an ArithmeticError if the factorization does not exist.

As an example, we solve

⌊                ⌋    ⌊      ⌋
  10    0 3    0         0  4
||  0    5 0  - 2 ||X = ||  1  5|| .
⌈  3    0 5    0 ⌉    ⌈  2  6⌉
   0  - 2 0    2         3  7
(6.2)

>>> from cvxopt import matrix, spmatrix, cholmod  
>>> A = spmatrix([10,3, 5,-2, 5, 2], [0,2, 1,3, 2, 3], [0,0, 1,1, 2, 3])  
>>> X = matrix(range(8), (4,2), ’d’)  
>>> cholmod.linsolve(A,X)  
>>> print X  
[-1.46e-01  4.88e-02]  
[ 1.33e+00  4.00e+00]  
[ 4.88e-01  1.17e+00]  
[ 2.83e+00  7.50e+00]

splinsolve(A, B[, p=None[, uplo=’L’]])

Similar to linsolve() except that B is a sparse matrix and that the solution is returned as an output argument (as a new sparse matrix). B is not modified.

The following code computes the inverse of the coefficient matrix in (6.2) as a sparse matrix.

>>> X = cholmod.splinsolve(A, spmatrix(1.0,range(4),range(4)))  
>>> print X  
[ 1.22e-01     0     -7.32e-02     0    ]  
[    0      3.33e-01     0      3.33e-01]  
[-7.32e-02     0      2.44e-01     0    ]  
[    0      3.33e-01     0      8.33e-01]

The functions linsolve() and splinsolve() are equivalent to symbolic() and numeric() called in sequence, followed by solve(), respectively, spsolve().

symbolic(A[, p=None[, uplo=’L’]])

Performs a symbolic analysis of a sparse real symmetric or complex Hermitian matrix A for one of the two factorizations:

PAP T = LLT ,   P AP T = LLH,
(6.3)

and

P AP T = LDLT ,   P AP T = LDLH,
(6.4)

where P is a permutation matrix, L is lower triangular (unit lower triangular in the second factorization), and D is nonsingular diagonal. The type of factorization depends on the value of options[’supernodal’] (see below).

If uplo is ’L’, only the lower triangular part of A is accessed and the upper triangular part is ignored. If uplo is ’U’, only the upper triangular part of A is accessed and the lower triangular part is ignored.

If p is not provided, a reordering from the AMD library is used.

The symbolic factorization is returned as an opaque C object that can be passed to cholmod.numeric().

numeric(A, F)

Performs a numeric factorization of a sparse symmetric matrix as (6.3) or (6.4). The argument F is the symbolic factorization computed by cholmod.symbolic() applied to the matrix A, or to another sparse matrix with the same sparsity pattern and typecode, or by cholmod.numeric() applied to a matrix with the same sparsity pattern and typecode as A.

If F was created by a cholmod.symbolic with uplo equal to ’L’, then only the lower triangular part of A is accessed and the upper triangular part is ignored. If it was created with uplo is ’U’, then only the upper triangular part of A is accessed and the lower triangular part is ignored.

On successful exit, the factorization is stored in F. Raises an ArithmeticError if the factorization does not exist.

solve(F, B[, sys=0])

Solves one of the following linear equations where B is a dense matrix and F is the numeric factorization (6.3) or (6.4) computed by cholmod_numeric(). sys is an integer with values between 0 and 8.



sys equation


0 AX = B


1 LDLTX = B


2 LDLX = B


3 DLTX = B


4 LX = B


5 LTX = B


6 DX = B


7 PTX = B


8 PX = B


(If F is a Cholesky factorization of the form (6.3), D is an identity matrix in this table. If A is complex, LT should be replaced by LH.)

The matrix B is a dense ’d’ or ’z’ matrix, with the same type as A. On exit it is overwritten by the solution.

spsolve(F, B[, sys=0])

Similar to solve(), except that B is a sparse matrix, and the solution is returned as an output argument (as a sparse matrix). B must have the same typecode as A.

For the same example as above:

>>> X = matrix(range(8), (4,2), ’d’)  
>>> F = cholmod.symbolic(A)  
>>> cholmod.numeric(A, F)  
>>> cholmod.solve(F, X)  
>>> print X  
[-1.46e-01  4.88e-02]  
[ 1.33e+00  4.00e+00]  
[ 4.88e-01  1.17e+00]  
[ 2.83e+00  7.50e+00]

diag(F)

Returns the diagonal elements of the Cholesky factor L in (6.3), as a dense matrix of the same type as A. Note that this only applies to Cholesky factorizations. The matrix D in an LDLT factorization can be retrieved via cholmod.solve() with sys equal to 6.

In the functions listed above, the default values of the control parameters described in the CHOLMOD user guide are used, except for Common->print which is set to 0 instead of 3 and Common->supernodal which is set to 2 instead of 1. These parameters (and a few others) can be modified by making an entry in the dictionary cholmod.options. The meaning of these parameters is as follows.

options[’supernodal’]
If equal to 0, a factorization (6.4) is computed using a simplicial algorithm. If equal to 2, a factorization (6.3) is computed using a supernodal algorithm. If equal to 1, the most efficient of the two factorizations is selected, based on the sparsity pattern. Default: 2.
options[’print’]
A nonnegative integer that controls the amount of output printed to the screen. Default: 0 (no output).

As an example that illustrates diag() and the use of cholmod.options, we compute the logarithm of the determinant of the coefficient matrix in (6.2) by two methods.

>>> import math  
>>> from cvxopt.cholmod import options  
>>> from cvxopt import log  
>>> F = cholmod.symbolic(A)  
>>> cholmod.numeric(A, F)  
>>> print 2.0 * sum(log(cholmod.diag(F)))  
5.50533153593

>>> options[’supernodal’] = 0  
>>> F = cholmod.symbolic(A)  
>>> cholmod.numeric(A, F)  
>>> Di = matrix(1.0, (4,1))  
>>> cholmod.solve(F, Di, sys=6)  
>>> print -sum(log(Di))  
5.50533153593