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

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[
{i,eqns}],
1],{eqns,vars}];

b=Table[Random[]-.5,{eqns}];