| MATLAB Function Reference | Search Help Desk |
| cholinc | Examples See Also |
Sparse Incomplete Cholesky and Cholesky-Infinity factorizations
R = cholinc(X,droptol) R = cholinc(X,options) R = cholinc(X,'0') [R,p] = cholinc(X,'0') R = cholinc(X,'inf') [R,p] = cholinc(X,'inf')
cholinc produces two different kinds of incomplete Cholesky factorizations: the drop tolerance and the 0 level of fill-in factorizations. These factors may be useful as preconditioners for a symmetric positive definite system of linear equations being solved by an iterative method such as pcg (Preconditioned Conjugate Gradients). cholinc works only for sparse matrices.
R = cholinc(X,droptol)
performs the incomplete Cholesky factorization of X, with drop tolerance droptol.
R = cholinc(X,options)
allows additional options to the incomplete Cholesky factorization. options is a structure with up to three fields:droptol |
drop tolerance of the incomplete factorization |
michol |
modified incomplete Cholesky |
rdiag |
replace zeros on the diagonal of R |
droptol is a non-negative scalar used as the drop tolerance for the incomplete Cholesky factorization. This factorization is computed by performing the incomplete LU factorization with the pivot threshold option set to 0 (which forces diagonal pivoting) and then scaling the rows of the incomplete upper triangular factor, U, by the square root of the diagonal entries in that column. Since the nonzero entries U(i,j) are bounded below by droptol*norm(X(:,j)) (see luinc), the nonzero entries R(i,j) are bounded below by the local drop tolerance droptol*norm(X(:,j))/R(i,i).
Setting droptol = 0 produces the complete Cholesky factorization, which is the default.
michol stands for modified incomplete Cholesky factorization. Its value is either 0 (unmodified, the default) or 1 (modified). This performs the modified incomplete LU factorization of X and scales the returned upper triangular factor as described above.
rdiag is either 0 or 1. If it is 1, any zero diagonal entries of the upper triangular factor R are replaced by the square root of the local drop tolerance in an attempt to avoid a singular factor. The default is 0.
R = cholinc(X,'0')
produces the incomplete Cholesky factor of a real symmetric positive definite sparse matrix with 0 level of fill-in (i.e. no fill-in). The upper triangular R has the same sparsity pattern as triu(X), although R may be zero in some positions where X is nonzero due to cancellation. The lower triangle of X is assumed to be the transpose of the upper. Note that the positive definiteness of X does not guarantee the existence of a factor with the required sparsity. An error message results if the factorization is not possible. If the factorization is successful, R'*R agrees with X over its sparsity pattern.
[R,p] = cholinc(X,'0')
with two output arguments, never produces an error message. If R exists, p is 0. But, if the incomplete factor does not exist, then p is a positive integer and R is an upper triangular matrix of size q-by-n where q = p-1 so that the sparsity pattern of R is that of the q-by-n upper triangle of X. R'*R agrees with X over the sparsity pattern of its first q rows and first q columns.
R = cholinc(X,'inf')
produces the Cholesky-Infinity factorization. This factorization is based on the Cholesky factorization, and additionally handles real positive semi-definite matrices. It may be useful for finding a solution to systems which arise in interior-point methods. When a zero pivot is encountered in the ordinary Cholesky factorization, the diagonal of the Cholesky-Infinity factor is set to Inf and the rest of that row is set to 0. This is designed to force a 0 in the corresponding entry of the solution vector in the associated system of linear equations. A negative pivot still results in an error.
[R,p] = cholinc(X,'inf')
with two output arguments, never produces an error message. If X is positive semi-definite, then p is 0. Otherwise, p is a positive integer and R is an upper triangular matrix of size(p-1)-by-n.
1. The 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.
2. The Cholesky-Infinity factorization is meant to be used within interior-point methods. Its use otherwise cannot be recommended.
Example 1.
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 = cholinc(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).
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.
norm(R'*R-S,1)/norm(S,1) in the next figure.
S = sparse([1 0 3 0;0 25 0 30;3 0 9 0;0 30 0 661])This symmetric sparse matrix is singular. The Cholesky factorization fails at the zero pivot in the third row, but
cholinc succeeds in computing all rows of the Cholesky-Infinity factorization.
[R,p] = chol(S);
full(R)
ans =
1 0 3 0
0 5 0 6
Rinf = cholinc(S,'inf');
full(Rinf)
ans =
1 0 3 0
0 5 0 6
0 0 Inf 0
0 0 0 25
cholinc works on square sparse matrices only. For cholinc(X,'0') and cholinc(X,'inf'), X must be real.
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.
R = 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.
R = cholinc(X,'inf') is based on the algorithm in Zhang.
chol Cholesky factorization
luinc Incomplete LU matrix factorizations
pcg Preconditioned Conjugate Gradients method