成对加权距离矢量化

Question

以下高效和矢量化的 Matlab 代码使用权重向量 WTS（每个维度 1 个权重；所有点的权重相同）计算两组点 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 个权重向量，而不是 1 个权重向量。

这个答案似乎满足了我的需要，但它在 Python 中，我不确定如何将其转换为 Matlab/R/Julia：

此外，不是的副本，因为该问题涉及单个权重向量的情况，而我明确要求 N 个权重向量的情况。

编辑：应用示例：RBF 网络和高斯混合模型，其中每个 neuron/component 您（可以）有 1 个权重向量。有效解决问题对于此类问题至关重要。

Answer 1

在 Julia 中，您不必对其进行矢量化以提高效率，只需编写循环，它无论如何都会比这些矢量化形式更快，因为它可以融合并消除临时变量。 Here's a pretty efficient implementation of pairwise applies in Julia 您可以从中工作。它有所有的花里胡哨，但如果你愿意，你可以把它配对。

请注意，矢量化不一定是 "fast"，它只是比 R/Python/MATLAB 中的循环更快，因为它只是对用低级语言编写的优化内核进行单个函数调用 (C/C++) 实际上是在循环。但是将向量化函数放在一起通常会有很多临时分配，因为每个向量化函数 return 数组。因此，如果你真的需要效率，你应该避免一般的向量化，并用允许低成本函数调用/循环的语言编写它。 This post explains more about issues with vectorization in high level languages.

这回答了您提出的三个问题之一。我对 MATLAB 或 R 没有很好的答案。

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)")

Answer 3

作为为未来读者提供的附加信息，包 Distances.jl 具有您能想到的大多数距离的有效实现。作为一般建议，如果一个操作在科学计算中很常见，那么会有一个包很好地实现它。

using Distances

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

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 的版本快大约 2 倍，并且使用相同的技巧，但可读性差。此外，对于首先误解问题的确切设置（并建议 Distance.jl 此处不直接适用），我深表歉意。

成对加权距离矢量化

Pairwise weighted distance vectorization

matlab

r

vectorization

euclidean-distance

julia