LinearSolve for underdetermined sparse systems
- To: mathgroup at smc.vnet.net
- Subject: [mg73153] LinearSolve for underdetermined sparse systems
- From: Veit Elser <ve10 at cornell.edu>
- Date: Sun, 4 Feb 2007 23:57:48 -0500 (EST)
The current version of Mathematica does not support the solution of
linear, underdetermined, sparse systems: A.x = b. The solution x is
unique when x is required to have the minimum Euclidean norm. For non-
sparse matrices one can use the Moore-Penrose inverse function: x =
PseudoInverse[A].b. As far as I know, Mathematica has no functions
yet for solving large sparse systems of this kind. Since this seems
like a common application, I provide below a conjugate-gradient based
solver I wrote.
Veit Elser
conjGradSolve[A_,b_]:=Module[{AA,n,x,r,p,a,d0,d1},
AA=Transpose[A].A;
n=Length[AA];
x=Table[0,{n}];
r=b.A;
p=r;
d1=r.r;
While[d0=d1;Log[10,d0/n]/2>-$MachinePrecision,
a=d0/(p.AA.p);
x+=a p;
r-=a AA.p;
d1=r.r;
p=r+(d1/d0)p;
];
x]
vars=20000;
eqns=10000;
sparseData=
Table[{Table[Random[Integer,{1,vars}],{3}],
Table[2Random[Integer]-1,{3}]},{eqns}];
A=SparseArray[
Flatten[Table[MapThread[{i,#1}\[Rule]#2&,sparseData[[i]]],
{i,eqns}],
1],{eqns,vars}];
b=Table[Random[]-.5,{eqns}];
Timing[x=conjGradSolve[A,b];]
{1.42867 Second,Null}
Norm[A.x-b]
\!\(2.4235873983871594`*^-14\)