Python 中的无限求和

Question

我有一个函数，我需要对（所有整数）数值进行无限求和。求和并不总是需要收敛，因为我可以更改内部参数。该函数看起来像，

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])

使用 Pythonic 伪代码

使用 Scipy 积分方法，我只是将 n 取整并像固定 x 一样积分，

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

这工作得很好，但我必须对 x 作为 x 的函数进行优化，然后对产生积分优化的积分进行另一个求和。几乎需要很长时间。

你知道有更好的求和方法吗？手写代码好像变慢了。

目前，我正在与

合作

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

但解决方案应该是通用的

本文来自 "The Wavelet Transform, Time-Frequency Localization and Signal Analysis" Daubechies (IEEE 1990)

谢谢

Answer 1

Numba 有可能显着提高速度 - http://numba.pydata.org

安装起来有点痛苦，但非常容易使用。看一下： https://jakevdp.github.io/blog/2015/02/24/optimizing-python-with-numpy-and-numba/

Answer 2

g(x) 几乎肯定是您的瓶颈。一个非常快速和肮脏的解决方案是将其矢量化以对整数数组进行操作，然后使用 np.trapz 使用梯形规则估计积分：

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)

除此之外，您可以使用 Cython 或 numba 重写 g(x) 以加快速度。

Answer 3

感谢所有有用的评论，我编写了自己的求和器，看起来运行相当快。如果有人有更好的建议，我会很乐意采纳。

我将针对我正在处理的问题对此进行测试，一旦证明成功，我将声明它可以正常工作。

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

Python 中的无限求和

Infinite Summation in Python

python

numpy

scipy

mathematical-optimization

integral