MATLAB Functions

Incomplete Cholesky factorizations

cholinc(X,'0')
R = cholinc(X,'0')
[R,p] = cholinc(X,'0')
R = cholinc(X,droptol)
R = cholinc(X,options)

Description

cholinc(X,'0')

produces the incomplete Cholesky factorization of a real symmetric positive definite sparse matrix with 0 level of fill-in. cholinc(X,'0') produces an upper triangular matrix. The lower triangle of X is assumed to be the transpose of the upper (X is symmetric).

R = cholinc(X,'0')

returns an upper triangular matrix which has the same sparsity pattern as the upper triangle of X. The product R'*R agrees with X over its sparsity pattern. The positive definiteness of X is not sufficient to guarantee the existence of the incomplete factor, and, in this case, an error message is printed.

[R,p] = cholinc(X,'0')

never produces an error message. If the incomplete factor exists, then p = 0 and R is the upper triangular factor. If the calculation of R breaks down due to a zero or negative pivot, then p is a positive integer and R is an upper triangular matrix of size q-by-n where q = p-1. The sparsity pattern of R is the same as the q-by-n upper triangle of X and the n-by-n product R'*R agrees with X over the sparsity pattern of its first q rows and columns X(1:q,:) and X(:,1:q).

R = cholinc(X,droptol)

computes the incomplete Cholesky factorization of any sparse matrix using a drop tolerance. droptol must be a non-negative scalar. cholinc(X,droptol) produces an approximation to the complete Cholesky factor returned by chol(X). For increasingly smaller values of the drop tolerance, this approximation improves, until the drop tolerance is 0, at which time the complete Cholesky factorization is produced, as in chol(X).

The off-diagonal entries R(i,j) which are smaller in magnitude than the local drop tolerance, which is given by droptol*norm(X(:,j))/R(i,i), are dropped from the factor. The diagonal entries are preserved even if they are too small in an attempt to avoid a singular factor.

R = cholinc(X,options)

specifies a structure with up to three fields which may be used in any combination: droptol, michol, rdiag. Additional fields are ignored.

droptol is the drop tolerance of the incomplete factorization.

If michol is 1, cholinc produces the modified incomplete Cholesky factorization which subtracts the dropped elements in any column from the diagonal element of the upper triangular factor. The default value is 0.

If rdiag is 1, any zeros on the diagonal of the upper triangular factor are replaced by the square root of the product of the drop tolerance and the norm of that column of X, sqrt(droptol*norm(X(:,j))). The default is 0. Note that the thresh option available in the incomplete LU factorization (see luinc) is not here as it is always set to 0.There are never any row interchanges during the formation of the incomplete Cholesky factor.

R = cholinc(X,droptol) and R = cholinc(X,options)

return an upper triangular matrix in R. The product R'*R is an approximation to the matrix X.

Remarks

These incomplete factorizations may be useful as preconditioners for solving large sparse systems of linear equations. A single 0 on the diagonal of the upper triangular factor makes it singular. The incomplete factorization with a drop tolerance prints a warning message if the upper triangular factor has zeros on the diagonal. Similarly, using the rdiag option to replace a zero diagonal only gets rid of the symptoms of the problem, but it does not solve it. The preconditioner may not be singular, but it probably is not useful, and a warning message is printed.

Examples

Start with a symmetric positive definite matrix, S.

S = delsq(numgrid('C',15));

S is the two-dimensional, five-point discrete negative Lapacian on the grid generated by numgrid('C',15).

Compute the Cholesky factorization and the incomplete Cholesky factorization of level 0 to compare the fill-in. Make S singular by zeroing out a diagonal entry and compute the (partial) incomplete Cholesky factorization of level 0.

C = chol(S);
R0 = chol(S,'0');
S2 = S; S2(101,101) = 0;
[R,p] = cholinc(S2,'0');

There is fill-in within the bands of S in the complete Cholesky factor, but none in the incomplete Cholesky factor. The incomplete factorization of the singular S2 stopped at row p = 101 resulting in a 100-by-139 partial factor.

D1 = (R0'*R0).*spones(S)-S;
D2 = (R'*R).*spones(S2)-S2;

D1 has elements of the order of eps, showing that R0'*R0 agrees with S over its sparsity pattern. D2 has elements of the order of eps over its first 100 rows and first 100 columns, D2(1:100,:) and D2(:,1:100).

The first subplot below shows that cholinc(S,0), the incomplete Cholesky factor with a drop tolerance of 0, is the same as the Cholesky factor of S. Increasing the drop tolerance increases the sparsity of the incomplete factors, as seen below.

Unfortunately, the sparser factors are poor approximations, as is seen by the plot of drop tolerance versus norm(R'*R-S,1)/norm(S,1) in the next figure.

Limitations

cholinc works on square matrices only. For cholinc(X,'0'), X must be real.

Algorithm

R = cholinc(X,droptol) is obtained from [L,U] = luinc(X,options), where options.droptol = droptol and options.thresh = 0. The rows of the uppertriangular U are scaled by the square root of the diagonal in that row, and this scaled factor becomes R.

R = cholinc(X,options) is produced in a similar manner, except the rdiag option translates into the udiag option and the milu option takes the value of the michol option.

cholinc(X,'0') is based on the "KJI" variant of the Cholesky factorization. Updates are made only to positions which are nonzero in the upper triangle of X.