Navigation

Operators and Keywords

Function List:

C++ API

: L = ichol (A)
: L = ichol (A, opts)

Compute the incomplete Cholesky factorization of the sparse square matrix A.

By default, ichol uses only the lower triangle of A and produces a lower triangular factor L such that L*L' approximates A.

The factor given by this routine may be useful as a preconditioner for a system of linear equations being solved by iterative methods such as PCG (Preconditioned Conjugate Gradient).

The factorization may be modified by passing options in a structure opts. The option name is a field of the structure and the setting is the value of field. Names and specifiers are case sensitive.

type

Type of factorization.

"nofill" (default)

Incomplete Cholesky factorization with no fill-in (IC(0)).

"ict"

Incomplete Cholesky factorization with threshold dropping (ICT).

diagcomp

A non-negative scalar alpha for incomplete Cholesky factorization of A + alpha * diag (diag (A)) instead of A. This can be useful when A is not positive definite. The default value is 0.

droptol

A non-negative scalar specifying the drop tolerance for factorization if performing ICT. The default value is 0 which produces the complete Cholesky factorization.

Non-diagonal entries of L are set to 0 unless

abs (L(i,j)) >= droptol * norm (A(j:end, j), 1).

michol

Modified incomplete Cholesky factorization:

"off" (default)

Row and column sums are not necessarily preserved.

"on"

The diagonal of L is modified so that row (and column) sums are preserved even when elements have been dropped during the factorization. The relationship preserved is: A * e = L * L' * e, where e is a vector of ones.

shape
"lower" (default)

Use only the lower triangle of A and return a lower triangular factor L such that L*L' approximates A.

"upper"

Use only the upper triangle of A and return an upper triangular factor U such that U'*U approximates A.

EXAMPLES

The following problem demonstrates how to factorize a sample symmetric positive definite matrix with the full Cholesky decomposition and with the incomplete one. 

A = [ 0.37, -0.05,  -0.05,  -0.07;
     -0.05,  0.116,  0.0,   -0.05;
     -0.05,  0.0,    0.116, -0.05;
     -0.07, -0.05,  -0.05,   0.202];
A = sparse (A);
nnz (tril (A))
ans =  9
L = chol (A, "lower");
nnz (L)
ans =  10
norm (A - L * L', "fro") / norm (A, "fro")
ans =  1.1993e-16
opts.type = "nofill";
L = ichol (A, opts);
nnz (L)
ans =  9
norm (A - L * L', "fro") / norm (A, "fro")
ans =  0.019736

Another example for decomposition is a finite difference matrix used to solve a boundary value problem on the unit square. 

nx = 400; ny = 200;
hx = 1 / (nx + 1); hy = 1 / (ny + 1);
Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)],
               [-1 0 1 ], nx, nx) / (hx ^ 2);
Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)],
               [-1 0 1 ], ny, ny) / (hy ^ 2);
A = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy);
nnz (tril (A))
ans =  239400
opts.type = "nofill";
L = ichol (A, opts);
nnz (tril (A))
ans =  239400
norm (A - L * L', "fro") / norm (A, "fro")
ans =  0.062327

References for implemented algorithms:

[1] Y. Saad. "Preconditioning Techniques." Iterative Methods for Sparse Linear Systems, PWS Publishing Company, 1996.

[2] M. Jones, P. Plassmann: An Improved Incomplete Cholesky Factorization, 1992.

See also: chol, ilu, pcg.

Package: octave