High dimensional integration by Metropolis

The classical method for high-dimensional integration is the Vegas algorithm, discussed before. However, Vegas can not resolve poles, which are not along one of the axis.

Problem in physics very commonly require integrals over space of many coordinates, and might have the form: \begin{eqnarray} I(\vec{r}_0) = \int f(\vec{r}_0,\vec{r}_1,\vec{r}_2,\cdots,\vec{r}_N) d^3\vec{r}_1 d^3\vec{r}_2\cdots d^3\vec{r}_N \end{eqnarray} Here $f$ might have poles for many combinations of $\vec{r}_i-\vec{r}_j$ or $\vec{r}_i-\vec{r}_j-\vec{r}_k$, etc.

For a concrete example, we will think of the following convolutions \begin{eqnarray} I(\vec{r}_0) = {\cal N} \int e^{-\frac{(\vec{r}_1-\vec{r}_0)^2}{w^2}} e^{-\frac{(\vec{r}_2-\vec{r}_1)^2}{w^2}} \cdots e^{-\frac{(\vec{r}_{N-1}-\vec{r}_{N})^2}{w^2}} e^{-\frac{{\vec{r}_N}^2}{w^2}} d^3\vec{r}_1 d^3\vec{r}_2\cdots d^3\vec{r}_N \end{eqnarray} with noralization ${\cal N}= \frac{(N+1)^{3/2}}{(w \sqrt{\pi})^{3 N}},$ the value of the integral is $I(\vec{r}_0)=e^{-\frac{\vec{r}_0^2}{(N+1) w^2}}$

If width $w$ is small, this is very difficult problem for MC, even for moderate values of $N$.


We could create a Markov's chain and take as the probability for a step proportionaly to $|f(x)|$, i.e, the $T(X\rightarrow X') = min(|f(X')|/|f(X)|,1)$. Here $X$ stands for configuration $(\vec{r}_0,\vec{r}_1,\cdots,\vec{r}_N)$.

We would meassure $f(X)/|f(X)|=\textrm{sign}(f(X))$. As a result we would get correct shape of the function $I(\vec{r}_0)$, however, the integral itself would be unknown.

Jumping between spaces

The most common way to determine the unknown constant is to define two spaces, the Physical space and the Alternative space. The function in the alternative space should be as similar as possible to the function we are integrating, but its integral must be known, and we should be allowed to jump between the two spaces. Suppose the function in the alternative space is $V_0 f_m(X)$, where $f_m(X)$ is normalized to unity ($\int f_m(X) dX=1$) and positive function, and $V_0$ is an arbitrary constant.

Probability density in the physical space $P(X)$ is then taken to be proportional to $|f(x)|$ in the physical Hilbert space, and $V_0 f_m(x)$ in the alternative space. The probability to visit a configuration $X$ in the physical and alternative space would be: \begin{eqnarray} & P_{physical}(X) = & |f(X)| \; {\cal C}\\ & P_{alternative}(X) = & V_0 f_m(X) \; {\cal C} \end{eqnarray}

In Metropolis we would have two types of moves:

If we are in the Physical Hilbert space, we would accept the step with $T(X\rightarrow X') = min(|f(X')|/|f(X)|,1)$. If we are in alternative space, we would accept the step with $T(X\rightarrow X') = min(f_m(X')/f_m(X),1)$, and if we are jumping from Physical to the alternative space, we would have $T(X_p\rightarrow X'_a) = min(V_0 f_m(X'_a)/|f(X)|,1)$.

