+ Advanced Sampling Algorithms + Mobashir Mohammad Hirak Sarkar - PowerPoint PPT Presentation

+ 2D Ising Model (3) 35 −𝛾 𝑡 𝑦𝑗 𝑡(𝑦𝑘) 𝑦𝑗𝑦𝑘 ∈𝐹 𝑓 −𝛾𝐼(𝑡) 𝑓  P( 𝑡 ) = = 𝑎 𝑎  If 𝛾 = 0 all spin configurations have same probability  If 𝛾 > 0 lower energy preferred

+ Exact Inference is Hard 36  Posterior distribution over 𝑡  If 𝑍 is any observed variable 𝑞 𝑡,𝑧  𝑄 𝑡 𝑍) = 𝑞(𝑡,𝑧)  Intractable computation  Joint probability distribution - 2 𝑂 possible combinations of 𝑡  Marginal probability distribution at a site  MAP estimation

+ The Big Question 37  Given 𝜌 on 𝑇 = 𝑡 1 , 𝑡 2 , … , 𝑡 𝑜 simulate a random object with distribution 𝜌

+ Methods 38  Generate random samples to estimate a quantity  Samples are generated “Markov - chain style”  Markov Chain Monte Carlo (MCMC)  Propp Wilson Simulation  Las Vegas variant  Sandwiching  Improvement on Propp Wilson

+ MARKOV CHAIN MONTE CARLO METHOD (MCMC) YAMILET R. SERRANO LL. 39

+ Recall 40 Given a probability distribution π on S = {s 1 , …, s k }, how do we simulate a random object with distribution π ?

+ Intuition 41 Given a probability distribution π on S = {s 1 , …, s k }, how do we simulate a random object with distribution π ? High Dimension Space

+ MCMC 42 1. Construct an irreducible and aperiodic Markov Chain [X 0 ,X 1 , … ], whose stationary distribution π . 2. If we run the chain with arbitrary initial distribution then 1. The Markov Chain Convergence Theorem guarantees that the distribution of the chain at time n converges to π . 3. Hence, if we run the chain for a sufficiently long time n, then the distribution of X n , will be very close to π . So it can be used as a sample.

+ MCMC(2) 43  Generally, two types of MCMC algorithm  Metropolis Hasting  Gibbs Sampling

+ Metropolis Hastings Algorithm 44  Original method:  Metropolis, Rosenbluth, Rosenbluth, Teller and Teller (1953).  Generalized by Hasting in 1970.  Rediscovered by Tanner and Wong (1987) and Gelfang and Smith (1990)  Is one way to implement MCMC. Nicholas Metropolis

+ Metropolis Hastings Algorithm 45 Basic Idea GIVEN: A probability distribution π on S = {s 1 , …, s k } GOAL : Approx. sample from π Start with a proposal distribution Q(x,y) x: current state • Q(x,y) specifies transition of Markov Chain y: new proposal • Q(x,y) plays the role of the transition matrix By accepting/rejecting the proposal, MH simulates a Markov Chain, whose stationary distribution is π

+ Metropolis Hastings Algorithm 46 Algorithm A probability distribution π on S = {s 1 , …, s k } GIVEN: Approx. sample from π GOAL : Given current sample,  Draw y from the proposal distribution,  Draw U ~ (0,1) and update where the acceptance probability is

+ Metropolis Hastings Algorithm 47 The Ising Model Consider m sites around a circle. Each site i can have one of two spins x i  {-1,1} The target distribution :

+ Metropolis Hastings Algorithm The Ising Model Target Distribution: Proposal Distribution: 1. Randomly pick one out of the m spins 2. Flip its sign Acceptance probability (say i-th spin flipped):

+ Metropolis Hastings Algorithm 49 The Ising Model Image has been taken from Montecarlo investigation of the Ising Model(2006) – Tobin Fricke

+ Disadvantages of MCMC 50  No matter how large n is taken to be in the MCMC algorithm, there will still be some discrepancy between the distribution of the output and the target distribution π .  In order to make the previous error small, we need to figure out how large n needs to be.

