优化现有的 3D Numpy 矩阵乘法
Optimize this existing 3D Numpy Matrix Multiplication
我有一些代码刚刚完成。它按预期工作。我选择在 Numpy 中使用 dot，因为根据我有限的经验，如果您的系统上安装了 BLAS，它比通常形式的矩阵乘法更快。但是，您会注意到我不得不转置很多东西。我注意到这一点实际上超过了使用 dot.
的好处
我提供论文中找到的数学方程式。注意是递归
这是我的代码实现：
我先提供必要的组件及其尺寸
P = array([[0.73105858, 0.26894142],
[0.26894142, 0.73105858]]) # shape (K,K)
B = array([[6.07061629e-09, 0.00000000e+00],
[0.00000000e+00, 2.57640371e-10]]) # shape (K,K)
dP = array([[[ 0.19661193, -0.19661193],
[ 0. , 0. ]],
[[ 0. , 0. ],
[ 0.19661193, -0.19661193]],
[[ 0. , 0. ],
[ 0. , 0. ]],
[[ 0. , 0. ],
[ 0. , 0. ]]]) # shape (L,K,K)
dB = array([[[ 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00]],
[[ 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00]],
[[-1.16721049e-09, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00]],
[[ 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, -1.27824683e-09]]]) # shape (L,K,K)
a = array([[[ 7.60485178e-08, 5.73923956e-07]],
[[-5.54100398e-09, -8.75213012e-08]],
[[-1.25878643e-08, -1.48361081e-07]],
[[-2.73494035e-08, -1.74585971e-07]]]) # shape (L,1,K)
alpha = array([[0.11594542, 0.88405458]]) # shape (1,K)
c = 1 # scalar
下面是实际的 Numpy 计算。注意所有的转置使用。
term1 = (a/c).dot(P).dot(B)
term2 = dP.transpose((0,2,1)).dot(alpha.T).transpose((0,2,1)).dot(B)
term3 = dB.dot( P.T.dot(alpha.T) ).transpose((0,2,1))
a = term1 + term2 + term3
然后应该得到：
>>> a
array([[[ 1.38388584e-10, -5.87312190e-12]],
[[ 1.05516813e-09, -4.47819530e-11]],
[[-3.76451117e-10, -2.88160549e-17]],
[[-4.06412069e-16, -8.65984406e-10]]])
请注意，alpha 和 a 的形状都是我选择的。如果您发现它可以提供卓越的性能，则可以更改这些设置。
我想指出，我认为现有代码很快。实际上，非常快。但是，我一直想知道是否可以做得更好。请试一试。我已经分析了我的代码（其中有很多 Numpy 广播和矢量化），这不一定是瓶颈，因为在我非常旧的机器上评估需要 23 微秒。但是，它是递归的单个步骤。这意味着它按顺序被评估 N 次。因此，即使是最小的收益也会对大序列产生很大的影响。
感谢您的宝贵时间。
EDIT/UPDATE:
感谢 @max9111，他建议我看一下这个问题 ，我管理了一些 Numba 代码，运行s 比 a 的 Numpy 计算快。与原来的 23 微秒相比，它需要 14 微秒。
这里是：
import numba as nb
@nb.njit(fastmath=True,parallel=True,boundscheck=False)
def get_next_a(a,alpha,P,dP,B,dB,c):
N,M,_ = dP.shape
new_a = np.zeros((N,1,M),dtype=np.float64)
new_a = np.zeros((N,1,M))
entry = 0
for idx in nb.prange(N):
for i in range(M):
for j in range(M):
term1 = a[idx,0,j]*P[j,i]*B[i,i]/c
term2 = alpha[0,j]*dP[idx,j,i]*B[i,i]
term3 = alpha[0,j]*P[j,i]*dB[idx,i,i]
entry += term1 + term2 + term3
new_a[idx,0,i] = entry
entry = 0
return new_a
人们会看到 get_next_a returns 想要的结果。但是，当我在包含 Numpy 的纯 python 函数中调用它时，它会抱怨。这是我的实际代码片段：
def forward_recursions(X,working_params):
# P,dP,B,dB,pi,dpi = setup(X,working_params)
# Dummy Data and Parameters instead of setup
working_params = np.random.uniform(0,2,size=100)
X = np.random.uniform(0,1,size=1000)
P = np.random.uniform(0,1,size=(10,10))
norm = P.sum(axis=1)
P = P/norm[:,None]
dP = np.random.uniform(-1,1,size=(100,10,10))
# We pretend that all 1000 of the 10 x 10 matrices
# only have diagonal entries
B = np.random.uniform(0,1,size=(1000,10,10))
dB = np.random.uniform(0,1,size=(1000,100,10,10))
pi = np.random.uniform(0,1,size=10)
norm = pi.sum()
pi = (pi/norm).reshape(1,10)
dpi = np.random.uniform(0,1,size=(100,1,10))
T = len(X)
N = len(working_params)
M = np.int(np.sqrt(N))
ones = np.ones((M,1))
alpha = pi.dot(B[0])
scale = alpha.dot(ones)
alpha = alpha/scale
ll = np.log(scale)
a = dpi.dot(B[0]) + dB[0].dot(pi.T).transpose((0,2,1))
for t in range(1,T):
old_scale = scale
alpha = alpha.dot(P).dot(B[t])
scale = alpha.dot(ones)
ll += np.log(scale)
alpha = alpha/scale
# HERE IS THE NUMBA FUNCTION
a = get_next_a(a,alpha,P,dP,B[t],dB[t],old_scale)
dll = a.dot(ones).reshape((N,1))/scale
return ll,dll,a
我知道包含我自己的代码取决于其他未包含的函数，因此意味着 forward_recursions 不会 运行。我只是希望它能提供一些观点。
我得到的错误是
TypingError: Invalid use of Function(<built-in function iadd>) with argument(s) of type(s): (Literal[int](0), array(float64, 2d, C))
Known signatures:
* (int64, int64) -> int64
* (int64, uint64) -> int64
* (uint64, int64) -> int64
* (uint64, uint64) -> uint64
* (float32, float32) -> float32
* (float64, float64) -> float64
* (complex64, complex64) -> complex64
* (complex128, complex128) -> complex128
* parameterized
In definition 0:
All templates rejected with literals.
In definition 1:
All templates rejected without literals.
In definition 2:
All templates rejected with literals.
In definition 3:
All templates rejected without literals.
In definition 4:
All templates rejected with literals.
In definition 5:
All templates rejected without literals.
In definition 6:
All templates rejected with literals.
In definition 7:
All templates rejected without literals.
In definition 8:
All templates rejected with literals.
In definition 9:
All templates rejected without literals.
In definition 10:
All templates rejected with literals.
In definition 11:
All templates rejected without literals.
In definition 12:
All templates rejected with literals.
In definition 13:
All templates rejected without literals.
In definition 14:
All templates rejected with literals.
In definition 15:
All templates rejected without literals.
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: typing of intrinsic-call at <ipython-input-251-50e636317ef8> (13)
我不明白这是什么意思。你也许知道我怎么能解决这样的问题。非常感谢您的宝贵时间。
Q : …if one could do better?
您的原样代码在我的 （看起来更老） 机器上执行，不在已发布的 ~ 23 [us] 上，但 ~ 45 [ms] 用于第一次调用和，享受 iCACHE 和 dCACHE 层次结构之间的所有预取 ~77..1xx [us]:
>>> from zmq import Stopwatch; aClk = Stopwatch()
>>> import numpy as np
>>>
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
44679
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
149
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
113
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
128
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
82
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
100
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
77
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
97
>>> a
array([[[ 1.38387304e-10, -5.87323502e-12]],
[[ 1.05516829e-09, -4.47819355e-11]],
[[-3.76450816e-10, -2.60843400e-20]],
[[-1.41384088e-18, -8.65984377e-10]]])
有趣的是，重新运行代码多次，将处理结果重新赋值回a实际上并没有更改 a:
中的值
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
97
>>> a
array([[[ 1.38387304e-10, -5.87323502e-12]],
[[ 1.05516829e-09, -4.47819355e-11]],
[[-3.76450816e-10, -2.60843400e-20]],
[[-1.41384088e-18, -8.65984377e-10]]])
这意味着，代码按原样做了很多工作，以便最终提供 a (一个重新制作的身份，代价是花费 ~ XY [我们] 这样做 - 你是唯一一个决定的人，你的目标应用程序是否可以）
关于希望space改进的评论：
嗯，考虑到 N ~ 1E(3..6) 和 K ~ 10 和 L ~ 100，对任何人都没有太多期望改进工作，赞助重新解决（到目前为止 a 的身份结果）希望提高性能。
寻求改进的目标处理将按顺序重复多于~1,000x …少于~ 1,000,000x，这意味着：
- RAM 绑定问题不是主要问题，因为缓存对静态部分的影响，所有大小都只有几个
[MB]，将以尽可能短的延迟从缓存中重用
- CPU 绑定问题已经在
numpy 工具的设计和工程中预先解决（利用 CPU 的 SIMD 资源和向量对齐的跨步技巧，在可行的情况下 - 所以对用户级编码没有太大期望）
最后，但同样重要的是，有人可能会评论“costs”-of-transposing - numpy 对于必须转置矩阵没有做任何其他事情，只是改变了索引的顺序 - 没有别的。如果这可能会产生一些影响，则可以通过查看 FORTRAN 类型的排序或 C 来预期- 将数据单元的底层存储排序到物理 RAM 中的语言类型，但规模为 1E1 x 1E1 ~ 1E2 x 1E1 x 1E1 最大，这使得这个方面变得可以忽略不计，并且被 in-cache 零回写或处理的性质很好地掩盖了其他与性能相关的影响。
结果：
鉴于上述所有事实和进一步观察，确实获得此处定义的计算的更高吞吐量的最便宜和最合理的步骤是移动到此处的线性工作自由度 - [GHz] CPU-chip，越好（这里性能线性增长）也有尽可能多的 AVX-512 寄存器和尽可能大的 L1i + L1d 缓存（亲和映射策略避免任何其他 O/S 噪声对于 HPC 级性能目标来说都是显而易见的）并且依赖于已经智能的 numpy 工具，针对 CPU-矩阵处理资源（如果需要超越float64 IEEE-754表示，另一个故事开始）。
不要指望用户级代码比这做得更好，numpy-原生 SIMD 对齐处理可以而且将会交付。
对于上述给定的规模，内联装配可以取得优势，但必须付出巨大的人力成本来制作和测试这样一个最终的，但在解决方案的概念上却是一个神秘的变化。如果市场确实需要这样做，请告诉我。
我决定 post 回答我自己的问题。我要感谢@max9111 和@user3666197 的帮助和建议。为了使用 Numba 生成优化的代码，我使用了他们教我的东西。
虽然我在使用 Numba 时遇到了一些问题。人们会注意到我优化的 Numba 版本在 Python 函数中 运行ning 有问题（正如我在 edit/update 部分中指出的那样）。我现在明白为什么了。
Numba 真的很棒，因为它非常快。与 Numpy 相比，Plus one 的计算非常明确。 Numpy 生成更短的代码，对我来说在实际阅读时在数学上更简洁。此外，Numpy 可与原生 Python 无缝协作，因此使用它编写代码非常轻松。就我个人而言，我认为 Numpy 是生成保证 运行 快速工作脚本的最快方法。 Numba 的表现不如 Numpy。事实上，它和 Numpy 相处得不是很好。例如，Numpy 的心爱 dot 不能在 Numba 中使用。
Numba 需要稍微改变一下自己的编码风格。我想认为人们使用 Numba 编码的功能范式：一个函数有输入，它输出输出可能通过其输入流入另一个函数。换句话说，如果您想将数据存放到另一个函数中，则不能将数据生成函数嵌套到将要处理它的函数中。您必须生成数据，然后通过其参数将其提供给处理函数。
然而，嵌套确实有效，但请注意您嵌套的函数也必须是 Numba 编译函数。我的代码不起作用，因为 setup 是一个嵌套的 Python 函数。 Numba 无法识别它，因此 TypingError: Invalid use of Function(<built-in function iadd>)。所以，如果你想嵌套函数，那么你必须确保你所有的函数都是 Numba 函数。这实际上是非常有限制的，并且会带走在 Python 中快速编写脚本的美感。因此，这就是为什么我提到在 Numba 中不嵌套函数，因为这给了你不必创建所有函数的自由 Numba。
既然 Numba 这么快，为什么不希望你的所有功能都在 Numba 中呢？简单的 Numpy 广播和矢量化速度快得吓人 在某些应用程序中比 Numba 更快，例如计算许多数据点和参数的概率密度函数。另外，您仍然想使用 numpy.linalg.solve 等
这是我的优化版本：
@nb.njit(fastmath=True,parallel=True,boundscheck=False)
def forward_recursions_numba(X,
P,
dP,
B,
dB,
pi,
dpi):
# P,dP,B,dB,pi,dpi = setup(X,working_params) # does not work with numba
# print(P)
T = len(X)
N,M,_ = dP.shape
# ones = np.ones((M,1))
# alpha = pi.dot(B[0])
# scale = alpha.dot(ones)
alpha = np.zeros((1,M))
scale = 0
for i in range(M):
alpha[0,i] = pi[0,i]*B[0,i,i]
scale += alpha[0,i]
alpha = alpha/scale
ll = np.log(scale)
# a = dpi.dot(B[0]) + dB[0].dot(pi.T).transpose((0,2,1))
a = np.zeros((N,1,M))
for idx in range(N):
for i in range(M):
a[idx,0,i] = dpi[idx,0,i]*B[0,i,i] + pi[0,i]*dB[0,idx,i,i]
for t in range(1,T):
old_scale = scale
# alpha = alpha.dot(P).dot(B[t])
# scale = alpha.dot(ones)
scale = 0
alpha_new = np.zeros(alpha.shape)
for i in range(M):
entry = 0
for j in range(M):
entry += alpha[0,j]*P[j,i]*B[t,i,i]
alpha_new[0,i] = entry
scale += alpha_new[0,i]
ll += np.log(scale)
alpha = alpha_new/scale
# term1 = (a/old_scale).dot(P).dot(B[t])
# term2 = dP.transpose((0,2,1)).dot(alpha.T).transpose((0,2,1)).dot(B[t])
# term3 = dB[t].dot( P.T.dot(alpha.T) ).transpose((0,2,1))
# a = term1 + term2 + term3
new_a = np.zeros((N,1,M))
for idx in nb.prange(N):
for i in range(M):
entry = 0
for j in range(M):
term1 = a[idx,0,j]*P[j,i]*B[t,i,i]/old_scale
term2 = alpha[0,j]*dP[idx,j,i]*B[t,i,i]
term3 = alpha[0,j]*P[j,i]*dB[t,idx,i,i]
entry += term1 + term2 + term3
new_a[idx,0,i] = entry
a = new_a
# dll = a.dot(ones).reshape((N,1))/scale
dll = np.zeros((N,1))
for idx in nb.prange(N):
dparam = 0
for i in range(M):
dparam += a[idx,0,i]
dll[idx] = dparam/scale
return ll,dll,a
如果我们运行将它放在一些虚拟数据上
X = np.random.uniform(0,1,size=1000)
P = np.random.uniform(0,1,size=(10,10))
norm = P.sum(axis=1)
P = P/norm[:,None]
dP = np.random.uniform(-1,1,size=(100,10,10))
# We pretend that all 1000 of the 10 x 10 matrices
# only have diagonal entries
B = np.random.uniform(0,1,size=(1000,10,10))
dB = np.random.uniform(0,1,size=(1000,100,10,10))
pi = np.random.uniform(0,1,size=10)
norm = pi.sum()
pi = (pi/norm).reshape(1,10)
dpi = np.random.uniform(0,1,size=(100,1,10))
让我们计时：
>>> %timeit forward_recursions_numba(X,P,dP,B,dB,pi,dpi)
51.3 ms ± 389 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
让我们将它与 Numpy 点版本进行比较。
>>> %timeit forward_recursions(X,P,dP,B,dB,pi,dpi)
271 ms ± 1.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
注意： 要获得 Numpy 版本，您只需重新注释 Numpy 代码，注释 Numba 循环并注释掉 @nb.njit 装饰器。
在我的 Lubuntu 机器上使用 lscpu 我们可以看到我的规格是 Intel(R) Core(TM)2 Quad CPU Q6700 @ 2.66GHz 和 6 GB DDR2 RAM (800 MHz)
所以我的代码得到了极大的优化。另外，我学到了很多东西。非常感谢大家的热心帮助和耐心等待。
我有一些代码刚刚完成。它按预期工作。我选择在 Numpy 中使用 dot，因为根据我有限的经验，如果您的系统上安装了 BLAS，它比通常形式的矩阵乘法更快。但是，您会注意到我不得不转置很多东西。我注意到这一点实际上超过了使用 dot.
我提供论文中找到的数学方程式。注意是递归
这是我的代码实现：
我先提供必要的组件及其尺寸
P = array([[0.73105858, 0.26894142],
[0.26894142, 0.73105858]]) # shape (K,K)
B = array([[6.07061629e-09, 0.00000000e+00],
[0.00000000e+00, 2.57640371e-10]]) # shape (K,K)
dP = array([[[ 0.19661193, -0.19661193],
[ 0. , 0. ]],
[[ 0. , 0. ],
[ 0.19661193, -0.19661193]],
[[ 0. , 0. ],
[ 0. , 0. ]],
[[ 0. , 0. ],
[ 0. , 0. ]]]) # shape (L,K,K)
dB = array([[[ 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00]],
[[ 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00]],
[[-1.16721049e-09, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00]],
[[ 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, -1.27824683e-09]]]) # shape (L,K,K)
a = array([[[ 7.60485178e-08, 5.73923956e-07]],
[[-5.54100398e-09, -8.75213012e-08]],
[[-1.25878643e-08, -1.48361081e-07]],
[[-2.73494035e-08, -1.74585971e-07]]]) # shape (L,1,K)
alpha = array([[0.11594542, 0.88405458]]) # shape (1,K)
c = 1 # scalar
下面是实际的 Numpy 计算。注意所有的转置使用。
term1 = (a/c).dot(P).dot(B)
term2 = dP.transpose((0,2,1)).dot(alpha.T).transpose((0,2,1)).dot(B)
term3 = dB.dot( P.T.dot(alpha.T) ).transpose((0,2,1))
a = term1 + term2 + term3
然后应该得到：
>>> a
array([[[ 1.38388584e-10, -5.87312190e-12]],
[[ 1.05516813e-09, -4.47819530e-11]],
[[-3.76451117e-10, -2.88160549e-17]],
[[-4.06412069e-16, -8.65984406e-10]]])
请注意，alpha 和 a 的形状都是我选择的。如果您发现它可以提供卓越的性能，则可以更改这些设置。
我想指出，我认为现有代码很快。实际上，非常快。但是，我一直想知道是否可以做得更好。请试一试。我已经分析了我的代码（其中有很多 Numpy 广播和矢量化），这不一定是瓶颈，因为在我非常旧的机器上评估需要 23 微秒。但是，它是递归的单个步骤。这意味着它按顺序被评估 N 次。因此，即使是最小的收益也会对大序列产生很大的影响。
感谢您的宝贵时间。
EDIT/UPDATE:
感谢 @max9111，他建议我看一下这个问题 a 的 Numpy 计算快。与原来的 23 微秒相比，它需要 14 微秒。
这里是：
import numba as nb
@nb.njit(fastmath=True,parallel=True,boundscheck=False)
def get_next_a(a,alpha,P,dP,B,dB,c):
N,M,_ = dP.shape
new_a = np.zeros((N,1,M),dtype=np.float64)
new_a = np.zeros((N,1,M))
entry = 0
for idx in nb.prange(N):
for i in range(M):
for j in range(M):
term1 = a[idx,0,j]*P[j,i]*B[i,i]/c
term2 = alpha[0,j]*dP[idx,j,i]*B[i,i]
term3 = alpha[0,j]*P[j,i]*dB[idx,i,i]
entry += term1 + term2 + term3
new_a[idx,0,i] = entry
entry = 0
return new_a
人们会看到 get_next_a returns 想要的结果。但是，当我在包含 Numpy 的纯 python 函数中调用它时，它会抱怨。这是我的实际代码片段：
def forward_recursions(X,working_params):
# P,dP,B,dB,pi,dpi = setup(X,working_params)
# Dummy Data and Parameters instead of setup
working_params = np.random.uniform(0,2,size=100)
X = np.random.uniform(0,1,size=1000)
P = np.random.uniform(0,1,size=(10,10))
norm = P.sum(axis=1)
P = P/norm[:,None]
dP = np.random.uniform(-1,1,size=(100,10,10))
# We pretend that all 1000 of the 10 x 10 matrices
# only have diagonal entries
B = np.random.uniform(0,1,size=(1000,10,10))
dB = np.random.uniform(0,1,size=(1000,100,10,10))
pi = np.random.uniform(0,1,size=10)
norm = pi.sum()
pi = (pi/norm).reshape(1,10)
dpi = np.random.uniform(0,1,size=(100,1,10))
T = len(X)
N = len(working_params)
M = np.int(np.sqrt(N))
ones = np.ones((M,1))
alpha = pi.dot(B[0])
scale = alpha.dot(ones)
alpha = alpha/scale
ll = np.log(scale)
a = dpi.dot(B[0]) + dB[0].dot(pi.T).transpose((0,2,1))
for t in range(1,T):
old_scale = scale
alpha = alpha.dot(P).dot(B[t])
scale = alpha.dot(ones)
ll += np.log(scale)
alpha = alpha/scale
# HERE IS THE NUMBA FUNCTION
a = get_next_a(a,alpha,P,dP,B[t],dB[t],old_scale)
dll = a.dot(ones).reshape((N,1))/scale
return ll,dll,a
我知道包含我自己的代码取决于其他未包含的函数，因此意味着 forward_recursions 不会 运行。我只是希望它能提供一些观点。
我得到的错误是
TypingError: Invalid use of Function(<built-in function iadd>) with argument(s) of type(s): (Literal[int](0), array(float64, 2d, C))
Known signatures:
* (int64, int64) -> int64
* (int64, uint64) -> int64
* (uint64, int64) -> int64
* (uint64, uint64) -> uint64
* (float32, float32) -> float32
* (float64, float64) -> float64
* (complex64, complex64) -> complex64
* (complex128, complex128) -> complex128
* parameterized
In definition 0:
All templates rejected with literals.
In definition 1:
All templates rejected without literals.
In definition 2:
All templates rejected with literals.
In definition 3:
All templates rejected without literals.
In definition 4:
All templates rejected with literals.
In definition 5:
All templates rejected without literals.
In definition 6:
All templates rejected with literals.
In definition 7:
All templates rejected without literals.
In definition 8:
All templates rejected with literals.
In definition 9:
All templates rejected without literals.
In definition 10:
All templates rejected with literals.
In definition 11:
All templates rejected without literals.
In definition 12:
All templates rejected with literals.
In definition 13:
All templates rejected without literals.
In definition 14:
All templates rejected with literals.
In definition 15:
All templates rejected without literals.
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: typing of intrinsic-call at <ipython-input-251-50e636317ef8> (13)
我不明白这是什么意思。你也许知道我怎么能解决这样的问题。非常感谢您的宝贵时间。
Q : …if one could do better?
您的原样代码在我的 （看起来更老） 机器上执行，不在已发布的 ~ 23 [us] 上，但 ~ 45 [ms] 用于第一次调用和，享受 iCACHE 和 dCACHE 层次结构之间的所有预取 ~77..1xx [us]:
>>> from zmq import Stopwatch; aClk = Stopwatch()
>>> import numpy as np
>>>
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
44679
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
149
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
113
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
128
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
82
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
100
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
77
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
97
>>> a
array([[[ 1.38387304e-10, -5.87323502e-12]],
[[ 1.05516829e-09, -4.47819355e-11]],
[[-3.76450816e-10, -2.60843400e-20]],
[[-1.41384088e-18, -8.65984377e-10]]])
有趣的是，重新运行代码多次，将处理结果重新赋值回a实际上并没有更改 a:
>>> aClk.start(); a = ( a / c ).dot( P ).dot( B ) + dP.transpose( ( 0, 2, 1) ).dot( alpha.T ).transpose( ( 0, 2, 1 ) ).dot( B ) + dB.dot( P.T.dot( alpha.T ) ).transpose( ( 0, 2, 1 ) ); aClk.stop()
97
>>> a
array([[[ 1.38387304e-10, -5.87323502e-12]],
[[ 1.05516829e-09, -4.47819355e-11]],
[[-3.76450816e-10, -2.60843400e-20]],
[[-1.41384088e-18, -8.65984377e-10]]])
这意味着，代码按原样做了很多工作，以便最终提供 a (一个重新制作的身份，代价是花费 ~ XY [我们] 这样做 - 你是唯一一个决定的人，你的目标应用程序是否可以）
关于希望space改进的评论：
嗯，考虑到 N ~ 1E(3..6) 和 K ~ 10 和 L ~ 100，对任何人都没有太多期望改进工作，赞助重新解决（到目前为止 a 的身份结果）希望提高性能。
寻求改进的目标处理将按顺序重复多于~1,000x …少于~ 1,000,000x，这意味着：
- RAM 绑定问题不是主要问题，因为缓存对静态部分的影响，所有大小都只有几个
[MB]，将以尽可能短的延迟从缓存中重用 - CPU 绑定问题已经在
numpy工具的设计和工程中预先解决（利用 CPU 的 SIMD 资源和向量对齐的跨步技巧，在可行的情况下 - 所以对用户级编码没有太大期望）
最后，但同样重要的是，有人可能会评论“costs”-of-transposing - numpy 对于必须转置矩阵没有做任何其他事情，只是改变了索引的顺序 - 没有别的。如果这可能会产生一些影响，则可以通过查看 FORTRAN 类型的排序或 C 来预期- 将数据单元的底层存储排序到物理 RAM 中的语言类型，但规模为 1E1 x 1E1 ~ 1E2 x 1E1 x 1E1 最大，这使得这个方面变得可以忽略不计，并且被 in-cache 零回写或处理的性质很好地掩盖了其他与性能相关的影响。
结果：
鉴于上述所有事实和进一步观察，确实获得此处定义的计算的更高吞吐量的最便宜和最合理的步骤是移动到此处的线性工作自由度 - [GHz] CPU-chip，越好（这里性能线性增长）也有尽可能多的 AVX-512 寄存器和尽可能大的 L1i + L1d 缓存（亲和映射策略避免任何其他 O/S 噪声对于 HPC 级性能目标来说都是显而易见的）并且依赖于已经智能的 numpy 工具，针对 CPU-矩阵处理资源（如果需要超越float64 IEEE-754表示，另一个故事开始）。
不要指望用户级代码比这做得更好，numpy-原生 SIMD 对齐处理可以而且将会交付。
对于上述给定的规模，内联装配可以取得优势，但必须付出巨大的人力成本来制作和测试这样一个最终的，但在解决方案的概念上却是一个神秘的变化。如果市场确实需要这样做，请告诉我。
我决定 post 回答我自己的问题。我要感谢@max9111 和@user3666197 的帮助和建议。为了使用 Numba 生成优化的代码，我使用了他们教我的东西。
虽然我在使用 Numba 时遇到了一些问题。人们会注意到我优化的 Numba 版本在 Python 函数中 运行ning 有问题（正如我在 edit/update 部分中指出的那样）。我现在明白为什么了。
Numba 真的很棒，因为它非常快。与 Numpy 相比，Plus one 的计算非常明确。 Numpy 生成更短的代码，对我来说在实际阅读时在数学上更简洁。此外，Numpy 可与原生 Python 无缝协作，因此使用它编写代码非常轻松。就我个人而言，我认为 Numpy 是生成保证 运行 快速工作脚本的最快方法。 Numba 的表现不如 Numpy。事实上，它和 Numpy 相处得不是很好。例如，Numpy 的心爱 dot 不能在 Numba 中使用。
Numba 需要稍微改变一下自己的编码风格。我想认为人们使用 Numba 编码的功能范式：一个函数有输入，它输出输出可能通过其输入流入另一个函数。换句话说，如果您想将数据存放到另一个函数中，则不能将数据生成函数嵌套到将要处理它的函数中。您必须生成数据，然后通过其参数将其提供给处理函数。
然而，嵌套确实有效，但请注意您嵌套的函数也必须是 Numba 编译函数。我的代码不起作用，因为 setup 是一个嵌套的 Python 函数。 Numba 无法识别它，因此 TypingError: Invalid use of Function(<built-in function iadd>)。所以，如果你想嵌套函数，那么你必须确保你所有的函数都是 Numba 函数。这实际上是非常有限制的，并且会带走在 Python 中快速编写脚本的美感。因此，这就是为什么我提到在 Numba 中不嵌套函数，因为这给了你不必创建所有函数的自由 Numba。
既然 Numba 这么快，为什么不希望你的所有功能都在 Numba 中呢？简单的 Numpy 广播和矢量化速度快得吓人 在某些应用程序中比 Numba 更快，例如计算许多数据点和参数的概率密度函数。另外，您仍然想使用 numpy.linalg.solve 等
这是我的优化版本：
@nb.njit(fastmath=True,parallel=True,boundscheck=False)
def forward_recursions_numba(X,
P,
dP,
B,
dB,
pi,
dpi):
# P,dP,B,dB,pi,dpi = setup(X,working_params) # does not work with numba
# print(P)
T = len(X)
N,M,_ = dP.shape
# ones = np.ones((M,1))
# alpha = pi.dot(B[0])
# scale = alpha.dot(ones)
alpha = np.zeros((1,M))
scale = 0
for i in range(M):
alpha[0,i] = pi[0,i]*B[0,i,i]
scale += alpha[0,i]
alpha = alpha/scale
ll = np.log(scale)
# a = dpi.dot(B[0]) + dB[0].dot(pi.T).transpose((0,2,1))
a = np.zeros((N,1,M))
for idx in range(N):
for i in range(M):
a[idx,0,i] = dpi[idx,0,i]*B[0,i,i] + pi[0,i]*dB[0,idx,i,i]
for t in range(1,T):
old_scale = scale
# alpha = alpha.dot(P).dot(B[t])
# scale = alpha.dot(ones)
scale = 0
alpha_new = np.zeros(alpha.shape)
for i in range(M):
entry = 0
for j in range(M):
entry += alpha[0,j]*P[j,i]*B[t,i,i]
alpha_new[0,i] = entry
scale += alpha_new[0,i]
ll += np.log(scale)
alpha = alpha_new/scale
# term1 = (a/old_scale).dot(P).dot(B[t])
# term2 = dP.transpose((0,2,1)).dot(alpha.T).transpose((0,2,1)).dot(B[t])
# term3 = dB[t].dot( P.T.dot(alpha.T) ).transpose((0,2,1))
# a = term1 + term2 + term3
new_a = np.zeros((N,1,M))
for idx in nb.prange(N):
for i in range(M):
entry = 0
for j in range(M):
term1 = a[idx,0,j]*P[j,i]*B[t,i,i]/old_scale
term2 = alpha[0,j]*dP[idx,j,i]*B[t,i,i]
term3 = alpha[0,j]*P[j,i]*dB[t,idx,i,i]
entry += term1 + term2 + term3
new_a[idx,0,i] = entry
a = new_a
# dll = a.dot(ones).reshape((N,1))/scale
dll = np.zeros((N,1))
for idx in nb.prange(N):
dparam = 0
for i in range(M):
dparam += a[idx,0,i]
dll[idx] = dparam/scale
return ll,dll,a
如果我们运行将它放在一些虚拟数据上
X = np.random.uniform(0,1,size=1000)
P = np.random.uniform(0,1,size=(10,10))
norm = P.sum(axis=1)
P = P/norm[:,None]
dP = np.random.uniform(-1,1,size=(100,10,10))
# We pretend that all 1000 of the 10 x 10 matrices
# only have diagonal entries
B = np.random.uniform(0,1,size=(1000,10,10))
dB = np.random.uniform(0,1,size=(1000,100,10,10))
pi = np.random.uniform(0,1,size=10)
norm = pi.sum()
pi = (pi/norm).reshape(1,10)
dpi = np.random.uniform(0,1,size=(100,1,10))
让我们计时：
>>> %timeit forward_recursions_numba(X,P,dP,B,dB,pi,dpi)
51.3 ms ± 389 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
让我们将它与 Numpy 点版本进行比较。
>>> %timeit forward_recursions(X,P,dP,B,dB,pi,dpi)
271 ms ± 1.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
注意： 要获得 Numpy 版本，您只需重新注释 Numpy 代码，注释 Numba 循环并注释掉 @nb.njit 装饰器。
在我的 Lubuntu 机器上使用 lscpu 我们可以看到我的规格是 Intel(R) Core(TM)2 Quad CPU Q6700 @ 2.66GHz 和 6 GB DDR2 RAM (800 MHz)
所以我的代码得到了极大的优化。另外，我学到了很多东西。非常感谢大家的热心帮助和耐心等待。