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');