More precisely, if there is non-symmetric trial step probability, the Metropolis acceptance probability would be $A(X\rightarrow X')=min\left(1,\frac{|f(X')|\omega_{X'\rightarrow X}}{|f(X)|\omega_{X\rightarrow X'}}\right)$, where $\omega_{X'\rightarrow X}$ is the trial step probability. But let's first think of steps which have symmetric trial step probability.

We would than sample the following quantities \begin{eqnarray} V_{physical}= \sum_{i=physical} \frac{f(X_i)}{|f(X_i)|} \end{eqnarray} when we are in the physical space, and when we are in alternative space we would just count steps in this space \begin{eqnarray} V_{alternative}= \sum_{i=alternative} 1 \end{eqnarray}

These quantities will converge to the following values \begin{eqnarray} V_{physical} &=& {\cal C} \int dX\; |f(X)| \frac{f(X)}{|f(X)|}\\ V_{alternative} &=& {\cal C} \int dX\; V_0 |f_m(X)| 1 \end{eqnarray} where $C$ is an arbitrary number, which we do not know. This is because the probability to reach a configuration $X$ in the Physical space is $\propto dX |f(X)|$ and in alternative space is $\propto dX V_0 |f_m(X)|$.

Because the function in alternative space is normalized $\int dX |f_m(x)|=1$, we conclude that \begin{eqnarray} \frac{V_{physical}}{V_{alternative}}=\frac{1}{V_0} \int dX f(X) \end{eqnarray} hence we know the integral \begin{eqnarray} \int dX f(X) = V_0 \frac{V_{physical}}{V_{alternative}}=V_0 \frac{\sum_{i\in physical} \textrm{sign}(f(X_i))}{\sum_{i\in alternative} 1} \end{eqnarray}

The art in this approach is to find a good function $f_m(X)$, which we know how to normalize, and adjust $V_0$ so that we spent some finite (but not too much time) in alternative space. A good rule is 90% in physical space, and 10% in alternative space.

The best approach is to self-consistently determine both $f_m(X)$ and the constant $V_0$ during the sampling.

But before we discuss how to determine both, we would sketch a simpler, and usually even more efficient, algorithm.

Evaluation of both functions on each configuration

Alternatively we can define the probability for visiting configuration $X$ to be \begin{eqnarray} P \propto |f(X)|+ V_0 f_m(X) \end{eqnarray} and we could have only the physical space of configurations $X$, on which we would simultaneously evaluate both $f(X)$ and $f_m(X)$.

We have freedom to adjust $V_0$. The idea is to adjust $V_0$ so that on average $|f(X)|$ is around 10-times larger than $V_0 f_m(X)$. As a result we would still visit most often those configurations in which function $|f(X)|$ is large, and if $f_m(X)$ is similar to $f(X)$ we would not visit configurations where $|f(X)|$ is very small. Even if the overlap between $f(X)$ and $f_m(X)$ is small, we will still visit configurations where $f(X)$ is large more often, because we will make sure that $V_0$ is such that on avergae $f(X)$ contributes more to the weight that $V_0 f_m(X)$.

The transition probability is therefore \begin{eqnarray} T(X\rightarrow X') = \frac{|f(X')| + V_0 f_m(X')}{|f(X)| + V_0 f_m(X)} \end{eqnarray}

We will sample two quantities \begin{eqnarray} V_{physical} = \sum_i \frac{f(X_i)}{|f(X_i)| + V_0 f_m(X_i)}\\ V_{alternative} = \sum_i \frac{V_0 f_m(X_i)}{|f(X_i)| + V_0 f_m(X_i)} \end{eqnarray}

The two quantities will converge towards \begin{eqnarray} V_{physical} &=& {\cal C}\int dX\; (|f(X)| + V_0 f_m(X)) \frac{f(X)}{|f(X)| + V_0 f_m(X)} = {\cal C} \int dX f(X) \\ V_{alternative} &=& {\cal C} \int dX\; (|f(X)| + V_0 f_m(X)) \frac{V_0 f_m(X)}{|f(X)| + V_0 f_m(X)}={\cal C}\; V_0 \end{eqnarray}

The desired integral will again be \begin{eqnarray} \int dX f(X) = V_0 \frac{V_{physical}}{V_{alternative}} \end{eqnarray}

During the simulations of the Markov chain we will adjust $V_0$ so that \begin{eqnarray} \frac{\widetilde{V}_{physical}}{V_{alternative}} \approx 10, \end{eqnarray} where \begin{eqnarray} \widetilde{V}_{physical}=\sum_i \frac{|f(X_i)|}{|f(X_i)| + V_0 f_m(X_i)} \end{eqnarray}

Suppose that during simulation we realize that currently $\frac{\widetilde{V}_{physical}}{V_{alternative}} \approx 1$. It means that $V_{alternative}$ is around 10-times too large, hence we would reduce $V_0$ for a factor (maybe 2 initially, and than recheck). We would also need to reduce with the same constant the current sum $V_{alternative}$, because it is proportional to $V_0$. We will keep adjusting $V_0$ untill the ratio of $\frac{\widetilde{V}_{physical}}{V_{alternative}}$ is close to the desired number.

Measuring function

We are left to determine the optimal function $f_m(X)$, which is as similar as possible to $|f(X)|$, but with added constrained that is normalizable, i.e., $\int dX f_m(X)=1$.

An obvious choice (think of Vegas) is a separable ansatz: \begin{eqnarray} f_m(X=(\vec{r}_1,\vec{r}_r,\cdots,\vec{r}_N))=g_1(r_1) g_2(r_2)\cdots g_N(r_N) \end{eqnarray}

To determine the projections $g_i(r_i)$, we would self-consistently project the sampled function $|f(X)|$ to all axis, just like in Vegas algorithm: \begin{eqnarray} g_i(x)\propto \int d\vec{r_1} d\vec{r}_2 \cdots d\vec{r}_N |f(X)|\delta(r_i-x) \end{eqnarray}

The integration of separable function is simple \begin{eqnarray} \int dX f_m(X) = \prod_{i} \int d\vec{r}_i g_i(r_i) \end{eqnarray} To avoid systematic error, we will not treat $g_i(r_i)$ function as continuous function, but rather as a step-wise constant function, with the following property \begin{eqnarray} g(x) = \sum_{l=0}^{M-1} g[l]\theta(x_l<x<x_{l+1}). \end{eqnarray} Unfortunately this makes analysis tedious, and integrals a bit more involved.

Newertheless, for separable function, we can readily compute \begin{eqnarray} &&1D : \int d\vec{r}_i g_i(r_i)= \Delta \sum_{l=0}^{M-1} g_i[l]\\ &&3D : \int d\vec{r}_i g_i(r_i)= \frac{4\pi \Delta^3}{3} \sum_{l=0}^{M-1} ((l+1)^3-l^3)g_i[l] \end{eqnarray}

For more complicated ansatz of convolutions, the analytic formulas can still be derived, and the equations are given below.


Here we will not show all details of implementing the measuring function. It is implemented in module mweight.py, which we will include here. We will rather implement Metropolis Markov's chain and wew will use 3D-convolutions as the measuring function.

To test the implementation, we will use Gaussian functions, with peaks along the diagonals (not along the axes). The Gausian's have the form \begin{eqnarray} f(k_0,k_1,\cdots k_{N-1}) = C e^{k_{N-1}^2/w^2} e^{-|k_{N-2}-k_{N-1}|^2/w^2} e^{-|k_{N-3}-k_{N-2}|^2/w^2}\cdots e^{-|k_1-k_0|^2/w^2} \end{eqnarray} where \begin{eqnarray} C = \frac{N^{3/2} }{ ( \sqrt{\pi} w )^{3(N-1)}} \end{eqnarray}

We will define several parameters for MC integration

This is the core of the algorithm, which implements the Metropolis random walk, and measures both functions $f$ and $f_m$ on each configuration, using $|f|+V_0 f_m$ as the weight of each configuration.

At the first iteration, the weight function $f_m$ will be

$$f_m(r) \propto \theta(r<k_F) + \left(\frac{k_F}{r}\right)^{dexp} \theta(r>k_F) $$

Later it is adjusted to

\begin{eqnarray} f_m(\vec{r}_1,\vec{r}_2,\cdots,\vec{r}_n)= g_1(r_1) h_1(|\vec{r}_n-\vec{r}_1|) g_2(r_2)h_2(|\vec{r}_n-\vec{r}_2|)\cdots g_n(r_n) \end{eqnarray}

where $g_i$ and $h_i$ are determined self-consistently from the histograms.

The external variable $r_0$ will be computed on the user defined mesh, which is named qx[i].

Since $r_0$ is only the length of the vector in 3D, we will sample (and integrate) over its angles $\theta$ and $\phi$. Here $$\vec{r}_0=\left(qx_i \sin(\theta)\cos(\phi), qx_i \sin(\theta)\sin(\phi), qx_i \cos(\theta)\right).$$

However, when we choose variables $(r,\theta,\phi)$ instead of $(x,y,z)$ the trial step probability is not symmetric, because we want to integrate uniformly over the entire volume $d(\cos(\theta)) d\phi\, d^3r_1 d^3 r_2\cdots d^3 r_N $.

The probability to jump into any configuration is proportional to its contribution to the volume, i.e, $\omega_{X\rightarrow X'} = \frac{d\Omega_{X'}}{d\Omega_{X}}= \frac{sin(\theta_{X'}) d\theta d\phi}{sin(\theta_{X}) d\theta d\phi} = \frac{sin(\theta_{X'})}{sin(\theta_{X})}$. If we were sampling the integral over $r$ as well, the trial step probability would be $\omega_{X\rightarrow X'} = \frac{dV_{X'}}{dV_{X}}= \frac{r_{X'}^2 sin(\theta_{X'})}{r_X^2 sin(\theta_{X})}$. In this example, $r_0$ is the external variable, which we do not integrate over, hence we should not add $r^2$ to trial step probability. The trial step will determine new $\theta_{new} = \xi_0 \pi$, and $\phi_{new}= \xi_1 2\pi$ with the trial step probability $\frac{\sin(\theta_{new})}{\sin(\theta_{old})}$

