Re: Speeding Up Numerical Matrix Calculations ???
- To: mathgroup at christensen.cybernetics.net
- Subject: [mg1124] Re: Speeding Up Numerical Matrix Calculations ???
- From: wagner at bullwinkle.cs.Colorado.EDU (Dave Wagner)
- Date: Mon, 15 May 1995 02:41:49 -0400
- Organization: University of Colorado, Boulder
In article <3ouung$995 at news0.cybernetics.net>,
Paul E Howland <PEHOWLAND at taz.dra.hmg.gb> wrote:
>My question is quite simple. Mathematica can be numerically very slow.
>How can I speed up my numerical calculations involving matrices that
>consist only of machine size numbers?
>
>It appears that Mathematica handles numerical matrices in much the same
>way as it does symbolic matrices. Here's a test I did using Mathematica
>V2.2.3 on a 486/33 PC running Windows 3.1.
>
>In[1]:= F[a_?MatrixQ, b_?MatrixQ, c_?MatrixQ] :=
> Transpose[a]a + 3 b c + Transpose[c]c
>In[2]:= NF[a_?MatrixQ, b_?MatrixQ, c_?MatrixQ] :=
> N[Transpose[a]a + 3 b c + Transpose[c]c]
>
>Now generate some machine sized numerical matrices ...
>
>And some symbolic matrices ...
>
>And compare timings. First compare F[] on both sets of data ...
>
>In[9] := Timing[Do[F[a,b,c], {i,10}]]
>Out[9] = {3.185 Second, Null}
>In[10] := Timing[Do[F[x,y,z], {i,10}]]
>Out[10] = {2.087 Second, Null}
>
>thus the symbolic processing is faster (!).
It certainly appears that ordinary multiplication is not optimized for
machine-precision arguments. Dot-product, however, is:
In[30]:=
F2[a_?MatrixQ, b_?MatrixQ, c_?MatrixQ] :=
Transpose[a].a + 3 b.c + Transpose[c].c
In[31]:=
Timing[Do[F2[a,b,c], {i,10}]]
Out[31]=
{0.583333 Second, Null}
In[32]:=
Timing[Do[F2[x,y,z], {i,10}]]
Out[32]=
{6.31667 Second, Null}
(Assuming what you really want to do is a dot-product.)
>It is not possible to use Compile[] on a function that accepts matrices
>as arguments because Compile[] only accepts _Integer, _Real and _Complex
>as its data types (and not _List).
You can sometimes get a speedup just by compiling very simple
functions, in this case Times:
In[43]:=
fasttimes = Compile[{{x, _Real}, {y, _Real}},
x * y
]
Out[43]=
CompiledFunction[{x, y}, x y, -CompiledCode-]
In[44]:=
fasttimes[3., 4.]
Out[44]=
12.
Now we need a fast way to apply this function pairwise to two matrices.
In[57]:=
matrixtimes[m1_, m2_] :=
MapThread[fasttimes, {m1, m2}, 2]
In[58]:=
matrixtimes[{{1.,2.},{3.,4.}},{{5.,6.},{7.,8.}}]
Out[58]=
{{5., 12.}, {21., 32.}}
Here's another version of F that uses matrixtimes rather than Times:
In[59]:=
F3[a_?MatrixQ, b_?MatrixQ, c_?MatrixQ] :=
matrixtimes[Transpose[a],a] +
3 matrixtimes[b,c] +
matrixtimes[Transpose[c],c]
It works:
In[67]:=
F[a,b,c] == F3[a,b,c]
Out[67]=
True
But it's apparently not much faster:
In[12]:=
Timing[Do[F[a,b,c], {i,10}]]
Out[12]=
{0.866667 Second, Null}
In[61]:=
Timing[Do[F3[a,b,c], {i,10}]]
Out[61]=
{0.75 Second, Null}
However, if the matrices are larger, the overheads of the initial function
call, etc., will not play so large a role in the timings:
In[72]:=
a2 = Table[Random[], {i,100}, {j,100}];
In[73]:=
b2 = Table[Random[], {i,100}, {j,100}];
In[74]:=
c2 = Table[Random[], {i,100}, {j,100}];
In[75]:=
Timing[F[a2,b2,c2];]
Out[75]=
{9.5 Second, Null}
In[76]:=
Timing[F3[a2,b2,c2];]
Out[76]=
{6.75 Second, Null}
Admittedly, this isn't terrific. Supposedly the compiler will be
substantially enhanced in the next version of Mma. Maybe someone
at WRI would care to comment on that?
>It would be nice if Mathematica was written such that
>when N[] was wrapped around simple numerical matrix operations, optimised
>numerical code was invoked.
N is the function that everyone uses but hardly anyone really understands.
The way N[expr] works is that first, expr is evaluated. If the result
is a normal expression (i.e., head[arg, arg, arg...]), *then* N looks
for rules of the form N[head[...]] := ... to apply; such rules are called
NValues (big surprise there!) Finally, N numericalizes every number it
can find in the result. So you see, since the expression inside the N
in your case evaluates before N can do anything, it's not possible for
"optimised numerical code" to be invoked. It is necessary for a
separate function to be written, say NTimes. In fact, this is the
reason that there are functions like NSum, NIntegrate, NRoots, etc: if
you call N[Integrate[expr, ...]] on an expression that can't be
integrated, first Integrate tries and fails, then N sees that the head
of the result is Integrate and so NIntegrate is called. It's much
faster if you call NIntegrate in the first place.
"Selected Tutorial Notes", available from WRI. Buy 'em, you'll be
happy you did.
Dave Wagner
Principia Consulting
(303) 786-8371
dbwagner at princon.com
http://www.princon.com/princon