2010 June « Journey into Randomness

Is this random ? Card shuffling, coins and their friends …

June 23, 2010 at 5:57 pm (probability, random permutation)
Tags: card shuffling, coins, persi diaconis, poker, randomness

A brilliant talk by Persi Diaconis !

1 Comment

Potts model and Monte Carlo Slow Down

June 22, 2010 at 12:24 am (markov chain, Monte Carlo, probability)
Tags: Ising, MCMC, Metropolis, Monte Carlo, Potts model

A simple model of interacting particles

The mean field Potts model is extremely simple: there are ${N}$ interacting particles ${x_1, \ldots, x_N}$ and each one of them can be in ${q}$ different states ${1,2, \ldots, q}$ . Define the Hamiltonian

$\displaystyle H_N(x) = -\frac{1}{N} \sum_{i,j} \delta(x_i, x_j)$

where ${x=(x_1, \ldots, x_N)}$ and ${\delta}$ is the Kronecker symbol. The normalization ${\frac{1}{N}}$ ensures that the energy is an extensive quantity so that the mean energy per particle ${h_N(x) = \frac{1}{N} H_N(x)}$ does no degenerate to ${0}$ or ${+\infty}$ for large values of ${N}$ . The sign minus is here to favorize configurations that have a lot of particles in the same state. The Boltzman distribution at inverse temperature ${\beta}$ on ${\{1, \ldots, q\}^N}$ is given by

$\displaystyle P_{N,\beta} = \frac{1}{Z_N(\beta)} e^{-\beta H_N(x)}$

where ${Z_N(\beta)}$ is a normalization constant. Notice that if we choose a configuration uniformly at random in ${\{1, \ldots, q\}^N}$ , with overwhelming probability the ratio of particles in state ${k}$ will be close to ${\frac{1}{q}}$ . Also it is obvious that if we define

$\displaystyle L^{(N)}_k(x) = \frac{1}{N} \, \Big( \textrm{Number of particles in state }k \Big)$

then ${L=(L^{(N)}_1, \ldots, L^{(N)}_q)}$ will be close to ${(\frac{1}{q}, \ldots, \frac{1}{q})}$ for a configuration taken uniformly at random. Stirling formula even says that the probability that ${L}$ is close to ${\nu = (\nu_1, \ldots, \nu_q)}$ is close to ${e^{-N \, R(\nu)}}$ where

$\displaystyle R(\nu) = \nu_1 \ln(q\nu_1) + \ldots + \nu_q \ln(q\nu_q).$

Indeed ${(\frac{1}{q}, \ldots, \frac{1}{q}) = \textrm{argmin} \, R(\nu)}$ . The situation is quite different under the Boltzman distribution since it favorizes the configurations that have a lot of particles in the same state: this is because the Hamiltonian ${H_N(x)}$ is minimized for configurations that have all the particles in the same state. In short there is a competition between the entropy (there are a lot of configurations close to the ratio ${(\frac{1}{q}, \ldots, \frac{1}{q})}$ ) and the energy that favorizes the configurations where all the particles are in the same state.

With a little more work, one can show that there is a critical inverse temperature ${\beta_c}$ such that:

for ${\beta < \beta_c}$ the entropy wins the battle: the most probable configurations are close to the ratio ${(\frac{1}{q}, \ldots, \frac{1}{q})}$
for ${\beta > \beta_c}$ the energy effect shows up: there are ${q}$ most probable configurations that are the permutations of ${(a_{\beta},b_{\beta},b_{\beta}, \ldots, b_{\beta})}$ where ${a_{\beta}}$ and ${b_{\beta}}$ are computable quantities.

The point is that above ${\beta_c}$ the system has more than one stable equilibrium point. Maybe more important, if we compute the energy of these most probable states

$\displaystyle h(\beta) = \lim \frac{1}{N} H_N(\textrm{most probable state})$

then this function has a discontinuity at ${\beta=\beta_c}$ . I will try to show in the weeks to come how this behaviour can dramatically slow down usual Monte-Carlo approach to the study of these kind of models.

