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