Creating a Matrix from Distances between Vectors

I start with a 3 x 46 array where each 3 x1 column in the array is viewed as a position vector of a point, let's say the first vector in the array is p_1, the next one is p_2.

I also have a 3 x 100 array, which I am viewing as the position vectors of 100 points (where these points are always of a lower 'radius' than the points in the 3 x 46 array). I basically need to create a matrix which looks as follows, where each distance is the row vector obtained by vector subtraction of the two vectors:

[Distance from p_1 to first vector in the 3 x 100 array Distance from p_1 to second vector etc ]

[Distance from p_2 to first vector in the 3 x 100 array Distance from p_2 to second vector etc ]

[ Distance from p_3 ..... ]

and so on, until all the positions in the 3 x 46 array are covered, let me know if my intention is not clear.

Answer 1

[Guillaume]

This can be achieved easily without a loop:

%demo data
m = 46
n = 100
A = rand(3,m)
B = rand(3,n)

distance = sqrt(sum((permute(A, [2, 3, 1]) - permute(B, [3, 2, 1])).^2, 3));

The result is a mxn matrix. It basically moves the (X, Y, Z) of each point into the 3rd dimension for both matrix, and move the 2nd dimension of the 1st matrix into the 1st. So points in A are along the rows, and points in B are along the column (and coordinates across the 3rd for both). We now have compatible array sizes, so can create a 3D matrix of (Xa-Xb, Ya-Yb, Za-Zb) for all the points. Square that, sum across (X, Y, Z) so it reduces to a 2D matrix, and take the square root for the distance.

Or, if you have the stats toolbox, you can use pdist2.

edit: Note that if you want to store the vectors between the points as suggested by Jon instead of the distance (norm of the vectors), then it's simply:

vectors = permute(A, [2, 3, 1]) - permute(B, [3, 2, 1])

vectors(p1, p2, :) is the vector between point A(:, p1) and point B(:, p2).

[Tom]

Hi Guillaume, I've realised that in fact it's not the distance vector which goes into the matrix each time, but a function into which you plug each distance to get a 3 x 3 matrix (where rij is the distance vector):

