有一个矩阵（numpy）等价于这些递归方程吗？

Question

我想用Python实现以下递归方程
（背景是金融市场）：

这很容易通过循环按顺序完成：

def compute(x, y, χ, γ, a0, b0):
    a = [a0]
    b = [b0]

    α = [a[0]]
    β = [b[0]]

    for i in range(len(x)):
        a.append(-α[i]*x[i] + β[i]*y[i]*γ[i])
        b.append(-β[i]*y[i] + α[i]*x[i]*χ[i])

        α.append(α[-1] + a[-1])
        β.append(β[-1] + b[-1])
    return α, β

一个小型阵列示例：

x = np.array([0.6, 0.4, 0., 0., 0.9])
y = np.array([0., 0., 0.3, 0.9, 0.])

χ = np.arange(100., 105., 1.)
γ = 1. / (χ - 1.)

print(np.array(compute(x, y, χ, γ, 1., 0.)))

会产生

>>> [[ 1.          0.4         0.24        0.46621782  0.93661782  0.09366178]  
     [ 0.         60.         76.16       53.312       5.3312     92.99862812]]

有没有办法在NumPy中实现这一点（我希望它能显著加快计算速度）

换句话说：要计算整个a和b向量而不使用循环，只需NumPy函数

Answer 1

如果您重新排列元素以获得以下形式，则可以使用矩阵乘法

之后，你可以通过简单的矩阵乘法来计算一切

请注意，2x2矩阵的k、j、l、m都可用，它们构造的矩阵可以预先计算

在这种情况下，它们将是：

k = 1-x  
l = y*γ  
m = x*χ  
n = 1-y 

此外，我建议预先分配任何可能使用的数组，因为大小始终可用（附加列表的成本非常高）

无论如何，for循环是不可避免的。但我想下面的内容会让它更整洁

from functools import reduce
import numpy as np
def compute(mat, inp):
    return reduce(np.dot, mat) @ inp

Answer 2

首先要明确的是：我不知道如何使用numpy使其更快。所以这不是你问题的答案

但是：您可以使用numba实现一些加速：

from numba import jit
import numpy as np

N = 1000000

x = np.random.rand(N)
y = np.random.rand(N)
gamma = np.random.rand(N)
chi = np.random.rand(N)

a = [.2]
b = [.3]
alpha = [a[0]]
beta = [b[0]]

def compute_in_python():
    for i in range(len(x)):
        a.append(-alpha[i]*x[i] + beta[i]*y[i]*gamma[i])
        b.append(-beta[i]*y[i] + alpha[i]*x[i]*chi[i])

        alpha.append(alpha[-1]+a[-1])
        beta.append(beta[-1]+b[-1])


@jit(nopython=True)
def compute_with_numba(x,y,gamma,chi,a0, b0, alpha0, beta0):
    N = len(x)
    a = np.empty(N+1)
    b = np.empty(N+1)
    alpha = np.empty(N+1)
    beta = np.empty(N+1)
    a[0] = a0
    b[0] = b0
    alpha[0] = alpha0
    beta[0] = beta0
    for i in range(N):
        a[i+1] = -alpha[i] * x[i] + beta[i] * y[i] * gamma[i]
        b[i+1] = -beta[i] * y[i] + alpha[i] * x[i] * chi[i]

        alpha[i+1] = alpha[i]+a[i+1]
        beta[i+1] = beta[i]+b[i+1]

    return a,b,alpha,beta

普通python循环：

In [23]: %time compute_in_python()
CPU times: user 1.6 s, sys: 24.8 ms, total: 1.63 s
Wall time: 1.63 s

麻木紧张：

In [42]: %time res = compute_with_numba(x,y,gamma,chi,a[0], b[0], alpha[0], beta[0])
CPU times: user 13 ms, sys: 3.36 ms, total: 16.4 ms
Wall time: 16.4 ms

请注意，对compute_with_numba的第一个调用将触发jit，因此您应该测量第二个调用的运行时

同样，这不是你问题的答案，但它仍然快了大约100倍