Possible To Calculate A Double Sum In Python Using List Comprehensions To Replace Both For-loops?

I have a double sum which reads basically sum = exp( x^2 + y^2 ) Of course, I could simply use two nested for loops but that tends to be time consuming for large numbers. I can use

Solution 1:

Why don't you do it the numpy way... without for loops:

x = np.arange(N_x)
y = np.arange(N_y)

xx, yy = np.meshgrid(x, y)
result = np.sum(np.exp((xx/N_x)**2 + (yy/N_y)**2))

Solution 2:

result = sum(np.exp(x**2 + y**2) for x in range(N_x) for y in range(N_y))

Solution 3:

You were almost there. The full sum can be written as the product of two 1D sums, i.e. (sum exp x^2) * (sum exp y^2):

>>> import numpy as np                      
>>>                       
>>> N_x = N_y = 100                                      
>>> 
# brute force                         
>>> result_1 = .0
>>> for x in xrange(N_x):
...     for y in xrange(N_y):
...         result_1 += np.exp( (float(x)/N_x)**2 + (float(y)/N_y)**2 )
... 
>>> result_1
21144.232143358553
>>> 
# single product method
>>> from __future__ import division
>>> 
>>> x, y = np.arange(N_x) / N_x, np.arange(N_y) / N_y
>>> np.exp(x*x).sum() * np.exp(y*y).sum()
21144.232143358469

My guess is this you can even do with list comp and beat the brute force numpy method:

>> rx, ry = 1.0 / (N_x*N_x), 1.0 / (N_y*N_y)
>>> sum([np.exp(rx*x*x) for x in xrange(N_x)]) * sum([np.exp(ry*y*y) for y in xrange(N_y)])
21144.232143358469

Indeed, timings done in Python3 because I don't know how to use timeit in Python2:

>>> from timeit import repeat
>>> 
>>> kwds = dict(globals=globals(), number=100)
>>> 
# single product - list comp
>>> repeat('sum(np.exp(rx*x*x) for x in range(N_x)) * sum(np.exp(ry*y*y) for y in range(N_y))', **kwds)
[0.0166887859813869, 0.016465034103021026, 0.016357041895389557]
>>>
# numpy brute force
>>> repeat('np.exp(np.add.outer(x*x, y*y)).sum()', **kwds)
[0.07063774298876524, 0.0348161740694195, 0.02283189189620316]

Obviously, numpy single product is even faster

>>> repeat('np.exp(x*x).sum() * np.exp(y*y).sum()', **kwds)
[0.0031406711786985397, 0.0031003099866211414, 0.0031157969497144222]