Hugo Touchette has a very nice review of statistical physics that I like a lot and a good survey of the Potts model. Also T. Tao has a very nice exposition of related models. The blog of Georg von Hippel is dedicated to similar models on lattices, which are far more complex that this mean field approximation presented here.

MCMC Simulations

These is extremely easy to simulate this mean field Potts model since we only need to keep track of the ratio ${L=(L_1, \ldots, L_q)}$ to have an accurate picture of the system. For example, a typical Markov Chain Monte Carlo approach would run as follows:

choose a particle ${x_i}$ uniformly at random in ${\{1,2, \ldots, N\}}$
try to switch its value uniformly in ${\{1,2, \ldots, q\} \setminus \{x_i\}}$
compute the Metropolis ratio
update accordingly.

If we do that ${10^5}$ times for ${q=3}$ states at inverse temperature ${\beta=1.5}$ and for ${100}$ particles (which is fine since we only need to keep track of the ${3}$ -dimensional ratio vector) and plot the result in barycentric coordinates we get a picture that looks like:

Here I started with a configuration where all the particles were in the same states i.e ratio vector equal to ${(1,0,0)}$ . We can see that even with ${10^5}$ steps, the algorithm struggles to go from one most probable position ${(a,b,b)}$ to the other two ${(b,a,b)}$ and ${(b,b,a)}$ – in this simulation, one of the most probable state has even not been visited! Indeed, this approach was extremely naive, and this is quite interesting to try to come up with better algorithms. Btw, Christian Robert’s blog has tons of interesting stuffs related to MCMC and how to boost up the naive approach presented here.

Challenge – Announcement

June 6, 2010 at 12:44 pm (17x17 challenge, Monte Carlo, probability)
Tags: 17x17 challenge, MCMC, Monte Carlo, optimization, simulated annealing, William Gasarch

Have you already heard about the 17×17 Challenge ? It all started 7 months ago, and still no solution has been found. I have the feeling though that no really serious attempt at solving the problem has been done, even if very encouraging results have been found. It could be very interesting to keep track of the progress made so that results will be published here. Send what you have reached so far!

In the next few weeks, I will try to write about the Monte Carlo approaches described in the excellent book “Monte Carlo Strategies in Scientific Computing” by Jun Liu: I do not think that one can easily solve the problem with a direct implementation of these methods, but I want to know more about it. These are all more or less Markov Chain Monte Carlo Methods: I have tried the “parallel tempering” method, which gives results that are not that bad, and was planning to check the “flat histogram” method. “Simulated annealing” methods have succesfully been used to solve the Sudoku problem though: see Christian Robert blog for example. I am also planning to continue my reading of the fascinating book “Information, Physics, and Computation” by Marc Mezard: methods described in this book (belief propagation, k-SAT problem, etc…) are more adapted to the 17×17 challenge I think, and I am quite optimistic about it.

Best Solution so far, by Rohan Puttagunta: only 4 monochromatic rectangles!

Read also here,here,here,here and there – and have a try with this on-line simulator !

Random permutations and giant cycles

June 3, 2010 at 11:07 pm (coagulation, fluid limit, fragmentation, markov chain, probability, random permutation)
Tags: Erdos Renyi, giant cycle, markov chain, phase transition, random graph, random permutation

The other day, Nathanael Berestycki presented some very nice results on random permutations. Start with the identity permutation of ${\mathfrak{S}_N}$ and compose with random transpositions: in short, at time ${k}$ we consider the random permutation ${\sigma_k = \tau_k \circ \ldots \tau_{2} \circ \tau_1}$ where ${\tau_i}$ is a transposition chosen uniformly at random among the ${\frac{N(N-1)}{2}}$ transpositions of ${\mathfrak{S}_N}$ . How long do we have to wait before seeing a permutation that has a cycle of length greater than ${\epsilon N}$ , where ${\epsilon>0}$ is a fixed constant ?

The answer has been known for quite a long time (O. Schramm,2005), but no simple proof was available: the highlight of Nathanael’s talk was a short and elegant proof of the fact that we only have to wait a little bit more than ${N/2}$ steps to see this giant cycle! The first half of the proof is such a beauty that I cannot resist in sharing it here!

1. Visual representation