The internal variables of integration $r_1,r_2\cdots r_N$ (are called momentum) are samples in cartesian coordinates, hence the trial step probability is symmetric. The trial step will be $$ \vec{r}_{new} \rightarrow \vec{r}_{old} + dk\;(2\xi_1-1, 2\xi_2-1, 2\xi_3-1),$$ where $dk$=p.dkF. But this is only provided that $|r_{new}|<= \Lambda_{cutoff}$. If $|r_{new}|>\Lambda_{cutoff}$, we instantly reject the step. The sampling is thus constrained to a set of spheres, so that none of particles distances (or momenta) can be too large.

As a test, we will use $6\times 3=18$-dimensional space, and Gaussians of width $1/\sqrt{6}$. The result should be Gaussian of width 1.

We expect the result to be gausian function with width of unity, i.e.,

Speed up by numba

We will speed up simulation using numba. All parts that are called every single step are being recast into numba functions.

The function, which is being evaluated, will also be speed-up by Numba

Checking that these Numba functions work.

Finally, using these Numba function in the Metropolis routine. We will also time steps, to understand if the code is optimized.

Next plotting the auxiliary function $f_m$ projections $g(r)$ to all $r$.

Here is the auxiliary function $f_m$ for convolutios $h(r)$ versus $r$.

