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