从Python到Matlab的张量乘积矢量化

Question

我正在将Python中的一些代码转换为Matlab。我有一些代码可以产生相同的结果，但我想知道是否有一种方法可以在Matlab代码中对我的一些for循环进行矢量化，因为它需要很长时间才能运行。在Nxd矩阵中的X，diff是NxNxd张量，kxy是NxN矩阵，gradK是NxNs2张量，并且sumkxy、dxkxy和obj都是Nxd阵。

以下是原始Python代码：

diff = x[:, None, :] - x[None, :, :]  # D_{ij, s}
kxy = np.exp(-np.sum(diff ** 2, axis=-1) / (2 * h ** 2)) / np.power(np.pi * 2.0 * h * h, d / 2) # -1 last dimension K_{ij]
gradK = -diff * kxy[:, :, None] / h ** 2  # N * N * 2
sumkxy = np.sum(kxy, axis=1)
dxkxy = np.sum(gradK, axis=1)  # N * 2    sum_{i} d_i K_{ij, s}
obj = np.sum(gradK / sumkxy[None, :, None], axis=1)  # N * 2

这是我最初的Matlab代码，包含所有for循环：

diff = zeros([n,n,d]);
for i = 1:n
for j = 1:n
for k = 1:d
diff(i,j,k) = x(i,k) - x(j,k);
end
end
end
kxy = exp(-sum(dif.^2, 3)/(2*h^2))/((2*pi*h^2)^(d/2));
sumkxy = sum(kxy,2);
gradK = zeros([n,n,d]);
for i = 1:n
for j = 1:n
for k = 1:d
gradK(i,j,k) = -diff(i,j,k)*kxy(i, j)/h^2;
end
end
end
dxkxy = squeeze(sum(gradK,2));
a = zeros([n,n,d]);
for i =1:n
for j = 1:n
for k = 1:d
a(i,j,k) = gradK(i,j,k)/sumkxy(i);
end
end
end
obj = squeeze(sum(a, 2));

我知道一种更快计算kxy项的方法是使用以下代码：

XY = x*x';
x2= sum(x.^2, 2);
X2e = repmat(x2, 1, n);

H = (X2e + X2e' - 2*XY); % calculate pairwise distance
Kxy = exp(-H/(2*h^2))/((2*pi*h*h)^(d/2));

但后来我很难在没有diff的情况下有效地计算gradK。如有任何帮助或建议，我们将不胜感激！

Answer 1

如果你的目标是计算obj，你甚至不需要计算gradK和a:

sx = sum(x.^2, 2);
H = sx - 2*x*x.' + sx.';
kxy = exp(-H/(2*h^2))/((2*pi*h^2)^(d/2));
kh = kxy / h^2; 
sumkxy = sum(kxy, 2);
khs = kh ./ sumkxy;
obj = khs * x - sum(khs, 2) .* x;

gradK和dif可以这样计算：

dif = reshape(x, n, 1, d) - reshape(x, 1, n, d);
gradK = -dif .* (kxy / h^2);.

Answer 2

我喜欢尝试通过将其分解为"；子部件"；使用一些将快速执行的伪数据，您可以使用这些伪数据来测试代码功能。您可能首先使用的第一个子组件是您的第一个嵌套循环计算diff:

n = 100;
d = 50;
x = round(100*rand(n,d));
tic
diff = zeros([n,n,d]);
for i = 1:n
for j = 1:n
for k = 1:d
diff(i,j,k) = x(i,k) - x(j,k);
end
end
end
toc

首先，考虑最内部的循环：

...
for k = 1:d
diff(i,j,k) = x(i,k) - x(j,k);
end
...

看看这个循环(至少对我来说！(可以大大简化事情。为了将这个"矢量化"；子部件"；我们可以写这样的东西：

diff(i,j,:) = x(i,:) - x(j,:);

既然唾手可得的果实已经过时了，让我们考虑下一层循环。做和工作前一样的把戏吗？

diff(i,:,:) = x(i,:) - x; % where x(:,:) can just be written as x.

如果您不确定，可以通过使用相同(强调相同(的伪数据运行嵌套循环版本和上面的版本来检查这一点，并使用isequal((检查它们是否相等。切入正题，结果应该是一样的，现在你的原始循环是：

tic
diff = zeros([n,n,d]);
for i = 1:n
diff(i,:,:) = x(i,:) - x;
end
toc

对于最后一位，您可以利用matlab的矩阵/数组整形/置换函数。有关更多详细信息，请查阅有关重新整形((或置换((的文档。简言之，如果你将x的一个副本的维数从Nxd重塑或更改为1xNxd，那么从另一个规则大小的矩阵中减去x将在matlab中执行元素运算。例如：

diff = permute(x,[1,3,2]) - permute(x,[3,1,2]); % this is Nx1xd - 1xNxd

应该有效地计算你在第一个循环中寻找的张量差！

如果你愿意的话，我可以扩展这个答案来展示其他循环是如何实现的，但首先用同样的逻辑来尝试其他循环。希望你能保持diff，然后更快地计算kxy。在不知道原始矩阵有多大的情况下，我不能说你应该期望多大的加速。

更新：

我应该加一个，以确保你正在进行元素乘法、除法和转置运算，确保加一个"。"在每个命令之前。例如

gradK(i,j,:) = -diff(i,j,:).*kxy(i, j)/h^2;

有关更多信息，请在Matlab 中查找元素操作