1D Convolutions

The power of Metropolis integration is that we have freedom to pick an appropriate normalized function. This can come from physical intuition into processes being modeled. In most physical process the convolutions (correlations between particles) turn out to be important, hence convolution type ansatz considerably improves the precision. The convolutions can still be worked out analytically, even though they are more tedious to derive.

To derive the integral, we will concentrate on a specific form of convolution,i.e., \begin{eqnarray} f_m(x_1,x_2,\cdots,x_n)= \sum_{j=1}^n g_j(x_j)\prod_{i=1,i\ne j}^n g_i(x_i) h_{ij}(x_j-x_i) \end{eqnarray} As before, the two functions can be determined self-consistently by projections \begin{eqnarray} &g_i(x)=&\int dy_1 dy_2\cdots dy_n |f(y_1,y_2,\cdots,y_n)|\delta(y_i-x)\\ &h_{ij}(x_j-x_i)=&\int dy_1 dy_2\cdots dy_n |f(y_1,y_2,\cdots,y_n)|\delta(x_j-x_i-(y_j-y_i)) \end{eqnarray}

To normalize $f_m$ we need to evaluate the integral \begin{eqnarray} \int dx_1\cdots dx_n f_m = \sum_{j=1}^n \int dx_j g(x_j) \prod_{i\ne j}^{n} \int dx_i g_i(x_i) h_{ij}(x_j-x_i) \end{eqnarray} Withouth loss of generality, we will work out one term in the sum, namely, \begin{eqnarray} I^N = \int dx_N g(x_N) \prod_{i\ne N} \int dx_i g_i(x_i) h_{ij}(x_j-x_i) \end{eqnarray} so that the end result is \begin{eqnarray} \int dx_1\cdots dx_n f_m = \sum_j I^j \end{eqnarray}

