Computing Mahalanobis Distance Between Set of Points and Set of Reference Points

Question

I have an n x p matrix - mX which is composed of n points in R^p.

I have another m x p matrix - mY which is composed of m reference points in R^p.

I would like to create an n x m matrix - mD which is the Mahalanobis Distance matrix.

D(i, j) means the Mahalanobis Distance between point i in mX, mX(i, :) and point j in mY, mY(j, :).

Namely, is computes the following:

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

Where mC is the given Mahalanobis Distance PSD Matrix.

It is easy to be done in a loop, is there a way to vectorize it?

Namely, is the a function which its inputs are mX, mY and mC and its output is mD and fully vectorized without using any MATLAB toolbox?

Thank You.

Could you add your loop-based solution and some realistic values for n,m and p to your question? At least the ideas I have to vectorize it are much less memory efficient as a temporary variable of size m x n x p would be created. Would this fit into the memory? — Daniel, Aug 09 '15 at 13:38
What values for n,m,p are realistic for your problem? Vectorizing is not the only technique to increaze the performance, maybe something else or a combination is the best solution. — Daniel, Aug 09 '15 at 21:04
Daniel, thank you for your interest. I want to have a vectorized code. Think there's unlimited memory so the size doesn't matter. I had the "Loop" code to begin with. I want a vectorized code and then I will test which one is superior in which cases. — Royi, Aug 10 '15 at 05:38
I'm not sure that I have properly understood the problem, but it seems to me that you can compute eigenvalue decomposition of `mC`, and then transform all the points into the space where `mC` is diagonal. As a result, computing distance between two points would take *O(p)* instead of *O(p^2)* time. Also, this code can surely be vectorized by SSE/AVX intrinsics, but I cannot say anything about MATLAB. — stgatilov, Aug 10 '15 at 11:00
@stgatilov, Go ahead, I would like to see that. Use "Pseudo" code. — Royi, Aug 10 '15 at 11:12

score 1 · Answer 1 · answered Aug 09 '15 at 20:10

1

Here is one solution that eliminates one loop

function d = mahalanobis(mX, mY)

    n = size(mX, 2);
    m = size(mY, 2);
    data = [mX, mY];
    mc = cov(transpose(data));

    dist = zeros(n,m);
    for i = 1 : n
        diff = repmat(mX(:,i), 1, m) - mY;
        dist(i,:) = sum((mc\diff).*diff , 1);
    end
    d = sqrt(dist);

end

You would invoke it as:

d = mahalanobis(transpose(X),transpose(Y))

answered Aug 09 '15 at 20:10

dpmcmlxxvi

1,292
1
9
13

Is there a way to do it without any loop at all? Thank You. – Royi Aug 10 '15 at 05:36

Divakar · Accepted Answer · 2015-08-12T08:19:14.350

Approach #1

Assuming infinite resources, here's one vectorized solution using bsxfun and 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);

Approach #2

Here's a comprise solution trying to reduce the loop complexity -

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 run -

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

Instead if you had :

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

The solutions would translate to the versions listed next.

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

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
    out(i) = A(i,:)*imC*A(i,:).';  %//'
end
out = reshape(out,m,n);

Sample run -

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

