Grabbing Overlapping Blocks and an Update Version of Jacobi

1) I do not need to use square blocks if I go the PseudoInverse route.  
2) I do not need to not overlap if I use the minimal update version

Here is Block Jacobi for a partitioned matrix A
	(Subscript[A, i,i ])(Subscript[x, i]^new)= (Subscript[b, i])-(∑_(j=1)^(i-1))(Subscript[A, i,j ])(Subscript[x, j]^old)-(∑_(j=i+1)^n)(Subscript[A, i,j ])(Subscript[x, j]^old)
Here is the update version for ∇(Subscript[x, i])=(Subscript[x, i]^new)-(Subscript[x, i]^old)
	(Subscript[A, i,i ])((Subscript[x, i]^old)+∇(Subscript[x, i]))= (Subscript[b, i])-(∑_(j=1)^(i-1))(Subscript[A, i,j ])(Subscript[x, j]^old)-(∑_(j=i+1)^n)(Subscript[A, i,j ])(Subscript[x, j]^old)
which gives
	(Subscript[A, i,i ])∇(Subscript[x, i])= (Subscript[b, i])-(∑_(j=1)^n)(Subscript[A, i,j ])(Subscript[x, j]^old)
which should be implemented as
	(Subscript[res, i])=(Subscript[b, i])-(∑_(j=1)^n)(Subscript[A, i,j ])(Subscript[x, j]^old)
	∇(Subscript[x, i])=LinearSolve[(Subscript[A, i,i]), (Subscript[res, i])]
or semi-equivalently
	(Subscript[res, i])=(Subscript[b, i])-(∑_(j=1)^n)(Subscript[A, i,j ])(Subscript[x, j]^old)
	(Subscript[x, i])+=LinearSolve[(Subscript[A, i,i]), (Subscript[res, i])]
Now as far as a pseudo inverse is involved the update version is different from the original version.  An undertermined pseudoinverse generates the minimal norm solution that gives a zero residual.  An overdetermined pseudoinverse generates a solution that minimizes the residual.

Test Matrices

I downloaded three FEA matrices (positive definate and probably symmetric) from Matrix Market.  I added stuff to the diagonals to make them diagonally dominant.

