Month: April 2022

What Is a Permutation Matrix?

A permutation matrix is a square matrix in which every row and every column contains a single 1 and all the other elements are zero. Such a matrix, P say, is orthogonal, that is, P^TP = PP^T = I_n, so it is nonsingular and has determinant \pm 1. The total number of n\times n permutation matrices is n!.

Premultiplying a matrix by P reorders the rows and postmultiplying by P reorders the columns. A permutation matrix P that has the desired reordering effect is constructed by doing the same operations on the identity matrix.

Examples of permutation matrices are the identity matrix I_n, the reverse identity matrix J_n, and the shift matrix K_n (also called the cyclic permutation matrix), illustrated for n = 4 by

\notag I_4 = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}, \qquad J_4 = \begin{bmatrix} 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \end{bmatrix}, \qquad K_4 = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \end{bmatrix}.

Pre- or postmultiplying a matrix by J_n reverses the order of the rows and columns, respectively. Pre- or postmultiplying a matrix by K_n shifts the rows or columns, respectively, one place forward and moves the first one to the last position—that is, it cyclically permutes the rows or columns. Note that J_n is a symmetric Hankel matrix and K_n is a circulant matrix.

An elementary permutation matrix P differs from I_n in just two rows and columns, i and j, say. It can be written P = I_n - (e_i-e_j)(e_i-e_j)^T, where e_i is the ith column of I_n. Such a matrix is symmetric and so satisfies P^2 = I_n, and it has determinant -1. A general permutation matrix can be written as a product of elementary permutation matrices P = P_1P_2\dots P_k, where k is such that \det(P) = (-1)^k.

It is easy to show that \det(\lambda I - K_n) = \lambda^n - 1, which means that the eigenvalues of K_n are 1, w, w^2, \dots, w^{n-1}, where w = \exp(2\pi\mathrm{i}/n) is the nth root of unity. The matrix K_n has two diagonals of 1s, which move up through the matrix as it is powered: K_n^i \ne I for i< n and K_n^n = I. The following animated gif superposes MATLAB spy plots of K_{64}, K_{64}^2, …, K_{64}^{64} = I_{64}.

shift_powers.gif

The shift matrix K_n plays a fundamental role in characterizing irreducible permutation matrices. Recall that a matrix A\in\mathbb{C}^{n\times n} is irreducible if there does not exist a permutation matrix P such that

\notag P^TAP = \begin{bmatrix} A_{11} & A_{12} \\ 0 & A_{22} \end{bmatrix},

where A_{11} and A_{22} are square, nonempty submatrices.

Theorem 1. For a permutation matrix P \in \mathbb{R}^{n \times n} the following conditions are equivalent.

  • P is irreducible.
  • There exists a permutation matrix Q such that Q^{-1}PQ = K_n
  • The eigenvalues of P are 1, w, w^2, \dots, w^{n-1}.

One consequence of Theorem 1 is that for any irreducible permutation matrix P, \mathrm{rank}(P - I) = \mathrm{rank}(K_n - I) = n-1.

The next result shows that a reducible permutation matrix can be expressed in terms of irreducible permutation matrices.

Theorem 2. Every reducible permutation matrix is permutation similar to a direct sum of irreducible permutation matrices.

Another notable permutation matrix is the vec-permutation matrix, which relates A\otimes B to B\otimes A, where \otimes is the Kronecker product.

A permutation matrix is an example of a doubly stochastic matrix: a nonnegative matrix whose row and column sums are all equal to 1. A classic result characterizes doubly stochastic matrices in terms of permutation matrices.

Theorem 3 (Birkhoff). A matrix is doubly stochastic if and only if it is a convex combination of permutation matrices.

In coding, memory can be saved by representing a permutation matrix P as an integer vector p, where p_i is the column index of the 1 within the ith row of P. MATLAB functions that return permutation matrices can also return the permutation in vector form. Here is an example with the MATLAB lu function that illustrates how permuting a matrix can be done using the vector permutation representation.

>> A = gallery('frank',4), [L,U,P] = lu(A); P
A =
     4     3     2     1
     3     3     2     1
     0     2     2     1
     0     0     1     1
P =
     1     0     0     0
     0     0     1     0
     0     0     0     1
     0     1     0     0
>> P*A
ans =
     4     3     2     1
     0     2     2     1
     0     0     1     1
     3     3     2     1
>> [L,U,p] = lu(A,'vector'); p
p =
     1     3     4     2
