vectorization of for loop code

prod1=zeros(1000,1000);

for n=1:1000;

prod1(n,(1:n))=temp*((n-(1:n)+1).^alpha-(n-(1:n)).^alpha);

end

[KSSV]

First tell me is that code working?

[Pooja Patil]

yes.This code creats lower triangular matrix of 1000*1000.

Answer 1

[Andrei Bobrov]

nn = 1:n;
z = temp*(nn.^alpha - (nn - 1).^alpha);
prod1 = tril(toeplitz(z));

[Stephen23]

+1 nice and tidy

Answer 2

[Jan]

I assume you want to vectorize this to gain more speed. Then start with simplifying the loop at first:

prod1 = zeros(m, m);   % Clean loop
v     = temp * diff((0:m) .^ alpha);  % Reduce number of calls to ^
for n = 1:m
  prod1(n:m, n) = v(1:m-n+1);
end

Some timings:

function SpeedTest
temp  = 0.2;  % Guessed parameters
alpha = 0.3;
m     = 1000;

tic;
for k = 1:10
  prod1 = zeros(m, m);   % Original loop
  for n = 1:m
    prod1(n,(1:n))=temp*((n-(1:n)+1).^alpha-(n-(1:n)).^alpha);
  end 
end
toc

tic;
for k = 1:10
  prod1 = zeros(m, m);   % Clean loop
  v     = temp * diff((0:m) .^ alpha);  % Reduce number of calls to ^
  for n = 1:m
    prod1(n:m, n) = v(1:m-n+1);
  end
end
toc

tic;
for k = 1:10
  nn = 1:n;              % Vectorized
  z = temp*(nn.^alpha - (nn - 1).^alpha);
  prod1 = tril(toeplitz(z));
end
toc

 Elapsed time is 1.687053 seconds.  % Original
 Elapsed time is 0.070619 seconds.  % Cleaned loop
 Elapsed time is 0.409178 seconds.  % Vectorized with TRIL/TOEPLITZ

If run time matters, using a clean loop can be more important then vectorization. The main rules for optimizing loops:

1. Avoid repeated calculations
2. Process array in column order
3. The power operation is expensive. Reduce the calls whenever possible.

Finally: Even if the lean loop is faster, Andrei's suggestion (with the leaner creation of the data) is nicer:

prod1 = tril(toeplitz(temp * diff((0:m) .^ alpha)));

The loop can save seconds of run time, the vectorized version minutes during debugging or expanding the code.

[Andrei Bobrov]

+1. Well!

cc = (0:n).^alpha;
k = temp*(cc(2:n+1) - cc(1:n));
prod1 = zeros(n);
for ii = 1:n
    prod1(ii:n,ii) = k(1:n - ii + 1);
end