如何知道一个算法是否可以实现

Question

我在做一个实验来计算信号的近似熵。详细信息（和实际代码）可以在它的Wikipedia page上找到。不幸的是，虽然算法本身可以工作，但对于大型数据集来说速度非常慢（例如，对于一个2000长的信号，大约需要25秒）。由于我需要对更长的信号进行计算，以这种速度，我希望我的实验至少持续1个月。我想知道有没有办法加快算法的速度。你知道吗

import numpy as np

def ApEn(U, m, r):

    def _maxdist(x_i, x_j):
        return max([abs(ua - va) for ua, va in zip(x_i, x_j)])

    def _phi(m):
        x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)]
        C = [len([1 for x_j in x if _maxdist(x_i, x_j) <= r]) / (N - m + 1.0) for x_i in x]
        return (N - m + 1.0)**(-1) * sum(np.log(C))

    N = len(U)

    return abs(_phi(m + 1) - _phi(m))

Answer 1

通常在优化时，你应该从减少算法复杂度的算法优化开始，而不是简单的常数。你知道吗

一个经验法则是查看最里面的循环—它包含执行次数最多的操作。你知道吗

我不确定我是否正确地阅读了代码，但它看起来像U是一个矩阵，_maxdist对它的列进行计算。在这种情况下，确保每列只执行一次计算是有意义的。你知道吗

例如，为每一列计算它的值，存储在数组中并在_phi中使用它。你知道吗

Answer 2

我没有看全部内容，但给你一个例子，说明如何使用向量计算优化函数：

def maxdist_opti(x_i,x_j):
    return max(abs(x_i-x_j))

当您的数据存储到numpy数组中时，您可以在它们上使用numpy操作符（有很多操作符，您可以在这里查看：https://docs.scipy.org/doc/numpy-1.13.0/user/index.html），这样会更快，在上面的例子中，我在numpy数组上使用了soustraction和np.max函数。你知道吗

这里，使用随机数据：

x_i = np.random.rand(10000)
x_j = np.random.rand(10000)

这里使用的数据并不长，但是您可以看到非常好的性能提升：

%timeit _maxdist(x_i,x_j)
100 loops, best of 3: 3.01 ms per loop

%timeit maxdist_opti(x_i,x_j)
10000 loops, best of 3: 28 µs per loop

您可以使用下面的逻辑只对整个公式进行向量计算，性能上的增益将是巨大的。你知道吗

请注意，数据越长，使用向量计算的优化程度就越高。你知道吗

Answer 3

如果您愿意将该函数移到cython并添加一些类型注释，那么将获得显著的性能提升。这是我的算法版本：

你知道吗阿彭.pyx地址：

cimport cython
from libc.math cimport fabs, log
import numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.initializedcheck(False)
@cython.cdivision(True)
cdef double max_dist(double[:] x_i, double[:] x_j, int m) nogil:
    #Performs the max function described in step 4 of ApEn algorithm
    cdef double out
    cdef double dist
    out = fabs(x_i[0] - x_j[0])
    for k in range(1, m - 1):
        dist = fabs(x_i[k] - x_j[k])
        if dist > out:
            out = dist
    return out

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.initializedcheck(False)
@cython.cdivision(True)
cdef double phi(double[:] Sn, int m, int r):
    cdef int N = len(Sn)
    cdef int i
    cdef int j
    cdef int k

    cdef int c_val
    cdef int counter
    cdef double phi_sum = 0
    cdef double phi
    cdef double m_dist

    #Performs step 3 of the ApEn algorithm
    cdef double[:, :] x = np.empty((N - m  + 1, m), dtype=np.float64)
    with nogil:
        for i in range(N - m + 1):
            for j in range(0, m):
                x[i, j] = Sn[j + i]

        #Performs a combined steps 4 & 5 of the ApEn algorithm
        for i in range(N - m + 1):
            counter = 0
            for j in range(N - m + 1):
                m_dist = max_dist(x[i], x[j], m)
                c_val = 1 if m_dist <= r else 0
                counter += c_val
            phi_sum += log(counter / (N - m + 1.0))
        phi = phi_sum / (N - m + 1.0)
    return phi

cpdef double approx_entropy(double[:] Sn, int m, int r):#Passing in steps 1 & 2 of the ApEn algorithm
    cdef double ApEn = abs(phi(Sn, m, r) - phi(Sn, m + 1, r))#Performs step 6 of the ApEn algorithm
    return ApEn

你知道吗apen.pxd公司地址：

cdef double max_dist(double[:] x_i, double[:] x_j, int m) nogil
cdef double phi(double[:] Sn, int m, int r)
cpdef double approx_entropy(double[:] Sn, int m, int r)

你知道吗设置.pxd地址：

from distutils.core import setup
from Cython.Build import cythonize
from distutils.core import Extension
import numpy as np

extensions = [
    Extension("apen", sources=["apen.pyx"], include_dirs=[np.get_include()], extra_compile_args=["-w"]),
]

setup(
    ext_modules = cythonize(extensions)
)

你知道吗主.py地址：

import time
import apen
import numpy as np

start = time.time()
data = np.random.rand(2000)
#data = np.array([85, 80, 89] * 17, dtype=np.float64)
answer = apen.approx_entropy(Sn=data, m=2, r=3)
print(answer)
end = time.time()
print(end - start)

用这个代码在我的笔记本电脑上计算2000个随机数据点，cython代码用0.36秒计算ApEn。相比之下，wikipedia代码用14.75秒。这相当于速度提升了40倍。希望这对你有帮助！你知道吗