Diagonal and Permutation Matrices
*********************************

Creating and Manipulating Diagonal and Permutation Matrices
===========================================================

A diagonal matrix is defined as a matrix that has zero entries outside the main diagonal; that is, D(i,j) == 0 if i != j. Most often, square diagonal matrices are considered; however, the definition can equally be applied to non-square matrices, in which case we usually speak of a rectangular 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) p such that P(i,j) == 1 if j == p(i) and P(i,j) == 0 otherwise.

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.

Creating Diagonal Matrices
--------------------------

The most common and easiest way to create a diagonal matrix is using the built-in function diag.  The expression diag (v), with V a vector, will create a square diagonal matrix with elements on the main diagonal given by the elements of V, and size equal to the length of V. diag (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.

Diagonal matrix with unit elements can be created using eye. Some other built-in functions can also return diagonal matrices.  Examples include balance or inv.

Example:

       diag (1:4)
     =>
     Diagonal Matrix

        1   0   0   0
        0   2   0   0
        0   0   3   0
        0   0   0   4

       diag (1:3,5,3)

     =>
     Diagonal Matrix

        1   0   0
        0   2   0
        0   0   3
        0   0   0
        0   0   0

Creating Permutation Matrices
-----------------------------

For creating permutation matrices, Octave does not introduce a new function, but rather overrides an existing syntax: permutation matrices can be conveniently created by indexing an identity matrix by permutation vectors. That is, if Q is a permutation vector of length N, the expression

       P = eye (n) (:, q);

will create a permutation matrix - a special matrix object.

     eye (n) (q, :)

will also work (and create a row permutation matrix), as well as

     eye (n) (q1, q2).

For example:

       eye (4) ([1,3,2,4],:)
     =>
     Permutation Matrix

        1   0   0   0
        0   0   1   0
        0   1   0   0
        0   0   0   1

       eye (4) (:,[1,3,2,4])
     =>
     Permutation Matrix

        1   0   0   0
        0   0   1   0
        0   1   0   0
        0   0   0   1

Mathematically, an identity matrix is both diagonal and permutation matrix. In Octave, eye (n) returns a diagonal matrix, because a matrix can only have one class.  You can convert this diagonal matrix to a permutation matrix by indexing it by an identity permutation, as shown below. This is a special property of the identity matrix; indexing other diagonal matrices generally produces a full matrix.

       eye (3)
     =>
     Diagonal Matrix

        1   0   0
        0   1   0
        0   0   1

       eye(3)(1:3,:)
     =>
     Permutation Matrix

        1   0   0
        0   1   0
        0   0   1

Some other built-in functions can also return permutation matrices.  Examples include inv or lu.

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 it specially, or implicitly convert themselves to full matrices.

Instances include matrix indexing, except for extracting a single element or a leading submatrix, indexed assignment, or applying most mapper functions, such as exp.

An explicit conversion to a full matrix can be requested using the built-in function 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.

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.

Expressions Involving Diagonal Matrices
---------------------------------------

Assume D is a diagonal matrix.  If M is a full matrix, then D*M will scale the rows of M.  That means, if S = D*M, then for each pair of indices i,j it holds

     S(i,j) = D(i,i) * M(i,j).

Similarly, M*D will do a column scaling.

The matrix D may also be rectangular, m-by-n where m != n. If m < n, then the expression D*M is equivalent to

     D(:,1:m) * M(1:m,:),

i.e., trailing n-m rows of M are ignored.  If m > n, then D*M is equivalent to

     [D(1:n,n) * M; zeros(m-n, columns (M))],

i.e., null rows are appended to the result. The situation for right-multiplication M*D is analogous.

The expressions D \ M and M / D perform inverse scaling. They are equivalent to solving a diagonal (or rectangular diagonal) in a least-squares minimum-norm sense.  In exact arithmetic, this is equivalent to multiplying by a pseudoinverse.  The pseudoinverse of a rectangular diagonal matrix is again a rectangular diagonal matrix with swapped 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.

Multiplication and division by diagonal matrices works efficiently also when combined with sparse matrices, i.e., D*S, where D is a diagonal matrix and S is a sparse matrix scales the rows of the sparse matrix and returns a sparse matrix.  The expressions S*D, D\S, S/D work analogically.

If D1 and D2 are both diagonal matrices, then the expressions

     D1 + D2
     D1 - D2
     D1 * D2
     D1 / D2
     D1 \ D2

again produce diagonal matrices, provided that normal dimension matching rules are obeyed.  The relations used are same as described above.

Also, a diagonal matrix 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., D(1:m,1:n), will produce a diagonal matrix, other indexing expressions will implicitly convert to full matrix.

Adding a diagonal matrix to a full matrix only operates on the diagonal elements.  Thus,

     A = A + eps * eye (n)

is an efficient method of augmenting the diagonal of a matrix.  Subtraction works analogically.

When involved in expressions with other element-by-element operators, .*, ./, .\ or .^, an implicit conversion to full matrix will take place.  This is not always strictly necessary but chosen to facilitate better consistency with MATLAB.

Expressions Involving Permutation Matrices
------------------------------------------

If P is a permutation matrix and M a matrix, the expression P*M will permute the rows of M.  Similarly, M*P will yield a column permutation. Matrix division P\M and 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 I = eye (n) is an identity matrix:

       I(p,:) * M = (I*M) (p,:) = M(p,:)

Similarly,

       M * I(:,p) = (M*I) (:,p) = M(:,p)

The expressions I(p,:) and 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 P\M is equivalent to 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.

Multiplication and division by permutation matrices works efficiently also when combined with sparse matrices, i.e., P*S, where P is a permutation matrix and S is a sparse matrix permutes the rows of the sparse matrix and returns a sparse matrix.  The expressions S*P, P\S, S/P work analogically.

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.

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

Diagonal Matrix Functions
-------------------------

inv and pinv can be applied to a diagonal matrix, yielding again a diagonal matrix.  det will use an efficient straightforward calculation when given a diagonal matrix, as well as cond. The following mapper functions can be applied to a diagonal matrix without converting it to a full one: abs, real, imag, conj, sqrt. A diagonal matrix can also be returned from the balance and svd functions. The sparse function will convert a diagonal matrix efficiently to a sparse matrix.

Permutation Matrix Functions
----------------------------

inv and pinv will invert a permutation matrix, preserving its specialness.  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 lu and qr, if a pivoted factorization is requested.

The sparse function will convert a permutation matrix efficiently to a sparse matrix. The find function will also work efficiently with a permutation matrix, making it possible to conveniently obtain the permutation indices.

Some Examples of Usage
======================

The following can be used to solve a linear system A*x = b using the pivoted LU factorization:

       [L, U, P] = lu (A); ## now L*U = P*A
       x = U \ L \ P*b;

This is one way to normalize columns of a matrix X to unit norm:

       s = norm (X, "columns");
       X /= diag (s);

The same can also be accomplished with broadcasting (see Broadcasting):

       s = norm (X, "columns");
       X ./= s;

The following expression is a way to efficiently calculate the sign of a permutation, given by a permutation vector P.  It will also work in earlier versions of Octave, but slowly.

       det (eye (length (p))(p, :))

Finally, here's how you solve a linear system A*x = b with Tikhonov regularization (ridge regression) using SVD (a skeleton only):

       m = rows (A); n = columns (A);
       [U, S, V] = svd (A);
       ## determine the regularization factor alpha
       ## alpha = ...
       ## transform to orthogonal basis
       b = U'*b;
       ## Use the standard formula, replacing A with S.
       ## S is diagonal, so the following will be very fast and accurate.
       x = (S'*S + alpha^2 * eye (n)) \ (S' * b);
       ## transform to solution basis
       x = V*x;

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 0*Inf and 0*NaN as 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 occurring 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 by any number, or divided by any nonzero number, yields *always* a zero, even when, e.g., multiplied by Inf or divided by 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:

  * scalar * diagonal matrix is a diagonal matrix

  * sparse matrix / scalar preserves the sparsity pattern

  * permutation matrix * matrix is equivalent to permuting rows

all of these natural mathematical truths would be invalidated by treating assumed zeros as numerical ones.

Note that MATLAB 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 mimic such behavior in Octave.

Examples of effects of assumed zeros vs. numerical zeros:

     Inf * eye (3)
     =>
        Inf     0     0
          0   Inf     0
          0     0   Inf

     Inf * speye (3)
     =>
     Compressed Column Sparse (rows = 3, cols = 3, nnz = 3 [33%])

       (1, 1) -> Inf
       (2, 2) -> Inf
       (3, 3) -> Inf

     Inf * full (eye (3))
     =>
        Inf   NaN   NaN
        NaN   Inf   NaN
        NaN   NaN   Inf

     diag (1:3) * [NaN; 1; 1]
     =>
        NaN
          2
          3

     sparse (1:3,1:3,1:3) * [NaN; 1; 1]
     =>
        NaN
          2
          3
     [1,0,0;0,2,0;0,0,3] * [NaN; 1; 1]
     =>
        NaN
        NaN
        NaN