# Path Integrals - Lévy Sampling¶

Doruk Efe Gökmen -- 19/08/2018 -- Ankara

## Quantum statistical mechanics - Density matrices¶

In a thermal ensemble, the probability of being in $n$th energy eigenstate is given by the Boltzmann factor $\pi(n)\propto e^{-\beta E_n}$, where $\beta=\frac{1}{k_BT}$. Hence, e.g the probability $\pi(x,n)$ to be in state $n$ and in position $x$ is proportional to $e^{-\beta E_n}|\psi_n(x)|^2$.

We can consider the diagonal density matrix $\rho(x,x,\beta)=\sum_n e^{\beta E_n}\psi_n(x)\psi_n^*(x)$, which is the probability $\pi(x)$ of being at position $x$. This is a special case of the more general density matrix $\rho(x,x',\beta)=\sum_n e^{\beta E_n}\psi_n(x)\psi_n^*(x')$, which is the central object of quantum statistical mechanics. The partition function is given by $Z(\beta)=\text{Tr}\rho_u=\int_{-\infty}^\infty \rho_u(x,x,\beta)\text{d}x$, where $\rho_u=e^{-\beta \mathcal{H}}$ is the unnormalised density matrix. It follows that $\rho(\beta)=\frac{e^{-\beta\mathcal{H}}}{\text{Tr}(e^{-\beta\mathcal{H}})}$.

Properties of the density matrix:

• The convolution property: $\int \rho(x,x',\beta_1) \rho(x',x'',\beta_2) \text{d}x' = \int \text{d}x' \sum_{n,m} \psi_n(x)e^{-\beta_1 E_n} \psi_n^*(x')\psi_m(x')e^{-\beta_2 E_m}\psi_m^*(x'')$ $= \sum_{n,m} \psi_n(x)e^{-\beta_1 E_n} \int \text{d}x' \psi_n^*(x')\psi_m(x')e^{-\beta_2 E_m}\psi_m^*(x'') = \sum_n \psi_n(x)e^{-(\beta_1+\beta_2)E_n}\psi_n^*(x'')=\rho(x,x'',\beta_1+\beta_2)$ $\implies \boxed{ \int \rho(x,x',\beta) \rho(x',x'',\beta) \text{d}x' = \rho(x,x'',2\beta)}$ (note that in the discrete case, this is just matrix squaring). So, if we have the density matrix at temperature $T=k_B/\beta$ this equation allows us to compute the density matrix at temperature $T/2$.

• The free density matrix for a system of infinte size is $\rho^\text{free}(x,x',\beta)=\frac{1}{\sqrt{2\pi\beta}}\exp{\left[-\frac{(x-x')^2}{2\beta}\right]}$. Notice that in the high temperature limit ($\beta\rightarrow 0$) the density matrix becomes classical: $\rho^\text{free}\rightarrow \delta(x-x')$. The quantum system exihibits its peculiar properties more visibly at low temperatures.

• High temperature limit and the Trotter decomposition. In general any Hamiltonian can be written as $\mathcal{H}=\mathcal{H}^\text{free}+V(x)$. At high temperatures ($\beta\rightarrow 0$) we can approximate the density matrix as $\rho(x,x',\beta)=e^{-\beta V(x)/2}\rho^\text{free}e^{-\beta V(x')/2}$ (Trotter expansion). Hence an explicit expression for the density matrix is available without solving the Schrödinger (or more preciesly Liouville) equation for any potential.

## Path integrals - Quantum Monte Carlo¶

### Path integral representation of the kernel¶

The kernel $K$ is the matrix element of the unitary time evolution operator $U(t_i-t_f)=e^{-i/\hbar(t_f-t_i)\mathcal{H}}$ in the position representation: $K(x_i,x_f;t_f-t_i)=\langle x_f \left| U(t_f-t_i) \right| x_i \rangle$. We can write $K(x_i,x_f;t_f-t_i)=\langle x_f \left| U^N((t_f-t_i)/N) \right| x_i \rangle$, that is, divide the time interval $[t_i,t_f]$ into $N$ equal intervals $[t_k,t_{k+1}]$ of length $\epsilon$, where $\epsilon=t_{k+1}-t_k=(t_f-t_i)/N$.

Then we can insert $N-1$ resolutions of identity ($\int_{-\infty}^\infty \text{d} x_k \left|x_k\rangle\langle x_k\right|$) to obtain

$K(x_i,x_f;t_f-t_i)= \left[\Pi_{k=1}^{N-1}\int_{-\infty}^\infty dx_k \right] \left[\Pi_{k=0}^{N-1} K(x_i,x_f;\epsilon = (t_f-t_i)/N)\right]$,

where $x_f=x_N$ and $x_i=x_0$. In the continuous limit, we would have

$K(x_i,x_f;t_f-t_i)= \lim_{N\rightarrow\infty} \left[\Pi_{k=1}^{N-1}\int_{-\infty}^\infty dx_k \right] \left[\Pi_{k=0}^{N-1} K(x_i,x_f;\epsilon = (t_f-t_i)/N)\right]$. (A)

Let us now consider the limit $\epsilon\rightarrow 0$ ($N\rightarrow \infty$) to obtain the short time kernel $K(x_i,x_f;\epsilon)$ and thereby switching from discrete to the continuous limit. It is known that for small $\epsilon$ the Trotter formula implies that to a very good approximation

