<html lang="en">
<head>
<title>Iterative Techniques - GNU Octave</title>
<meta http-equiv="Content-Type" content="text/html">
<meta name="description" content="GNU Octave">
<meta name="generator" content="makeinfo 4.13">
<link title="Top" rel="start" href="index.html#Top">
<link rel="up" href="Sparse-Matrices.html#Sparse-Matrices" title="Sparse Matrices">
<link rel="prev" href="Sparse-Linear-Algebra.html#Sparse-Linear-Algebra" title="Sparse Linear Algebra">
<link rel="next" href="Real-Life-Example.html#Real-Life-Example" title="Real Life Example">
<link href="http://www.gnu.org/software/texinfo/" rel="generator-home" title="Texinfo Homepage">
<meta http-equiv="Content-Style-Type" content="text/css">
<style type="text/css"><!--
pre.display { font-family:inherit }
pre.format { font-family:inherit }
pre.smalldisplay { font-family:inherit; font-size:smaller }
pre.smallformat { font-family:inherit; font-size:smaller }
pre.smallexample { font-size:smaller }
pre.smalllisp { font-size:smaller }
span.sc { font-variant:small-caps }
span.roman { font-family:serif; font-weight:normal; }
span.sansserif { font-family:sans-serif; font-weight:normal; }
--></style>
</head>
<body>
<div class="node">
<a name="Iterative-Techniques"></a>
<p>
Next: <a rel="next" accesskey="n" href="Real-Life-Example.html#Real-Life-Example">Real Life Example</a>,
Previous: <a rel="previous" accesskey="p" href="Sparse-Linear-Algebra.html#Sparse-Linear-Algebra">Sparse Linear Algebra</a>,
Up: <a rel="up" accesskey="u" href="Sparse-Matrices.html#Sparse-Matrices">Sparse Matrices</a>
<hr>
</div>
<h3 class="section">22.3 Iterative Techniques applied to sparse matrices</h3>
<p>The left division <code>\</code> and right division <code>/</code> operators,
discussed in the previous section, use direct solvers to resolve a
linear equation of the form <var>x</var><code> = </code><var>A</var><code> \ </code><var>b</var> or
<var>x</var><code> = </code><var>b</var><code> / </code><var>A</var>. Octave equally includes a number of
functions to solve sparse linear equations using iterative techniques.
<!-- pcg scripts/sparse/pcg.m -->
<p><a name="doc_002dpcg"></a>
<div class="defun">
— Function File: <var>x</var> = <b>pcg</b> (<var>A, b, tol, maxit, m1, m2, x0, <small class="dots">...</small></var>)<var><a name="index-pcg-2295"></a></var><br>
— Function File: [<var>x</var>, <var>flag</var>, <var>relres</var>, <var>iter</var>, <var>resvec</var>, <var>eigest</var>] = <b>pcg</b> (<var><small class="dots">...</small></var>)<var><a name="index-pcg-2296"></a></var><br>
<blockquote>
<p>Solve the linear system of equations <var>A</var><code> * </code><var>x</var><code> = </code><var>b</var>
by means of the Preconditioned Conjugate Gradient iterative
method. The input arguments are
<ul>
<li><var>A</var> can be either a square (preferably sparse) matrix or a
function handle, inline function or string containing the name
of a function which computes <var>A</var><code> * </code><var>x</var>. In principle
<var>A</var> should be symmetric and positive definite; if <code>pcg</code>
finds <var>A</var> to not be positive definite, you will get a warning
message and the <var>flag</var> output parameter will be set.
<li><var>b</var> is the right hand side vector.
<li><var>tol</var> is the required relative tolerance for the residual error,
<var>b</var><code> - </code><var>A</var><code> * </code><var>x</var>. The iteration stops if
<code>norm (</code><var>b</var><code> - </code><var>A</var><code> * </code><var>x</var><code>) <=
</code><var>tol</var><code> * norm (</code><var>b</var><code> - </code><var>A</var><code> * </code><var>x0</var><code>)</code>.
If <var>tol</var> is empty or is omitted, the function sets
<var>tol</var><code> = 1e-6</code> by default.
<li><var>maxit</var> is the maximum allowable number of iterations; if
<code>[]</code> is supplied for <code>maxit</code>, or <code>pcg</code> has less
arguments, a default value equal to 20 is used.
<li><var>m</var> = <var>m1</var> * <var>m2</var> is the (left) preconditioning matrix, so that
the iteration is (theoretically) equivalent to solving by <code>pcg</code>
<var>P</var><code> *
</code><var>x</var><code> = </code><var>m</var><code> \ </code><var>b</var>, with <var>P</var><code> = </code><var>m</var><code> \ </code><var>A</var>.
Note that a proper choice of the preconditioner may dramatically
improve the overall performance of the method. Instead of matrices
<var>m1</var> and <var>m2</var>, the user may pass two functions which return
the results of applying the inverse of <var>m1</var> and <var>m2</var> to
a vector (usually this is the preferred way of using the preconditioner).
If <code>[]</code> is supplied for <var>m1</var>, or <var>m1</var> is omitted, no
preconditioning is applied. If <var>m2</var> is omitted, <var>m</var> = <var>m1</var>
will be used as preconditioner.
<li><var>x0</var> is the initial guess. If <var>x0</var> is empty or omitted, the
function sets <var>x0</var> to a zero vector by default.
</ul>
<p>The arguments which follow <var>x0</var> are treated as parameters, and
passed in a proper way to any of the functions (<var>A</var> or <var>m</var>)
which are passed to <code>pcg</code>. See the examples below for further
details. The output arguments are
<ul>
<li><var>x</var> is the computed approximation to the solution of
<var>A</var><code> * </code><var>x</var><code> = </code><var>b</var>.
<li><var>flag</var> reports on the convergence. <var>flag</var><code> = 0</code> means
the solution converged and the tolerance criterion given by <var>tol</var>
is satisfied. <var>flag</var><code> = 1</code> means that the <var>maxit</var> limit
for the iteration count was reached. <var>flag</var><code> = 3</code> reports that
the (preconditioned) matrix was found not positive definite.
<li><var>relres</var> is the ratio of the final residual to its initial value,
measured in the Euclidean norm.
<li><var>iter</var> is the actual number of iterations performed.
<li><var>resvec</var> describes the convergence history of the method.
<var>resvec</var><code> (i,1)</code> is the Euclidean norm of the residual, and
<var>resvec</var><code> (i,2)</code> is the preconditioned residual norm,
after the (<var>i</var>-1)-th iteration, <var>i</var><code> =
1, 2, ..., </code><var>iter</var><code>+1</code>. The preconditioned residual norm
is defined as
<code>norm (</code><var>r</var><code>) ^ 2 = </code><var>r</var><code>' * (</code><var>m</var><code> \ </code><var>r</var><code>)</code> where
<var>r</var><code> = </code><var>b</var><code> - </code><var>A</var><code> * </code><var>x</var>, see also the
description of <var>m</var>. If <var>eigest</var> is not required, only
<var>resvec</var><code> (:,1)</code> is returned.
<li><var>eigest</var> returns the estimate for the smallest <var>eigest</var><code>
(1)</code> and largest <var>eigest</var><code> (2)</code> eigenvalues of the
preconditioned matrix <var>P</var><code> = </code><var>m</var><code> \ </code><var>A</var>. In
particular, if no preconditioning is used, the estimates for the
extreme eigenvalues of <var>A</var> are returned. <var>eigest</var><code> (1)</code>
is an overestimate and <var>eigest</var><code> (2)</code> is an underestimate,
so that <var>eigest</var><code> (2) / </code><var>eigest</var><code> (1)</code> is a lower bound
for <code>cond (</code><var>P</var><code>, 2)</code>, which nevertheless in the limit should
theoretically be equal to the actual value of the condition number.
The method which computes <var>eigest</var> works only for symmetric positive
definite <var>A</var> and <var>m</var>, and the user is responsible for
verifying this assumption.
</ul>
<p>Let us consider a trivial problem with a diagonal matrix (we exploit the
sparsity of A)
<pre class="example"> n = 10;
A = diag (sparse (1:n));
b = rand (n, 1);
[l, u, p, q] = luinc (A, 1.e-3);
</pre>
<p><span class="sc">Example 1:</span> Simplest use of <code>pcg</code>
<pre class="example"> x = pcg (A,b)
</pre>
<p><span class="sc">Example 2:</span> <code>pcg</code> with a function which computes
<var>A</var><code> * </code><var>x</var>
<pre class="example"> function y = apply_a (x)
y = [1:N]' .* x;
endfunction
x = pcg ("apply_a", b)
</pre>
<p><span class="sc">Example 3:</span> <code>pcg</code> with a preconditioner: <var>l</var> * <var>u</var>
<pre class="example"> x = pcg (A, b, 1.e-6, 500, l*u)
</pre>
<p><span class="sc">Example 4:</span> <code>pcg</code> with a preconditioner: <var>l</var> * <var>u</var>.
Faster than <span class="sc">Example 3</span> since lower and upper triangular matrices
are easier to invert
<pre class="example"> x = pcg (A, b, 1.e-6, 500, l, u)
</pre>
<p><span class="sc">Example 5:</span> Preconditioned iteration, with full diagnostics. The
preconditioner (quite strange, because even the original matrix
<var>A</var> is trivial) is defined as a function
<pre class="example"> 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, eigest] = ...
pcg (A, b, [], [], "apply_m");
semilogy (1:iter+1, resvec);
</pre>
<p><span class="sc">Example 6:</span> Finally, a preconditioner which depends on a
parameter <var>k</var>.
<pre class="example"> 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, eigest] = ...
pcg (A, b, [], [], "apply_m", [], [], 3)
</pre>
<p>References:
<ol type=1 start=1>
<li>C.T. Kelley, <cite>Iterative Methods for Linear and Nonlinear Equations</cite>,
SIAM, 1995. (the base PCG algorithm)
<li>Y. Saad, <cite>Iterative Methods for Sparse Linear Systems</cite>, PWS 1996.
(condition number estimate from PCG) Revised version of this book is
available online at <a href="http://www-users.cs.umn.edu/~saad/books.html">http://www-users.cs.umn.edu/~saad/books.html</a>
</ol>
<!-- Texinfo @sp should work but in practice produces ugly results for HTML. -->
<!-- A simple blank line produces the correct behavior. -->
<!-- @sp 1 -->
<p class="noindent"><strong>See also:</strong> <a href="doc_002dsparse.html#doc_002dsparse">sparse</a>, <a href="doc_002dpcr.html#doc_002dpcr">pcr</a>.
</p></blockquote></div>
<!-- pcr scripts/sparse/pcr.m -->
<p><a name="doc_002dpcr"></a>
<div class="defun">
— Function File: <var>x</var> = <b>pcr</b> (<var>A, b, tol, maxit, m, x0, <small class="dots">...</small></var>)<var><a name="index-pcr-2297"></a></var><br>
— Function File: [<var>x</var>, <var>flag</var>, <var>relres</var>, <var>iter</var>, <var>resvec</var>] = <b>pcr</b> (<var><small class="dots">...</small></var>)<var><a name="index-pcr-2298"></a></var><br>
<blockquote>
<p>Solve the linear system of equations <var>A</var><code> * </code><var>x</var><code> = </code><var>b</var>
by means of the Preconditioned Conjugate Residuals iterative
method. The input arguments are
<ul>
<li><var>A</var> can be either a square (preferably sparse) matrix or a
function handle, inline function or string containing the name
of a function which computes <var>A</var><code> * </code><var>x</var>. In principle
<var>A</var> should be symmetric and non-singular; if <code>pcr</code>
finds <var>A</var> to be numerically singular, you will get a warning
message and the <var>flag</var> output parameter will be set.
<li><var>b</var> is the right hand side vector.
<li><var>tol</var> is the required relative tolerance for the residual error,
<var>b</var><code> - </code><var>A</var><code> * </code><var>x</var>. The iteration stops if
<code>norm (</code><var>b</var><code> - </code><var>A</var><code> * </code><var>x</var><code>) <=
</code><var>tol</var><code> * norm (</code><var>b</var><code> - </code><var>A</var><code> * </code><var>x0</var><code>)</code>.
If <var>tol</var> is empty or is omitted, the function sets
<var>tol</var><code> = 1e-6</code> by default.
<li><var>maxit</var> is the maximum allowable number of iterations; if
<code>[]</code> is supplied for <code>maxit</code>, or <code>pcr</code> has less
arguments, a default value equal to 20 is used.
<li><var>m</var> is the (left) preconditioning matrix, so that the iteration is
(theoretically) equivalent to solving by <code>pcr</code> <var>P</var><code> *
</code><var>x</var><code> = </code><var>m</var><code> \ </code><var>b</var>, with <var>P</var><code> = </code><var>m</var><code> \ </code><var>A</var>.
Note that a proper choice of the preconditioner may dramatically
improve the overall performance of the method. Instead of matrix
<var>m</var>, the user may pass a function which returns the results of
applying the inverse of <var>m</var> to a vector (usually this is the
preferred way of using the preconditioner). If <code>[]</code> is supplied
for <var>m</var>, or <var>m</var> is omitted, no preconditioning is applied.
<li><var>x0</var> is the initial guess. If <var>x0</var> is empty or omitted, the
function sets <var>x0</var> to a zero vector by default.
</ul>
<p>The arguments which follow <var>x0</var> are treated as parameters, and
passed in a proper way to any of the functions (<var>A</var> or <var>m</var>)
which are passed to <code>pcr</code>. See the examples below for further
details. The output arguments are
<ul>
<li><var>x</var> is the computed approximation to the solution of
<var>A</var><code> * </code><var>x</var><code> = </code><var>b</var>.
<li><var>flag</var> reports on the convergence. <var>flag</var><code> = 0</code> means
the solution converged and the tolerance criterion given by <var>tol</var>
is satisfied. <var>flag</var><code> = 1</code> means that the <var>maxit</var> limit
for the iteration count was reached. <var>flag</var><code> = 3</code> reports t
<code>pcr</code> breakdown, see [1] for details.
<li><var>relres</var> is the ratio of the final residual to its initial value,
measured in the Euclidean norm.
<li><var>iter</var> is the actual number of iterations performed.
<li><var>resvec</var> describes the convergence history of the method,
so that <var>resvec</var><code> (i)</code> contains the Euclidean norms of the
residual after the (<var>i</var>-1)-th iteration, <var>i</var><code> =
1,2, ..., </code><var>iter</var><code>+1</code>.
</ul>
<p>Let us consider a trivial problem with a diagonal matrix (we exploit the
sparsity of A)
<pre class="example"> n = 10;
A = sparse (diag (1:n));
b = rand (N, 1);
</pre>
<p><span class="sc">Example 1:</span> Simplest use of <code>pcr</code>
<pre class="example"> x = pcr (A, b)
</pre>
<p><span class="sc">Example 2:</span> <code>pcr</code> with a function which computes
<var>A</var><code> * </code><var>x</var>.
<pre class="example"> function y = apply_a (x)
y = [1:10]' .* x;
endfunction
x = pcr ("apply_a", b)
</pre>
<p><span class="sc">Example 3:</span> Preconditioned iteration, with full diagnostics. The
preconditioner (quite strange, because even the original matrix
<var>A</var> is trivial) is defined as a function
<pre class="example"> 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);
</pre>
<p><span class="sc">Example 4:</span> Finally, a preconditioner which depends on a
parameter <var>k</var>.
<pre class="example"> 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)
</pre>
<p>References:
<p>[1] W. Hackbusch, <cite>Iterative Solution of Large Sparse Systems of
Equations</cite>, section 9.5.4; Springer, 1994
<!-- Texinfo @sp should work but in practice produces ugly results for HTML. -->
<!-- A simple blank line produces the correct behavior. -->
<!-- @sp 1 -->
<p class="noindent"><strong>See also:</strong> <a href="doc_002dsparse.html#doc_002dsparse">sparse</a>, <a href="doc_002dpcg.html#doc_002dpcg">pcg</a>.
</p></blockquote></div>
<p>The speed with which an iterative solver converges to a solution can be
accelerated with the use of a pre-conditioning matrix <var>M</var>. In this
case the linear equation <var>M</var><code>^-1 * </code><var>x</var><code> = </code><var>M</var><code>^-1 *
</code><var>A</var><code> \ </code><var>b</var> is solved instead. Typical pre-conditioning matrices
are partial factorizations of the original matrix.
<!-- luinc src/DLD-FUNCTIONS/luinc.cc -->
<p><a name="doc_002dluinc"></a>
<div class="defun">
— Loadable Function: [<var>L</var>, <var>U</var>, <var>P</var>, <var>Q</var>] = <b>luinc</b> (<var>A, '0'</var>)<var><a name="index-luinc-2299"></a></var><br>
— Loadable Function: [<var>L</var>, <var>U</var>, <var>P</var>, <var>Q</var>] = <b>luinc</b> (<var>A, droptol</var>)<var><a name="index-luinc-2300"></a></var><br>
— Loadable Function: [<var>L</var>, <var>U</var>, <var>P</var>, <var>Q</var>] = <b>luinc</b> (<var>A, opts</var>)<var><a name="index-luinc-2301"></a></var><br>
<blockquote><p><a name="index-LU-decomposition-2302"></a>Produce the incomplete LU factorization of the sparse matrix <var>A</var>.
Two types of incomplete factorization are possible, and the type
is determined by the second argument to <code>luinc</code>.
<p>Called with a second argument of '0', the zero-level incomplete
LU factorization is produced. This creates a factorization of <var>A</var>
where the position of the non-zero arguments correspond to the same
positions as in the matrix <var>A</var>.
<p>Alternatively, the fill-in of the incomplete LU factorization can
be controlled through the variable <var>droptol</var> or the structure
<var>opts</var>. The <span class="sc">umfpack</span> multifrontal factorization code by Tim A.
Davis is used for the incomplete LU factorization, (availability
<a href="http://www.cise.ufl.edu/research/sparse/umfpack/">http://www.cise.ufl.edu/research/sparse/umfpack/</a>)
<p><var>droptol</var> determines the values below which the values in the
LU factorization are dropped and replaced by zero. It must be a
positive scalar, and any values in the factorization whose absolute value
are less than this value are dropped, expect if leaving them increase the
sparsity of the matrix. Setting <var>droptol</var> to zero results in a complete
LU factorization which is the default.
<p><var>opts</var> is a structure containing one or more of the fields
<dl>
<dt><code>droptol</code><dd>The drop tolerance as above. If <var>opts</var> only contains <code>droptol</code>
then this is equivalent to using the variable <var>droptol</var>.
<br><dt><code>milu</code><dd>A logical variable flagging whether to use the modified incomplete
LU factorization. In the case that <code>milu</code> is true, the dropped
values are subtracted from the diagonal of the matrix <var>U</var> of the
factorization. The default is <code>false</code>.
<br><dt><code>udiag</code><dd>A logical variable that flags whether zero elements on the diagonal of
<var>U</var> should be replaced with <var>droptol</var> to attempt to avoid singular
factors. The default is <code>false</code>.
<br><dt><code>thresh</code><dd>Defines the pivot threshold in the interval [0,1]. Values outside that
range are ignored.
</dl>
<p>All other fields in <var>opts</var> are ignored. The outputs from <code>luinc</code>
are the same as for <code>lu</code>.
<p>Given the string argument 'vector', <code>luinc</code> returns the values of
<var>p</var> <var>q</var> as vector values.
<!-- Texinfo @sp should work but in practice produces ugly results for HTML. -->
<!-- A simple blank line produces the correct behavior. -->
<!-- @sp 1 -->
<p class="noindent"><strong>See also:</strong> <a href="doc_002dsparse.html#doc_002dsparse">sparse</a>, <a href="doc_002dlu.html#doc_002dlu">lu</a>.
</p></blockquote></div>
</body></html>