+ Bounds on Convergence of MCMC Aditya Kulkarni 51

+ Seminar on probabilities 52 (Strasunbourge - 1983)  If it is difficult to obtain asymptotic bounds on the convergence time of MCMC algorithms for Ising models, use quantitative bounds  Use a characteristics of Ising model MCMC algorithm to obtain quantitative bounds David Aldous

+ What we need to ensure 53  As we go along the time We approach to the stationary distribution in a monotonically decreasing fashion. 𝑒 𝑢 → 0 as 𝑢 → ∞  When we stop The sample should follow a distribution which is not further apart from a stationary distribution by some factor 𝜻 . 𝑒 𝑢 ≤ 𝜁

+ Total Variation Distance || ⋅ || 𝑈𝑊 54  Given two distributions 𝑞 and 𝑟 over a finite set of states 𝑇 ,  The total variation distance ‖ ⋅ ‖ 𝑈𝑊 is 1 1 | 𝑞 − 𝑟 | 𝑈𝑊 = 2 ||𝑞 − 𝑟|| 1 = 2 𝑦∈𝑇 |𝑞 𝑦 − 𝑟 𝑦 |

+ Convergence Time 55  Let 𝑌 = 𝑌 0 , 𝑌 1 , … be a Markov Chain of a state space with stationary distribution 𝜌  We define 𝑒 𝑢 as the worst total variation distance at time 𝑢 . 𝑒 𝑢 = 𝑛𝑏𝑦 𝑦∈𝑊 ||𝑄 𝑌 𝑢 𝑌 0 = 𝑦 − 𝜌|| 𝑈𝑊  The mixing time 𝑢 is the minimum time 𝑢 such that 𝑒 𝑢 is at most 𝜗 . τ (𝜁) = min {𝑢 ∶ 𝑒 𝑢 ≤ 𝜁} 1 2𝑓 = min {𝑢 ∶ 𝑒 𝑢 ≤ 1 2𝑓 } τ = τ  Define

+ Quantitative results 56 1 𝒆 𝒖 𝟐 𝟑𝒇 𝜻 τ τ (𝜻) 0 𝒖

+ Lemma 57  Consider two random walks started from state i and state j Define 𝜍 𝑢 as the worst total variation distance between their respective probability distributions at time 𝑢 a. 𝜍 𝑢 ≤ 2𝑒 𝑢 b. 𝑒(𝑢) is decreasing

