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

Question

我有一个大型数据集(约 300,000 个数据点)，我从中抽取约 300,000 个数字。我首先形成一个经验 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

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

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

我还查看了 interp1d source code .是不是 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

这会导致以下错误:AttributeError: 'int' object has no attribute 'clip'。我做错了什么？

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)。但这里不值得头疼。 )