>> A(p,:)
ans =
     4     3     2     1
     0     2     2     1
     0     0     1     1
     3     3     2     1

For more on handling permutations in MATLAB see section 24.3 of MATLAB Guide.

Notes

For proofs of Theorems 1–3 see Zhang (2011, Sec. 5.6). Theorem 3 is also proved in Horn and Johnson (2013, Thm. 8.7.2).

Permutations play a key role in the fast Fourier transform and its efficient implementation; see Van Loan (1992).

References

Related Blog Posts

This article is part of the “What Is” series, available from https://nhigham.com/index-of-what-is-articles/ and in PDF form from the GitHub repository https://github.com/higham/what-is.

What’s New in MATLAB R2022a?

In this post I discuss some of the new features in MATLAB R2022a, focusing on ones that relate to my particular interests. See the release notes for a detailed list of the many changes in MATLAB and its toolboxes. For my articles about new features in earlier releases, see here.

Themes

MATLAB Online now has themes, including a dark theme (which is my preference). We will have to wait for a future release for themes to be supported on desktop MATLAB.

Economy Factorizations

One can now write qr(A,'econ') instead of qr(A,0) and gsvd(A,B,'econ') instead of gsvd(A,B) for the “economy size” decompositions. This is useful as the 'econ' form is more descriptive. The svd function already supported the 'econ' argument. The economy-size QR factorization is sometimes called the thin QR factorization.

Tie Breaking in the round Function

The round function, which rounds to the nearest integer, now breaks ties by rounding away from zero by default and has several other tie-breaking options (albeit not stochastic rounding). See a sequence of four blog posts on this topic by Cleve Moler starting with this one from February 2021.

Tolerances for null and orth

The null (nullspace) and orth (orthonormal basis for the range) functions now accept a tolerance as a second argument, and any singular values less than that tolerance are treated as zero. The default tolerance is max(size(A)) * eps(norm(A)). This change brings the two functions into line with rank, which already accepted the tolerance. If you are working in double precision (the MATLAB default) and your matrix has inherent errors of order 10^{-8} (for example), you might set the tolerance to 10^{-8}, since singular values smaller than this are indistinguishable from zero.

Unit Testing Reports

The unit testing framework can now generate docx, html, and pdf reports after test execution, by using the function generatePDFReport in the latter case. This is useful for keeping a record of test results and for printing them. We use unit testing in Anymatrix and have now added an option to return the results in a variable so that the user can call one of these new functions.

Checking Arrays for Special Values

Previously, if you wanted to check whether a matrix had all finite values you would need to use a construction such as all(all(isfinite(A))) or all(isfinite(A),'all'). The new allfinite function does this in one go: allfinite(A) returns true or false according as all the elements of A are finite or not, and it works for arrays of any dimension.

Similarly, anynan and anymissing check for NaNs or missing values. A missing value is a NaN for numerical arrays, but is indicated in other ways for other data types.

Linear Algebra on Multidimensional Arrays

The new pagemldivide, pagemrdivide, and pageinv functions solve linear equations and calculate matrix inverses using pages of d-dimensional arrays, while tensorprod calculates tensor products (inner products, outer products, or a combination of the two) between two d-dimensional arrays.

Animated GIFs

The append option of the exportgraphics function now supports the GIF format, enabling one to create animated GIFs (previously only multipage PDF files were supported). The key command is exportgraphics(gca,file_name,"Append",true). There are other ways of creating animated GIFs in MATLAB, but this one is particularly easy. Here is an example M-file (based on cheb3plot in MATLAB Guide) with its output below.

%CHEB_GIF  Animated GIF of Chebyshev polynomials.
%   Based on cheb3plot in MATLAB Guide.
x = linspace(-1,1,1500)';
p = 49
Y = ones(length(x),p);

Y(:,2) = x;
for k = 3:p
  Y(:,k) = 2*x.*Y(:,k-1) - Y(:,k-2);
end

delete cheby_animated.gif
a = get(groot,'defaultAxesColorOrder'); m = length(a);

for j = 1:p-1 % length(k)
    plot(x,Y(:,j),'LineWidth',1.5,'color',a(1+mod(j-1,m),:));
    xlim([-1 1]), ylim([-1 1])  % Must freeze axes.
    title(sprintf('%2.0f', j),'FontWeight','normal')
    exportgraphics(gca,"cheby_animated.gif","Append",true)
end

cheby_animated.gif