LCM CA PRNG

TestImage=ImageTake[
Import["C:\\Users\\Public\\Pictures\\Sample Pictures\\Tree.jpg"],
{257,512},
{257,512}];
NMax=1024;
NRun=45;
A0=Round[NMax*ImageData[ColorConvert[TestImage,"Grayscale"]]];
Ker=RandomInteger[{0,NMax},{4,4}];
A=A0;
Data=Table[
A=Mod[ListConvolve[Ker,A,1],NMax],
{NRun}];
(*Manipulate[
MatrixPlot[Data[[i]]],
{i,1, Length[Data],1}]*)

p=Length[IntegerDigits[NMax,2]];
Length[Fred=Flatten[Map[Drop[IntegerDigits[#,2,p],1]&, Flatten[Data]]]]

Might be a good idea to drop the most siginifcant digit!

Histogram[Map[Length,Split[Take[Fred, {12000, 24000}]]]]

ListPlot[Take[Fred,240]]

Manipulate[
Histogram[Flatten[Data[[i]]]],
{i,1, Length[Data],1}]

Histogram[Flatten[A]]

Jacobi Like Iterations.

Making a Sparse Array

n=5;
A= SparseArray[
Band[{1,1}]:>RandomReal[{1,2}],
{n,n}];
{i,j}=RandomInteger[{1,n},{2}];
Do[
{i,j}=RandomInteger[{1,n},{2}];
If[i≠j,A[[i,j]]+=RandomReal[{-1/n,0}]],
{n}];
A=A+Transpose[A]
MatrixPlot[A];
ListLogPlot[Abs[Eigenvalues[A]],
PlotRange->All]

RandomSparseMatrix[n_,Sparsity_]:= Module[
{A},
 A= SparseArray[
Band[{1,1}]:>RandomReal[{1,2}],
{n,n}];
If[Sparsity<1/n, Return[Print["Idiot"]]];
Do[
{i,j}=RandomInteger[{1,n},{2}];
If[i≠j,a=RandomReal[{0,1/n}];A[[i,j]]+=a;A[[j,i]]+=a],
{Floor[n*(1/n - Sparsity)]}
];
While[
((Length[ArrayRules[A]]-1)/(n^2))≤Sparsity,
{i,j}=RandomInteger[{1,n},{2}];
If[i≠j,a=RandomReal[{0,1/n}];A[[i,j]]+=a;A[[j,i]]+=a]
];
A
]

JacobiIterate[A_,x_,b_]:=Module[
{n=Length[A], Dia,LPlusU},
Dia=Diagonal[A];
LPlusU=A-DiagonalMatrix[Dia];
Table[
(b[[i]]-LPlusU[[i]].x)/Dia[[i]],
{i,n}] 
]
JacobiIterateRandomOrder[A_,x_,b_]:=Module[
{n=Length[A],order, Dia,y=0*x},
order=RandomSample[Range[n]];
Dia=Diagonal[A];
LPlusU=A-DiagonalMatrix[Dia];
Table[
i=order[[j]];
y[[i]]=(b[[i]]-LPlusU[[i]].x)/Dia[[i]],
{j,n}] ;
y
]
SetAttributes[JacobiIterateInPlace, HoldAll]
JacobiIterateInPlace[A_,x_,b_]:=Module[
{n=Length[A], Dia},
Dia=Diagonal[A];
LPlusU=A-DiagonalMatrix[Dia];
Do[
x[[i]]=(b[[i]]-LPlusU[[i]].x)/Dia[[i]],
{i,n}];
x
];
SetAttributes[JacobiIterateInPlaceRO, HoldAll]
JacobiIterateInPlaceRO[A_,x_,b_]:=Module[
{n=Length[A], Dia,order,LPlusU},
order=RandomSample[Range[n]];
Dia=Diagonal[A];
LPlusU=A-DiagonalMatrix[Dia];
Do[
i=order[[j]];
x[[i]]=(b[[i]]-LPlusU[[i]].x)/Dia[[i]],
{j,n,1,-1}];
x
]
GaussSeidelIterate[A_,x_,b_]:=Module[
{n=Length[A], DiaPlusU,L},
L=LowerTriangularize[A,-1];
DiaPlusU=A-L;
LinearSolve[DiaPlusU,b-L.x]
]

Testing the Convergence

n=400;
A=RandomSparseMatrix[n,0.1];
b=RandomReal[{-1,1},{n}];
x=RandomReal[{-1,1},{n}];

x1=JacobiIterate[A,x,b];
x2=JacobiIterateRandomOrder[A,x,b];
x1==x2
x3=GaussSeidelIterate[A,x,b];

NIter=5;
xx=LinearSolve[A,b];
x=0*b;
DataJ=Table[
x=JacobiIterate[A,x,b],
{NIter}];
x=0*b;
DataJIP=Table[
JacobiIterateInPlace[A,x,b],
{NIter}];
x=0*b;
DataJIPRO=Table[
JacobiIterateInPlaceRO[A,x,b],
{NIter}];
x=0*b;
DataGS=Table[
x=GaussSeidelIterate[A,x,b],
{NIter}];
(*order=RandomSample[Range[n]];
x=0*b;
Ap=Permute[A,order];
bp=Permute[b,order];
DataJIPRandomizedOnce=Table[
JacobiIterateInPlaceRO[Ap,x,bp],
{12}];*)

ColorList={Red, Blue,Green,Purple};
Res=Transpose[
Table[
{
Norm[DataJ[[i]]-xx],
Norm[DataJIP[[i]]-xx], 
Norm[DataJIPRO[[i]]-xx],
Norm[DataGS[[i]]-xx](*,
Norm[DataJIPRandomizedOnce[[i]]-Permute[xx,order]]*)
},{i,Length[DataJ]}]
];
Res=Res[[All,3;;-1]];

Res=Map[Transpose[{Range[3,NIter],#}]&, Res];
Dots=ListLogLogPlot[
Res,
PlotStyle->ColorList,
PlotRange->All,
Frame->True]

Res

Clear[x]
f[x_]=Map[Fit[#,{1,x},x]&,Log[ Res]]
Fits=LogLogPlot[Exp[f[Log[x]]],{x,1, NIter},
AspectRatio->1,
PlotStyle->ColorList]
Show[Dots,Fits]

f[1.2]

Block Computations

BlockSize=4;
NumBlock=12;
MatrixSize=BlockSize*NumBlock;
A=RandomReal[{-1,1},{MatrixSize,MatrixSize}];
Do[A[[i,i]]+=RandomReal[{12, 16}],{i,1, MatrixSize}]
APartitioned=Partition[A,{BlockSize,BlockSize}];

A block Jacobi looks like.

NIter=12;
xPartitioned=Table[0,{NumBlock},{BlockSize}];
bPartitioned=RandomReal[{-1,1},{NumBlock,BlockSize}];
b=Flatten[bPartitioned];
i=2;
Do[
Do[
xPartitioned[[i]]=LinearSolve[APartitioned[[i,i]],bPartitioned[[i]]-Sum[APartitioned[[i,j]].xPartitioned[[j]],{j,1,i-1}]
-Sum[APartitioned[[i,j]].xPartitioned[[j]],{j,i+1,NumBlock}]
],
{i,1,NumBlock}],
{NIter}
]

SparsePartition[A_,{nSub_,mSub_}]:= Module[
{i,j, nA,mA},
{nA,mA}=Dimensions[A];
Table[
A[[1+i*nSub;;(i+1)*nSub, 1+j*mSub;;(j+1)*mSub]],
{i,0, Floor[nA/nSub]-1},
{j,0, Floor[mA/mSub]-1}
]
]

Ap=SparsePartition[A,{4,4}];
AA=ArrayFlatten[Ap]
Norm[A-AA]

SetAttributes[BlockJacobiIterateInPlace, HoldAll]
BlockJacobiIterateInPlace[Ap_,xp_,bp_]:=Module[
{n=Length[Ap]},
Do[
xp[[i]]=LinearSolve[Ap[[i,i]],bp[[i]]-Sum[Ap[[i,j]].xp[[j]],{j,1,i-1}]-Sum[Ap[[i,j]].xp[[j]],{j,i+1,n}]
],
{i,n}];
xp
];

400/16

NIter=7;
nSub=8;
bp=Partition[b,nSub];
Ap=SparsePartition[A,{nSub, nSub}];
xSol=LinearSolve[A,b];
xp=ConstantArray[0,Dimensions[bp]];
Dimensions[xp]
ListLogLogPlot[
Table[
Norm[xSol-Flatten[BlockJacobiIterateInPlace[Ap,xp,bp]]],
{NIter}]
]

Dimensions[bp]

Partitioning Sparse Arrays

I am pretty sure that partition will work just fine on a SparseArray. Bad news (well sort of) is that the sub matrices are dense stored.

n=24;
A=RandomSparseMatrix[n,0.1];
ASplit=Partition[A,{6,6}];
AA=ArrayFlatten[ASplit];
Norm[A-AA]

I can fix the sparsity issue.

SparseArray[ASplit[[1,1]]]

Rolling PseudoInverse

n=5;
A=RandomReal[{-1,1},{n,n}];
MatrixForm[APInv=PseudoInverse[A]];
B=A;
B[[1]]=ConstantArray[0,n];
B[[All,1]]=B[[1]];
B[[1,1]]=0.0;
MatrixForm[BPInv=PseudoInverse[B]];
MatrixForm[Hmm=Chop[BPInv.A]]

Rank One and Rank Two updates.

Do[Hmm[[i,i]]=0,{i,n}];
MatrixRank[Hmm]

Understanding the packed sparse storage scheme

B={{0, 5, 13, 22, 30, 37, 45, 53, 61, 68, 74}, {{1}, 
   {2}, {6}, {8}, {9}, {1}, {2}, {3}, {4}, {5}, {7}, 
   {9}, {10}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, 
   {10}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {10}, {2}, 
   {3}, {4}, {5}, {6}, {7}, {8}, {1}, {3}, {4}, {5}, 
   {6}, {7}, {8}, {10}, {2}, {3}, {4}, {5}, {6}, {7}, 
   {8}, {9}, {1}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, 
   {1}, {2}, {3}, {7}, {8}, {9}, {10}, {2}, {3}, {4}, 
   {6}, {9}, {10}}};
Dimensions[A]
MatrixForm[A]

Block

A1=Import[ "ExampleData/bcsstk08.mtx.gz"]
A2=Import["ExampleData/mcca.rua.gz"]
ListLogPlot[Eigenvalues[A2],
PlotRange->All
]

$HistoryLength=0;
SetDirectory["C:\\Program Files\\Wolfram Research\\Mathematica\\8.0\\AddOns\\ExtraPackages\\LinearAlgebraExamples\\Data"]
FileNames[]
TestMatFiles={"cavity01.rua","dwg961b.mtx","gr_30_30.rsa","nos6.mtx","west0381.mtx"};
Mats=Map[Import,TestMatFiles]
Map[Chop,Map[Det, Mats]]

A=Mats[[-1]];
ASplit=Map[Sparse,Partition[A, {16,16},1],2];
ReconsA=ArrayFlatten[ASplit];
Map[ArrayPlot,{A,ReconsA}]

ListPlot[Sort[Table[Det[ASplit[[i,i]]],{i,1,Length[ASplit]}]]]