Exercise: Vectorization of MC Simulation for π

Exercise: Vectorization of MC Simulation for π

1 clear; clc;

2

3 n = 1e5;

4 x = rand(n, 1);

5 y = rand(n, 1);

6 m = sum(x .ˆ 2 + y .ˆ 2 < 1);

7 result = 4 * m / n

• More clear and faster!!!

while Loops

• The whileloops are used to repeat the instructionsuntil the continuation criterion is not satisfied.

1 while criterion

2 % body

3 end

• Be aware that theif statement executes only once; you should use the whileloop if you want to repeat some actions.

Zheng-Liang Lu 82

Example: Compounding

• Let balance be the initial amount of some investment, and r be the annualized return rate.

• Write a program which calculates the holding years when this investment doubles it value.

Zheng-Liang Lu 84

Solution

• In this case, we don’t know how many iterations we need before the loop.

1 clear; clc;

2

3 balance = 100;

4 r = 0.01;

5 goal = 200;

6

7 holding years = 0;

8 while balance < goal

9 balance = balance * (1 + r);

10 holding years = holding years + 1;

11 end

12 holding years

• Note that the criterion is evaluated to continue the loop.

Infinite Loops

1 while true

2 disp("Press ctrl+c to stop me!!!");

3 end

• Note that your program can terminate the program by pressing ctrl+c.

Zheng-Liang Lu 86

More Exercises (Optional)

• Let a > b be two any positive integers.

• Write a program which calculates the remainder of a divided by b.

• Do not use mod(a, b).

• Write a program which determines the greatest common divisor (GCD) of a and b.

• Do not use gcd(a, b).

Numerical Example: Bisection Method for Root-Finding

Zheng-Liang Lu 88

Problem Formulation

Input

- Target function f (x ) = x³− x − 2.

- Initial search interval [a, b] = [1, 2].

- Error tolerance  = 1e − 9.

Output

- The approximate root ˆr .

Solution

1 clear; clc;

2

3 a = 1; b = 2; eps = 1e-9;

4

5 while b - a > eps

6

7 c = (a + b) / 2;

8 fa = a * a * a - a - 2;

9 fc = c * c * c - c - 2;

10

11 if fa * fc < 0

12 b = c;

13 else

14 a = c;

15 end

16 17 end

18 root = c

19 residual = fc

Zheng-Liang Lu 90

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -8

-6 -4 -2 0 2 4

x³ - x - 2 Root

c = 1.52137970691547

“All science is dominated by the idea of approximation.”

– Bertrand Russell(1872–1970)

Zheng-Liang Lu 92

Jump Statements

• Abreakstatement terminates afor orwhileloop immediately.

• Akaearly termination.

• A continuestatement skips instructions behind it and start the next iteration.

• Directly jump to the very beginning of the loop; still in the loop.

• Notice that the breakandcontinuestatements must be conditional.

Example: Primality Test

• Let x be any positive integer larger than 2 as input.

• Then x is aprimenumber if ∀y ∈ {2, 3, . . . , x − 1}, y is not a divisor of x , denoted by y - x.

• In other words, x is called a compositenumber if

∃y ∈ {2, 3, . . . , x − 1}, y | x.

• Now write a program which determines if x is a prime number.

1Also see Manindra Agrawal, Neeraj Kayal, Nitin Saxena (2002).

Zheng-Liang Lu 94

1 clear; clc;

2

3 x = input('Enter x > 2? ');

4 isPrime = true; % a flag, true if the number is prime

5 for y = 2 : sqrt(x)

6 if mod(x, y) == 0

7 isPrime = false;

8 break;

9 end

10 end

11

12 if isPrime

13 disp([num2str(x) ' is a prime number.']);

14 else

15 disp([num2str(x) ' is a composite number.']);

16 end

Equivalence: for and while Loops

• Whatever you can do with afor loop can be done with awhile loop, and vice versa.

1 clear; clc;

2

3 balance = 100; goal = 200; r = 0.01;

4

5 for years = 1 : inf % inf: a huge but finite integer

6 balance = balance * (1 + r);

7 if balance >= goal

8 break;

9 end

10 end

11 years

Zheng-Liang Lu 96

• For another example,

1 clear; clc;

2

3 x = input("Enter x > 2? ");

4

5 isPrime = true; y = 2;

6 while isPrime && y < x

7 isPrime = mod(x, y);

8 y = y + 1;

9 end

10

11 if isPrime

12 disp(num2str(x) + " is a prime number.");

13 else

14 disp(num2str(x) + " is a composite number.");

15 end

Nested Loops

• Write a program which outputs the following patterns:

Zheng-Liang Lu 98

• You may use fprintf(”*”) and fprintf(”\n”) to print a single star and break a new line, respectively.

1 clear; clc;

2

3 % case (a)

4 for i = 1 : 5

5 for j = 1 : i

6 fprintf("*");

7 end

8 fprintf("\n");

9 end

Exercise: e ∼ 2.7183

• Write a program to estimate the Euler constant by Monte Carlo simulation.

• It can be done as follows.

• Let N be the number of iterations.

• For each iteration, find the minimal number n so that Pn

i =1r_i > 1 where r_i is the random variable following the standard uniform distribution (you can simply use rand).

• Then e is the average of n.

Zheng-Liang Lu 100

Special Issue: Sort

1 >> stocks = {"GOOG", 15;

2 "TSMC", 12;

3 "AAPL", 18};

4 >> [~, idx] = sort([stocks{:, 2}], "descend")

5

6 idx =

7

8 3 1 2