First, this is more visual to represent a permutation by its cycle decomposition. Consider for example the permutation ${\sigma = (1,3,4)(5,6,8,9,10,11)(2)(7) \in \mathfrak{S}_{10}}$ that sends ${1 \rightarrow 3 \rightarrow 4 \rightarrow 1}$ etc… It can be graphically represented by the following pictures:

cycle structure - visual representation

If this permutation is composed with the transposition ${(4,6)}$ , this gives the following cycle structure:

two cycles merge

If it is composed with ${(8,11)}$ this gives,

a cycles breaks down into two

It is thus pretty clear that the cycle structure of ${\sigma_{k+1}}$ can be obtained from the cycle structure of ${\sigma_k}$ in only 2 ways:

or two cycles of ${\sigma_k}$ merge into one,
or one cycle of ${\sigma_k}$ breaks into two cycles.

In short, the process ${\{\sigma_k\}_k}$ behaves like a fragmentation-coalescence processes. At the beginning, the cycle structure of ${\sigma_0=\text{Id}}$ is trivial: there are ${N}$ trivial cycles of length one! Then, as times goes by, they begin to merge, and sometimes they also break down. At each step two integers ${i}$ and ${j}$ are chosen uniformly at random (with ${i \neq j}$ ): if ${i}$ and ${j}$ belong to the same cycle, this cycle breaks into two parts. Otherwise the cycle that contains ${i}$ merges with the cycle that contains ${j}$ . Of course, since at the beginning all the cycles are extremely small, the probability that a cycles breaks into two part is extremely small: a cycle of size ${k}$ breaks into two part with probability roughly equal to ${\Big(\frac{k}{N}\Big)^2}$ .

2. Comparison with the Erdos-Renyi random Graph

Since the fragmentation part of the process can be neglected at the beginning, this is natural to compare the cycle structure at time ${k}$ with the rangom graph ${G(N,p_k)}$ where ${\frac{N(N-1)}{2}p_k = k}$ : this graph has ${N}$ vertices and roughly ${k}$ edges since each vertex is open with probability ${p_k}$ . A coupling argument makes this idea more precise:

start with ${\sigma_0=\textrm{Id}}$ and the graph ${G_0}$ that has ${N}$ disconnected vertices ${\{1,2, \ldots, N\}}$
choose two different integers ${i \neq j}$
define ${\sigma_{k+1}=(i,j) \circ \sigma_k}$
to define ${G_{k+1}}$ , take ${G_k}$ and join the vertices ${i}$ and ${j}$ (never mind if they were already joined)
go to step ${2}$ .

The fundamental remark is that each cycle of ${\sigma_k}$ is included in a connected component of ${G_k}$ : this is true at the beginning, and this property obviously remains true at each step. The Erdos-Renyi random graph ${G(n,p)}$ is one of the most studied object: it is well-known that for ${c<1}$ , with probability tending to ${1}$ , the size of the largest connected component of ${G(N,\frac{c}{N})}$ is of order ${\ln(N)}$ . Since ${p_k \sim \frac{2k}{N^2}}$ , this immediately shows that for ${k \approx \alpha N}$ , with ${\alpha < \frac{1}{2}}$ the length of the longuest cycle in ${\sigma_k}$ is of order ${\ln(N)}$ or less. This is one half of the result. It remains to work a little bit to show that if we wait a little bit more than ${\frac{N}{2}}$ , the first cycle of order ${N}$ appears.

For a permutation ${\sigma \in \mathfrak{S}_N}$ with cycle decomposition ${\sigma = \prod_{j=1}^r c_i}$ (cycles of length ${1}$ included), one can consider the quantity

$\displaystyle T(\sigma) = \sum_{j=1}^R |c_j|^2.$

This is a quantity between ${N}$ (for the identity) and ${N^2}$ that measures in some extends the sizes of the cycles. In order to normalize everything and to take into account that the cycle structure abruptly changes near the critical value ${\frac{N}{2}}$ , one can also observe this quantity on a logarithmic scale:

$\displaystyle s(\sigma) = \frac{ \ln\Big( T(\sigma)/N^2 \Big) }{\ln(N)}.$

This quantity ${s(\sigma)}$ lies between ${0}$ and ${1}$ , and evolves quite smoothly with time: look at what it looks like, once the time has been rescaled by a factor ${N}$ :

