Search code examples
pythonnumpyscipymathematical-optimizationintegral

Infinite Summation in Python


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


Solution

  • 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