Mathematics

Permutation and Reordering

A permutation of the rows and columns of a sparse matrix S can be represented in two ways:

• A permutation matrix P acts on the rows of S as P*S or on the columns as S*P'.
• A permutation vector p, which is a full vector containing a permutation of 1:n, acts on the rows of S as S(p,:), or on the columns as S(:,p).

For example, the statements

• p = [1 3 4 2 5]
  I = eye(5,5);
  P = I(p,:);
  e = ones(4,1);
  S = diag(11:11:55) + diag(e,1) + diag(e,-1)
  

produce

• p =
  
       1     3     4     2     5
  
  P =
  
       1     0     0     0     0
       0     0     1     0     0
       0     0     0     1     0
       0     1     0     0     0
       0     0     0     0     1
  
  S =
  
      11     1     0     0     0
       1    22     1     0     0
       0     1    33     1     0
       0     0     1    44     1
       0     0     0     1    55
  

You can now try some permutations using the permutation vector p and the permutation matrix P. For example, the statements S(p,:) and P*S produce

• ans =
  
      11     1     0     0     0
       0     1    33     1     0
       0     0     1    44     1
       1    22     1     0     0
       0     0     0     1    55
  

Similarly, S(:,p) and S*P' produce

• ans =
  
      11     0     0     1     0
       1     1     0    22     0
       0    33     1     1     0
       0     1    44     0     1
       0     0     1     0    55
  

If P is a sparse matrix, then both representations use storage proportional to n and you can apply either to S in time proportional to nnz(S). The vector representation is slightly more compact and efficient, so the various sparse matrix permutation routines all return full row vectors with the exception of the pivoting permutation in LU (triangular) factorization, which returns a matrix compatible with earlier versions of MATLAB.

To convert between the two representations, let I = speye(n) be an identity matrix of the appropriate size. Then,

• P = I(p,:)
  P' = I(:,p)
  p = (1:n)*P'
  p = (P*(1:n)')'
  

The inverse of P is simply R = P'. You can compute the inverse of p with r(p) = 1:n.

• r(p) = 1:5
  
  r =
  
       1     4     2     3     5
  

Reordering for Sparsity

Reordering the columns of a matrix can often make its LU or QR factors sparser. Reordering the rows and columns can often make its Cholesky, factors sparser. The simplest such reordering is to sort the columns by nonzero count. This is sometimes a good reordering for matrices with very irregular structures, especially if there is great variation in the nonzero counts of rows or columns.

The function p = colperm(S) computes this column-count permutation. The colperm M-file has only a single line.

• [ignore,p] = sort(full(sum(spones(S))));
  

This line performs these steps:

1. The inner call to spones creates a sparse matrix with ones at the location of every nonzero element in S.
2. The sum function sums down the columns of the matrix, producing a vector that contains the count of nonzeros in each column.
3. full converts this vector to full storage format.
4. sort sorts the values in ascending order. The second output argument from sort is the permutation that sorts this vector.

Reordering to Reduce Bandwidth

The reverse Cuthill-McKee ordering is intended to reduce the profile or bandwidth of the matrix. It is not guaranteed to find the smallest possible bandwidth, but it usually does. The function symrcm(A) actually operates on the nonzero structure of the symmetric matrix A + A', but the result is also useful for asymmetric matrices. This ordering is useful for matrices that come from one-dimensional problems or problems that are in some sense "long and thin."

Minimum Degree Ordering

The degree of a node in a graph is the number of connections to that node. This is the same as the number of off-diagonal nonzero elements in the corresponding row of the adjacency matrix. The minimum degree algorithm generates an ordering based on how these degrees are altered during Gaussian elimination or Cholesky factorization. It is a complicated and powerful algorithm that usually leads to sparser factors than most other orderings, including column count and reverse Cuthill-McKee.

MATLAB functions implement two methods for each of two types of matrices: symamd and symmmd for symmetric matrices, and colamd and colmmd for nonsymmetric matrices. colamd and colmmd also work for symmetric matrices of the form A*A' or A'*A.

Because the most time-consuming part of a minimum degree ordering algorithm is keeping track of the degree of each node, all four functions use an approximation to the degree, rather the exact degree. As a result:

• Factorizations obtained using colmmd and symmmd tend to have more nonzero elements than if the implementation used exact degrees.
• colamd and symamd use a tighter approximation than colmmd and symmmd. They generate orderings that are as good as could be obtained using exact degrees.
• colamd and symamd are faster than colmmd and symmmd, respectively. This is true particularly for very large matrices.

You can change various parameters associated with details of the algorithms using the spparms function.

For details on the algorithms used by colmmd and symmmd, see [4]. For details on the algorithms used by colamd and symamd, see [5]. The approximate degree used in colamd and symamd is based on [1].

Standard Mathematical Operations

Factorization