Prev: Robot design
Next: Exception in thread "AWT-EventQueue-0" java.lang.NullPointerException
From: Quentin CAUDRON on 6 Oct 2009 08:03
Hi everyone,

As an example for an introductory computing class I'm giving, I've written some code to calculate pi using the "darts" idea ( Monte Carlo derivative ).

I've written the following function, which runs in about eight to nine seconds on my system :

function [mypi, time] = pi_approx(iterations)

% Start the clock
tic;

counter = 0;

for i = 1 : iterations
x(1) = rand;
x(2) = rand;

if norm(x) <= 1
counter = counter + 1;
end
end

mypi = counter * 4 / iterations;

time = toc;



However, when I replace the vector variable x with two doubles, x and y, and replace the norm() function with the square root of the sum of x^2 and y^2, the code ( below ) runs about two seconds.

function [mypi, time] = pi_approx2(iterations)

% Start the clock
tic;

counter = 0;

for i = 1 : iterations
x = rand;
y = rand;

if sqrt(x^2 + y^2) <= 1
counter = counter + 1;
end
end

mypi = counter * 4 / iterations;

time = toc;




Any ideas as to why this is so much faster ?

Thank you,
Quentin
From: Sebastiaan on 6 Oct 2009 08:44
"Quentin CAUDRON" <q.caudron(a)warwick.ac.uk> wrote in message <hafblo$fkc$1(a)fred.mathworks.com>...
> Hi everyone,
>
> As an example for an introductory computing class I'm giving, I've written some code to calculate pi using the "darts" idea ( Monte Carlo derivative ).
>
> I've written the following function, which runs in about eight to nine seconds on my system :
>
> function [mypi, time] = pi_approx(iterations)
>
> % Start the clock
> tic;
>
> counter = 0;
>
> for i = 1 : iterations
> x(1) = rand;
> x(2) = rand;
>
> if norm(x) <= 1
> counter = counter + 1;
> end
> end
>
> mypi = counter * 4 / iterations;
>
> time = toc;
>
>
>
> However, when I replace the vector variable x with two doubles, x and y, and replace the norm() function with the square root of the sum of x^2 and y^2, the code ( below ) runs about two seconds.
>
> function [mypi, time] = pi_approx2(iterations)
>
> % Start the clock
> tic;
>
> counter = 0;
>
> for i = 1 : iterations
> x = rand;
> y = rand;
>
> if sqrt(x^2 + y^2) <= 1
> counter = counter + 1;
> end
> end
>
> mypi = counter * 4 / iterations;
>
> time = toc;
>
>
>
>
> Any ideas as to why this is so much faster ?
>
Function overhead. Calling the 'norm' function in stead of 'sqrt' increases time. Some timings from me:

1) Your 1st approach: x(1), x(2), norm
6.77 secs

2) Your 2nd approach: x, y, sqrt
1.98 secs

3) Use x(1) = rand, x(2) = rand, AND sqrt(x(1)^2+x(2)^2)
1.90 secs

So, try to minimize the number of function calls in a loop. Note that you are also calling the rand function in every iteration. Call this once:

4)
x = rand(iterations, 2);
for i = 1 : iterations
if sqrt(x(i,1)^2+x(i,2)^2) <= 1
counter = counter + 1;
end
end
0.87 secs

Note that for a large number of iterations the memory might become a burden, so you might want to split the randomizer in blocks, e.g.
N = 5;
for j=1:N
x = rand(iterations/N, 2);
for i=1:iterations/N
....

Now that I am busy with your code, let's do something more challenging by vectorizing it completely ;)
5)
x = rand(iterations, 2);
y = sqrt(x(:,1).^2+x(:,2).^2);
counter = sum(y<=1);
0.68 secs

6) Or not use y:
x = rand(iterations, 2);
counter = sum(sqrt(x(:,1).^2+x(:,2).^2)<=1);
0.49 secs

