[Numpy-discussion] how to deal with large arrays

[Darren Dale]

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