成对加权距离向量化

Question

下面这段高效的向量化 Matlab 代码使用一个权重向量 WTS（每个维度一个权重；所有点使用相同的权重）计算两组点 A 和 B 之间的加权欧几里得距离：

    WTS = sqrt(WTS); 

    % modify A and B against weight values
    A = WTS(ones(1,size(A,1)),:).*A;
    B = WTS(ones(1,size(B,1)),:).*B; 

    % calculate distance
    AA = sum(A.*A,2);  
    BB = sum(B.*B,2)'; 
    D = sqrt(AA(:,ones(1,size(B,1))) + BB(ones(1,size(A,1)),:) - 2*A*B'); 

(来源：https://github.com/nolanbconaway/pairdist/blob/master/pairdist.m)

我的问题是：是否有一种有效的向量化形式（Matlab、R或Julia都可以），用于类似的计算，但WTS是与A相同大小的权重向量集合？换句话说，我需要为A中的每个点提供1个权重向量，而不是一个权重向量。

这个答案似乎可以满足我的需求，但它是用Python编写的，我不确定如何将其转换为Matlab/R/Julia：https://dev59.com/JHfZa4cB1Zd3GeqPV8qN#19285289

此外，这不是Efficiently calculating weighted distance in MATLAB的重复，因为那个问题涉及单个权重向量情况，而我明确要求N个权重向量情况。

编辑：示例应用程序：RBF网络和高斯混合模型，其中您可以为每个神经元/组件拥有1个权重向量。对于这些类型的问题，有效的解决方案至关重要。

- rcpinto

---

如果有人想要复制代码，请访问此网址：https://pastebin.com/0bnms3gE - rcpinto

@rahnema1 那个截图已经是第二次运行的结果了。 - rcpinto

不要使用 Float32。使用默认的 Float64 进行测试的结果是 1.73、81.22、23.50 和 2.81 秒。 - rahnema1

@DanGetz，正如我在Distance.jl解决方案下评论的那样，它要求A和B具有相同的大小，W是一个向量，但OP希望A和B具有不同的大小，并且W是一个与A大小相同的矩阵。你能展示一下你的代码吗？ - rahnema1

@DanGetz 我评论回复了你之前提到我的名字的评论。你认为这是一件坏事吗？ - rahnema1

Answer 1

在Julia中，你不必将其向量化以提高效率，只需编写循环即可。实际上，它比这些向量化形式更快，因为它可以融合并消除临时变量。您可以使用这个相当高效的Julia中的成对应用程序的实现进行工作。它具有所有需要的功能，但您可以将其简化。

请注意，向量化并不一定意味着“快”，它只是比在R/Python/MATLAB中循环更快，因为它只是调用一个优化的内核（使用较低级别语言（C/C++）编写），该内核实际上正在循环。但是，组合向量化函数通常会产生许多临时分配，因为每个向量化函数都返回数组。因此，如果您真的需要效率，请避免通常的向量化，并使用允许低成本函数调用/循环的语言进行编写。此帖子详细介绍了高级语言中向量化的问题。

这回答了您其中的一个问题。我没有关于MATLAB或R的好答案。

- Chris Rackauckas

---

我知道Julia具有良好的非向量化性能，也认为它可能是最佳解决方案，但我尝试了一个非加权版本的非向量化解决方案，速度要慢得多。然而，这只是一个天真的2-for实现，也许你提供的链接建议可以给我更好的结果。我会尝试着去改进它们，谢谢 :) - rcpinto

你是否将它放在一个函数中，使用@code_warntype检查其类型稳定性，并消除所有临时分配？ - Chris Rackauckas

我刚刚在问题的评论中添加了我的函数。如果您能看一下就好了。我用@code_warntype运行了它们，但没能找到可以改进的地方（自从上次使用Julia以来已经很长时间了，我还不太记得如何正确使用这些日志）。 - rcpinto

D[i,j] = sum(((X[i, :] - Y[j, :]).*W[i,:]).^2) 这段代码有两个主要问题。首先，每个 X[i,:] 切片都会创建临时数组。其次，这是使用行编写的。在列序语言中始终沿列计算（这也适用于 MATLAB）。这两个问题都非常不好。您应该在其中再加入另一个循环来切割临时数组并更改顺序。 - Chris Rackauckas

Answer 2

这里是MATLAB（R2016b及更高版本）中的矢量化版本：

W2 = 1./W.^2;
D = sqrt(sum((A./W).^2 ,2) - 2 * (A .* W2) * B.' +W2 * (B.^2).');

在R2016b之前的版本中，您可以使用以下代码：

W2 = 1./W.^2;
D = sqrt(bsxfun(@plus,sum((A./W).^2 ,2) , -2 * (A .* W2) * B.' +W2 * (B.^2).'));

MATLAB翻译成Julia：

W2 = 1./W.^2;
z=sqrt.(broadcast(+,sum((A./W).^2 ,2) , -2 * (A .* W2) * B.' .+W2 * (B.^2).'));

在这里，我提出的方法矢量化(Vectorization)与@DanGetz提供的循环(Loop)方法进行比较。其他解决方案在这里不适用。

我们可以看到，对于小于128的维度，循环版本比向量化版本更快。随着维度数量的增加，循环版本的性能会变得更差。

以下代码用于生成图表：

function pdist_vectorized (A::Matrix, B::Matrix, W::Matrix)
    W2 = 1./W.^2;
    return sqrt.(broadcast(+,sum((A./W).^2 ,2) , -2 * (A .* W2) * B.' .+W2 * (B.^2).'));
end

result = zeros(10,2);
for i = 1:10
    A = rand( 3000, 2^i);
    B = rand( 2000, 2^i);
    W = ones(size(A));
    result[i,1]=(@timed pdist_1alloc(A,B,W))[2];
    result[i,2]=(@timed pdist_vectorized(A,B,W))[2];
end

using Plots
pyplot()
plot(2.^(1:10), result, title="Pairwise Weighted Distance",
    label=["Loop" "Vectorization"], lw=3,
    xlabel = "Dimension", ylabel = "Time Elapsed(seconds)")

- rahnema1

---

谢谢，@rahnema1。请看一下我在问题中的新评论，也许你有什么要补充的。 - rcpinto

Answer 3

作为未来读者的附加信息，Distances.jl软件包拥有你能想到的大多数距离的高效实现。一般来说，如果在科学计算中某个操作非常常见，就会有一个良好的工具包来实现它。

using Distances

D = pairwise(WeightedEuclidean(weights), A, B)

- juliohm

---

你的解决方案有两个问题：1. A 和 B 应该是相同大小的矩阵。2. 权重应该是一个向量。 - rahnema1

不需要，A和B的大小可以不同。关于权重，你自己想出来不是很简单吗？我的回答重点是有一个包可以使用，而不是重新发明轮子。其余的只是细节。 - juliohm

矩阵A和B（或X和Y）应该将点的坐标作为列包含在内。权重向量应该包含空间每个维度的一个权重。如果您在R ^ 2中，则这些是2个权重。这是计算成对加权欧几里得距离的最小表示形式。任何偏离此表示形式的内容都可以证明是次优的。 - juliohm

@juliohm，“细节”正是我的问题所在。此外，使用更多的权重并不是次优的，它们是完全不同的问题。在RBF网络或高斯混合模型中只有一个权重向量是没有意义的。 - rcpinto

谢谢@rcpinto，很高兴你解决了问题。未来的用户可能会遇到类似的问题，每个答案都提供了不同的有趣变化。 - juliohm

Answer 4

另一个版本进行了优化，只分配结果矩阵而不分配其他内容：

function pdist_1alloc(A::Matrix, B::Matrix, W::Matrix)
    LA, LD = size(A) ; LB = size(B,1)
    res = zeros(LB, LA)
    indA = 0 ; indB = 0 ; indres = 0
    @inbounds for i=1:LD
        for j=1:LA
            a = A[indA+j] ; w = W[indA+j] ; a2w = a^2*w ; awtmp = -2.0*a*w
            for k=1:LB
                indres += 1
                b = B[indB+k] ; b2w = b^2*w
                res[indres] += a2w+awtmp*b+b2w
            end
        end
        indA += LA ; indB += LB ; indres = 0
    end
    res .= sqrt.(res)
    return res
end

这个版本比 @rahnema1 的版本快约两倍，使用了相同的技巧，但可读性不如前者。此外，我对于一开始误解问题的确切设置（并建议 Distance.jl，这在这里并不直接适用）深感歉意。

- Dan Getz

---

@rahnema1 或许我们对“加权距离”有不同的定义。rahnema1的计算得到这个距离的值为sqrt(2)。我的版本则是4.0（考虑向量A和向量B都为1时，得到单一的值）。根据我的理解，每个维度的平方差被权重乘以，即(3-1)^22+(3-1)^22 = 16，开方后得到4.0。正确答案应该是什么？pdist([1.0 1.0], [3.0 3.0], [2.0 2.0]) - Dan Getz

定义来自于这个python的帖子。我认为你应该转换一下结果。 - rahnema1

@rahnema1 转置是正确的（并且无关紧要），但我们版本给出的值是不同的。对于 pdist([1.0 1.0], [3.0 3.0], [2.0 2.0])，一个是 sqrt(2.0)，另一个是 4.0。 - Dan Getz

我正在处理的定义是 ||(Xi - Yj)/Wi||（或者 ||(Xi - Yj)*Wi||），我不确定是否有被认为是正确的规范数学形式。 - rcpinto

谢谢，丹，我明白了。对我来说，两个定义都可以，因为这些权重最终都会被学习到。顺便说一下，我在JuliaBox上进行了新的测试，使用了A 2000x1000和B 3000x1000的Float32（左）和Float64（右）：https://imgur.com/a/S57YU（第4行是rahnema1的，第5行是你的）。 - rcpinto