fluid limit ?

The plot represented the evolution of the quantity ${s\Big(\sigma(\alpha N)\Big)}$ between time $k=0$ and $k=\frac{N}{2}$ for dimensions $N=10000,20000,30000$ .

Question: does the random quantity $s(\alpha,N)={s\Big(\sigma(\alpha N)\Big)}$ converge (in probability) towards a non trivial constant ${S(\alpha) \in (0,1)}$ ? It must be the analogue quantity for the random graph ${G(N,\frac{2\alpha}{N})}$ , so I bet it is a known result …

1 Comment

Doob H-transforms

June 2, 2010 at 11:21 am (Brownian motion, Markov process, probability, Stochastic analysis)
Tags: brownian motion, conditioning, Doob, h-transform, Markov process, probability

I read today about Doob h-transforms in the Rogers-Williams … It is done quite quickly in the book so that I decided to practice on some simple examples to see how this works.

So we have an Markov process ${X_t}$ living in ${S}$ , and we want to see how this process looks like if we condition on the event ${X_T \in A}$ , say. To fix the notations we define ${p(t,t+s,x,y) = P(X_{t+s}=y|X_t=x)}$ and ${h(t,x)=P(X_T \in A \, | X_t=x)}$ . The conditioned semi-group ${\hat{p}(t,t+s,x,y)=P(X_{t+s}=y|X_t=x, X_T \in A)}$ is quite easily computed from ${p}$ and ${h}$ . Indeed, this is also equal to

$\displaystyle \hat{p}(t,t+s,x,y) = \frac{P(X_{t+s}=y; X_T \in A\;|X_t=x)}{P(X_T \in A \,|X_t=x)} = p(t,t+s,x,y) \frac{h(t+s,y)}{h(t,x)}.$

Notice also that ${\hat{p}(t,t+s,x,y) = p(t,t+s,x,y) \frac{h(t+s,y)}{h(t,x)}}$ is indeed a Markov kernel in the sense that ${\int_{y} \hat{p}(t,t+s,x,y) \, dy = 1}$ : the only property needed for that is

$\displaystyle h(t,x) = \int_{y} p(t,t+s,x,y)h(t+s,y)\,dy = E\left[ h(t+s,X_{t+s}) \, |X_t=x\right].$

In fact, we could take any function ${h}$ that satisfies this equality and define a new Markovian kernel ${\hat{p}}$ and study the associated Markov process. That’s what people usually do by the way.

Remark 1 we almost never know explicitly the quantity ${h(t,x)}$ , except in some extremely simple cases !

Before trying these ideas on some simple examples, let us see what this says on the generator of the process:

continuous time Markov chains, finite state space:let us suppose that the intensity matrix is ${Q}$ and that we want to know the dynamic on ${[0,T]}$ of this Markov chain conditioned on the event ${X_T=z}$ . Indeed ${p(t,t+s,i,j) = [\exp(sQ)]_{i,j}}$ so that ${\hat{p}(t,t+s,i,j) = [\exp(sQ)]_{i,j} \frac{p(t+s,T,j,z)}{p(t,T,i,z)}}$ so that in the limit we see that at time ${t}$ , the intensity of the jump from ${i}$ to ${j}$ of the conditioned Markov chain is
$\displaystyle Q(i,j) \frac{p(t+s,T,j,z)}{p(t,T,i,z)}.$

Notice how this behaves while ${t \rightarrow T}$ : if at ${t=T-\epsilon}$ the Markov chain is in state ${i \neq z}$ then the intensity of jump from ${i}$ to ${z}$ is equivalent to ${\approx \frac{1}{\epsilon}}$ .
diffusion processes:this time consider a ${1}$ -dimensional diffusion ${dX_t = \mu(X_t) \, dt + \sigma(X_t) \, dW_t}$ on ${[0,T]}$ conditioned on the event ${X_T \in A}$ and define as before ${h(t,x)=P(X_T \in A \,|X_t=x)}$ . The generator of the (non-homogeneous) conditioned diffusion is defined at time ${t}$ by
$\displaystyle \begin{array}{rcl} \mathcal{G}^{(t)} f(x) &=& \lim_{s \rightarrow 0} \frac{1}{s} \Big( E\left[f(X_{t+s}) \,| X_t=x, X_T \in A\right]-f(x) \Big)\\ &=& \lim_{s \rightarrow 0} \frac{E\left[ f(X_{t+s}) h(t+s, X_{t+s}) \,| X_t=x\right]-f(x) }{s\,h(t,x)} \end{array}$

