Scipy leastsq() 函数开销答案

【问题标题】：Scipy leastsq() function overheadScipy leastsq() 函数开销
【发布时间】：2013-01-27 17:39:38
【问题描述】：

我正在开发一个图像分析程序，我已将瓶颈缩小到尝试多次将 2D 高斯拟合到小窗口 (20x20) 像素。 90% 的执行时间都花在了这段代码上。

我正在使用 scipy 食谱中给出的代码来解决这个问题：

   def gaussian(height, center_x, center_y, width_x, width_y):
           """Returns a gaussian function with the given parameters"""
        width_x = float(width_x)
        width_y = float(width_y)
        return lambda x,y: height*exp(
                    -(((center_x-x)/width_x)**2+((center_y-y)/width_y)**2)/2)


def moments(data):
    """Returns (height, x, y, width_x, width_y)
    the gaussian parameters of a 2D distribution by calculating its
    moments """
    total = data.sum()
    X, Y = indices(data.shape)
    x = (X*data).sum()/total
    y = (Y*data).sum()/total
    col = data[:, int(y)]
    width_x = sqrt(abs((arange(col.size)-y)**2*col).sum()/col.sum())
    row = data[int(x), :]
    width_y = sqrt(abs((arange(row.size)-x)**2*row).sum()/row.sum())
    height = data.max()
    return height, x, y, width_x, width_y


def fitgaussian(data):
    """Returns (height, x, y, width_x, width_y)
    the gaussian parameters of a 2D distribution found by a fit"""
    params = moments(data)
    errorfunction = lambda p: ravel(gaussian(*p)(*indices(data.shape)) -
                                 data)
    p, success = optimize.leastsq(errorfunction, params, maxfev=50, ftol=1.49012e-05)
    return p

通过结合 errorfunction() 和 gaussian() 函数，我能够将执行时间减半，因此每次至少最小平方 () 调用 errorfunction() 时都会调用一个函数而不是两个。

这让我相信大部分剩余的执行时间都花在了函数调用开销上，因为 leastsq() 算法调用了 errorfunction()。

有没有办法减少这个函数调用的开销？我不知所措，因为 leastsq() 将函数作为输入。

如果我的描述令人困惑，我提前道歉，我是一名受过培训的机械工程师，并且我正在学习 Python。如果还有其他有用的信息，请告诉我。

【问题讨论】：

如果你适合的真的是高斯，你真的需要打电话给leastsq---你不能只计算矩吗，你已经做了。
我不敢相信我没有检查过这个，谢谢你的建议！快速浏览一下，大多数情况下计算出的矩和最终拟合的高斯相差不到 10%。我将尝试在我的算法中使用它们，看看它的效果如何，然后我会回来查看。
好吧，我回来了。似乎使用时刻已经足够了。谢谢振亚！
Zhenya 可能应该发布他们的建议作为答案，因为它似乎确实回答了这个问题......

标签： python numpy scipy least-squares

【解决方案1】：

由于 exp 是单调的，你可以使用高斯的对数作为你的误差函数，例如

def log_gaussian(height, center_x, center_y, width_x, width_y):
    """Returns a gaussian function with the given parameters"""
    width_x = float(width_x)
    width_y = float(width_y)
    log_height = log(height)
    return lambda x,y: (log_height - 
             (((center_x-x)/width_x)**2 - ((center_y-y)/width_y)**2)/2)

这将导致每次迭代调用 1 次 log，而不是每次迭代调用每个数据集行 1 次 exp。

【讨论】：