Loop unrolling question

Hi, I'm trying to generate the matrix A usually used for least-square fitting, namely columns of sine and cosine of increasing frequency. The number of desired frequencies (columns) is passed as a parameter. I'm trying to figure out a 'non-for-loop' way to generate the matrix. I managed to vectorize the column length using the time vector, but I can't find how to vectorize the number of columns!

For example,

Given: npoints = 1024; w = 1; fs = 55000 t = 0:(1/fs):((npoints-1)/fs);

Now if ncolumns=1, we want: A=[ones(npoints,1) cos(t'*w/fs) sin(t'*w/fs)];

Now if ncolumns=2, we want A=[ones(npoints,1) cos(t'*w/fs) sin(t'*w/fs) cos(2*t'*w/fs) sin(2*t'*w/fs)];

Now if ncolumns=3, we want A=[ones(npoints,1) cos(t'*w/fs) sin(t'*w/fs) cos(2*t'*w/fs) sin(2*t'*w/fs) cos(3*t'*w/fs) sin(3*t'*w/fs)];

Is there a nice way to do this without a for loop or a bunch of 'if'?

Thanks! David

Answer 1

[Andrei Bobrov]

variant

npoints = 1024;
w = 1;
fs = 55000;
t = (0:npoints-1)'/fs;
p1 = t*w/fs*(1:ncls);
A = [ones(npoints,1) reshape(permute(cat(3,sin(p1),cos(p1)),[1 3 2]),npoints,[])];

variant 2

npoints = 1024;
w = 1;
fs = 55000;
t = (0:npoints-1)'/fs;
p1 = t*w/fs*(1:ncls);

data = sortrows([1:2:2*ncls,2:2:2*ncls;cos(p1),sin(p1)]')';
A = [ones(npoints,1) data(2:end,:)];

variant 3

p = (0:npoints-1)'/fs*w/fs;

for i1 = ncls:-1:1
    ii = i1*2;
    ang = i1*p;
    A1(:,ii+1) = sin(ang);
    A1(:,ii) = cos(ang);
end
A1(:,1) = ones(npoints,1);

comparing the run-time of variants 1 and 3

>> t1=zeros(20,2);
npoints = 1024;
w = 1;
fs = 1;
ncls =1500;
for j1 = 1:20
%variant 3 loop
tic,
p = (0:npoints-1)'/fs*w/fs;
for i1 = ncls:-1:1
    ii = i1*2;
    ang = p*i1;
    A1(:,ii+1) = sin(ang);
    A1(:,ii) = cos(ang);
end
A1(:,1) = ones(npoints,1);
t1(j1,1) = toc;clear p A1
%variant 1 vectorized
tic
p1 = (0:npoints-1)'/fs*w/fs*(1:ncls);
A2 = [ones(npoints,1) reshape(permute(cat(3,cos(p1),sin(p1)),[1 3 2]),npoints,[])];
t1(j1,2) = toc;clear p1 A2
end
runtime = [min(t1);mean(t1);max(t1)];
>> runtime
runtime =
    0.1639    0.1510
    0.1720    0.1720
    0.1868    0.2067

[David]

Woahh...had my head spinning for a while. This works amazing, I just had to add a transpose of the reshaped data.

A = [ones(npoints,1) reshape(permute(cat(3,sin(p1),cos(p1)),[1 3 2]),npoints,[])'];

So in the end is there any speed advantage to this compared to a 'for' loop? It's not the easiest code to understand so unless there is a good reason to do it...

Here was my for loop:

for i=2:2:(2*n+1)

A(:,i) = cos(w*t*(i/2));

A(:,i+1) = sin(w*t*(i/2));

end

Thanks a lot!

[David]

Spoke too fast...this added transpose only worked for the case where 2*ncls is the same as npoints.

Things get twisted around too much, I can't keep up!

[Andrei Bobrov]

more add 'variant 2'

[Andrei Bobrov]

more add 'variant 3' and comparing the run-time

Answer 2

[David]

Hi, Now both variants and the for loop work great. However, I tried calculating the time taken to create the matrix for large value of [t] and ncls using tic/toc. Turned out the for loop was substantially much faster (about 3x)! I'm now confused with loop unrolling!

Answer 3

[Derek O'Connor]

David,

What you call loop unrolling is not. Here is a rolled-up loop:

    for i = 1:n
       a(i) = b(i)*c(i)
    end

Here is the same loop unrolled with a stride of 3:

    for i = 1:3:n
       a(i)   = b(i)*c(i);
       a(i+1) = b(i+1)*c(i+1);
       a(i+2) = b(i+2)*c(i+2);
    end

Loop unrolling may speed up the loop because there are now n/3 rather than n loop control tests. Optimizing compilers may automatically do loop unrolling, even when you don't want it.

I suspect you are interested in vectorization rather than loop unrolling.

Your latest comment suggests you have discovered some of the "joys" of vectorization.

Answer 4

[David]

Hi Derek, I understand what you're saying. I meant vectorization.

But I guess my question still holds. One of the first thing than any Matlab manual says is to use vectorization as much as possible to speed up the program. Doesn't seem to be the case here...

Thanks, David

Answer 5

[Derek O'Connor]

Most of those manuals or books are out of date, or are unaware of the JIT compiler, or do not understand what the JIT compiler does (who does?).

Vectorization is not always a good thing. I have many examples where it both slows down the code and increases the memory.

My advice is to first write your program in the simplest possibe way, using loops. Time it. Vectorize and time it. Compare.

Here is a "tip from the horse's mouth", by Matlab's Doug Hull:

http://www.mathworks.com/matlabcentral/answers/7308-when-to-break-the-matlab-rules

"I have seen some tortured logic in the name of vectorizing code. These cases could easily be replaced with a clear, readable for loop. The coder gained nothing in speed with this, yet lost readability.

Vectorizing is the right thing to do in many cases, taking advantage of MATLAB's matrix centric view of the world. It is not always the right thing though. Since the JIT (Just In Time) compiler was introduced, the former 'for loop penalty' is largely gone. Vectorizing is not always important for speed now."