随机过程的循环优化

Question

因为我调用这个代码超过1000次，有没有一种方法可以优化这个代码，这需要1.73秒？你知道吗

def generate():
    S0    = 0
    T     = 1.
    nt    = 100000
    lbd   = 500.
    mu    = 0
    sigma = 1.

    dt = T/nt
    St  = [S0]  * nt
    sqrtdt = np.sqrt(dt)
    dBt = np.random.normal(0, sqrtdt, nt)

    for k in xrange(1, nt):
        dSt = lbd * (mu - St[k-1]) * dt + sigma * dBt[k]
        St[k] = St[k-1] + dSt
    return St

Answer 1

您可以从for-loop中挤出更多的工作，但是同时生成所有路径（假设您有足够的内存）：

import numpy as np

def generate_orig(T=1., nt=100000, lbd=500., mu=0, sigma=1., S0=0):
    dt = T/nt
    St  = [S0]  * nt
    sqrtdt = np.sqrt(dt)
    dBt = np.random.normal(0, sqrtdt, nt)

    for k in xrange(1, nt):
        dSt = lbd * (mu - St[k-1]) * dt + sigma * dBt[k]
        St[k] = St[k-1] + dSt
    return St

def generate(T=1., nt=100000, lbd=500., mu=0, sigma=1., S0=0, npaths=1):
    dt = T/nt
    St  = np.full((nt, npaths), S0)
    sqrtdt = np.sqrt(dt)
    dBt = np.random.normal(0, sqrtdt, size=(nt, npaths))

    for k in xrange(1, nt):
        dSt = lbd * (mu - St[k-1]) * dt + sigma * dBt[k]
        St[k] = St[k-1] + dSt
    return St

下面是100条路径上的timeit基准测试。你知道吗

In [55]: %timeit [generate_orig() for i in xrange(100)]
1 loops, best of 3: 23.6 s per loop

In [56]: %timeit generate(npaths=100)
1 loops, best of 3: 1.97 s per loop

您还可以通过使用Cython来改进for-loop的性能。你知道吗

Answer 2

我想为unutbu's answer提供两种备选方案。他写的是对的™ 如果您不想依赖于Cython或JIT编译器，那么批量生成输出St是合适的。你知道吗

---

Numba

我从他的答案中抓取generate_orig()，并将Python列表St转换为numpy数组：

import numpy as np

def generate_orig(T=1., nt=100000, lbd=500., mu=0, sigma=1., S0=0):
    dt = T/nt
    St  = np.full(nt, fill_value=S0, dtype=np.float64)
    sqrtdt = np.sqrt(dt)
    dBt = np.random.normal(0, sqrtdt, nt)

    for k in xrange(1, nt):
        dSt = lbd * (mu - St[k-1]) * dt + sigma * dBt[k]
        St[k] = St[k-1] + dSt
    return St

时间安排：

%timeit [generate_orig() for i in xrange(100)]
1 loops, best of 3: 25.4 s per loop

到目前为止没有改进，和以前一样。但是，对于Numba，只需添加@autojit：

import numpy as np
from numba import autojit

@autojit
def generate_orig(T=1., nt=100000, lbd=500., mu=0, sigma=1., S0=0):
    # The rest is exactly the same as before

计时下降：

%timeit [generate_orig(1., 100000, 500., 0, 1., 0) for i in xrange(100)]
1 loops, best of 3: 642 ms per loop

我觉得太棒了！只需添加@autojit即可加速40倍！你知道吗

---

赛顿

下面是带有typed memoryviews的Cython版本：

%%cython
# cython: infer_types=True
# cython: boundscheck=False
# cython: wraparound=False
import numpy as np
cimport numpy as np

def generate_cython(double T=1., int nt=100000, double lbd=500., double mu=0, double sigma=1., double S0=0):
    cdef int k
    cdef double dt, dSt
    cdef double[:] vSt, vdBt
    dt = T/nt
    St  = np.full(nt, fill_value=S0, dtype=np.float64)
    vSt = St
    vdBt = np.random.normal(0.0, np.sqrt(dt), nt)
    for k in xrange(1, nt):
        dSt = lbd * (mu - vSt[k-1]) * dt + sigma * vdBt[k]
        vSt[k] = vSt[k-1] + dSt
    return St

时间安排：

%timeit [generate_cython(1., 100000, 500., 0, 1., 0) for i in xrange(100)]
1 loops, best of 3: 638 ms per loop

代码和Numba版本一样快（微小的差别只是噪音）。然而，代码变得丑陋，所有这些类型声明使它笨拙。：（好吧，不是灾难，但仍然是。你知道吗

---

与unutbu的答案相比，这两种解决方案的速度都提高了3倍，而unutbu的答案在我的机器上运行时间为1.97秒。然而，正如我在一开始所说的，如果你不想依赖Cython或Nuba，他的解决方案是可行的。（两者都有缺点；如果有人想避免这种依赖，这是可以理解的。）

如果我们对unutbu的解应用Numba或Cython，会发生什么？这会导致更快的代码吗？不，和麻木没有区别。赛顿让事情变得更糟。嗯，也许一个赛昂大师能想出更好的解决办法。。。你知道吗