7) To avoid the problem of allocation memory for a large number of iterations, split it up in blocks:
block = 1e5;
counter = 0;
for i=1:ceil(iterations/block)
% Capture remainder of blocks
if block*i<=iterations
NoIter = block;
else
NoIter = mod(iterations, block);
end

x = rand(NoIter, 2);
counter = counter + sum(sqrt(x(:,1).^2+x(:,2).^2)<=1);
end


Sorry if I took away the fun of the excersise :s

Sebastiaan

> Thank you,
> Quentin
From: John D'Errico on 6 Oct 2009 08:50
"Quentin CAUDRON" <q.caudron(a)warwick.ac.uk> wrote in message <hafblo$fkc$1(a)fred.mathworks.com>...
> Hi everyone,
>
> As an example for an introductory computing class I'm giving, I've written some code to calculate pi using the "darts" idea ( Monte Carlo derivative ).
>
> I've written the following function, which runs in about eight to nine seconds on my system :
>
> function [mypi, time] = pi_approx(iterations)
>
> % Start the clock
> tic;
>
> counter = 0;
>
> for i = 1 : iterations
> x(1) = rand;
> x(2) = rand;
>
> if norm(x) <= 1
> counter = counter + 1;
> end
> end
>
> mypi = counter * 4 / iterations;
>
> time = toc;
>

First of all, use the capabilities of MATLAB when
you teach it!!!!!!! Don't assign the elements of x
independently, as distinct statements. Teach your
students to work with vectors of elements from
the very beginning! So a looped version might
look like this (after I corrected your code, because
you used the wrong direction for the inequality.)

tic;
counter = 0;
for iter = 1 : iterations
x = rand(1,2);
counter = counter + (norm(x) <= 1);
end
mypi = counter * 4 / iterations;
time = toc;

Even better, I might try this:

tic;
counter = 0;
for i = 1 : iterations
counter = counter + (abs(rand + j*rand) <= 1);
end
mypi = counter * 4 / iterations;
time = toc;

But really, I would NEVER do it using a loop at all!
The MATLAB way might be like this:

tic
count = sum(abs(rand(1,iterations) + j*rand(1,iterations)) <= 1);
mypi = count*4/iterations;
toc

Or, if you prefer, use sqrt.

tic
count = sum(sqrt(rand(1,iterations).^2 + rand(1,iterations).^2) <= 1);
mypi = count*4/iterations;
toc

Note that for 10000 iterations, either of these methods
took a fraction of the time required for your loop.

Use the capability of MATLAB. Don't write MATLAB as
if you are writing C or basic.

John
From: Sebastiaan on 6 Oct 2009 08:59
More is possible: remove the sqrt call: sqrt(x)<=1 is the same as x<=1:

0.35 secs
counter = counter + sum((x(:,1).^2+x(:,2).^2)<=1);

Not bad, coming from 6.77 to 0.35 secs: speedup of 20.
From: Rune Allnor on 6 Oct 2009 09:00
On 6 Okt, 14:50, "John D'Errico" <woodch...(a)rochester.rr.com> wrote:

> First of all, use the capabilities of MATLAB when
> you teach it!!!!!!!

Or rather, decide if you want to teach programming
or if you want to teach matlab.

The first approach about separating the algorithm
into small, easily understood, easily tested separate
functions is what is considered 'best practices' when
programming the compiled languages: The focus is on
the human reader and maintainer of the code, the
compiler is left to optimize away function overheads
and so on.

Matlab is not a compiled language, and thus has a
lot of inherent bottlenecks.

Just about anything that makes matlab go fast, breaks
with what is otherwise considered 'best practice'
programming. Which is a very good reason for *not*
teaching 'the capacities of matlab' for newbie
programmers.

Rune
 |  Next  |  Last
Pages: 1 2
Prev: Robot design
Next: Exception in thread "AWT-EventQueue-0" java.lang.NullPointerException