这里通常给出的答案是基于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.

对于更详细的评估w.r.t.数据点的维数和数量遵循下面的讨论(@comments).事实证明,在不同的设置中应该优先选择不同的algos.在非时间紧迫的情况下,只需使用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维的三维矩阵可能是有用的.特别是对于d = 1,它提供了一种通过简单矩阵乘法计算成对差的方法：

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

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