Question? Leave a message!




Stochastic algorithms

Stochastic algorithms
ECE 250 Algorithms and Data Structures Stochastic algorithms Douglas Wilhelm Harder, M.Math. LEL Department of Electrical and Computer Engineering University of Waterloo Waterloo, Ontario, Canada ece.uwaterloo.ca dwharderalumni.uwaterloo.ca © 20062013 by Douglas Wilhelm Harder. Some rights reserved.Stochastic algorithms 2 Outline In this topic, we will cover: – Concepts about random numbers – Random number generators – Monte Carlo methods • Monte Carlo integration – Randomized quicksort – Skip listsStochastic algorithms 3 Stochastic The word stochastic comes from Greek: στόχος, – The relative translation is a word describing a target stick for archery practices – The pattern of arrow hits represents a stochastic process http://marshfieldrodandgunclub.com/ArcheryStochastic algorithms 4 Random Numbers I flipped a coin 30 times: HTHTHHHTTTHHHTHTTTTTTHHHHHHTTH One may almost suggest such a sequence is not random: what is the probability of getting 6 tails and then 6 heads – Unfortunately, this is exactly what I got with the first attempt to generate such a sequence, and is as random a sequence as I can generate – Humans see patterns where none are…Stochastic algorithms 5 Random Numbers Generating random numbers is difficult… A small radioactive source is believed to have been used to generate the Rand Corporation's book A Million Random Digits with 100,000 Normal Deviates – Called an electronic “roulette wheel” • uses a ―random frequency pulse source‖ – 20,000 lines with 50 decimal digits per line divided into five columns of ten digits each – Available at the low price of USD 90 http://www.wps.com/projects/million/index.htmlStochastic algorithms 6 Random Numbers Problem: humans will not open the book to a random page... 73735 45963 78134 63873 02965 58303 90708 20025 The instructions include: 98859 23851 27965 62394 ―...open the book to an unselected page... 33666 62570 64775 78428 blindly choose a 5digit number; this 81666 26440 20422 05720 number reduced modulo 2 determines the 15838 47174 76866 14330 starting line; the two digits to the right... 89793 34378 08730 56522 determine the starting column...‖ 78155 22466 81978 57323 16381 66207 11698 99314 75002 80827 53867 37797 Column: 81978 / 20000 = 4 99982 27601 62686 44711 84543 87442 50033 14021 Line: 81978 20000 = 01978 77757 54043 46176 42391 80871 32792 87989 72248 30500 28220 12444 71840Stochastic algorithms 7 Random Numbers Strategy 0: – Roll a die Not very useful, especially for computers http://xkcd.com/221/Stochastic algorithms 8 Random Numbers Strategy 1: – Use the system clock This is very useful for picking one 5decimaldigit number about once a day It is also very poor for picking many random numbers, as any program doing so will most likely do so periodicallyStochastic algorithms 9 Random Numbers Strategy 2: – Use the keyboard Certain encryption programs use intervals between keystrokes to generate random numbers The user is asked to generate some text which is then analyzed – Random, but slow and tediousStochastic algorithms 10 Random Numbers Our best hope is to generate (using mathematics) a sequence of numbers which appears to be random A sequence of numbers which appears to be random but is not is said to be pseudorandom – We will look at one random number generatorStochastic algorithms 11 Linear Congruential Generators The first pseudorandom number generator which we will look as described by Lehmer in 1951 It is called a linear congruential generator because: – it uses a linear multiplier – it produces objects of the same shapeStochastic algorithms 12 Linear Congruential Generators Select an initial value x (called the seed) 0 Choose two fixed values A 0 and M 1 Given x , calculate i x = Ax mod M i + 1 i There are some restrictions on these values: – M should be prime – A and x should both be relatively prime to M 0Stochastic algorithms 13 Linear Congruential Generators For example, if M = 17, A = 6, and x = 3, the first 17 numbers are: 0 3, 1, 6, 2, 12, 4, 7, 8, 14, 16, 11, 15, 5, 13, 10, 9, 3 th The 17 number equals the first, so after this, the sequence will repeat itself A different seed generates a different sequence, e.g., with x = 4 we 0 get 4, 7, 8, 14, 16, 11, 15, 5, 13, 10, 9, 3, 1, 6, 2, 12, 4 however, it is only a shift...Stochastic algorithms 14 Linear Congruential Generators The choice of A will affect quality of the sequence Given the previous seed and modulus, but if we choose A = 2, we get the sequence 3, 6, 12, 7, 14, 11, 5, 10, 3, 6, 12, 7, 14, 11, 5, 10, 3 which you may note repeats itself earlier than when we used A = 2Stochastic algorithms 15 Linear Congruential Generators 31 If we choose the prime M = 2 – 1, then one value of A recommended by Lehmer is A = 48271 31 This generates a sequence of 2 – 2 unique numbers before it repeats We can then use modulus to map this onto a smaller range, or divisionStochastic algorithms 16 Linear Congruential Generators A word on linear congruential generators and the like: ―Any one who considers arithmetical methods of producing random digits is, of course, in a state of sin...‖ John von Neumann The point here is, these methods produce pseudorandom numbers and not actually random numbersStochastic algorithms 17 Linear Congruential Generators The C++ standard library comes equipped with the int rand() function – Don’t use it... Most C++ standard libraries come with: Function Name Range 31 long lrand48() 0, 1, ..., 2 – 1 31 31 long mrand48() –2 , ..., 2 – 1 double drand48() 0, 1)Stochastic algorithms 18 Linear Congruential Generators By default, these use a 48bit seed with x = 0 and 0 x = (Ax + C) mod M k + 1 k where A = 25214903917 10111011110111011001110011001101101 2 C = 11 1011 2 48 M = 2 and where the most significant bits of x are used k + 1Stochastic algorithms 19 Linear Congruential Generators All three generates use the same value – Suppose that x is the value is the following 48 bits: k + 1Stochastic algorithms 20 Linear Congruential Generators lrand48() uses the most significant 31 bits – Recall the avalanche effectStochastic algorithms 21 Linear Congruential Generators mrand48() uses the most significant 32 bits st – This is a simple cast: the 1 bit is the sign bit – This number is positive—use 2s complement to get the valueStochastic algorithms 22 Linear Congruential Generators drand48() uses all 48 bits and treats it as a fractional value in 0, 1) – The mantissa stores 52 bits with an implied leading 1. – The last five bits will be 0Stochastic algorithms 23 Linear Congruential Generators These can be easily mapped onto different ranges: long lrand48ab( long a, long b ) return a + lrand48() (b – a + 1); double drand48ab( double a, double b ) return a + drand48()(b – a); Stochastic algorithms 24 Other PseudoRandom Number Generators Other random number generators: Nonlinear additive feedback long random() RANDMAX varies from platform to platform Matrix linear recurrence over a finite binary field F 2 Mersenne Twister Excellent for simulations, poor for cryptography, long periods By Matsumoto and Nishimura in 1997 Slow but good for cryptography 2 BlumBlumShub x = x mod pq where p and q are large prime numbers k + 1 k By Blum, Blum and Shub in 1986 In 1995, George Marsaglia published his diehard battery of statistical tests for pseudorandom number generatorsStochastic algorithms 25 Linear Congruential Generators The course text discusses a few further features which are of note: – Do not attempt to come up with your own random number generators: use what is in the literature – Using the given value of M 00111111111111111111111111111111 = 3FFF 2 16 calculating Ax may cause an overflow i – It is possible to program Ax mod M to avoid this iStochastic algorithms 26 Normal Distributions If you take 1000 similar inverters and measure their delays, you will find that the delay are not all the same: – All will be close to some expected value (the mean m) – Most will be close to m while some will be further away David Rennie, ECE 437Stochastic algorithms 27 Normal Distributions The spread is defined by the standard deviation s – If the mean and standard deviation of a manufacturing process is known, then • 99.7 will fall within 3s of the mean • 99.993 will fall within 4s of the mean Many natural systems follow such distributions – The variability in a resistor or capacitor, etc. David Rennie, ECE 437Stochastic algorithms 28 Normal Distributions A standard normal has a mean m = 0 and standard deviation s = 1 David Rennie, ECE 437Stochastic algorithms 29 Normal Distributions There are more efficient algorithms: – The Marsaglia polar method is efficient, is to choose two pseudo random numbers u and u from the uniform distribution on (–1, 1) and if 1 2 2 2 s = u + u 1, calculate 1 2 2ln s   xu = 11 s 2ln s   xu = 22 s both of which are independent (one gives no information about the other) and normally distributed – If these functions are not available, calculate twelve pseudorandom numbers from the uniform distribution on 0, 1 and subtract 6 • The distribution approximates a normal distribution via the central limit theoremStochastic algorithms 30 Normal Distributions Recall the Greek origin of the word stochastic – Here are two archery targets with ten hits generated using the Marsaglia polar method for ( int i = 0; i 10; ++i ) double x, y, s; do x = 2.0drand48() 1; y = 2.0drand48() 1; s = xx + yy; while ( s = 1.0 ); s = std::sqrt( 2.0std::log(s)/s ); // log(x) is the natural logarithm cout "(" xs ", " ys ")" std::endl; Stochastic algorithms 31 Monte Carlo Methods Consider the following operations: – integration, equality testing, simulations In some cases, it may be impossible to accurately calculate correct answers – integrals in higher dimensions over irregular regions – simulating the behaviour of the arrival of clients in a clientserver modelStochastic algorithms 32 Monte Carlo Methods An algorithm which attempts to approximate an answer using repeated sampling is said to be a Monte Carlo method Conor Ogle http://www.gouv.mc/Stochastic algorithms 33 Monte Carlo Methods Consider the problem of determining if p is a prime number 2 – If p is not prime, then given a random number 1 x ≤ 2 ln (p), the Miller Rabin primality test has at least a 75 chance of showing that p is not prime n n n • Uses Fermat’s little theorem (not x + y = z ) – Therefore, test n different integers, and the probability of not being able n to show that a composite p is not prime is 1/4Stochastic algorithms 34 Monte Carlo Integration We will look at a technique for estimating an integral on irregular region V, i.e., f (vv )d  V n Here, V may be a region in R where the dimension could be reasonably high – It may either be unreasonable or impossible to perform an exact approximation using calculus – If the dimension is too high, it may be unreasonable to even us a uniform grid of pointsStochastic algorithms 35 Monte Carlo Integration To solve this, the average value of a function is b 1 f = f (x)d x ab ,  ba  a Solving for the integral, we get: b f (x)d x= f b a   ab ,  a Thus, if we can approximate the average value f , we can ab , approximate the integralStochastic algorithms 36 Monte Carlo Integration Similarly, suppose we wish to integrate a scalarvalued function f(v) over a region – The formula for the average value is the same: 1 ff = (vv )d V  V V – Solving for the integral, we get f (vv )d = f V V  V f – If we can approximate both the average value and the volume of the V region V, we can approximate the integralStochastic algorithms 37 Monte Carlo Integration First, we will look at the easier problem of determining the volume of this region by estimating: 1d v  V n Suppose we know that V⊆ 0, 1 = S for some dimension n – Let V be the volume of the region In this case, given our previous algorithm, we can easily pick a T random point in S by creating a vector v = (x , ..., x ) where each 1 n x∈ 0, 1 k – Create M random vectors in S, and count the number m that fall inside VStochastic algorithms 38 Monte Carlo Integration If the distribution of the vectors v∈ S is random, then the proportion m/M of these vectors in V will be approximately the ratio V/S, that is, the volume of V over the volume of S m VS = Therefore MStochastic algorithms 39 Monte Carlo Integration For example, let us find the area of a circle – Choose 1000 points in 0, 1 × 0, 1:Stochastic algorithms 40 Monte Carlo Integration If we only consider those points which are less than 0.5 distance from the center, we get:Stochastic algorithms 41 Monte Carlo Integration By coincidence, the number of points within the given circle was exactly 800 (in this case) – By our algorithm, this would mean that an approximation of the area of that circle is 800/1000 = 0.8 – The correct area is 0.25p ≈ 0.7853981635 – 0.8 is not that poor an approximation... http://xkcd.com/10/Stochastic algorithms 42 Monte Carlo Integration Suppose we wanted to integrate the function 11 f v==   2 22 xy  v 2 over the circle given in the previous region – For each point (x , x ) that is in the circle, evaluate the function and 1 2 average these – In this example, we have: 800 1 1 2203.972409 f  = = 2.754965511 V  22 800 xy  800 k=1 kkStochastic algorithms 43 Monte Carlo Integration In this case, the approximation of the integral is f (vv )d = f V  0.82.754965511= 2.203972409 V  V To determine how well this algorithm works, use geometry, calculus and MapleStochastic algorithms 44 Monte Carlo Integration The region we are integrating over has: – For any value of x such that 0 ≤ x ≤ 1, 2 2 1 1 – We have that y runs from  xx to xx 2 2 y xStochastic algorithms 45 Monte Carlo Integration lcirc := 1/2 – sqrt(x x2); ucirc := 1/2 + sqrt(x – x2); int( int( 1/(x2 + y2), y = lcirc..ucirc ), x = 0..1 ); evalf( ); Maple got one integral, we need numerical algorithms for the second 2.1775860903 Our approximation: f (vv )d  2.203972409  VStochastic algorithms 46 Monte Carlo Integration Our approximation was, again, reasonably close considering we did no calculus: f (xx )d= 2.1775860903 2.203972409  V The relative error is 1.2 Stochastic algorithms 47 Monte Carlo Testing Suppose you have created a reasonably complex IC – Your design will assume ―ideal‖ characteristics – Reality is different: Consider the propagation delay of an inverter David Rennie, ECE 437Stochastic algorithms 48 Monte Carlo Testing Suppose you have created a reasonably complex IC – Accurate shapes will impossible for wiring and transistors David Rennie, ECE 437Stochastic algorithms 49 Monte Carlo Testing Suppose you have created a reasonably complex IC – There will be variability in doping Friedberg and Spanos, SPIE05Stochastic algorithms 50 Monte Carlo Testing Variability in digital gates will require designers to compensate for uncertainties Possible solutions: – Use and design for the worstcase delay for each gate • Consequences: unnecessary, pessimistic and results in overengineering – Use statistics • In ECE 316, you will learn about probability distributions • This, however, is computationally prohibitiveStochastic algorithms 51 Monte Carlo Testing A reasonable solution is to run a few hundred simulations where the parameter on the circuit elements, the delays, etc. are perturbed randomly within the expected distribution – Suppose one resistor is 120 W ± 5 W – For each simulation, create a standard normal variable x and then calculate 120 + 5x • This number will be a random sample with mean 120 and a standard deviation of 5 – Run the simulation with this resistor at that particular random valueStochastic algorithms 52 Randomized quicksort Recall from quicksort that we select a pivot and then separate the elements based on whether they are greater than or less than the pivot If we always choose the median of the first, middle, and last entries an array when performing quicksort, this gives us an algorithm for 2 devising a list which will run in Q(n ) time – The medianofthree algorithm, however, ensures that quicksort will not perform poorly if we pass it an already sorted listStochastic algorithms 53 Randomized quicksort Building this worstcase would be difficult: for each sublist, the median of the three must be either: – the second smallest number, or – the second largest number For example, this will give us one iteration of the worst case: 1, 3, 4, 5, 6, 7, 8, 2, 9, 10, 11, 12, 13, 14, 15Stochastic algorithms 54 Randomized quicksort However, further if we want three steps to fail, we must have the following form: 1, 3, 5, 7, 8, 9, 10, 2, 4, 6, 11, 12, 13, 14, 15 1, 2, 3, 5, 7, 8, 9, 10, 4, 6, 11, 12, 13, 14, 15 1, 2, 3, 4, 5, 7, 8, 9, 10, 6, 11, 12, 13, 14, 15 Watching this pattern, we may hypothesize that the worst case scenario is the input 1, 3, 5, 7, 9, 11, 13, 2, 4, 6, 8, 10, 12, 14, 15Stochastic algorithms 55 Randomized quicksort If instead of choose three random numbers and take the median of those three, then we would still partition the points into two sublists which are, on average, of a ratio of size 2:1 – This will still give us, on average, the expected run timeStochastic algorithms 56 Randomized quicksort However, because the choice of these three changes with each application of quicksort, we cannot create one sequence which must 2 result in Q(n ) behaviour Thus, for a given list, it may run more slow than with the deterministic algorithm, but it will also speed up other cases, and we cannot construct one example which must failStochastic algorithms 57 Skip lists A linked list has desirable insertion and deletion characteristics if we have a pointer to one of the nodes – Searching a linked list with n elements, however is O(n) – Can we achieve some of the characteristics of an array for searchingStochastic algorithms 58 Skip lists In an array, we have a means of finding the central element in O(1) time – To achieve this in a linked list, we require a pointer to the central elementStochastic algorithms 59 Skip lists Consider the following linked list... – searching is O(n):Stochastic algorithms 60 Skip lists Suppose, however, if we had an array of head pointersStochastic algorithms 61 Skip lists First, we could point to every second node – but this is still O(n/2) = O(n)Stochastic algorithms 62 Skip lists th We continue by pointing to every 4 entryStochastic algorithms 63 Skip lists th And then every 8 entry, and so on...Stochastic algorithms 64 Skip lists Suppose we search for 47 th – Following the 4 level pointer, first 32 47 but the next pointer is 0Stochastic algorithms 65 Skip lists Suppose we search for 47 rd – Continuing with the 3 level pointer from 32, 53 47Stochastic algorithms 66 Skip lists Suppose we search for 47 nd – Continuing with the 2 level pointer from 32, 47 42 but 53 47Stochastic algorithms 67 Skip lists Suppose we search for 47 st – Continuing with the 1 level pointer from 42, we find 47Stochastic algorithms 68 Skip lists Suppose we search for 24 th – With the 4 level pointer, 32 24Stochastic algorithms 69 Skip lists Suppose we search for 24 rd – Following the 3 level pointers, 4 24 but 32 24Stochastic algorithms 70 Skip lists Suppose we search for 24 nd – Following the 2 level pointer from 14, 23 24 but 32 24Stochastic algorithms 71 Skip lists Suppose we search for 24 st – Following the 1 level pointer from 23, 25 24 – Thus, 24 is not in the skip listStochastic algorithms 72 Skip lists Problem: in order to maintain this structure, we would be forced to perform a number of expensive operations – E.g., suppose we insert the number 8Stochastic algorithms 73 Skip lists Suppose that with each insertion, we don’t require that the exact shape is kept, but rather, only the distribution of how many next pointers a node may have Instead, note that: 1 node has 4 next pointers 2 nodes have 3 next pointers 4 nodes have 2 next pointers, and 8 nodes have only 1 next pointerStochastic algorithms 74 Skip lists h h + 1 If we are inserting into a skip list with 2 – 1 ≤ n 2 , we will choose a height between 1 and h such that the probability of having the k 1 2 height k is h 21  – We will choose a value between 1 and ⌊lg(n + 1)⌋ For example: Table 1. For 7 ≤ n 15 Table 2. For 15 ≤ n 31 Height Probability Height Probability 1 4/7 1 8/15 2 2/7 2 4/15 3 1/7 3 2/15 4 1/15Stochastic algorithms 75 Table 2. For 15 ≤ n 31 Skip lists Value k 0001 1 0010 2 To get such a distribution, we need a 0011 1 h number from 1 to 2 – 1 and determine 0100 3 the location of the leastsignificant bit 0101 1 that is equal to 1 0110 2 Table 1. For 7 ≤ n 15 0111 1 Value k 1000 4 001 1 1001 1 010 2 1010 2 011 1 1011 1 100 3 1100 3 101 1 1101 1 110 2 1110 2 111 1 1111 1Stochastic algorithms 76 Skip List Example This is relatively easy to program: int newheight() const int n = 1 height; // Generate a random number between 1 and 2height 1, inclusive int rand = ( lrand48() (n 1) ) + 1; // Starting with the leastsignificant bit, check if it equals 1 for ( int mask = 1, k = 1; true; mask = 1, ++k ) // If the bit is set to 1, return k if ( rand mask ) return k; // This should never be executed: requires include cassert assert( false ); Stochastic algorithms 77 Skip List Example For example, suppose we insert the following numbers into a skip list: 94, 28, 15, 45, 31, 52, 88, 76Stochastic algorithms 78 Skip List Example The first number, 94, is placed into the list as if it were a sorted linked listStochastic algorithms 79 Skip List Example The second number, 28, is also simply inserted into the correct locationStochastic algorithms 80 Skip List Example Now ⌊lg(3 + 1)⌋ = 2; hence for the third number, 15, we pick a random number between 1 and 3: 3 = 11 2Stochastic algorithms 81 Skip List Example Similarly with for the next number, 45, we pick a random number between 1 and 3: 1 = 01 2Stochastic algorithms 82 Skip List Example When we insert 31, the random number was 2 = 10 : 2Stochastic algorithms 83 Skip List Example For 52, the random number is 3 = 11 : 2Stochastic algorithms 84 Skip List Example Inserting 88, we find ⌊lg(7 + 1)⌋ = 3: 5 = 101 : 2Stochastic algorithms 85 Skip List Example Inserting 76, the random number is 4 = 100 : 2Stochastic algorithms 86 Skip List Example The problem with skip lists is that it doesn’t seem initially intuitive that it works – We need significantly larger examples – There is an implementation of skip lists at http://ece.uwaterloo.ca/ece250/Algorithms/Skiplists/Stochastic algorithms 87 Skip List Example The following are five examples of skip lists of randomly generated listsStochastic algorithms 88 Skip List Search Time This plot shows: – The average number of searches for a binary search (blue) – The average number of searches in a skip list with a best fitting line of the form a + b ln(n)Stochastic algorithms 89 References Wikipedia, http://en.wikipedia.org/wiki/Topologicalsorting 1 Cormen, Leiserson, and Rivest, Introduction to Algorithms, McGraw Hill, 1990, §11.1, p.200. rd 2 Weiss, Data Structures and Algorithm Analysis in C++, 3 Ed., Addison Wesley, §9.2, p.3425. These slides are provided for the ECE 250 Algorithms and Data Structures course. The material in it reflects Douglas W. Harder’s best judgment in light of the information available to him at the time of preparation. Any reliance on these course slides by any party for any other purpose are the responsibility of such parties. Douglas W. Harder accepts no responsibility for damages, if any, suffered by any party as a result of decisions made or actions based on these course slides for any other purpose than that for which it was intended.
Website URL
Comment