MathGroup Archive 1997

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

Search the Archive

More on: Newbie question: big matrix calculations

  • To: mathgroup at smc.vnet.net
  • Subject: [mg9298] More on: [mg9253] Newbie question: big matrix calculations
  • From: Luci Ellis <elisha at dot.net.au>
  • Date: Mon, 27 Oct 1997 02:47:08 -0500
  • Sender: owner-wri-mathgroup at wolfram.com

Dennis I have looked at your code and have realised there is a lot
quicker way to do it. See below:

Your original code:

In[1]:=
size =5;
In[2]:=
matrxD=Table[Switch[i-j,0,5,_,0],{i,size},{j,size}];

In[3]:=
matrxU=Table[Switch[i-j,-1,-3,_,0],{i,size},{j,size}];

In[4]:=
matrxL=Table[Switch[i-j,1,-2,_,0],{i,size},{j,size}];

In[5]:=
b=Table[1,{i,size}];
In[6]:=
x1=Table[1,{i,size}];
In[9]:=
x2 = -Inverse[matrxD].(matrxL+matrxU).x1+Inverse[matrxD].b;//Timing
Out[9]=
{0.716667 Second,Null}

My original suggestion about not inverting D twice doesn't buy you much
for small sizes, if anything you'll go backwards (-;

In[11]:=
x2a =
With[{dinv=Inverse[matrxD]},-dinv.(matrxL+matrxU).x1+dinv.b];//Timing

Out[11]=
{0.95 Second,Null}

Let's look at what's going on here:

In[10]:=
x2
Out[10]=
{4/5, 6/5, 6/5, 6/5, 3/5}

This is the first issue -- you are working in integers, with infinite
precision, but you are implementing what is in essence a numerical
procedure. So you let's try it in Real numbers not rationals.

In[19]:=
size = 5;

In[22]:=
matrxD=Table[Switch[i-j,0,5.,_,0.],{i,size},{j,size}]; In[23]:=
matrxU=Table[Switch[i-j,-1,-3.,_,0.],{i,size},{j,size}]; In[24]:=
matrxL=Table[Switch[i-j,1.,-2.,_,0.],{i,size},{j,size}]; In[16]:=
b=Table[1.,{i,size}];
In[17]:=
x1=Table[1.,{i,size}];
In[25]:=
x2 = -Inverse[matrxD].(matrxL+matrxU).x1+Inverse[matrxD].b;//Timing
Out[25]=
{0.05 Second,Null}

Much quicker! (But I got a different answer, I don't know why...)

In[26]:=
x2
Out[26]=
{0.8,0.8,0.8,0.8,0.2}

But there's more to it than that.  Looking at your code I see that what
you are trying to do is solve for the matrix X that satisfies the
equation: X = -inv(D).(L+U).x +inv(D).b
or if you like:  inv(D)-inv(D).(L+U).x +inv(D).b-X=0.

So why don't you Solve for the equation explicitly instead of using
repeated evaluation to find the solution

In[32]:=
X= Table[x[i],{i,size}]
Out[32]=
{x[1],x[2],x[3],x[4],x[5]}

In[36]:=
eqns=Thread[ -Inverse[matrxD].(matrxL+matrxU).X+Inverse[matrxD].b-X==0]
Out[36]=
{0.2 -1. x[1]+0.6 x[2]+0. x[3]+0. x[4]+0. x[5]==0,
  0.2 +0. x[1]-1. x[2]+0.6 x[3]+0. x[4]+0. x[5]==0,
  0.2 +0. x[1]+0. x[2]-1. x[3]+0.6 x[4]+0. x[5]==0,
  0.2 +0. x[1]+0. x[2]+0. x[3]-1. x[4]+0.6 x[5]==0,
  0.2 +0. x[1]+0. x[2]+0. x[3]+0. x[4]-1. x[5]==0} In[40]:=
solution=Solve[eqns,X]//Timing
Out[40]=
{0.0666667
    Second,{{x[1]->0.46112,x[2]->0.4352,x[3]->0.392,
      x[4]->0.32,x[5]->0.2}}}

Pretty fast considering I'm doing this on a three-year-old Mac Quadra
630! So you could write a function:

In[45]:=
solvesparsething[size_Integer]:=
Module[{matrxD,matrxU,matrxL,b,X},
	b=Table[1.,{i,size}];
	matrxD=Table[Switch[i-j,0,5.,_,0.],{i,size},{j,size}];
 	matrxU=Table[Switch[i-j,-1,-3.,_,0.],{i,size},{j,size}];
 	matrxL=Table[Switch[i-j,1.,-2.,_,0.],{i,size},{j,size}];
	X= Table[x[i],{i,size}];
	eqns=Thread[
-Inverse[matrxD].(matrxL+matrxU).X+Inverse[matrxD].b-X==0];
	Flatten[X/.Solve[eqns,X]]]

In[46]:=
solvesparsething[6]//Timing
Out[46]=
{0.766667 Second,{0.476672,0.46112,0.4352,0.392,0.32,0.2}} In[48]:=
solvesparsething[10]//Timing
Out[48]=
{1.
Second,{0.496977,0.494961,0.491602,0.486003,0.476672,0.46112,0.4352,0.392,
    0.32,0.2}}
In[49]:=
solvesparsething[50]//Timing
Out[49]=
{31.2833
Second,{0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,
    0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.499999,0.499999,0.499999,0.499998,
   
0.499996,0.499993,0.499989,0.499982,0.49997,0.499949,0.499915,0.499859,
   
0.499765,0.499608,0.499347,0.498912,0.498186,0.496977,0.494961,0.491602,
    0.486003,0.476672,0.46112,0.4352,0.392,0.32,0.2}} In[50]:=
solvesparsething[100]//Timing

Hope this helps,

Regards,
Luci

--------------
Luci Ellis:  elisha at dot.net.au
http://www.dot.net.au/~elisha



  • Prev by Date: Re: Newbie question: big matrix calculations
  • Next by Date: Re: Finding "start directory" with Mathematica 3.0 ?
  • Previous by thread: graphing
  • Next by thread: How to display *.bmp image in notebook?