This integral is combersome because all the functions are pice-wise constant, which leads to somewhat lengthly derivations. We can write $x_N= \Delta (j_N + t)$, $x_i = \Delta(j_i + u)$ with $t\in[0,1]$ and $u\in[0,1]$ and $j_i$ integers. Note that $g(\Delta(j+u))=g[j]$ for any $u\in [0,1]$. For simplicity, below we will drop $\Delta$ inside function arguments, i.e., $g(\Delta(j+u))\equiv g(j+u)$

The integral can thus be turned into \begin{eqnarray} \int f_m = \Delta \sum_{j_N} g[j_N] \int_0^1 dt \prod_{i=1}^{N-1} \Delta \sum_{j_i} \int_0^1 du g_i(j_i+u) h_i(j_N+t-j_i-u) \\ =\Delta^N \sum_{j_N} g[j_N] \int_0^1 dt \prod_{i=1}^{N-1} \sum_{j_i} g_i[j_i] \int_0^1 du\; h_i(j_N-j_i+t-u) \end{eqnarray} Since both $t$ and $u$ are in the interval $[0,1]$, $t-u$ is in the interval $[-1,1]$, and there are two possibilities,

we thus recognize \begin{eqnarray} \int_0^1 du\; h_i(j_N-j_i+t-u) = \int_0^1 du\; \left(\begin{array}{c} h_i[j_N-j_i] \Theta(0<u<t)) \\ h_i[j_N-j_i-1] \Theta(t<u<1) \end{array}\right) \\ =h_i[j_N-j_i-1] (1-t) + h_i[j_N-j_i] t %\\=h_i[j_N-j_i] + (h_i[j_N-j_i-1]-h_i[j_N-j_i])(1-t) \end{eqnarray}