I think your fully vectorized solution gives the transpose of mD. Am I right? Thank You. Could you explain the fully vectorized solution? — Royi, Aug 10 '15 at 11:14
The output from these solutions would be of size `n x m`, i.e. number of rows in `mX` x number of rows in `mY`. So, it should be same as `mD`. Here's a sample run - http://pastebin.com/V9K6Mjdf — Divakar, Aug 10 '15 at 11:37
Yes. I made a mistake in my code, it should have been `mD(j, i) = (mX(i, :) - mY(j, :)) * inv(mC) * (mX(i, :) - mY(j, :)).';`. If you could add what would you change (Don't delete how it is now, just add) I would be grateful. I +1 your answer. I want more people to answer. since you are the first I will mark you in few days. Thank You!!! — Royi, Aug 10 '15 at 13:29
@Drazick Could you edit your question to reflect that. I would edit my answer accordingly after your edit? — Divakar, Aug 10 '15 at 13:38
@Divkar, that's the thing, I'd rather have both of them online. Hence if you could add "In the case of..." It would be great. Thank You. — Royi, Aug 10 '15 at 15:39
Regarding the alternative solution #Approach 1, I don't think that reshaping with m instead of n will do. The Transpose is the solution, yet doing it that way won't yield the transpose. Unless I'm missing another change. — Royi, Aug 12 '15 at 06:33
@Drazick Think there was a typo. Check out the just added sample runs? — Divakar, Aug 12 '15 at 08:19

score 1 · Answer 3 · answered Aug 10 '15 at 12:23

Reduce to L2

It seems that Mahalanobis Distance can be reduced to ordinary L2 distance if you are allowed to preprocess matrix mC and you are not afraid of numerical differences.

First of all, compute Cholesky decomposition of mC:

mR = chol(mC)      % C = R^t * R, where R is upper-triangular

Now we can use these factors to reformulate Mahalanobis Distance:

(Xi-Yj) * inv(C) * (Xi-Yj)^t = || (Xi-Yj) inv(R) ||^2 = ||TXi - TYj||^2
where:  TXi = Xi * inv(R)
        TYj = Yj * inv(R)

So the idea is to transform points Xi, Yj to TXi, TYj first, and then compute euclidean distances between them. Here is the algorithm outline:

Compute mR - Cholesky factor of covariance matrix mC (takes O(p^3) time).
Invert triangular matrix mR (takes O(p^3) time).
Multiply both mX and mY by inv(mR) on the right (takes O(p^2 (m+n)) time).
Compute squared L2 distances between pairs of points (takes O(m n p) time).

Total time is O(m n p + (m + n) p^2 + p^3) versus original O(m n p^2). It should work faster when 1 << p << n,m. In such case step 4 would takes most of the time and should be vectorized.

Vectorization

I have little experience of MATLAB, but quite a lot of SIMD vectorization on x86 CPUs. In raw computations, it would be enough to vectorize along one sufficiently large array dimension, and make trivial loops for the other dimensions.

If you expect p to be large enough, it may probably be OK to vectorize along coordinates of points, and make two nested loops for i <= n and j <= m. That's similar to what @Daniel posted.

If p is not sufficiently large, you can vectorize along one of the point sequences instead. This would be similar to solution posted by @dpmcmlxxvi: you have to subtract single row of one matrix from all the rows of the second matrix, then compute squared norms of the resulting rows. Repeat n times ( or m times).

As for me, full vectorization (which means rewriting with matrix operations instead of loops in MATLAB) does not sound like a clever performance goal. Most likely partially vectorized solutions would be optimally fast.

Now all we need is the most efficient way to calculate the L2 Distance between each item of 2 sets of vectors. — Royi, Aug 10 '15 at 16:55

score 0 · Answer 4 · answered Aug 09 '15 at 14:09

I came to the conclusion that vectorizing this problem is not efficient. My best idea for vectorizing this problem would require m x n x p x p working memory, at least if everything is processed at once. This means with n=m=p=152 the code would already require 4GB Ram. At these dimensions, my system can run the loop in less than a second:

mD=zeros(size(mX,1),size(mY,1));
ImC=inv(mC);
for i=1:size(mX,1)
    for j=1:size(mY,1)
        d=mX(i, :) - mY(j, :);
        mD(i, j) = (d) * ImC * (d).';
    end
end

OK, Yet could you offer the vectorized form? Thank You. – Royi Aug 09 '15 at 14:26 — Royi, Aug 09 '15 at 14:26

Computing Mahalanobis Distance Between Set of Points and Set of Reference Points

4 Answers4

Reduce to L2

Vectorization