Curve_fit 到 apply_along_axis。如何加快速度？答案

【问题标题】：Curve_fit to apply_along_axis. How to speed it up?Curve_fit 到 apply_along_axis。如何加快速度？
【发布时间】：2016-11-05 15:21:14
【问题描述】：

我有一些大型数据集，我想在其中拟合单指数时间衰减。

数据由在不同时间获取的多个 4D 数据集组成，因此拟合应该沿着第 5 维（通过数据集）运行。

我目前使用的代码如下：

import numpy as np
import scipy.optimize as opt

[... load 4D datasets ....]
data = (dataset1, dataset2, dataset3)
times = (10, 20, 30)

def monoexponential(t, M0, t_const):
    return M0*np.exp(-t/t_const)

# Starting guesses to initiate  descent.
M0_init = 80.0
t_const_init = 50.0
init_guess = (M0_init, t_const_init)

def fit(vector):
    try:
        nlfit, nlpcov = opt.curve_fit(monoexponential, times, vector,
                                      p0=init_guess,
                                      sigma=None,
                                      check_finite=False,
                                      maxfev=100, ftol=0.5, xtol=1,
                                      bounds=([0, 2000], [0, 800]))
        M0, t_const = nlfit
    except:
        t_const = 0

    return t_const

# Concatenate datasets in data into a single 5D array.
concat5D = np.concatenate([block[..., np.newaxis] for block in data],
                     axis=len(data[0].shape))

# And apply the curve fitting along the last dimension.
decay_map = np.apply_along_axis(fit, len(concat5D.shape) - 1, concat5D)

代码运行良好，但需要很长时间（例如，dataset1.shape == (100,100,50,500)）。我读过一些其他主题提到apply_along_axis 非常慢，所以我猜这是罪魁祸首。不幸的是，我真的不知道这里可以使用什么替代方案（除了可能是显式的 for 循环？）。

有没有人知道我可以做些什么来避免apply_along_axis 并加快curve_fit 被多次调用的速度？

【问题讨论】：

标签： python performance python-2.7 numpy scipy

【解决方案1】：

所以您将 fit 操作 100*100*50*500 次应用于一维数组（示例中为 3 个值，现实生活中更多？）？

apply_along_axis 会遍历输入数组的所有维度，除了一个。 fit 不会同时在多个轴上编译或执行此操作。

如果没有apply_along_axis，最简单的方法是将数组重塑为二维数组，将 (100,100,50,500) 压缩为一维 (250...,)，然后对其进行迭代。然后重塑结果。

我在想，在最后一个轴上连接 datasets 可能比在第一个轴上连接要慢，但时间表明并非如此。

np.stack 是concatenate 的新版本，可以轻松地在任何位置添加新轴。

In [319]: x=np.ones((2,3,4,5),int)
In [320]: d=[x,x,x,x,x,x]

In [321]: np.stack(d,axis=0).shape   # same as np.array(d)
Out[321]: (6, 2, 3, 4, 5)

In [322]: np.stack(d,axis=-1).shape
Out[322]: (2, 3, 4, 5, 6)

对于更大的列表（带有一个微不足道的sum 函数）：

In [295]: d1=[x]*1000       # make a big list

In [296]: timeit np.apply_along_axis(sum,-1,np.stack(d1,-1)).shape
10 loops, best of 3: 39.7 ms per loop

In [297]: timeit np.apply_along_axis(sum,0,np.stack(d1,0)).shape
10 loops, best of 3: 39.2 ms per loop

一个使用数组重塑时间大致相同的显式循环

In [312]: %%timeit 
   .....: d2=np.stack(d1,-1)
   .....: d2=d2.reshape(-1,1000)
   .....: res=np.stack([sum(i) for i in d2],0).reshape(d1[0].shape)
   .....: 
10 loops, best of 3: 39.1 ms per loop

但是像sum 这样的函数可以在整个数组上工作，而且速度要快得多

In [315]: timeit np.stack(d1,-1).sum(-1).shape
100 loops, best of 3: 3.52 ms per loop

因此，更改堆叠和迭代方法不会对速度产生太大影响。但是改变“拟合”使其可以在多个维度上工作可能会有很大帮助。我对optimize.fit 了解不多，不知道这是否可能。

=====================

我刚刚研究了apply_along_axis 的代码。它基本上构造了一个看起来像ind=(0,1,slice(None),2,1) 的索引，然后执行func(arr[ind])，然后递增它，就像带进位的长算术一样。所以它只是系统地遍历所有元素，同时保持一个轴为: 切片。

【讨论】：

感谢您的回答。我查看了 scipy 的 curve_fit 和 least_squares 代码。我发现有数以千计的（在我的情况下是不必要的）检查、函数调用和辅助变量赋值，在多次调用函数时可能会造成巨大的开销。我想这里最好的办法是编写自己的函数来直接调用求解器trf_no_bounds。

【解决方案2】：

在这种特殊情况下，您要拟合单个指数，最好记录您的数据。然后拟合变成线性的，比非线性最小二乘要快得多，并且很可能会被矢量化，因为它几乎变成了线性代数问题。

（当然，如果你知道如何改进 least_squares，scipy 开发人员可能会感激你。）

【讨论】：