深度學習與量子蒙特卡洛（1）

08-11

深度學習與量子蒙特卡洛（1）

10 人贊了文章

1）讓我們從一個簡單的高斯積分開始. 一維情況下：

$I_1 = frac{1}{sqrt{2pi}} int_{-infty}^infty ext{e}^{-frac{x^2}{2}} dx = 1. ~~~~~~( ext{Eq}~1.1)$

考慮 $x^2$ 的期望：

$<x^2> = frac{1}{sqrt{2pi}} int_{-infty}^infty x^2 ext{e}^{-frac{x^2}{2}} dx = 1. ~~~~~~( ext{Eq}~1.2)$

多維情況下：

$I_n = frac{1}{(sqrt{2pi})^n} int_{-infty}^infty ext{e}^{-frac{x_1^2+...+x_n^2}{2}} dx_1 ...dx_n = 1. ~~~~~~( ext{Eq}~1.3)$

現在我們嘗試用數值的方法計算方程 $( ext{Eq}~1.1)$ 。這裡我們使用蒙特卡洛方法，而不是通常的離散求和的方法。原因是對於以後將要討論的高維積分，通常的方法將會失效，而蒙特卡洛的方法依舊適用。

2）為了計算 $x^2$ 的期望，我們需要去按照一下分布隨機生成 $x$ 的值

$p(x)= frac{1}{sqrt{2pi}} ext{e}^{-frac{x^2}{2}} . ~~~~~~( ext{Eq}~2.1)$

例如產生了100000個 $x$ , 那麼 $x^2$ 的期望w值為：

$<x^2> = frac{1}{100000}sum_{i=1}^{100000} x_i^2 ~~~~~~( ext{Eq}~2.2)$

因而重點是如何按照分布 $p(x)$ 快速抽樣出100000個點。

這裡我們使用 Hybrid Monte-Carlo 方法來實現這個目的。（2017/12/28/0：27）

代碼詳見附錄B。可以看出無論是1維積分還是5維積分，蒙特卡洛的方法都可以在幾乎相同的時間內給出很好的結果。這是傳統的均勻畫網格的方法所不可及的。

3）總結

在物理中，我們經常遇到如下積分

$Z = int e^{-S(x)} dx ~~~~ ~~~~~~~~~~~~ ext{Eq 3.1}\ <O> = frac{1}{Z} int O(x)e^{-S(x)}dx ~~~~~~~~~~~~~~~~~~~~~ ext{Eq 3.2}$

當 $S(x)$ 為實數時，可以賦予因子 $e^{-S(x)}$ 概率解釋。利用蒙特卡洛方法產生分布 $p(x) sim e^{-S(x)}$ , 此時

$<O> = frac{1}{N}sum_{i=1}^N O(x_i)$

當 $S(x)$ 是複數時，概率解釋不成立，蒙特卡洛方法失效。這個問題被稱為符號問題。

如何解決符號問題？待續（2017/12/28/12：34）

附錄）

這裡我們給出上文提到的一些結果的代碼。

A. 對於高斯積分 (Mathematica code)：

s1 = Integrate[Exp[-x^2/2], {x, -Infinity, Infinity}]s2 = Integrate[x^2 Exp[-x^2/2], {x, -Infinity, Infinity}]I1 = s1/Sqrt[2 Pi]x2 = s2/s1

B. Hybrid Monte-Carlo 數值積分 (python3 code)：

# K.J. Sun 2017/12/28 copyright reserved.from __future__ import print_functionimport numpy as npimport matplotlib.pyplot as pltimport matplotlib.image as mpingfrom pyhmc import hmc # one can install hybrid monte-carlo by pip3 install pyhmcdef logprob(x, ivar): logp = -0.5 * np.sum(ivar * x**2) # log of probability function grad = -ivar * x # first derivative oflog of probability function return logp, gradivar = 1. # np.random.rand(5)dimx = 1 # dimension of xsamples = hmc(logprob, x0=np.random.randn(dimx), args=(ivar,), n_samples=100000)# generating 100000*dimx array of xx2 = np.mean(np.power(samples,2),0) # expectation value of x^2print(x2)

把該代碼保存為 hmc_integral.py, 然後運行 python3 hmc_integral.py,

dimx = 1 時輸出結果為 0.966

dimx = 2 時輸出結果為 [ 1.017 0.967]

dimx = 5 時輸出結果為 [ 1.02 1.036 0.983 0.989 1.022]