# 経路積分 Monte Carlo：Feynman–Kac 公式・時間分割の誤差・推定量と符号問題 > 量子系の有限温度の計算を Monte Carlo 法で行うときの誤差を証明つきで扱う。対称分割の誤差 O(β³/M²) の交換子による評価、調和振動子の離散化誤差の閉じた形、Brown 橋の二分法、熱力学的推定量とビリアル推定量の分散、局所更新の緩和の M² 則、符号問題の指数則。 > https://rikai.mugen-giken.com/computer-science/stochastic-computation/path-integral-monte-carlo ## 0. この記事の要点 量子力学で従う粒子の温度を持つ状態を計算機で調べたいとき、状態の数が多すぎて直接は扱えないことがよくあります。この記事の方法は、1 個の量子的な粒子を、小さなばねで輪につないだ多数の古典的な粒子（ビーズ）に置き換え、輪の形をサイコロで次々に変えながら平均を取るというものです。読み終えると、ビーズの数を増やしたときに誤差がどれだけ減るか、どの量の平均を取ると計算が安定するか、どんな場合に方法が破綻するかを説明できるようになります。 - 量子系の分配関数 $Z = \mathrm{Tr}\,e^{-\beta H}$ は、$e^{-\beta H}$ を $M$ 個の短い時間発展の積に分けることで、$M$ 個の配位からなる閉じた経路（環状の鎖）についての古典的な積分になります。経路積分 Monte Carlo 法は、その積分を Markov 連鎖 Monte Carlo 法で求めます。 - 分割の誤差は交換子（$[X,Y]=XY-YX$、掛ける順序を入れ替えたときの差）で書けます。$\|e^{-sA}\| \le 1$、$\|e^{-sB}\| \le 1$、$\|e^{-s(A+B)}\| \le 1$ のとき、対称分割の誤差は $(\beta^3/M^2)(\|[A,[A,B]]\|/24 + \|[B,[B,A]]\|/12)$ 以下、対称でない分割では $(\beta^2/2M)\|[A,B]\|$ 以下です（）。分配関数はトレースの巡回性からどちらの分割でも同じ値で、誤差は $O(\beta^3/M^2)$ です（）。 - 調和振動子では $M$ 分割の分配関数が閉じた形で求まり、エネルギーの離散化誤差の先頭項は $(\beta\omega)^2/M^2$ に比例します（）。 - 両端を固定した自由粒子の経路（Brown 橋。始点と終点を決めたブラウン運動の経路のこと）は、中点を正規分布から選ぶことを繰り返す二分法で、棄却なしに生成できます（）。ポテンシャルは Metropolis の採択確率（新しい配位を受け入れるかどうかを決める確率）にだけ現れます。 - エネルギーの熱力学的推定量（$\ln Z$ の微分から素直に作った平均の取り方）の分散は $M$ に比例して増えます。ビリアル推定量（部分積分で書き換えた別の平均の取り方）は同じ期待値をもち、分散が $M$ について有界です（）。隣り合う 1 点だけを動かす更新は、緩和に $M^2$ に比例する掃引回数を要します（）。 - 重みが符号を変えるとき、平均符号は $e^{-\beta\Delta F}$ で、相対誤差を一定に保つのに要る標本数は $e^{2\beta\Delta F}$ に比例します（）。 ## 1. 動機 温度 $1/\beta$ の量子系の平衡状態は、ハミルトニアン $H$ から $e^{-\beta H}/Z$ で決まります。$H$ を対角化できれば話は終わりですが、粒子数とともに状態空間の次元が指数的に増えるので、対角化できるのは小さな系に限られます。一方、$H$ が運動エネルギーとポテンシャルエネルギーの和 $A + B$ で、それぞれの $e^{-\tau A}$、$e^{-\tau B}$ なら具体的に書ける、という状況はよくあります。$A$ と $B$ が交換しないので $e^{-\beta(A+B)} \ne e^{-\beta A}e^{-\beta B}$ ですが、$\beta$ を $M$ 個に分けて $$ e^{-\beta(A+B)} \approx \bigl(e^{-\tau A}e^{-\tau B}\bigr)^M, \qquad \tau = \frac\beta M $$ とすれば、$M \to \infty$ で誤差が消えます（Lie–Trotter の積公式、）。積の間に完全系を $M$ 回はさむと、トレースは $M$ 個の配位 $x_1, \ldots, x_M$ についての積分になり、被積分関数が正なら確率分布と見なせます。これが経路積分 Monte Carlo 法（path integral Monte Carlo）です。ここで「完全系をはさむ」とは、各積の間に恒等作用素（位置の固有状態についての和 $\int|x\rangle\langle x|\,dx$）を入れることです。すると各因子は、位置 $x_k$ から $x_{k+1}$ への行列要素という普通の数になり、トレースはそれらの積を $x_1,\ldots,x_M$ について積分したものになります。 この方法の誤差には、性質の違う 2 種類があります。第 1 は、有限の $M$ で打ち切ることによる**離散化誤差**（偏り）で、標本をいくら増やしても消えません。第 2 は Monte Carlo 法の**統計誤差**で、こちらは推定量の選び方と連鎖の更新の仕方で大きく変わります。厄介なのは、$M$ を増やして離散化誤差を減らすと、素朴な推定量の分散が増え、素朴な更新の緩和が遅くなることです。この章では、離散化誤差の次数（§3・§4）、経路の効率のよい生成（§5）、推定量の分散（§6）、更新の緩和（§7）を順に定理として評価し、最後に重みが正でない場合に起きる符号問題（§8）を扱います。 収束そのもの（次数のつかない極限）は既存の章にあります。Euclid 時間分割公式と転送行列は[Euclid 経路積分と転送行列](/physics/lattice-field-theory/euclidean-path-integral)、Feynman–Kac 公式は[熱方程式](/mathematics/pde/heat-equation)（）と [Feller 過程](/mathematics/probability/markov-processes-feller)（）、非有界な作用素の Trotter 積公式は です。この章はそれらを再証明せず、誤差の大きさを扱います。単位は $\hbar = 1$、粒子の質量は $1$ とします。 ## 2. 準備 ### 2.1. 分配関数と熱平均 $H$ を $n$ 次のエルミート行列、$\beta > 0$ とする。 $$ Z(\beta) = \mathrm{Tr}\,e^{-\beta H}, \qquad \langle O\rangle_\beta = \frac{\mathrm{Tr}(O\,e^{-\beta H})}{Z(\beta)} \quad (O \in M_n(\mathbb{C})) $$ を**分配関数**、$O$ の**熱平均**という。$E(\beta) = \langle H\rangle_\beta = -\dfrac{d}{d\beta}\ln Z(\beta)$ を内部エネルギーという。 行列指数関数 $e^{X} = \sum_kX^k/k!$ とその性質（指数法則、$\frac{d}{dt}e^{tX} = Xe^{tX} = e^{tX}X$、相似変換との交換）は と を使います。$H$ の固有値を $E_1 \le \cdots \le E_n$ とすると $Z = \sum_je^{-\beta E_j} > 0$ です。 行列のノルムは、$\mathbb{C}^n$ のユークリッドノルムから誘導した作用素ノルム $\|X\| = \max_{|v|=1}|Xv|$（）とします。$\|XY\| \le \|X\|\|Y\|$、$\|I\| = 1$ です。対角成分は $|X_{ii}| = |\langle e_i, Xe_i\rangle| \le \|X\|$ を満たすので $|\mathrm{Tr}\,X| \le n\|X\|$ です。行列値の連続関数の積分は成分ごとに定め、$\|\int_a^bX(s)\,ds\| \le \int_a^b\|X(s)\|\,ds$ が成り立ちます（Riemann 和についての三角不等式の極限）。交換子を $[X, Y] = XY - YX$ と書きます。 ### 2.2. 縮小半群を生成する行列 誤差評価で必要なのは「短い時間発展を何回掛けても誤差が増幅されない」ことです。それを次の条件で表します。 $A \in M_n(\mathbb{C})$ が次のどちらかを満たせば、すべての $s \ge 0$ で $\|e^{-sA}\| \le 1$ である。 1. $A$ はエルミートで半正定値（固有値がすべて $0$ 以上）。 2. $A$ は反エルミート（$A^* = -A$）。このとき $e^{-sA}$ はユニタリで $\|e^{-sA}\| = 1$。 どちらの性質も和で保たれる。したがって $A, B$ がともに 1 を満たす（またはともに 2 を満たす）とき、$A$、$B$、$A + B$ のすべてについて結論が成り立つ。 **1.** スペクトル定理（）から $A = U\Lambda U^*$（$U$ はユニタリ、$\Lambda = \mathrm{diag}(\lambda_j)$、$\lambda_j \ge 0$）です。 の 4 から $e^{-sA} = Ue^{-s\Lambda}U^*$ で、$e^{-s\Lambda} = \mathrm{diag}(e^{-s\lambda_j})$ です。$|v| = 1$ のとき $w = U^*v$ も $|w| = 1$ で、$|e^{-sA}v|^2 = \sum_je^{-2s\lambda_j}|w_j|^2 \le \sum_j|w_j|^2 = 1$ です。 **2.** 級数の各項の随伴を取ると $(e^{-sA})^* = e^{-sA^*} = e^{sA}$ で、 の 2 から $(e^{-sA})^*e^{-sA} = e^{sA}e^{-sA} = I$ です。よって $|e^{-sA}v|^2 = \langle v, (e^{-sA})^*e^{-sA}v\rangle = |v|^2$ です。 半正定値エルミート行列の和は半正定値エルミート（$\langle v, (A+B)v\rangle = \langle v, Av\rangle + \langle v, Bv\rangle \ge 0$）、反エルミート行列の和は反エルミートです。 1 は虚時間の発展 $e^{-\tau H}$（この章の主題）、2 は実時間の発展 $e^{-itH}$（$A = iH_1$、$B = iH_2$）に当たります。次の定理は両方を同時に扱います。 $X, Y \in M_n(\mathbb{C})$、$s \ge 0$ について $$ [Y, e^{-sX}] = \int_0^se^{-(s-u)X}\,[X, Y]\,e^{-uX}\,du. $$ したがって、すべての $u \ge 0$ で $\|e^{-uX}\| \le 1$ ならば $\|[Y, e^{-sX}]\| \le s\,\|[X, Y]\|$。 $G(s) = Ye^{-sX} - e^{-sX}Y$ とおきます。$G(0) = 0$ で、$\frac{d}{ds}e^{-sX} = -Xe^{-sX}$（ の 3）から $$ G'(s) = -YXe^{-sX} + Xe^{-sX}Y = -X\bigl(Ye^{-sX} - e^{-sX}Y\bigr) + (XY - YX)e^{-sX} = -XG(s) + [X, Y]e^{-sX} $$ です。これは $G$ の各列についての非斉次線形方程式 $g' = -Xg + b(s)$ で、斉次方程式の基本解行列は $e^{-sX}$ です。定数変化法（）と $e^{-sX}(e^{-uX})^{-1} = e^{-(s-u)X}$ から $$ G(s) = \int_0^se^{-(s-u)X}[X, Y]e^{-uX}\,du $$ です。後半は、$0 \le u \le s$ で $\|e^{-(s-u)X}\| \le 1$、$\|e^{-uX}\| \le 1$ であることと、積分のノルムの評価から従います。 ## 3. 時間分割の誤差 ### 3.1. 対称分割と非対称分割 $A, B \in M_n(\mathbb{C})$、$\beta \ge 0$、$M \in \mathbb{N}$、$\tau = \beta/M$ とする。 $$ S_\tau = e^{-\tau A/2}e^{-\tau B}e^{-\tau A/2}, \qquad L_\tau = e^{-\tau A}e^{-\tau B} $$ とおき、$S_\tau^M$ を $e^{-\beta(A+B)}$ の**対称分割**（$M$ 分割の原始近似、Strang 分割）、$L_\tau^M$ を**非対称分割**（Lie–Trotter 分割）という。$\tau$ を**時間刻み**という。 2 つの分割は、$A$ と $B$ の短い時間発展を交互に並べる点は同じで、並べ方だけが違います。下の図は $M=3$ の場合の並びです。
ABA BAB A/2BA BABA/2
$A, B \in M_n(\mathbb{C})$ が、すべての $s \ge 0$ で $$ \|e^{-sA}\| \le 1, \qquad \|e^{-sB}\| \le 1, \qquad \|e^{-s(A+B)}\| \le 1 $$ を満たすとする。$\beta \ge 0$、$M \in \mathbb{N}$、$\tau = \beta/M$ について次が成り立つ。 1. （対称分割） $$ \bigl\|e^{-\beta(A+B)} - \bigl(e^{-\tau A/2}e^{-\tau B}e^{-\tau A/2}\bigr)^M\bigr\| \le \frac{\beta^3}{M^2}\Bigl(\frac{\|[A,[A,B]]\|}{24} + \frac{\|[B,[B,A]]\|}{12}\Bigr). $$ 2. （非対称分割） $$ \bigl\|e^{-\beta(A+B)} - \bigl(e^{-\tau A}e^{-\tau B}\bigr)^M\bigr\| \le \frac{\beta^2}{2M}\,\|[A, B]\|. $$ 仮定は、$A, B$ がともに半正定値エルミート行列のとき、およびともに反エルミート行列のときに満たされる（）。 証明の見通しです。全体の誤差は、1 ステップの誤差（局所誤差）を $M$ 個足したもので抑えられます。1 ステップの誤差は、「厳密な発展 $e^{-\tau(A+B)}$ が満たす方程式」からのずれを時間で積分して求めます。対称分割では、ずれの先頭の項（$A$ と $B$ の交換子に比例する項）が左右で打ち消し合い、残りが 2 重交換子になります。それで次数がひとつ上がります。 $H = A + B$、$C = [A, B]$ とおきます。 **テレスコープの和。** $F$ を $\|F\| \le 1$ である行列とすると $$ F^M - e^{-\beta H} = F^M - (e^{-\tau H})^M = \sum_{k=0}^{M-1}F^k\bigl(F - e^{-\tau H}\bigr)e^{-(M-1-k)\tau H} $$ です（右辺は隣り合う項が打ち消し合います）。$\|F^k\| \le 1$、$\|e^{-(M-1-k)\tau H}\| \le 1$ なので $$ \|F^M - e^{-\beta H}\| \le M\,\|F - e^{-\tau H}\| $$ です。$S_\tau$ と $L_\tau$ は縮小の積なのでノルムが $1$ 以下です。したがって 1 ステップの誤差（局所誤差）を評価すればよいことになります。 **局所誤差の表示。** $t \ge 0$ の行列値関数 $F(t)$ が $F(0) = I$ で微分可能とし、$R(t) = F'(t) + HF(t)$（方程式 $F' = -HF$ からのずれ）とおきます。$F' = -HF + R$ に定数変化法（。基本解行列は $e^{-tH}$）を列ごとに使うと $$ F(t) - e^{-tH} = \int_0^te^{-(t-s)H}R(s)\,ds, \qquad \|F(t) - e^{-tH}\| \le \int_0^t\|R(s)\|\,ds $$ です（Duhamel の公式）。 **2 の証明。** $F(t) = e^{-tA}e^{-tB}$ とします。$F'(t) = -Ae^{-tA}e^{-tB} - e^{-tA}Be^{-tB}$ なので $$ R(t) = F' + (A + B)F = Be^{-tA}e^{-tB} - e^{-tA}Be^{-tB} = [B, e^{-tA}]\,e^{-tB} $$ です。 から $\|R(t)\| \le t\|[A, B]\|$ で、$\|L_\tau - e^{-\tau H}\| \le \int_0^\tau s\,ds\,\|[A,B]\| = \tau^2\|[A,B]\|/2$ です。$M$ 倍して $M\tau^2/2 = \beta^2/(2M)$ です。 **1 の証明。** $F(t) = e^{-tA/2}e^{-tB}e^{-tA/2}$ とします。 $$ F'(t) = -\frac A2F(t) - e^{-tA/2}Be^{-tB}e^{-tA/2} - F(t)\frac A2 $$ なので $$ R(t) = \frac A2F - F\frac A2 + BF - e^{-tA/2}Be^{-tB}e^{-tA/2} = \underbrace{[B, e^{-tA/2}]\,e^{-tB}e^{-tA/2}}_{R_1(t)} + \underbrace{\tfrac12e^{-tA/2}\,[A, e^{-tB}]\,e^{-tA/2}}_{R_2(t)} $$ です（$A$ は $e^{-tA/2}$ と交換するので、$\frac A2F - F\frac A2 = \frac12e^{-tA/2}(Ae^{-tB} - e^{-tB}A)e^{-tA/2}$）。$R_1$ と $R_2$ はそれぞれ $t$ の 1 次の大きさですが、先頭の項が打ち消し合います。それを、逆行列 $e^{+tA}$ を使わずに取り出します。 $R_1$ について。（$X = A$、$Y = B$）から $[B, e^{-tA/2}] = \int_0^{t/2}e^{-(t/2-u)A}Ce^{-uA}\,du$ です。$Ce^{-uA} = e^{-uA}C + [C, e^{-uA}]$ と書くと、$e^{-(t/2-u)A}e^{-uA} = e^{-tA/2}$ なので $$ R_1(t) = \frac t2\,e^{-tA/2}Ce^{-tB}e^{-tA/2} + E_1(t), \qquad E_1(t) = \int_0^{t/2}e^{-(t/2-u)A}\,[C, e^{-uA}]\,du\;e^{-tB}e^{-tA/2} $$ です。（$X = A$、$Y = C$）から $\|[C, e^{-uA}]\| \le u\|[A, C]\|$ なので $$ \|E_1(t)\| \le \int_0^{t/2}u\,du\,\|[A, C]\| = \frac{t^2}{8}\|[A,[A,B]]\|. $$ $R_2$ について。（$X = B$、$Y = A$）から $[A, e^{-tB}] = \int_0^te^{-(t-v)B}[B, A]e^{-vB}\,dv = -\int_0^te^{-(t-v)B}Ce^{-vB}\,dv$ です。$e^{-(t-v)B}C = Ce^{-(t-v)B} - [C, e^{-(t-v)B}]$ と書くと $$ R_2(t) = -\frac t2\,e^{-tA/2}Ce^{-tB}e^{-tA/2} + E_2(t), \qquad E_2(t) = \frac12e^{-tA/2}\int_0^t[C, e^{-(t-v)B}]\,e^{-vB}\,dv\;e^{-tA/2} $$ です。（$X = B$、$Y = C$）から $\|[C, e^{-(t-v)B}]\| \le (t-v)\|[B, C]\|$ で、$[B, C] = [B,[A,B]] = -[B,[B,A]]$ なので $$ \|E_2(t)\| \le \frac12\int_0^t(t - v)\,dv\,\|[B, C]\| = \frac{t^2}{4}\|[B,[B,A]]\|. $$ $R_1$ と $R_2$ の第 1 項は打ち消し合い、$R(t) = E_1(t) + E_2(t)$ です。したがって $$ \|S_\tau - e^{-\tau H}\| \le \int_0^\tau\Bigl(\frac{s^2}{8}\|[A,[A,B]]\| + \frac{s^2}{4}\|[B,[B,A]]\|\Bigr)ds = \tau^3\Bigl(\frac{\|[A,[A,B]]\|}{24} + \frac{\|[B,[B,A]]\|}{12}\Bigr) $$ で、$M$ 倍して $M\tau^3 = \beta^3/M^2$ です。 対称に並べるだけで次数が 1 つ上がります。手間はほとんど変わりません。$S_\tau^M = e^{-\tau A/2}(e^{-\tau B}e^{-\tau A})^{M-1}e^{-\tau B}e^{-\tau A/2}$ なので、対称分割は非対称分割の両端を半分ずつにしたものです。$A$ と $B$ が交換すれば誤差は $0$ で、上界もそうなっています。誤差を決めるのは $\beta$ と $M$ だけではなく交換子の大きさで、たとえばポテンシャルが急に変わる系では $[A,[A,B]]$（運動エネルギーとポテンシャルの 2 重交換子で、ポテンシャルの高階導関数を含みます）が大きく、細かい分割が要ります。 (a) **半正定値でないエルミート行列。** $A$ の最小固有値を $a_0$、$B$ の最小固有値を $b_0$ とすると、$A - a_0I$ と $B - b_0I$ は半正定値です。定数を引いても交換子は変わらず、$e^{-s(A - a_0I)} = e^{sa_0}e^{-sA}$ なので、定理を $A - a_0I$、$B - b_0I$ に使って両辺に $e^{-\beta(a_0+b_0)}$ を掛けると、任意のエルミート行列 $A, B$ について右辺に因子 $e^{-\beta(a_0 + b_0)}$ を付けた評価が得られます。 (b) **実時間の発展。** $H_1, H_2$ がエルミートのとき、$A = iH_1$、$B = iH_2$ は反エルミートで、$\beta$ を時刻 $t$ と読めば $\|e^{-it(H_1+H_2)} - (e^{-i\tau H_1/2}e^{-i\tau H_2}e^{-i\tau H_1/2})^M\| \le t\tau^2(\|[H_1,[H_1,H_2]]\|/24 + \|[H_2,[H_2,H_1]]\|/12)$ です（$i$ の 3 乗は絶対値 $1$）。 (c) **非有界な作用素。** 連続空間の $H = -\frac12\Delta + V$ では、運動エネルギーもポテンシャルも有界な行列ではなく、交換子のノルムは有限ではありません。このとき $M \to \infty$ での収束（次数なし）は Trotter 積公式（）と Euclid 時間分割公式（：$V$ についての適当な仮定のもとで、時間分割した積分核が $e^{-\tau H}$ に強収束する）が保証します。次数つきの評価は、作用させるベクトルの滑らかさ（交換子がそのベクトルに作用した結果が有限であること）を仮定すれば得られます。Schrödinger 方程式の Strang 分割についての T. Jahnke, C. Lubich, *BIT* 40 (2000) がその形の定理です。この記事では、格子や基底で切断して有限次元にした系だけを定理にします。§9 の実験で、切断の次元を上げると交換子のノルムが大きくなり、上界が実際の誤差からかけ離れることを見ます。 (d) **高次の分割。** $[B,[A,B]]$ がポテンシャルの勾配の 2 乗になることを使ってポテンシャルを補正する方法（M. Takahashi, M. Imada, *J. Phys. Soc. Jpn.* 53 (1984)、S. A. Chin, *Phys. Lett. A* 226 (1997)）は、トレースの誤差を $O(\tau^4)$ にします。 ### 3.2. 分配関数と熱平均 $A, B$ を $n$ 次の半正定値エルミート行列、$H = A + B$、$Z = \mathrm{Tr}\,e^{-\beta H}$ とし、 $$ Z_M = \mathrm{Tr}\,S_\tau^M, \qquad \varepsilon_M = \frac{\beta^3}{M^2}\Bigl(\frac{\|[A,[A,B]]\|}{24} + \frac{\|[B,[B,A]]\|}{12}\Bigr) $$ とおく（$\tau = \beta/M$）。 1. $\mathrm{Tr}\,L_\tau^M = \mathrm{Tr}\,S_\tau^M = Z_M$。すなわち、対称分割と非対称分割の分配関数は一致する。 2. $|Z_M - Z| \le n\,\varepsilon_M$。 3. $O \in M_n(\mathbb{C})$ について $\langle O\rangle_M = \mathrm{Tr}(O\,S_\tau^M)/Z_M$ とおく。$n\,\varepsilon_M \le Z/2$ ならば $$ \bigl|\langle O\rangle_M - \langle O\rangle_\beta\bigr| \le \frac{4n\,\|O\|}{Z}\,\varepsilon_M. $$ **1.** $S_\tau^M = e^{-\tau A/2}(e^{-\tau B}e^{-\tau A})^{M-1}e^{-\tau B}e^{-\tau A/2}$ です。トレースの巡回性 $\mathrm{Tr}(XY) = \mathrm{Tr}(YX)$ で右端の $e^{-\tau A/2}$ を左端に回すと $\mathrm{Tr}\,S_\tau^M = \mathrm{Tr}\bigl(e^{-\tau A}(e^{-\tau B}e^{-\tau A})^{M-1}e^{-\tau B}\bigr) = \mathrm{Tr}(e^{-\tau A}e^{-\tau B})^M$ です。 **2.** の 1 と $|\mathrm{Tr}\,X| \le n\|X\|$ から $|Z_M - Z| = |\mathrm{Tr}(S_\tau^M - e^{-\beta H})| \le n\,\varepsilon_M$ です。 **3.** $S_\tau = (e^{-\tau B/2}e^{-\tau A/2})^*(e^{-\tau B/2}e^{-\tau A/2})$ は半正定値エルミートなので $S_\tau^M$ もそうで、$Z_M$ は実数です。2 と仮定から $Z_M \ge Z - n\varepsilon_M \ge Z/2 > 0$ です。次に $|\langle O\rangle_\beta| \le \|O\|$ を示します。$e^{-\beta H}$ の固有分解を $\sum_jp_jv_jv_j^*$（$p_j = e^{-\beta E_j} > 0$、$v_j$ は正規直交）とすると $\mathrm{Tr}(Oe^{-\beta H}) = \sum_jp_j\langle v_j, Ov_j\rangle$ で、絶対値は $\|O\|\sum_jp_j = \|O\|Z$ 以下です。そこで $$ \langle O\rangle_M - \langle O\rangle_\beta = \frac{\mathrm{Tr}\bigl(O(S_\tau^M - e^{-\beta H})\bigr)}{Z_M} + \langle O\rangle_\beta\,\frac{Z - Z_M}{Z_M} $$ と分けると、第 1 項の絶対値は $n\|O\|\varepsilon_M/Z_M$ 以下、第 2 項は $\|O\|\cdot n\varepsilon_M/Z_M$ 以下で、合わせて $2n\|O\|\varepsilon_M/Z_M \le 4n\|O\|\varepsilon_M/Z$ です。 1 は、経路積分 Monte Carlo 法で分配関数やエネルギーの誤差が $\tau^2$ の次数になる理由です。実際の計算では $e^{-\tau A}e^{-\tau B}$ の形の重みを使いますが、トレースを取る限り対称分割と同じものです。観測量についても、$O$ が $A$ と交換する場合（$A$ をポテンシャル、$O$ を位置の関数とする場合がこれに当たります）は、$\mathrm{Tr}(OL_\tau^M) = \mathrm{Tr}(OS_\tau^M)$ が同じ巡回の議論で成り立ち、誤差は $O(\tau^2)$ です。$A$ と交換しない観測量を非対称分割のまま測ると、誤差は $O(\tau)$ に落ちます。 誤差が $\tau^2$ の偶数冪で展開できることを使えば、2 つの $M$ での結果から $\tau \to 0$ への外挿ができます。外挿の誤差の扱いは [Richardson 補外](/mathematics/numerical-analysis/extrapolation-and-romberg)と同じです。 $\sigma_x = \begin{pmatrix}0&1\\1&0\end{pmatrix}$、$\Gamma > 0$ とし、$H = -\Gamma\sigma_x$（横磁場の中のスピン $1/2$）を考えます。$\sigma_z$ の固有状態 $s = \pm1$ を基底に取ると、$\sigma_x^2 = I$ から $$ e^{\tau\Gamma\sigma_x} = \cosh(\tau\Gamma)\,I + \sinh(\tau\Gamma)\,\sigma_x $$ で、行列要素は「同じ状態にとどまる」が $\cosh(\tau\Gamma)$、「反転する」が $\sinh(\tau\Gamma)$ です。$Z = \mathrm{Tr}(e^{\tau\Gamma\sigma_x})^M$ は、周期的な列 $s_1, \ldots, s_M$ についての和になり、反転（キンク）が $k$ 回ある列の重みは $\sinh^k(\tau\Gamma)\cosh^{M-k}(\tau\Gamma) > 0$ です。周期的なので $k$ は偶数で、キンクの位置の選び方が $\binom Mk$ 通り、出発点が 2 通りです。$M \to \infty$（$\tau = \beta/M$）で $$ \binom Mk\sinh^k\Bigl(\frac{\beta\Gamma}{M}\Bigr)\cosh^{M-k}\Bigl(\frac{\beta\Gamma}{M}\Bigr) \to \frac{(\beta\Gamma)^k}{k!} $$ なので、キンクの個数は平均 $\beta\Gamma$ の Poisson 分布を偶数に制限したものに従い、$Z = 2\sum_{k\text{ 偶数}}(\beta\Gamma)^k/k! = 2\cosh(\beta\Gamma)$ が再現されます。虚時間の経路は、強度 $\Gamma$ の Poisson 過程（）の時刻で反転する階段関数です。この系では分割する相手がないので離散化誤差はありません。縦磁場 $-h\sigma_z$ を加えると $A = -h\sigma_z$、$B = -\Gamma\sigma_x$ の分割になり、重みに $\prod_ke^{\tau hs_k}$ が掛かって、誤差は の形になります。連続時間の極限を直接標本にする方法（連続時間 Monte Carlo 法）では、この誤差もなくなります。 ## 4. 環状の鎖としての経路 ### 4.1. 連続空間の離散化した作用 1 次元の粒子 $H = -\frac12\frac{d^2}{dx^2} + V(x)$ を考えます。運動エネルギーの部分の $e^{-\tau H_0}$ は熱核 $$ p(\tau, x - y) = \frac{1}{\sqrt{2\pi\tau}}\exp\Bigl(-\frac{(x-y)^2}{2\tau}\Bigr) $$ を積分核とする作用素で、$p(\tau, \cdot)$ は時間 $\tau$ の間の Brown 運動（）の変位の密度です。熱方程式の章の表示（）は、生成子が $\Delta$ になる規約（$X_t = x + \sqrt2B_t$、熱核 $\Phi(y, t) = (4\pi t)^{-1/2}e^{-y^2/(4t)}$）で書かれているので、$p(\tau, y) = \Phi(y, \tau/2)$、$e^{-\tau(-\frac12\Delta + V)} = e^{-(\tau/2)(-\Delta + 2V)}$ と読み替えます（そこでの時刻がここの $\tau/2$、ポテンシャルがここの $2V$ に当たります）。ポテンシャルの部分 $e^{-\tau V}$ は掛け算です。分割 $(e^{-\tau V/2}e^{-\tau H_0}e^{-\tau V/2})^M$ のトレースを積分核で書くと、次の量になります。以下、$V$ について次を仮定します。 - **仮定 (C)**：$V : \mathbb{R} \to \mathbb{R}$ は $C^1$ 級で、定数 $c_0 > 0$、$c_1, c_2 \ge 0$、$r \ge 1$ があって $V(x) \ge c_0x^2 - c_1$、$|V'(x)| \le c_2(1 + |x|^r)$。 $\beta > 0$、$M \in \mathbb{N}$、$\tau = \beta/M$ とする。$x = (x_1, \ldots, x_M) \in \mathbb{R}^M$（$x_{M+1} = x_1$ とする）について $$ S_M(x) = \sum_{k=1}^M\Bigl[\frac{(x_{k+1} - x_k)^2}{2\tau} + \tau V(x_k)\Bigr], \qquad Z_M = \frac{1}{(2\pi\tau)^{M/2}}\int_{\mathbb{R}^M}e^{-S_M(x)}\,dx $$ を**離散化した作用**、$M$ 分割の分配関数といい、確率密度 $$ \pi_M(x) = \frac{e^{-S_M(x)}}{(2\pi\tau)^{M/2}Z_M} $$ に従う $x$ を**環状の鎖**（ring polymer）、各 $x_k$ を**ビーズ**（時間切片）という。$\pi_M$ についての期待値を $\langle\cdot\rangle_M$ と書く。 仮定 (C) から $S_M(x) \ge \tau c_0|x|^2 - \beta c_1$ なので、積分は収束します。$\pi_M$ は、隣どうしがばね定数 $1/\tau$ のばねでつながれ、各ビーズが外場 $\tau V$ を受ける古典的な環状の鎖の Boltzmann 分布です。1 個の量子的な粒子が、$M$ 個の古典的な粒子の環に化けたことになります。$M$ を増やすとばねは硬くなり、外場は弱くなります。 これは[格子上の Euclid 経路積分](/physics/lattice-field-theory/euclidean-path-integral)の周期的な作用と同じものです。そこでの転送行列 $T_a = e^{-aV/2}e^{-aH_0}e^{-aV/2}$（）を使うと、$Z_M = \mathrm{Tr}\,T_\tau^M$ です。$M \to \infty$ の極限では、$\pi_M$ のもとでの経路は、重み $\exp(-\int_0^\beta V(X_s)\,ds)$ を付けた周期的な Brown 運動に近づきます。これが Feynman–Kac 公式（、）の内容です。ただし、引用した 2 つの定理は有界なポテンシャルについてのものです。仮定 (C) のような非有界なポテンシャルでの極限は、この章の定理では使いません（以下の主張はすべて有限の $M$ についてのものです）。 ### 4.2. 調和振動子の離散化誤差 調和振動子では $Z_M$ が閉じた形で求まるので、離散化誤差の次数と係数を具体的に確かめられます。 $V(x) = \frac12\omega^2x^2$（$\omega > 0$）とし、$\theta > 0$ を $\sinh(\theta/2) = \tau\omega/2$ で定める（$\tau = \beta/M$）。 1. $Z_M = \dfrac{1}{2\sinh(M\theta/2)}$。 2. $E_M = -\dfrac{\partial}{\partial\beta}\ln Z_M$ は $$ E_M = \frac{\omega}{2\cosh(\theta/2)}\coth\frac{M\theta}{2} $$ である。また $\langle x_k^2\rangle_M = E_M/\omega^2$。 3. $\beta, \omega$ を固定して $M \to \infty$ とするとき、$E = \frac\omega2\coth\frac{\beta\omega}{2}$（厳密な内部エネルギー）、$x = \beta\omega/2$ として $$ E_M - E = -\frac\omega2\Bigl(\frac{\beta\omega}{M}\Bigr)^2\Bigl[\frac{\coth x}{8} - \frac{x}{24\sinh^2x}\Bigr] + O(M^{-4}) $$ で、角括弧の中は正である。 **1.** $S_M(x) = \frac{1}{2\tau}\sum_k[(x_{k+1} - x_k)^2 + \tau^2\omega^2x_k^2]$ は、格子間隔 $a = \tau$、質量 $1$ の格子調和振動子の周期的な作用で、$\cosh\theta = 1 + \tau^2\omega^2/2$ は $\sinh(\theta/2) = \tau\omega/2$ と同値です（$\cosh\theta = 1 + 2\sinh^2(\theta/2)$）。 から $Z_M = 1/(2\sinh(M\theta/2))$ です。（証明の筋は、$S_M = \frac{1}{2\tau}x^{\mathsf T}Qx$ の行列 $Q$ が巡回行列で固有値が $2\cosh\theta - 2\cos(2\pi j/M)$、Gauss 積分から $Z_M = (\det Q)^{-1/2}$、$\det Q = \prod_j(2\cosh\theta - 2\cos(2\pi j/M)) = (2\sinh(M\theta/2))^2$ です。） **2.** $\ln Z_M = -\ln(2\sinh(M\theta/2))$ で、$\theta$ は $\tau = \beta/M$ を通じて $\beta$ によります。$\sinh(\theta/2) = \tau\omega/2$ を $\tau$ で微分すると $\frac12\cosh(\theta/2)\,\frac{d\theta}{d\tau} = \frac\omega2$、すなわち $\frac{d\theta}{d\tau} = \omega/\cosh(\theta/2)$ で、$\frac{\partial\theta}{\partial\beta} = \frac1M\frac{d\theta}{d\tau}$ です。よって $$ E_M = \frac M2\coth\frac{M\theta}{2}\cdot\frac{\partial\theta}{\partial\beta} = \frac{\omega}{2\cosh(\theta/2)}\coth\frac{M\theta}{2}. $$ $\langle x_k^2\rangle_M$ について。$\ln Z_M$ を $\omega$ で微分します。定義の積分の側では、積分記号の下で微分して（。被積分関数の $\omega$ 微分は $\omega$ の有界な範囲で可積分な関数で抑えられます）$\frac{\partial}{\partial\omega}\ln Z_M = -\tau\omega\langle\sum_kx_k^2\rangle_M = -\beta\omega\langle x_1^2\rangle_M$ です（鎖の巡回対称性から $\langle x_k^2\rangle_M$ は $k$ によりません）。閉じた形の側では、$\frac{\partial\theta}{\partial\omega} = \tau/\cosh(\theta/2)$ から $\frac{\partial}{\partial\omega}\ln Z_M = -\frac M2\coth\frac{M\theta}{2}\cdot\frac{\tau}{\cosh(\theta/2)} = -\frac{\beta}{\omega}E_M$ です。比べて $\langle x_1^2\rangle_M = E_M/\omega^2$ です。 **3.** $u = \tau\omega = \beta\omega/M$ とおきます。$\theta = 2\,\mathrm{arsinh}(u/2) = u - u^3/24 + O(u^5)$、$\cosh(\theta/2) = \sqrt{1 + u^2/4} = 1 + u^2/8 + O(u^4)$ です。$M\theta/2 = x\,\theta/u = x(1 - u^2/24 + O(u^4))$ で、$\frac{d}{dx}\coth x = -1/\sinh^2x$ から $$ \coth\frac{M\theta}{2} = \coth x + \frac{x}{\sinh^2x}\cdot\frac{u^2}{24} + O(u^4) $$ です。掛け合わせて $$ E_M = \frac\omega2\Bigl(1 - \frac{u^2}{8}\Bigr)\Bigl(\coth x + \frac{xu^2}{24\sinh^2x}\Bigr) + O(u^4) = E - \frac\omega2u^2\Bigl[\frac{\coth x}{8} - \frac{x}{24\sinh^2x}\Bigr] + O(u^4) $$ です。$x \le \sinh x\cosh x$（$x \ge 0$）から $x/\sinh^2x \le \coth x$ なので、角括弧の中は $\coth x/8 - \coth x/24 = \coth x/12 > 0$ 以上です。 誤差は $(\beta\omega/M)^2 = (\tau\omega)^2$ に比例し、 の次数 $\tau^2$ と合っています。符号は負で、離散化した系のエネルギーは真の値より小さく出ます。低温（$x \to \infty$）では角括弧が $1/8$ に近づき、相対誤差は $(\tau\omega)^2/8$ です。高温（$x \to 0$）では角括弧が $1/(12x)$ に近づいて $E_M - E \approx -\frac1\beta\cdot\frac{(\beta\omega)^2}{12M^2}$ となります。どちらの場合も、必要な分割数を決めるのは $\tau\omega = \beta\omega/M$、つまり「1 刻みの間に系の最も速い振動がどれだけ進むか」です。$\beta\omega = 4$ で相対誤差を $10^{-3}$ にするには $M \approx 44$ が要ります（$x = 2$ で相対誤差は $0.119\,(\beta\omega/M)^2$）。 ## 5. Brown 橋と経路の生成 ### 5.1. 熱核の積の公式 $\pi_M$ から標本を作るとき、運動エネルギーの部分 $\prod_kp(\tau, x_{k+1} - x_k)$ は Gauss 型で、厳密に標本抽出できます。そのための計算が次の補題です。平均 $m$、分散 $v$ の正規分布の密度を $\phi(c; m, v) = (2\pi v)^{-1/2}\exp(-(c-m)^2/(2v))$ と書きます。 $\tau_1, \tau_2 > 0$、$a, b, c \in \mathbb{R}$ について $$ p(\tau_1, c - a)\,p(\tau_2, b - c) = p(\tau_1 + \tau_2, b - a)\;\phi(c; m, v), \qquad m = \frac{\tau_2a + \tau_1b}{\tau_1 + \tau_2}, \quad v = \frac{\tau_1\tau_2}{\tau_1 + \tau_2}. $$ 特に $\int_{\mathbb{R}}p(\tau_1, c - a)\,p(\tau_2, b - c)\,dc = p(\tau_1 + \tau_2, b - a)$（Chapman–Kolmogorov の等式）。 指数の中身について、$c$ の 2 次式として平方完成します。 $$ \frac{(c-a)^2}{\tau_1} + \frac{(b-c)^2}{\tau_2} = \frac{(c - m)^2}{v} + \frac{(b-a)^2}{\tau_1 + \tau_2} $$ を確かめます。$c^2$ の係数は左辺で $1/\tau_1 + 1/\tau_2 = (\tau_1+\tau_2)/(\tau_1\tau_2) = 1/v$。$c$ の係数は左辺で $-2(a/\tau_1 + b/\tau_2)$、右辺で $-2m/v = -2(\tau_2a + \tau_1b)/(\tau_1\tau_2)$ で一致します。定数項は、左辺の $a^2/\tau_1 + b^2/\tau_2 = (a^2\tau_2 + b^2\tau_1)(\tau_1+\tau_2)/(\tau_1\tau_2(\tau_1+\tau_2))$ から右辺の $m^2/v = (\tau_2a + \tau_1b)^2/(\tau_1\tau_2(\tau_1+\tau_2))$ を引くと、分子が $\tau_1\tau_2(a^2 + b^2 - 2ab)$ となり、$(b-a)^2/(\tau_1+\tau_2)$ です。係数は $(2\pi\tau_1)^{-1/2}(2\pi\tau_2)^{-1/2} = (2\pi(\tau_1+\tau_2))^{-1/2}(2\pi v)^{-1/2}$ です（$\tau_1\tau_2 = (\tau_1+\tau_2)v$）。後半は $\int\phi(c; m, v)\,dc = 1$ から従います。 Brown 橋は、直観的には「時刻 $t_1$ に $a$ から出発したブラウン運動の経路のうち、時刻 $t_2$ にちょうど $b$ に着くものだけを集めたもの」です。両端を釘で留めた揺れるひもを思い浮かべてください。次の定義は、それを有限個の時刻での密度で書いたものです。 $t_1 < t_2$、$a, b \in \mathbb{R}$ とする。確率過程 $(X_s)_{t_1 < s < t_2}$ が、任意の時刻 $t_1 < s_1 < \cdots < s_k < t_2$ について $(X_{s_1}, \ldots, X_{s_k})$ の密度が $$ q_{a\to b}^{[t_1, t_2]}(y_1, \ldots, y_k) = \frac{\prod_{i=0}^{k}p(s_{i+1} - s_i,\ y_{i+1} - y_i)}{p(t_2 - t_1,\ b - a)} \qquad (s_0 = t_1,\ y_0 = a,\ s_{k+1} = t_2,\ y_{k+1} = b) $$ であるとき、$X$ を $(t_1, a)$ から $(t_2, b)$ への **Brown 橋**という。 の後半を繰り返し使うと、$q$ を全変数で積分すると $1$ になり、1 つの変数について積分するとその時刻を除いた密度になる（時刻の取り方について矛盾がない）ことが分かります。$q$ は、$a$ から出発した Brown 運動（）の時刻 $s_1, \ldots, s_k, t_2$ での同時密度 $\prod_ip$ を、終点 $b$ の密度 $p(t_2 - t_1, b - a)$ で割ったもの、すなわち終点が $b$ であるという条件のもとでの条件付き密度です。$q$ は $y$ の 2 次式の指数関数なので、Brown 橋は Gauss 過程（）です。 環状の鎖との関係は次のとおりです。ビーズ $x_i$ と $x_{i+L}$ を固定すると、$e^{-S_M}$ のうちその間のビーズ $x_{i+1}, \ldots, x_{i+L-1}$ に依存する部分は $\prod_{k=i}^{i+L-1}p(\tau, x_{k+1} - x_k)\cdot\exp(-\tau\sum_{k=i+1}^{i+L-1}V(x_k))$ に比例します。したがって、$\pi_M$ のもとでの内側のビーズの条件付き密度は $$ q_{x_i\to x_{i+L}}^{[0, L\tau]}(x_{i+1}, \ldots, x_{i+L-1})\,\exp\Bigl(-\tau\sum_{k=i+1}^{i+L-1}V(x_k)\Bigr) $$ に比例します（時刻は $s_j = j\tau$）。自由粒子（$V = 0$）なら Brown 橋そのものです。 ### 5.2. 二分法 $X$ を $(t_1, a)$ から $(t_2, b)$ への Brown 橋、$T = t_2 - t_1$、$t_{\mathrm m} = (t_1 + t_2)/2$ とする。 1. $t_1 < s < t_2$ について、$X_s$ は平均 $a + \dfrac{s - t_1}{T}(b - a)$、分散 $\dfrac{(s - t_1)(t_2 - s)}{T}$ の正規分布に従う。特に中点 $X_{t_{\mathrm m}}$ は平均 $(a+b)/2$、分散 $T/4$ の正規分布に従う。 2. 時刻 $t_1 < s_1 < \cdots < s_k < t_{\mathrm m} < s_1' < \cdots < s_l' < t_2$ について、$(X_{s_1}, \ldots, X_{s_k}, X_{t_{\mathrm m}}, X_{s_1'}, \ldots, X_{s_l'})$ の $(y, c, y')$ での密度は $$ \phi\Bigl(c;\ \frac{a+b}{2},\ \frac T4\Bigr)\;q_{a\to c}^{[t_1, t_{\mathrm m}]}(y)\;q_{c\to b}^{[t_{\mathrm m}, t_2]}(y') $$ に等しい。すなわち、中点の値が $c$ であるという条件のもとで、左半分と右半分は独立で、それぞれ $(t_1, a)$ から $(t_{\mathrm m}, c)$ へ、$(t_{\mathrm m}, c)$ から $(t_2, b)$ への Brown 橋である。 **1.** で $k = 1$ とすると、$X_s$ の密度は $p(s - t_1, y - a)\,p(t_2 - s, b - y)/p(T, b - a)$ です。 を $\tau_1 = s - t_1$、$\tau_2 = t_2 - s$ で使うと、これは $\phi(y; m, v)$、$m = ((t_2 - s)a + (s - t_1)b)/T = a + (s - t_1)(b - a)/T$、$v = (s - t_1)(t_2 - s)/T$ です。$s = t_{\mathrm m}$ では $m = (a+b)/2$、$v = (T/2)^2/T = T/4$ です。 **2.** 定義から、密度は $$ \frac{\Pi_{\mathrm L}(y, c)\,\Pi_{\mathrm R}(c, y')}{p(T, b - a)} $$ です。ここで $\Pi_{\mathrm L}$ は $a$ から $y_1, \ldots, y_k$ を経て $c$ に至る $k + 1$ 個の熱核の積、$\Pi_{\mathrm R}$ は $c$ から $y_1', \ldots, y_l'$ を経て $b$ に至る $l + 1$ 個の熱核の積です。分子と分母に $p(T/2, c - a)\,p(T/2, b - c)$ を掛けて $$ \frac{\Pi_{\mathrm L}(y, c)}{p(T/2, c - a)}\cdot\frac{\Pi_{\mathrm R}(c, y')}{p(T/2, b - c)}\cdot\frac{p(T/2, c - a)\,p(T/2, b - c)}{p(T, b - a)} $$ と書くと、第 1 の因子は $q_{a\to c}^{[t_1, t_{\mathrm m}]}(y)$、第 2 の因子は $q_{c\to b}^{[t_{\mathrm m}, t_2]}(y')$ で、第 3 の因子は （$\tau_1 = \tau_2 = T/2$）から $\phi(c; (a+b)/2, T/4)$ です。$c$ を固定すると密度が $y$ の関数と $y'$ の関数の積なので、条件付きで独立です。 この定理を繰り返すと、両端を固定した区間の内点を、粗いほうから順に生成できます（P. Lévy による Brown 運動の構成と同じ手順です）。区間の長さを $L = 2^l$ 刻みとします。まず中央のビーズを両端の平均を中心とする分散 $L\tau/4$ の正規分布から選びます。次に、できた 2 つの半区間のそれぞれの中央を、その両端の平均を中心とする分散 $L\tau/8$ の正規分布から選びます。これを $l$ 段繰り返すと、$2^l - 1$ 個の内点がすべて決まり、その同時密度は の 2 から、ちょうど $q_{a\to b}^{[0, L\tau]}$ です。棄却はなく、手間は内点の個数に比例します。
段 1段 2段 2 段 3段 3段 3段 3 0Lτ/2Lτ
仮定 (C) のもとで、環状の鎖の区間 $x_i, \ldots, x_{i+L}$（$L = 2^l < M$）について次の更新を考える。両端 $x_i$、$x_{i+L}$ とその他のビーズを固定し、内点の候補 $y' = (x_{i+1}', \ldots, x_{i+L-1}')$ を二分法で（$q_{x_i\to x_{i+L}}^{[0, L\tau]}$ から）生成し、現在の内点 $y$ から確率 $$ \alpha(y, y') = \min\Bigl(1,\ \exp\Bigl(-\tau\sum_{k=i+1}^{i+L-1}\bigl(V(x_k') - V(x_k)\bigr)\Bigr)\Bigr) $$ で $y'$ に移り、そうでなければとどまる。この更新は $\pi_M$ を不変に保つ。 固定したビーズのもとでの内点の条件付き密度を $\pi^{\mathrm c}(y)$ と書くと、§5.1 から $\pi^{\mathrm c}(y) = c\,q(y)e^{-\tau\sum V(y_k)}$（$q = q_{x_i\to x_{i+L}}^{[0,L\tau]}$、$c$ は $y$ によらない定数）です。提案の密度は現在の $y$ によらず $q(y')$ です。Metropolis–Hastings 法（：密度 $\pi$、提案核 $g(y, y')$ に対して採択確率を $\min(1, \pi(y')g(y', y)/(\pi(y)g(y, y')))$ とする方法）の採択確率は $$ \min\Bigl(1,\ \frac{\pi^{\mathrm c}(y')\,q(y)}{\pi^{\mathrm c}(y)\,q(y')}\Bigr) = \min\Bigl(1,\ \frac{e^{-\tau\sum V(y_k')}}{e^{-\tau\sum V(y_k)}}\Bigr) = \alpha(y, y') $$ で、$q$ が打ち消し合ってポテンシャルだけが残ります。この遷移は $\pi^{\mathrm c}$ について詳細つり合いを満たし、$\pi^{\mathrm c}$ を不変に保ちます（）。固定したビーズの周辺分布は変わらず、内点の条件付き分布が保たれるので、同時分布 $\pi_M$ が保たれます。 運動エネルギーの部分（硬いばね）は提案が厳密に処理し、採択の判定にはなめらかで弱い外場 $\tau V$ しか現れません。$M$ を増やすとばねは硬くなりますが、採択率は $\tau\sum_kV$ の変化、つまり虚時間 $L\tau$ の区間でのポテンシャルの積分の変化で決まるので、$L\tau$ を固定すれば $M$ によらずほぼ一定です（§9 の実験では、区間の長さ $\beta/2$ で採択率が $M = 8$ から $64$ まで $0.79$〜$0.80$ です）。採択率が低いときは、粗い段で先に判定して見込みのない候補を早く捨てる多段の二分法が使われます。段階分割（staging）は、同じ Brown 橋を左から順に「次のビーズを、現在のビーズと終点で決まる正規分布から選ぶ」形で生成する方法で、 の 1 を $s$ が隣の時刻の場合に使います。 ## 6. エネルギーの推定量 内部エネルギー $E_M = -\partial_\beta\ln Z_M$ を、$\pi_M$ からの標本の平均として求めます。同じ期待値をもつ関数（推定量）はいくつもあり、分散がまったく違います。 $x \in \mathbb{R}^M$ について $$ E_{\mathrm T}(x) = \frac{M}{2\beta} - \frac{M}{2\beta^2}\sum_{k=1}^M(x_{k+1} - x_k)^2 + \frac1M\sum_{k=1}^MV(x_k), \qquad E_{\mathrm V}(x) = \frac1M\sum_{k=1}^M\Bigl[\frac12x_kV'(x_k) + V(x_k)\Bigr] $$ を、エネルギーの**熱力学的推定量**、**ビリアル推定量**という。 $E_{\mathrm T}$ は $\ln Z_M$ の被積分関数を $\beta$ で微分したものです。運動エネルギーに当たる部分は、$M/(2\beta)$ という大きな定数から、ほぼ同じ大きさの揺らぐ量を引いた形をしています。 1. 仮定 (C) のもとで $\langle E_{\mathrm T}\rangle_M = \langle E_{\mathrm V}\rangle_M = -\dfrac{\partial}{\partial\beta}\ln Z_M$。 2. $V(x) = \frac12\omega^2x^2$ とし、$\kappa_j = \dfrac4\tau\sin^2\dfrac{\pi j}{M}$（$j = 0, \ldots, M-1$）、$\nu = \tau\omega^2$ とおく。$\pi_M$ のもとでの分散は $$ \mathrm{Var}(E_{\mathrm T}) = \frac{1}{2\beta^2}\sum_{j=0}^{M-1}\Bigl(\frac{\kappa_j - \nu}{\kappa_j + \nu}\Bigr)^2, \qquad \mathrm{Var}(E_{\mathrm V}) = \frac{2\omega^4}{M^2}\sum_{j=0}^{M-1}\frac{1}{(\kappa_j + \nu)^2} $$ で、 $$ \frac{M - \beta^2\omega^2}{2\beta^2} \le \mathrm{Var}(E_{\mathrm T}) \le \frac{M}{2\beta^2}, \qquad \mathrm{Var}(E_{\mathrm V}) \le \frac{2}{\beta^2} + \frac{\beta^2\omega^4}{32} $$ が成り立つ。熱力学的推定量の分散は $M$ に比例して増え、ビリアル推定量の分散は $M$ について有界である。 証明の見通しです。1 の前半は、$\ln Z_M$ の微分が $E_{\mathrm T}$ の期待値になるという計算です。後半は、「全部のビーズを同じ割合で伸び縮みさせる」変数変換に対応する部分積分から、運動エネルギーの部分を $V'$ を含む式に書き換えて $E_{\mathrm V}$ を得ます（Euler の関係式は、2 次の斉次式 $K$ について $\sum x_k\partial_kK=2K$ となる性質のことです）。2 の分散は、鎖を揺れの独立な成分に分けて 1 成分ずつ計算します。 **1.** $w(x, \beta) = (2\pi\beta/M)^{-M/2}\exp\bigl(-\frac{M}{2\beta}\sum_k(x_{k+1}-x_k)^2 - \frac\beta M\sum_kV(x_k)\bigr)$ とおくと $Z_M = \int w\,dx$ で、 $$ \frac{\partial}{\partial\beta}\ln w = -\frac{M}{2\beta} + \frac{M}{2\beta^2}\sum_k(x_{k+1} - x_k)^2 - \frac1M\sum_kV(x_k) = -E_{\mathrm T}(x) $$ です。仮定 (C) から $|V(x)| \le |V(0)| + c_2(|x| + |x|^{r+1})$ で、$\beta$ が $(0,\infty)$ のコンパクトな区間を動くとき $|\partial_\beta w| = |E_{\mathrm T}|w$ は「多項式 $\times\ e^{-c|x|^2}$」の形の可積分な関数で一様に抑えられます。積分記号の下で微分して（）$\partial_\beta Z_M = -\int E_{\mathrm T}w\,dx$、すなわち $-\partial_\beta\ln Z_M = \langle E_{\mathrm T}\rangle_M$ です。 次に $\langle E_{\mathrm V}\rangle_M = \langle E_{\mathrm T}\rangle_M$ を部分積分で示します。各 $k$ について、他の変数を固定すると $x_ke^{-S_M(x)}$ は $|x_k| \to \infty$ で $0$ に収束し（$S_M \ge \tau c_0|x|^2 - \beta c_1$）、その $x_k$ についての導関数 $(1 - x_k\partial_kS_M)e^{-S_M}$ は $\mathbb{R}^M$ 上で可積分です。微分積分学の基本定理と Fubini の定理（）から $\int_{\mathbb{R}^M}\partial_k(x_ke^{-S_M})\,dx = 0$、すなわち $$ \langle x_k\,\partial_kS_M\rangle_M = 1 \qquad (k = 1, \ldots, M) $$ です。$S_M = K + \tau\sum_kV(x_k)$、$K(x) = \frac{1}{2\tau}\sum_k(x_{k+1} - x_k)^2$ と分けます。$K$ は 2 次の斉次式なので $\sum_kx_k\partial_kK = 2K$（Euler の関係式。$K(\lambda x) = \lambda^2K(x)$ を $\lambda$ で微分して $\lambda = 1$ とする）です。$k$ について足すと $$ M = \Bigl\langle\sum_kx_k\partial_kS_M\Bigr\rangle_M = \frac1\tau\Bigl\langle\sum_k(x_{k+1} - x_k)^2\Bigr\rangle_M + \tau\Bigl\langle\sum_kx_kV'(x_k)\Bigr\rangle_M $$ です。これを $E_{\mathrm T}$ の第 2 項に代入すると、$\tau = \beta/M$ から $$ \frac{M}{2\beta^2}\Bigl\langle\sum_k(x_{k+1}-x_k)^2\Bigr\rangle_M = \frac{M\tau}{2\beta^2}\Bigl(M - \tau\Bigl\langle\sum_kx_kV'(x_k)\Bigr\rangle_M\Bigr) = \frac{M}{2\beta} - \frac{1}{2M}\Bigl\langle\sum_kx_kV'(x_k)\Bigr\rangle_M $$ で、$\langle E_{\mathrm T}\rangle_M = \frac{1}{2M}\langle\sum_kx_kV'(x_k)\rangle_M + \frac1M\langle\sum_kV(x_k)\rangle_M = \langle E_{\mathrm V}\rangle_M$ です。 **2.** $S_M = \frac12x^{\mathsf T}Px$、$P = \frac1\tau\Lambda + \tau\omega^2I$ です。ここで $\Lambda$ は $\sum_k(x_{k+1} - x_k)^2 = x^{\mathsf T}\Lambda x$ となる実対称な巡回行列（対角成分 $2$、巡回的に隣り合う成分 $-1$）で、固有値は $4\sin^2(\pi j/M)$（$j = 0, \ldots, M-1$。固有ベクトルは $(e^{2\pi ijk/M})_k$）です。スペクトル定理（）から実直交行列 $U$ で $U^{\mathsf T}\Lambda U = \mathrm{diag}(4\sin^2(\pi j/M))$ となり、$P$ の固有値は $\kappa_j + \nu$ です。$\pi_M$ は平均 $0$、共分散 $P^{-1}$ の正規分布なので（ の 3）、$y = U^{\mathsf T}x$ は共分散 $\mathrm{diag}((\kappa_j + \nu)^{-1})$ の正規ベクトル（同じ定理の 4）で、成分 $y_j$ は独立です（）。$x^{\mathsf T}\Lambda x = \tau\sum_j\kappa_jy_j^2$、$|x|^2 = \sum_jy_j^2$ なので $$ E_{\mathrm T} = \frac{M}{2\beta} + \sum_j\Bigl(-\frac{M\tau\kappa_j}{2\beta^2} + \frac{\omega^2}{2M}\Bigr)y_j^2 = \frac{M}{2\beta} + \sum_j\frac{\nu - \kappa_j}{2\beta}\,y_j^2, \qquad E_{\mathrm V} = \frac{\omega^2}{M}\sum_jy_j^2 $$ です（$M\tau = \beta$、$\omega^2/M = \nu/\beta$）。平均 $0$、分散 $\sigma^2$ の正規確率変数 $y$ について $\mathbb{E}[y^4] = 3\sigma^4$（）なので $\mathrm{Var}(y^2) = 2\sigma^4$ です。独立な項の分散を足して $$ \mathrm{Var}(E_{\mathrm T}) = \sum_j\frac{(\nu - \kappa_j)^2}{4\beta^2}\cdot\frac{2}{(\kappa_j+\nu)^2}, \qquad \mathrm{Var}(E_{\mathrm V}) = \frac{\omega^4}{M^2}\sum_j\frac{2}{(\kappa_j + \nu)^2} $$ で、主張の式です。 不等式を示します。$\kappa_j \ge 0$、$\nu > 0$ から各項 $((\kappa_j - \nu)/(\kappa_j+\nu))^2 \le 1$ で、上界 $M/(2\beta^2)$ が出ます。下界について、$1 - ((\kappa - \nu)/(\kappa+\nu))^2 = 4\kappa\nu/(\kappa+\nu)^2 \le 4\nu/\kappa$（$\kappa > 0$）です。$1 \le j \le M - 1$ について $j' = \min(j, M - j) \le M/2$ とおくと、$\sin y \ge 2y/\pi$（$0 \le y \le \pi/2$。$\sin$ がこの区間で上に凸で、両端で等号）から $\sin(\pi j/M) = \sin(\pi j'/M) \ge 2j'/M$ で、 $$ \kappa_j \ge \frac{16j'^2}{\tau M^2}, \qquad \frac{4\nu}{\kappa_j} \le \frac{\tau^2M^2\omega^2}{4j'^2} = \frac{\beta^2\omega^2}{4j'^2} $$ です。各 $j'$ は高々 2 つの $j$ に対応し、$\sum_{i\ge1}i^{-2} \le 1 + \sum_{i\ge2}\frac{1}{i(i-1)} = 2$ なので、$\sum_{j=1}^{M-1}4\nu/\kappa_j \le 2\cdot\frac{\beta^2\omega^2}{4}\cdot2 = \beta^2\omega^2$ です。$j = 0$ の項は $\kappa_0 = 0$ から $1$ です。よって $\sum_j((\kappa_j-\nu)/(\kappa_j+\nu))^2 \ge M - \beta^2\omega^2$ です。 $E_{\mathrm V}$ について。$j = 0$ の項は $\frac{2\omega^4}{M^2\nu^2} = \frac{2}{M^2\tau^2} = \frac{2}{\beta^2}$ です。$j \ge 1$ では $(\kappa_j + \nu)^{-2} \le \kappa_j^{-2} \le \tau^2M^4/(256j'^4)$ で、$\sum_{i\ge1}i^{-4} \le 2$ から $\sum_{j=1}^{M-1}(\kappa_j+\nu)^{-2} \le 2\cdot\frac{\tau^2M^4}{256}\cdot2 = \frac{\tau^2M^4}{64}$、これに $2\omega^4/M^2$ を掛けて $\omega^4\tau^2M^2/32 = \beta^2\omega^4/32$ です。 $N$ 個の独立な標本で $E_M$ を求めるときの統計誤差は、熱力学的推定量では $\sqrt{M/(2\beta^2N)}$ 程度です。離散化誤差を半分にするために $M$ を $\sqrt2$ 倍にすると、統計誤差が $2^{1/4}$ 倍になります。ビリアル推定量にはこの問題がありません。証明の 1 が示すとおり、2 つの推定量の差は「全体を一様に伸縮する変数変換」についての部分積分（恒等式 $\langle x_k\partial_kS_M\rangle = 1$）で、期待値には寄与せず分散だけを変えます。同じ期待値をもつ推定量の中から分散の小さいものを選ぶという点で、[制御変量](/computer-science/stochastic-computation/monte-carlo-integration)と同じ考え方です。 $\mathrm{Var}(E_{\mathrm V})$ の上界の第 1 項 $2/\beta^2$ は、鎖の重心（$j = 0$ の成分）の揺らぎです。重心からのずれ $x_k - \bar x$ を使う重心ビリアル推定量 $\frac{1}{2\beta} + \frac{1}{2M}\sum_k(x_k - \bar x)V'(x_k) + \frac1M\sum_kV(x_k)$ は、高温でこの項を取り除きます。導出は 1 と同じで、伸縮の中心を重心に取ります。 ## 7. 局所更新の緩和 ビーズを 1 個ずつ動かす更新は実装が簡単ですが、$M$ が大きいと緩和が遅くなります。運動エネルギーの部分だけを取り出した模型で、緩和の速さを厳密に計算します。 熱浴更新とは、1 個の変数を選び、他の変数を固定したときのその変数の条件付き分布から新しい値を選び直す更新のことです。棄却がなく、選び直す相手が常に正しい分布になっています。次の命題は、これを 1 点ずつ繰り返しても鎖全体のゆるやかなうねりが変わるのに時間がかかる、ということを固有値で示します。 $M \ge 2$ とし、両端 $x_0 = x_M = 0$ を固定した自由な鎖の内点 $x = (x_1, \ldots, x_{M-1})$ の分布 $$ \pi(x) \propto \exp\Bigl(-\frac{1}{2\tau}\sum_{k=0}^{M-1}(x_{k+1} - x_k)^2\Bigr) $$ を考える。$k \in \{1, \ldots, M-1\}$ を一様に選び、$x_k$ を平均 $(x_{k-1} + x_{k+1})/2$、分散 $\tau/2$ の正規分布から選び直す更新（無作為走査の熱浴更新）の遷移作用素を $P$ とする。 1. $\pi$ は $P$ の不変分布である。 2. $v_j = (\sin(jk\pi/M))_{k=1}^{M-1}$、$f_j(x) = \langle v_j, x\rangle$（$j = 1, \ldots, M-1$）について $$ Pf_j = \mu_jf_j, \qquad \mu_j = 1 - \frac{1 - \cos(j\pi/M)}{M - 1}. $$ 3. 定常な連鎖（$X_0 \sim \pi$）で $A_t = f_1(X_t)$ とすると、$\rho(t) = \mu_1^t$ で、積分自己相関時間は $$ \tau_{\mathrm{int}} = \frac{1 + \mu_1}{2(1 - \mu_1)} \ge \frac{2M^2(M-1)}{\pi^2} - \frac12 $$ である（単位は 1 点の更新 1 回）。$M - 1$ 回の更新を 1 掃引とすると、$\tau_{\mathrm{int}}$ は $2M^2/\pi^2 - 1/2$ 掃引以上である。 **1.** 指数の中で $x_k$ を含む項は $$ \frac{(x_k - x_{k-1})^2 + (x_{k+1} - x_k)^2}{2\tau} = \frac1\tau\Bigl(x_k - \frac{x_{k-1} + x_{k+1}}{2}\Bigr)^2 + (x_k\text{ によらない項}) $$ なので、他の成分を固定したときの $x_k$ の条件付き分布は、平均 $(x_{k-1}+x_{k+1})/2$、分散 $\tau/2$ の正規分布です。$X \sim \pi$ のとき、$X_k$ をその条件付き分布から選び直したベクトルの密度は「他の成分の周辺密度 $\times$ 条件付き密度」で、$\pi$ に等しくなります。$k$ を無作為に選んでも同じです。 **2.** $T$ を、対角成分 $2$、すぐ上とすぐ下の成分 $-1$ の $M - 1$ 次の三重対角行列とします。$x_0 = x_M = 0$ なので $(Tx)_k = 2x_k - x_{k-1} - x_{k+1}$ で、$k$ を更新したあとの $x_k'$ の条件付き期待値は $(x_{k-1} + x_{k+1})/2 = x_k - (Tx)_k/2$ です。$f(x) = \langle v, x\rangle$ について、$k$ が選ばれる確率が $1/(M-1)$ で、そのとき変わるのは第 $k$ 成分だけなので $$ (Pf)(x) = \mathbb{E}[f(x')\mid x] = f(x) - \frac{1}{M-1}\sum_kv_k\frac{(Tx)_k}{2} = f(x) - \frac{\langle Tv, x\rangle}{2(M-1)} $$ です（$T$ は対称）。三重対角 Toeplitz 行列の固有値（ を $n = M - 1$、$\alpha = 2$、$\beta = -1$ で使う）から $Tv_j = (2 - 2\cos(j\pi/M))v_j$ なので、$Pf_j = (1 - (1 - \cos(j\pi/M))/(M-1))f_j$ です。 **3.** $\pi$ は平均 $0$ の正規分布なので $\mathbb{E}_\pi[f_1] = 0$、$0 < \mathbb{E}_\pi[f_1^2] < \infty$ です。Markov 性から $\mathbb{E}[f_1(X_t)\mid X_0] = (P^tf_1)(X_0) = \mu_1^tf_1(X_0)$ なので $$ \mathrm{Cov}(A_0, A_t) = \mathbb{E}\bigl[f_1(X_0)\,\mathbb{E}[f_1(X_t)\mid X_0]\bigr] = \mu_1^t\,\mathbb{E}_\pi[f_1^2] $$ で、$\rho(t) = \mu_1^t$ です。$0 < \mu_1 < 1$ なので $\tau_{\mathrm{int}} = \frac12 + \sum_{t\ge1}\mu_1^t = \frac{1 + \mu_1}{2(1 - \mu_1)} = \frac{1}{1 - \mu_1} - \frac12$ です。$1 - \cos y \le y^2/2$ から $1 - \mu_1 \le \pi^2/(2M^2(M-1))$ で、下界を得ます。 最も遅いのは、鎖全体がゆっくりうねる波長の長い変形（$j = 1$）です。1 点の更新は隣の平均のまわりに $\sqrt{\tau/2}$ しか動けず、長い波長の変形は拡散的に、$M^2$ 掃引かけて緩和します。1 掃引の手間が $M$ に比例するので、独立な標本 1 個あたりの手間は $M^3$ に比例します。これは格子場の理論で連続極限に近づくときの臨界減速と同じ現象です。 二分法や段階分割がこれを避ける理由は、 から明らかです。この命題の模型（自由な鎖）では、二分法で内点を全部生成し直した結果は直前の配位と独立で、$\tau_{\mathrm{int}} = 1/2$、手間は $M$ に比例します。ポテンシャルがあるときも、長さ $L\tau$ の区間の中の長波長の変形が 1 回の更新で入れ替わるので、緩和に要る掃引数は $M$ ではなく、採択率を保てる $L\tau$ と $\beta$ の比で決まります。 ## 8. 符号問題 ### 8.1. 符号つきの重みと平均符号 ここまでは、経路の重み $e^{-S_M}$ が正でした。重みが負（または複素数）になる系では、重みを確率と見なせません。同種の Fermi 粒子がその代表です。$N_{\mathrm p}$ 個の Fermi 粒子では、虚時間の終わりに粒子の置換 $\sigma$ を許し、奇置換の経路に符号 $-1$ を付けて足します。フラストレーションのあるスピン系、実時間の経路積分（重みが $e^{iS}$）でも同じことが起きます。 $(\Omega, \mu)$ を測度空間、$w : \Omega \to \mathbb{R}$ を可測関数で $0 < Z_{|w|} = \int|w|\,d\mu < \infty$、$Z = \int w\,d\mu \ne 0$ とする。$O\,w$ が可積分な関数 $O$ について $$ \langle O\rangle_w = \frac1Z\int O\,w\,d\mu $$ と書く。$p = |w|/Z_{|w|}$ は $\mu$ についての確率密度で、$p$ についての期待値を $\langle\cdot\rangle_{|w|}$ と書く。$s = \mathrm{sgn}\,w$（$w$ の符号）について $\langle s\rangle_{|w|}$ を**平均符号**という。 の設定で次が成り立つ。 1. $\langle s\rangle_{|w|} = \dfrac{Z}{Z_{|w|}}$、$|\langle s\rangle_{|w|}| \le 1$、 $$ \langle O\rangle_w = \frac{\langle Os\rangle_{|w|}}{\langle s\rangle_{|w|}}. $$ $Z > 0$ のとき $Z = e^{-\beta F}$、$Z_{|w|} = e^{-\beta F_{|w|}}$ で $F$、$F_{|w|}$ を定めると、$\langle s\rangle_{|w|} = e^{-\beta\Delta F}$、$\Delta F = F - F_{|w|} \ge 0$。 2. $X_1, X_2, \ldots$ を密度 $p$ に従う独立な確率変数とし、$O^2|w|$ が可積分とする。 $$ \hat R_N = \frac{\sum_{i=1}^NO(X_i)s(X_i)}{\sum_{i=1}^Ns(X_i)} $$ （分母が $0$ のときは $0$ とする）は確率 $1$ で $\langle O\rangle_w$ に収束し、$\sqrt N(\hat R_N - \langle O\rangle_w)$ は平均 $0$、分散 $$ \sigma_R^2 = \frac{\langle(O - \langle O\rangle_w)^2\rangle_{|w|}}{\langle s\rangle_{|w|}^2} = e^{2\beta\Delta F}\,\langle(O - \langle O\rangle_w)^2\rangle_{|w|} $$ の正規分布に分布収束する（最後の等号は $Z > 0$ のとき）。 **1.** $w = s|w|$ なので $\int Ow\,d\mu = Z_{|w|}\int Os\,p\,d\mu = Z_{|w|}\langle Os\rangle_{|w|}$ です。$O = 1$ とすると $Z = Z_{|w|}\langle s\rangle_{|w|}$ で、割り算して $\langle O\rangle_w = \langle Os\rangle_{|w|}/\langle s\rangle_{|w|}$ です。$|Z| = |\int w\,d\mu| \le \int|w|\,d\mu = Z_{|w|}$ から $|\langle s\rangle_{|w|}| \le 1$ で、$Z > 0$ のとき $e^{-\beta\Delta F} = Z/Z_{|w|} \le 1$ から $\Delta F \ge 0$ です。 **2.** $r = \langle O\rangle_w$、$Y_i = (O(X_i) - r)s(X_i)$ とおきます。$p$ のもとで $w \ne 0$ が確率 $1$ で成り立つので $s(X_i)^2 = 1$ で、 $$ \mathbb{E}[Y_i] = \langle Os\rangle_{|w|} - r\langle s\rangle_{|w|} = 0, \qquad \mathbb{E}[Y_i^2] = \langle(O - r)^2\rangle_{|w|} < \infty $$ です（1 と、$O^2|w|$ の可積分性）。$\bar s_N = N^{-1}\sum_is(X_i)$、$\bar Y_N = N^{-1}\sum_iY_i$ とすると、$\bar s_N \ne 0$ のとき $$ \hat R_N - r = \frac{\sum_i(O(X_i) - r)s(X_i)}{\sum_is(X_i)} = \frac{\bar Y_N}{\bar s_N} $$ です。大数の強法則（）から確率 $1$ で $\bar s_N \to \langle s\rangle_{|w|} \ne 0$、$\bar Y_N \to 0$ なので、$\hat R_N \to r$ です。中心極限定理（）から $\sqrt N\bar Y_N$ は $N(0, \mathbb{E}[Y_1^2])$ に分布収束し（$\mathbb{E}[Y_1^2] = 0$ のときは $Y_i = 0$ で自明）、$\bar s_N$ は定数 $\langle s\rangle_{|w|}$ に確率収束するので、Slutsky の定理（）から $\sqrt N(\hat R_N - r) = \sqrt N\bar Y_N/\bar s_N$ は $N(0, \mathbb{E}[Y_1^2]/\langle s\rangle_{|w|}^2)$ に分布収束します（$\bar s_N = 0$ となる事象の確率は $0$ に向かうので、極限に影響しません）。 2 の分散は、比の推定量に 2 変数の delta 法（ の多変数版）を使って得られる式と同じです。標準誤差を $\varepsilon$ にするのに要る標本数は $N \approx \sigma_R^2/\varepsilon^2$ で、$e^{2\beta\Delta F}$ に比例します。平均符号そのものについても、$\mathrm{Var}(s) = 1 - \langle s\rangle^2$ なので相対誤差は $\sqrt{(1 - \langle s\rangle^2)/N}/\langle s\rangle \approx e^{\beta\Delta F}/\sqrt N$ です。 $\Delta F$ は、符号つきの系（Fermi 粒子）と、符号を無視した系（多くの場合は同じハミルトニアンの Bose 粒子）の自由エネルギーの差です。2 つの系がどちらも熱力学極限をもち、自由エネルギーが示量的（粒子数にほぼ比例する）で、1 粒子あたりの値に差があるとき、$\Delta F$ は粒子数 $N_{\mathrm p}$ にほぼ比例します（これは個々の系について確かめるべき仮定で、 が証明しているのは $\Delta F$ を与えたときの関係だけです）。そのとき必要な標本数は $e^{2\beta N_{\mathrm p}\Delta f}$（$\Delta f$ は 1 粒子あたりの差）で、温度を下げても粒子を増やしても指数的に増えます。Monte Carlo 法の誤差が次元によらず $N^{-1/2}$ であるという利点が、ここで失われます。失われるのは収束の次数ではなく、$N^{-1/2}$ の前の係数です。 振動数 $\omega$ の 1 次元調和振動子の中の、スピンのない同種粒子 2 個を考えます。1 粒子の分配関数は $Z_1(\beta) = 1/(2\sinh(\beta\omega/2)) = q^{1/2}/(1-q)$（$q = e^{-\beta\omega}$）です。経路は、置換が恒等のとき 2 本の別々の閉じた経路（重みの和 $Z_1(\beta)^2/2$）、互換のとき長さ $2\beta$ の 1 本の閉じた経路（重みの和 $Z_1(2\beta)/2$）です。Bose 粒子と Fermi 粒子の分配関数は $$ Z_{\mathrm B} = \frac{Z_1(\beta)^2 + Z_1(2\beta)}{2}, \qquad Z_{\mathrm F} = \frac{Z_1(\beta)^2 - Z_1(2\beta)}{2} $$ で、Fermi 粒子の計算で符号を無視した系が Bose 粒子です。$Z_1(2\beta)/Z_1(\beta)^2 = (1-q)^2/(1-q^2) = (1-q)/(1+q)$ から $$ \langle s\rangle = \frac{Z_{\mathrm F}}{Z_{\mathrm B}} = \frac{1 - \frac{1-q}{1+q}}{1 + \frac{1-q}{1+q}} = q = e^{-\beta\omega} $$ で、$\Delta F = \omega$ がすべての温度で厳密に成り立ちます。基底エネルギーの差（Fermi 粒子は $\omega/2 + 3\omega/2 = 2\omega$、Bose 粒子は $\omega$）に等しい値です。$\beta\omega = 8$ では $\langle s\rangle = 3.4\times10^{-4}$ で、恒等置換と互換が $0.50017 : 0.49983$ の割合で現れ、その差の中に答えがあります。§9 でこの模型を数値的に調べます。 符号問題は基底の取り方に依存します。$H$ を対角化する基底では重みはすべて正ですが、その基底を見つけることが元の問題です。符号問題を一般に（多項式時間で）解く方法があれば、計算量の理論で NP と呼ばれるクラスの問題がすべて確率的な多項式時間で解けることになる、という結果があります（M. Troyer, U.-J. Wiese, *Phys. Rev. Lett.* 94 (2005) 170201）。したがって一般的な解決は期待されておらず、個々の系の構造を使った方法（節を固定する近似、符号のない基底への変換、複素化した積分路の変形）が研究されています。Bose 粒子の系には符号問題がなく、置換を効率よく標本にするワーム・アルゴリズム（M. Boninsegni, N. V. Prokof'ev, B. V. Svistunov, *Phys. Rev. Lett.* 96 (2006) 070601）によって、数千粒子の超流動ヘリウムの計算ができます。経路積分 Monte Carlo 法の全体像は D. M. Ceperley の総説（*Rev. Mod. Phys.* 67 (1995)）にあります。 ## 9. 数値実験 コードは Python 3.9.6・NumPy 2.0.2・SciPy 1.13.1・mpmath 1.3.0 で実行し、出力をそのまま載せています。乱数は `numpy.random.default_rng(20261008)` です。 ### 9.1. 分割の誤差の次数 $\beta = 1$ とし、2 組の半正定値エルミート行列で と を確かめます。1 つは $2\times2$ の $A = \mathrm{diag}(1, 0)$、$B = \frac12\begin{pmatrix}1&-1\\-1&1\end{pmatrix}$。もう 1 つは、調和振動子の固有基底を 20 次で切断した非調和振動子で、$A = p^2/2$、$B = x^4/10$（$x$、$p$ は切断した行列で、$x^4$ は行列 $x^2$ の 2 乗）です。参照値 $e^{-\beta H}$ は `scipy.linalg.expm` で計算し、`mpmath.expm`（30 桁）と比べます。 ```python import numpy as np from scipy.linalg import expm import mpmath as mp np.seterr(divide="ignore", over="ignore", invalid="ignore") # macOS の Accelerate BLAS が matmul で出す偽の警告を抑える def comm(x, y): return x @ y - y @ x def study(name, A, B, beta=1.0): n = len(A) H = A + B ref = expm(-beta * H) mp.mp.dps = 30 ref_mp = np.array(mp.expm(mp.matrix((-beta * H).tolist())).tolist(), dtype=float) nrm = lambda x: np.linalg.norm(x, 2) c2 = nrm(comm(A, comm(A, B))) / 24 + nrm(comm(B, comm(B, A))) / 12 c1 = nrm(comm(A, B)) / 2 print(f"{name}: n = {n}, |expm(scipy) - expm(mpmath)| = {nrm(ref - ref_mp):.1e}, " f"c_sym = {c2:.4g}, c_lie = {c1:.4g}") print(" M err(sym) bound err(Lie) bound |Z_M - Z|") rows = [] for k in range(9): M = 2**k t = beta / M ea, eb, eh = expm(-t * A), expm(-t * B), expm(-t * A / 2) sym = np.linalg.matrix_power(eh @ eb @ eh, M) lie = np.linalg.matrix_power(ea @ eb, M) e = (nrm(sym - ref), nrm(lie - ref), abs(np.trace(sym) - np.trace(ref)), abs(np.trace(lie) - np.trace(sym))) rows.append(e) print(f" {M:4d} {e[0]:.3e} {beta**3 * c2 / M**2:.3e} {e[1]:.3e} {beta**2 * c1 / M:.3e} {e[2]:.3e}") r = np.array(rows) ks = np.arange(4, 9) sl = [np.polyfit(ks, np.log2(r[4:, j]), 1)[0] for j in range(3)] print(f" slope (M >= 16): sym {sl[0]:.2f}, Lie {sl[1]:.2f}, Z {sl[2]:.2f}; max |Tr Lie^M - Tr sym^M| = {r[:, 3].max():.1e}") A2 = np.array([[1.0, 0.0], [0.0, 0.0]]) B2 = 0.5 * np.array([[1.0, -1.0], [-1.0, 1.0]]) study("2x2", A2, B2) n = 20 a = np.diag(np.sqrt(np.arange(1, n)), 1) # 消滅演算子（調和振動子の基底で 20 次に切断） x = (a + a.T) / np.sqrt(2) p2 = -((a.T - a) @ (a.T - a)) / 2 # p^2, p = i (a^T - a)/sqrt(2) x2 = x @ x study("20x20 (p^2/2 + x^4/10)", p2 / 2, x2 @ x2 / 10) ``` ```text 2x2: n = 2, |expm(scipy) - expm(mpmath)| = 6.7e-17, c_sym = 0.0625, c_lie = 0.25 M err(sym) bound err(Lie) bound |Z_M - Z| 1 2.235e-02 6.250e-02 1.233e-01 2.500e-01 8.055e-03 2 5.726e-03 1.562e-02 5.601e-02 1.250e-01 2.063e-03 4 1.440e-03 3.906e-03 2.649e-02 6.250e-02 5.189e-04 8 3.605e-04 9.766e-04 1.286e-02 3.125e-02 1.299e-04 16 9.016e-05 2.441e-04 6.335e-03 1.562e-02 3.249e-05 32 2.254e-05 6.104e-05 3.144e-03 7.812e-03 8.124e-06 64 5.635e-06 1.526e-05 1.566e-03 3.906e-03 2.031e-06 128 1.409e-06 3.815e-06 7.814e-04 1.953e-03 5.078e-07 256 3.522e-07 9.537e-07 3.903e-04 9.766e-04 1.269e-07 slope (M >= 16): sym -2.00, Lie -1.00, Z -2.00; max |Tr Lie^M - Tr sym^M| = 7.5e-15 20x20 (p^2/2 + x^4/10): n = 20, |expm(scipy) - expm(mpmath)| = 2.7e-15, c_sym = 1461, c_lie = 107.5 M err(sym) bound err(Lie) bound |Z_M - Z| 1 5.278e-02 1.461e+03 1.253e-01 1.075e+02 8.766e-02 2 1.612e-02 3.652e+02 5.817e-02 5.376e+01 2.048e-02 4 4.344e-03 9.130e+01 2.755e-02 2.688e+01 5.041e-03 8 1.110e-03 2.282e+01 1.334e-02 1.344e+01 1.269e-03 16 2.792e-04 5.706e+00 6.556e-03 6.720e+00 3.183e-04 32 6.991e-05 1.427e+00 3.250e-03 3.360e+00 7.966e-05 64 1.748e-05 3.566e-01 1.618e-03 1.680e+00 1.992e-05 128 4.372e-06 8.916e-02 8.071e-04 8.400e-01 4.980e-06 256 1.093e-06 2.229e-02 4.031e-04 4.200e-01 1.245e-06 slope (M >= 16): sym -2.00, Lie -1.01, Z -2.00; max |Tr Lie^M - Tr sym^M| = 3.7e-14 ``` - 2 通りの参照値の差は丸めの大きさ（$10^{-15}$ 以下）で、参照値は信頼できます。 - 誤差の傾きは、対称分割で $-2.00$、非対称分割で $-1.00$、分配関数で $-2.00$ です。対称分割と非対称分割のトレースの差は $10^{-14}$ 以下で、 の 1 のとおり丸めを除いて一致します。 - $2\times2$ の例では、誤差は上界の $0.36$〜$0.37$ 倍（対称）と $0.4$〜$0.5$ 倍（非対称）で、上界は定数まで含めて良い評価です（定数 $c_{\mathrm{sym}} = 1/16$、$c_{\mathrm{Lie}} = 1/4$ の手計算は ）。 - $20\times20$ の例では、上界は成り立っていますが、実際の誤差の $2\times10^4$〜$3\times10^4$ 倍（対称）です。交換子のノルムは切断した空間の高いエネルギーの状態で決まり、$e^{-\beta H}$ の主な部分である低いエネルギーの状態での誤差を表していません。切断の次元を上げるとノルムはさらに大きくなり、連続空間では無限大です。 の (c) で述べた、非有界な作用素では状態の滑らかさを使った評価が要るという事情が、有限次元でもすでに見えています。次数 $-2$ のほうは切断によらず成り立っています。 ### 9.2. 調和振動子の経路積分 Monte Carlo $\omega = 1$、$\beta = 4$ の調和振動子を、二分法の更新（。区間の長さは $M/2$ 刻みで、始点を無作為に選ぶ）で計算します。独立な 4000 本の鎖を同時に走らせ、各鎖で 300 回の更新を捨てたあと 3000 回の更新について 2 つの推定量を平均し、鎖ごとの平均のばらつきから標準誤差を求めます（鎖が独立なので、自己相関の解析なしで標準誤差が出ます）。 ```python import numpy as np beta, omega = 4.0, 1.0 V = lambda x: 0.5 * omega**2 * x * x def exact(M): # 離散化した系の厳密な値: E_M, Var(E_T), Var(E_V) tau = beta / M theta = 2 * np.arcsinh(tau * omega / 2) e_m = omega / (2 * np.cosh(theta / 2)) / np.tanh(M * theta / 2) ell = 4 * np.sin(np.pi * np.arange(M) / M) ** 2 kap, nu = ell / tau, tau * omega**2 var_t = ((kap - nu) / (kap + nu)) ** 2 return e_m, var_t.sum() / (2 * beta**2), 2 * omega**4 / M**2 * (1.0 / (kap + nu) ** 2).sum() def bisect(rng, left, right, tau, levels): # 両端 left, right を固定した 2^levels - 1 個の内点を Brown 橋から生成する n = 2**levels path = np.empty((left.shape[0], n + 1)); path[:, 0], path[:, n] = left, right step = n while step > 1: half = step // 2 for i in range(0, n, step): mean = 0.5 * (path[:, i] + path[:, i + step]) path[:, i + half] = mean + np.sqrt(step * tau / 4) * rng.standard_normal(left.shape[0]) step = half return path[:, 1:n] rng = np.random.default_rng(20261008) walkers, burn, meas = 4000, 300, 3000 print(f"beta = {beta}, exact E = {omega / 2 / np.tanh(beta * omega / 2):.6f}") print(" M E_M(exact) E_T (se) E_V (se) Var E_T (exact) Var E_V (exact) acc") for M in (8, 16, 32, 64): tau, levels = beta / M, int(np.log2(M)) - 1 # 区間の長さ M/2 seg = 2**levels x = np.zeros((walkers, M)) rows = np.arange(walkers)[:, None] sum_t = np.zeros(walkers); sum_v = np.zeros(walkers); sq_t = sq_v = 0.0; acc = 0.0 for it in range(burn + meas): start = rng.integers(0, M, size=walkers) idx = (start[:, None] + np.arange(1, seg)[None, :]) % M new = bisect(rng, x[np.arange(walkers), start], x[np.arange(walkers), (start + seg) % M], tau, levels) old = x[rows, idx] ok = rng.random(walkers) < np.exp(-tau * (V(new) - V(old)).sum(axis=1)) x[rows[ok], idx[ok]] = new[ok] if it >= burn: d = np.roll(x, -1, axis=1) - x e_t = M / (2 * beta) - M / (2 * beta**2) * (d * d).sum(axis=1) + V(x).mean(axis=1) e_v = (0.5 * x * omega**2 * x + V(x)).mean(axis=1) sum_t += e_t; sum_v += e_v; sq_t += (e_t**2).mean(); sq_v += (e_v**2).mean(); acc += ok.mean() mt, mv = sum_t / meas, sum_v / meas se = lambda m: m.std(ddof=1) / np.sqrt(walkers) e_m, vt, vv = exact(M) print(f" {M:3d} {e_m:.6f} {mt.mean():.5f} ({se(mt):.5f}) {mv.mean():.5f} ({se(mv):.5f}) " f"{sq_t / meas - mt.mean()**2:7.3f} ({vt:7.3f}) {sq_v / meas - mv.mean()**2:.4f} ({vv:.4f}) {acc / meas:.2f}") ``` ```text beta = 4.0, exact E = 0.518657 M E_M(exact) E_T (se) E_V (se) Var E_T (exact) Var E_V (exact) acc 8 0.503935 0.50416 (0.00027) 0.50407 (0.00041) 0.150 ( 0.150) 0.1521 (0.1521) 0.80 16 0.514848 0.51442 (0.00040) 0.51421 (0.00040) 0.392 ( 0.392) 0.1487 (0.1496) 0.79 32 0.517697 0.51809 (0.00056) 0.51793 (0.00039) 0.890 ( 0.890) 0.1486 (0.1489) 0.79 64 0.518417 0.51824 (0.00081) 0.51856 (0.00040) 1.891 ( 1.890) 0.1491 (0.1487) 0.79 ``` - **離散化誤差**。2 つの推定量の平均は、連続極限の値 $0.518657$ ではなく、離散化した系の値 $E_M$（ の 2）に標準誤差の範囲で一致します。$E_M$ との差は標準誤差の $+0.8, -1.1, +0.7, -0.2$ 倍（$E_{\mathrm T}$）、$+0.3, -1.6, +0.6, +0.4$ 倍（$E_{\mathrm V}$）です。$M = 8$ では真の値との差 $-0.0147$ が標準誤差の 50 倍あり、標本を増やしても消えない偏りです。先頭項の式（ の 3）の値は $M = 8$ で $-0.0154$、$M = 32$ で $-0.00096$ で、厳密な差 $-0.0147$、$-0.00096$ と合います。 - **推定量の分散**。$\mathrm{Var}(E_{\mathrm T})$ は $0.150, 0.392, 0.890, 1.891$ と $M$ にほぼ比例して増え、 の閉じた形と 3 桁一致します（$M = 64$ で上界 $M/(2\beta^2) = 2$、下界 $1.5$）。$\mathrm{Var}(E_{\mathrm V})$ は $0.149$ 前後で変わりません。標準誤差も、$E_{\mathrm T}$ では $M$ とともに増え、$E_{\mathrm V}$ では一定です。 - **採択率**は $M$ によらず約 $0.79$ です。更新する虚時間の長さ $\beta/2$ を固定しているためです。 ### 9.3. 符号問題の模型 の 2 個の Fermi 粒子（$\omega = 1$）で、置換の種類（恒等か互換か）だけを標本にします。恒等は確率 $(1+q)/2$ で符号 $+1$、互換は確率 $(1-q)/2$ で符号 $-1$ です。観測量 $O$ は、それぞれの重みの $-\partial_\beta\ln$、すなわち $2E_1(\beta)$ と $2E_1(2\beta)$（$E_1(\beta) = \frac12\coth\frac\beta2$）とします。このとき $\langle Os\rangle/\langle s\rangle = -\partial_\beta\ln Z_{\mathrm F} = 2 + \frac{2q}{1-q} - \frac{q}{1+q}$ が厳密に成り立ちます。経路の揺らぎを省いた、符号だけの模型です。 ```python import numpy as np # 調和振動子（omega = 1）の中の同種 Fermi 粒子 2 個。置換の扇（恒等・互換）だけを標本にする模型。 rng = np.random.default_rng(20261008) N = 1_000_000 e1 = lambda b: 0.5 / np.tanh(b / 2) # 1 粒子のエネルギー -d ln Z_1/d beta print(" beta exact est (se) E_F exact E_F est (se) N for 1% (E_F)") for beta in (1.0, 2.0, 4.0, 6.0, 8.0): q = np.exp(-beta) s = np.where(rng.random(N) < (1 + q) / 2, 1.0, -1.0) # 恒等: 確率 (1+q)/2、符号 + o = np.where(s > 0, 2 * e1(beta), 2 * e1(2 * beta)) # 扇ごとのエネルギー -d ln(重み)/d beta r = (o * s).mean() / s.mean() y = (o - r) * s # 比の推定量の標準誤差（線形化） se_r = y.std(ddof=1) / np.sqrt(N) / abs(s.mean()) e_f = 2 + 2 * q / (1 - q) - q / (1 + q) v = ((1 + q) / 2 * (2 * e1(beta) - e_f) ** 2 + (1 - q) / 2 * (2 * e1(2 * beta) - e_f) ** 2) / q**2 print(f" {beta:4.1f} {q:.3e} {s.mean():+.3e} ({s.std(ddof=1) / np.sqrt(N):.1e}) " f"{e_f:.5f} {r:+9.4f} ({se_r:.4f}) {v / (0.01 * e_f) ** 2:.2e}") ``` ```text beta exact est (se) E_F exact E_F est (se) N for 1% (E_F) 1.0 3.679e-01 +3.680e-01 (9.3e-04) 2.89501 +2.8947 (0.0029) 1.02e+04 2.0 1.353e-01 +1.340e-01 (9.9e-04) 2.19383 +2.2039 (0.0076) 1.16e+05 4.0 1.832e-02 +1.843e-02 (1.0e-03) 2.01933 +2.0133 (0.0540) 7.31e+06 6.0 2.479e-03 +3.446e-03 (1.0e-03) 2.00250 +1.7218 (0.2087) 4.06e+08 8.0 3.355e-04 +1.068e-03 (1.0e-03) 2.00034 +1.3144 (0.2941) 2.22e+10 ``` 標本数は $N = 10^6$ で固定です。平均符号の標準誤差は $\beta$ によらず $10^{-3}$（$\approx 1/\sqrt N$）ですが、平均符号そのものが $e^{-\beta}$ で小さくなるので、$\beta = 6$ では真の値 $2.5\times10^{-3}$ が標準誤差の $2.5$ 倍しかなく（推定値は $3.4\times10^{-3}$）、$\beta = 8$ では真の値 $3.4\times10^{-4}$ が標準誤差の $1/3$ です。エネルギーの推定値は $\beta \le 4$ では厳密な値と標準誤差の範囲で合いますが、$\beta = 6$ で誤差が $0.2$、$\beta = 8$ では推定値 $1.31$ が厳密な値 $2.00$ から標準誤差の $2.3$ 倍ずれています。分母が誤差に埋もれると、線形化で求めた標準誤差そのものも信用できません。最後の列は の $\sigma_R^2$ から求めた、相対誤差 $1\%$ に要る標本数で、$\beta$ が $2$ 増えるごとに約 $e^4 \approx 55$ 倍になっています。粒子が 2 個でこの有様で、粒子数に比例して指数が増えます。 ## 10. 演習 §9.1 の $2\times2$ の例 $A = \mathrm{diag}(1, 0)$、$B = \frac12\begin{pmatrix}1&-1\\-1&1\end{pmatrix}$ について、$[A, B]$、$[A,[A,B]]$、$[B,[B,A]]$ とそれらの作用素ノルムを計算し、 の定数が $1/16$（対称）と $1/4$（非対称）であることを確かめよ。 $AB = \frac12\begin{pmatrix}1&-1\\0&0\end{pmatrix}$、$BA = \frac12\begin{pmatrix}1&0\\-1&0\end{pmatrix}$ から $C = [A, B] = \frac12\begin{pmatrix}0&-1\\1&0\end{pmatrix}$。これは回転行列の $1/2$ 倍で $\|C\| = 1/2$ です。非対称分割の定数は $\|C\|/2 = 1/4$。 $[A, C] = AC - CA = \frac12\begin{pmatrix}0&-1\\0&0\end{pmatrix} - \frac12\begin{pmatrix}0&0\\1&0\end{pmatrix} = -\frac12\begin{pmatrix}0&1\\1&0\end{pmatrix}$ で、ノルムは $1/2$。$BC = \frac14\begin{pmatrix}-1&-1\\1&1\end{pmatrix}$、$CB = \frac14\begin{pmatrix}1&-1\\1&-1\end{pmatrix}$ から $[B, C] = \frac12\begin{pmatrix}-1&0\\0&1\end{pmatrix}$ で、$[B,[B,A]] = -[B, C]$ のノルムも $1/2$。対称分割の定数は $\frac{1/2}{24} + \frac{1/2}{12} = \frac1{48} + \frac1{24} = \frac1{16}$ で、出力の `c_sym = 0.0625`、`c_lie = 0.25` と一致します。 $(0, a)$ から $(T, b)$ への Brown 橋を二分法で生成する。段 1 で $X_{T/2}$、段 2 で $X_{T/4}$ を作るとき、$X_{T/4}$ の（$X_{T/2}$ について平均した）分布の平均と分散を求め、 の 1 の式と一致することを確かめよ。 段 1：$X_{T/2} = (a+b)/2 + \sqrt{T/4}\,\xi_1$。段 2：区間 $[0, T/2]$（長さ $T/2$）の中点なので $X_{T/4} = (a + X_{T/2})/2 + \sqrt{T/8}\,\xi_2$（$\xi_1, \xi_2$ は独立な標準正規）。代入すると $$ X_{T/4} = \frac{3a + b}{4} + \frac12\sqrt{\frac T4}\,\xi_1 + \sqrt{\frac T8}\,\xi_2 $$ で、独立な正規確率変数の一次結合なので正規分布に従い、平均は $(3a+b)/4 = a + \frac14(b-a)$、分散は $\frac14\cdot\frac T4 + \frac T8 = \frac{3T}{16}$ です。定理の式は平均 $a + \frac{T/4}{T}(b-a)$、分散 $\frac{(T/4)(3T/4)}{T} = \frac{3T}{16}$ で、一致します。 調和振動子で $M = 1$（ビーズ 1 個）のとき、$\pi_1$、$\langle E_{\mathrm T}\rangle$、$\langle E_{\mathrm V}\rangle$、$\mathrm{Var}(E_{\mathrm T})$、$\mathrm{Var}(E_{\mathrm V})$ を直接計算し、 と の式と比べよ。 $M = 1$ では $x_2 = x_1$ でばねの項が消え、$S_1 = \beta\omega^2x^2/2$、$\pi_1$ は平均 $0$、分散 $\sigma^2 = 1/(\beta\omega^2)$ の正規分布（古典的な Boltzmann 分布）です。$E_{\mathrm T} = \frac{1}{2\beta} + \frac12\omega^2x^2$、$E_{\mathrm V} = \omega^2x^2$ で、$\langle x^2\rangle = \sigma^2$ から $\langle E_{\mathrm T}\rangle = \frac1{2\beta} + \frac1{2\beta} = \frac1\beta$、$\langle E_{\mathrm V}\rangle = \frac1\beta$（古典的な等分配則）。$\mathrm{Var}(x^2) = 2\sigma^4$ から $\mathrm{Var}(E_{\mathrm T}) = \frac{\omega^4}{4}\cdot\frac{2}{\beta^2\omega^4} = \frac{1}{2\beta^2}$、$\mathrm{Var}(E_{\mathrm V}) = \frac{2}{\beta^2}$。 命題の式では、$\sinh(\theta/2) = \beta\omega/2$ で $E_1 = \frac{\omega}{2\cosh(\theta/2)}\coth\frac\theta2 = \frac{\omega}{2\sinh(\theta/2)} = \frac1\beta$。分散は $j = 0$（$\kappa_0 = 0$）の項だけで $\frac{1}{2\beta^2}\cdot1$ と $\frac{2\omega^4}{\nu^2} = \frac{2}{\beta^2}$ で、一致します。$M = 1$ では熱力学的推定量のほうが分散が小さく、$M$ が大きくなると逆転します。 の系で、平均符号 $\langle s\rangle = e^{-\beta\omega}$ を相対誤差 $10\%$ で求めるのに要る独立な標本の数を、$\beta\omega = 5$ と $\beta\omega = 10$ について求めよ。 $s = \pm1$ なので $\mathrm{Var}(s) = 1 - q^2$（$q = e^{-\beta\omega}$）で、$N$ 個の平均の相対誤差は $\sqrt{(1-q^2)/N}/q$ です。これを $0.1$ にするには $N = 100(1 - q^2)/q^2 \approx 100\,e^{2\beta\omega}$。$\beta\omega = 5$ で $2.2\times10^6$、$\beta\omega = 10$ で $4.9\times10^{10}$ です。温度を半分にすると必要な標本数が 2 万倍になります。 ## 参考文献 1. D. M. Ceperley, "Path integrals in the theory of condensed helium", *Reviews of Modern Physics* 67 (1995), 279–355 — 原始近似、二分法と多段の更新、エネルギーの推定量、置換の標本抽出。 2. H. F. Trotter, "On the product of semi-groups of operators", *Proceedings of the American Mathematical Society* 10 (1959), 545–551 — 積公式の原論文。 3. A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, S. Zhu, "Theory of Trotter error with commutator scaling", *Physical Review X* 11 (2021), 011020 — 分割の誤差を交換子で評価する一般論。2 次の分割の定数 $1/24$、$1/12$ はここにある形。[arXiv:1912.08854](https://arxiv.org/abs/1912.08854) 4. T. Jahnke, C. Lubich, "Error bounds for exponential operator splittings", *BIT Numerical Mathematics* 40 (2000), 735–744 — 非有界な作用素に対する Strang 分割の誤差。 5. M. F. Herman, E. J. Bruskin, B. J. Berne, "On path integral Monte Carlo simulations", *The Journal of Chemical Physics* 76 (1982), 5150–5155 — ビリアル推定量と、熱力学的推定量の分散が分割数とともに増えること。 6. M. Troyer, U.-J. Wiese, "Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations", *Physical Review Letters* 94 (2005), 170201 — 符号問題の計算量。[arXiv:cond-mat/0408370](https://arxiv.org/abs/cond-mat/0408370)