Next: Real Life Example, Previous: Sparse Linear Algebra, Up: Sparse Matrices [Contents][Index]
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.
Solve the linear system of equations A * x = b
by means of the Preconditioned Conjugate Gradient iterative method.
The input arguments are:
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 - 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.
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.
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:
A * x = b
. If the algorithm did not converge,
then x is the iteration which has the minimum residual.
eps * norm (x,2)
.
length(resvec) - 1
.
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(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:
Solve the linear system of equations A * x = b
by
means of the Preconditioned Conjugate Residuals iterative method.
The input arguments are
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 - 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.
[]
is
supplied for maxit, or pcr
has less arguments, a default
value equal to 20 is used.
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.
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
A * x = b
.
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.
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
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.
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 of factorization.
"nofill"
(default)Incomplete Cholesky factorization with no fill-in (IC(0)).
"ict"
Incomplete Cholesky factorization with threshold dropping (ICT).
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.
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)
.
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.
"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.
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: https://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.
Next: Real Life Example, Previous: Sparse Linear Algebra, Up: Sparse Matrices [Contents][Index]