在Matlab中有效地计算成对平方欧几里德距离

Question

给出两组 - d维点.如何在Matlab中最有效地计算成对平方欧氏距离矩阵？

符号: 集合1由(numA,d)-matrix 给出A,集合2由(numB,d)-matrix 给出B.得到的距离矩阵应为格式(numA,numB).

示例点:

d = 4;            % dimension
numA = 100;       % number of set 1 points
numB = 200;       % number of set 2 points
A = rand(numA,d); % set 1 given as matrix A
B = rand(numB,d); % set 2 given as matrix B

Answer 1

这里通常给出的答案是基于bsxfun(参见例如[1]).我提出的方法基于矩阵乘法,结果比我能找到的任何类似算法快得多:

helpA = zeros(numA,3*d);
helpB = zeros(numB,3*d);
for idx = 1:d
    helpA(:,3*idx-2:3*idx) = [ones(numA,1), -2*A(:,idx), A(:,idx).^2 ];
    helpB(:,3*idx-2:3*idx) = [B(:,idx).^2 ,    B(:,idx), ones(numB,1)];
end
distMat = helpA * helpB';

请注意: 对于常量,d可以for通过硬编码实现替换-loop,例如

helpA(:,3*idx-2:3*idx) = [ones(numA,1), -2*A(:,1), A(:,1).^2, ... % d == 2
                          ones(numA,1), -2*A(:,2), A(:,2).^2 ];   % etc.

评价:

%% create some points
d = 2; % dimension
numA = 20000;
numB = 20000;
A = rand(numA,d);
B = rand(numB,d);

%% pairwise distance matrix
% proposed method:
tic;
helpA = zeros(numA,3*d);
helpB = zeros(numB,3*d);
for idx = 1:d
    helpA(:,3*idx-2:3*idx) = [ones(numA,1), -2*A(:,idx), A(:,idx).^2 ];
    helpB(:,3*idx-2:3*idx) = [B(:,idx).^2 ,    B(:,idx), ones(numB,1)];
end
distMat = helpA * helpB';
toc;

% compare to pdist2:
tic;
pdist2(A,B).^2;
toc;

% compare to [1]:
tic;
bsxfun(@plus,dot(A,A,2),dot(B,B,2)')-2*(A*B');
toc;

% Another method: added 07/2014
% compare to ndgrid method (cf. Dan's comment)
tic;
[idxA,idxB] = ndgrid(1:numA,1:numB);
distMat = zeros(numA,numB);
distMat(:) = sum((A(idxA,:) - B(idxB,:)).^2,2);
toc;

结果:

Elapsed time is 1.796201 seconds.
Elapsed time is 5.653246 seconds.
Elapsed time is 3.551636 seconds.
Elapsed time is 22.461185 seconds.

有关数据点的维度和数量的更详细评估,请参阅下面的讨论(@comments).事实证明,在不同的环境中应该首选不同的算法.在非时间紧急情况下,只需使用该pdist2版本.

进一步发展: 人们可以考虑用基于相同原理的任何其他指标替换平方欧几里得:

help = zeros(numA,numB,d);
for idx = 1:d
    help(:,:,idx) = [ones(numA,1), A(:,idx)     ] * ...
                    [B(:,idx)'   ; -ones(1,numB)];
end
distMat = sum(ANYFUNCTION(help),3);

然而,这非常耗时.这可能是替换为更小的有用d的三维基质help由d2维矩阵.特别是d = 1它提供了一种通过简单矩阵乘法计算成对差异的方法:

pairDiffs = [ones(numA,1), A ] * [B'; -ones(1,numB)];

你有什么进一步的想法吗？