We can define new arrays \begin{eqnarray} a_i[j_N] &=& \Delta \sum_j g_i[j] h_i[j_N-j] \\ b_i[j_N] &=& \Delta \sum_j g_i[j] (h_i[j_N-j_i-1]-h_i[j_N-j_i])\\ \end{eqnarray} to obtain the final integral \begin{eqnarray} I^N =\Delta \sum_{j_N} g[j_N] \int_0^1 dt \prod_{i=1}^{N-1} \left( a_i[j_N] + b_i[j_N] t \right) \end{eqnarray} which can be expanded to a polynomial of $N-1$-th order, and evaluated \begin{eqnarray} I^N =\Delta \sum_{j_N} g[jN] \left( \prod{i=1}^{N-1} a_i[j_N]+ \cdots

3D convolutions

Multiple 3D-type ingerals are very common, and deserve special attention. It is many times physically more appropriate to write the ansatz in terms of radial functions $g(|\vec{r}_i|)$ rather than use separate $g(x)$, $g(y)$, $g(z)$. We might want to use a convolution ansatz of the form \begin{eqnarray} f_m(\vec{r}_1,\vec{r}_2,\cdots,\vec{r}_n)= g_1(r_1) h_1(|\vec{r}_n-\vec{r}_1|) g_2(r_2)h_2(|\vec{r}_n-\vec{r}_2|)\cdots g_n(r_n) \end{eqnarray} The histogram for projections to $|\vec{r}_i|$ and $|\vec{r}_n-\vec{r}_i|$ is straighforward to obtain during simulation.

The analytic integral $\int f_m$ is somewhat more challenging for 3D case. We derive it below: \begin{eqnarray} \int f_m = \int d^3r_n g_n(r_n)\prod_i \int d^3 r_i g_i(r_i) h_i(|\vec{r}_n-\vec{r}_i|) \\= 4\pi \int d r_n r_n^2 g_n(r_n)\prod_{i=1}^{n-1}\int dr_i 2\pi r_i^2 \\ \int_{-1}^1 d(\cos\theta)g_i(r_i) h_i(\sqrt{r_n^2+r_i^2-2 r_n r_i \cos(\theta)}) \end{eqnarray}

\begin{eqnarray} \int f_m = 2 (2\pi)^n \int d r_n r_n^{3-n} g_n(r_n)\prod_{i=1}^{n-1}\int dr_i r_i g_i(r_i) \int_{|r_n-r_i|}^{r_n+r_i} dx\, x\, h_i(x) \end{eqnarray}

Next we integrate over pice-wise constant functions \begin{eqnarray} r_n = j_n + t\\ r_i = j + u \end{eqnarray} with $t\in[0,1]$ and $u\in[0,1]$, and we get \begin{eqnarray} &&\int f_m = 2 (2\pi\Delta^3)^n \sum_{j_n} g_n[j_n]\int_0^1 dt (j_n+t)^{3-n}\times \nonumber\\ && \times \prod_{i=1}^{n-1}\sum_j g_i[j] \int_0^1 du (j+u) \int_{|j_n-j+t-u|}^{j_n+j+t+u} dx\, x\, h_i(x) \end{eqnarray} Lets define \begin{eqnarray} F_i(j_n+t) \equiv \sum_j g_i[j] \int_0^1 du (j+u) \int_{|j_n-j+t-u|}^{j_n+j+t+u} dx\, x\, h_i(x) \end{eqnarray} and by power counting we can see that $F_i(t)$ is polynomial in $t$ of the forth order.

The derivation of polynomials $F_i[j_n](t)$ is somewhat tedious and will be given below in Sec.(Integrals $K_n$). But before that, lets us show how to compute the integral $\int f_m$ if the polynomials $F_i[j_n](t)$ are given: \begin{eqnarray} \int f_m = 2 (2\pi\Delta^3)^n \sum_{j_n} g_n[j_n]\int_0^1 dt (j_n+t)^{3-n} \prod_{i=1}^{n-1} F_i[j_n](t) \label{Eq:31} \end{eqnarray} Note that $\prod_{i=1}^{n-1} F_i[j_n](t)$ is a polynomial of $4(n-1)$-th order, and we write it as \begin{eqnarray} \prod_{i=1}^{n-1} F_i[j_n](t) \equiv \sum_{k=0}^{4(n-1)} a_k t^k \end{eqnarray} Such multiplication of polynomials can be straighforwardly implemented on the computer, and if polynomials $F_i[j_n](t)$ are known, the coefficients $a_k$ are easy to determine.

The integral of the function $f_m$ than requires one to evaluate \begin{eqnarray} K_n[j_n]= \int_0^1 dt (j_n+t)^{3-n} \sum_{k=0}^{4(n-1)} a_k t^k \end{eqnarray} so that the final result becomes \begin{eqnarray} \int f_m =2 (2\pi\Delta^3)^n \sum_{j_n} g_n[j_n] K_n[j_n] \end{eqnarray}

Integrals $K_n$

For $n \ge 3$ all powers in the integral $K_n$ are positive, and we can write \begin{eqnarray} K_n[j]= \sum_{p=0}^{3-n} {3-n\choose p}(j)^p \sum_{k=0}^{4(n-1)} a_k \int_0^1 dt\, t^{3-n-p+k} \\= \sum_{p=0}^{3-n} {3-n\choose p} (j)^p \sum_{k=0}^{4(n-1)} \frac{a_k}{4+k-n-p} \end{eqnarray} For $n=4$ we have \begin{eqnarray} K_n[j]= \sum_{k=0}^{4(n-1)} a_k \int_0^1 dt \frac{t^k}{j+t} \end{eqnarray} which can be computed by recursion. Let's define \begin{eqnarray} I_k[j]\equiv \int_0^1 dt \frac{t^k}{j+t} \end{eqnarray} so that \begin{eqnarray} K_n[j]= \sum_{k=0}^{4(n-1)} a_k I_k[j] \end{eqnarray} It is easy to see that $I_0[j]=\log(\frac{j+1}{j})$ and that $$I_{k+1}[j]=\int_0^1 t^k dt - j I_k[j]= \frac{1}{k+1} - j I_k[j]$$ which gives $$I_k[j] = \sum_{p=0}^{k-1}\frac{ (-j)^p}{k-p} + (-j)^k \log(\frac{j+1}{j})$$

Finally, for $n>4$, we define $$ I^p_k[j] = \int_0^1 \frac{t^k}{(j+t)^{p+1}} dt$$ so that $K_n$ can be simply computed by \begin{eqnarray} K_n[j_n]= \sum_{k=0}^{4(n-1)} a_k I_k^{n-4}[j_n] \end{eqnarray} and we derive the recursion relation for integrals $I_k^p[j]$. Notice that the previous case $n=4$ corresponds to $p=0$, but now we will derive recurison for any $n \ge 4$.

We notice that \begin{eqnarray} I_{k+1}^p= I_k^{p-1}- j\, I_k^p, \label{Eq:recursion} \end{eqnarray} which follows from the definition \begin{eqnarray} I^p_{k+1} = \int_0^1 \frac{t^k(t+j-j)}{(j+t)^{p+1}} dt = \int_0^1 \frac{t^k}{(j+t)^{p}} dt -j\int_0^1 \frac{t^k}{(j+t)^{p+1}} dt \end{eqnarray} This can be used to compute any $I_k^p$ if we know $I^{-1}_k$, which is simply given by $$I^{-1}_k=\frac{1}{k+1}$$ and $I_0^p$ for arbitrary $p$. We notice that \begin{eqnarray} I^{p>0}_{0} = \int_0^1 \frac{1}{(j+t)^{p+1}} dt = \frac{1}{p}\left(\frac{1}{j^p}-\frac{1}{(j+1)^p}\right) \end{eqnarray} and $I^0_0 = \log\left(\frac{j+1}{j}\right)$ hence the above recursion gives an arbitrary $I_k^p$.

To show that all $I_k^p$ are now determined, lets show for a first few $p$. For $p=0$ we get \begin{eqnarray} I^0_{k+1}= I^{-1}_k - j I^0_k = \frac{1}{k+1} - j I^0_k \end{eqnarray} which is exactly the recursion relation we derived earlier. Once all $I^0_k$ are determined, we can proceed to compute $I^1_{k}$ by the recursion \begin{eqnarray} I^1_{k+1}= I^{0}_k - j I^1_k \end{eqnarray} Since we know $I^1_0$, all $I^1_k$ can readily be obtained from this recursion relation. And once $I^1_k$ are known, $I^2_k$ are determined by the recursion relation specified above $I_{k+1}^p= I_k^{p-1}- j\, I_k^p$


Using Metropolis integration evaluate the following Linhardt function:

\begin{eqnarray} P(\Omega,\vec{q}) = -2 \int \frac{d^3k}{(2\pi)^3} \frac{f(\varepsilon_{\vec{k}+\vec{q}})-f(\varepsilon_{\vec{k}})}{\Omega-\varepsilon_{\vec{k}+\vec{q}}+\varepsilon_\vec{k}+i\delta} \end{eqnarray}

Here $\varepsilon_\vec{k}=k^2-k_F^2$.

We will use $k_F=\frac{(\frac{9\pi}{4})^{1/3}}{rs}$ with $r_s=2$, $f(x) = \frac{1}{\exp(x/T)+1}$, $T=0.02 k_F^2$, $\delta = 0.002 k_F^2$, and cutoff is $3 k_F$.

We will use linear mesh of frequency $\Omega=[0,k_F^2]$ with 100 points, and $q$ mesh $q=[0.1 k_F,0.2 k_F,0.3 k_F,0.4 k_F]$. We will plot the family of curves for $P(\Omega,q)$.

We will run for p.Nitt=5000000 steps.

The result for $P(\Omega=0,q<< k_F)$ should be close to $-n_F$, where $n_F = \frac{k_F}{2\pi^2}$.

We will use Linhardt function from previous homework, but slightly modify it, so that the input momentum contains $\vec{q}$ is the first vector and $\vec{k}$ as the second vector. (The dimension Ndim will thus be set to 2.) The function will return an array of values, corresponding to all frequencies.

The Metropolis algorithm has to be modified so that it takes an array of values. The probability will be proportional to the value at zero frequency $P(\Omega=0)$, but we will measure all frequencies at the same time.

Here we summarize all the necessary changes to the IntegrateByMetropolis3 routine

We marked all lines that require change by

# Next line needs CORRECTION for homework