J=1/(8*pi*0.7*sqrt(2/pi))*(eye(3)/norm(rij)+(rij'*rij)/((norm(rij))^3));

So instead of a row each time for the distance vector, each time you get a 3 x 3 matrix and the final size of the matrix ends up being 300 x 138, rather than 300 x 46. Let me know if this is not clear what I mean.

[Jon]

@Guillaume Thank you for explaining how you can eliminate the loop using permute. This is a very elegant approach, that I have definitely learned something from! Occassionally I have stored data in multidimensional arrays when it was natural to do so. For actual calcuation I have generally confined myself to thinking in terms of 1d and 2d arrays. I can see now that there is a whole other, quite powerful way of thinking about this.

[Guillaume]

@Tom.

If it involves matrix multiplication then I'm afraid you'll have to use loops.

Note that instead of a 300x138 matrix, I'd create a 100x46x3x3 matrix, or a 100x46 cell array of 3x3 matrices (The 4D matrix is more efficient).

@Jon,

[Guillaume]

Here's an easy way to create the 100x46 cell array of 3x3 matrices:

%A: 3xm matrix of points
%B: 3xn matrix of points

%distance formula, as an anonymous function
J = @(rij) 1/(8*pi*0.7*sqrt(2/pi))*(eye(3)/norm(rij)+(rij'*rij)/((norm(rij))^3));
[idA, idB] = ndgrid(1:size(A, 2), 1:size(B, 2));
distance = arrayfun(@(cola, colb) J(B(:, colb) - A(:, cola)), idA, idB, 'UniformOutput', false)

[Jon]

Also very nice! It took me awhile to understand this example but I now have a much deeper appreciation for the power of arrayfun to solve this kind of problem.

@Tom, as you finally wanted a 3xn by 3xm array, you could add one more line to turn the cell array into a matrix, and transpose it

distance = cell2mat(distance)'

[Tom]

Yes, I was just about to ask how to get the cell back into a matrix so thanks for that. I've got my solution, thanks for explaining these things for me, I didn't realise you could do so much in MATLAB with only a few lines of code.

Also is there a way that I can double-check that the distance formula script between all the positions is doing what I think it should be doing, as obviously the size of everything is right but the entries might not be right if I failed to explain my intention properly.

[Tom]

I now have a slightly harder version of the problem, I now need to do the same thing but create a 'distance' matrix which is 400 x 184 (so the matrices in the cell array are now 4 x 4).

A B C D

E F G H

I J K L

M N O P

There are now 4 functions of the distances between points rij, call them f1, f2, f3 and f4. The matrix formed by A B C E F G I J K corresponded to the case to our first example, where we plugged in a distance to the function each time and got a 3 x 3 matrix.

Now in this example, the same must be done, but the function f1 returns the matrix formed by

A B C E F G I J K when you plug in each distance between points, f2 returns the column vector formed by D H L when you plug in the same distance, f3 returns the row vector formed by

M N O when you plug in the distance, and f4 returns the scalar P when you plug in the distance, so this way you are repeating what was done before but now making a 4 x 4 matrix each time. Let me know if this is not clear, as I need to be clear when trying to explain this stuff.

[Guillaume]

It should be pretty much the same:

%A, B: matrices of points
f{1} = @(rij) ...  %returns a 3x3 matrix
f{2} = @(rij) ...  %returns a 3x1 column vector
f{3} = @(rij) ...  %returns a 1x3 row vector
f{4} = @(rij) ...  %returns a scalar
fall = @(rij) [f{1}(rij), f{2}(rij); f{3}(rij), f{4}(rij)];

[idA, idB] = ndgrid(1:size(A, 2), 1:size(B, 2));
distance = arrayfun(@(cola, colb) fall(B(:, colb) - A(:, cola)), idA, idB, 'UniformOutput', false)
distasmat = cell2mat(distance);   

[Tom]

Another slight subtlety is that some of the functions depend on the normal vector of the point on the radial line, if the position vector of a point on the radial line is called p, this as defined as p divided by the modulus of p.

That means f{1} might not just be a function of rij (the distance from a point on the radial line to the inner points), but also a function of the position vector of the radial point which is being considered, I suppose this can also be added in to the arrayfun?

[Tom]

@Guillame: I've afraid I can't get the above to work, I will try later but I'm running out of time and would be willing to just loops now. Basically I have these functions, where rij is the distance we discussed and normal is just the position vector for one of the points on the radial line divided by the modulus.

function [GV]=gradletSlip(rij,normal)
global Kn;   
GV=((1/(8*pi*Kn)*(eye(3)/norm(rij)+(rij'*rij)/((norm(rij))^3)))+sqrt(pi/2)*((-1/(4*pi))*(rij'*normal* ...
  (eye(3)/norm(rij)^3 -3*(rij'*rij)/norm(rij)^5)*(eye(3)-normal'*normal)))+(3*Kn/(8*pi)*(eye(3)/norm(rij)^3-3* ...
   (rij'*rij)/norm(rij)^5)*(eye(3)-normal'*normal))/5+(3*Kn/(pi*20)*(-eye(3)/norm(rij)^3+3*(rij'*rij)/norm(rij)^5)));  
end

function [FV]=thermletSlip(rij,normal)
global Kn; 
FV=((-Kn/(5*pi)*normal*(eye(3)/norm(rij)^3-3*(rij'*rij)/norm(rij)^5)*(eye(3)-normal'*normal))') ...
   *sqrt(pi/2)+((1/(4*pi*norm(rij)^3)*rij*(eye(3)-normal'*normal))')/5;  
end

function [GT]=gradletJump(rij,normal)
global Kn;
GT=((3*Kn/(8*pi)*(eye(3)/norm(rij)^3-3*(rij'*rij)/norm(rij)^5)*normal')*sqrt(pi/8)+(-1/(4*pi)*rij'* ...
    ((normal*(eye(3)/norm(rij)^3-3*(rij'*rij)/norm(rij)^5))*normal'))/4);    
end

function [FT]=thermletJump(rij,normal)
global Kn;
FT=((1/(15*pi*Kn*norm(rij)))+(normal*rij'/(4*pi*norm(rij)^3))*sqrt(pi/8)+(-Kn/(5*pi) ...
    *((normal*(eye(3)/norm(rij)^3-3*(rij'*rij)/norm(rij)^5))*normal'))/4);  
end

which evaluate to give matrices, columns in the order you gave. I basically just want to evaluate these in the way I explained above with loops if necessary to get a 400 x 138 matrix. I am trying to create the matrix with this code where p is 3 x 46 and sing is 3 x 100 (where totalNodes = 100), but I think something is wrong with the indices and I am getting a 400 x 400 matrix. Bear in mind some functions have 'normal' or 'in' which is just the vector for the position divided by the modulus.

for II=1:46 
    for IJ=1:totalNodes   
        in=(p(1:3,II)/(norm(p(1:3,II))))'; 
        r=(p(1:3,II)-sing(1:3,IJ))'; 
 A((1+(II-1)*3):(3+(II-1)*3),(1+(IJ-1)*3):(3+(IJ-1)*3))=gradletSlip(r,in)';
    end
    for IJ=3*totalNodes+1:4*totalNodes  
        in=(p(1:3,II)/(norm(p(1:3,II))))';  
        r=(p(1:3,II)-sing(1:3,IJ-3*totalNodes))'; 
        A((1+(II-1)*3):(3+(II-1)*3),IJ)=thermletSlip(r,in);
    end
end
for II=1:46 
     for IJ=1:totalNodes   
        in=(p(1:3,II)/(norm(p(1:3,II))))'; 
        r=(p(1:3,II)-sing(1:3,IJ))'; 
    A(3*totalNodes+II,(1+(IJ-1)*3):(3+(IJ-1)*3))=gradletJump(r,in);
     end
      for IJ=3*totalNodes+1:4*totalNodes  
     in=(p(1:3,II)/(norm(p(1:3,II))))';  
        r=(p(1:3,II)-sing(1:3,IJ-3*totalNodes))'; 
      A(3*totalNodes+II,IJ)=thermletJump(r,in);
      end
end

[Guillaume]

arrayfun is just a for loop in disguise. It just saves you on having to compute the bounds yourself (more on that later). With either method, if speed is not an issue, I'd create a cell array rather than index into the destination. The code is a lot easier to read this way.

So, I'd do:

matrixbuilder = @(rij, normal) [gradletSlip(rij, normal), thermletSlip(rij, normal);
                                gradletJump(rij, normal), thermletJump(rij, normal)];
[colp, colsing] = ndgrid(1:size(colp, 2), 1:size(colsing, 2))                            
distance = arrayfun(@(colp, colsing) matrixbuilder((p(:, colp) - sing(:, colsing))', p(:, colp)' / norm(p(:, colp))), colp, colsing, 'UniformOutput', false)
distmat = cell2mat(distance);

Answer 2

[Jon]

So from your description you actually want to find and store the vectors that go from each point in the first matix to each point in the second matrix. Here is a way you could do that

m = 46
n = 100
A = rand(3,m)
B = rand(3,n)

% preallocate 3d matrix to hold results
C = zeros(m,3,n);

% loop through columns of  A calculating "distance" vectors and storing
% them
for k = 1:m
    C(k,:,:) = B - A(:,k);
end

% to look, for example, at "distance vectors" from the 3rd column of A to all of the
% columns of B you can use
disp(squeeze(C(3,:,:)))

Note, the squeeze function gets rid of extraneous single dimension in multidimensional arrays, which would otherwise make it difficult to see the result

If instead you just want the scalar Euclidean distance from each point in A to each point in B, you could do the following

D = zeros(m,n);
for k = 1:m
    D(k,:) = vecnorm(B - A(:,k));
end

If you wanted both you could put both calculations in the same loop.

It may be possible to vectorize this further and eliminate the loop, but it is not immediately obvious to me how to do this. Maybe someone else has a suggestion.

[Tom]

I've realised that in fact it's not the distance vector which goes into the matrix each time, but a function into which you plug each distance to get a 3 x 3 matrix (where rij is the distance vector):

J=1/(8*pi*0.7*sqrt(2/pi))*(eye(3)/norm(rij)+(rij'*rij)/((norm(rij))^3));

So instead of a row each time for the distance vector, each time you get a 3 x 3 matrix and the final size of the matrix ends up being 300 x 138, rather than 300 x 46. Let me know if this is not clear what I mean, I don't mind using a loop if necessary.

[Jon]

I think this should be at least close to what you want. You should proof read the actual formulae to make sure.

% with 3d result for each pairing

% preallocate 4d array to hold results
J = zeros(m,n,3,3);
for i = 1:m
    r = B - A(:,k);
    rnorm = vecnorm(r);
    for j = 1:n
        J(i,j,:,:) =1/(8*pi*0.7*sqrt(2/pi))*(eye(3)/rnorm(j)+(r(:,j)'*r(:,j))/rnorm(j)^3);
        
    end
end

So now to see the 3d result from the pairing between the ith column of A and jth column of B you would look at J(i,j,:,:)

We will see if Guillaume can get rid of one or both of these loops using his approach, would be a good challenge. Unless you are doing huge matrices or performance is otherwise crucial I think the looping will be OK. Still nice if it could be vectorized.

[Tom]

I've read through and given it a go, but it does not seem to like the index k for the 3 x 100 array A.

[Jon]

Sorry there's a typo there. I cut and pasted from my earlier example, but now the index for the elements of A is i not k so instead you should have

m = 4
n = 5
A = rand(3,m)
B = rand(3,n)



% with 3d result for each pairing

% preallocate 4d array to hold results
J = zeros(m,n,3,3);
for i = 1:m
    r = B - A(:,i);
    rnorm = vecnorm(r);
    for j = 1:n
        J(i,j,:,:) =1/(8*pi*0.7*sqrt(2/pi))*(eye(3)/rnorm(j)+(r(:,j)'*r(:,j))/rnorm(j)^3);
        
    end
end

Note you will have to adjust m and n for your case. I just wanted a small example here so you could see the results

[Tom]

Is this right, as in my case if we call the 3 x 46 array P and the 3 x 100 array S, I need the distance between points with position vecotrs in P and points with position vectors in S (which are closer in radially), so it should be P - S?

Also, is there a way to avoid having a 4D array as I have not used these before? I just want a 300 x 138 matrix which I can multiply with a 300 x 1 array which I already have prepared, and that will give me my solution. (I don't mind if it is not efficient, it just needs to work).

[Jon]

I don't know what your application is so I can't give you any advice on whether you want to define the "distance" vectors as P-S or S-P. I don't think it changes the main logic that I outlined. If I'm understanding correctly it just changes some algebraic signs.

You could probably work out how to permute the 4d array that I have obtained above to get the 3n x 3m array that you want. Otherwise, here is how you can store it this way in the first place. You should check the details, to make sure it is right, but this should give you enough of an idea to go on.

m = 4
n = 5
A = rand(3,m)
B = rand(3,n)



% with 3d result for each pairing

% preallocate 2d array to hold results
J = zeros(n*3,m*3);
for i = 1:m
    r = B - A(:,i);
    rnorm = vecnorm(r);
    for j = 1:n
        startRow = (j-1)*3 + 1; % starting row where this 3x3 result is stored
        startCol = (i-1)*3 + 1; % starting column where this 3x3 result is stored
        J(startRow:startRow+2,startCol:startCol+2) = ...
            1/(8*pi*0.7*sqrt(2/pi))*(eye(3)/rnorm(j)+(r(:,j)'*r(:,j))/rnorm(j)^3);
        
    end
end

[Tom]

If it is a 4D array will I still be multiply it into the 300 x 1 column vector which I have ready to go?

[Guillaume]

Working in N-D is no more harder than working in 2-D. And in your case it would make it much easier to find the matrix related to a pair of point. You can either have the indices of the points as the first 2 indices, so:

distmat = squeeze(J(p1, p2, :, :));  %3x3 distance matrix between A(:, p1) and B(:, p2)

or as the last 2 indices:

distmat = J(:, :, p1, p2); %3x3 distance matrix between A(:, p1) and B(:, p2)

Storing the same in a 2D array, it's a lot less clear:

distmat = J((p1-1)*3 + (1:3), (p2-1)*3 + (1:3));

Another option, as I suggested is to store the matrices in a 2D cell array, in which case:

distmat = J{p1, p2};

The disadvantage of cell arrays is the memory and speed overhead. For an easy way to do this, see comment to my answer.