Ytukyg

Consider the following calculation of the tangent tangent correlation which is performed in a for loop

v1=rand(25,1);

v2=rand(25,1);

n=25;

nSteps=10;

mean_theta = zeros(nSteps,1);

for j=1:nSteps

    theta=;

    for i=1:(n-j)

        d = dot([v1(i) v2(i)],[v1(i+j) v2(i+j)]);

        n1 = norm([v1(i) v2(i)]);

        n2 = norm([v1(i+j) v2(i+j)]);

        theta = [theta acosd(d/n1/n2)];  



    end

    mean_theta(j)=mean(theta);

end

plot(mean_theta)

How can matlab matrix calculations be utilized to make this performance better?

Here is a full vectorized solution:

i = 1:n-1;

j = (1:nSteps).';

ij= min(i+j,n);

a = cat(3, v1(i).', v2(i).');

b = cat(3, v1(ij), v2(ij));



d = sum(a .* b, 3);

n1 = sum(a .^ 2, 3);

n2 = sum(b .^ 2, 3);

theta = acosd(d./sqrt(n1.*n2));

idx = (1:nSteps).' <= (n-1:-1:1);



mean_theta = sum(theta .* idx ,2) ./ sum(idx,2);

Result of Octave timings for my method,method4 from the answer provided by @CrisLuengo and the original method (n=250):

Full vectorized    : 0.000864983 seconds

Method4(Vectorize) : 0.002774 seconds

Original(loop)     : 0.340693 seconds

There are several things you can do to speed up your code. First, always preallocate. This converts:

theta = ;

for i = 1:(n-j)

    %...

    theta = [theta acosd(d/n1/n2)];

end

into:

theta = zeros(1,n-j);

for i = 1:(n-j)

    %...

    theta(i) = acosd(d/n1/n2);

end

Next, move the normalization out of the loops. There is no need to normalize over and over again, just normalize the input:

v = [v1,v2];

v = v./sqrt(sum(v.^2,2)); % Can use VECNORM in newest MATLAB

%...

    theta(i) = acosd(dot(v(i,:),v(i+j,:)));

This does change the output very slightly, within numerical precision, because the different order of operations leads to different floating-point rounding error.

Finally, you can remove the inner loop by vectorizing that computation:

i = 1:(n-j);

theta = acosd(dot(v(i,:),v(i+j,:),2));

Timings (n=25):

• Original: 0.0019 s


• Preallocate: 0.0013 s


• Normalize once: 0.0011 s


• Vectorize: 1.4176e-04 s

Timings (n=250):

• Original: 0.0185 s


• Preallocate: 0.0146 s


• Normalize once: 0.0118 s


• Vectorize: 2.5694e-04 s

Note how the vectorized code is the only one whose timing doesn't grow linearly with n.

Timing code:

function so

n = 25;

v1 = rand(n,1);

v2 = rand(n,1);

nSteps = 10;

mean_theta1 = method1(v1,v2,nSteps);

mean_theta2 = method2(v1,v2,nSteps);

fprintf('diff method1 vs method2: %gn',max(abs(mean_theta1(:)-mean_theta2(:))));

mean_theta3 = method3(v1,v2,nSteps);

fprintf('diff method1 vs method3: %gn',max(abs(mean_theta1(:)-mean_theta3(:))));

mean_theta4 = method4(v1,v2,nSteps);

fprintf('diff method1 vs method4: %gn',max(abs(mean_theta1(:)-mean_theta4(:))));



timeit(@()method1(v1,v2,nSteps))

timeit(@()method2(v1,v2,nSteps))

timeit(@()method3(v1,v2,nSteps))

timeit(@()method4(v1,v2,nSteps))



function mean_theta = method1(v1,v2,nSteps)

n = numel(v1);

mean_theta = zeros(nSteps,1);

for j = 1:nSteps

    theta=;

    for i=1:(n-j)

        d = dot([v1(i) v2(i)],[v1(i+j) v2(i+j)]);

        n1 = norm([v1(i) v2(i)]);

        n2 = norm([v1(i+j) v2(i+j)]);

        theta = [theta acosd(d/n1/n2)];

    end

    mean_theta(j) = mean(theta);

end



function mean_theta = method2(v1,v2,nSteps)

n = numel(v1);

mean_theta = zeros(nSteps,1);

for j = 1:nSteps

    theta = zeros(1,n-j);

    for i = 1:(n-j)

        d = dot([v1(i) v2(i)],[v1(i+j) v2(i+j)]);

        n1 = norm([v1(i) v2(i)]);

        n2 = norm([v1(i+j) v2(i+j)]);

        theta(i) = acosd(d/n1/n2);

    end

    mean_theta(j) = mean(theta);

end



function mean_theta = method3(v1,v2,nSteps)

v = [v1,v2];

v = v./sqrt(sum(v.^2,2)); % Can use VECNORM in newest MATLAB

n = numel(v1);

mean_theta = zeros(nSteps,1);

for j = 1:nSteps

    theta = zeros(1,n-j);

    for i = 1:(n-j)

        theta(i) = acosd(dot(v(i,:),v(i+j,:)));

    end

    mean_theta(j) = mean(theta);

end



function mean_theta = method4(v1,v2,nSteps)

v = [v1,v2];

v = v./sqrt(sum(v.^2,2)); % Can use VECNORM in newest MATLAB

n = numel(v1);

mean_theta = zeros(nSteps,1);

for j = 1:nSteps

    i = 1:(n-j);

    theta = acosd(dot(v(i,:),v(i+j,:),2));

    mean_theta(j) = mean(theta);

end