Incomplete Cholesky factorization
L = ichol(A)
L = ichol(A,options)
L = ichol( performs the incomplete Cholesky
A with zero-fill.
ichol references the lower triangle of
and produces lower triangular factors.
Structure with up to five fields:
Incomplete Cholesky Factorization
This example generates an incomplete Cholesky factorization.
Start with a symmetric positive definite matrix,
N = 100; A = delsq(numgrid('S',N));
A is the two-dimensional, five-point discrete negative Laplacian on a 100-by-100 square grid with Dirichlet boundary conditions. The size of
98*98 = 9604 (not 10000 as the borders of the grid are used to impose the Dirichlet conditions).
The no-fill incomplete Cholesky factorization is a factorization which contains only nonzeros in the same position as
A contains nonzeros. This factorization is extremely cheap to compute. Although the product
L*L' is typically very different from
A, the product
L*L' will match
A on its pattern up to round-off.
L = ichol(A); norm(A-L*L','fro')./norm(A,'fro')
ans = 0.0916
ans = 4.9606e-17
ichol may also be used to generate incomplete Cholesky factorizations with threshold dropping. As the drop tolerance decreases, the factor tends to get more dense and the product
L*L' tends to be a better approximation of
A. The following plots show the relative error of the incomplete factorization plotted against the drop tolerance as well as the ratio of the density of the incomplete factors to the density of the complete Cholesky factor.
n = size(A,1); ntols = 20; droptol = logspace(-8,0,ntols); nrm = zeros(1,ntols); nz = zeros(1,ntols); nzComplete = nnz(chol(A,'lower')); for k = 1:ntols L = ichol(A,struct('type','ict','droptol',droptol(k))); nz(k) = nnz(L); nrm(k) = norm(A-L*L','fro')./norm(A,'fro'); end figure loglog(droptol,nrm,'LineWidth',2) title('Drop tolerance vs norm(A-L*L'',''fro'')./norm(A,''fro'')')
figure semilogx(droptol,nz./nzComplete,'LineWidth',2) title('Drop tolerance vs fill ratio ichol/chol')
The relative error is typically on the same order as the drop tolerance, although this is not guaranteed.
Using ichol as a Preconditioner
This example shows how to use an incomplete Cholesky factorization as a preconditioner to improve convergence.
Create a symmetric positive definite matrix,
N = 100; A = delsq(numgrid('S',N));
Create an incomplete Cholesky factorization as a preconditioner for
pcg. Use a constant vector as the right hand side. As a baseline, execute
pcg without a preconditioner.
b = ones(size(A,1),1); tol = 1e-6; maxit = 100; [x0,fl0,rr0,it0,rv0] = pcg(A,b,tol,maxit);
fl0 = 1 indicating that
pcg did not drive the relative residual to the requested tolerance in the maximum allowed iterations. Try the no-fill incomplete Cholesky factorization as a preconditioner.
L1 = ichol(A); [x1,fl1,rr1,it1,rv1] = pcg(A,b,tol,maxit,L1,L1');
fl1 = 0, indicating that
pcg converged to the requested tolerance and did so in 59 iterations (the value of
it1). Since this matrix is a discretized Laplacian, however, using modified incomplete Cholesky can create a better preconditioner. A modified incomplete Cholesky factorization constructs an approximate factorization that preserves the action of the operator on the constant vector. That is,
norm(A*e-L*(L'*e)) will be approximately zero for
e = ones(size(A,2),1) even though
norm(A-L*L','fro')/norm(A,'fro') is not close to zero. It is not necessary to specify type for this syntax since
nofill is the default, but it is good practice.
options.type = 'nofill'; options.michol = 'on'; L2 = ichol(A,options); e = ones(size(A,2),1); norm(A*e-L2*(L2'*e))
ans = 3.7983e-14
[x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,L2,L2');
pcg converges (
fl2 = 0) but in only 38 iterations. Plotting all three convergence histories shows the convergence.
semilogy(0:maxit,rv0./norm(b),'b.'); hold on semilogy(0:it1,rv1./norm(b),'r.'); semilogy(0:it2,rv2./norm(b),'k.'); legend('No Preconditioner','IC(0)','MIC(0)');
The plot shows that the modified incomplete Cholesky preconditioner creates a much faster convergence.
You can also try incomplete Cholesky factorizations with threshold dropping. The following plot shows convergence of
pcg with preconditioners constructed with various drop tolerances.
L3 = ichol(A, struct('type','ict','droptol',1e-1)); [x3,fl3,rr3,it3,rv3] = pcg(A,b,tol,maxit,L3,L3'); L4 = ichol(A, struct('type','ict','droptol',1e-2)); [x4,fl4,rr4,it4,rv4] = pcg(A,b,tol,maxit,L4,L4'); L5 = ichol(A, struct('type','ict','droptol',1e-3)); [x5,fl5,rr5,it5,rv5] = pcg(A,b,tol,maxit,L5,L5'); figure semilogy(0:maxit,rv0./norm(b),'b-','linewidth',2); hold on semilogy(0:it3,rv3./norm(b),'b-.','linewidth',2); semilogy(0:it4,rv4./norm(b),'b--','linewidth',2); semilogy(0:it5,rv5./norm(b),'b:','linewidth',2); legend('No Preconditioner','ICT(1e-1)','ICT(1e-2)', ... 'ICT(1e-3)','Location','SouthEast');
Note the incomplete Cholesky preconditioner constructed with drop tolerance
1e-2 is denoted as
As with the zero-fill incomplete Cholesky, the threshold dropping factorization can benefit from modification (i.e.
options.michol = 'on') since the matrix arises from an elliptic partial differential equation. As with
MIC(0), the modified threshold based dropping incomplete Cholesky will preserve the action of the preconditioner on constant vectors, that is
norm(A*e-L*(L'*e)) will be approximately zero.
Using the diagcomp Option
This example illustrates the use of the
diagcomp option of
Incomplete Cholesky factorizations of positive definite matrices do not always exist. The following code constructs a random symmetric positive definite matrix and attempts to solve a linear system using
S = rng('default'); A = sprandsym(1000,1e-2,1e-4,1); rng(S); b = full(sum(A,2)); [x0,fl0,rr0,it0,rv0] = pcg(A,b,1e-6,100);
Since convergence is not attained, try to construct an incomplete Cholesky preconditioner.
L = ichol(A,struct('type','ict','droptol',1e-3));
Error using ichol Encountered nonpositive pivot.
ichol breaks down as above, you can use the
diagcomp option to construct a shifted incomplete Cholesky factorization. That is, instead of constructing
L such that
ichol with diagonal compensation constructs
L such that
M = A + alpha*diag(diag(A)) without explicitly forming
M. As incomplete factorizations always exist for diagonally dominant matrices,
alpha can be found to make
M diagonally dominant.
alpha = max(sum(abs(A),2)./diag(A))-2
alpha = 62.9341
L1 = ichol(A, struct('type','ict','droptol',1e-3,'diagcomp',alpha)); [x1,fl1,rr1,it1,rv1] = pcg(A,b,1e-6,100,L1,L1');
pcg still fails to converge to the desired tolerance within the desired number of iterations, but as the plot below shows, convergence is better for
pcg with this preconditioner than with no preconditioner. Choosing a smaller
alpha may help. With some experimentation, we can settle on an appropriate value for
alpha = .1; L2 = ichol(A, struct('type','ict','droptol',1e-3,'diagcomp',alpha)); [x2,fl2,rr2,it2,rv2] = pcg(A,b,1e-6,100,L2,L2');
pcg converges and a plot can show the convergence histories with each preconditioner.
semilogy(0:100,rv0./norm(b),'b.'); hold on; semilogy(0:100,rv1./norm(b),'r.'); semilogy(0:it2,rv2./norm(b),'k.'); legend('No Preconditioner','\alpha \approx 63','\alpha = .1'); xlabel('Iteration Number'); ylabel('Relative Residual');
 Saad, Yousef. “Preconditioning Techniques.” Iterative Methods for Sparse Linear Systems. PWS Publishing Company, 1996.
 Manteuffel, T.A. “An incomplete factorization technique for positive definite linear systems.” Math. Comput. 34, 473–497, 1980.
Run code in the background using MATLAB®
backgroundPool or accelerate code with Parallel Computing Toolbox™
This function fully supports thread-based environments. For more information, see Run MATLAB Functions in Thread-Based Environment.
Partition large arrays across the combined memory of your cluster using Parallel Computing Toolbox™.
Usage notes and limitations:
If you specify the
options, it must be
Amust be real symmetric or complex hermitian.
For more information, see Run MATLAB Functions with Distributed Arrays (Parallel Computing Toolbox).