LU Decomposition

First here is a complete LU decomposition set up in the notation and language of the link.

LU[A_,0]:=Module[{LU=A, n=Length[A],d,e,r,i,j},
Do[
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
Do[
LU[[i,j]]-= e* LU[[r,j]],
{j,r+1,n}
],
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]

n=3;
A=RandomReal[{-1,1},{n,n}];
MatrixForm[A]
{L,U}=LU[A,0];
{MatrixForm[L],MatrixForm[U]}
MatrixForm[L.U]

Here is a version vectorized on the interior (column vector) loop.

LU[A_,1]:=Module[{LU=A, n=Length[A],d,e,r,i,j},
Do[
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
LU[[i,r+1;;n]]-= e* LU[[r,r+1;;n]],
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]

Here is a version vectorized on the next inner most loop.  This is now a block sub-matrix update with an outer product.

LU[A_,2]:=Module[{LU=A, n=Length[A],d,e,r,i,j},
Do[
d=1/LU[[r,r]];
e = d LU[[r+1;;n,r]];
LU[[r+1;;n,r]]=e;
LU[[r+1;;n,r+1;;n]]-=Outer[Times,e,LU[[r,r+1;;n]]],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]

n=3;
A=RandomReal[{-1,1},{n,n}];
MatrixForm[A]
{L,U}=LU[A,0];
{MatrixForm[L],MatrixForm[U]}
MatrixForm[L.U]

Lets see which one is fastest.  We will include the built-in LUDecompostion in the timing.

Data=Table[
n=Floor[50 (1.5^p)];
A=RandomReal[{-1,1},{n,n}];
{
{n,Timing[LU[A,0];][[1]]}, 
{n,Timing[LU[A,1];][[1]]},
{n,Timing[LU[A,2];][[1]]}
},
{p,0, 5}]
ListLogLogPlot[
Transpose[Data],
Joined->True,
Mesh->All,
FrameLabel->{"n","Sec"},
GridLines->Automatic,
Frame->True,
PlotStyle->{Red, Green, Blue}
]

Data=Table[
n=Floor[50 (1.5^p)];
A=RandomReal[{-1,1},{n,n}];
{ 
{n,Timing[LU[A,1];][[1]]},
{n,Timing[LU[A,2];][[1]]},
{n,Timing[LUDecomposition[A];][[1]]}
},
{p,0, 9}]
ListLogLogPlot[
Transpose[Data],
Joined->True,
Mesh->All,
FrameLabel->{"n","Sec"},
GridLines->Automatic,
Frame->True,
PlotStyle->{Red, Green, Blue}
]

Data=Table[
n=Floor[50 (1.5^p)];
A=RandomReal[{-1,1},{n,n}];
{ 
{n,Timing[LU[A,1];][[1]]},
{n,Timing[LUDecomposition[A];][[1]]}
},
{p,5, 10}]

ListLogLogPlot[
Join[ {{{300, 0.1}, {4000,0.1 ((((4000/300)))^2)}}},Transpose[Data]],
Joined->True,
Mesh->All,
FrameLabel->{"n","Sec"},
GridLines->Automatic,
Frame->True,
PlotStyle->{Red, Green, Blue}
]

Incomplete LU Decomposition

The ILU description (which is ILU(0)) at  http://www.cfd-online.com/Wiki/Incomplete_LU_factorization_-_ILU is simple to understand and fairly simple implement.  You only update entries in the positions where A has non-zeros. The easiest to alter is the block updating version "2".   All we need to do is fix the updating to only do the

Now lets work on the ILU part.  To do this in the form on the web we need to find out where a given sparse matrix has entries. This is pretty simple but if we start working from the LU[A,0]  version we test for A[[i,j]]==0 a lot.  To keep the code simple I want to build a matrix that has ones in all the entries that A has non-zeroes. I am going to call such a matrix a “pattern”.

SparsePattern[AIn_SparseArray]:= Module[
{RulesForA},
RulesForA=ArrayRules[AIn];
SparseArray[
Append[Table[
RulesForA[[i,1]]-> 1,
{i,1, Length[RulesForA]-1}],
RulesForA[[-1]]
]
]
]

Clear[ILU,A]
ILU[A_,2]:=Module[
{LU=A, Patt,n=Length[A],d,e,r, Skeleton},
Patt=SparsePattern[A];
Do[
d=1/LU[[r,r]];
e = d LU[[r+1;;n,r]];
LU[[r+1;;n,r]]=e;
LU[[r+1;;n,r+1;;n]]-=Patt[[r+1;;n,r+1;;n]]*Outer[Times,e,LU[[r,r+1;;n]]];,
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]

A=Import["http://math.nist.gov/MatrixMarket/data/Harwell-Boeing/lanpro/nos5.mtx.gz"]
MatrixPlot[A]

Timing[{L,U}=LU[A,2]]
Timing[{IL,IU}=ILU[A,2]]

ADense=Normal[A];
Timing[{LDense,UDense}=LU[ADense,2];]
{MatrixPlot[A],MatrixPlot[L],MatrixPlot[U]}

{MatrixPlot[IL],MatrixPlot[IU]}

The hope is that IL.IU is close to A in an eigenvalue sense.

Id=Inverse[A].L.U;
ListLogPlot[Eigenvalues[Id],
PlotRange->All]

AppId=Inverse[A].IL.IU;
ListLogPlot[Eigenvalues[AppId],
PlotRange->All]

A simpler but much less good preconditioner is a “diagonal” preconditioner

DiaAppId=DiagonalMatrix[1/Normal[Diagonal[A]]].A;
ListLogPlot[Eigenvalues[DiaAppId],
PlotRange->All]

The hope is that Hyper power algorithms converge rapidly for

SparseIdentity[n_]:=SparseArray[ Band[{1,1}]->1,{n,n}]
HyperPowerIteration[X_,B_]:= 2 X - X.B.X
Sparsity[A_]:=(Length[ArrayRules[A]]-1.0)/Apply[Times,Dimensions[A]]

σ=SingularValueList[A,1][[1]];
α=1.9/(σ^2)
X=α Transpose[A];
ListLogPlot[
Table[
X=HyperPower[X,A];
Abs[Eigenvalues[A.X]],
{4}],
PlotStyle->{Red,Green,Blue,Cyan}
]

B=Inverse[IL.IU].A;
σ=SingularValueList[B,1][[1]];
α=2.0/(σ^2)
X=α Transpose[B];
ListLogPlot[
Table[
X=HyperPower[X,B];
Abs[Eigenvalues[B.X]],
{12}],
PlotStyle->{Red,Green,Blue,Cyan}
]

B=Inverse[IL.IU].A;
X= SparseIdentity[Length[B]]
λs[0]=Eigenvalues[X.Inverse[IL.IU].A];
Data=Table[
X=HyperPower[X,B];
λs[i]=Eigenvalues[X.Inverse[IL.IU].A];
1-Abs[λs[i]],
{i,1,12}];

ListLogPlot[
Data,
PlotStyle->{Red,Green,Blue,Cyan},
PlotRange->All,
Frame->True,
GridLines->Automatic
]

ListPlot[λs[3],
PlotRange->All]

Dimensions[Data]

Li

468^2

Break

Sparsity[A_SparseArray,f]:=Module[
{mn,nz},
mn=Apply[Times,Dimensions[A]];
nz=Length[ArrayRules[A]]-1;
N[(nz/mn)]
]
Sparsity[A_SparseArray,c]:=Module[
{mn,nz},
mn=Apply[Times,Dimensions[A]];
nz=Length[ArrayRules[A]]-1;
nz
]
RowSparsity[A_SparseArray]:=Module[
{nz,n},
n=Length[A];
nz=Table[ Sparsity[A[[i]],c],{i,1,n}]
]

ListLogPlot[{RowSparsity[A],RowSparsity[A.A]},
PlotRange->All,
Frame->True,
GridLines->Automatic
]

Sparsity[A]
Sparsity[A.A]

Head[A][[0]]

{IL,IU}=ILU[A,0];
IA=IL.IU;
ApproxId=Inverse[IU].Inverse[IL].A;
ListLogPlot[
Eigenvalues[ApproxId],
PlotRange->All]

SparseCombine[A_,

NZAs=Drop[ArrayRules[A],-1];
NZIAs=Drop[ArrayRules[IA],-1];

ΔA=SparseArray[
Table[ NZs[[i,1]]-> 1-((NZIAs[[i,2]])/(NZAs[[i,2]])),{i,1, Length[NZAs]}]
]
MatrixPlot[ΔA]
Max[Abs[ΔA]]
Norm[ΔA]

MatrixPlot[
IA-A
]