大型数据集python(scipy)存在interp1d性能问题的线性插值

Question

我有一个大数据集(约300000个数据点(，我从中采样了约300000个数字。我首先形成一个经验CDF，然后使用intrep1d为CDF的逆创建一个插值对象。然后，我从均匀分布中生成随机数，并得到插值函数的值，即采样数：

def sampleFromDistr(data, sampleSize):
    t0 = time.time()
    # forming empirical CDF
    sortedData = np.sort(np.array(data))
    yvals = np.arange(len(sortedData)) / float(len(sortedData))
    # linear interpolation object for the inverse of the cdf
    f = interp1d(yvals, sortedData)
    # create the sample set
    sample = []
    with click.progressbar(range(sampleSize), label='sampling') as bar:
        for i in bar:
            # sampling one by one
            sample.append(float(f(random.uniform(0, max(yvals)))))
    t1 = time.time()
    print t1 - t0
    return sample

问题是，这段代码的工作速度非常慢。它的工作速度似乎因数据而异。

因此，我使用均匀分布的数字作为数据集进行了一些测试：

>>> test = [random.random() for s in range(1000)]
>>> sample = sampleFromDistr(test, 10)
sampling  [####################################]  100%
0.00515699386597
>>> sample = sampleFromDistr(test, 100)
sampling  [####################################]  100%
0.0200171470642
>>> sample = sampleFromDistr(test, 1000)
sampling  [####################################]  100%
0.16183590889
>>> sample = sampleFromDistr(test, 10000)
sampling  [####################################]  100%          
1.56129717827
>>> sample = sampleFromDistr(test, 100000)
sampling  [####################################]  100%          
16.2284870148
>>> sample = sampleFromDistr(test, 1000000)
sampling  [####################################]  100%          
174.504947901

这令人惊讶，因为对于约300000个元素的数据集，估计采样时间约为2小时。所以我尝试增加数据集的大小：

>>> test = [random.random() for s in range(10000)]
>>> sample = sampleFromDistr(test, 1000000)
sampling  [#####-------------------------------]   15%  00:09:42

我还查看了interp1d的源代码。是intrep1d通过调用numpy.searchsorted()来查找最近邻居的部分导致代码变慢了吗？如果是，我该怎么做才能使代码更快？

编辑：我发现bisect.bisect()比numpy.searchsorted()快10倍。是否可以在不修改原始文件的情况下更改原始interp1d方法的这一部分？

第2版：我尝试的解决方案：

import numpy as np
from scipy.interpolate import interp1d
import random
import pdb
import click
import time
import bisect
clip = np.clip

class interp1dMod(interp1d):
    def _call_linear(self, x_new):
        x_new_indices = bisect.bisect_left(self.x, x_new)
        x_new_indices = x_new_indices.clip(1, len(self.x) - 1).astype(int)
        lo = x_new_indices - 1
        hi = x_new_indices
        x_lo = self.x[lo]
        x_hi = self.x[hi]
        y_lo = self._y[lo]
        y_hi = self._y[hi]
        slope = (y_hi - y_lo) / (x_hi - x_lo)[:, None]
        y_new = slope * (x_new - x_lo)[:, None] + y_lo
        return y_new

def sampleFromDistr(data, sampleSize):
    t0 = time.time()
    sortedData = np.sort(np.array(data))
    yvals = np.arange(len(sortedData)) / float(len(sortedData))
    f = interp1dMod(yvals, sortedData)
    sample = []
    with click.progressbar(range(sampleSize), label='sampling') as bar:
        for i in bar:
            sample.append(float(f(random.uniform(0, max(yvals)))))
    t1 = time.time()
    print t1 - t0
    return sample

导致以下错误：CCD_ 7。我做错了什么？

Answer 1

你在这里做了很多奇怪的事情。：-(

您要为每个点再次计算max(yvals)，这意味着您每次都必须循环使用len(sortedData)数字，并且您要使用Python函数来执行此操作；您没有利用矢量化，而是使用缓慢的Python级别循环；甚至你的进度条似乎也会放慢速度。在新代码中，您使用的是bisect.bisect，但它只会返回一个Python整数，因此调用结果x_new_indices似乎很奇怪。

无论如何，如果我们把自己限制在numpy(而不是scipy.stats.rv_continuous的子类(，我会做一些类似的事情

def sampleFromDistr_vectorized(data, sampleSize):
    t0 = time.time()
    # forming empirical CDF
    sortedData = np.sort(np.array(data))
    yvals = np.arange(len(sortedData)) / float(len(sortedData))
    # linear interpolation object for the inverse of the cdf
    f = interp1d(yvals, sortedData)
    # get the random numbers
    r = np.random.uniform(0, yvals.max(), sampleSize)
    # interpolate
    sample = f(r)
    t1 = time.time()
    print(t1 - t0)
    return sample

这给了我

>>> test = np.random.random(10**3)
>>> sample = sampleFromDistr(test, 10**4)
sampling  [####################################]  100%          
1.4801428318023682
>>> sample = sampleFromDistr_onemax_noprogressbar(test, 10**4)
0.26591944694519043
>>> sample = sampleFromDistr_vectorized(test, 10**4)
0.00497126579284668

因此

>>> test = np.random.random(10**6)
>>> sample = sampleFromDistr_vectorized(test, 10**6)
0.3583641052246094

与

>>> sample = sampleFromDistr(test, 10**6)
sampling  [------------------------------------]    0%  12:23:25

(在不到一秒钟的时间里，我会用它运行，但如果这开始占用太多时间，我会使用别名方法，预处理后为O(1(。但这里不值得头疼。(