加快矩阵计算(在子阵列上循环)[numpy]
问题描述:
我正在尝试为我的学生作业编写一个算法,它运行良好。但是,计算需要很长时间,尤其是对于大数组。 这部分代码会减慢所有程序。加快矩阵计算(在子阵列上循环)[numpy]
Shapes: X.shape = mask.shape = logBN.shape = (500,500,1000),
F.shape = (20,20),
A.shape = (481,481),
s2 -- scalar.
我应该如何改变这种代码,使其更快?
h = F.shape[0]
w = F.shape[1]
q = np.zeros((A.shape[0], A.shape[1], X.shape[2]))
for i in range(A.shape[0]):
for j in range(A.shape[1]):
mask[:,:,:] = 0
mask[i:i + h,j:j + w,:] = 1
q[i,j,:] = ((logBN*(1 - mask)).sum(axis=(0,1)) +
(np.log(norm._pdf((X[i:i + h,j:j + w,:]-F[:,:,np.newaxis])/s2)/s2)).sum(axis=(0,1))
答
只是试图让你的内环感
mask[:,:,:] = 0
mask[i:i + h,j:j + w,:] = 1
q[i,j,:] = ((logBN*(1 - mask)).sum(axis=(0,1)) +
(np.log(norm._pdf((X[i:i + h,j:j + w,:]-F[:,:,np.newaxis])/s2)/s2)).sum(axis=(0,1))
看起来像
idx = (slice(i,i+h), slice(j,j_w), slice(None))
mask = np.zeros(X.shape)
mask(idx) = 1
mask = 1 - mask
# alt mask=np.ones(X.shape);mask[idx]=0
term1 = (logBN*mask).sum(axis=(0,1))
term2 = np.log(norm._pdf((X[idx] - F[...,None])/s2)/s2).sum(axis=(0,1))
q[i,j,:] = term1 + term2
所以idx
和mask
在A
定义一个子阵列。您正在阵列外部使用logBN
;里面还有term
。您正在对第1个2个dim进行求和,因此term1
和term2
的形状为X.shape[2]
,您可以将其保存在q
中。
该面具/窗口是20x20。
作为第一次切割,我试图一次计算所有i
,j
的term2
。这看起来像一个典型的滑动窗口问题。我也尝试将term1
表示为减法 - 整个logBN
减去此窗口。
答
后通过log
,exp
,power
代数运算重杂耍,这一切来到这个 -
# Params
m,n = F.shape[:2]
k1 = 1.0/(s2*np.sqrt(2*np.pi))
k2 = -0.5/s2**2
k3 = np.log(k1)*m*n
out = np.zeros((A.shape[0], A.shape[1], X.shape[2]))
for i in range(A.shape[0]):
for j in range(A.shape[1]):
mask[:] = 1
mask[i:i + h,j:j + w,:] = 0
XF = (X[i:i + h,j:j + w,:]-F[:,:,np.newaxis])
p1 = np.einsum('ijk,ijk->k',logBN,mask)
p2 = k2*np.einsum('ijk,ijk->k',XF,XF)
out[i,j,:] = p1 + p2
out += k3
使用几件事情是 -
1] norm._pdf基本上是:norm.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
。所以,我们可以内嵌实现并在脚本级别对其进行优化。
2]由标量划分将不会有效,所以这些被它们的倒数乘法代替。所以,在进入循环之前,作为一个预处理存储他们的倒数。
这是不完整的 - 将所有变量(F,A,X)放在一起,以便人们可以使用某些东西。如果在数组上迭代,通常最好转换为python列表,因为它非常慢 - 使用向量操作最快。 – kabanus
@kabanus我不能在程序工作期间生成它们。 – Acapello
我建议打印一次,并在开始时粘贴结果。 – kabanus