[Numpy-discussion] how to deal with large arrays
Darren Dale
dd55 at cornell.edu
Fri Oct 15 16:29:03 CDT 2004
Thank you for your response, Tim,
On Friday 15 October 2004 06:07 pm, Tim Hochberg wrote:
> Darren Dale wrote:
> >Hello,
> >
> >I have two 2D arrays, Q is 3-by-q and R is r-by-3. At each q, I need to
> > sum q(r) over R, so I take a dot product RQ and then sum along one axis
> > to get a 1-by-q result.
> >
> >I'm doing this with dot products because it is much faster than the
> > equivalent for or while loop. The intermediate r-by-q array can get very
> > large though (200MB in my case), so I was wondering if there is a better
> > way to go about it?
>
> I'm fairly certain (but I urge you to double check), that this reduces to:
>
> result_2 = na.dot(na.sum(R, 0), Q)
>
Yes. As usual, I left out a bit of information that turned out to be
important. See below
A modified test:
from numarray import *
from numarray import random_array
q = 10
r = 12
R = random_array.random((r,3))
Q = random_array.random((3,q))
x1 = sum( exp(1j*dot(R,Q)), 0) #note complex argument to exp()
x2 = exp(1j*dot(sum(R, 0), Q))
print allclose(x1, x2)
The complex arithmetic changes things. I am still learning how to keep my code
efficient. The following code is actually almost as fast as using the large
dot product, apparently I had some other sinks in my original tests:
phase = zeros(len(Q[0]),'d')
for i in range(len(Q[0])):
phase[i] = phase[i] + sum(exp(1j*dot(R,Q[:,i])), 0)
If q=1000 and r=2500, the for loop takes about 13% longer than the dot product
method. Incredibly, if q=10,000 and r=2500, the for loop is 17% faster. So I
am going to use it instead. Apparently I had some other time sink in my
original test.
from numarray import *
from numarray import random_array
from time import clock
q = 10000
r = 2500
R = random_array.random((r,3))
Q = random_array.random((3,q))
t0 = clock()
x1 = sum(exp(1j*dot(R,Q)), 0) #note complex argument to exp()
t1 = clock()
dt1 = t1-t0
phase = zeros(len(Q[0]),'d')
for i in range(len(Q[0])):
phase[i] = phase[i] + sum(exp(1j*dot(R,Q[:,i])), 0)
t2 = clock()
dt2 = t2-t1
print (dt2-dt1)/dt1
--
Darren
More information about the Numpy-discussion
mailing list