Next: , Previous: , Up: Sparse Matrices   [Contents][Index]


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:

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:

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: sparse, pcr, gmres, bicg, bicgstab, 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

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

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)

Reference:

W. Hackbusch, Iterative Solution of Large Sparse Systems of Equations, section 9.5.4; Springer, 1994

See also: sparse, 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: chol, ilu, pcg.

: 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: 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.

See also: lu, ichol, bicg, gmres.


Next: , Previous: , Up: Sparse Matrices   [Contents][Index]