1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
|
<html lang="en">
<head>
<title>Iterative Techniques - Untitled</title>
<meta http-equiv="Content-Type" content="text/html">
<meta name="description" content="Untitled">
<meta name="generator" content="makeinfo 4.11">
<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">
<p>
<a name="Iterative-Techniques"></a>
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">21.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.
<!-- ./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-1752"></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-1753"></a></var><br>
<blockquote>
<p>Solves 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><span class="sc">References</span>
<p> [1] C.T.Kelley, 'Iterative methods for linear and nonlinear equations',
SIAM, 1995 (the base PCG algorithm)
<p> [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 http://www-users.cs.umn.edu/~saad/books.html
<!-- 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>
<!-- ./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-1754"></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-1755"></a></var><br>
<blockquote>
<p>Solves 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><span class="sc">References</span>
<p> [1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of
Equations", 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.
<!-- ./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-1756"></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-1757"></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-1758"></a></var><br>
<blockquote><p><a name="index-LU-decomposition-1759"></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 <dfn>luinc</dfn>.
<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 UMFPACK 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 U 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 U
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 <dfn>luinc</dfn>
are the same as for <dfn>lu</dfn>.
<p>Given the string argument 'vector', <dfn>luinc</dfn> 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>
|
