(octave.info)Iterative Techniques
22.3 Iterative Techniques Applied to Sparse Matrices
====================================================
The left division ‘\’ and right division ‘/’ operators, discussed in the
previous section, use direct solvers to resolve a linear equation of the
form ‘X = A \ B’ or ‘X = B / A’. Octave also includes a number of
functions to solve sparse linear equations using iterative techniques.
-- : X = pcg (A, B, TOL, MAXIT, M1, M2, X0, ...)
-- : X = pcg (A, B, TOL, MAXIT, M, [], X0, ...)
-- : [X, FLAG, RELRES, ITER, RESVEC, EIGEST] = pcg (A, B, ...)
Solve the linear system of equations ‘A * X = B’ by means of the
Preconditioned Conjugate Gradient iterative method.
The input arguments are:
• A is the matrix of the linear system and it must be square. A
can be passed as a matrix, function handle, or inline function
‘Afun’ such that ‘Afun(x) = A * x’. Additional parameters to
‘Afun’ may be passed after X0.
A has to be Hermitian and Positive Definite (HPD). If ‘pcg’
detects A not to be positive definite, a warning is printed
and the FLAG output is set.
• B is the right-hand side vector.
• TOL is the required relative tolerance for the residual error,
‘B - A * X’. The iteration stops if
‘norm (B - A * X)’ ≤ ‘TOL * norm (B)’. If TOL is omitted or
empty, then a tolerance of 1e-6 is used.
• MAXIT is the maximum allowed number of iterations; if MAXIT is
omitted or empty then a value of 20 is used.
• M is a HPD preconditioning matrix. For any decomposition ‘M =
P1 * P2’ such that ‘inv (P1) * A * inv (P2)’ is HPD, the
conjugate gradient method is formally applied to the linear
system ‘inv (P1) * A * inv (P2) * Y = inv (P1) * B’, with ‘X =
inv (P2) * Y’ (split preconditioning). In practice, at each
iteration of the conjugate gradient method a linear system
with matrix M is solved with ‘mldivide’. If a particular
factorization ‘M = M1 * M2’ is available (for instance, an
incomplete Cholesky factorization of A), the two matrices M1
and M2 can be passed and the relative linear systems are
solved with the ‘mldivide’ operator. Note that a proper
choice of the preconditioner may dramatically improve the
overall performance of the method. Instead of matrices M1 and
M2, the user may pass two functions which return the results
of applying the inverse of M1 and M2 to a vector. If M1 is
omitted or empty ‘[]’, then no preconditioning is applied. If
no factorization of M is available, M2 can be omitted or left
[], and the input variable M1 can be used to pass the
preconditioner M.
• X0 is the initial guess. If X0 is omitted or empty then the
function sets X0 to a zero vector by default.
The arguments which follow X0 are treated as parameters, and passed
in an appropriate manner to any of the functions (A or M1 or M2)
that have been given to ‘pcg’. See the examples below for further
details.
The output arguments are:
• X is the computed approximation to the solution of
‘A * X = B’. If the algorithm did not converge, then X is the
iteration which has the minimum residual.
• FLAG reports on the convergence:
• 0: The algorithm converged to within the prescribed
tolerance.
• 1: The algorithm did not converge and it reached the
maximum number of iterations.
• 2: The preconditioner matrix is singular.
• 3: The algorithm stagnated, i.e., the absolute value of
the difference between the current iteration X and the
previous is less than ‘EPS * norm (X,2)’.
• 4: The algorithm detects that the input (preconditioned)
matrix is not HPD.
• RELRES is the ratio of the final residual to its initial
value, measured in the Euclidean norm.
• ITER indicates the iteration of X which it was computed.
Since the output X corresponds to the minimal residual
solution, the total number of iterations that the method
performed is given by ‘length(resvec) - 1’.
• RESVEC describes the convergence history of the method.
‘RESVEC (I, 1)’ is the Euclidean norm of the residual, and
‘RESVEC (I, 2)’ is the preconditioned residual norm, after the
(I-1)-th iteration, ‘I = 1, 2, ..., ITER+1’. The
preconditioned residual norm is defined as ‘R' * (M \ R)’
where ‘R = B - A * X’, see also the description of M. If
EIGEST is not required, only ‘RESVEC (:, 1)’ is returned.
• EIGEST returns the estimate for the smallest ‘EIGEST(1)’ and
largest ‘EIGEST(2)’ eigenvalues of the preconditioned matrix
‘P = M \ A’. In particular, if no preconditioning is used,
the estimates for the extreme eigenvalues of A are returned.
‘EIGEST(1)’ is an overestimate and ‘EIGEST(2)’ is an
underestimate, so that ‘EIGEST(2) / EIGEST(1)’ is a lower
bound for ‘cond (P, 2)’, which nevertheless in the limit
should theoretically be equal to the actual value of the
condition number.
Let us consider a trivial problem with a tridiagonal matrix
n = 10;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
b = A * ones (n, 1);
M1 = ichol (A); # in this tridiagonal case it corresponds to chol (A)'
M2 = M1';
M = M1 * M2;
Afun = @(x) A * x;
Mfun = @(x) M \ x;
M1fun = @(x) M1 \ x;
M2fun = @(x) M2 \ x;
EXAMPLE 1: Simplest use of ‘pcg’
x = pcg (A, b)
EXAMPLE 2: ‘pcg’ with a function which computes ‘A * X’
x = pcg (Afun, b)
EXAMPLE 3: ‘pcg’ with a preconditioner matrix M
x = pcg (A, b, 1e-06, 100, M)
EXAMPLE 4: ‘pcg’ with a function as preconditioner
x = pcg (Afun, b, 1e-6, 100, Mfun)
EXAMPLE 5: ‘pcg’ with preconditioner matrices M1 and M2
x = pcg (A, b, 1e-6, 100, M1, M2)
EXAMPLE 6: ‘pcg’ with functions as preconditioners
x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun)
EXAMPLE 7: ‘pcg’ with as input a function requiring an argument
function y = Ap (A, x, p) # compute A^p * x
y = x;
for i = 1:p
y = A * y;
endfor
endfunction
Apfun = @(x, p) Ap (A, x, p);
x = pcg (Apfun, b, [], [], [], [], [], 2);
EXAMPLE 8: explicit example to show that ‘pcg’ uses a split
preconditioner
M1 = ichol (A + 0.1 * eye (n)); # factorization of A perturbed
M2 = M1';
M = M1 * M2;
## reference solution computed by pcg after two iterations
[x_ref, fl] = pcg (A, b, [], 2, M)
## split preconditioning
[y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2)
x = M2 \ y # compare x and x_ref
References:
1. C.T. Kelley, ‘Iterative Methods for Linear and Nonlinear
Equations’, SIAM, 1995. (the base PCG algorithm)
2. Y. Saad, ‘Iterative Methods for Sparse Linear Systems’, PWS
1996. (condition number estimate from PCG) Revised version of
this book is available online at
<https://www-users.cs.umn.edu/~saad/books.html>
See also: Note: sparse, Note: pcr, Note:
gmres, Note: bicg, *note bicgstab:
XREFbicgstab, Note: cgs.
-- : X = pcr (A, B, TOL, MAXIT, M, X0, ...)
-- : [X, FLAG, RELRES, ITER, RESVEC] = pcr (...)
Solve the linear system of equations ‘A * X = B’ by means of the
Preconditioned Conjugate Residuals iterative method.
The input arguments are
• A can be either a square (preferably sparse) matrix or a
function handle, inline function or string containing the name
of a function which computes ‘A * X’. In principle A should
be symmetric and non-singular; if ‘pcr’ finds A to be
numerically singular, you will get a warning message and the
FLAG output parameter will be set.
• B is the right hand side vector.
• TOL is the required relative tolerance for the residual error,
‘B - A * X’. The iteration stops if ‘norm (B - A * X) <= TOL
* norm (B - A * X0)’. If TOL is empty or is omitted, the
function sets ‘TOL = 1e-6’ by default.
• MAXIT is the maximum allowable number of iterations; if ‘[]’
is supplied for MAXIT, or ‘pcr’ has less arguments, a default
value equal to 20 is used.
• M is the (left) preconditioning matrix, so that the iteration
is (theoretically) equivalent to solving by ‘pcr’ ‘P * X = M \
B’, with ‘P = M \ A’. Note that a proper choice of the
preconditioner may dramatically improve the overall
performance of the method. Instead of matrix M, the user may
pass a function which returns the results of applying the
inverse of M to a vector (usually this is the preferred way of
using the preconditioner). If ‘[]’ is supplied for M, or M is
omitted, no preconditioning is applied.
• X0 is the initial guess. If X0 is empty or omitted, the
function sets X0 to a zero vector by default.
The arguments which follow X0 are treated as parameters, and passed
in a proper way to any of the functions (A or M) which are passed
to ‘pcr’. See the examples below for further details.
The output arguments are
• X is the computed approximation to the solution of ‘A * X =
B’.
• FLAG reports on the convergence. ‘FLAG = 0’ means the
solution converged and the tolerance criterion given by TOL is
satisfied. ‘FLAG = 1’ means that the MAXIT limit for the
iteration count was reached. ‘FLAG = 3’ reports a ‘pcr’
breakdown, see [1] for details.
• RELRES is the ratio of the final residual to its initial
value, measured in the Euclidean norm.
• ITER is the actual number of iterations performed.
• RESVEC describes the convergence history of the method, so
that ‘RESVEC (i)’ contains the Euclidean norms of the residual
after the (I-1)-th iteration, ‘I = 1,2, ..., ITER+1’.
Let us consider a trivial problem with a diagonal matrix (we
exploit the sparsity of A)
n = 10;
A = sparse (diag (1:n));
b = rand (N, 1);
EXAMPLE 1: Simplest use of ‘pcr’
x = pcr (A, b)
EXAMPLE 2: ‘pcr’ with a function which computes ‘A * X’.
function y = apply_a (x)
y = [1:10]' .* x;
endfunction
x = pcr ("apply_a", b)
EXAMPLE 3: Preconditioned iteration, with full diagnostics. The
preconditioner (quite strange, because even the original matrix A
is trivial) is defined as a function
function y = apply_m (x)
k = floor (length (x) - 2);
y = x;
y(1:k) = x(1:k) ./ [1:k]';
endfunction
[x, flag, relres, iter, resvec] = ...
pcr (A, b, [], [], "apply_m")
semilogy ([1:iter+1], resvec);
EXAMPLE 4: Finally, a preconditioner which depends on a parameter
K.
function y = apply_m (x, varargin)
k = varargin{1};
y = x;
y(1:k) = x(1:k) ./ [1:k]';
endfunction
[x, flag, relres, iter, resvec] = ...
pcr (A, b, [], [], "apply_m"', [], 3)
References:
[1] W. Hackbusch, ‘Iterative Solution of Large Sparse Systems of
Equations’, section 9.5.4; Springer, 1994
See also: Note: sparse, Note: pcg.
The speed with which an iterative solver converges to a solution can
be accelerated with the use of a pre-conditioning matrix M. In this
case the linear equation ‘M^-1 * X = M^-1 * A \ B’ is solved instead.
Typical pre-conditioning matrices are partial factorizations of the
original matrix.
-- : 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: Note: chol, Note: ilu, *note pcg:
XREFpcg.
-- : ilu (A)
-- : ilu (A, OPTS)
-- : [L, U] = ilu (...)
-- : [L, U, P] = ilu (...)
Compute the incomplete LU factorization of the sparse square matrix
A.
‘ilu’ returns a unit lower triangular matrix L, an upper triangular
matrix U, and optionally a permutation matrix P, such that ‘L*U’
approximates ‘P*A’.
The factors given by this routine may be useful as preconditioners
for a system of linear equations being solved by iterative methods
such as BICG (BiConjugate Gradients) or GMRES (Generalized Minimum
Residual Method).
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)
ILU factorization with no fill-in (ILU(0)).
Additional supported options: ‘milu’.
"crout"
Crout version of ILU factorization (ILUC).
Additional supported options: ‘milu’, ‘droptol’.
"ilutp"
ILU factorization with threshold and pivoting.
Additional supported options: ‘milu’, ‘droptol’, ‘udiag’,
‘thresh’.
‘droptol’
A non-negative scalar specifying the drop tolerance for
factorization. The default value is 0 which produces the
complete LU factorization.
Non-diagonal entries of U are set to 0 unless
‘abs (U(i,j)) >= droptol * norm (A(:,j))’.
Non-diagonal entries of L are set to 0 unless
‘abs (L(i,j)) >= droptol * norm (A(:,j))/U(j,j)’.
‘milu’
Modified incomplete LU factorization:
"row"
Row-sum modified incomplete LU factorization. The
factorization preserves row sums: ‘A * e = L * U * e’,
where e is a vector of ones.
"col"
Column-sum modified incomplete LU factorization. The
factorization preserves column sums: ‘e' * A = e' * L *
U’.
"off" (default)
Row and column sums are not necessarily preserved.
‘udiag’
If true, any zeros on the diagonal of the upper triangular
factor are replaced by the local drop tolerance ‘droptol *
norm (A(:,j))/U(j,j)’. The default is false.
‘thresh’
Pivot threshold for factorization. It can range between 0
(diagonal pivoting) and 1 (default), where the maximum
magnitude entry in the column is chosen to be the pivot.
If ‘ilu’ is called with just one output, the returned matrix is ‘L
+ U - speye (size (A))’, where L is unit lower triangular and U is
upper triangular.
With two outputs, ‘ilu’ returns a unit lower triangular matrix L
and an upper triangular matrix U. For OPTS.type == "ilutp", one of
the factors is permuted based on the value of OPTS.milu. When
OPTS.milu == "row", U is a column permuted upper triangular factor.
Otherwise, L is a row-permuted unit lower triangular factor.
If there are three named outputs and OPTS.milu != "row", P is
returned such that L and U are incomplete factors of ‘P*A’. When
OPTS.milu == "row", P is returned such that L and U are incomplete
factors of ‘A*P’.
EXAMPLES
A = gallery ("neumann", 1600) + speye (1600);
opts.type = "nofill";
nnz (A)
ans = 7840
nnz (lu (A))
ans = 126478
nnz (ilu (A, opts))
ans = 7840
This shows that A has 7,840 nonzeros, the complete LU factorization
has 126,478 nonzeros, and the incomplete LU factorization, with 0
level of fill-in, has 7,840 nonzeros, the same amount as A. Taken
from: <http://www.mathworks.com/help/matlab/ref/ilu.html>
A = gallery ("wathen", 10, 10);
b = sum (A, 2);
tol = 1e-8;
maxit = 50;
opts.type = "crout";
opts.droptol = 1e-4;
[L, U] = ilu (A, opts);
x = bicg (A, b, tol, maxit, L, U);
norm (A * x - b, inf)
This example uses ILU as preconditioner for a random FEM-Matrix,
which has a large condition number. Without L and U BICG would not
converge.
See also: Note: lu, Note: ichol, *note bicg:
XREFbicg, Note: gmres.
automatically generated by info2www version 1.2.2.9