Mathematics |
Permutation and Reordering
A permutation of the rows and columns of a sparse matrix S
can be represented in two ways:
P
acts on the rows of S
as P*S
or on the columns as S*P'
.
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)
.
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)
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
Similarly, S(:,p)
and S*P'
produce
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,
The inverse of P
is simply R = P'
. You can compute the inverse of p
with r(p) = 1:n
.
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.
This line performs these steps:
spones
creates a sparse matrix with ones at the location of every nonzero element in S
.
sum
function
sums down the columns of the matrix, producing a vector that contains the count of nonzeros in each column.
full
converts this vector to full storage format.
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."
Approximate 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 approximate 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. Because the keeping track of the degree of each node is very time-consuming, the approximate minimum degree algorithm uses an approximation to the degree, rather than the exact degree.
The following MATLAB functions implement the approximate minimum degree algorithm:
symamd
-- Use with symmetric matrices
colamd
-- Use with nonsymmetric matrices and symmetric matrices of the form A*A'
or A'*A
.
See Reordering and Factorization for an example using symamd
.
You can change various parameters associated with details of the algorithms using the spparms
function.
For details on the algorithms used by colamd
and symamd
, see [5]. The approximate degree the algorithms use is based on [1].
Standard Mathematical Operations | Factorization |
© 1994-2005 The MathWorks, Inc.