How to Vectorize this Code?

I have this 147000 x 3 array, where each index represents a 3D vector. I want to use the vector at each index to find the angle between it and a predefined vector "b0" [ 0 1 0], this results in a 147000 x 1 array where each index has a single angle. I have been able to do this by using a normal for loop, however I was wondering if there was a way to vectorize this so it will be faster. I have been trying, but it either says that it can't do the dot product because arrays are of incompatible size, or it will only calculate the angle of the first index, then copy that 147000 times into the Nx1 array.

Thanks for any help.

[James Tursa]

You can use a matrix multiply with b0 as a column vector to get the dot products. E.g.,

d = n3Array * b0(:); % the (:) syntax results in a column vector

[Krishan Bansal]

@James Tursa when I do that the results aren't right, I'm probably doing something wrong still. But now b0angles is filled with values around 89. The final result should look like this.

[Matt J]

@Krishan Bansal Please post your code and output as text, not images, to facilitate copy/pasting.

[Krishan Bansal]

@Matt J ok the code is just:

n3Array = squeeze(Orientation4D(:, 1, 1,1:3));

N = length(n3Array);

b0angles = zeros(N,1);

for i = 1:N

b0angles(i) = rad2deg(acos(dot(n3Array(i,:), b0)/(norm(n3Array(i,:)) * norm(b0))));

end

Answer 1

[Matt J]

n3Array=normalize(n3Array,2,'n');
b0=b0/norm(b0);

b0angles=acosd(n3Array*b0);

[Krishan Bansal]

Wow somehow this worked, all I had to do was squeeze b0 angles again because acosd(n3Array*b0) produced an Nx3 matrix where the 1st and 3rd columns all have 90, but the 2nd column has the correct angles.

Thanks so much!

[James Tursa]

You might be interested in this post also:

https://www.mathworks.com/matlabcentral/answers/101590-how-can-i-determine-the-angle-between-two-vectors-in-matlab?s_tid=srchtitle

[Krishan Bansal]

@James Tursa Thanks. I actually found the formula from that post in the first place. The reason I used the cosine one is because when I used the tangent one the person I'm developing this script for told me he wanted me to use the cosine formula instead. I'm not sure why.

Now I'm just trying to understand why Matt's method works. Because I have never used normalize(). Not sure what the 'n' parameter does.

[James Tursa]

'n' picks the 'norm' option for the method to use for normalization. By normalizing everything ahead of time you don't have to divide by the vector norms later. The 2 argument means the normalization will happen along the 2nd dimension (i.e., effectively normalizing the rows).

You should ask the person you are developing the script for why he wants to use a numerically less robust method ...

[Krishan Bansal]

@James Tursa ok that makes much more sense, I guess trying to do everything in one step was not possible using the vectorized version of the loop. It is much faster than what I wrote originally though! I wonder if I put the dot product back into the formula I won't have to manually squeeze out the 1st and 3rd rows after doing acosd(n3array*b0).

[James Tursa]

If n3Array is Nx3 and b0 is 1x3, why aren't you doing n3Array*b0(:) as suggested earlier to calculate the dot products between the rows and b0? Why would there be any further manipulation needed?

[Krishan Bansal]

@James Tursa When I tried that it didn't work. If I do:

d = n3Array * b0(:);

i = 1:N;

b0angles(i) = rad2deg(acos(d(i)/(norm(n3Array(i,:)) * norm(b0))));

I end up with an Nx1 matrix, but every value is just some number close to 89, which is not right. Not sure if this is what you meant to do.

[James Tursa]

I meant using Matt's approach without the loop but doing this

b0angles = acosd(n3Array*b0(:));

[Krishan Bansal]

@James Tursa b0anglesMethod2=acosd(n3Array*b0(:)); creates an Nx1 complex double where every element is 0 for some reason.

[James Tursa]

This creates complex doubles?

n3Array = squeeze(Orientation4D(:, 1, 1,1:3));
b0 = [0 1 0];
n3Array = normalize(n3Array,2,'n');
b0 = b0/norm(b0);
b0angles = acosd(n3Array*b0(:));

[Krishan Bansal]

@James Tursa That worked perfectly, I forgot to use the normalized version of n3 I think. Thanks so much! This is super concise and efficient code.