Infinite Summation in Python

Question

I have a function for with i need to do an infinite summation on (over all the integers) numerically. The summation doesn't always need to converge as I can change internal parameters. The function looks like,

m(g, x, q0) = sum(abs(g(x - n*q0))^2 for n in Integers)
m(g, q0) = minimize(m(g, x, q0) for x in [0, q0])

using a Pythonic pseudo-code

Using Scipy integration methods, I was just flooring the n and integrating like for a fixed x,

m(g, z, q0) = integrate.quad(lambda n:
                             abs(g(x - int(n)*q0))**2,
                             -inf, +inf)[0]

This works pretty well, but then I have to do optimization on the x as a function of x, and then do another summation on that which yields a integral of a optimization of an integral. Pretty much it takes a really long time.

Do you know of a better way to do the summation that is faster? Hand coding it seemed to go slower.

Currently, I am working with

g(x) = (2/sqrt(3))*pi**(-0.25)*(1 - x**2)*exp(-x**2/2)

but the solution should be general

The paper this comes from is "The Wavelet Transform, Time-Frequency Localization and Signal Analysis" by Daubechies (IEEE 1990)

Thank you

Answer 1

Thanks to all the useful comment, I wrote my own summator that seems to run pretty fast. It anyone has any recommendations to make it better, I will gladly take them.

I will test this on the problem I am working on and once it demonstrates success, I will claim it functional.

def integers(blk_size=100):
    x = arange(0, blk_size)
    while True:
        yield x
        yield -x -1
        x += blk_size

#                                                                                                                                                                                                            
# For convergent summation                                                                                                                                                                                   
# on not necessarily finite sequences                                                                                                                                                                        
# processes in blocks which can be any size                                                                                                                                                                  
# shape that the function can handle                                                                                                                                                                         
#                                                                                                                                                                                                            
def converge_sum(f, x_strm, eps=1e-5, axis=0):
    total = sum(f(x_strm.next()), axis=axis)
    for x_blk in x_strm:
        diff = sum(f(x_blk), axis=axis)
        if abs(linalg.norm(diff)) <= eps:
            # Converged                                                                                                                                                                                      
            return total + diff
        else:
            total += diff

Answer 2

g(x) is almost certainly your bottleneck. A very quick-and-dirty solution would be to vectorize it to operate on an array of integers, then use np.trapz to estimate the integral using the trapezoid rule:

import numpy as np

# appropriate range and step size depends on how accurate you need to be and how
# quickly the sum converges
xmin = -1000000
xmax = 1000000
dx = 1

x = np.arange(xmin, xmax + dx, dx)
gx = (2 / np.sqrt(3)) * np.pi**(-0.25)*(1 - x**2) * np.exp(-x**2 / 2)
sum_gx = np.trapz(gx, x, dx)

Aside from that, you could re-write g(x) using Cython or numba to speed it up.