Mercurial > octave-libgccjit

comparison doc/interpreter/diagperm.txi @ 8917:d707aa3bbc36

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
manual improvements
author Jaroslav Hajek <highegg@gmail.com>
date Fri, 06 Mar 2009 13:42:24 +0100
parents 46fdf8714acf
children 349616d9c38e
comparison
equal deleted inserted replaced
8916:a2878ba31a9e 8917:d707aa3bbc36
26 * Example Codes:: Some Examples of Usage 26 * Example Codes:: Some Examples of Usage
27 * Zeros Treatment:: The Differences in Treatment of Zero Elements 27 * Zeros Treatment:: The Differences in Treatment of Zero Elements
28 @end menu 28 @end menu
29 29
30 @node Basic Usage 30 @node Basic Usage
31 @section The Creation and Manipulation of Diagonal and Permutation Matrices 31 @section Creating and Manipulating Diagonal and Permutation Matrices
32 32
33 A diagonal matrix is defined as a matrix that has zero entries outside the main diagonal; 33 A diagonal matrix is defined as a matrix that has zero entries outside the main diagonal;
34 that is, 34 that is,
35 @iftex 35 @iftex
36 @tex
36 $D_{ij} = 0$ if $i \neq j$ 37 $D_{ij} = 0$ if $i \neq j$
38 @end tex
37 @end iftex 39 @end iftex
38 @ifnottex 40 @ifnottex
39 @code{D(i,j) == 0} if @code{i != j}. 41 @code{D(i,j) == 0} if @code{i != j}.
40 @end ifnottex 42 @end ifnottex
41 Most often, square diagonal matrices are considered; however, the definition can equally 43 Most often, square diagonal matrices are considered; however, the definition can equally
44 46
45 A permutation matrix is defined as a square matrix that has a single element equal to unity 47 A permutation matrix is defined as a square matrix that has a single element equal to unity
46 in each row and each column; all other elements are zero. That is, there exists a 48 in each row and each column; all other elements are zero. That is, there exists a
47 permutation (vector) 49 permutation (vector)
48 @iftex 50 @iftex
51 @tex
49 $p$ such that $P_{ij}=1$ if $j = p_i$ and 52 $p$ such that $P_{ij}=1$ if $j = p_i$ and
50 $P_{ij}=0$ otherwise. 53 $P_{ij}=0$ otherwise.
54 @end tex
51 @end iftex 55 @end iftex
52 @ifnottex 56 @ifnottex
53 @code{p} such that @code{P(i,j) == 1} if @code{j == p(i)} and 57 @code{p} such that @code{P(i,j) == 1} if @code{j == p(i)} and
54 @code{P(i,j) == 0} otherwise. 58 @code{P(i,j) == 0} otherwise.
55 @end ifnottex 59 @end ifnottex
67 71
68 @node Creating Diagonal Matrices 72 @node Creating Diagonal Matrices
69 @subsection Creating Diagonal Matrices 73 @subsection Creating Diagonal Matrices
70 74
71 The most common and easiest way to create a diagonal matrix is using the built-in 75 The most common and easiest way to create a diagonal matrix is using the built-in
72 function @dfn{diag}. The expression @code{diag (@var{v})}, with @var{v} a vector, 76 function @dfn{diag}. The expression @code{diag (v)}, with @var{v} a vector,
73 will create a square diagonal matrix with elements on the main diagonal given 77 will create a square diagonal matrix with elements on the main diagonal given
74 by the elements of @var{v}, and size equal to the length of @var{v}. 78 by the elements of @var{v}, and size equal to the length of @var{v}.
75 @code{diag (@var{v}, m, n)} can be used to construct a rectangular diagonal matrix. 79 @code{diag (v, m, n)} can be used to construct a rectangular diagonal matrix.
76 The result of these expressions will be a special diagonal matrix object, rather 80 The result of these expressions will be a special diagonal matrix object, rather
77 than a general matrix object. 81 than a general matrix object.
78 82
83 Diagonal matrix with unit elements can be created using @dfn{eye}.
79 Some other built-in functions can also return diagonal matrices. Examples include 84 Some other built-in functions can also return diagonal matrices. Examples include
80 @dfn{balance} or @dfn{inv}. 85 @dfn{balance} or @dfn{inv}.
86
87 Example:
88 @example
89 diag (1:4)
90 @result{}
91 Diagonal Matrix
92
93 1 0 0 0
94 0 2 0 0
95 0 0 3 0
96 0 0 0 4
97
98 diag(1:3,5,3)
99
100 @result{}
101 Diagonal Matrix
102
103 1 0 0
104 0 2 0
105 0 0 3
106 0 0 0
107 0 0 0
108 @end example
81 109
82 @node Creating Permutation Matrices 110 @node Creating Permutation Matrices
83 @subsection Creating Permutation Matrices 111 @subsection Creating Permutation Matrices
84 112
85 For creating permutation matrices, Octave does not introduce a new function, but 113 For creating permutation matrices, Octave does not introduce a new function, but
86 rather overrides an existing syntax: permutation matrices can be conveniently 114 rather overrides an existing syntax: permutation matrices can be conveniently
87 created by indexing an identity matrix by permutation vectors. 115 created by indexing an identity matrix by permutation vectors.
88 That is, if @var{q} is a permutation vector of length @var{n}, the expression 116 That is, if @var{q} is a permutation vector of length @var{n}, the expression
89 @example 117 @example
90 @var{P} = eye (@var{n}) (:, @var{q}); 118 P = eye (n) (:, q);
91 @end example 119 @end example
92 will create a permutation matrix - a special matrix object. 120 will create a permutation matrix - a special matrix object.
93 @example 121 @example
94 eye (@var{n}) (@var{q}, :) 122 eye (n) (q, :)
95 @end example 123 @end example
96 will also work (and create a row permutation matrix), as well as 124 will also work (and create a row permutation matrix), as well as
97 @example 125 @example
98 eye (@var{n}) (@var{q1}, @var{q2}). 126 eye (n) (q1, q2).
127 @end example
128
129 For example:
130 @example
131 eye (4) ([1,3,2,4],:)
132 @result{}
133 Permutation Matrix
134
135 1 0 0 0
136 0 0 1 0
137 0 1 0 0
138 0 0 0 1
139
140 eye (4) (:,[1,3,2,4])
141 @result{}
142 Permutation Matrix
143
144 1 0 0 0
145 0 0 1 0
146 0 1 0 0
147 0 0 0 1
99 @end example 148 @end example
100 149
101 Mathematically, an identity matrix is both diagonal and permutation matrix. 150 Mathematically, an identity matrix is both diagonal and permutation matrix.
102 In Octave, @code{eye (@var{n})} returns a diagonal matrix, because a matrix 151 In Octave, @code{eye (n)} returns a diagonal matrix, because a matrix
103 can only have one class. You can convert this diagonal matrix to a permutation 152 can only have one class. You can convert this diagonal matrix to a permutation
104 matrix by indexing it by an identity permutation, as shown above. 153 matrix by indexing it by an identity permutation, as shown below.
105 This is a special property of the identity matrix; indexing other diagonal 154 This is a special property of the identity matrix; indexing other diagonal
106 matrices generall produces a full matrix. 155 matrices generally produces a full matrix.
156
157 @example
158 eye (3)
159 @result{}
160 Diagonal Matrix
161
162 1 0 0
163 0 1 0
164 0 0 1
165
166 eye(3)(1:3,:)
167 @result{}
168 Permutation Matrix
169
170 1 0 0
171 0 1 0
172 0 0 1
173 @end example
107 174
108 Some other built-in functions can also return permutation matrices. Examples include 175 Some other built-in functions can also return permutation matrices. Examples include
109 @dfn{inv} or @dfn{lu}. 176 @dfn{inv} or @dfn{lu}.
110 177
111 @node Explicit and Implicit Conversions 178 @node Explicit and Implicit Conversions
112 @subsection Explicit and Implicit Conversions 179 @subsection Explicit and Implicit Conversions
113 180
114 The diagonal and permutation matrices are special objects in their own right. A number 181 The diagonal and permutation matrices are special objects in their own right. A number
115 of operations and built-in functions, are defined for these matrices to use special, 182 of operations and built-in functions are defined for these matrices to use special,
116 more efficient code than would be used for a full matrix in the same place. Examples 183 more efficient code than would be used for a full matrix in the same place. Examples
117 are given in further sections. 184 are given in further sections.
118 185
119 To facilitate smooth mixing with full matrices, backward compatibility, and 186 To facilitate smooth mixing with full matrices, backward compatibility, and
120 compatibility with Matlab, the diagonal and permutation matrices should allow 187 compatibility with Matlab, the diagonal and permutation matrices should allow
121 any operation that works on full matrices, and will either treat is specially, 188 any operation that works on full matrices, and will either treat it specially,
122 or implicitly convert themselves to full matrices. 189 or implicitly convert themselves to full matrices.
123 190
124 Instances include matrix indexing (except for extracting a single element), 191 Instances include matrix indexing, except for extracting a single element or
125 indexed assignment, or applying most mapper functions, such as @dfn{exp}. 192 a leading submatrix, indexed assignment, or applying most mapper functions,
193 such as @dfn{exp}.
126 194
127 An explicit conversion to a full matrix can be requested using the built-in 195 An explicit conversion to a full matrix can be requested using the built-in
128 function @dfn{full}. It should also be noted that the diagonal and permutation 196 function @dfn{full}. It should also be noted that the diagonal and permutation
129 matrix objects will cache the result of the conversion after it is first 197 matrix objects will cache the result of the conversion after it is first
130 requested (explicitly or implicitly), so that subsequent conversions will 198 requested (explicitly or implicitly), so that subsequent conversions will
145 213
146 @node Expressions Involving Diagonal Matrices 214 @node Expressions Involving Diagonal Matrices
147 @subsection Expressions Involving Diagonal Matrices 215 @subsection Expressions Involving Diagonal Matrices
148 216
149 Assume @var{D} is a diagonal matrix. If @var{M} is a full matrix, 217 Assume @var{D} is a diagonal matrix. If @var{M} is a full matrix,
150 then @code{@var{D}*@var{M}} will scale the rows of @var{M}. That means, 218 then @code{D*M} will scale the rows of @var{M}. That means,
151 if @code{@var{S} = @var{D}*@var{M}}, then for each pair of indices 219 if @code{S = D*M}, then for each pair of indices
152 i,j it holds 220 i,j it holds
153 @iftex 221 @iftex
222 @tex
154 $$S_{ij} = D_{ii} M_{ij}$$ 223 $$S_{ij} = D_{ii} M_{ij}$$
224 @end tex
155 @end iftex 225 @end iftex
156 @ifnottex 226 @ifnottex
157 @example 227 @example
158 S(i,j) = D(i,i) * M(i,j). 228 S(i,j) = D(i,i) * M(i,j).
159 @end example 229 @end example
176 The expressions @code{D \ M} and @code{M / D} perform inverse scaling. 246 The expressions @code{D \ M} and @code{M / D} perform inverse scaling.
177 They are equivalent to solving a diagonal (or rectangular diagonal) 247 They are equivalent to solving a diagonal (or rectangular diagonal)
178 in a least-squares minimum-norm sense. In exact arithmetics, this is 248 in a least-squares minimum-norm sense. In exact arithmetics, this is
179 equivalent to multiplying by a pseudoinverse. The pseudoinverse of 249 equivalent to multiplying by a pseudoinverse. The pseudoinverse of
180 a rectangular diagonal matrix is again a rectangular diagonal matrix 250 a rectangular diagonal matrix is again a rectangular diagonal matrix
181 of the same dimensions, where each nonzero diagonal element is replaced 251 with swapped dimensions, where each nonzero diagonal element is replaced
182 by its reciprocal. 252 by its reciprocal.
183 The matrix division algorithms do, in fact, use division rather than 253 The matrix division algorithms do, in fact, use division rather than
184 multiplication by reciprocals for better numerical accuracy; otherwise, they 254 multiplication by reciprocals for better numerical accuracy; otherwise, they
185 honor the above definition. Note that a diagonal matrix is never truncated due 255 honor the above definition. Note that a diagonal matrix is never truncated due
186 to ill-conditioning; otherwise, it would not be much useful for scaling. This 256 to ill-conditioning; otherwise, it would not be much useful for scaling. This
188 happens to be diagonal (an is thus not a special object) is of course treated 258 happens to be diagonal (an is thus not a special object) is of course treated
189 normally. 259 normally.
190 260
191 If @var{D1} and @var{D2} are both diagonal matrices, then the expressions 261 If @var{D1} and @var{D2} are both diagonal matrices, then the expressions
192 @example 262 @example
193 @var{D1} + @var{D2} 263 D1 + D2
194 @var{D1} - @var{D2} 264 D1 - D2
195 @var{D1} * @var{D2} 265 D1 * D2
196 @var{D1} / @var{D2} 266 D1 / D2
197 @var{D1} \ @var{D2} 267 D1 \ D2
198 @end example 268 @end example
199 again produce diagonal matrices, provided that normal 269 again produce diagonal matrices, provided that normal
200 dimension matching rules are obeyed. The relations used are same as described above. 270 dimension matching rules are obeyed. The relations used are same as described above.
201 271
202 Also, a diagonal matrix @var{D} can be multiplied or divided by a scalar, or raised 272 Also, a diagonal matrix @var{D} can be multiplied or divided by a scalar, or raised
203 to a scalar power if it is square, producing diagonal matrix result in all cases. 273 to a scalar power if it is square, producing diagonal matrix result in all cases.
204 274
205 A diagonal matrix can also be transposed or conjugate-transposed, giving the expected 275 A diagonal matrix can also be transposed or conjugate-transposed, giving the expected
206 result. Extracting a leading submatrix of a diagonal matrix, i.e. @code{@var{D}(1:m,1:n)}, 276 result. Extracting a leading submatrix of a diagonal matrix, i.e. @code{D(1:m,1:n)},
207 will produce a diagonal matrix, other indexing expressions will implicitly convert to 277 will produce a diagonal matrix, other indexing expressions will implicitly convert to
208 full matrix. 278 full matrix.
209 279
210 When involved in expressions with the element-by-element operators @code{.*}, 280 Adding a diagonal matrix to a full matrix only operates on the diagonal elements. Thus,
281 @example
282 A = A + eps * eye (n)
283 @end example
284 is an efficient method of augmenting the diagonal of a matrix. Subtraction works
285 analogically.
286
287 When involved in expressions with other element-by-element operators, @code{.*},
211 @code{./}, @code{.\} or @code{.^}, an implicit conversion to full matrix will 288 @code{./}, @code{.\} or @code{.^}, an implicit conversion to full matrix will
212 also take place. This is not always strictly necessary but chosen to facilitate 289 take place. This is not always strictly necessary but chosen to facilitate
213 better consistency with Matlab. 290 better consistency with Matlab.
214 291
215 @node Expressions Involving Permutation Matrices 292 @node Expressions Involving Permutation Matrices
216 @subsection Expressions Involving Permutation Matrices 293 @subsection Expressions Involving Permutation Matrices
217 294
293 @example 370 @example
294 [L, U, P] = lu (A); ## now L*U = P*A 371 [L, U, P] = lu (A); ## now L*U = P*A
295 x = U \ L \ P*b; 372 x = U \ L \ P*b;
296 @end example 373 @end example
297 374
375 @noindent
298 This is how you normalize columns of a matrix @var{X} to unit norm: 376 This is how you normalize columns of a matrix @var{X} to unit norm:
299 @example 377 @example
300 s = norm (X, "columns"); 378 s = norm (X, "columns");
301 X = diag (s) \ X; 379 X = diag (s) \ X;
302 @end example 380 @end example
303 381
382 @noindent
304 The following expression is a way to efficiently calculate the sign of a 383 The following expression is a way to efficiently calculate the sign of a
305 permutation, given by a permutation vector @var{p}. It will also work 384 permutation, given by a permutation vector @var{p}. It will also work
306 in earlier versions of Octave, but slowly. 385 in earlier versions of Octave, but slowly.
307 @example 386 @example
308 det (eye (length (p))(p, :)) 387 det (eye (length (p))(p, :))
309 @end example 388 @end example
310 389
390 @noindent
311 Finally, here's how you solve a linear system @code{A*x = b} 391 Finally, here's how you solve a linear system @code{A*x = b}
312 with Tikhonov regularization (ridge regression) using SVD (a skeleton only): 392 with Tikhonov regularization (ridge regression) using SVD (a skeleton only):
313 @example 393 @example
314 m = rows (A); n = columns (A); 394 m = rows (A); n = columns (A);
315 [U, S, V] = svd (A); 395 [U, S, V] = svd (A);
316 ## determine the regularization factor alpha 396 ## determine the regularization factor alpha
317 ## alpha = ... 397 ## alpha = ...
318 ## transform to orthogonal basis 398 ## transform to orthogonal basis
319 x = U'*b; 399 b = U'*b;
320 ## Use the standard formula, replacing A with S. We work with 400 ## Use the standard formula, replacing A with S.
321 ## diagonal matrices, so this will be plausibly efficient. 401 ## S is diagonal, so the following will be very fast and accurate.
322 x = (S'*S + alpha^2 * eye (n)) \ (S' * x); 402 x = (S'*S + alpha^2 * eye (n)) \ (S' * b);
323 ## transform to solution basis 403 ## transform to solution basis
324 x = V*x; 404 x = V*x;
325 @end example 405 @end example
406
326 407
327 @node Zeros Treatment 408 @node Zeros Treatment
328 @section The Differences in Treatment of Zero Elements 409 @section The Differences in Treatment of Zero Elements
329 410
330 Making diagonal and permutation matrices special matrix objects in their own 411 Making diagonal and permutation matrices special matrix objects in their own
345 An "assumed zero", on the contrary, is a zero matrix element implied by the 426 An "assumed zero", on the contrary, is a zero matrix element implied by the
346 matrix structure (diagonal, triangular) or a sparsity pattern; its value is 427 matrix structure (diagonal, triangular) or a sparsity pattern; its value is
347 usually not stored explicitly anywhere, but is implied by the underlying 428 usually not stored explicitly anywhere, but is implied by the underlying
348 data structure. 429 data structure.
349 430
350 The primary distinction is that an assumed zero, when multiplied or divided 431 The primary distinction is that an assumed zero, when multiplied
351 by any number, yields *always* a zero, even when, e.g., multiplied by @code{Inf} 432 by any number, or divided by any nonzero number,
433 yields *always* a zero, even when, e.g., multiplied by @code{Inf}
352 or divided by @code{NaN}. 434 or divided by @code{NaN}.
353 The reason for this behavior is that the numerical multiplication is not 435 The reason for this behavior is that the numerical multiplication is not
354 actually performed anywhere by the underlying algorithm; the result is 436 actually performed anywhere by the underlying algorithm; the result is
355 just assumed to be zero. Equivalently, one can say that the part of the 437 just assumed to be zero. Equivalently, one can say that the part of the
356 computation involving assumed zeros is performed symbolically, not numerically. 438 computation involving assumed zeros is performed symbolically, not numerically.