This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English version of the page.

Note: This page has been translated by MathWorks. Click here to see
To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

Sparse Matrix Reordering

This example shows how reordering the rows and columns of a sparse matrix can influence the speed and storage requirements of a matrix operation.

Visualizing a Sparse Matrix

A spy plot shows the nonzero elements in a matrix.

This spy plot shows a sparse symmetric positive definite matrix derived from a portion of the Harwell-Boeing test matrix west0479. This matrix describes connections in a model of a diffraction column in a chemical plant.

load west0479.mat
A = west0479;
S = A * A' + speye(size(A));
pct = 100 / numel(A);

title('A Sparse Symmetric Matrix')
nz = nnz(S);
xlabel(sprintf('Nonzeros = %d (%.3f%%)',nz,nz*pct));

Computing the Cholesky Factor

Compute the Cholesky factor L, where S = L*L'. Notice that L contains many more nonzero elements than the unfactored S, because the computation of the Cholesky factorization creates fill-in nonzeros. These fill-in values slow down the algorithm and increase storage cost.

L = chol(S,'lower');
title('Cholesky Decomposition of S')
nc(1) = nnz(L);
xlabel(sprintf('Nonzeros = %d (%.2f%%)',nc(1),nc(1)*pct));

Reordering to Speed Up Calculation

By reordering the rows and columns of a matrix, it is possible to reduce the amount of fill-in that factorization creates, thereby reducing time and storage cost.

Three different reorderings supported by MATLAB® are:

  • Reverse Cuthill-McKee

  • Column count

  • Minimum degree

Test the effects of these sparse matrix reorderings on the west0479 matrix.

Reordering 1: Reverse Cuthill-McKee

The symrcm command uses the reverse Cuthill-McKee reordering algorithm to move all nonzero elements closer to the diagonal, reducing the bandwidth of the original matrix.

p = symrcm(S);
title('S(p,p) After Cuthill-McKee Ordering')
nz = nnz(S);
xlabel(sprintf('Nonzeros = %d (%.3f%%)',nz,nz*pct));

The fill-in produced by Cholesky factorization is confined to the band, so factorizing the reordered matrix takes less time and less storage.

L = chol(S(p,p),'lower');
title('chol(S(p,p)) After Cuthill-McKee Ordering')
nc(2) = nnz(L);
xlabel(sprintf('Nonzeros = %d (%.2f%%)', nc(2),nc(2)*pct));

Reordering 2: Column Count

The colperm command uses the column count reordering algorithm to move rows and columns with higher nonzero count towards the end of the matrix.

q = colperm(S);
title('S(q,q) After Column Count Ordering')
nz = nnz(S);
xlabel(sprintf('Nonzeros = %d (%.3f%%)',nz,nz*pct));

For this matrix, the column count ordering happens to reduce the time and storage for Cholesky factorization, but this behavior is not guaranteed in general.

L = chol(S(q,q),'lower');
title('chol(S(q,q)) After Column Count Ordering')
nc(3) = nnz(L);
xlabel(sprintf('Nonzeros = %d (%.2f%%)',nc(3),nc(3)*pct));

Reordering 3: Minimum Degree

The symamd command uses the approximate minimum degree algorithm (a powerful graph-theoretic technique) to produce large blocks of zeros in the matrix.

r = symamd(S);
title('S(r,r) After Minimum Degree Ordering')
nz = nnz(S);
xlabel(sprintf('Nonzeros = %d (%.3f%%)',nz,nz*pct));

The Cholesky factorization preserves the blocks of zeros produced by the minimum degree algorithm. This structure can significantly reduce time and storage costs.

L = chol(S(r,r),'lower');
title('chol(S(r,r)) After Minimum Degree Ordering')
nc(4) = nnz(L);
xlabel(sprintf('Nonzeros = %d (%.2f%%)',nc(4),nc(4)*pct));

Summarizing Results

This bar chart summarizes the effects of reordering the matrix before performing the Cholesky factorization. While the Cholesky factorization of the original matrix had about 13% of its elements as nonzeros, using symamd reduces that density to only about 4%.

labels = {'Original','Cuthill-McKee','Column Count','Min Degree'};
title('Nonzeros After Cholesky Factorization')
ax = gca;
ax.XTickLabel = labels;
ax.XTickLabelRotation = -45;

See Also

| | | | |