matlab - 计算点集和引用点集之间的马氏距离

Question

我有一个 n x p 矩阵 - mX，它由 R^p 中的 n 个点组成。

我有另一个 m x p 矩阵 - mY，它由 R^p 中的 m 个引用点组成。

我想创建一个 n x m 矩阵 - mD，即 Mahalanobis Distance矩阵。

D(i, j) 表示 Mahalanobis Distance在 mX, mX(i, :) 中的点 i 和 mY, mY(j, :) 中的点 j 之间。

即，计算如下:

mD(i, j) = (mX(i, :) - mY(j, :)) * inv(mC) * (mX(i, :) - mY(j, :)).';

其中 mC 是给定的马氏距离 PSD 矩阵。

在循环中很容易完成，有没有办法对其进行矢量化？

即，是一个函数，它的输入是 mX、mY 和 mC，它的输出是 mD 并且完全矢量化而不使用任何 MATLAB 工具箱？

谢谢你。

Answer 1

方法#1

假设资源无限，这是一个使用 bsxfun 的矢量化解决方案和 matrix-multiplication -

A = reshape(bsxfun(@minus,permute(mX,[1 3 2]),permute(mY,[3 1 2])),[],p);
out = reshape(diag(A*inv(mC)*A.'),n,m);

方法#2

这是一个尝试降低循环复杂性的组合解决方案 -

A = reshape(bsxfun(@minus,permute(mX,[1 3 2]),permute(mY,[3 1 2])),[],p);
imC = inv(mC);
out = zeros(n*m,1);
for ii = 1:n*m
    out(ii) = A(ii,:)*imC*A(ii,:).';
end
out = reshape(out,n,m);

sample 运行 -

>> n = 3;  m = 4;   p = 5;
mX = rand(n,p);
mY = rand(m,p);
mC = rand(p,p);
imC = inv(mC);
>> %// Original solution
for i = 1:n
    for j = 1:m
        mD(i, j) = (mX(i, :) - mY(j, :)) * inv(mC) * (mX(i, :) - mY(j, :)).'; %//'
    end
end
>> mD
mD =
      -8.4256       10.032       2.8929       7.1762
      -44.748      -4.3851      -13.645      -9.6702
      -4.5297       3.2928      0.11132       2.5998
>> %// Approach #1
A = reshape(bsxfun(@minus,permute(mX,[1 3 2]),permute(mY,[3 1 2])),[],p);
out = reshape(diag(A*inv(mC)*A.'),n,m);  %//'
>> out
out =
      -8.4256       10.032       2.8929       7.1762
      -44.748      -4.3851      -13.645      -9.6702
      -4.5297       3.2928      0.11132       2.5998
>> %// Approach #2
A = reshape(bsxfun(@minus,permute(mX,[1 3 2]),permute(mY,[3 1 2])),[],p);
imC = inv(mC);
out1 = zeros(n*m,1);
for ii = 1:n*m
    out1(ii) = A(ii,:)*imC*A(ii,:).';  %//'
end
out1 = reshape(out1,n,m);
>> out1
out1 =
      -8.4256       10.032       2.8929       7.1762
      -44.748      -4.3851      -13.645      -9.6702
      -4.5297       3.2928      0.11132       2.5998

相反，如果您有:

mD(j, i) = (mX(i, :) - mY(j, :)) * inv(mC) * (mX(i, :) - mY(j, :)).';

解决方案将转换为接下来列出的版本。

方法#1

A = reshape(bsxfun(@minus,permute(mY,[1 3 2]),permute(mX,[3 1 2])),[],p);
out = reshape(diag(A*inv(mC)*A.'),m,n);

方法#2

A = reshape(bsxfun(@minus,permute(mY,[1 3 2]),permute(mX,[3 1 2])),[],p);
imC = inv(mC);
out1 = zeros(m*n,1);
for i = 1:n*m
    out(i) = A(i,:)*imC*A(i,:).';  %//'
end
out = reshape(out,m,n);

sample 运行 -

>> n = 3; m = 4; p = 5;
mX = rand(n,p);    mY = rand(m,p);     mC = rand(p,p);  imC = inv(mC);
>> %// Original solution
for i = 1:n
    for j = 1:m
        mD(j, i) = (mX(i, :) - mY(j, :)) * inv(mC) * (mX(i, :) - mY(j, :)).'; %//'
    end
end
>> mD
mD =
      0.81755      0.33205      0.82254
       1.7086       1.3363       2.4209
      0.36495      0.78394     -0.33097
      0.17359       0.3889      -1.0624
>> %// Approach #1
A = reshape(bsxfun(@minus,permute(mY,[1 3 2]),permute(mX,[3 1 2])),[],p);
out = reshape(diag(A*inv(mC)*A.'),m,n);  %//'
>> out
out =
      0.81755      0.33205      0.82254
       1.7086       1.3363       2.4209
      0.36495      0.78394     -0.33097
      0.17359       0.3889      -1.0624
>> %// Approach #2
A = reshape(bsxfun(@minus,permute(mY,[1 3 2]),permute(mX,[3 1 2])),[],p);
imC = inv(mC);
out1 = zeros(m*n,1);
for i = 1:n*m
    out1(i) = A(i,:)*imC*A(i,:).';  %//'
end
out1 = reshape(out1,m,n);
>> out1
out1 =
      0.81755      0.33205      0.82254
       1.7086       1.3363       2.4209
      0.36495      0.78394     -0.33097
      0.17359       0.3889      -1.0624