increase compute speed compute angle between 2 vectors

hi, I'm student in a space compagny and I have built a matlab soft to compute orbit , but the run take more than 48h, in fact, one function was call more than billion time, and also I search to win some milliseconde in my funtion.

this funtion compute angle between 2 vectors input : 2 vectors v & u (3 by n) output : angle between u and v in rad (n by 1)

here after the code , but I don't find more solution to optimize it :

function angle = searchAngle(u, v)
norm = @(v) sqrt(sum(v.^2, 2));
dot = @(u, v) sum(u .* v, 2);

cross = @(a, b) [ a(:,2) .* b(:,3) - a(:,3) .* b(:,2), ...
      a(:,3) .* b(:,1) - a(:,1) .* b(:,3), ...
      a(:,1) .* b(:,2) - a(:,2) .* b(:,1) ];

normVect = norm(u) .* norm(v);

dotVect = dot(u, v); 
threshold = normVect * 0.9999;

idx1 = dotVect > threshold;
axis = cross(v(idx1,:), u(idx1,:));
angle(idx1) = asin(norm(axis) ./ normVect(idx1));

idx2 = dotVect < -threshold;
axis = cross(v(idx2,:), u(idx2,:));
angle(idx2) = pi - asin(norm(axis) ./ normVect(idx2));

idx = ~(idx1 | idx2);
angle(idx) = acos(dotVect(idx) ./ normVect(idx));

end

thx for any help

[Titus Edelhofer]

Hi,

one question: is your function searchAngle called billions of times? The function itself looks good, even for large u and v. Maybe you can arrange optimizing the calling function to call searchAngle less often but with larger blocks ...?

Titus

[franck lepapaix]

it was my first idea to reduce the time. but I can't do that , in fact I make computation of many objects in space and for position compution (orbit) I compute many point each 10 minutes during one year, and for that I can't increase blocks, my vectors are tri-dimensionnal !(thk mister Kepler)

[Rena Berman]

(Answers dev) Restored question.

Answer 1

[Jan]

The indirection of anonymous functions costs time. So either use the built-in functions with the same names cross, norm and dot, or hard code the functions directly.

Instead of the expensive trick to determine the positions of instabilities in the ASIN and ACOS methods, use a stable method directly:

atan2(norm(cross(N1 x N2)), dot(N1, N2))

Where N1 and N2 are the normalized input vectors.

N1      = bsxfun(@rdivide, a, sqrt(sum(a .* a ,1)))
N2      = bsxfun(@rdivide, b, sqrt(sum(b .* b ,1)))
N1dotN2 = N1(:, 1) .* N2(:, 1) + N1(:, 2) .* N2(:, 2) + N1(:, 3) .* N2(:, 3);
N1xN2   = [(N1(:, 2) .* N2(:, 3) - N1(:, 3) .* N2(:, 2)), ...
         (N1(:, 3) .* N2(:, 1) - N1(:, 1) .* N2(:, 3)), ...
         (N1(:, 1) .* N2(:, 2) - N1(:, 2) .* N2(:, 1))];
Angle   = atan2(sqrt(sum(N1xN2 .* N1xN2, 1)), N1dotN2);

[franck lepapaix]

sorry, but....

on one small test of computation, with my function the compute time is : 48 sec with your function : 52 sec... more than 4sec...

I will check your idea about hard code funtion use directly...

[Titus Edelhofer]

Which version of MATLAB are you using? R2015b or an older one? That could make quite some difference for your code ...

[franck lepapaix]

matlab 2015a

[Jan]

@franck: Do I understand correctly, that a and b are [1 x 3] vectors? Then try:

N1      = a ./ sqrt(a * a.');
N2      = b ./ sqrt(b * b.');
N1dotN2 = N1(1) .* N2(1) + N1(2) .* N2(2) + N1(3) .* N2(3);
N1xN2   = [(N1(2) .* N2(3) - N1(3) .* N2(2)), ...
           (N1(3) .* N2(1) - N1(1) .* N2(3)), ...
           (N1(1) .* N2(2) - N1(2) .* N2(1))];
Angle   = atan2(sqrt(N1xN2 * N1xN2.')), N1dotN2);

I have this a C-Mex function also, but without the initial normalization. But the main drawback remains, that teh code is called with wingle vectors instead of a vectorized calling with thousands of elements. A vectorization would be the best improvement.

[franck lepapaix]

thx for your help, my vectors are [222651 x 3] size !! that's why I use n in my post ! I have no idea about the use of vectorization in matlab, I think this is the next step for me to optimize the speed..

[Jan]

And you provide this [222651 x 3] matrix as input, or do you call the function in a loop for each [1 x 3] vector? The command norm(u) in your code seems to imply, that you call it for vectors. The code in my answer can process the complete matrix in one call, which should be substantially faster. Even a fast C-Mex function, which avoids the creation of large temporary arrays, would suffer from beeing called hundret thousands of times due to the calling overhead.