# HG changeset patch # User Jaroslav Hajek # Date 1235393744 -3600 # Node ID fcba62cc454939ea068fbfc8a581640d85bb66e3 # Parent dea5bd01e6d744179b7a5e9f38d44f09d5b7d336 add chapter about diagonal and permutation matrices to manual diff -r dea5bd01e6d7 -r fcba62cc4549 doc/ChangeLog --- a/doc/ChangeLog Mon Feb 23 12:47:30 2009 +0100 +++ b/doc/ChangeLog Mon Feb 23 13:55:44 2009 +0100 @@ -1,3 +1,11 @@ +2009-02-23 Jaroslav Hajek + + * interpreter/diagperm.txi: New file. + * interpreter/octave.texi: Reference it. + * intepreter/Makefile.in: Include it. + * interpreter/matrix.txi: Move @DOCSTRING(diag) to diagperm.txi. + + 2009-02-19 John W. Eaton * doc/interperter: Include @DOCSTRING commands for the following diff -r dea5bd01e6d7 -r fcba62cc4549 doc/interpreter/Makefile.in --- a/doc/interpreter/Makefile.in Mon Feb 23 12:47:30 2009 +0100 +++ b/doc/interpreter/Makefile.in Mon Feb 23 13:55:44 2009 +0100 @@ -110,9 +110,9 @@ SUB_SOURCE := arith.txi audio.txi basics.txi bugs.txi \ container.txi contrib.txi cp-idx.txi data.txi \ - debug.txi diffeq.txi dynamic.txi emacs.txi errors.txi eval.txi \ - expr.txi fn-idx.txi func.txi geometry.txi gpl.txi \ - grammar.txi image.txi install.txi interp.txi \ + debug.txi diffeq.txi diagperm.txi dynamic.txi emacs.txi \ + errors.txi eval.txi expr.txi fn-idx.txi func.txi geometry.txi \ + gpl.txi grammar.txi image.txi install.txi interp.txi \ intro.txi io.txi linalg.txi matrix.txi nonlin.txi numbers.txi \ oop.txi op-idx.txi optim.txi package.txi plot.txi poly.txi preface.txi \ quad.txi set.txi signal.txi sparse.txi stats.txi \ diff -r dea5bd01e6d7 -r fcba62cc4549 doc/interpreter/diagperm.txi --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/doc/interpreter/diagperm.txi Mon Feb 23 13:55:44 2009 +0100 @@ -0,0 +1,399 @@ +@c Copyright (C) 2009 Jaroslav Hajek +@c +@c This file is part of Octave. +@c +@c Octave is free software; you can redistribute it and/or modify it +@c under the terms of the GNU General Public License as published by the +@c Free Software Foundation; either version 3 of the License, or (at +@c your option) any later version. +@c +@c Octave is distributed in the hope that it will be useful, but WITHOUT +@c ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +@c FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +@c for more details. +@c +@c You should have received a copy of the GNU General Public License +@c along with Octave; see the file COPYING. If not, see +@c . + +@ifhtml +@set htmltex +@end ifhtml +@iftex +@set htmltex +@end iftex + +@node Diagonal and Permutation Matrices +@chapter Diagonal and Permutation Matrices + +@menu +* Basic Usage:: Creation and Manipulation of Diagonal and Permutation Matrices +* Matrix Algebra:: Linear Algebra with Diagonal and Permutation Matrices +* Function Support:: Functions That Are Aware of These Matrices +* Example Codes:: Some Examples of Usage +* Zeros Treatment:: The Differences in Treatment of Zero Elements +@end menu + +@node Basic Usage +@section The Creation and Manipulation of Diagonal and Permutation Matrices + +A diagonal matrix is defined as a matrix that has zero entries outside the main diagonal; +that is, @code{D(i,j) == 0} if @code{i != j}. +Most often, square diagonal matrices are considered; however, the definition can equally +be applied to nonsquare matrices, in which case we usually speak of a rectangular diagonal +matrix. For more information, see @url{http://en.wikipedia.org/wiki/Diagonal_matrix}. + +A permutation matrix is defined as a square matrix that has a single element equal to unity +in each row and each column; all other elements are zero. That is, there exists a +permutation (vector) @code{p} such that @code{P(i,j) == 1} if @code{j == p(i)} and +@code{P(i,j) == 0} otherwise. For more information, see +@url{http://en.wikipedia.org/wiki/Permutation_matrix}. + +Octave provides special treatment of real and complex rectangular diagonal matrices, +as well as permutation matrices. They are stored as special objects, using efficient +storage and algorithms, facilitating writing both readable and efficient matrix algebra +expressions in the Octave language. + +@menu +* Creating Diagonal Matrices:: +* Creating Permutation Matrices:: +* Explicit and Implicit Conversions:: +@end menu + +@node Creating Diagonal Matrices +@subsection Creating Diagonal Matrices + +The most common and easiest way to create a diagonal matrix is using the built-in +function @dfn{diag}. The expression @code{diag (@var{v})}, with @var{v} a vector, +will create a square diagonal matrix with elements on the main diagonal given +by the elements of @var{v}, and size equal to the length of @var{v}. +@code{diag (@var{v}, m, n)} can be used to construct a rectangular diagonal matrix. +The result of these expressions will be a special diagonal matrix object, rather +than a general matrix object. + +Some other built-in functions can also return diagonal matrices. Examples include +@dfn{balance} or @dfn{inv}. + +@DOCSTRING(diag) + +@node Creating Permutation Matrices +@subsection Creating Permutation Matrices + +For creating permutation matrices, Octave does not introduce a new function, but +rather overrides an existing syntax. That is, if @var{q} is a permutation vector +of length @var{n}, the expression +@example + @var{P} = eye (@var{n}) (:, @var{q}); +@end example +will create a permutation matrix - a special matrix object. +@example +eye (@var{n}) (@var{q}, :) +@end example +will also work (a row permutation matrix), as well as +@example +eye (@var{n}) (@var{q1}, @var{q2}). +@end example +Note that the identity, @code{eye (@var{n})}, is a diagonal matrix by definition, +but should work in any place where a permutation matrix is requested. + +Some other built-in functions can also return permutation matrices. Examples include +@dfn{inv} or @dfn{lu}. + +@node Explicit and Implicit Conversions +@subsection Explicit and Implicit Conversions + +The diagonal and permutation matrices are special objects in their own right. A number +of operations and built-in functions, are defined for these matrices to use special, +more efficient code than would be used for a full matrix in the same place. Examples +are given in further sections. + +To facilitate smooth mixing with full matrices, backward compatibility, and +compatibility with Matlab, the diagonal and permutation matrices should allow +any operation that works on full matrices, and will either treat is specially, +or implicitly convert themselves to full matrices. + +Instances include matrix indexing (except for extracting a single element), +indexed assignment, or applying most mapper functions, such as @dfn{exp}. + +An explicit conversion to a full matrix can be requested using the built-in +function @dfn{full}. It should also be noted that the diagonal and permutation +matrix objects will cache the result of the conversion after it is first +requested (explicitly or implicitly), so that subsequent conversions will +be very cheap. + +@node Matrix Algebra +@section Linear Algebra with Diagonal and Permutation Matrices + +As has been already said, diagonal and permutation matrices make it +possible to use efficient algorithms while preserving natural linear +algebra syntax. This section describes in detail the operations that +are treated specially when performed on these special matrix objects. + +@menu +* Expressions Involving Diagonal Matrices:: +* Expressions Involving Permutation Matrices:: +@end menu + +@node Expressions Involving Diagonal Matrices +@subsection Expressions Involving Diagonal Matrices + +Assume @var{D} is a diagonal matrix. If @var{M} is a full matrix, +then @code{@var{D}*@var{M}} will scale the rows of @var{M}. That means, +if @code{@var{S} = @var{D}*@var{M}}, then for each pair of indices +i,j it holds +@example +S(i,j) = D(i,i) * M(i,j). +@end example +Similarly, @code{M*D} will do a column scaling. + +The matrix @var{D} may also be rectangular, m-by-n where @code{m != n}. +If @code{m < n}, then the expression @code{D*M} is equivalent to +@example +D(:,1:m) * M(1:m,:), +@end example +i.e. trailing @code{n-m} rows of @var{M} are ignored. If @code{m > n}, +then @code{D*M} is equivalent to +@example +[D(1:n,n) * M; zeros(m-n, columns (M))], +@end example +i.e. null rows are appended to the result. +The situation for right-multiplication @code{M*D} is analogous. + +The expressions @code{D \ M} and @code{M / D} perform inverse scaling. +They are equivalent to solving a diagonal (or rectangular diagonal) +in a least-squares minimum-norm sense. In exact arithmetics, this is +equivalent to multiplying by a pseudoinverse. The pseudoinverse of +a rectangular diagonal matrix is again a rectangular diagonal matrix +of the same dimensions, where each nonzero diagonal element is replaced +by its reciprocal. +The matrix division algorithms do, in fact, use division rather than +multiplication by reciprocals for better numerical accuracy; otherwise, they +honor the above definition. Note that a diagonal matrix is never truncated due +to ill-conditioning; otherwise, it would not be much useful for scaling. This +is typically consistent with linear algebra needs. A full matrix that only +happens to be diagonal (an is thus not a special object) is of course treated +normally. + +If @var{D1} and @var{D2} are both diagonal matrices, then the expressions +@example +@var{D1} + @var{D2} +@var{D1} - @var{D2} +@var{D1} * @var{D2} +@var{D1} / @var{D2} +@var{D1} \ @var{D2} +@end example +again produce diagonal matrices, provided that normal +dimension matching rules are obeyed. The relations used are same as described above. + +Also, a diagonal matrix @var{D} can be multiplied or divided by a scalar, or raised +to a scalar power if it is square, producing diagonal matrix result in all cases. + +A diagonal matrix can also be transposed or conjugate-transposed, giving the expected +result. Extracting a leading submatrix of a diagonal matrix, i.e. @code{@var{D}(1:m,1:n)}, +will produce a diagonal matrix, other indexing expressions will implicitly convert to +full matrix. + +When involved in expressions with the element-by-element operators @code{.*}, +@code{./}, @code{.\} or @code{.^}, an implicit conversion to full matrix will +also take place. This is not always strictly necessary but chosen to facilitate +better consistency with Matlab. + +@node Expressions Involving Permutation Matrices +@subsection Expressions Involving Permutation Matrices + +If @var{P} is a permutation matrix and @var{M} a matrix, the expression +@code{P*M} will permute the rows of @var{M}. Similarly, @code{M*P} will +yield a column permutation. +Matrix division @code{P\M} and @code{M/P} can be used to do inverse permutation. + +The previously described syntax for creating permutation matrices can actually +help an user to understand the connection between a permutation matrix and +a permuting vector. Namely, the following holds, where @code{I = eye (n)} +is an identity matrix: +@example + I(p,:) * M = (I*M) (p,:) = M(p,:) +@end example +Similarly, +@example + M * I(:,p) = (M*I) (:,p) = M(:,p) +@end example + +The expressions @code{I(p,:)} and @code{I(:,p)} are permutation matrices. + +A permutation matrix can be transposed (or conjugate-transposed, which is the +same, because a permutation matrix is never complex), inverting the +permutation, or equivalently, turning a row-permutation matrix into a +column-permutation one. For permutation matrices, transpose is equivalent to +inversion, thus @code{P\M} is equivalent to @code{P'*M}. Transpose of a +permutation matrix (or inverse) is a constant-time operation, flipping only a +flag internally, and thus the choice between the two above equivalent +expressions for inverse permuting is completely up to the user's taste. + +Two permutation matrices can be multiplied or divided (if their sizes match), performing +a composition of permutations. Also a permutation matrix can be indexed by a permutation +vector (or two vectors), giving again a permutation matrix. +Any other operations do not generally yield a permutation matrix and will thus +trigger the implicit conversion. + +@node Function Support +@section Functions That Are Aware of These Matrices + +This section lists the built-in functions that are aware of diagonal and +permutation matrices on input, or can return them as output. Passed to other +functions, these matrices will in general trigger an implicit conversion. +(Of course, user-defined dynamically linked functions may also work with +diagonal or permutation matrices). + +@menu +* Diagonal Matrix Functions:: +* Permutation Matrix Functions:: +@end menu + +@node Diagonal Matrix Functions +@subsection Diagonal Matrix Functions + +@dfn{inv} and @dfn{pinv} can be applied to a diagonal matrix, yielding again +a diagonal matrix. @dfn{det} will use an efficient straightforward calculation +when given a diagonal matrix, as well as @dfn{cond}. +The following mapper functions can be applied to a diagonal matrix +without converting it to a full one: +@dfn{abs}, @dfn{real}, @dfn{imag}, @dfn{conj}, @dfn{sqrt}. +A diagonal matrix can also be returned from the @dfn{balance} +and @dfn{svd} functions. + +@node Permutation Matrix Functions +@subsection Permutation Matrix Functions + +@dfn{inv} and @dfn{pinv} will invert a permutation matrix, preserving its +specialness. @dfn{det} can be applied to a permutation matrix, efficiently +calculating the sign of the permutation (which is equal to the determinant). + +A permutation matrix can also be returned from the built-in functions +@dfn{lu} and @dfn{qr}, if a pivoted factorization is requested. + +@node Example Codes +@section Some Examples of Usage + +The following can be used to solve a linear system @code{A*x = b} +using the pivoted LU factorization: +@example + [L, U, P] = lu (A); ## now L*U = P*A + x = U \ L \ P*b; +@end example + +This is how you normalize columns of a matrix @var{X} to unit norm: +@example + s = norm (X, "columns"); + X = diag (s) \ X; +@end example + +The following expression is a way to efficiently calculate the sign of a +permutation, given by a permutation vector @var{p}. It will also work +in earlier versions of Octave, but slowly. +@example + det (eye (length (p))(p, :)) +@end example + +Finally, here's how you solve a linear system @code{A*x = b} +with Tikhonov regularization using SVD (a skeleton only): +@example + m = rows (A); n = columns (A); + [U, S, V] = svd (A); + ## determine the regularization factor alpha + ## ... + ## and solve + x = U'*b; + x = (S + alpha*eye (m, n)) \ x; ## this will be very efficient + x = V*x; +@end example + +@node Zeros Treatment +@section The Differences in Treatment of Zero Elements + +Making diagonal and permutation matrices special matrix objects in their own +right and the consequent usage of smarter algorithms for certain operations +implies, as a side effect, small differences in treating zeros. +The contents of this section applies also to sparse matrices, discussed in +the following chapter. + +The IEEE standard defines the result of the expressions @code{0*Inf} and +@code{0*NaN} as @code{NaN}, as it has been generally agreed that this is the +best compromise. +Numerical software dealing with structured and sparse matrices (including +Octave) however, almost always makes a distinction between a "numerical zero" +and an "assumed zero". +A "numerical zero" is a zero value occuring in a place where any floating-point +value could occur. It is normally stored somewhere in memory as an explicit +value. +An "assumed zero", on the contrary, is a zero matrix element implied by the +matrix structure (diagonal, triangular) or a sparsity pattern; its value is +usually not stored explicitly anywhere, but is implied by the underlying +data structure. + +The primary distinction is that an assumed zero, when multiplied or divided +by any number, yields *always* a zero, even when, e.g., multiplied by @code{Inf} +or divided by @code{NaN}. +The reason for this behavior is that the numerical multiplication is not +actually performed anywhere by the underlying algorithm; the result is +just assumed to be zero. Equivalently, one can say that the part of the +computation involving assumed zeros is performed symbolically, not numerically. + +This behavior not only facilitates the most straightforward and efficient +implementation of algorithms, but also preserves certain useful invariants, +like: +@itemize +@item scalar * diagonal matrix is a diagonal matrix +@item sparse matrix / scalar preserves the sparsity pattern +@item permutation matrix * matrix is equivalent to permuting rows +@end itemize +all of these natural mathematical truths would be invalidated by treating +assumed zeros as numerical ones. + +Note that certain competing software does not strictly follow this principle +and converts assumed zeros to numerical zeros in certain cases, while not doing +so in other cases. As of today, there are no intentions to mimick such behavior +in Octave. + +Examples of effects of assumed zeros vs. numerical zeros: +@example +Inf * eye (3) +@result{} + Inf 0 0 + 0 Inf 0 + 0 0 Inf + +Inf * speye (3) +@result{} +Compressed Column Sparse (rows = 3, cols = 3, nnz = 3 [33%]) + + (1, 1) -> Inf + (2, 2) -> Inf + (3, 3) -> Inf + +Inf * full (eye (3)) +@result{} + Inf NaN NaN + NaN Inf NaN + NaN NaN Inf + +@end example + +@example +diag(1:3) * [NaN; 1; 1] +@result{} + NaN + 2 + 3 + +sparse(1:3,1:3,1:3) * [NaN; 1; 1] +@result{} + NaN + 2 + 3 +[1,0,0;0,2,0;0,0,3] * [NaN; 1; 1] +@result{} + NaN + NaN + NaN +@end example + diff -r dea5bd01e6d7 -r fcba62cc4549 doc/interpreter/matrix.txi --- a/doc/interpreter/matrix.txi Mon Feb 23 12:47:30 2009 +0100 +++ b/doc/interpreter/matrix.txi Mon Feb 23 13:55:44 2009 +0100 @@ -239,8 +239,6 @@ @DOCSTRING(randperm) -@DOCSTRING(diag) - The functions @code{linspace} and @code{logspace} make it very easy to create vectors with evenly or logarithmically spaced elements. @xref{Ranges}. diff -r dea5bd01e6d7 -r fcba62cc4549 doc/interpreter/octave.texi --- a/doc/interpreter/octave.texi Mon Feb 23 12:47:30 2009 +0100 +++ b/doc/interpreter/octave.texi Mon Feb 23 13:55:44 2009 +0100 @@ -165,6 +165,7 @@ * Arithmetic:: * Linear Algebra:: * Nonlinear Equations:: +* Diagonal and Permutation Matrices:: * Sparse Matrices:: * Numerical Integration:: * Differential Equations:: @@ -419,6 +420,14 @@ * Matrix Factorizations:: * Functions of a Matrix:: +Diagonal and Permutation Matrices + +* Basic Usage:: +* Matrix Algebra:: +* Function Support:: +* Example Codes:: +* Zeros Treatment:: + Sparse Matrices * Basics:: @@ -594,6 +603,7 @@ @include arith.texi @include linalg.texi @include nonlin.texi +@include diagperm.texi @include sparse.texi @include quad.texi @include diffeq.texi