9

10 >> stocks = stocks(idx, :)

11

12 stocks =

13

14 "AAPL" [18]

15 "GOOG" [15]

16 "TSMC" [12]

Programming Exercise: Sorting Algorithm

• Let A be any array.

• Write a program which outputs the sorted array of A (in ascending order).

• For example, A = [5, 4, 1, 2, 3].

• Then the sorted array is [1, 2, 3, 4, 5].

2See https://visualgo.net/sorting.

Zheng-Liang Lu 102

Special Issue: Random Permutation

• Use randperm to generate an index array with arandom order.

1 >> A = ["Matlab", "Python", "Java", "C++"];

2 >> idx = randperm(length(A))

3

4 idx =

5

6 3 1 2 4

7

8 >> A(idx)

9

10 ans =

11

12 1x4 string array

13

14 "Java" "Matlab" "Python" "C++"

“Exploring the unknown requires tolerating uncertainty.”

– Brian Greene

“I can live with doubt, and uncertainty, and not knowing.

I think it is much more interesting to live not knowing than have answers which might be wrong.”

– Richard Feynman

Zheng-Liang Lu 104

Speedup: Vectorization (Revisited)

• Vector in, vector out.

1 >> x = randi(100, 1, 5)

2

3 x =

4

5 88 30 90 73 82

6

7 >> dx = diff(x)

8 9 dx =

10

11 -58 60 -17 9

Advantages from Vectorization

• Appearance: vectorized mathematical code appears more like the mathematical expressions found in textbooks, making the code easier to understand.

• Less error prone: without loops, vectorized code is often shorter.

• Fewer lines of code mean fewer opportunities to introduce programming errors.

• Performance: vectorized code often runs much fasterthan the corresponding code containing loops.

Zheng-Liang Lu 106

Performance Analysis: Profiling

• Use a timerto measure your performance.⁴

• In newer version, press the button Run and Time.

• Identify which functions are consuming the most time.

• Know why you are calling them and then look for alternatives to improve the overall performance.

4Note that the results may differ depending on the difference of run-time environments, so make sure that you benchmark the algorithms on thesame

tic & toc

• The command tic makes a stopwatch timer start.

• The command toc returns the elapsed time from the stopwatch timer started by tic.

1 >> tic

2 >> toc

3 Elapsed time is 0.786635 seconds.

4 >> toc

5 Elapsed time is 1.609685 seconds.

6 >> toc

7 Elapsed time is 2.417677 seconds.

Zheng-Liang Lu 108

Selected Performance Suggestions

• Preallocate arrays.

• Instead of continuously resizing arrays, consider preallocating the maximum amount of space required for an array.

• Vectorizeyour code.

• Create new variables if data type changes.

• Use functions instead of scripts.

• Avoid overloading Matlab built-in functions.

Programming Exercise: A Benchmark

• Let N = 1e1, 1e2, 1e3, 1e4, 1e5.

• Write a program which produces a benchmark for the following three cases:

• Generate an array of 1 : N by dynamically resizing the array.

• Generate an array of 1 : N by allocating an array of size N and filling up sequentially.

• Generate an array of 1 : N by vectorization.

Zheng-Liang Lu 110

Analysis of Algorithms (Optional)

• For one problem, there exist various algorithms (solutions).

• We then compare these algorithms for various considerations and choose the most appropriate one.

• In general, we want efficient algorithms.

• Except for real-time performance analysis, could we predict before the program is completed?

• Definitely yes.

Growth Rate

• Now we use f (n) to denote thegrowth rate of time cost as a function of n.

• In general, n refers to the data size.

• For simplicity, assume that every instruction (e.g. + − ×÷) takes 1 unit of computation time.

• Find f (n) for the following problem.

• Sum(n): ?

• Triangle(n): ?

Zheng-Liang Lu 112

O-notation

• In math, O-notation describes the limiting behaviorof a function, usually in terms of simple functions.

• We say that

f (n) ∈ O(g (n)) as n → ∞ if and only if ∃c > 0, n₀ > 0 such that

| f (n) | ≤ c| g (n) | ∀n ≥ n₀.

• So O(g (n)) is a collection featured by a simple function g (n).

• We use f (n) ∈ O(g (n)) to denote that f (n) is one instance of O(g (n)).

• Big-O is used for the asymptotic upper boundof time complexity of algorithm.

• In layman’s term, Big-O describes the worstcase of this algorithm.

Zheng-Liang Lu 114

• For example, 8n²− 3n + 4 ∈ O(n²).

• For large n, you could ignore the last two terms. (Why?)

• It is easy to find a constant c > 0 so that cn²> 8n², say c = 9.

• Hence the statement is proved.

• Also, 8n²− 3n + 4 ∈ O(n³) but we seldom say this. (Why?)

• However, 8n²− 3n + 4 /∈ O(n). (Why?)

• What is this analysis related to the algorithm?

• Any insight?

Common Simple Functions

7See Table 4.1 and Figure 4.2 in Goodrich and etc, p. 161.

Zheng-Liang Lu 116

Remarks

• We often make a trade-offbetween time and space.

• Unlike time, we can reuse memory.

• Users are sensitive to time.

• Playing game well is hard.⁸

• Solve the problem P ?= NP, which is one of Millennium Prize Problems.⁹

8See https://en.wikipedia.org/wiki/Game_complexity.

“All roads lead to Rome.”

– Anonymous

“但如你根本並無招式，敵人如何來破你的招式？”

–風清揚。笑傲江湖。第十回。傳劍

Zheng-Liang Lu 118