tangent tangent correlation calculation in matrix form MATLAB

Question

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?

Do you have the Econometrics toolbox? If so you can have a look at `autocorr`, which calculates the auto correlation. You should be able to adapt it to calculate the tangent tangent correlation instead. When I say adapt, I mean create a new similar function, not changing the existing one. — Nicky Mattsson, Nov 26 '18 at 13:28
For starters, growing arrays, like you do with `theta` is very, very, very bad for MATLAB performance. Just preallocate it to the correct size (that's known in this case, although variable per iteration) and you'll save a lot of time already. — Adriaan, Nov 26 '18 at 13:28
@NickyMattsson, I'd like to stick with what matlab has to offer without the toolboxes. — jarhead, Nov 26 '18 at 13:54
@Adriaan, if you can post an answer with specifying how my example is bad for matlab and how you propose to write it differently, it would be great and will also benefit others. — jarhead, Nov 26 '18 at 13:55
@jarhead see [this question](https://stackoverflow.com/q/48351041/5211833). I'm not going to write an answer, as preallocating is too evident imo, as it should always be done. No value there. — Adriaan, Nov 26 '18 at 15:20
Simply write `theta=zeros(n-j,1);` where you now have `theta=[];`, and then do `theta(i) = acosd(d/n1/n2);` That should speed up your code significantly. This is indeed [standard practice](https://www.mathworks.com/help/matlab/matlab_prog/preallocating-arrays.html). — Cris Luengo, Nov 26 '18 at 15:55

score 2 · Answer 1 · answered Nov 26 '18 at 16:32

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: %g\n',max(abs(mean_theta1(:)-mean_theta2(:))));
mean_theta3 = method3(v1,v2,nSteps);
fprintf('diff method1 vs method3: %g\n',max(abs(mean_theta1(:)-mean_theta3(:))));
mean_theta4 = method4(v1,v2,nSteps);
fprintf('diff method1 vs method4: %g\n',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

score 1 · Accepted Answer · answered Nov 26 '18 at 19:12

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

tangent tangent correlation calculation in matrix form MATLAB

2 Answers2