so that if ${\mathcal{L} = \mu \partial_x + \frac{1}{2} \sigma^2 \partial^2_{xx}}$ is the generator of the original diffusion we get

$\displaystyle \mathcal{G}^{(t)} f = \frac{1}{h} \Big(\partial_t + \mathcal{L})(hf).$

Because ${(\partial_t + \mathcal{L})h=0}$ , this also reads

$\displaystyle \mathcal{G}^{(t)} f = \mathcal{L}f + \sigma^2 \frac{\partial_x h}{h} \partial_x f.$

This means that the conditioned diffusion ${Z}$ follows the SDE:

$\displaystyle dZ_t = \Big( \mu(Z_T) + \sigma(Z_t)^2 \frac{\partial_x h(t,Z_t)}{h(t,Z_t)} \Big) \, dt+ \sigma(Z_t) dW_t.$

The volatility function remains the same while an additional drift shows up.

We will try these ideas on some examples where the probability densities are extremely simple. Notice that in the case of diffusions, if we take ${A=\{x^+\}}$ , the function ${(t,x) \mapsto P(X_T = x^+ \, |X_t=x)}$ is identically equal to ${0}$ (except degenerate cases): to condition on the event ${X_T=x}$ we need instead to take ${h(t,x)}$ to be the transition probability ${p(t,T,x,x^+)}$ . This follows from the approximation ${P(X_T \in (x^+,x^+ + dx) \,|X_t=x ) \approx p(t,T,x,x^+) \, dx + o(dx)}$ . Let’s do it:

Brownian Bridge on ${[0,T]}$ :in this case ${p(t,T,x,0) \propto e^{-\frac{|x-y|^2}{2(T-t)}}}$ so that the additional drift reads ${\frac{-x}{T-t}}$ : a Brownian bridge follows the SDE
$\displaystyle dX_t = -\frac{X_t}{T-t} \, dt + dW_T.$

This might not be the best way to simulate a Brownian bridge though!
Poisson Bridge on ${[0,T]}$ :we condition a Poisson process of rate ${\lambda}$ on the event ${X_T=N}$ . The intensity matrix is simply ${Q(k,k+1)=\lambda=-Q(k,k)}$ and ${0}$ everywhere else while the transition probabilities are given by ${p(t,T,k,N) = e^{-\lambda (T-t)} \frac{(\lambda (T-t) )^{N-k}}{(N-k)!}}$ . This is why at time ${t}$ , the intensity from ${k}$ to ${k+1}$ is given by
$\displaystyle \lambda(t,k,k+1) = \frac{N-k}{T-t}.$

Again, that might not be the most efficient way to simulate a Poisson Bridge ! Notice how the intensity ${\lambda}$ has disappeared …
Ornstein-Uhlenbeck Bridge:Let’s consider the usual OU process given by the dynamic ${dX_t = -X_t + \sqrt{2}dW_t}$ : the invariant probability is the usual centred Gaussian distribution. Say that we want to know how does such an OU process behave if we condition on the event ${X_T = z}$ . Because ${p(t,T,x,z) \propto \exp(-\frac{|z-e^{-(T-t)}x|^2}{2(1-e^{-2(T-t)})^2})}$ we find that the conditioned O-U process follows the SDE
$\displaystyle dX_t = \Big(\frac{z-e^{-(T-t)x}}{e^{T-t}(1-e^{-2(T-t)})^2} - X_t \Big)\, dt+ \sqrt{2} \, dW_t.$

If we Taylor expand the additonal drift, it can be seen that this term behaves exactly as in the case of the Brownian bridge. Below is a plot of an O-U process conditioned on the event ${X_{10} = 10}$ , starting from ${X_0=0}$ .

conditioned O-U process