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

Is `beta(s)` just a scalar constant? The `s` parameter does not seem to be doing anything. — ali_m, Jul 15 '15 at 20:04
Oh, sorry, I copied half of one equation and the other half of the one above. Let me fix that. g(x) is arbitrary, currently in my code it is the second derivative of the gaussian. — , Jul 15 '15 at 20:24

score 4 · Accepted Answer · edited Jul 15 '15 at 23:16

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

Very nice. I just fixed a couple of errors in your indentation. — ali_m, Jul 15 '15 at 23:16
`np.linalg.norm` is rather slow - you could probably do substantially better by taking the norm yourself, e.g. `np.sqrt(diff.dot(diff))` for 1D — ali_m, Jul 15 '15 at 23:25
Thanks. I decided to stick with linalg.norm because I am using higher dimensional structures (3,4,5 ndarrays). It may be slow, but everything is running remarkably fast. For the record, the code above converges normally, but it not fully robust. 90%-ish percent of the numbers I get are equal to the numbers in the paper, but for certain case, I get different values as though I am missing areas of high density in the approximation — , Jul 17 '15 at 16:47
For fast multidimensional L2 norms you could do `np.sqrt(diff.ravel().dot(diff.ravel()))` — ali_m, Jul 17 '15 at 17:37
So this method converges most of the time. I posted my code on my website, links following. I used it for wavelets. http://acsweb.ucsd.edu/~amacdona/notes/Wavelets/ http://acsweb.ucsd.edu/~amacdona/notes/Wavelets/daubechies-1990-affine-bounds.py — , Jul 22 '15 at 14:43
I left UCSD, so the above link might go stale. More permanent link, https://github.com/aidan-plenert-macdonald/website/tree/master/notes/Wavelets — , Jul 24 '17 at 02:11

score 1 · Answer 2 · answered Jul 15 '15 at 20:53

1

There's a chance Numba improves speed significantly - http://numba.pydata.org

It's slightly painful to install but very easy to use. Have a look at: https://jakevdp.github.io/blog/2015/02/24/optimizing-python-with-numpy-and-numba/

answered Jul 15 '15 at 20:53

alex314159

3,159
2
20
28

ali_m · Answer 3 · 2015-07-15T21:16:46.460

1

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.

edited Jul 15 '15 at 21:16

answered Jul 15 '15 at 21:09

ali_m

71,714
23
223
298

Thanks. I was using a vectorized version of g(x), but the optimization step destroy the vector nature. I am looking into using another optimizer that can do bounded vector optimization (scipy.optimize.differential_evolution). Then I can vectorize everything. Also, I think numerically, it is better to just sum the gx vector. Trapazoidal integration might cut edges too much – Jul 15 '15 at 21:27
Also, I know differential_evolution is slow, but my problem is not necessarily convex, so I need bounded global minimization – Jul 15 '15 at 21:30
If you are interested, I posted an answer with a solution that does the summation for vectorized functions. – Jul 15 '15 at 23:08

Infinite Summation in Python

3 Answers3