|
3294
|
1 @c Copyright (C) 1996, 1997 John W. Eaton |
|
|
2 @c This is part of the Octave manual. |
|
|
3 @c For copying conditions, see the file gpl.texi. |
|
|
4 |
|
|
5 @node Linear Algebra, Nonlinear Equations, Arithmetic, Top |
|
|
6 @chapter Linear Algebra |
|
|
7 |
|
|
8 This chapter documents the linear algebra functions of Octave. |
|
|
9 Reference material for many of these functions may be found in |
|
|
10 Golub and Van Loan, @cite{Matrix Computations, 2nd Ed.}, Johns Hopkins, |
|
|
11 1989, and in @cite{@sc{Lapack} Users' Guide}, SIAM, 1992. |
|
|
12 |
|
|
13 @menu |
|
|
14 * Basic Matrix Functions:: |
|
|
15 * Matrix Factorizations:: |
|
|
16 * Functions of a Matrix:: |
|
|
17 @end menu |
|
|
18 |
|
|
19 @node Basic Matrix Functions, Matrix Factorizations, Linear Algebra, Linear Algebra |
|
|
20 @section Basic Matrix Functions |
|
|
21 |
|
|
22 @deftypefn {Loadable Function} {@var{aa} =} balance (@var{a}, @var{opt}) |
|
|
23 @deftypefnx {Loadable Function} {[@var{dd}, @var{aa}] =} balance (@var{a}, @var{opt}) |
|
|
24 @deftypefnx {Loadable Function} {[@var{cc}, @var{dd}, @var{aa}, @var{bb]} =} balance (@var{a}, @var{b}, @var{opt}) |
|
|
25 |
|
|
26 @code{[dd, aa] = balance (a)} returns @code{aa = dd \ a * dd}. |
|
|
27 @code{aa} is a matrix whose row and column norms are roughly equal in |
|
|
28 magnitude, and @code{dd} = @code{p * d}, where @code{p} is a permutation |
|
|
29 matrix and @code{d} is a diagonal matrix of powers of two. This allows |
|
|
30 the equilibration to be computed without roundoff. Results of |
|
|
31 eigenvalue calculation are typically improved by balancing first. |
|
|
32 |
|
|
33 @code{[cc, dd, aa, bb] = balance (a, b)} returns @code{aa = cc*a*dd} and |
|
|
34 @code{bb = cc*b*dd)}, where @code{aa} and @code{bb} have non-zero |
|
|
35 elements of approximately the same magnitude and @code{cc} and @code{dd} |
|
|
36 are permuted diagonal matrices as in @code{dd} for the algebraic |
|
|
37 eigenvalue problem. |
|
|
38 |
|
|
39 The eigenvalue balancing option @code{opt} is selected as follows: |
|
|
40 |
|
|
41 @table @asis |
|
|
42 @item @code{"N"}, @code{"n"} |
|
|
43 No balancing; arguments copied, transformation(s) set to identity. |
|
|
44 |
|
|
45 @item @code{"P"}, @code{"p"} |
|
|
46 Permute argument(s) to isolate eigenvalues where possible. |
|
|
47 |
|
|
48 @item @code{"S"}, @code{"s"} |
|
|
49 Scale to improve accuracy of computed eigenvalues. |
|
|
50 |
|
|
51 @item @code{"B"}, @code{"b"} |
|
|
52 Permute and scale, in that order. Rows/columns of a (and b) |
|
|
53 that are isolated by permutation are not scaled. This is the default |
|
|
54 behavior. |
|
|
55 @end table |
|
|
56 |
|
|
57 Algebraic eigenvalue balancing uses standard @sc{Lapack} routines. |
|
|
58 |
|
|
59 Generalized eigenvalue problem balancing uses Ward's algorithm |
|
|
60 (SIAM Journal on Scientific and Statistical Computing, 1981). |
|
|
61 @end deftypefn |
|
|
62 |
|
|
63 @deftypefn {} {} cond (@var{a}) |
|
|
64 Compute the (two-norm) condition number of a matrix. @code{cond (a)} is |
|
|
65 defined as @code{norm (a) * norm (inv (a))}, and is computed via a |
|
|
66 singular value decomposition. |
|
|
67 @end deftypefn |
|
|
68 |
|
|
69 @deftypefn {Loadable Function} {} det (@var{a}) |
|
|
70 Compute the determinant of @var{a} using @sc{Linpack}. |
|
|
71 @end deftypefn |
|
|
72 |
|
|
73 @deftypefn {Loadable Function} {@var{lambda} =} eig (@var{a}) |
|
|
74 @deftypefnx {Loadable Function} {[@var{v}, @var{lambda}] =} eig (@var{a}) |
|
|
75 The eigenvalues (and eigenvectors) of a matrix are computed in a several |
|
|
76 step process which begins with a Hessenberg decomposition, followed by a |
|
|
77 Schur decomposition, from which the eigenvalues are apparent. The |
|
|
78 eigenvectors, when desired, are computed by further manipulations of the |
|
|
79 Schur decomposition. |
|
|
80 @end deftypefn |
|
|
81 |
|
|
82 @deftypefn {Loadable Function} {@var{G} =} givens (@var{x}, @var{y}) |
|
|
83 @deftypefnx {Loadable Function} {[@var{c}, @var{s}] =} givens (@var{x}, @var{y}) |
|
|
84 @iftex |
|
|
85 @tex |
|
|
86 Return a $2\times 2$ orthogonal matrix |
|
|
87 $$ |
|
|
88 G = \left[\matrix{c & s\cr -s'& c\cr}\right] |
|
|
89 $$ |
|
|
90 such that |
|
|
91 $$ |
|
|
92 G \left[\matrix{x\cr y}\right] = \left[\matrix{\ast\cr 0}\right] |
|
|
93 $$ |
|
|
94 with $x$ and $y$ scalars. |
|
|
95 @end tex |
|
|
96 @end iftex |
|
|
97 @ifinfo |
|
|
98 Return a 2 by 2 orthogonal matrix |
|
|
99 @code{@var{G} = [@var{c} @var{s}; -@var{s}' @var{c}]} such that |
|
|
100 @code{@var{G} [@var{x}; @var{y}] = [*; 0]} with @var{x} and @var{y} scalars. |
|
|
101 @end ifinfo |
|
|
102 |
|
|
103 For example, |
|
|
104 |
|
|
105 @example |
|
|
106 @group |
|
|
107 givens (1, 1) |
|
|
108 @result{} 0.70711 0.70711 |
|
|
109 -0.70711 0.70711 |
|
|
110 @end group |
|
|
111 @end example |
|
|
112 @end deftypefn |
|
|
113 |
|
|
114 @deftypefn {Loadable Function} {} inv (@var{a}) |
|
|
115 @deftypefnx {Loadable Function} {} inverse (@var{a}) |
|
|
116 Compute the inverse of the square matrix @var{a}. |
|
|
117 @end deftypefn |
|
|
118 |
|
|
119 @deftypefn {Function File} {} norm (@var{a}, @var{p}) |
|
|
120 Compute the p-norm of the matrix @var{a}. If the second argument is |
|
|
121 missing, @code{p = 2} is assumed. |
|
|
122 |
|
|
123 If @var{a} is a matrix: |
|
|
124 |
|
|
125 @table @asis |
|
|
126 @item @var{p} = @code{1} |
|
|
127 1-norm, the largest column sum of @var{a}. |
|
|
128 |
|
|
129 @item @var{p} = @code{2} |
|
|
130 Largest singular value of @var{a}. |
|
|
131 |
|
|
132 @item @var{p} = @code{Inf} |
|
|
133 @cindex infinity norm |
|
|
134 Infinity norm, the largest row sum of @var{a}. |
|
|
135 |
|
|
136 @item @var{p} = @code{"fro"} |
|
|
137 @cindex Frobenius norm |
|
|
138 Frobenius norm of @var{a}, @code{sqrt (sum (diag (@var{a}' * @var{a})))}. |
|
|
139 @end table |
|
|
140 |
|
|
141 If @var{a} is a vector or a scalar: |
|
|
142 |
|
|
143 @table @asis |
|
|
144 @item @var{p} = @code{Inf} |
|
|
145 @code{max (abs (@var{a}))}. |
|
|
146 |
|
|
147 @item @var{p} = @code{-Inf} |
|
|
148 @code{min (abs (@var{a}))}. |
|
|
149 |
|
|
150 @item other |
|
|
151 p-norm of @var{a}, @code{(sum (abs (@var{a}) .^ @var{p})) ^ (1/@var{p})}. |
|
|
152 @end table |
|
|
153 @end deftypefn |
|
|
154 |
|
|
155 @deftypefn {Function File} {} null (@var{a}, @var{tol}) |
|
|
156 Return an orthonormal basis of the null space of @var{a}. |
|
|
157 |
|
|
158 The dimension of the null space is taken as the number of singular |
|
|
159 values of @var{a} not greater than @var{tol}. If the argument @var{tol} |
|
|
160 is missing, it is computed as |
|
|
161 |
|
|
162 @example |
|
|
163 max (size (@var{a})) * max (svd (@var{a})) * eps |
|
|
164 @end example |
|
|
165 @end deftypefn |
|
|
166 |
|
|
167 @deftypefn {Function File} {} orth (@var{a}, @var{tol}) |
|
|
168 Return an orthonormal basis of the range space of @var{a}. |
|
|
169 |
|
|
170 The dimension of the range space is taken as the number of singular |
|
|
171 values of @var{a} greater than @var{tol}. If the argument @var{tol} is |
|
|
172 missing, it is computed as |
|
|
173 |
|
|
174 @example |
|
|
175 max (size (@var{a})) * max (svd (@var{a})) * eps |
|
|
176 @end example |
|
|
177 @end deftypefn |
|
|
178 |
|
|
179 @deftypefn {Function File} {} pinv (@var{x}, @var{tol}) |
|
|
180 Return the pseudoinverse of @var{x}. Singular values less than |
|
|
181 @var{tol} are ignored. |
|
|
182 |
|
|
183 If the second argument is omitted, it is assumed that |
|
|
184 |
|
|
185 @example |
|
|
186 tol = max (size (@var{x})) * sigma_max (@var{x}) * eps, |
|
|
187 @end example |
|
|
188 |
|
|
189 @noindent |
|
|
190 where @code{sigma_max (@var{x})} is the maximal singular value of @var{x}. |
|
|
191 @end deftypefn |
|
|
192 |
|
|
193 @deftypefn {Function File} {} rank (@var{a}, @var{tol}) |
|
|
194 Compute the rank of @var{a}, using the singular value decomposition. |
|
|
195 The rank is taken to be the number of singular values of @var{a} that |
|
|
196 are greater than the specified tolerance @var{tol}. If the second |
|
|
197 argument is omitted, it is taken to be |
|
|
198 |
|
|
199 @example |
|
|
200 tol = max (size (@var{a})) * sigma (1) * eps; |
|
|
201 @end example |
|
|
202 |
|
|
203 @noindent |
|
|
204 where @code{eps} is machine precision and @code{sigma} is the largest |
|
|
205 singular value of @var{a}. |
|
|
206 @end deftypefn |
|
|
207 |
|
|
208 @deftypefn {Function File} {} trace (@var{a}) |
|
|
209 Compute the trace of @var{a}, @code{sum (diag (@var{a}))}. |
|
|
210 @end deftypefn |
|
|
211 |
|
|
212 @node Matrix Factorizations, Functions of a Matrix, Basic Matrix Functions, Linear Algebra |
|
|
213 @section Matrix Factorizations |
|
|
214 |
|
|
215 @deftypefn {Loadable Function} {} chol (@var{a}) |
|
|
216 @cindex Cholesky factorization |
|
|
217 Compute the Cholesky factor, @var{r}, of the symmetric positive definite |
|
|
218 matrix @var{a}, where |
|
|
219 @iftex |
|
|
220 @tex |
|
|
221 $ R^T R = A $. |
|
|
222 @end tex |
|
|
223 @end iftex |
|
|
224 @ifinfo |
|
|
225 |
|
|
226 @example |
|
|
227 r' * r = a. |
|
|
228 @end example |
|
|
229 @end ifinfo |
|
|
230 @end deftypefn |
|
|
231 |
|
|
232 @deftypefn {Loadable Function} {@var{h} =} hess (@var{a}) |
|
|
233 @deftypefnx {Loadable Function} {[@var{p}, @var{h}] =} hess (@var{a}) |
|
|
234 @cindex Hessenberg decomposition |
|
|
235 Compute the Hessenberg decomposition of the matrix @var{a}. |
|
|
236 |
|
|
237 The Hessenberg decomposition is usually used as the first step in an |
|
|
238 eigenvalue computation, but has other applications as well (see Golub, |
|
|
239 Nash, and Van Loan, IEEE Transactions on Automatic Control, 1979. The |
|
|
240 Hessenberg decomposition is |
|
|
241 @iftex |
|
|
242 @tex |
|
|
243 $$ |
|
|
244 A = PHP^T |
|
|
245 $$ |
|
|
246 where $P$ is a square unitary matrix ($P^HP = I$), and $H$ |
|
|
247 is upper Hessenberg ($H_{i,j} = 0, \forall i \ge j+1$). |
|
|
248 @end tex |
|
|
249 @end iftex |
|
|
250 @ifinfo |
|
|
251 @code{p * h * p' = a} where @code{p} is a square unitary matrix |
|
|
252 (@code{p' * p = I}, using complex-conjugate transposition) and @code{h} |
|
|
253 is upper Hessenberg (@code{i >= j+1 => h (i, j) = 0}). |
|
|
254 @end ifinfo |
|
|
255 @end deftypefn |
|
|
256 |
|
|
257 @deftypefn {Loadable Function} {[@var{l}, @var{u}, @var{p}] =} lu (@var{a}) |
|
|
258 @cindex LU decomposition |
|
|
259 Compute the LU decomposition of @var{a}, using subroutines from |
|
|
260 @sc{Lapack}. The result is returned in a permuted form, according to |
|
|
261 the optional return value @var{p}. For example, given the matrix |
|
|
262 @code{a = [1, 2; 3, 4]}, |
|
|
263 |
|
|
264 @example |
|
|
265 [l, u, p] = lu (a) |
|
|
266 @end example |
|
|
267 |
|
|
268 @noindent |
|
|
269 returns |
|
|
270 |
|
|
271 @example |
|
|
272 l = |
|
|
273 |
|
|
274 1.00000 0.00000 |
|
|
275 0.33333 1.00000 |
|
|
276 |
|
|
277 u = |
|
|
278 |
|
|
279 3.00000 4.00000 |
|
|
280 0.00000 0.66667 |
|
|
281 |
|
|
282 p = |
|
|
283 |
|
|
284 0 1 |
|
|
285 1 0 |
|
|
286 @end example |
|
|
287 @end deftypefn |
|
|
288 |
|
|
289 @deftypefn {Loadable Function} {[@var{q}, @var{r}, @var{p}] =} qr (@var{a}) |
|
|
290 @cindex QR factorization |
|
|
291 Compute the QR factorization of @var{a}, using standard @sc{Lapack} |
|
|
292 subroutines. For example, given the matrix @code{a = [1, 2; 3, 4]}, |
|
|
293 |
|
|
294 @example |
|
|
295 [q, r] = qr (a) |
|
|
296 @end example |
|
|
297 |
|
|
298 @noindent |
|
|
299 returns |
|
|
300 |
|
|
301 @example |
|
|
302 q = |
|
|
303 |
|
|
304 -0.31623 -0.94868 |
|
|
305 -0.94868 0.31623 |
|
|
306 |
|
|
307 r = |
|
|
308 |
|
|
309 -3.16228 -4.42719 |
|
|
310 0.00000 -0.63246 |
|
|
311 @end example |
|
|
312 |
|
|
313 The @code{qr} factorization has applications in the solution of least |
|
|
314 squares problems |
|
|
315 @iftex |
|
|
316 @tex |
|
|
317 $$ |
|
|
318 \min_x \left\Vert A x - b \right\Vert_2 |
|
|
319 $$ |
|
|
320 @end tex |
|
|
321 @end iftex |
|
|
322 @ifinfo |
|
|
323 |
|
|
324 @example |
|
|
325 @code{min norm(A x - b)} |
|
|
326 @end example |
|
|
327 |
|
|
328 @end ifinfo |
|
|
329 for overdetermined systems of equations (i.e., |
|
|
330 @iftex |
|
|
331 @tex |
|
|
332 $A$ |
|
|
333 @end tex |
|
|
334 @end iftex |
|
|
335 @ifinfo |
|
|
336 @code{a} |
|
|
337 @end ifinfo |
|
|
338 is a tall, thin matrix). The QR factorization is |
|
|
339 @iftex |
|
|
340 @tex |
|
|
341 $QR = A$ where $Q$ is an orthogonal matrix and $R$ is upper triangular. |
|
|
342 @end tex |
|
|
343 @end iftex |
|
|
344 @ifinfo |
|
|
345 @code{q * r = a} where @code{q} is an orthogonal matrix and @code{r} is |
|
|
346 upper triangular. |
|
|
347 @end ifinfo |
|
|
348 |
|
|
349 The permuted QR factorization @code{[@var{q}, @var{r}, @var{p}] = |
|
|
350 qr (@var{a})} forms the QR factorization such that the diagonal |
|
|
351 entries of @code{r} are decreasing in magnitude order. For example, |
|
|
352 given the matrix @code{a = [1, 2; 3, 4]}, |
|
|
353 |
|
|
354 @example |
|
|
355 [q, r, pi] = qr(a) |
|
|
356 @end example |
|
|
357 |
|
|
358 @noindent |
|
|
359 returns |
|
|
360 |
|
|
361 @example |
|
|
362 q = |
|
|
363 |
|
|
364 -0.44721 -0.89443 |
|
|
365 -0.89443 0.44721 |
|
|
366 |
|
|
367 r = |
|
|
368 |
|
|
369 -4.47214 -3.13050 |
|
|
370 0.00000 0.44721 |
|
|
371 |
|
|
372 p = |
|
|
373 |
|
|
374 0 1 |
|
|
375 1 0 |
|
|
376 @end example |
|
|
377 |
|
|
378 The permuted @code{qr} factorization @code{[q, r, p] = qr (a)} |
|
|
379 factorization allows the construction of an orthogonal basis of |
|
|
380 @code{span (a)}. |
|
|
381 @end deftypefn |
|
|
382 |
|
|
383 |
|
|
384 @deftypefn {Function File} {@var{lambda} =} qz (@var{a}, @var{b}) |
|
|
385 Generalized eigenvalue problem @math{A x = s B x}, |
|
|
386 @var{QZ} decomposition. Three ways to call: |
|
|
387 @enumerate |
|
|
388 @item @code{lambda = qz(A,B)} |
|
|
389 |
|
|
390 Computes the generalized eigenvalues @var{lambda} of @math{(A - sB)}. |
|
|
391 |
|
|
392 @item @code{[AA, BB, Q, Z @{, V, W, lambda@}] = qz (A, B)} |
|
|
393 |
|
|
394 Computes qz decomposition, generalized eigenvectors, and |
|
|
395 generalized eigenvalues of @math{(A - sB)} |
|
|
396 @example |
|
|
397 @group |
|
|
398 A V = B V diag(lambda) |
|
|
399 W' A = diag(lambda) W' B |
|
|
400 AA = Q'*A*Z, BB = Q'*B*Z with Q, Z orthogonal (unitary)= I |
|
|
401 @end group |
|
|
402 @end example |
|
|
403 |
|
|
404 @item @code{[AA,BB,Z@{,lambda@}] = qz(A,B,opt)} |
|
|
405 |
|
|
406 As in form [2], but allows ordering of generalized eigenpairs |
|
|
407 for (e.g.) solution of discrete time algebraic Riccati equations. |
|
|
408 Form 3 is not available for complex matrices and does not compute |
|
|
409 the generalized eigenvectors V, W, nor the orthogonal matrix Q. |
|
|
410 @table @var |
|
|
411 @item opt |
|
|
412 for ordering eigenvalues of the GEP pencil. The leading block |
|
|
413 of the revised pencil contains all eigenvalues that satisfy: |
|
|
414 @table @code |
|
|
415 @item "N" |
|
|
416 = unordered (default) |
|
|
417 |
|
|
418 @item "S" |
|
|
419 = small: leading block has all |lambda| <=1 |
|
|
420 |
|
|
421 @item "B" |
|
|
422 = big: leading block has all |lambda >= 1 |
|
|
423 |
|
|
424 @item "-" |
|
|
425 = negative real part: leading block has all eigenvalues |
|
|
426 in the open left half-plant |
|
|
427 |
|
|
428 @item "+" |
|
|
429 = nonnegative real part: leading block has all eigenvalues |
|
|
430 in the closed right half-plane |
|
|
431 @end table |
|
|
432 @end table |
|
|
433 @end enumerate |
|
|
434 |
|
|
435 Note: qz performs permutation balancing, but not scaling (see balance). |
|
|
436 Order of output arguments was selected for compatibility with MATLAB |
|
|
437 |
|
|
438 See also: balance, dare, eig, schur |
|
|
439 @end deftypefn |
|
|
440 |
|
|
441 @deftypefn {Function File} {[@var{aa}, @var{bb}, @var{q}, @var{z}] =} qzhess (@var{a}, @var{b}) |
|
|
442 Compute the Hessenberg-triangular decomposition of the matrix pencil |
|
|
443 @code{(@var{a}, @var{b})}, returning |
|
|
444 @code{@var{aa} = @var{q} * @var{a} * @var{z}}, |
|
|
445 @code{@var{bb} = @var{q} * @var{b} * @var{z}}, with @var{q} and @var{z} |
|
|
446 orthogonal. For example, |
|
|
447 |
|
|
448 @example |
|
|
449 @group |
|
|
450 [aa, bb, q, z] = qzhess ([1, 2; 3, 4], [5, 6; 7, 8]) |
|
|
451 @result{} aa = [ -3.02244, -4.41741; 0.92998, 0.69749 ] |
|
|
452 @result{} bb = [ -8.60233, -9.99730; 0.00000, -0.23250 ] |
|
|
453 @result{} q = [ -0.58124, -0.81373; -0.81373, 0.58124 ] |
|
|
454 @result{} z = [ 1, 0; 0, 1 ] |
|
|
455 @end group |
|
|
456 @end example |
|
|
457 |
|
|
458 The Hessenberg-triangular decomposition is the first step in |
|
|
459 Moler and Stewart's QZ decomposition algorithm. |
|
|
460 |
|
|
461 Algorithm taken from Golub and Van Loan, @cite{Matrix Computations, 2nd |
|
|
462 edition}. |
|
|
463 @end deftypefn |
|
|
464 |
|
|
465 @deftypefn {Loadable Function} {} qzval (@var{a}, @var{b}) |
|
|
466 Compute generalized eigenvalues of the matrix pencil |
|
|
467 @iftex |
|
|
468 @tex |
|
|
469 $a - \lambda b$. |
|
|
470 @end tex |
|
|
471 @end iftex |
|
|
472 @ifinfo |
|
|
473 @code{@var{a} - lambda @var{b}}. |
|
|
474 @end ifinfo |
|
|
475 |
|
|
476 The arguments @var{a} and @var{b} must be real matrices. |
|
|
477 @end deftypefn |
|
|
478 |
|
|
479 @deftypefn {Loadable Function} {@var{s} =} schur (@var{a}) |
|
|
480 @deftypefnx {Loadable Function} {[@var{u}, @var{s}] =} schur (@var{a}, @var{opt}) |
|
|
481 @cindex Schur decomposition |
|
|
482 The Schur decomposition is used to compute eigenvalues of a |
|
|
483 square matrix, and has applications in the solution of algebraic |
|
|
484 Riccati equations in control (see @code{are} and @code{dare}). |
|
|
485 @code{schur} always returns |
|
|
486 @iftex |
|
|
487 @tex |
|
|
488 $S = U^T A U$ |
|
|
489 @end tex |
|
|
490 @end iftex |
|
|
491 @ifinfo |
|
|
492 @code{s = u' * a * u} |
|
|
493 @end ifinfo |
|
|
494 where |
|
|
495 @iftex |
|
|
496 @tex |
|
|
497 $U$ |
|
|
498 @end tex |
|
|
499 @end iftex |
|
|
500 @ifinfo |
|
|
501 @code{u} |
|
|
502 @end ifinfo |
|
|
503 is a unitary matrix |
|
|
504 @iftex |
|
|
505 @tex |
|
|
506 ($U^T U$ is identity) |
|
|
507 @end tex |
|
|
508 @end iftex |
|
|
509 @ifinfo |
|
|
510 (@code{u'* u} is identity) |
|
|
511 @end ifinfo |
|
|
512 and |
|
|
513 @iftex |
|
|
514 @tex |
|
|
515 $S$ |
|
|
516 @end tex |
|
|
517 @end iftex |
|
|
518 @ifinfo |
|
|
519 @code{s} |
|
|
520 @end ifinfo |
|
|
521 is upper triangular. The eigenvalues of |
|
|
522 @iftex |
|
|
523 @tex |
|
|
524 $A$ (and $S$) |
|
|
525 @end tex |
|
|
526 @end iftex |
|
|
527 @ifinfo |
|
|
528 @code{a} (and @code{s}) |
|
|
529 @end ifinfo |
|
|
530 are the diagonal elements of |
|
|
531 @iftex |
|
|
532 @tex |
|
|
533 $S$ |
|
|
534 @end tex |
|
|
535 @end iftex |
|
|
536 @ifinfo |
|
|
537 @code{s} |
|
|
538 @end ifinfo |
|
|
539 If the matrix |
|
|
540 @iftex |
|
|
541 @tex |
|
|
542 $A$ |
|
|
543 @end tex |
|
|
544 @end iftex |
|
|
545 @ifinfo |
|
|
546 @code{a} |
|
|
547 @end ifinfo |
|
|
548 is real, then the real Schur decomposition is computed, in which the |
|
|
549 matrix |
|
|
550 @iftex |
|
|
551 @tex |
|
|
552 $U$ |
|
|
553 @end tex |
|
|
554 @end iftex |
|
|
555 @ifinfo |
|
|
556 @code{u} |
|
|
557 @end ifinfo |
|
|
558 is orthogonal and |
|
|
559 @iftex |
|
|
560 @tex |
|
|
561 $S$ |
|
|
562 @end tex |
|
|
563 @end iftex |
|
|
564 @ifinfo |
|
|
565 @code{s} |
|
|
566 @end ifinfo |
|
|
567 is block upper triangular |
|
|
568 with blocks of size at most |
|
|
569 @iftex |
|
|
570 @tex |
|
|
571 $2\times 2$ |
|
|
572 @end tex |
|
|
573 @end iftex |
|
|
574 @ifinfo |
|
|
575 @code{2 x 2} |
|
|
576 @end ifinfo |
|
|
577 blocks along the diagonal. The diagonal elements of |
|
|
578 @iftex |
|
|
579 @tex |
|
|
580 $S$ |
|
|
581 @end tex |
|
|
582 @end iftex |
|
|
583 @ifinfo |
|
|
584 @code{s} |
|
|
585 @end ifinfo |
|
|
586 (or the eigenvalues of the |
|
|
587 @iftex |
|
|
588 @tex |
|
|
589 $2\times 2$ |
|
|
590 @end tex |
|
|
591 @end iftex |
|
|
592 @ifinfo |
|
|
593 @code{2 x 2} |
|
|
594 @end ifinfo |
|
|
595 blocks, when |
|
|
596 appropriate) are the eigenvalues of |
|
|
597 @iftex |
|
|
598 @tex |
|
|
599 $A$ |
|
|
600 @end tex |
|
|
601 @end iftex |
|
|
602 @ifinfo |
|
|
603 @code{a} |
|
|
604 @end ifinfo |
|
|
605 and |
|
|
606 @iftex |
|
|
607 @tex |
|
|
608 $S$. |
|
|
609 @end tex |
|
|
610 @end iftex |
|
|
611 @ifinfo |
|
|
612 @code{s}. |
|
|
613 @end ifinfo |
|
|
614 |
|
|
615 The eigenvalues are optionally ordered along the diagonal according to |
|
|
616 the value of @code{opt}. @code{opt = "a"} indicates that all |
|
|
617 eigenvalues with negative real parts should be moved to the leading |
|
|
618 block of |
|
|
619 @iftex |
|
|
620 @tex |
|
|
621 $S$ |
|
|
622 @end tex |
|
|
623 @end iftex |
|
|
624 @ifinfo |
|
|
625 @code{s} |
|
|
626 @end ifinfo |
|
|
627 (used in @code{are}), @code{opt = "d"} indicates that all eigenvalues |
|
|
628 with magnitude less than one should be moved to the leading block of |
|
|
629 @iftex |
|
|
630 @tex |
|
|
631 $S$ |
|
|
632 @end tex |
|
|
633 @end iftex |
|
|
634 @ifinfo |
|
|
635 @code{s} |
|
|
636 @end ifinfo |
|
|
637 (used in @code{dare}), and @code{opt = "u"}, the default, indicates that |
|
|
638 no ordering of eigenvalues should occur. The leading |
|
|
639 @iftex |
|
|
640 @tex |
|
|
641 $k$ |
|
|
642 @end tex |
|
|
643 @end iftex |
|
|
644 @ifinfo |
|
|
645 @code{k} |
|
|
646 @end ifinfo |
|
|
647 columns of |
|
|
648 @iftex |
|
|
649 @tex |
|
|
650 $U$ |
|
|
651 @end tex |
|
|
652 @end iftex |
|
|
653 @ifinfo |
|
|
654 @code{u} |
|
|
655 @end ifinfo |
|
|
656 always span the |
|
|
657 @iftex |
|
|
658 @tex |
|
|
659 $A$-invariant |
|
|
660 @end tex |
|
|
661 @end iftex |
|
|
662 @ifinfo |
|
|
663 @code{a}-invariant |
|
|
664 @end ifinfo |
|
|
665 subspace corresponding to the |
|
|
666 @iftex |
|
|
667 @tex |
|
|
668 $k$ |
|
|
669 @end tex |
|
|
670 @end iftex |
|
|
671 @ifinfo |
|
|
672 @code{k} |
|
|
673 @end ifinfo |
|
|
674 leading eigenvalues of |
|
|
675 @iftex |
|
|
676 @tex |
|
|
677 $S$. |
|
|
678 @end tex |
|
|
679 @end iftex |
|
|
680 @ifinfo |
|
|
681 @code{s}. |
|
|
682 @end ifinfo |
|
|
683 @end deftypefn |
|
|
684 |
|
|
685 @deftypefn {Loadable Function} {@var{s} =} svd (@var{a}) |
|
|
686 @deftypefnx {Loadable Function} {[@var{u}, @var{s}, @var{v}] =} svd (@var{a}) |
|
|
687 @cindex singular value decomposition |
|
|
688 Compute the singular value decomposition of @var{a} |
|
|
689 @iftex |
|
|
690 @tex |
|
|
691 $$ |
|
|
692 A = U\Sigma V^H |
|
|
693 $$ |
|
|
694 @end tex |
|
|
695 @end iftex |
|
|
696 @ifinfo |
|
|
697 |
|
|
698 @example |
|
|
699 a = u * sigma * v' |
|
|
700 @end example |
|
|
701 @end ifinfo |
|
|
702 |
|
|
703 The function @code{svd} normally returns the vector of singular values. |
|
|
704 If asked for three return values, it computes |
|
|
705 @iftex |
|
|
706 @tex |
|
|
707 $U$, $S$, and $V$. |
|
|
708 @end tex |
|
|
709 @end iftex |
|
|
710 @ifinfo |
|
|
711 U, S, and V. |
|
|
712 @end ifinfo |
|
|
713 For example, |
|
|
714 |
|
|
715 @example |
|
|
716 svd (hilb (3)) |
|
|
717 @end example |
|
|
718 |
|
|
719 @noindent |
|
|
720 returns |
|
|
721 |
|
|
722 @example |
|
|
723 ans = |
|
|
724 |
|
|
725 1.4083189 |
|
|
726 0.1223271 |
|
|
727 0.0026873 |
|
|
728 @end example |
|
|
729 |
|
|
730 @noindent |
|
|
731 and |
|
|
732 |
|
|
733 @example |
|
|
734 [u, s, v] = svd (hilb (3)) |
|
|
735 @end example |
|
|
736 |
|
|
737 @noindent |
|
|
738 returns |
|
|
739 |
|
|
740 @example |
|
|
741 u = |
|
|
742 |
|
|
743 -0.82704 0.54745 0.12766 |
|
|
744 -0.45986 -0.52829 -0.71375 |
|
|
745 -0.32330 -0.64901 0.68867 |
|
|
746 |
|
|
747 s = |
|
|
748 |
|
|
749 1.40832 0.00000 0.00000 |
|
|
750 0.00000 0.12233 0.00000 |
|
|
751 0.00000 0.00000 0.00269 |
|
|
752 |
|
|
753 v = |
|
|
754 |
|
|
755 -0.82704 0.54745 0.12766 |
|
|
756 -0.45986 -0.52829 -0.71375 |
|
|
757 -0.32330 -0.64901 0.68867 |
|
|
758 @end example |
|
|
759 |
|
|
760 If given a second argument, @code{svd} returns an economy-sized |
|
|
761 decomposition, eliminating the unnecessary rows or columns of @var{u} or |
|
|
762 @var{v}. |
|
|
763 @end deftypefn |
|
|
764 |
|
|
765 @node Functions of a Matrix, , Matrix Factorizations, Linear Algebra |
|
|
766 @section Functions of a Matrix |
|
|
767 |
|
|
768 @deftypefn {Loadable Function} {} expm (@var{a}) |
|
|
769 Return the exponential of a matrix, defined as the infinite Taylor |
|
|
770 series |
|
|
771 @iftex |
|
|
772 @tex |
|
|
773 $$ |
|
|
774 \exp (A) = I + A + {A^2 \over 2!} + {A^3 \over 3!} + \cdots |
|
|
775 $$ |
|
|
776 @end tex |
|
|
777 @end iftex |
|
|
778 @ifinfo |
|
|
779 |
|
|
780 @example |
|
|
781 expm(a) = I + a + a^2/2! + a^3/3! + ... |
|
|
782 @end example |
|
|
783 |
|
|
784 @end ifinfo |
|
|
785 The Taylor series is @emph{not} the way to compute the matrix |
|
|
786 exponential; see Moler and Van Loan, @cite{Nineteen Dubious Ways to |
|
|
787 Compute the Exponential of a Matrix}, SIAM Review, 1978. This routine |
|
|
788 uses Ward's diagonal |
|
|
789 @iftex |
|
|
790 @tex |
|
|
791 Pad\'e |
|
|
792 @end tex |
|
|
793 @end iftex |
|
|
794 @ifinfo |
|
|
795 Pade' |
|
|
796 @end ifinfo |
|
|
797 approximation method with three step preconditioning (SIAM Journal on |
|
|
798 Numerical Analysis, 1977). Diagonal |
|
|
799 @iftex |
|
|
800 @tex |
|
|
801 Pad\'e |
|
|
802 @end tex |
|
|
803 @end iftex |
|
|
804 @ifinfo |
|
|
805 Pade' |
|
|
806 @end ifinfo |
|
|
807 approximations are rational polynomials of matrices |
|
|
808 @iftex |
|
|
809 @tex |
|
|
810 $D_q(a)^{-1}N_q(a)$ |
|
|
811 @end tex |
|
|
812 @end iftex |
|
|
813 @ifinfo |
|
|
814 |
|
|
815 @example |
|
|
816 -1 |
|
|
817 D (a) N (a) |
|
|
818 @end example |
|
|
819 |
|
|
820 @end ifinfo |
|
|
821 whose Taylor series matches the first |
|
|
822 @iftex |
|
|
823 @tex |
|
|
824 $2 q + 1 $ |
|
|
825 @end tex |
|
|
826 @end iftex |
|
|
827 @ifinfo |
|
|
828 @code{2q+1} |
|
|
829 @end ifinfo |
|
|
830 terms of the Taylor series above; direct evaluation of the Taylor series |
|
|
831 (with the same preconditioning steps) may be desirable in lieu of the |
|
|
832 @iftex |
|
|
833 @tex |
|
|
834 Pad\'e |
|
|
835 @end tex |
|
|
836 @end iftex |
|
|
837 @ifinfo |
|
|
838 Pade' |
|
|
839 @end ifinfo |
|
|
840 approximation when |
|
|
841 @iftex |
|
|
842 @tex |
|
|
843 $D_q(a)$ |
|
|
844 @end tex |
|
|
845 @end iftex |
|
|
846 @ifinfo |
|
|
847 @code{Dq(a)} |
|
|
848 @end ifinfo |
|
|
849 is ill-conditioned. |
|
|
850 @end deftypefn |
|
|
851 |
|
|
852 @deftypefn {Loadable Function} {} logm (@var{a}) |
|
|
853 Compute the matrix logarithm of the square matrix @var{a}. Note that |
|
|
854 this is currently implemented in terms of an eigenvalue expansion and |
|
|
855 needs to be improved to be more robust. |
|
|
856 @end deftypefn |
|
|
857 |
|
|
858 @deftypefn {Loadable Function} {} sqrtm (@var{a}) |
|
|
859 Compute the matrix square root of the square matrix @var{a}. Note that |
|
|
860 this is currently implemented in terms of an eigenvalue expansion and |
|
|
861 needs to be improved to be more robust. |
|
|
862 @end deftypefn |
|
|
863 |
|
|
864 @deftypefn {Function File} {} kron (@var{a}, @var{b}) |
|
|
865 Form the kronecker product of two matrices, defined block by block as |
|
|
866 |
|
|
867 @example |
|
|
868 x = [a(i, j) b] |
|
|
869 @end example |
|
|
870 |
|
|
871 For example, |
|
|
872 |
|
|
873 @example |
|
|
874 @group |
|
|
875 kron (1:4, ones (3, 1)) |
|
|
876 @result{} 1 2 3 4 |
|
|
877 1 2 3 4 |
|
|
878 1 2 3 4 |
|
|
879 @end group |
|
|
880 @end example |
|
|
881 @end deftypefn |
|
|
882 |
|
|
883 @deftypefn {Loadable Function} {@var{x} =} syl (@var{a}, @var{b}, @var{c}) |
|
|
884 Solve the Sylvester equation |
|
|
885 @iftex |
|
|
886 @tex |
|
|
887 $$ |
|
|
888 A X + X B + C = 0 |
|
|
889 $$ |
|
|
890 @end tex |
|
|
891 @end iftex |
|
|
892 @ifinfo |
|
|
893 |
|
|
894 @example |
|
|
895 A X + X B + C = 0 |
|
|
896 @end example |
|
|
897 @end ifinfo |
|
|
898 using standard @sc{Lapack} subroutines. For example, |
|
|
899 |
|
|
900 @example |
|
|
901 @group |
|
|
902 syl ([1, 2; 3, 4], [5, 6; 7, 8], [9, 10; 11, 12]) |
|
|
903 @result{} [ -0.50000, -0.66667; -0.66667, -0.50000 ] |
|
|
904 @end group |
|
|
905 @end example |
|
|
906 @end deftypefn |
