Matrix Ordering and LU Decomposition

In this demonstration, we see that the ordering of a sparse matrix can make a big difference in the sparsity of the LU decomposition of the matrix. We first define a sparse matrix:

ArrowMatrix = SparseArray[{Band[{1,1}]-> 1,Band[{1,10},Automatic,{1,0}]->0.1,Band[{10,1},Automatic,{0,1}]->0.1},{10,10}]

ArrowMatrix // MatrixForm

MatrixPlot[ArrowMatrix]

This code recovers the explicit LU decomposition from the combined form provided by LUDecomposition.

ExplicitLUDecomposition[M_] := Module[{lu,p,c,L,U,i,j,n},
n = Dimensions[M][[1]];
{lu,p,c} = LUDecomposition[M];
L=lu SparseArray[{i_,j_}/;j<i->1,{n,n}]+IdentityMatrix[n];
U=lu SparseArray[{i_,j_}/;j≥i->1,{n,n}];
{L,U,p,c}
]

{L,U,p,c} = ExplicitLUDecomposition[ArrowMatrix] ;

And here’s the LU decomposition in graphical form (so that we can see how sparse the L and U factors are):

GraphicsGrid[{{MatrixPlot[L],MatrixPlot[U],MatrixPlot[ArrowMatrix]}}]

We now try reordering the ArrowMatrix by applying permutations to see what effect that has on the LU decomposition.

PermuteMatrix[RowPerm_,ColPerm_,M_] := 
Transpose[Permute[Transpose[Permute[M,RowPerm]],ColPerm]];

ReorderedArrowMatrix = PermuteMatrix[Cycles[{{1,10}}],Cycles[{{1,10}}],ArrowMatrix]

MatrixPlot[ReorderedArrowMatrix]

{L,U,p,c} = ExplicitLUDecomposition[ReorderedArrowMatrix] ;
GraphicsGrid[{{MatrixPlot[L],MatrixPlot[U],MatrixPlot[ReorderedArrowMatrix]}}]

L // MatrixForm

Grim! If we measure the sparsity of the L and U entries of a matrix by counting nonzero entries in the L and U factors and dividing by n^2 (the maximum number of nonzero entries), as follows:

