MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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\)


  • Prev by Date: Re: Text Function with Graphic
  • Next by Date: Re: Text Function with Graphic
  • Previous by thread: Falloon's SpheroidalHarmonics package
  • Next by thread: Re: Questions about Discrete Fourier Transform