PorkDiagonals[A_]:= Module[
{m=Max[Abs[Diagonal[A]]]},
A+SparseArray[Band[{1,1}]->m, {Length[A],Length[A]}]
]
SetDirectory["/home/campus03/struther/Desktop/Link to Courses/5580/CodeFragments/LinearAlgebra/"];
names=FileNames["*mtx.gz"];
{A1,A2,A3}=Sort[Map[Import,names], Length[#1]<Length[#2]&];
{A1,A2,A3}=Map[PorkDiagonals, {A1,A2,A3}];
Map[MatrixPlot,{A1,A2,A3}]
{n1,n2,n3}=Map[Length, {A1,A2,A3}];
{u1,u2,u3}=Map[ Sin[Range[0.0,1,(1/(#-1))]]&, {n1,n2,n3}];
{b1,b2,b3}={A1.u1,A2.u2,A3.u3};
{Norm[LinearSolve[A1,b1]-u1],Norm[LinearSolve[A2,b2]-u2],Norm[LinearSolve[A3,b3]-u3]}

Apparently the condition number of some of the matrices is relatively large.

Jacobi and Related Iterations

The following defines a Jacobi Iteration with a specified order, an inplace Jacobi like Iterate and Gauss-Seidel iteration.  The standard order is just 1,2,3,4,…

JacobiIterate[xIn_,{{Dia_,LPlusU_},b_},order_]:=Module[
{i,j},
Table[
i=order[[j]];
(b[[i]]-LPlusU[[i]].x)/Dia[[i]],
{j,Length[Dia]}] 
]
SetAttributes[JacobiIterateInPlace, HoldAll]
JacobiIterateInPlace[x_,{{Dia_,LPlusU_},b_},order_]:=Module[
{i,j},
Do[
i=order[[j]];
x[[i]]=(b[[i]]-LPlusU[[i]].x)/Dia[[i]],
{j,Length[Dia]}] 
]
GaussSeidelIterate[xIn_,{{L_,DiaPlusU_},b_},order_]:=Module[
{},
LinearSolve[DiaPlusU,b-L.xIn]
]

All three of these should converge on the test problems.

NIter=23;
{TestMatrix,TestRHS,TestSol}={A1,b1,u1};
order=Range[1,Length[TestMatrix]];
Dia=Diagonal[TestMatrix];
LPlusU=TestMatrix-DiagonalMatrix[Dia];
L=LowerTriangularize[TestMatrix,-1];
DiaPlusU=TestMatrix-L;
x=ConstantArray[0,Length[TestMatrix]];
DataJ=Table[
x=JacobiIterate[x,{{Dia,LPlusU},TestRHS},order];
Norm[x-TestSol],
{NIter}];
x=ConstantArray[0,Length[TestMatrix]];
DataJIP=Table[
JacobiIterateInPlace[x,{{Dia,LPlusU},TestRHS},order];
Norm[x-TestSol],
{NIter}];
x=ConstantArray[0,Length[TestMatrix]];
DataJIPRandom=Table[
JacobiIterateInPlace[x,{{Dia,LPlusU},TestRHS},RandomSample[order,Length[TestMatrix]]];
Norm[x-TestSol],
{NIter}];
x=ConstantArray[0,Length[TestMatrix]];
FixedRandomOrder1=RandomSample[order,Length[TestMatrix]];
DataJIPFixedRandom1=Table[
JacobiIterateInPlace[x,{{Dia,LPlusU},TestRHS},FixedRandomOrder1];
Norm[x-TestSol],
{NIter}];
x=ConstantArray[0,Length[TestMatrix]];
FixedRandomOrder2=RandomSample[order,Length[TestMatrix]];
DataJIPFixedRandom2=Table[
JacobiIterateInPlace[x,{{Dia,LPlusU},TestRHS},FixedRandomOrder2];
Norm[x-TestSol],
{NIter}];
x=ConstantArray[0,Length[TestMatrix]];
DataGS=Table[
x=GaussSeidelIterate[x,{{L,DiaPlusU},TestRHS},order];
Norm[x-TestSol],
{NIter}];
ListLogPlot[{DataJ, DataJIP,DataGS,DataJIPRandom,DataJIPFixedRandom1,DataJIPFixedRandom2 },
PlotStyle->{Red, Blue,Green,Pink, Cyan, Black}]

ListLogPlot[{DataGS,DataJIPFixedRandom1,DataJIPFixedRandom2 },
PlotStyle->{Green,Cyan, Black}]

The order does not seem to matter too much provided it is consistent from run to run. I wonder howmuach it matters how consistent.

Block Ideas

I can implement block versions of these as well. I need to be able to block a sparse array and keeps the pieces sparse.  Unfortunately, Partition makes dense arrays but there is no real problem.

SparsePartition[A_SparseArray,{nSub_,mSub_}]:= Module[
{n,m},
{n,m}=Dimensions[A];
Table[
A[[i;;i-1+nSub,j;;j-1+mSub]],
{i, 1, n+1-nSub, nSub},
{j,1, m+1-mSub, mSub}
]
]

The thing above merely takes what fits.  I need to make my partitioning fit my problem.  Well actually I need to fix that.  This is just a test that the thing is working correctly.   It also shows that there are a lot of zero blocks.

SplitA1=SparsePartition[A1,{13,13}];
{Dimensions[A1],Dimensions[SplitA1]}
Norm[A1-ArrayFlatten[SplitA1]]
MatrixPlot[Map[Norm, SplitA1,{2}]]

To make things work in block land I need to block my right hand side and vectors.  I should also check that it is working correctly

Splitu1=Partition[u1,13];
Splitb1=Partition[b1,13];
test=Table[
Sum[SplitA1[[i,j]].Splitu1[[j]],{j,1,Length[Splitu1]}],
{i,1, Length[SplitA1]}
];
Norm[Flatten[test]-b1]

SetAttributes[JacobiIPBlock, HoldAll]
JacobiIPBlock[Splitx_,{{SplitADia_,SplitLPlusU_},Splitb_},order_]:=Module[
{i,j,n,m,RHS},
{n,m}=Take[Dimensions[SplitLPlusU],2];
Do[
i=order[[k]];
RHS= Sum[SplitLPlusU[[i,j]].Splitx[[j]],{j,1,i-1}]+ Sum[SplitLPlusU[[i,j]].Splitx[[j]],{j,i+1,m}];
RHS=Splitb[[i]]-RHS;
Splitx[[i]]=LinearSolve[SplitADia[[i]],RHS],
{k,n}
] 
]

SplitA1Dia=Diagonal[SplitA1];
SplitA1LPlusU = SplitA1;
Do[SplitA1LPlusU[[i,i]]=SparseArray[{},Dimensions[SplitA1LPlusU[[i,i]]]],
{i,1, Length[SplitA1LPlusU]}]

NIter=15;
Splitx=0*Splitu1;
order=Range[1,36];
DataJacobiIPBlock=Table[
JacobiIPBlock[Splitx,{{SplitA1Dia,SplitA1LPlusU},Splitb1},order];
Norm[Flatten[Splitx]-u1],
{NIter}];
Splitx=0*Splitu1;
DataJacobiIPBlockRO=Table[
JacobiIPBlock[Splitx,{{SplitA1Dia,SplitA1LPlusU},Splitb1},RandomSample[order,Length[order]]];
Norm[Flatten[Splitx]-u1],
{NIter}];
Clear[FixedRO, DataJacobiIPBlockFRO,p]
Do[
Splitx=0*Splitu1;
FixedRO[p]=RandomSample[order,Length[order]];
DataJacobiIPBlockFRO[p]=Table[
JacobiIPBlock[Splitx,{{SplitA1Dia,SplitA1LPlusU},Splitb1},FixedRO[p]];
Norm[Flatten[Splitx]-u1],
{NIter}],
{p,1,100}]
ListLogPlot[{DataJacobiIPBlock, DataJacobiIPBlockRO,DataJacobiIPBlockFRO[1],DataGS},
PlotStyle->{Red, Green,Pink, Blue}]

ListLogPlot[Table[DataJacobiIPBlockFRO[p],{p,1,100}],
PlotStyle->{Red, Green,Pink, Blue}]

Overlaping Blocks

This can now have overlapping blocks for the solving portion of things.  I need to do an updating version with a relaxation parameter.  First just making the update form.

SetAttributes[JacobiIPBlockUpdate, HoldAll]
JacobiIPBlockUpdate[Splitx_,{{SplitADia_,SplitA_},Splitb_},order_]:=Module[
{i,j,n,m,RHS},
{n,m}=Take[Dimensions[SplitA],2];
Do[
i=order[[k]];
RHS= Splitb[[i]]-Sum[SplitA[[i,j]].Splitx[[j]],{j,1,m}];
Splitx[[i]]+=LinearSolve[SplitADia[[i]],RHS],
{k,n}
] 
]

First done using the same blcoks as before.

NIter=15;
Splitx=0*Splitu1;
order=Range[1,36];
DataJacobiIPBlock=Table[
JacobiIPBlockUpdate[Splitx,{{SplitA1Dia,SplitA1},Splitb1},order];
Norm[Flatten[Splitx]-u1],
{NIter}];
Splitx=0*Splitu1;
DataJacobiIPBlockRO=Table[
JacobiIPBlockUpdate[Splitx,{{SplitA1Dia,SplitA1},Splitb1},RandomSample[order,Length[order]]];
Norm[Flatten[Splitx]-u1],
{NIter}];
Clear[FixedRO, DataJacobiIPBlockFRO,p]
Do[
Splitx=0*Splitu1;
FixedRO[p]=RandomSample[order,Length[order]];
DataJacobiIPBlockFRO[p]=Table[
JacobiIPBlockUpdate[Splitx,{{SplitA1Dia,SplitA1},Splitb1},FixedRO[p]];
Norm[Flatten[Splitx]-u1],
{NIter}],
{p,1,100}]
ListLogPlot[{DataJacobiIPBlock, DataJacobiIPBlockRO,DataJacobiIPBlockFRO[1],DataGS},
PlotStyle->{Red, Green,Pink, Blue}]

Here is the spread in the orders of the updates.

ListLogPlot[
Table[DataJacobiIPBlockFRO[i],{i,100}],
Joined->False
]

here is a twiddler with the resiudals at each stage.

Manipulate[
Histogram[ Table[DataJacobiIPBlockFRO[i][[j]],{i,100}],
PlotLabel->j,
PlotRange->All],
{j, 1, 15,1}]

Now if I go and make the blocks overlap I might see something cool!

Note: I need to change what I do with the vectors.

Clear[JacobiIPUpdatePInv]
SetAttributes[JacobiIPUpdatePInv, HoldAll]
JacobiIPUpdatePInv[x_,{A_,b_},ω_,{EqOrder_,VarOrder_}]:=Module[
{Eqs,Vars,k,ASub,RHS},
Do[
Eqs=EqOrder[[k]];
Vars=VarOrder[[k]];
ASub=A[[Eqs,Vars]];
RHS=( b-A.x)[[Eqs]];
x[[Vars]]+=ω PseudoInverse[ASub].RHS,
{k,Length[EqOrder]}
] 
]

Lets test the damm thing!

Map[Dimensions,{A1,b1,u1}]
Norm[A1.u1-b1]

NIter=15;
BlockSize=36;
Dim=Length[TestSol];
x=0*u1;
EqOrder=Table[Range[1+i*BlockSize,BlockSize(i+1)],{i, 0, Dim/BlockSize-1}];
VarOrder=EqOrder;
(*VarOrder=Table[Range[1,Length[TestSol]],{13}];*)
order=Range[1, Length[TestSol]];
order=RandomSample[order, Length[TestSol]];
Clear[DataJacobiIPUpdatePInv,ω]
{ωMin,ωMax,Δω}={0.95,1.05,0.001};
DataJacobiIPUpdatePInv=Table[
x=0*u1;
Table[
JacobiIPUpdatePInv[x,{TestMatrix,TestRHS},ω,{EqOrder,VarOrder}];
Norm[x-TestSol],
{NIter}],
{ω, ωMin,ωMax,Δω}];
x=0*u1;
DataJIP=Table[
JacobiIterateInPlace[x,{{Dia,LPlusU},TestRHS},order];
Norm[x-TestSol],
{NIter}];
ListLogPlot[ Join[{DataJIP},DataJacobiIPUpdatePInv],
PlotStyle->Join[
{{PointSize[0.03],Black}},
Table[Hue[i/Length[DataJacobiIPUpdatePInv]],{i,1,Length[DataJacobiIPUpdatePInv]}]]
]

ListLogPlot[ DataJacobiIPUpdatePInv[[All,-6]],
DataRange->{ωMin,ωMax},
AxesLabel->{"ω","Res"},
PlotLabel->BlockSize]

Break

Hyper Power algorithms

Two ideas.  
1) As an approximate sub roblem solver for a block algorithm.
2) As simple preconditioner.  Including sparse vector idea.

I am implementing the q=2 versions of Stanmirovic [1].

Clear[HyperPower,A,k,q]
HyperPower[A_,k_,q_]:= Module[
{m,n,Id,α,Y,T,M,j},
{m,n}=Dimensions[A];
Id =IdentityMatrix[n];
α=2/Tr[Transpose[A].A];
Y=α Transpose[A];
Do[
T[0]=Id - Y.A;
Do[T[j]=T[0].T[j-1],{j,1,q-1}];
M=Id+Sum[T[j],{j,0,q-1}];
Y=Normal[M.Y],
{k}];
Y
]

Clear[q]
n=32;
NIter=12;
A=RandomReal[{0,1},{n,n}]+ DiagonalMatrix[RandomReal[{2,3},{n}]];
PInvA=PseudoInverse[A];
ListLogPlot[
Table[
B=HyperPower[A,k,q];
Norm[B-PInvA],
{q,2,12},
{k, 1,NIter}
],
Mesh->All,
Joined->True,
PlotRange->All
]

Clear[q]
n=32;
NIter=12;
A=RandomReal[{0,1},{n,n}]+ DiagonalMatrix[RandomReal[{2,3},{n}]];
A=A+Transpose[A];
PInvA=PseudoInverse[A];
evals=Eigenvalues[PInvA]
ListLogPlot[
Table[
B=HyperPower[A,k,q];
Norm[B-PInvA],
{q,2,4},
{k, 1,NIter}
],
Mesh->All,
Joined->True
]
Table[
B=HyperPower[A,k,q];
ListPlot[{Eigenvalues[B]/evals},
PlotRange->{-3,3}],
{q,2,12},
{k, 1,NIter}
]

I have a random selection of mat market stuff here.

SetDirectory["/home/campus03/struther/Desktop/Link to Courses/5580/Stuff/"]
FileNames[]

A= Import["nos5.mtx.gz"]
{m,n}=Dimensions[A];
Timing[PInvA=PseudoInverse[A];]
Norm[A.PInvA-IdentityMatrix[m]]
evals=Eigenvalues[PInvA,50];
NIter=21;
ListLogPlot[
Table[
B=HyperPower[A,k,q];
Norm[B-PInvA],
{q,2,2},
{k, 1,NIter}
],
Mesh->All,
Joined->True,
PlotRange->All
]

Norm[A-Transpose[A]]

ListLogPlot[{
Eigenvalues[PInvA],
Eigenvalues[B]
}]

SparsifyVector[Vec_, nz_]:= Module[
{rules},
Sort[

n=1000;
nz=50;
NZ=10;
coords=RandomSample[Range[1,n],nz];
vals=RandomReal[NormalDistribution[ 0,23],{nz}];
b=SparseArray[ Table[coords[[i]]->vals[[i]],{i, nz}]];
Barney=ArrayRules[b];
ListLogPlot[
Table[
((Norm[SparseArray[Barney[[Ordering[ArrayRules[b],NZ,Abs[#1[[2]]]>Abs[#2[[2]]]&]]]]])/(Norm[b])),
{NZ,1,nz}],
PlotRange->All
]

Barney[[24]]

Order

FullForm[{246}->0.5232769192040148`]

FullForm[{698}->0.9454930310081582`]

A
B=A.Transpose[A];
B

195654/(1806.^2)
195654/63454.