matlab simplification of code

Last week I showed my MatLab code to my major professor and when he saw the following loop he went crazy and went into a 15 min long lecture about how I am not using the full capabilities of MatLab and that I’m creating ineffective and computationally heavy code, etc. Anyhow, what I got out of the lecture is that there is a better way of writing this loop (even though it works as it is). Can you help me figure out what that is? I am trying to find:

Thanks Thoroughly rebuked.

>> a = [1:1000];
>> s=0;
>> for n = 1:1000
     b(n) = 1000-a(n);
     s = s + a(n)*b(n);
   end

[Jan]

I've formatted your code to increase the readablility.

Answer 1

[Joseph Cheng]

There is no need to perform the for loop and this problem can be solved through matrix arithmetic. to simplify it lets say

 a = [1 2 3];

if we go

 b= 1000 - a;

b will return as [999 998 997]

matlab will perform the matrix subtraction and there is no need to perform the each index individually.

the next thing is you're computing s which is the sum of the element by element multiplication of a and b. the element by element multiplication can be done using the .* notation or the function times().

after that you can perform the sum using the function sum(). which if you read the documentation on sum() you'll see it performs the summation of an array you give it.

in the end using what matlab can give you is

 a = 1:1000
 b = 1000 -a; 
 s = sum(a.*b);

[Stephen23]

Hers is a timing comparison using tic and toc, for one hundred repetitions:

Elapsed time is 5.39954 seconds.
Elapsed time is 0.00300002 seconds.

Joseph Cheng's solution is more than a thousand times faster than using the original slow loop. MATLAB is not a low-level programming language, so why limit it to low-level constructions (like loops)?

http://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

http://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html

[Thorsten]

Hm. On my machine it is just about 10 times faster. And if you pre-allocate b it is only twice as fast, at best.

Another way would be to write

a = 1:1000; 
s = a*(1000-a');

So here are the results for my 100 repetitions:

Elapsed time is 0.010911 seconds. % original
Elapsed time is 0.001583 seconds. % original with b preallocated
Elapsed time is 0.000954 seconds. % Chen
Elapsed time is 0.000818 seconds. % my solution

Generated with

clear
tic
for i = 1:100
   a = 1:1000;
   s= 0; for n = 1:1000, 
   b(n) = 1000-a(n); s = s + a(n)*b(n); end, 
end
toc

clear
tic
for i = 1:100
   a = 1:1000;
   b = zeros(1, 1000);
   s= 0; for n = 1:1000, 
   b(n) = 1000-a(n); s = s + a(n)*b(n); end, 
end
toc

clear
tic
for i = 1:100,
  a = 1:1000;
  b = 1000 - a; 
  s1 = sum(a.*b);
end
toc

clear 
tic, for i = 1:100, 
  x = 1:1000; 
  s2 = x*(1000-x'); 
end
toc

[Jan]

The code is not fair: The "original" version has a pre-allocated b after the first iteration. I've added a lean loop method using the index instead of creating a vector.

function test
iter = 1000;

tic
for i = 1:iter
   clear('b');
   a = 1:1000;
   s = 0;
   for n = 1:1000,
      b(n) = 1000-a(n);
      s = s + a(n)*b(n);
   end
end
toc

tic
for i = 1:iter
   a = 1:1000;
   b = zeros(1, 1000);
   s = 0;
   for n = 1:1000,
      b(n) = 1000-a(n);
      s = s + a(n)*b(n);
   end
end
toc

tic
for i = 1:iter
   s = 0;
   for n = 1:1000,
      s = s + n * (1000 - n);
   end
end
toc

tic
for i = 1:iter
   a = 1:1000;
   b = 1000 - a;
   s1 = sum(a.*b);
end
toc

tic
for i = 1:iter
   x = 1:1000;
   s2 = x*(1000-x');
end
toc

 Matlab 2011b/64/Win:
Elapsed time is 0.438072 seconds.  % Original
Elapsed time is 0.016922 seconds.  % pre-allocated 'b'
Elapsed time is 0.006091 seconds.  % direct index usage
Elapsed time is 0.009741 seconds.  % SUM
Elapsed time is 0.007666 seconds.  % BLAS

Note that Matlab's JIT acceleration requires, that a line does not contain more than one command. And the JIT has more power in functions than in scripts.

My conclusion: The benefit of pre-allocation is obvious even for this very artificial example. The power of Matlab's vectorized commands can be measured here also, but in addition the SUM and the BLAS methods are nicer: You have less chances for typos.

By the way: b(1000) is zero, so stop the sum at n=999 ;-)

Answer 2

[Jan]

Your code:

a = 1:1000;  % Square brackets not required: This is a vector already
s = 0;
for n = 1:1000
   b(n) = 1000-a(n);  % b is growing iteratively
   s = s + a(n)*b(n);
 end

The iterative growing of b is a standard mistake. In the first iteration, Matlab reserves 8 bytes for a scalar double (plus about 100 bytes overhead for the internat organization). In the second iteration Matlab reserves 16 bytes for a [1 x 2] double vector and copies the first double from the former variable. The the formerly used memory is released. And so on.

Finally Matlab (or the operating system to be exact) did not allocate memory for 1000 double values, but sum(1:1000) doubles, which means 4'004'000 bytes and almost the same numer of bytes have been copied in addition. This is expensive! Much more expensive the additions and the multiplications.

The solution is "pre-allocation". Create b as a [1 x 1000] array initially:

b = zeros(1, 1000);

Another solution: If b is not required, storing it in an array should be avoided, and the same for a:

s = 0;
for n = 1:1000
  s = s + n * (1000 - n);
end

If the index n equals the contents of a(n), why not using the index directly?