LUSparsity[M_] := Module[{lu,p,c},
{lu,p,c} = LUDecomposition[M];
Length[Select[Flatten[lu],Abs[#] > 0 &]]/(Dimensions[M][[1]]^2)
];

We can then measure the sparsity of these two decompositions:

LUSparsity[ArrowMatrix]

LUSparsity[ReorderedArrowMatrix]

What happens if we take a random permutation of the matrix?

RandomArrow := PermuteMatrix[RandomPermutation[10],RandomPermutation[10],ArrowMatrix];

R = RandomArrow

MatrixPlot[R]

{L,U,p,c} = ExplicitLUDecomposition[R] ;
GraphicsGrid[{{MatrixPlot[L],MatrixPlot[U],MatrixPlot[R]}}]

LUSparsity[R]

We now collect some data on the sparsity of LU decompositions of random permutations of our matrix.

SparsityData = Table[LUSparsity[RandomArrow],{i,1,1000}];

Histogram[SparsityData]

Algorithms for optimizing the sparsity of the LU decomposition.

Mathematica has some demonstration algorithms for reducing the fill in matrices using various heuristics. Two of the variants in most common use are the “Reverse Cuthill-McKee Algorithm”, and “Sloan’s algorithm”. This shows the result of each of these on our demonstration matrix.

Needs["GraphUtilities`"];

{r,c} = MinimumBandwidthOrdering[ArrowMatrix,Method->"Sloan"]
MBOSloan = PermuteMatrix[FindPermutation[Table[i,{i,1,10}],r],FindPermutation[Table[i,{i,1,10}],c],ArrowMatrix]

{r,c} = MinimumBandwidthOrdering[ArrowMatrix,Method->"RCMD"]
MBORCMD = PermuteMatrix[FindPermutation[Table[i,{i,1,10}],r],FindPermutation[Table[i,{i,1,10}],c],ArrowMatrix]

GraphicsGrid[{{MatrixPlot[ArrowMatrix],MatrixPlot[MBOSloan],MatrixPlot[MBORCMD]}}]

{L,U,p,c} = ExplicitLUDecomposition[ArrowMatrix] ;
GraphicsGrid[{{MatrixPlot[L],MatrixPlot[U],MatrixPlot[ArrowMatrix]}}]
LUSparsity[ArrowMatrix]

{L,U,p,c} = ExplicitLUDecomposition[MBO] ;
GraphicsGrid[{{MatrixPlot[L],MatrixPlot[U],MatrixPlot[MBO]}}]
LUSparsity[MBOSloan]

{L,U,p,c} = ExplicitLUDecomposition[MBORCMD] ;
GraphicsGrid[{{MatrixPlot[L],MatrixPlot[U],MatrixPlot[MBORCMD]}}]
LUSparsity[MBORCMD]

A more realistic example

This is a matrix from the Harwell-Boeing example collection, described as coming from the design of a nuclear power plant. I can’t comment on the fact that the matrix is 666 x 666. It came that way from the collection, I promise!

nnc666 = Import["~/courses/math4500/nnc666.rua"]

MatrixPlot[nnc666]

N[LUSparsity[nnc666]]

We can try out the RCMD and Sloan algorithms on this matrix.

{r,c} = MinimumBandwidthOrdering[nnc666,Method->"Sloan"];
MBOSloan = nnc666[[r,c]];
{L,U,p,cond} = ExplicitLUDecomposition[MBOSloan] ;
GraphicsGrid[{{MatrixPlot[L],MatrixPlot[U],MatrixPlot[MBOSloan]}}]
N[LUSparsity[MBOSloan]]

{r,c} = MinimumBandwidthOrdering[nnc666,Method->"RCMD"];
MBORCMD= nnc666[[r,c]];
{L,U,p,cond} = ExplicitLUDecomposition[MBORCMD] ;
GraphicsGrid[{{MatrixPlot[L],MatrixPlot[U],MatrixPlot[MBORCMD]}}]
N[LUSparsity[MBORCMD]]

As you can see, Sloan gives a little less fill on the LU matrices than RCMD (0.0669) vs (0.0791) and both beat the standard ordering (0.0835). We see that this is actually considerably better than we usually get from a random permutation:

RandomNNC := PermuteMatrix[RandomPermutation[666],RandomPermutation[666],nnc666];

MatrixPlot[RandomNNC]

SparsityData = Table[LUSparsity[RandomNNC],{i,1,100}];

Histogram[SparsityData]

Another example matrix

We need to define a little import cleanup to import our next matrix:

ImportExport`RegisterImport[
"RRMAT",
RRMAT`RRMATImport
]
RRMAT`RRMATImport[filename_String, options___]:=
Module[{strm,m,n,header,data,rows,cols},
strm = OpenRead[filename];
header = ReadList[strm,"String",5];
rows = ToExpression[StringCases[header[[4]],NumberString]][[1]];
cols = ToExpression[StringCases[header[[5]],NumberString]][[1]];
Print["Matrix of size: ",rows," x ",cols, " loaded.\n"];
data = Partition[ReadList[strm,Real],cols];
(* the converter returns a list of rules,
 in the form of ElementName -> ElementValue *)
{
"data" -> data
}
]

We now import this (a non-square matrix from an optimization problem with sparse constraints).

rrMatrix = "data" /. Import["~/courses/math4500/rrmatrix.mat","RRMat"];

MatrixPlot[rrMatrix]

We will see that when solving least squares problems involving matrices, a quick and dirty solution is to multiply A by A^T and solve a problem based on this.

sqMat = rrMatrix.Transpose[rrMatrix];
MatrixPlot[sqMat]

If we just take the standard LU decomposition (equivalent to the Cholesky decomposition), we get quite a bit of fill:

{L,U,p,cond} = ExplicitLUDecomposition[sqMat] ;
GraphicsGrid[{{MatrixPlot[L],MatrixPlot[U],MatrixPlot[sqMat]}}]
N[LUSparsity[sqMat]]

RCMD and Sloan each result in fill of about 0.22, which is half the fill from the standard ordering.

{r,c} = MinimumBandwidthOrdering[sqMat,Method->"RCMD"];
MBORCMD = sqMat[[r,c]];
{L,U,p,cond} = ExplicitLUDecomposition[MBORCMD] ;
GraphicsGrid[{{MatrixPlot[L],MatrixPlot[U],MatrixPlot[MBORCMD]}}]
N[LUSparsity[MBORCMD]]

{r,c} = MinimumBandwidthOrdering[sqMat,Method->"Sloan"];
MBOSloan= sqMat[[r,c]];
{L,U,p,cond} = ExplicitLUDecomposition[MBOSloan] ;
GraphicsGrid[{{MatrixPlot[L],MatrixPlot[U],MatrixPlot[MBOSloan]}}]
N[LUSparsity[MBOSloan]]

As before, Sloan is a little better at reducing fill in the LU decomposition, but both are better than a random permutation.

RandomSq := PermuteMatrix[RandomPermutation[498],RandomPermutation[498],sqMat];

MatrixPlot[RandomSq]

SparsityData = Table[LUSparsity[RandomSq],{i,1,100}];

Histogram[SparsityData]

You can see from this that a random permutation gives you something that’s very close to full fill.