vectorizing nested loops

I'm wondering if anyone can suggest a way of vectorizing these nested loops:

for i = 1:n
  for j = 1:n
      for k = 1:n
          for l = 1:n
              if (i+k)<=n && (j+l)<=n
                  r(i,j, k, l) = a(i, j)*a(k, l)*a(i+k, j+l);
              end
          end
      end
  end
end

Thanks a lot!

Answer 1

[Kye Taylor]

Sorry, but I don't have a suggestion to vectorize this thing. However, as the following code demonstrates, it looks like you could get away with only computing half of the values as n gets large :)

Someone explain these curves to me... I guess the blue curve has something to do with volumes of the corners of hypercubes in high dimensions. The red curve is due to symmetry in the assignment?

I know I'm not preallocating... that's just to avoid questions about me getting the size wrong.

 ns = 3:40;
nonZeros = zeros(1,length(ns));
uniqueVals = zeros(1,length(ns));
numElements = zeros(1,length(ns));
for n = ns;
a = rand(n,n,n,n);
for i = 1:n
  for j = 1:n
      for k = 1:n
          for l = 1:n
              if (i+k)<=n && (j+l)<=n
                  r(i,j, k, l) = a(i, j)*a(k, l)*a(i+k, j+l);
              end
          end
      end
  end
end 
nonZeros(n-min(ns)+1) = nnz(r(:));
uniqueVals(n-min(ns)+1) = length(unique(r(:)));
numElements(n-min(ns)+1) = numel(r(:));
end
figure,plot(ns,nonZeros./numElements,'.',ns,uniqueVals./nonZeros,'r.')
xlabel(n)
legend('Density of nonzeros in r','Ratio: unique to nonzeros');

[Max]

Thanks, Kye. This is really interesting and helpful.

Answer 2

[Teja Muppirala]

Vectorization would probably involve calling NDGRID and then doing a bunch of logical indexing. It would take a lot more memory.

But as a small step you can rewrite the indices a little better, and get rid of the IF statement:

for i = 1:n-1
    for k = 1:n-i
        for j = 1:n-1
            for l = 1:n-j
                r(i,j, k, l) = a(i, j)*a(k, l)*a(i+k, j+l);
            end
        end
    end
end

[Max]

It does save quite a bit of time. I got a factor of about 5 which makes a big difference for me. Thanks, Teja.

Answer 3

[Robert Cumming]

preallocate your variable r

That will improve the performance.