+ Upper bound on d ( 𝑢 ) proof 58 From part a and the definition of τ = τ 1 2𝑓 = min {𝑢 ∶ 𝑒 𝑢 ≤ 1 2𝑓 } we get 𝜍 𝜐 ≤ 𝑓 −1 Also from part b, we deduce that upper bound of 𝜍 at a particular time 𝑢 0 gives upper bound for later times 𝜍 𝑢 ≤ (𝜍(𝑢 0 )) 𝑜 ; 𝑜𝑢 0 ≤ 𝑢 ≤ (𝑜 + 1)𝑢 0 𝑢 𝑢 0 −1) 𝜍 𝑢 ≤ (𝜍(𝑢 0 )) ( Substitute 𝜐 instead of 𝑢 0 𝑢 𝜐 −1) 𝜍 𝑢 ≤ (𝜍(𝜐)) ( 𝑢 τ , 𝑢 ≥ 0 𝑒 𝑢 ≤ exp 1 −

+ Upper bound on τ ( 𝜁 ) proof 59  Algebraic calculations: 𝑢 τ 𝑒 𝑢 ≤ exp 1 − log 𝑒(𝑢) ≤ 1 − 𝑢 𝜐 𝑢 1 𝑒(𝑢) ) 𝜐 ≤ 1 − log 𝑒 𝑢 ≤ 1 + log( 𝑢 ≤ 𝜐 ∗ (1 + log( 1 𝜁 )) 𝜐 𝜁 ≤ 2𝑓𝑢 ∗ (1 + log( 1 𝜁 ))

+ Upper bounds 60 𝑢 τ ), 𝑢 ≥ 0 d 𝑢 ≤ min(1, exp 1 − 1 𝜁 ), 0 < 𝜁 < 1 τ ( 𝜁 ) ≤ τ ∗ (1 + log τ ( 𝜁 ) ≤ 2𝑓𝑢 ∗ (1 + log 1 𝜁 ), 0 < 𝜁 < 1

+ Entropy of initial distribution 61 Measure of randomness: 𝑓𝑜𝑢 𝜈 = − 𝜈 𝑦 log 𝜈 𝑦 𝑦∈𝑊 𝑦 is initial state 𝜈 is initial probability distribution

+ Few more lemma’s 62 1. Let (𝑌 𝑢 ) is a random walk associated with 𝜈 , and 𝜈 𝑢 be the distribution of 𝑌 𝑢 then 𝑓𝑜𝑢 𝜈 𝑢 ≤ 𝑢 ⋅ 𝑓𝑜𝑢(𝜈) 2. If 𝑤 is a distribution of 𝑇 such that 𝑤 − 𝜌 ≤ 𝜁 then 𝑓𝑜𝑢 𝑤 ≥ (1 − 𝜁) log |𝑇|

+ Lower bound on 𝑒(𝑢) proof 63  From lemma 2, 𝑓𝑜𝑢 𝑤 ≥ 1 − 𝜁 log 𝑇 𝑓𝑜𝑢 𝑤 log |𝑇| ≥ 1 − 𝜁 ≥ (1 − 𝑒(𝑢)) 𝑒 𝑢 ≥ 1 − 𝑓𝑜𝑢 𝑤 log |𝑇| 𝑢⋅𝑓𝑜𝑢 𝜈 𝑒 𝑢 ≥ 1 − log |𝑇|

+ Lower bound on τ ( ε ) proof 64  From lemma 1, 𝑓𝑜𝑢 𝜈 𝑢 ≤ 𝑢 𝑓𝑜𝑢(𝜈) 𝑓𝑜𝑢 𝜈 𝑢 𝑢 ≥ 𝑓𝑜𝑢(𝜈)  From lemma 2, 𝑓𝑜𝑢 𝜈 𝑢 ≥ (1 − 𝜁) log |𝑇| 𝑢 ≥ 1 − 𝜁 log |𝑇| 𝑓𝑜𝑢(𝜈) τ ( 𝜁 ) ≥ 1 − 𝜁 log |𝑇| 𝑓𝑜𝑢(𝜈)

+ 65 Lower bounds 𝑒 𝑢 ≥ 1 − 𝑢 ⋅ 𝑓𝑜𝑢(𝜈) log 𝑇 1 − 𝜁 log |𝑇| τ 𝜁 ≥ 𝑓𝑜𝑢(𝜈)

+ Propp Wilson Tobias Bertelsen 66

+ An exact version of MCMC 67  Problems with MCMC A. Have accuracy error, which depends on starting state We must know number of iterations B.  James Propp and David Wilson propos Coupling from the past (1996)  A.k.a. the Propp-Wilson algorithm  Idea:  Solve problems by running chain infinitely

+ An exact version of MCMC 68 Theoretical • Runs all configurations infinitely • Literarily takes ∞ time • Impossible Coupling from the past • Runs all configurations for finite time • Might take 1000’s of years • Infeasible Sandwiching • Run few configurations for finite time • Takes seconds • Practicable

+ Theoretical exact sampling 69  Recall the convergence theorem:  We will approach the stationary distribution as the number of steps goes to infinity  Intuitive approach:  To sample perfectly we start a chain and run for infinity  Start at 𝑢 = 0 , sample at 𝑢 = ∞  Problem: We never get a sample  Alternative approach:  To sample perfectly we take a chain that have already been running for an infinite amount of time  Start at 𝑢 = −∞ , sample at 𝑢 = 0

+ Theoretical independence of 70 starting state  Sample from a Markov chain in MCMC depends solely on  Starting state 𝑌 −∞  Sequence of random numbers 𝑉  We want to be independent of the starting state.  For a given sequence of random numbers 𝑉 −∞ , … , 𝑉 −1 we want to ensure that the starting state 𝑌 −∞ has no effect on 𝑌 0

+ Theoretical independence of 71 starting state Collisions:  For a given 𝑉 −∞ , … , 𝑉 −1 if two Markov chains is at the same state at some 𝑢′ the will continue on together:

+ Coupling from the past 72  At some finite past time 𝑢 = −𝑂  All past chains has already run infinitively and has coupled into one ∞ − 𝑂 = ∞  We want to continue that coupled chain to 𝑢 = 0  But we don’t know which at state they will be at 𝑢 = −𝑂  Run all states from −𝑂 instead of −∞

+ Coupling from the past 73 Let 𝑉 = 𝑉 −1 , 𝑉 −2 , 𝑉 −3 , … be the sequence of independnent 1. uniformly random numbers For 𝑂 𝑘 ∈ 1,2,4,8, … 2. Extend 𝑉 to length 𝑜 𝑘 , keeping 𝑉 −1 , … , 𝑉 𝑜 𝑘−1 the same 1. Start one chain from each state at 𝑢 = −𝑂 2. 𝑘 For 𝑢 from −𝑂 𝑘 to zero: Simulate the chains using 𝑉 𝑢 3. If all chains has converged at 𝑇 𝑗 at 𝑢 = 0 , return 𝑇 𝑗 4. Else repeat loop 5.

+ Coupling from the past 74

+ Questions 75  Why do we double the lengths?  Worst case < 4𝑂 𝑝𝑞𝑢 steps, where 𝑂 𝑝𝑞𝑢 is the minimal 𝑂 at which we can achieve convergence  Compare to N ∈ [1,2,3,4,5, … ] , 𝑃 𝑂 𝑝𝑞𝑢 2 steps  Why do we have to use the same random numbers?  Different samples might take longer or shorter to converge.  We must evaluate the same sample in each iteration.  The sample should only be dependent on 𝑉 not the different 𝑂

+ Using the same 𝑉 76  We have 𝑙 states with the following update function 1 1 𝑇 𝑗 , 𝑉 = 𝑇 1 𝑉 < 𝜚 𝑇 𝑙 , 𝑉 = 𝑇 1 𝑉 < 2 2 𝑇 𝑗+1 𝑝𝑢ℎ𝑓𝑠𝑥𝑗𝑡𝑓 𝑇 𝑙 𝑝𝑢ℎ𝑓𝑠𝑥𝑗𝑡𝑓 1 2  𝜌 𝑇 𝑗 = 2 𝑗 , 𝜌 𝑙 = 2 𝑙

+ 77 Using the same 𝑉  The probability of 𝑇 1 only depends on the last random number: 𝑄 𝑇 1 = 𝑄 𝑉 −1 < 1 = 1 2 2  Lets assume we generate a new 𝑉 for each run: 𝑉 1 , 𝑉 2 , 𝑉 3 , … 𝑸 𝑻 𝟐 acc. 𝒐 𝑽 𝑸 𝑻 𝟐 1 1 < 1 1 50 % 𝑉 −1 𝑄 𝑉 −1 2 2 , 𝑉 −1 2 2 < 1 1 < 1 2 𝑉 −2 75 % 𝑄 𝑉 −1 2 ∨ 𝑉 −1 2 3 , … , 𝑉 −1 3 3 < 1 2 < 1 1 < 1 4 81.25 % 𝑉 −4 𝑄 𝑉 −1 2 ∨ 𝑉 −1 2 ∨ 𝑉 −1 2 4 , … , 𝑉 −1 4 ] 3 < 1 2 < 1 1 < 1 4 < 1 8 [𝑉 −8 81.64 % 𝑄 𝑉 −1 2 ∨ 𝑉 −1 2 ∨ 𝑉 −1 2 ∨ 𝑉 1− 2

+ 78 Using the same 𝑉  The probability of 𝑇 1 only depends on the last random number: 𝑄 𝑇 1 = 𝑄 𝑉 −1 < 1 = 1 2 2  Lets instead use the same 𝑉 for each run: 𝑉 1 𝑸 𝑻 𝟐 acc. 𝒐 𝑽 𝑸 𝑻 𝟐 1 1 < 1 1 50 % 𝑉 −1 𝑄 𝑉 −1 2 1 , 𝑉 −1 1 1 < 1 2 𝑉 −2 50 % 𝑄 𝑉 −1 2 1 , … , 𝑉 −1 1 1 < 1 4 50 % 𝑉 −4 𝑄 𝑉 −1 2 1 , … , 𝑉 −1 1 ] 1 < 1 8 [𝑉 −8 50 % 𝑄 𝑉 −1 2

+ Problem 79  In each step we update up to 𝑙 chains  Total execution time is 𝑃 𝑂 𝑝𝑞𝑢 𝑙  BUT 𝑙 = 2 𝑀 2 in the Ising model  Worse than the naïve approach 𝑃(𝑙)

+ Sandwiching Nirandika Wanigasekara 80

+ Sandwiching 81 • Many vertices  running 𝑙 Propp Wilson chains will take time algorithm • Impractical for large 𝑙 Choose a • Can we still get the same relatively small results? state space • Try sandwiching

+ Sandwiching 82  Idea  Find two chains bounding all other chains  If we have such two boundary chains  Check if those two chains converge  Then all other chains have also converged

+ Sandwiching 83  To come up with the boundary chains  Need a way to order the states  𝑇 1 ≤ 𝑇 2 ≤ 𝑇 3 …  Chain in higher state does not cross a chain in lower state 2 2 2 2 2 2 1 1 1 1 1 1  Results in a Markov Chain obeying certain monotonicity properties

+ Sandwiching 84  Lets consider  A fixed set of states  𝑙  State space 𝑇 = {1, … . , 𝑙}  A transition matrix 12 = 1  𝑄 11 = 𝑄 2  𝑄 𝑙𝑙 = 𝑄 𝑙,𝑙−1 = 1 2  𝑔𝑝𝑠 𝑗 =2, ….. 𝑙 − 1, 𝑄 𝑗,𝑗−1 = 𝑄 𝑗,𝑗+1 = 1 2  All the other entries are 0

+ Sandwiching 85  What is this Markov Chain doing?  Take one step up and one step down the ladder, at each integer 1 time, with probability 2  If at the top or the bottom(state 𝑙 )  It will stays where it is  Ladder Walk on 𝒍 vertices  The stationary distribution 𝜌 of this Markov Chain 1  𝜌 𝑗 = 𝑙 𝑔𝑝𝑠 𝑗 = 1, … . , 𝑙

+ Sandwiching 86  Propp-Wilson for this Markov Chain with  Valid update function 𝜚  if u < 1 2 then step down  if 𝑣 ≥ 1 2 then step up  negative starting times 𝑂 1 , 𝑂 2 , … = 1, 2, 4, 8, …   States 𝑙 = 5

+ Sandwiching 87

+ Sandwiching 88  Update function preserves ordering between states  for all 𝑉 ∈ 0, 1 𝑏𝑜𝑒 𝑏𝑚𝑚 𝑗, 𝑘 ∈ 1, … . , 𝑙 𝑡𝑣𝑑ℎ 𝑢ℎ𝑏𝑢 𝑗 ≤ 𝑘 𝑥𝑓 ℎ𝑏𝑤𝑓 𝜚 𝑗, 𝑉 ≤ 𝜚(𝑘, 𝑉)  It is sufficient to run only 2 chains rather than k

+ Sandwiching 89  Are these conditions always met  No, not always  But, there are frequent instances where these conditions are met  Especially useful when k is large  Ising model is a good example for this

+ Ising Model Malay Singh 90

+ 2D Ising Model 91  A grid with sites numbered from 1 to 𝑀 2 where 𝑀 is grid size.  Each site 𝑗 can have spin 𝑦 𝑗 ∈ −1, +1  𝑊 is set of sites 𝑇 = −1,1 𝑊 defines all possible configurations (states)  Magnetisation 𝑛 for a state 𝑡 is 𝑗 𝑡 𝑦 𝑗  𝑛 𝑡 = 𝑀 2  The energy of a state 𝑡  H 𝑡 = − 𝑗𝑘 𝑡 𝑦 𝑗 𝑡 𝑦 𝑘

+ Ordering of states 92  For two states 𝑡, 𝑡 ′ we say s ≤ 𝑡 ′ if 𝑡 𝑦 < 𝑡 ′ 𝑦 ∀ 𝑦 ∈ 𝑊 Maxima 𝑛 = 1 Minima 𝑛 = −1 𝑡 min x = −1 ∀𝑦 ∈ 𝑊 𝑡max x = +1 ∀𝑦 ∈ 𝑊 Hence 𝑡 min ≤ 𝑡 ≤ 𝑡 max for all 𝑡

+ The update function 93 We use the sequence of random numbers 𝑉 𝑜 , 𝑉 𝑜−1 , … , 𝑉 0 We are updating the state at 𝑌 𝑜 . We choose a site 𝑦 in 0, 𝑀 2 uniformly. exp 2𝛾 𝐿 + 𝑦, 𝑡 − 𝐿 − 𝑦, 𝑡 +1, 𝑗𝑔 𝑉 𝑜+1 < 𝑌 𝑜+1 𝑦 = exp 2𝛾 𝐿 + 𝑦, 𝑡 − 𝐿 − 𝑦, 𝑡 + 1 −1, 𝑝𝑢ℎ𝑓𝑠𝑥𝑗𝑡𝑓

+ Maintaining ordering 94  Ordering after update (From 𝑌 𝑜 to 𝑌 𝑜+1 )  We choose the same site 𝑦 to update in both chains  The spin of 𝑦 at 𝑌 𝑜+1 depends on update function  We want to check exp 2𝛾 𝐿 + 𝑦, 𝑡 ′ − 𝐿 − 𝑦, 𝑡 ′ exp 2𝛾 𝐿 + 𝑦, 𝑡 − 𝐿 − 𝑦, 𝑡 ≤ exp 2𝛾 𝐿 + 𝑦, 𝑡 ′ − 𝐿 − 𝑦, 𝑡 ′ exp 2𝛾 𝐿 + 𝑦, 𝑡 − 𝐿 − 𝑦, 𝑡 + 1 + 1  That is equivalent to checking 𝐿 + 𝑦, 𝑡 − 𝐿 − 𝑦, 𝑡 ≤ 𝐿 + 𝑦, 𝑡 ′ − 𝐿 − 𝑦, 𝑡 ′

+ Maintaining ordering 95  As 𝑡 ≤ 𝑡 ′ we have  𝐿 + 𝑦, 𝑡 ≤ 𝐿 + 𝑦, 𝑡′  𝐿 − 𝑦, 𝑡 ≥ 𝐿 − 𝑦, 𝑡′  First equation minus second equation  𝐿 + 𝑦, 𝑡 − 𝐿 − 𝑦, 𝑡 ≤ 𝐿 + 𝑦, 𝑡′ − 𝐿 − 𝑦, 𝑡′

+ Ising Model L = 4 96 𝑈 = 3.5 and 𝑂 = 512 𝑈 = 4.8 and 𝑂 = 128

+ Ising Model 𝑀 = 8 97 𝑈 = 5.9 and 𝑂 = 512

+ Ising Model 𝑀 = 16 98 𝑈 = 5.3 and 𝑂 = 16 384

+ Summary 99  Exact sampling  Markov chains monte carlo but when we converge  Propp Wilson with Sandwiching to rescue

+ Questions? 100