cos() between vectors in matrix Optimization

Hello all,

I'm writing to ask for some help optimizing the following code. I'm taking a matrix, S, into my function, and I want to collect the cosine of the angle between each pair of columns in S into Dc, as you can see, I have naively done this with two for loops, what I want to know is if there is a way of doing this that takes advantage of the matrix operations of MATLAB, so that it runs faster. I am currently testing on a small(ish) matrix, Dc ends up on the order of 1300x1300, but I hope to expand the code to operate on larger matrices, between 10 and 15 times this size.

Any pointers would be much appreciated.

Dc = zeros(size(S,2));
for i = 1:length(S)
    for j = 1:length(S)
        Dc(i,j) = (S(:,i)' * S(:,j))/(norm(S(:,i))*norm(S(:,j)));
    end
end

Answer 1

[Roger Stafford]

This might be faster. It avoids repetition in computing norms.

 Dc = S'*S;
 N = sqrt(diag(Dc));
 Dc = Dc./(N*N');

[Matthew]

I appreciate the response, but your code doesn't do the same thing that mine does.

Run them on a 3x3 S matrix, and you'll see that the Dc that are created are not the same.

[Roger Stafford]

I have done that and get the same answer for both methods within round-off error. What result would you get on the following?

 format long
 S = magic(3);

 Dc = zeros(size(S,2));
 for i = 1:length(S)
   for j = 1:length(S)
     Dc(i,j) = (S(:,i)'*S(:,j))/(norm(S(:,i))*norm(S(:,j)));
   end
 end
 Dc1 = Dc;

 Dc = S'*S;
 N  = sqrt(diag(Dc));
 Dc = Dc./(N*N');
 Dc2 = Dc;

 Dc1 =
   1.00000000000000   0.60459579105879   0.86516853932584
   0.60459579105879   1.00000000000000   0.60459579105879
   0.86516853932584   0.60459579105879   1.00000000000000

 Dc2 =
   1.00000000000000   0.60459579105879   0.86516853932584
   0.60459579105879   1.00000000000000   0.60459579105879
   0.86516853932584   0.60459579105879   1.00000000000000

I repeated the test for S = randn(25,32) and again the maximum difference between the two results was 1.332267629550188e-15, well within expected round-off error.

There is one way you can get significantly different results. You have written "for i = 1:length(S)" instead of "for i = 1:size(S,2)" and similarly for the j-loop, which will give you an error message if there are more rows than columns in S. Remember, 'length' is the maximum of the number of rows and of the number of columns.

(Corrected)

[Matt J]

Using DNorm2, you can eliminate the outer product N*N',

 Dc = bsxfun(@rdivide,S, DNorm2(S,1));
 Dc = Dc'*Dc;

[Matthew]

Thank you Rodger.

The size v length was my main issue. I was testing with some square matricies, and some rectangular matricies.