$K(x_i,x_f;\epsilon = (t_f-t_i)/N) \simeq \langle x_{k+1} \left| e^{-i(\hbar\epsilon T} e^{-i/\hbar \epsilon V} \right| x_k\rangle$,

which becomes exact as $\epsilon\rightarrow 0$. If we insert resolution of identity $\int \text{d}p_k \left| p_k \rangle\langle p_k \right|$, we get

$K(x_i,x_f;\epsilon) = \int_{-\infty}^\infty \text{d}p_k \langle x_{k+1} \left| e^{-i(\hbar\epsilon T} \left| p_k \rangle\langle p_k \right| e^{-i/\hbar \epsilon V} \right| x_k\rangle = \int_{-\infty}^\infty \text{d}p_k \langle x_{k+1} \left| p_k \rangle\langle p_k \right| x_k\rangle e^{-i/\hbar \epsilon \left(\frac{p_k}{2m} + V(x)\right)}$

$\implies K(x_i,x_f;\epsilon) = \frac{1}{2\pi \hbar}\int_{-\infty}^\infty \text{d}p_k e^{i/\hbar \epsilon \left[p_k\frac{x_{k+1}-x_k}{\epsilon}-\mathcal{H}(p_k,x_k) \right]}$. (B)

Hence, inserting (B) into (A) we get

$K(x_i,x_f;t_f-t_i) = \lim_{N\rightarrow \infty}\left[\Pi_{k=1}^{N-1}\int_{-\infty}^\infty dx_k \right] \left \{ \Pi_{k=0}^{N-1} \int_{-\infty}^\infty \text{d}p_k e^{i/\hbar \epsilon \left[p_k\frac{x_{k+1}-x_k}{\epsilon}-\mathcal{H}(p_k,x_k) \right]} \right\}$. (C)

We can simplify the exponent of the integrand in the limiting case $N\rightarrow \infty$,

$\lim_{N\rightarrow \infty} \epsilon \sum_{k=0}^{N-1}\left[p_k\frac{x_{k+1}-x_k}{\epsilon}-\mathcal{H}(p_k,x_k) \right] =\int_{t_1}^{t_2}\text{d}t[p(t)\dot{x}(t)-\mathcal{H}[p(t),x(t)]]$

$=\int_{t_1}^{t_2}\text{d}t \mathcal{L}[x(t),\dot{x}(t)] = \mathcal{S}[x(t);t_f,t_i]$, (D)

where $\mathcal{L}[x(t),\dot{x}(t)] = \frac{m}{2}\dot{x}(t)^2-V[x(t)]$ is the Lagrangian and $\mathcal{S}[x(t);t_f,t_i]$ is the action between times $t_f$ and $t_i$.

Furthermore we can introduce the following notation for the integrals over paths:

$\lim_{N\rightarrow \infty}\left(\Pi_{k=1}^{N-1} \int_{-\infty}^\infty \text{d}x_k\right)=\int_{x(t_i)=x_i}^{x(t_f)=x_f}\mathcal{D}[x(t)]$, (E.1)

$\lim_{N\rightarrow \infty}\left(\Pi_{k=1}^{N-1}\int_{-\infty}^\infty\frac{\text{d}p_k}{2\pi\hbar}\right) =\int \mathcal{D}\left[\frac{p(t)}{2\pi\hbar}\right]$. (E.2)

Using (D) and (E) in (C), we get the path integral representation of the kernel

$K(x_i,x_f;t_f-t_i)= \int_{x(t_i)=x_i}^{x(t_f)=x_f}\mathcal{D}[x(t)] \int \mathcal{D}\left[\frac{p(t)}{2\pi\hbar}\right] e^{i/\hbar \mathcal{S}[x(t)]}$

$\implies \boxed{K(x_i,x_f;t_f-t_i)= \mathcal{N} \int_{x(t_i)=x_i}^{x(t_f)=x_f}\mathcal{D}[x(t)] e^{i/\hbar \mathcal{S}[x(t)]}}$, (F)

where $\mathcal{N}$ is the normalisation factor.

Here we see that each path has a phase proportional to the action. The equation (F) implies that we sum over all paths, which in fact interfere with one another. The true quantum mechanical amplitude is determined by the constructive and destructive interferences between these paths. For example, actions that are very large compared to $\hbar$, lead to very different phases even between nearby paths that differ only slightly, and that causes destructive interference between them. Only in the extremely close vicinity of the classical path $\bar x(t)$, where the action changes little when the phase varies, will neighbouring paths contirbute to the interference constructively. This leads to a classical deterministic path $\bar x(t)$, and this is why the classical approximation is valid when the action is very large compared to $\hbar$. Hence we see how the classical laws of motion arise from quantum mechanics.

### Path integral representation of the partition function¶

Heuristic derivation of the discrete case: Recall the convolution property of the density matrix, we can apply it repeatedly:

$\rho(x_0,x_2,\beta) = \int \rho(x_0,x_2,\beta/2) \rho(x_2,x_1,\beta/2) \text{d}x_2 = \int \int \int \rho(x_0,x_3,\beta/4) \rho(x_3, x_2,\beta/4) \rho(x_2,x_4,\beta/4) \rho(x_4,x_1 ,\beta/4) \text{d}x_2 \text{d}x_3 \text{d}x_4 = \cdots$

In other words: $\rho(x_0,x_N,\beta) = \int\int \cdots \int \text{d}x_1 \text{d}x_2 \cdots \text{d}x_{N-1}\rho(x_0,x_1,\beta/N)\rho(x_1,x_2,\beta/N)\cdots\rho(x_{N-1},x_N,\beta/N)$. The variables $x_k$ in this integral is called a path. We can imagine the variable $x_k$ to be at position $x_k$ at given slice $k\beta/N$ of an imaginary time variable $\tau$ that goes from $0$ to $\beta$ in steps of $\Delta\tau=\beta/N$. Density matrices and partition functions can thus be expressed as multiple integrals over path variables, which are none other than the path integrals that were introduced in the previous subsection.

Given the unnormalised density matrix $\rho_u$, the discrete partition $Z_d(\beta)$ function can be written as a path integral for all closed paths (because of the trace property), i.e., paths with the same beginning and end points ($x_0=x_N$), over a “time” interval $−i\hbar\beta$.

$Z_d(\beta)= \text{Tr}(e^{-\beta \mathcal{H}}) = \text{Tr}(\rho_u(x_0,x_N,\beta) )=\int \text{d}x_0 \rho_u (x_0,x_N=x_0,\beta)$ $= \int \int\int \cdots \int \text{d}x_0 \text{d}x_1 \text{d}x_2 \cdots \text{d}x_{N-1}\rho_u(x_0,x_1,\beta/N)\rho_u(x_1,x_2,\beta/N)\cdots\rho_u(x_{N-1},x_N,\beta/N)\rho_u(x_{N-1},x_0,\beta/N)$.

The integrand is the probabilistic weight $\Phi\left[\{x_i\}\right]$ of the discrete path consisting of points $\{x_i\}$. The continuous case can be obtained by taking the limit $N\rightarrow \infty$. By defining

$\Phi[x(\tau)] = \lim_{N\rightarrow \infty} \rho_u(x_0,x_1,\beta/N)\cdots \rho_u(x_{N-1},x_N,\beta/N)$, (G)

(note that this is the probability weight of a particular continuous path), and by using (E.1), we can express the continuous partition function $Z(\beta)$ as

$Z(\beta) = \int_{x(0)}^{x(\hbar \beta)=x(0)}\mathcal{D}[x(\tau)] \Phi[x(\tau)]$. (H)

But what is $\Phi[x(\tau)]$?

Derivation of the continuous case: Again we start from $Z(\beta)= \text{Tr}(e^{-\beta \mathcal{H}})$. The main point of the argument that follows is the operational resemblance between the unitary time evolution operator $U(t)=e^{-(i/\hbar) t\mathcal{H}}$ and the unnormalised density matrix $e^{-\beta \mathcal{H}}$: the former is used to define the kernel which reads $K(x,x';t)=\langle x \left| e^{-(i/\hbar) t\mathcal{H}} \right| x' \rangle$; and the latter is used in defining the density matrix which reads $\rho(x,x';\beta)=\langle x \left| e^{-\beta \mathcal{H}} \right| x' \rangle$. If we regard $\beta$ as the analytic continuation of the real time $t$ to the imaginary values: $t\rightarrow i \tau \rightarrow i \hbar \beta$, and $t=t_i-t_f$, we get the cousin of the partition function that lives in the imaginary spacetime (i.e. Minkowskian rather than Euclidian, but this has nothing to do with relativity, or has it?)

$Z\left[\beta\rightarrow -\frac{i}{\hbar}(t_f-t_i)\right]=\text{Tr}\left[U(t_f-t_i)\right]=\int_{-\infty}^\infty \text{d}x \langle x \left| U(t_f-t_i) \right| x \rangle$

$=\int_{-\infty}^\infty \text{d}x K(x,x;t_f-t_i)$

$=\int_{-\infty}^\infty \text{d}x \mathcal{N} \int_{x(t_i)=x}^{x(t_f)=x}\mathcal{D}[x(t)] e^{i/\hbar \int_{t_i}^{t_f}\text{d}t \mathcal{L}[x(t),\dot{x}(t)]}$ (using (F))

$=\mathcal{N} \int_{x(t_f)=x(t_i)}\mathcal{D}[x(t)] e^{i/\hbar \int_{t_i}^{t_f}\text{d}t \mathcal{L}[x(t),\dot{x}(t)]} = \mathcal{N} \int_{x(t_f)=x(t_i)}\mathcal{D}[x(t)] e^{i/\hbar \int_{t_i}^{t_f}\text{d}t \left[\frac{m}{2}\dot{x}(t)^2-V[x(t)]\right]}$,

which means that one is integrating not over all paths but over all closed paths (loops) at $x$. We are now ready to get the path integral representation of the real partition function by making the transformations $t\rightarrow i\tau$ so that $t_i\rightarrow 0$ and $t_f\rightarrow -i\hbar \beta$ (also note that $\dot{x(t)}=\frac{\partial x(t)}{\partial t}\rightarrow -i \frac{\partial x(\tau)} {\partial \tau} = -i x'(\tau) \implies \dot{x}(t)^2 \rightarrow -x'(\tau)^2$):

$\implies Z(\beta)=\mathcal{N} \int_{x(\hbar \beta)=x(0)}\mathcal{D}[x(\tau)] e^{-\frac{1}{\hbar} \int_{0}^{\beta \hbar}\text{d}\tau\left( \frac{m}{2}x'(\tau)^2+V[x(\tau)]\right)}$

$\implies \boxed{ Z(\beta)=\mathcal{N} \int_{x(\hbar \beta)=x(0)}\mathcal{D}[x(\tau)] e^{-\frac{1}{\hbar} \int_{0}^{\beta \hbar}\text{d}\tau \mathcal{H}[p(\tau),x(\tau)]} }$. (I)

Notice that by comparing (H) and (I) we get an expression for the probabilistic weight $\Phi[x(\tau)]$ of a particular path $x(\tau)$, that is

$\Phi[x(\tau)] = \lim_{N\rightarrow \infty} \rho_u(x_0,x_1;\beta/N)\cdots \rho_u(x_{N-1},x_N;\beta/N) = \exp{\left\{ e^{-\frac{1}{\hbar} \int_{0}^{\beta \hbar}\text{d}\tau \mathcal{H}[p(\tau),x(\tau)]}\right\}}$ (J), which is very intuitive, considering the definition of the unnormalised density matrix $\rho_u$. This is an intriguing result, since we were able to obtain the complete statistical description of a quantum mechanical system without the appearance of complex numbers.

Because of this reason, using (J) it is easy to see why some paths contribute very little to the path integral: those are paths for which the exponent is very large due to high energy, and thus the integrand is negligibly small. Furthermore, it is unnecessary to consider whether or not nearby paths cancel each other's contributions, for in the present case they do not interfere (since no complex numbers involved) i.e. all contributions add together with some being large and others small.

#### Path integral Monte Carlo¶

In the algorithm, so called the naïve path integral (Markov-chain) Monte Carlo, we move from one path configuration consisting of $\{x_i\}$ to another one consisting of $\{x'_i\}$ by choosing a single position $x_k$ and by making a little displacement $\Delta x$ that can be positive or negative. We compute the weight before ($\Phi[\{x_i\}]$) this move and after ($\Phi[\{x'_i\}]$) the move and accept the move with the Metropolis acceptance rate (reject with certainty if the new weight is greater than the old one, smaller the new weight is, the higher the acceptance rate). Defining $\epsilon \equiv \beta/N$, we can approximate $\Phi[\{x_i\}]$ by making a Trotter decomposition only around the point $x_k$:

$\Phi\left[\{x_i\}\right]\approx \cdots \rho^\text{free}(x_{k-1},x_k;\epsilon) e^{-\frac{1}{2}\epsilon V(x_k)} e^{-\frac{1}{2}\epsilon V(x_k)} \rho^\text{free}(x_{k},x_{k+1};\epsilon)\cdots$.

Therefore, the acceptance ratio $\frac{\Phi\left[\{x'_i\}\right]}{\Phi\left[\{x_i\}\right]}$ can be approximated as

$\frac{\Phi\left[\{x'_i\}\right]}{\Phi\left[\{x_i\}\right]}\approx\frac{\rho^\text{free}(x_{k-1},x'_k;\epsilon) e^{-\epsilon V(x'_k)}\rho^\text{free}(x'_k,x_{k+1};\epsilon)}{\rho^\text{free}(x_{k-1},x_k;\epsilon) e^{-\epsilon V(x_k)} \rho^\text{free}(x_k,x_{k+1};\epsilon)}$.

This is implemented in the following program.

In [4]:
%pylab qt
import math, random, pylab, os

# Exact quantum position distribution:
def p_quant(x, beta):
p_q = sqrt(tanh(beta / 2.0) / pi) * exp(- x**2.0 * tanh(beta / 2.0))
return p_q

def rho_free(x, y, beta):        # free off-diagonal density matrix
return math.exp(-(x - y) ** 2 / (2.0 * beta))

output_dir = 'snapshots_naive_harmonic_path'
if not os.path.exists(output_dir): os.makedirs(output_dir)

fig = pylab.figure(figsize=(6, 10))
def show_path(x, k, x_old, Accepted, hist_data, step, fig):
pylab.clf()
path = x + [x[0]] #Final position is the same as the initial position.
#Note that this notation appends the first element of x as a new element to x
y_axis = range(len(x) + 1) #construct the imaginary time axis

#Plot the paths
if Accepted:
old_path = x[:] #save the updated path as the old path
old_path[k] = x_old #revert the update to get the actual old path
old_path = old_path + [old_path[0]] #final position is the initial position
ax.plot(old_path, y_axis, 'ko--', label='old path')
if not Accepted and step !=0:
old_path = x[:]
old_path[k] = x_old
old_path = old_path + [old_path[0]]
ax.plot(old_path, y_axis, 'ro-', label='rejection', linewidth=3)
ax.plot(path, y_axis, 'bo-', label='new path') #plot the new path
ax.legend()
ax.set_xlim(-2.5, 2.5)
ax.set_ylabel('$\\tau$', fontsize=14)
ax.set_title('Naive path integral Monte Carlo, step %i' % step)
ax.grid()

#Plot the histogram
x = [a / 10.0 for a in range(-100, 100)]
y = [p_quant(a, beta) for a in x]
ax.plot(x, y, c='gray', linewidth=1.0, label='Exact quantum distribution')
ax.hist(hist_data, 10, histtype='step', normed = 'True', label='Path integral Monte Carlo') #histogram of the sample
ax.set_title('Position distribution at $T=%.2f$' % T)
ax.set_xlim(-2.5, 2.5) #restrict the range over which the histogram is shown
ax.set_xlabel('$x$', fontsize = 14)
ax.set_ylabel('$\pi(x)=e^{-\\beta E_n}|\psi_n(x)|^2$', fontsize = 14)
ax.legend(fontsize = 6)
ax.grid()

pylab.pause(0.2)
pylab.savefig(output_dir + '/snapshot_%05i.png' % step)

beta = 4.0                                           # inverse temperature
T = 1 / beta
N = 8                                                # number of (imagimary time) slices
dtau = beta / N
delta = 1.0                                          # maximum displacement on one slice
n_steps = 4                                        # number of Monte Carlo steps
hist_data = []
x = [random.uniform(-1.0, 1.0) for k in range(N)]    # initial path (a position for each time)
show_path(x, 0, 0.0, False, hist_data, 0, fig)                       #show the initial path

for step in range(n_steps):
#print 'step',step
k = random.randint(0, N - 1)                     # randomly choose slice
knext, kprev = (k + 1) % N, (k - 1) % N          # next/previous slices
x_old = x[k]
x_new = x[k] + random.uniform(-delta, delta)     # new position at slice k
#calculate the weight before and after the move
old_weight  = (rho_free(x[knext], x_old, dtau) *
rho_free(x_old, x[kprev], dtau) *
math.exp(-0.5 * dtau * x_old ** 2))
new_weight  = (rho_free(x[knext], x_new, dtau) *
rho_free(x_new, x[kprev], dtau) *
math.exp(-0.5 * dtau * x_new ** 2))
if random.uniform(0.0, 1.0) < new_weight / old_weight: #accept with metropolis acceptance rate
x[k] = x_new
Accepted = True
else:
Accepted = False
show_path(x, k, x_old, Accepted, hist_data, step + 1, fig)
hist_data.append(x[k])

Populating the interactive namespace from numpy and matplotlib


Note that the above program satisfies the detailed balance condition, is aperiodic and irreducible. However it is very slow, as it takes very long to explore all of the available configuration space.

To make progress, let us now move back by two steps: consider the free particle instead of the harmonic potential, and instead of sampling $x_k$ with random $k$ at each step, keep $k$ fixed: we modify the position $x_k$ between the fixed positions $x_{k\pm 1}$ with the Metropolis algorithm.

The distribution of the position $x_k$ between fixed values of $x_{k\pm 1}$ is a gaussian. This is not a coincidence. Because $\pi^\text{free}(x_k|x_{k-1},x_{k+1})\propto \rho^\text{free}(x_{k-1},x_k;\epsilon) \rho^\text{free}(x_{k},x_{k+1};\epsilon)\propto \exp{\left[\frac{(x_k-\langle x_k \rangle)^2}{2\sigma^2}\right]}$, i.e. the product of two gaussians is also a gaussian. Here we have the mean $\langle x_k \rangle = \frac{x_{k-1} + x_{k+1}}{2}$, and the variance $\sigma=\frac{\epsilon}{2}$.

On the other hand, since we can derive that the distribution is a gaussian $\pi(x) = \sqrt{\tanh(\beta / 2) / \pi} \exp[-x^2 * \tanh(\beta/2)]$, we do not need to use the Metropolis algorithm. We can use direct sampling of the gaussian with the above mean and variance instead! However, instead of doing that, let us generalise the problem further.

Consider a position $x_k$ between $x'$ and $x''$ with corresponding intervals $\epsilon'$ and $\epsilon''$ (not necessarily equal), respectively. Then $\pi^\text{free}(x_k|x',x'')\propto \rho^\text{free}(x',x_k;\epsilon') \rho^\text{free}(x_{k},x'';\epsilon'')\propto \exp{\left[\frac{(x_k-\langle x_k \rangle)^2}{2\sigma^2}\right]}$ with the mean $\langle x_k \rangle = \frac{\epsilon''x' + \epsilon'x''}{\epsilon'+\epsilon''}$, and the variance $\sigma=\left(\frac{1}{\epsilon'} + \frac{1}{\epsilon''}\right)^{-1}$.

In [4]:
%pylab inline
import math, random, pylab

def rho_free(x, y, beta):
return math.exp(-(x - y) ** 2 / (2.0 * beta))

dtau_prime  = 0.1
dtau_dprime = 0.2
x_prime  = 0.0
x_dprime = 1.0
delta = 1.0                 # maximum displacement of xk
n_steps = 100000            # number of Monte Carlo steps
data_hist = []
xk = 0.0                    # initial value of xk
for step in xrange(n_steps):
xk_new = xk + random.uniform(-delta, delta)
old_weight  = (rho_free(x_dprime, xk, dtau_dprime) *
rho_free(xk, x_prime, dtau_prime))
new_weight  = (rho_free(x_dprime, xk_new, dtau_dprime) *
rho_free(xk_new, x_prime, dtau_prime))
if random.random() < new_weight / old_weight:
xk = xk_new
data_hist.append(xk)

def pi_analytic(xk, x_prime, x_dprime, dtau_prime, dtau_dprime):
mean = (dtau_dprime * x_prime + dtau_prime * x_dprime) / (dtau_prime + dtau_dprime)
sigma = 1.0 / math.sqrt(1.0 / dtau_prime + 1.0 / dtau_dprime)
return math.exp(-(xk - mean) ** 2 / (2.0 * sigma ** 2)) / math.sqrt(2.0 * math.pi) / sigma

pylab.title('Distribution on slice k', fontsize=18)
histo, bin_edges, dummy = pylab.hist(data_hist, bins=100, normed=True)
bin_centers = 0.5 * (bin_edges[1:] + bin_edges[:-1])
pylab.plot(bin_centers, [pi_analytic(x, x_prime, x_dprime, dtau_prime, dtau_dprime) for x in bin_centers], 'r-', lw=3)
pylab.xlabel('$x_k$', fontsize=18)
pylab.ylabel('$\pi(x_k)$', fontsize=18)
pylab.savefig('plot-path_slice.png')

Populating the interactive namespace from numpy and matplotlib


#### Lévy sampling of quantum paths¶

$\pi(x_k | x', x'')$ gives the statistical weight of all paths that start at x' that pass through $x_k$ and that end up at $x''$. Let's apply this idea for $x_k = x_1 x' = x_0$ and $x'' = x_N$. We can sample the position $x_1$, given $x_0$ and $x_N$. Between the freshly sampled value of $x_1$ and $x_N$ we can then sample the value of $x_2$ and then of $x_3\cdots$ and eventually the entire path. Here we use $\pi^\text{free}(x_k|x',x'')\propto \exp{\left[\frac{(x_k-\langle x_k \rangle)^2}{2\sigma^2}\right]}$ with the mean $\langle x_k \rangle = \frac{\epsilon''x' + \epsilon'x''}{\epsilon'+\epsilon''}$, and the variance $\sigma=\left(\frac{1}{\epsilon'} + \frac{1}{\epsilon''}\right)^{-1}$ in order to apply this condition.

In [7]:
import math, random, pylab, os

fig = pylab.figure(figsize=(6, 10))
def show_path(x,step,title):
path = x[:] #Final position is the same as the initial position.
#Note that this notation appends the first element of x as a new element to x
y_axis = range(len(x)) #construct the imaginary time axis
#Plot the paths
pylab.plot(path, y_axis, 'ko-')
pylab.xlabel('$x$', fontsize = 14)
pylab.ylabel('$\\tau$', fontsize=14)
pylab.title(title, fontsize=14)
pylab.xlim([-1,2])
pylab.grid()
pylab.savefig(output_dir + '/snapshot_%05i.png' % step)
pylab.clf()

output_dir = 'snapshots_Levy_harmonic_path'
if not os.path.exists(output_dir): os.makedirs(output_dir)

beta = 1.0
N = 4
dtau = beta / N
nsteps = 100                     # number of paths to be generated
xstart, xend = 0.0, 1.0          # initial and final points
for step in range(nsteps):
title = 'Levy sampling quantum path, step %i' % step
x = [xstart] #initialise the path
for k in range(1, N):        # loop over internal slices
dtau_prime = (N - k) * dtau #current "time" step
x_mean = (dtau_prime * x[k - 1] + dtau * xend) / \
(dtau + dtau_prime) #current mean is calculated from the above formula
sigma = math.sqrt(1.0 / (1.0 / dtau + 1.0 / dtau_prime)) #current variance calculated similarly
x.append(random.gauss(x_mean, sigma)) #generate a gaussian random variable using above parameters
x.append(xend)
show_path(x,step,title)

<Figure size 432x720 with 0 Axes>

Note that without the condition that the path will finish at $x_N$, this is simply a random walk.

In [98]:
import math, random

beta = 4.0
N = 8
sigma = math.sqrt(beta / N)
x = [0.0]
title = 'Random walk'
for k in range(N - 1):
x.append(random.gauss(x[-1], sigma))
fig = pylab.figure(figsize=(6, 10))
show_path(x,title)

In [99]:
import math, random

beta = 1.0
N = 8
sigma = math.sqrt(beta / N)
xend = 1.0
Upsilon = [0.0]
title = 'Trivial quantum path algorithm'
for k in range(N):
Upsilon.append(random.gauss(Upsilon[-1], sigma))
x = [0.0] + [Upsilon[k] + (xend - Upsilon[-1]) * \
k / float(N) for k in range(1, N + 1)]
fig = pylab.figure(figsize=(6, 10))
show_path(x,title)


Since the harmonic oscillator density matrix is also a gaussian (see e.g. Feynman), we can use the Lévy sampling here.

In [4]:
%pylab inline
import math, random, pylab

beta = 2.0
N = 20
dtau = beta / N
nsteps = 1
xstart, xend = 2.0, 1.0
fig = pylab.figure(figsize=(6, 10))
for step in range(nsteps):
pylab.clf()
x = [xstart]
for k in range(1, N):
dtau_prime = (N - k) * dtau
Ups1 = 1.0 / math.tanh(dtau) + \
1.0 / math.tanh(dtau_prime)
Ups2 = x[k - 1] / math.sinh(dtau) + \
xend / math.sinh(dtau_prime)
x.append(random.gauss(Ups2 / Ups1, \
1.0 / math.sqrt(Ups1)))
x.append(xend)
# graphics
pylab.plot(x, [j * dtau for j in range(N + 1)], 'k-o')
pylab.xlabel('$x$', fontsize=18)
pylab.ylabel('$\\tau$', fontsize=18)
pylab.title('Harmonic Levy path %i' % step)
#pylab.xlim(-2.0, 4.0)
pylab.grid()
pylab.show()
pylab.pause(0.1)

Populating the interactive namespace from numpy and matplotlib


We can easily generalise this to the multidimensional case. We have no problem understanding that a three-dimensional harmonic oscillator quantum path can be sampled by independent Lévy constructions in x, y and z. The following is a 2D implementation.

In [91]:
%matplotlib qt
import math, random, os
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

def levy_harmonic_1d(start, end, dtau):
x = [start]
for k in range(1, N):
dtau_prime = (N - k) * dtau
Ups1 = 1.0 / math.tanh(dtau) + \
1.0 / math.tanh(dtau_prime)
Ups2 = x[k - 1] / math.sinh(dtau) + \
end / math.sinh(dtau_prime)
x.append(random.gauss(Ups2 / Ups1, \
1.0 / math.sqrt(Ups1)))
x.append(end)
return x

def show_path_2d(x,y,title):
plt.gca(projection='3d')
plt.plot(x, y, [j * dtau for j in range(N + 1)], 'k-o')
plt.title(title)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.gca(projection='3d').set_zlabel('$\\tau$')
plt.grid()
plt.show()
plt.pause(0.1)
plt.clf()

output_dir = 'snapshots_Levy_2d_harmonic_path'
if not os.path.exists(output_dir): os.makedirs(output_dir)

beta = 1.0
N = 20
nsteps = 50
dtau = beta / float(N)
[xstart, ystart, zstart] = [1.0, -2.0, 1.5]
[xend, yend, zend] = [-2.5, 0.0, -0.5]

for step in range(nsteps):
plt.clf()
title = '2D quantum path %.i th step' % step
x = levy_harmonic_1d(xstart, xend, dtau)
y = levy_harmonic_1d(ystart, yend, dtau)
#z = levy_harmonic_1d(zstart, zend, dtau)
#for i in range(N + 1):
#    print 'slice %2i:  ' % i, x[i], y[i], z[i]
show_path_2d(x,y,title)
plt.savefig(output_dir + '/snapshot_%05i.png' % step)

<Figure size 864x576 with 0 Axes>

### Path sampling methods implemented¶

#### Preliminaries¶

P.1.a The following program proposes uniformly distributed configurations in the heliport square (-1 to 1 in $x$, and -1 to 1 in $y$), then accepts these configurations with the mixed harmonic and quartic weight $\pi(x,y)=\exp[-0.5 (x^2 + y^2) - \alpha (x^4 + y^4)]$, for $\alpha = 0.5$. Also produces a two-dimensional histogram to visualise the data. Note that this is a pure rejection-sampling algorithm.

In [56]:
import random, math, pylab

alpha = 0.5
nsamples = 1000000
samples_x = []
samples_y = []
for sample in xrange(nsamples):
while True:
x = random.uniform(-1.0, 1.0)
y = random.uniform(-1.0, 1.0)
p = math.exp(-0.5 * (x ** 2 + y ** 2) - alpha * (x ** 4 + y ** 4))
if random.uniform(0.0, 1.0) < p:
break
samples_x.append(x)
samples_y.append(y)

pylab.hexbin(samples_x, samples_y, gridsize=50, bins=1000)
pylab.axis([-1.0, 1.0, -1.0, 1.0])
cb = pylab.colorbar()
pylab.xlabel('x')
pylab.ylabel('y')
pylab.title('A1_1')
pylab.savefig('plot_A1_1.png')
pylab.show()


P.1.b The same algorithm can be implemented by proposing the random variables with gaussian distribution with a cutoff to restrict into the region $[-1,1]$. Then the algorithm becomes partly rejection-free for the harmonic part of $\pi(x,y)$. Note that the two programs do not agree for negative values of alpha.

In [59]:
import random, math, pylab

def gauss_cut():
while True:
x = random.gauss(0.0, 1.0)
if abs(x) <= 1.0:
return x

alpha = 0.5
nsamples = 1000000
samples_x = []
samples_y = []
for sample in xrange(nsamples):
while True:
x = gauss_cut()
y = gauss_cut()
p = math.exp(- alpha * (x ** 4 + y ** 4))
if random.uniform(0.0, 1.0) < p:
break
samples_x.append(x)
samples_y.append(y)

pylab.hexbin(samples_x, samples_y, gridsize=50, bins=1000)
pylab.axis([-1.0, 1.0, -1.0, 1.0])
cb = pylab.colorbar()
pylab.xlabel('x')
pylab.ylabel('y')
pylab.title('A1_2')
pylab.savefig('plot_A1_2.png')
pylab.show()


P.2.a We can break $\pi(x,y)$ into a product of one dimensional distributions, propose uniform random variables between -1, 1. Note that although a direct-sampling approach is used in the coordinate that is changed, this is now a Markov-chain sampling algorithm because the outcome of the move depends on the present configuration.

In [45]:
import random, math, pylab

alpha = 0.5
nsteps = 1000000
samples_x = []
samples_y = []
x, y = 0.0, 0.0
for step in range(nsteps):
if step % 2 == 0:
while True:
x = random.uniform(-1.0, 1.0)
p = math.exp(-0.5 * x ** 2 - alpha * x ** 4)
if random.uniform(0.0, 1.0) < p:
break
else:
while True:
y = random.uniform(-1.0, 1.0)
p = math.exp(-0.5 * y ** 2 - alpha * y ** 4)
if random.uniform(0.0, 1.0) < p:
break
samples_x.append(x)
samples_y.append(y)

pylab.hexbin(samples_x, samples_y, gridsize=50, bins=1000)
pylab.axis([-1.0, 1.0, -1.0, 1.0])
cb = pylab.colorbar()
pylab.xlabel('x')
pylab.ylabel('y')
pylab.title('A2_1')
pylab.savefig('plot_A2_1.png')
pylab.show()


P.2.a Propose gaussian variables as in P.1.b.

In [46]:
import random, math, pylab

def gauss_cut():
while True:
x = random.gauss(0.0, 1.0)
if abs(x) <= 1.0:
return x

alpha = 0.5
nsteps = 1000000
samples_x = []
samples_y = []
x, y = 0.0, 0.0
for step in range(nsteps):
if step % 2 == 0:
while True:
x = gauss_cut()
p = math.exp(-alpha * x ** 4)
if random.uniform(0.0, 1.0) < p:
break
else:
while True:
y = gauss_cut()
p = math.exp(-alpha * y ** 4)
if random.uniform(0.0, 1.0) < p:
break
samples_x.append(x)
samples_y.append(y)

pylab.hexbin(samples_x, samples_y, gridsize=50, bins=1000)
pylab.axis([-1.0, 1.0, -1.0, 1.0])
cb = pylab.colorbar()
pylab.xlabel('x')
pylab.ylabel('y')
pylab.title('A2_2')
pylab.savefig('plot_A2_2.png')
pylab.show()


P.3.a Full-fledged Markov-chain Monte Carlo sampling using Metropolis algorithm:

In [74]:
import random, math, pylab

alpha = 0.5
nsteps = 1000000
samples_x = []
samples_y = []
x, y = 0.0, 0.0
for step in range(nsteps):
xnew = random.uniform(-1.0, 1.0)
ynew = random.uniform(-1.0, 1.0)
exp_new = - 0.5 * (xnew ** 2 + ynew ** 2) - alpha * (xnew ** 4 + ynew ** 4)
exp_old = - 0.5 * (x ** 2 + y ** 2) - alpha * (x ** 4 + y ** 4)
#Use the Metropolis algorithm:
if random.uniform(0.0, 1.0) < math.exp(exp_new - exp_old): #the transition rate pi(x_old->x_new)
x = xnew
y = ynew
samples_x.append(x)
samples_y.append(y)

pylab.hexbin(samples_x, samples_y, gridsize=50, bins=1000)
pylab.axis([-1.0, 1.0, -1.0, 1.0])
cb = pylab.colorbar()
pylab.xlabel('x')
pylab.ylabel('y')
pylab.title('A3_1')
pylab.savefig('plot_A3_1.png')
pylab.show()


P.3.b Again, we can improve the above program by proposing gaussian random variables instead of uniform random variables, hence the Metropolis acceptance rate simplifies.

In [77]:
import random, math, pylab

def gauss_cut():
while True:
x = random.gauss(0.0, 1.0)
if abs(x) <= 1.0:
return x

alpha = 0.5
nsteps = 1000000
samples_x = []
samples_y = []
x, y = 0.0, 0.0
for step in range(nsteps):
xnew = gauss_cut()
ynew = gauss_cut()
exp_new = - alpha * (xnew ** 4 + ynew ** 4)
exp_old = - alpha * (x ** 4 + y ** 4)
#Use the Metropolis algorithm:
if random.uniform(0.0, 1.0) < math.exp(exp_new - exp_old): #the transition rate pi(x_old->x_new)
x = xnew
y = ynew
samples_x.append(x)
samples_y.append(y)

pylab.hexbin(samples_x, samples_y, gridsize=50, bins=1000)
pylab.axis([-1.0, 1.0, -1.0, 1.0])
cb = pylab.colorbar()
pylab.xlabel('x')
pylab.ylabel('y')
pylab.title('A3_2')
pylab.savefig('plot_A3_2.png')
pylab.show()


#### Harmonic oscillator¶

##### Naive sampling¶

The naive algorithm uses the Trotter decomposition formula, hence it is only correct in the limit $\beta/N \rightarrow 0$. Although it is slow, it provides a good approximation and constitutes an intuitive starting point.

In [66]:
import math, random, pylab, os

def show_path(x, step, title):
fig = pylab.figure(figsize=(6, 10))
pylab.clf()
path = x[:] #Final position is the same as the initial position.
#Note that this notation appends the first element of x as a new element to x
y_axis = range(len(x)) #construct the imaginary time axis
#Plot the paths
pylab.plot(path, y_axis, 'ko-')
pylab.xlabel('$x$', fontsize = 14)
pylab.ylabel('$\\tau$', fontsize=14)
pylab.title(title, fontsize=14)
pylab.grid()
pylab.savefig('snapshot_%i.png' % step)

def rho_free(x, y, beta):
return math.exp(-(x - y) ** 2 / (2.0 * beta))

beta = 20.0
N = 80
dtau = beta / N
delta = 1.0
n_steps = 4000000
x = [5.0] * N
#show_path(x, 1, 'initial configuration')
data = []
for step in range(n_steps):
k = random.randint(0, N - 1)
knext, kprev = (k + 1) % N, (k - 1) % N
x_new = x[k] + random.uniform(-delta, delta)
old_weight  = (rho_free(x[knext], x[k], dtau) *
rho_free(x[k], x[kprev], dtau) *
math.exp(-0.5 * dtau * x[k] ** 2))
new_weight  = (rho_free(x[knext], x_new, dtau) *
rho_free(x_new, x[kprev], dtau) *
math.exp(-0.5 * dtau * x_new ** 2))
if random.uniform(0.0, 1.0) < new_weight / old_weight:
x[k] = x_new
if step % N == 0:
k = random.randint(0, N - 1)
data.append(x[k])

#Write the final path configuration on a file:
filename = 'naive_harmonic_path_configuration_N%i.txt' % N
f = open(filename, 'w')
for i in range(N):
f.write(str(x[i])+ '\n')
f.close()

show_path(x, step, 'Naive harmonic path final configuration, N = %i' %N)
pylab.show()

pylab.clf()
fig = pylab.figure(figsize=(7, 5))
pylab.hist(data, normed=True, bins=100, label='QMC')
list_x = [0.1 * a for a in range (-30, 31)]
list_y = [math.sqrt(math.tanh(beta / 2.0)) / math.sqrt(math.pi) * \
math.exp(-x ** 2 * math.tanh(beta / 2.0)) for x in list_x]
pylab.plot(list_x, list_y, label='analytic')
pylab.legend()
pylab.xlabel('$x$')
pylab.ylabel('$\\pi(x)$ (normalized)')
pylab.title('naive_harmonic_path (beta=%s, N=%i)' % (beta, N))
pylab.xlim(-2, 2)
pylab.savefig('plot_B1_beta%s.png' % beta)
pylab.show()

<Figure size 432x288 with 0 Axes>
##### Lévy sampling¶

The Lévy algorithm for the harmonic oscillator paths is exact, it uses no approximation.

In [67]:
import math, random, pylab, os

def show_path(x, step, title):
fig = pylab.figure(figsize=(6, 10))
pylab.clf()
path = x[:] #Final position is the same as the initial position.
#Note that this notation appends the first element of x as a new element to x
y_axis = range(len(x)) #construct the imaginary time axis
#Plot the paths
pylab.plot(path, y_axis, 'ko-')
pylab.xlabel('$x$', fontsize = 14)
pylab.ylabel('$\\tau$', fontsize=14)
pylab.title(title, fontsize=14)
pylab.grid()
pylab.savefig('levy_snapshot_%i.png' % step)

def levy_harmonic_path(xstart, xend, dtau, N):
x = [xstart]
for k in range(1, N):
dtau_prime = (N - k) * dtau
Ups1 = 1.0 / math.tanh(dtau) + \
1.0 / math.tanh(dtau_prime)
Ups2 = x[k - 1] / math.sinh(dtau) + \
xend / math.sinh(dtau_prime)
x.append(random.gauss(Ups2 / Ups1, \
1.0 / math.sqrt(Ups1)))
return x

beta = 20.0
N = 80
dtau = beta / N
delta = 1.0
n_steps = 40000
x = [5.0] * N
data = []
for step in range(n_steps):
Ncut = N/2
x = x[Ncut:] + x[:Ncut] #slices the x data from the middle and swaps the two parts, hence x[0] changes
x = levy_harmonic_path(x[0], x[0], dtau, N)
k = random.randint(0, N - 1)
data.append(x[k])

#Write the final path configuration on a file:
filename = 'levy_harmonic_path_configuration_levy_N%i.txt' % N
f = open(filename, 'w')
for i in range(N):
f.write(str(x[i])+ '\n')
f.close()

show_path(x, step, 'Levy harmonic final configuration, N = %i' %N)
pylab.show()

pylab.clf()
fig = pylab.figure(figsize=(7, 5))
pylab.hist(data, normed=True, bins=100, label='QMC')
list_x = [0.1 * a for a in range (-30, 31)]
list_y = [math.sqrt(math.tanh(beta / 2.0)) / math.sqrt(math.pi) * \
math.exp(-x ** 2 * math.tanh(beta / 2.0)) for x in list_x]
pylab.plot(list_x, list_y, label='analytic')
pylab.legend()
pylab.xlabel('$x$')
pylab.ylabel('$\\pi(x)$ (normalized)')
pylab.title('levy_harmonic_path (beta=%s, N=%i)' % (beta, N))
pylab.xlim(-2, 2)
pylab.savefig('plot_B2_beta%s.png' % beta)
pylab.show()

<Figure size 432x288 with 0 Axes>