How to calculate the pairwise distince between two dataset without any loop

I have two datasets, say, and , where m and n is the number of observations. if I want to calculate the following distinance:

where and , and . and is symmetric

--------------------------------------------------------------------------------------------------------------------------------

update:

for example:

X = [1 2 1;3 4 2;5 6 3]
Y = [5 6 1;7 8 2]

and the final results should be , where each element of Matrix C is calculated by

for example,

h = X(2, 1:2) - Y(1, 1:2);
u = X(2, 3) - Y(1, 3);
C(2,1) = (h-u)*inv(eye(2) + sigma*u^2))*(h-u)'

and I need to calculate the Matrix C

How can I implement this in a vectorized way (without any loops)?

Thanks!

[Image Analyst]

Can you say it in English instead of mathematical jargon? And attach your "datasets" so we know what they really are.

TUTORIAL: How to ask a question (on Answers) and get a fast answer

Meanwhile the sum symbol is done by the sum() function, and the T transpose is done by the apostrophe symbol. The -1 can be done with inv(), so give it a try, like

h = X(i, 1:2) - Y(j, 1:2);
u = X(i, 3) - Y(j, 3);
distance = (h-u) * inv(I2 + sum(u.^2)) * (h-u)'

Or something like that.

[Edward]

Hello！

Thanks for answering, I have updated the problem and make it more clear. Could you help me implement it without loop?

Answer 1

[Bruno Luong]

If you have latest MATLAB release (for pagemtimes) and download this Multiple same size solver FEX

X = [1 2 1;3 4 2;5 6 3]
Y = [5 6 1;7 8 2]
sigma = rand(2);

m = size(X,1);
n = size(Y,1);
XX = permute(X, [3 2 1 4]);
YY = permute(Y, [3 2 4 1]);
XmY = reshape(XX-YY,1,[],m*n);
h = XmY(:,1:2,:);
u = XmY(:,3,:);
hmu = h-u;
I2 = eye(2)
A = I2 + sigma.*(u.^2);
% https://www.mathworks.com/matlabcentral/fileexchange/24260-multiple-same-size-linear-solver?s_tid=srchtitle
C = SliceMultiSolver(A, permute(hmu,[2 1 3])); % WARNING: conjugate needed for complex?
C = pagemtimes(hmu, C);
C = reshape(C,[m n])

[Edward]

Hello! Bruno, Thank you for your solutions. I have updated my Matlab and download your package to implement this code. It works!

Frankly, this equation is one part of my final calculation. At first, I want to learn how to code it without loops. and later try code the left part by myself. However, I still couldn't understand how it works. (I am trying...

Actually, my whole calculation is

where and , we now can set randomly, say [2,1] (Identity matrix) and is positive symmetric (e.g. [2 1;1 2]). is the determinant of a Matrix.

Could you help me code this? Any help would be greatly appreciated.

Thanks!

[Bruno Luong]

To compute the determinant of multiple matrices once, you need another FEX MultipleQR

But frankly if the vectorized code is too obscure for you, you are off better using for-loop for the maintainability.

X = rand(10,3);
Y = rand(8,3);
sigma = rand(2);
sigma = sigma'*sigma;
I2 = eye(2);

m = size(X,1);
n = size(Y,1);
XX = permute(X, [3 2 1 4]);
YY = permute(Y, [3 2 4 1]);
XmY = reshape(XX-YY,1,[],m*n);
h = XmY(:,1:2,:);
u = XmY(:,3,:);
hmu = h-u;
A = I2 + sigma.*(u.^2);
% https://www.mathworks.com/matlabcentral/fileexchange/24260-multiple-same-size-linear-solver?s_tid=srchtitle
C = SliceMultiSolver(A, permute(hmu,[2 1 3])); % WARNING: conjugate needed for complex?
C = pagemtimes(hmu, C);

% Compute det(A(:,:,k)) for k=1,2,... without loop
% https://www.mathworks.com/matlabcentral/fileexchange/68976-multipleqr
[~,R] = MultipleQR(A);
p = size(R,1);
i = 1:p;
i = sub2ind(size(R),i,i,ones(size(i)));
ilin = i + (p*p)*reshape((0:m*n-1),1,1,[]);
detA = (-1)^p*prod(R(ilin),2);

C = exp(-C)./sqrt(detA); % EDIT correct bug of missing minus sign
C = reshape(C,[m n])

[Edward]

Hello! Bruno, Thank you for your reply and help! I wrote a loop code to implement it. I found our results is not the same. Could you help check my code? Thanks!

my_C = zeros(m,n);
for i = 1:m
    for j = 1:n
        h = X(i,1:2)-Y(j,1:2);
        u = X(i,3)-Y(j,3);
        d = h-u;
%         d = h-miu*u;
        my_C(i,j) = exp(-d*inv((eye(2)+sigma*u^2))*d')/sqrt(det(eye(2)+sigma*u^2));
    end
end
my_C

Beside, if I wanna to change be (see the equations), Can I only change your code like the following?

hmu = h-miu.*u;

[Bruno Luong]

The bug is on my side I just edit it.

No I don't change anymore my code. I hate question that is modified. I don't have free time to support questions that are constantly modified.

And beside if you don't feel comfortable of maintaining the vectorizing code, then I strongly suggest you NOT using it.

[Edward]

Thank you! Bruno.