1 00:00:02,910 --> 00:00:08,620 Hello, bonjour, konnitchiwa, god dag, 2 00:00:08,620 --> 00:00:15,390 welcome to the first tutorial of Statistical Mechanics: Algorithms and Computations, 3 00:00:15,690 --> 00:00:18,660 from the Physics Department of Ecole Normale Superieure. 4 00:00:19,290 --> 00:00:26,780 Our subject today will be the convergence of Markov chain Monte Carlo algorithms 5 00:00:27,440 --> 00:00:33,650 and our main tool will be this stone, this pebble, 6 00:00:34,060 --> 00:00:39,180 moving around the 3x3 pebble game. 7 00:00:39,180 --> 00:00:46,820 It will allow us to derive and illustrate a profound mathematical result 8 00:00:46,820 --> 00:00:49,700 that is of great practical importance. 9 00:00:50,060 --> 00:00:58,230 We will show that the convergence of Markov chain Monte Carlo algorithms is exponential. 10 00:00:58,710 --> 00:01:03,760 This means that a Markov chain Monte Carlo algorithm 11 00:01:03,760 --> 00:01:06,950 is characterized by a time scale, 12 00:01:06,950 --> 00:01:09,500 the convergence time tau. 13 00:01:10,040 --> 00:01:14,470 It is on this scale that we approach equilibrium 14 00:01:15,030 --> 00:01:20,200 and that the pebble positions on the 3x3 pebble game 15 00:01:20,200 --> 00:01:22,550 become equally probable. 16 00:01:23,120 --> 00:01:25,880 After a few times tau, 17 00:01:26,360 --> 00:01:29,540 a Markov Chain Monte Carlo algorithm 18 00:01:29,540 --> 00:01:31,140 is in equilibrium 19 00:01:31,440 --> 00:01:34,060 for all intents and purposes. 20 00:01:34,920 --> 00:01:40,780 In the same way that after a few nuclear half-lives 21 00:01:41,130 --> 00:01:46,540 a radioactive compound has completely decayed. 22 00:01:47,170 --> 00:01:51,450 Let us now return to the 3x3 pebble game 23 00:01:51,450 --> 00:01:55,610 and understand its exponential convergence. 24 00:01:56,120 --> 00:02:00,130 On the 3x3 grid, we label the sites 25 00:02:00,550 --> 00:02:06,520 from 0 in the lower left corner to 8 in the upper right corner. 26 00:02:07,650 --> 00:02:14,040 The Markov chain Monte Carlo algorithm, as discussed in Lecture 1, 27 00:02:14,370 --> 00:02:17,610 consists in moving the pebble 28 00:02:17,610 --> 00:02:23,470 with probability 1/4 to the right and to the left, 29 00:02:23,470 --> 00:02:27,040 up, and down. 30 00:02:27,840 --> 00:02:34,470 When a move is not possible, it is rejected 31 00:02:34,470 --> 00:02:39,700 and the pebble remains in the same position for the next iteration. 32 00:02:40,680 --> 00:02:44,960 In Python, this gives the following program. 33 00:02:45,440 --> 00:02:52,100 The list 'neighbor' contains nine lists with four sites each, 34 00:02:52,840 --> 00:03:00,540 corresponding to the moves to the right, up, to the left and down. 35 00:03:01,550 --> 00:03:06,720 So here is the configuration at time 0 36 00:03:06,720 --> 00:03:09,640 with the pebble on site 8. 37 00:03:10,170 --> 00:03:14,660 The function randint(0,3) 38 00:03:15,520 --> 00:03:21,890 samples random numbers 0, 1, 2 and 3 with equal probability. 39 00:03:23,130 --> 00:03:29,820 If we choose 0, we'll try to move to the right, but this move is rejected. 40 00:03:30,750 --> 00:03:36,720 If we have randint equal to one, we'll move up, but again this move will be rejected. 41 00:03:37,260 --> 00:03:44,180 If randint equal to 2, we'll move to the left; and if it's equal to 3, we'll move down. 42 00:03:45,880 --> 00:03:50,570 Before running this algorithm, let us ask a few questions. 43 00:03:52,540 --> 00:03:53,810 First question: 44 00:03:55,300 --> 00:04:04,220 if at time t=0 the pebble is in the upper right corner, on site 8, 45 00:04:06,210 --> 00:04:11,740 on which sites can it be at time t=1? 46 00:04:15,840 --> 00:04:24,820 Of course, it can be on the sites 5, 7 and 8. 47 00:04:25,510 --> 00:04:31,690 In fact, the pebble will remain on site 8 with probability 1/2, 48 00:04:32,730 --> 00:04:38,970 and it will move to sites 5 and 7 with probability 1/4 each. 49 00:04:39,810 --> 00:04:40,970 Question 2: 50 00:04:41,620 --> 00:04:49,290 at what time does the probability to be at site 0 become finite? 51 00:04:52,800 --> 00:04:59,720 The answer is that it becomes finite at time t=4. 52 00:05:00,350 --> 00:05:10,400 at t=0, 1, 2 and 3 because it takes at least 4 steps to connect the site 8 53 00:05:10,400 --> 00:05:12,400 with the site 0. 54 00:05:12,400 --> 00:05:20,360 The probability to be at site 0 is non-zero for time t=4 and larger. 55 00:05:21,160 --> 00:05:25,230 Now, final question in this sequence. Question 3: 56 00:05:27,310 --> 00:05:35,280 what is the probability to be on site 0 at time t=4, 57 00:05:35,670 --> 00:05:40,840 given that the pebble started at site 8 at t=0? 58 00:05:42,830 --> 00:05:48,180 Attention, the probability is a fraction, 59 00:05:48,180 --> 00:05:50,920 and it is not easy to find. 60 00:05:54,750 --> 00:06:05,370 The answer is that this probability is equal to 6/256, or 3/128. 61 00:06:07,360 --> 00:06:13,510 The denominator (256) is equal to 4 to the power of 4 62 00:06:14,610 --> 00:06:19,610 because there are four moves and in each move we have four possibilities. 63 00:06:21,600 --> 00:06:32,870 The numerator (6) corresponds to the 6 ways to connect site 8 to site 0 in four steps. 64 00:06:33,940 --> 00:06:42,560 From our stone age 3x3 pebble game we now move on to a program capable of doing graphics. 65 00:06:43,240 --> 00:06:48,030 The red pebble that you see here starts on site 8, 66 00:06:48,030 --> 00:06:54,950 then moves on to neighboring sites, hops back, goes on, touches site 0 67 00:06:54,950 --> 00:06:56,740 and so on and so on, 68 00:06:56,740 --> 00:07:02,540 until - after long time - all the positions become equally likely. 69 00:07:03,040 --> 00:07:08,120 Next in this tutorial, Michael will discuss the simulation 70 00:07:08,120 --> 00:07:13,270 without watching pebbles hopping around the 3x3 pebble game. 71 00:07:14,490 --> 00:07:19,480 But before moving on, please take a moment 72 00:07:19,480 --> 00:07:24,560 to download, run and modify the programs we just discussed. 73 00:07:25,720 --> 00:07:32,590 On the Coursera website, you will find the program pebble_basic.py 74 00:07:32,650 --> 00:07:34,350 that we just discussed. 75 00:07:34,350 --> 00:07:40,170 You will also see pebble_basic_movie.py 76 00:07:40,170 --> 00:07:45,050 that produces the nice graphics output. 77 00:07:50,490 --> 00:07:55,750 Now let us modify a bit the simulation program for the 3x3 pebble game. 78 00:07:55,750 --> 00:08:00,720 Rather than having one long run, we will now consider many runs 79 00:08:00,720 --> 00:08:05,320 and each of the run will start at site 8 in the upper right corner of the grid. 80 00:08:05,830 --> 00:08:09,460 Now this program contains an additional loop 81 00:08:09,460 --> 00:08:14,270 over the runs, and we have a total of n_runs runs in the program. 82 00:08:14,540 --> 00:08:19,560 The program finally outputs the final position of the pebble 83 00:08:19,560 --> 00:08:22,240 at time t_max. 84 00:08:22,240 --> 00:08:26,700 Output for t_max=4 is shown here. 85 00:08:27,740 --> 00:08:33,090 Let us now represent the output using two-dimensional histograms, 86 00:08:33,090 --> 00:08:37,070 first for the case t_max=0. 87 00:08:37,700 --> 00:08:43,790 For this time, all the runs are still on site 8, 88 00:08:43,790 --> 00:08:47,270 and there are no pebbles on other sites. 89 00:08:47,570 --> 00:08:51,020 As all the pebbles are on site 8, 90 00:08:51,020 --> 00:08:54,940 their proportion on site 8 is equal to one, 91 00:08:55,120 --> 00:08:56,940 here represented in white. 92 00:08:57,530 --> 00:09:04,630 All the other sites are shown in black, as they do not appear in the output of the program. 93 00:09:05,550 --> 00:09:13,200 At t_max=1, that is after one iteration, the output of our program looks as follows, 94 00:09:14,240 --> 00:09:17,650 and the corresponding histogram is shown here. 95 00:09:18,100 --> 00:09:24,460 At this time, about half of the pebbles is still on site 8, 96 00:09:25,210 --> 00:09:32,580 while there is about a quarter of the pebbles on site 5 and site 7 each. 97 00:09:33,050 --> 00:09:38,940 The other sites are still black, because they cannot be reached within one move. 98 00:09:39,300 --> 00:09:44,290 At t_max=2, we have the following histogram. 99 00:09:44,740 --> 00:09:51,640 we see that now some pebbles also reach the sites 2, 4 and 6. 100 00:09:51,640 --> 00:09:58,770 For larger times we see how the pebble diffuses over the entire 3x3 pebble game 101 00:09:59,160 --> 00:10:02,810 till the color becomes more and more uniform. 102 00:10:03,200 --> 00:10:08,460 Please take a moment to make sure that you understand the color code in this figure, 103 00:10:08,460 --> 00:10:12,330 and test your understanding in the following little quiz. 104 00:10:12,800 --> 00:10:17,350 We consider the program with n_runs=100. 105 00:10:18,240 --> 00:10:22,260 As is indicated by the color bar on the right of the diagram 106 00:10:22,460 --> 00:10:28,620 the color of site 0 represents the number 0.09, 107 00:10:28,620 --> 00:10:30,620 which is output of the program. 108 00:10:31,000 --> 00:10:37,390 Question 1: what is the meaning of this number 0.09? 109 00:10:41,640 --> 00:10:44,730 The correct answer is answer 2. 110 00:10:44,850 --> 00:10:48,060 The histogram simply analyzes the output 111 00:10:48,060 --> 00:10:52,070 of a simulation with n_runs=100 112 00:10:53,090 --> 00:10:54,630 Question 2: 113 00:10:54,990 --> 00:11:01,800 on the starting site 8, the probability falls from 1 at t= 0 114 00:11:02,120 --> 00:11:04,590 to one half at t=1 115 00:11:04,860 --> 00:11:08,960 and is always larger than the probability at any other site. 116 00:11:09,290 --> 00:11:10,960 Is this always true? 117 00:11:11,190 --> 00:11:15,470 Even if we choose another starting site at t=0? 118 00:11:17,460 --> 00:11:24,900 It is false, as can be seen from the case of the pebble game starting at site 4, the center. 119 00:11:25,280 --> 00:11:29,680 The probability to be found at the starting point at t=1 120 00:11:30,040 --> 00:11:32,600 is zero, in this case. 121 00:11:33,520 --> 00:11:35,150 Question 3: 122 00:11:35,150 --> 00:11:42,170 would our results change, if instead of playing 100000 consecutive pebble games 123 00:11:42,170 --> 00:11:47,490 we would let 100000 pebbles move across the grid at the same time? 124 00:11:49,120 --> 00:11:51,590 The result would not change at all. 125 00:11:52,130 --> 00:11:55,610 The probabilities represented by the color code in the figure 126 00:11:55,610 --> 00:12:02,030 can be correctly interpreted as the normed density of an ensemble of a large number of pebbles 127 00:12:02,270 --> 00:12:06,040 moving across the sites, without interactions between them. 128 00:12:06,840 --> 00:12:14,160 Next in this Tutorial 1, Alberto will discuss the simulation using a compact analytic method, 129 00:12:14,160 --> 00:12:16,160 the transfer matrix. 130 00:12:16,160 --> 00:12:22,510 It will give us the exact value of the probabilities that we just discussed. 131 00:12:22,510 --> 00:12:24,510 Before moving on, 132 00:12:24,770 --> 00:12:30,510 take your time to download, run and modify the programs. 133 00:12:30,510 --> 00:12:36,310 On the Coursera website, you will find the program pebble_multirun.py 134 00:12:36,310 --> 00:12:43,680 and also the program that allows you to generate the histogram at time t_max. 135 00:12:44,210 --> 00:12:48,820 There's also a program which generates all the histograms 136 00:12:48,820 --> 00:12:54,440 at times t=0, t=1 and so on, in one run. 137 00:12:59,490 --> 00:13:04,610 As indicated by Michael, we now discuss the transfer matrix method, 138 00:13:04,610 --> 00:13:10,320 which allows to solve exactly the Monte Carlo dynamics of the 3x3 pebble game. 139 00:13:10,580 --> 00:13:15,840 For larger system, this method does not apply, but its lesson remains valid. 140 00:13:16,260 --> 00:13:19,500 If we start at the right upper corner at t=0, 141 00:13:19,500 --> 00:13:23,540 we know that the probability to find the pebble at site 8 is one. 142 00:13:24,530 --> 00:13:31,630 At time t=1, we know that with probability 1/4 we will be at the site 5 or 7 143 00:13:31,630 --> 00:13:34,430 and with probability 1/2 at the site 8. 144 00:13:35,140 --> 00:13:38,830 To compute the exact probabilities at all time steps, 145 00:13:38,830 --> 00:13:42,100 we need to introduce a vector: pi_t, 146 00:13:42,100 --> 00:13:52,410 which contains the probabilities to find the pebble at time t at the site 0, 1, 2 and so on. 147 00:13:52,920 --> 00:13:56,250 As explained by Werner in Lecture 1, 148 00:13:56,250 --> 00:14:01,120 the Monte Carlo algorithm is nothing but a matrix, the transfer matrix 149 00:14:01,210 --> 00:14:05,670 of transition probabilities p(a->b) 150 00:14:05,670 --> 00:14:09,540 for moving from the site a to the site b. 151 00:14:10,370 --> 00:14:16,940 In the 3x3 pebble game, this matrix is 9x9, because there are 9 configurations. 152 00:14:18,390 --> 00:14:20,650 Let's fill in this matrix. 153 00:14:21,310 --> 00:14:25,530 The probability to go from the site 0 to 0 is 1/2, 154 00:14:25,530 --> 00:14:29,810 the probability to go from 0 to 1 is 1/4, 155 00:14:29,810 --> 00:14:34,840 the probability to go from 1 to 2 is 1/4, 156 00:14:34,840 --> 00:14:38,080 so let's fill in the entire matrix. 157 00:14:38,940 --> 00:14:43,610 Note the factor 1/4 in front of the matrix. 158 00:14:44,650 --> 00:14:47,170 It is now crucial to realize 159 00:14:47,170 --> 00:14:50,560 that the probability vector at time t+1 160 00:14:50,890 --> 00:14:53,770 is the product between the transfer matrix 161 00:14:53,770 --> 00:14:56,420 and the probability vector at time t. 162 00:14:56,420 --> 00:14:59,270 In Python, this gives the following program, 163 00:14:59,750 --> 00:15:04,450 where thanks to numpy we can write in one line the matrix vector product. 164 00:15:04,830 --> 00:15:12,260 The output is the probability vector at time 0, 1, 2 and so on, as shown here. 165 00:15:13,450 --> 00:15:16,630 This is the entire Monte Carlo dynamics, 166 00:15:16,630 --> 00:15:20,410 for an infinite number of pebble games, and it is exact. 167 00:15:20,740 --> 00:15:24,070 For example the probability on the site 0 168 00:15:24,070 --> 00:15:28,620 is zero at time 0, 1, 2 and 3 169 00:15:28,820 --> 00:15:39,140 and is 6/256 - which means 0.023 - at time 4, as you proved earlier. 170 00:15:39,820 --> 00:15:41,140 We also see 171 00:15:41,580 --> 00:15:46,570 that the equilibrium value of 1/9, which is 0.111, 172 00:15:46,570 --> 00:15:50,260 is reached very fast, for all sites in the system. 173 00:15:51,300 --> 00:15:57,720 At small time, the matrix-vector multiplication modifies the probability vector, 174 00:15:57,720 --> 00:16:04,910 while at longer times this multiplication does not change the vector. 175 00:16:04,910 --> 00:16:08,100 Does this ring a bell with you? 176 00:16:08,100 --> 00:16:15,320 It means that the equilibrium probability vector is an eigenvector of the transfer matrix, 177 00:16:15,320 --> 00:16:17,320 with eigenvalue 1. 178 00:16:18,830 --> 00:16:23,640 Eigenvalues and eigenvectors can be computed in one line, in Python, 179 00:16:24,000 --> 00:16:28,850 and this is the modified program pebble_transfer.py 180 00:16:29,170 --> 00:16:33,450 and the nine eigenvalues of the transfer matrix. 181 00:16:33,840 --> 00:16:37,560 We see that 1 is the largest eigenvalue, 182 00:16:37,770 --> 00:16:42,080 the second largest eigenvalue is 0.75, 183 00:16:42,080 --> 00:16:45,290 then we have 0.5 and -0.5. 184 00:16:46,120 --> 00:16:49,890 We have already understood the role of the eigenvalue 1, 185 00:16:49,890 --> 00:16:53,310 it is associated to the equilibrium eigenvector. 186 00:16:53,880 --> 00:17:01,490 Now we will study what the other eigenvalues, and in particular 0.75, are good for. 187 00:17:01,960 --> 00:17:07,020 To do so, we consider again the output of pebble_transfer.py. 188 00:17:07,380 --> 00:17:11,510 Here we find the probability to reach any site 189 00:17:11,510 --> 00:17:14,600 starting from the upper right corner. 190 00:17:15,080 --> 00:17:21,080 Let us now subtract the equilibrium value, 0.1111, 191 00:17:21,080 --> 00:17:23,100 and take the absolute value. 192 00:17:23,370 --> 00:17:28,390 This is done in pebble_transfer_sub.py, 193 00:17:28,390 --> 00:17:30,390 and this is the output 194 00:17:32,290 --> 00:17:37,640 Let us now plot the output for the site 0, in the semilog scale. 195 00:17:38,200 --> 00:17:44,330 The straight line indicates an exponential convergence to equilibrium 196 00:17:44,330 --> 00:17:47,390 as announced at the beginning of the tutorial. 197 00:17:47,390 --> 00:17:55,510 The slope of this line is the same as the function (0.75)^t 198 00:17:55,830 --> 00:17:58,960 the second eigenvalue power t. 199 00:17:59,190 --> 00:18:04,490 Let us now put 0.75 into an exponential 200 00:18:04,660 --> 00:18:11,470 this gives exp(-t ln(0.75)), 201 00:18:11,470 --> 00:18:16,470 which can also be written as exp(-t/tau). 202 00:18:17,180 --> 00:18:23,210 tau, which is more or less 3.4, is the correlation time 203 00:18:23,300 --> 00:18:26,570 introduced by Werner at the beginning of the tutorial. 204 00:18:27,020 --> 00:18:32,910 Let's go back to the original evolution of the probability vector pi_t. 205 00:18:34,100 --> 00:18:41,880 It is after few correlation times tau - let's say 4 or 5 - that the equilibrium is reached 206 00:18:41,880 --> 00:18:48,510 and the corrections to the equilibrium disappear for all intents and purposes 207 00:18:49,080 --> 00:18:55,260 To understand our observations, we remark that the initial probability vector 208 00:18:55,260 --> 00:18:59,660 can be decomposed on the eigenvectors of the transfer matrix. 209 00:18:59,780 --> 00:19:04,810 The component associated with the eigenvalue 1 is conserved in time, 210 00:19:04,840 --> 00:19:08,460 because 1^t is equal to 1. 211 00:19:08,940 --> 00:19:14,940 The other components associated with the other eigenvalues lambda decay exponentially 212 00:19:15,270 --> 00:19:20,210 and the slowest decay is given by the second largest eigenvalue. 213 00:19:20,210 --> 00:19:27,550 Before moving on to Vivien's part about rigorous solutions, 214 00:19:27,940 --> 00:19:35,190 let's take a moment to download, run and modify the programs that we have discussed in this section. 215 00:19:35,930 --> 00:19:41,190 On the Coursera website, you will find pebble_transfer.py 216 00:19:41,400 --> 00:19:45,210 which contains the construction of the transfer matrix 217 00:19:45,210 --> 00:19:50,140 and the matrix-vector multiplication using numpy. 218 00:19:50,860 --> 00:19:58,500 The eigenvalues of the transfer matrix are computed in transfer_pebble_eigen.py. 219 00:19:59,450 --> 00:20:05,250 The graphics and the subtraction of the equilibrium value 1/9 220 00:20:05,250 --> 00:20:11,610 is contained in the program transfer_pebble_sub.py. 221 00:20:17,550 --> 00:20:23,050 Alberto has indicated that the evolution of the Monte Carlo dynamics 222 00:20:23,410 --> 00:20:29,950 is determined by the decomposition of the probability vector at initial time 223 00:20:30,160 --> 00:20:34,650 into the eigenvectors of the transfer matrix. 224 00:20:35,240 --> 00:20:41,600 All eigenvalues have modulus less than or equal to one, 225 00:20:41,990 --> 00:20:49,010 an eigenvalue larger than one would lead to an explosion of the probability 226 00:20:49,010 --> 00:20:51,010 at large time. 227 00:20:51,410 --> 00:20:54,920 At large times, the dynamics converges 228 00:20:54,920 --> 00:21:00,240 to the eigenvector with the largest eigenvalue equal to one. 229 00:21:00,990 --> 00:21:08,330 this is what we will call the probability vector of the equilibrium state, 230 00:21:08,330 --> 00:21:10,330 or the steady state. 231 00:21:11,100 --> 00:21:14,750 When does this picture really apply? 232 00:21:15,200 --> 00:21:19,090 What are the rigorous mathematical conditions 233 00:21:19,090 --> 00:21:26,230 ensuring that one converges at large times to a unique steady state? 234 00:21:27,000 --> 00:21:33,480 This is what we will find out in the final part of this tutorial. 235 00:21:33,480 --> 00:21:39,750 There is a simple physical example in which a strange phenomenon occurs. 236 00:21:40,530 --> 00:21:47,840 Consider a system in which you have two copies of the 3x3 pebble game, 237 00:21:47,840 --> 00:21:52,030 one in red and one in blue. 238 00:21:52,510 --> 00:21:58,600 In each of these games, the pebble jumps from site to site 239 00:21:58,600 --> 00:22:02,410 with the same rules as we have seen before. 240 00:22:03,000 --> 00:22:09,280 If it is started in the red one, it will remain there forever 241 00:22:09,390 --> 00:22:16,650 and if it is started in the blue one the same thing will happen. 242 00:22:16,650 --> 00:22:24,410 The combined system is described by a 18x18 transfer matrix 243 00:22:24,410 --> 00:22:28,990 because we now have 18 sites. 244 00:22:29,760 --> 00:22:33,120 There are lot of zeros in this matrix 245 00:22:33,120 --> 00:22:41,320 indicating that the pebble cannot jump from the red game to the blue game and vice versa. 246 00:22:41,320 --> 00:22:48,040 In fact, we are in an annoying and physically unacceptable situation, 247 00:22:48,040 --> 00:22:54,880 where the outcome of the simulation depends on the starting configuration 248 00:22:54,880 --> 00:22:58,540 even at infinite times. 249 00:22:58,950 --> 00:23:02,970 We are also in a mathematically annoying situation, 250 00:23:02,970 --> 00:23:10,760 because the transfer matrix of our dual pebble game has two eigenvalues one. 251 00:23:11,230 --> 00:23:18,900 This is described in the program pebble_dual_eigen.py. 252 00:23:19,880 --> 00:23:30,610 Mathematicians describe this situation that we must avoid as stemming from the reducibility of the transfer matrix 253 00:23:30,610 --> 00:23:33,350 of our dual game. 254 00:23:33,770 --> 00:23:43,130 Reducible means that we can split apart the transfer matrix of our system 255 00:23:43,130 --> 00:23:49,190 into two subsystems without affecting the rest. 256 00:23:49,190 --> 00:23:58,140 One of the two mathematically rigorous conditions, for a Monte Carlo Markov chain algorithm 257 00:23:58,140 --> 00:24:03,460 is that this pulling apart is not possible, 258 00:24:03,460 --> 00:24:06,110 that would be irreducible. 259 00:24:06,350 --> 00:24:11,400 We can in fact render our dual pebble game irreducible, 260 00:24:11,400 --> 00:24:15,680 by adding a small probability of transition 261 00:24:15,680 --> 00:24:17,760 between the two games 262 00:24:17,760 --> 00:24:20,320 as depicted here. 263 00:24:20,650 --> 00:24:27,130 We now observe that the pebble can change its game. 264 00:24:27,130 --> 00:24:34,710 This small probability can also be observed in the transfer matrix 265 00:24:34,710 --> 00:24:40,210 in the program pebble_dual_eigen.py, 266 00:24:40,210 --> 00:24:48,800 where it allows to recover a unique eigenvector of eigenvalue one 267 00:24:48,800 --> 00:24:53,050 that is to say, a unique steady state. 268 00:24:53,410 --> 00:24:59,890 We realize we must have a unique eigenvector of eigenvalue one, 269 00:24:59,890 --> 00:25:05,030 ensuring that the matrix is irreducible. 270 00:25:05,600 --> 00:25:08,300 In fact, this is not enough. 271 00:25:09,020 --> 00:25:11,460 We need a second condition, 272 00:25:11,460 --> 00:25:16,180 ensuring that the probability converges 273 00:25:16,180 --> 00:25:20,910 to the steady state, in the large time limit. 274 00:25:22,280 --> 00:25:27,840 We could face the situation where the second eigenvalue lambda_2 275 00:25:27,840 --> 00:25:31,840 (different from one and possibly complex) 276 00:25:31,840 --> 00:25:35,150 is of unit modulus. 277 00:25:35,150 --> 00:25:39,700 In this situation, lambda_2 to the power of t 278 00:25:39,700 --> 00:25:44,630 does not converge to zero as t increases. 279 00:25:45,110 --> 00:25:49,780 This means that the probability vector oscillates 280 00:25:49,870 --> 00:25:54,770 and never converges at infinite times. 281 00:25:55,370 --> 00:25:59,350 To illustrate this behavior in a simple example, 282 00:25:59,350 --> 00:26:03,750 consider the following 2x1 pebble game 283 00:26:03,750 --> 00:26:06,990 and adopt the following jump rule. 284 00:26:06,990 --> 00:26:09,800 The poor pebble is unlucky. 285 00:26:09,800 --> 00:26:15,000 Not only it has only two sites where to be, 286 00:26:15,540 --> 00:26:23,820 but also at each time step it must go to the other accessible position. 287 00:26:25,630 --> 00:26:28,610 Starting from the left at time 0, 288 00:26:28,610 --> 00:26:32,110 it has to go to the right at time 1, 289 00:26:32,110 --> 00:26:38,130 to the left at time 2, and so on, in a deterministic way. 290 00:26:40,930 --> 00:26:45,500 The transfer matrix of this dynamics is simple to write. 291 00:26:46,740 --> 00:26:51,610 The ones in the matrix indicate that the particles jumps 292 00:26:51,610 --> 00:26:55,210 to its neighboring site at each time step, 293 00:26:55,210 --> 00:27:00,550 and the zeros indicate that it cannot stay in the same position. 294 00:27:01,830 --> 00:27:05,950 The eigenvectors and eigenvalues take a simple form. 295 00:27:07,060 --> 00:27:09,320 The steady state is uniform 296 00:27:09,320 --> 00:27:13,960 and the other eigenvector has eigenvalue -1. 297 00:27:14,540 --> 00:27:16,160 This explains why 298 00:27:16,160 --> 00:27:21,370 the probability vector does not converge in the infinite time limit. 299 00:27:22,590 --> 00:27:25,370 How is this phenomenon general? 300 00:27:26,090 --> 00:27:29,370 Let us remark that every two time steps 301 00:27:29,810 --> 00:27:34,250 the pebble comes back to its original state. 302 00:27:35,280 --> 00:27:37,380 From a general point of view 303 00:27:37,660 --> 00:27:42,740 we now arrived at the second of the two mathematical conditions 304 00:27:42,950 --> 00:27:48,890 that we require for our Monte Carlo Markov chain algorithm. 305 00:27:50,010 --> 00:27:53,390 We must avoid recurring states: 306 00:27:53,590 --> 00:27:58,340 situations where when starting from one state 307 00:27:58,340 --> 00:28:04,500 we are sure to find it back after a fixed number of time steps. 308 00:28:05,450 --> 00:28:10,770 In our 2x1 pebble game, we have two recurring states 309 00:28:10,990 --> 00:28:15,280 which we are sure to find back after two time steps. 310 00:28:15,280 --> 00:28:19,920 They correspond to starting either from the right 311 00:28:20,560 --> 00:28:23,280 or from the left states. 312 00:28:24,520 --> 00:28:28,380 Fortunately recurring states can be easily avoided 313 00:28:28,680 --> 00:28:35,260 In our example, imagine that we allow the pebble to stay in the same configuration 314 00:28:35,260 --> 00:28:38,460 with a small probability epsilon. 315 00:28:38,620 --> 00:28:44,100 Every now and then, the pebble will stop and stay in its position 316 00:28:45,020 --> 00:28:50,140 In the program pebble_recurrent_eigen.py, 317 00:28:50,660 --> 00:28:56,030 you will see how with a small finite epsilon 318 00:28:56,320 --> 00:29:03,530 the steady state becomes the unique eigenvector of unit modulus. 319 00:29:03,760 --> 00:29:08,090 What we have seen here has been cast by mathematicians 320 00:29:08,090 --> 00:29:12,090 into two conditions that we must supplement 321 00:29:12,090 --> 00:29:15,670 to the global or detailed balance condition 322 00:29:15,670 --> 00:29:18,610 that Werner discussed during the lecture. 323 00:29:19,710 --> 00:29:22,850 The Markov chain must satisfy 324 00:29:22,850 --> 00:29:25,160 global balance and 325 00:29:25,160 --> 00:29:29,160 condition 1: it is irreducible 326 00:29:29,550 --> 00:29:33,600 and condition 2: it is aperiodic, 327 00:29:33,710 --> 00:29:37,390 that is: there are no recurrent states. 328 00:29:38,650 --> 00:29:43,540 Those are in total the only three conditions 329 00:29:43,540 --> 00:29:48,770 that we require the Markov chain algorithm to verify 330 00:29:48,770 --> 00:29:52,770 for it to behave in a proper way. 331 00:29:53,590 --> 00:29:58,450 Please take a moment to download, modify and run 332 00:29:58,450 --> 00:30:01,640 the programs that we have discussed in this session. 333 00:30:02,400 --> 00:30:07,000 On the Coursera website, you will find the two following programs 334 00:30:07,600 --> 00:30:15,650 pebble_dual_eigen.py and pebble_dual_movie.py 335 00:30:15,650 --> 00:30:21,430 In these two programs, you can add a small probability epsilon 336 00:30:21,760 --> 00:30:26,760 for the pebble to jump from the red game to the blue game. 337 00:30:27,170 --> 00:30:34,550 You will see how this allows the eigenstate of eigenvalue one 338 00:30:34,550 --> 00:30:39,320 that is to say the steady state, to be the unique one. 339 00:30:40,080 --> 00:30:47,320 You will also find the associate programs pebble_recurrent_eigen.py 340 00:30:47,320 --> 00:30:51,320 and pebble_recurrent_movie.py 341 00:30:51,320 --> 00:30:55,670 In these two programs, you can add a small probability epsilon 342 00:30:55,890 --> 00:30:59,670 allowing the pebble to stay in the same position 343 00:31:00,200 --> 00:31:02,560 an you will discover by yourself 344 00:31:02,560 --> 00:31:05,800 how the steady state becomes unique. 345 00:31:05,800 --> 00:31:13,470 In conclusion, we have seen that Markov chain Monte Carlo calculations 346 00:31:13,650 --> 00:31:17,780 are equivalent to transfer matrix multiplications. 347 00:31:18,870 --> 00:31:24,120 They converge to equilibrium after an infinite number of iterations 348 00:31:25,250 --> 00:31:29,020 But this infinite number is not a serious restriction 349 00:31:29,280 --> 00:31:36,790 because of the existence of a time scale, set by the second largest eigenvalue of the transfer matrix. 350 00:31:37,320 --> 00:31:45,630 In the 3x3 pebble game, the second eigenvalue gave a time scale tau=3.4 351 00:31:46,600 --> 00:31:54,200 and convergence towards equilibrium was reached after few times tau, like 15 or 20 iterations. 352 00:31:55,240 --> 00:31:59,170 The concept of equilibrium is far reaching 353 00:31:59,860 --> 00:32:06,550 and the interest in Monte Carlo calculations is very strong because of the existence of this time scale 354 00:32:07,830 --> 00:32:10,770 that separates fast and slow things. 355 00:32:11,720 --> 00:32:17,730 There is a considerable tragedy attached to this question of the timescale tau, 356 00:32:18,330 --> 00:32:25,990 because in usual Monte Carlo calculations we have a hard time estimating tau. 357 00:32:27,540 --> 00:32:30,730 But this will be food for thought, 358 00:32:30,730 --> 00:32:36,340 for another session of Statistical Mechanics: Algorithms and Computations.