物性論における数値計算について, 初貝、2001 年夏

1

物性論における数値計算について
東京大学大学院工学系研究科 物理工学 初貝 安弘1

2004 年補足 版
目次
1

物理的問題：局在と遍歴

22

多体問題における数値的対角化
2.1 冪乗法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Lanczos 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 固有ベクトルについて . . . . . . . . . . . . . . . . . . . . . . . . .

4
4
5
8

3

多粒子系におけるモンテカルロ法
9
3.1 全積分の困難と確率的手法 . . . . . . . . . . . . . . . . . . . . . . . 9
3.2 Importance サンプリング . . . . . . . . . . . . . . . . . . . . . . . 10
3.3 負符号の問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

密度行列繰り込み群
4.1 (純粋状態に対する) 密度行列 . . . . . . . . . . . . . . . . . . . . .
4.2 部分的なトレースと (混合状態に対する) 密度行列 . . . . . . . . . .
4.3 近似関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 特異値分解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

陥りがちな数値計算の非常識
19
5.1 計算資源は有限である . . . . . . . . . . . . . . . . . . . . . . . . . 19
5.2 計算精度は有限である . . . . . . . . . . . . . . . . . . . . . . . . . 23

6

ふれられなかった話題

13
13
13
15
17

28

物性論における数値計算について, 初貝、2001 年夏1

2

物理的問題：局在と遍歴

物性論、特に凝縮系でいたるところにでてくる大問題として遍歴性と局在性の競合がある。これは、量子力学の基礎である粒子の波動性と粒子性の２面性に起因をもつとも言える。遍歴性とはある原子上にある電子が近くにある原子上の電子との波動関数との重なりのため電子波として拡がる性質である。この遍歴性の最も顕著な例が次のような自由電子系である。
 †H0 = −tci cj + c†j ci
i,j

ここで i は原子の場所を指定し、ci は原子 i 上の電子の消滅演算子である ({ci , c†j } =δjj  )。−t が飛び移り積分で、これにより電子は系全体に拡がる傾向を持つ。特に原子が規則的にならんでいる場合電子はブロッホ状態として完全に拡がる。
では具体的に N 個の原子が鎖状にならんでいる場合を周期的境界条件で考えよう。この場合ハミルトニアンは次のように書ける。
⎞
⎛
0 1
1
⎛
⎞
⎟
⎜ 1 0 1c1⎟⎜⎟⎜
⎜ c2 ⎟
⎟
⎜
1 0 1
⎟
† †
†
†
†
⎟, c = ⎜H0 = c H 0 c, c = (c1 , c2 , · · · cN ), H 0 = −t ⎜⎜ .. ⎟ ,
.. ..
⎟
⎜
⎝
.
.
. ⎠
⎟
⎜
⎟
⎜
⎝
1 0 1 ⎠
cN
1
1 0
この系でグランドカノニカル集団による統計力学の議論をしよう。まず、分配関
数 Ξ を求めることを考える。2

Ξ = Tr2L e−β(H0 −μNt ) = e−βJ†

ここで Nt = c† c = j=1 cj cj である。3 このトレースを具体的に計算するために
ここでいわゆる波数表示に移ろう。つまり N × N のユニタリ行列 U を使って4N

12π
d = U c, {U }mn = √ exp(i mn)
NNm
dk = {d}m , k = 2πNとすると、この dk もフェルミ演算子の反交換関係をみたし5
H0 =k nk , Nt =nk , nk = d†k dk ,kkとなる。d†k は定義からわかるように空間すべての原子上に有限の振幅を持つ電子を生成し、この系の電子の遍歴性を示している。6 また nk は遍歴状態 k の占有数

物性論における数値計算について, 初貝、2001 年夏

3

演算子となる。この表示ではハミルトニアンも粒子数も対角的なので容易にトレースが計算できてNΞ=(1 + e−β(k −μ) ) = detN (I N + e−β(H0 −μI N ) )

(1)

j=1

となる。7 ここでトレースをとる空間の次元が 2N から N へ著しく減少したことに注意しよう。8 この系の固有状態と固有値 E は各一粒子状態 k の占有数 nk を指定することで次のように書ける。
|{nk } =(d†k )nk |0, E = E({nk }) =
k nk
k

k

特に基底状態は k の小さい方から粒子数だけ詰めた状態となり、系が周期的な場
合には Fermi Sea という。9一方、電子系は粒子として考えたとき負に帯電した粒子であるため静電エネル
ギーによる粒子間の相互作用としてのクーロン反発力をもつ。この効果は最も簡
単には隣り合った原子上の電子間の反発力としてモデル化できハミルトニアンとしてはHint = V
ni nj
<i,j>

と書ける。ここで ni = c†i ci は i 番目の原子上にいる電子の個数演算子である。この項により、電子系は隣り合った原子上にくるとエネルギーを損することとなり可能ならば一つおき以上の間隔をあけて原子上に局在した状態がその基底状態となる。電子の局在性はこのクーロン相互作用により最も典型的に表現されるのである。
この遍歴性をもたらす H0 と局在性をもたらす Hint はお互いに競合し興味深い効果をもたらす。ここで H0 を対角的にする nk は波数表示の粒子数演算子であり、Hint を対角的にする ni は実空間のある原子上の電子数の粒子数演算子であり両立しないところに困難と逆に興味があるわけである。


厳密解の知られていない問題

d 次元格子上有限系 Ld サイト上に M 個の粒子をおいたとき次のハミルトニ
アンを考え
HLd = H0 + Hint
その全固有値と固有ベクトルを求め、固有ベクトルを用いて相関関数を求めよ。



この問題が解ければ遍歴と局在に関する重要な問題がかなりの部分理解できることとなる。
もし V = 0 であれば波数表示を用いることでこの問題に答えることができる。
10

物性論における数値計算について, 初貝、2001 年夏

24

多体問題における数値的対角化

前節で与えた解けていない問題に対して数値的になにがしか取り組むことを考
えよう。
 †

HLd = −tci cj + c†j ci + V
ni nj
i,j

<i,j>

具体的には系の大きさ L を有限とし粒子数 M も有限として具体的に固有値固有ベクトルを数値的に求めることを考えよう。まずヒルベルト空間の次元 D を数えてみよう。並進対称性等の対称性を利用すれば必要な次元は減らせるが、ここでは簡単のため粒子数が保存されることだけを用いよう。すると例えば 1 次元 20 サ
イトの系に 10 個粒子がある場合のヒルベルト空間の次元は

D=

Ld
M

=

20
10

= 184756となるが非常に大きな数となりハミルトニアン行列をそのままの行列として扱うとなるとこれを格納するメモリーは実数倍精度でD ∗ 2 ∗ 8/(1024 ∗ 1024 ∗ 1024) ≈ 254GBとなり現時点でのワークステーションのメモリーが約 2GB 程度であることを考えれば格納不可能である。この様に多粒子系のヒルベルト空間のの次元は L, M が少し大きくなると非常に大きな数となることに注意しよう。
ただし問題によってはこの程度の大きさの系から学ぶことも多く、物理的には低温の物理に限れば特に基底状態とその近傍の固有空間の情報で十分なことも多い。
この事実に注意すれば行列自体を扱うのは不可能でも状態ベクトルはぎりぎり扱えることが多い。そこで基底状態近傍の状態のみを扱う数値的手法が有用となる。この典型的な例として次の冪乗法と Lanczos 法を説明しよう。

2.1

冪乗法

ある状態ベクトル |ΨT  が与えられたときハミルトニアン行列を何度も作用させることを考えよう。
H k |ΨT 
ここで |ΨT  を次のように固有空間分解すれば

ΨT  =

D

cj |j

j=1

H|j = Ej |j

物性論における数値計算について, 初貝、2001 年夏

5Ejmax を絶対値最大の固有値として
H k |ΨT  ≈ Ejkmax |jmax 
このときの相対誤差はおよそ

Ejmax

k

Ejnext max)と見積もれる。すなわち行列を十分な回数書ければ絶対値最大の固有値と固有ベクトルを得ることができることとなる。ここで重要なことはハミルトニアン行列を全体で格納することはここでは必要でなく、ある状態ベクトルにハミルトニア
ンを作用させることができれば十分であることにある。
よってある物理的ハミルトニアン H に対して C を十分大きな数としてH  = H − CIをハミルトニアンと考えれば基底状態の情報が得られることとなる。

2.2

Lanczos 法

理論的には冪乗法と類似するがより精度よく効率的であるためよく使われる方法がここで説明する Lanczos 法である。この方法によると次の手続きにより、ハ
ミルトニアン行列を変形すると考える。
まず、任意の規格化された初期ベクトル|Ψ1 , Ψ1 |Ψ1  = 1
から始めよう。これから|Ψ2  = c2 H|Ψ1  − a1 |Ψ1 と作る。c2 は |Ψ2  の規格化から
c2 = ||big(H|Ψ1 − a1 |Ψ1  ||と定める。ここで係数 a1 は

Ψ1 |Ψ2  = 0となるように決めよう。つまり、|Ψ2  は |Ψ1  に垂直すなわち |Ψ1  方向の成分がないようにするわけである。よって

a1 = Ψ1 |H|Ψ1 = H11

物性論における数値計算について, 初貝、2001 年夏

6

ここで

Hab = Ψa |H|Ψbとした。次に |Ψ3  を|Ψ3  = c3 H|Ψ2  − a2 |Ψ2  − b2 |Ψ1 と作る。ここでもできた |Ψ3  はこれまで作った |Ψj  (j = 1, 2) にすべて垂直すな
わち |Ψj  方向の成分をすべてもたないようにしよう。（また規格化から c3 は定める。）よってまず、

Ψ2 |Ψ3  = 0= Ψ2 |H|Ψ2 − a2 Ψ2 |Ψ2  − b2 Ψ2 |Ψ1 = Ψ2 |H|Ψ2 − a2a2 = H22
つぎにΨ1 |Ψ3  = 0= Ψ1 |H|Ψ2 − a2 Ψ1 |Ψ2  − b2 Ψ1 |Ψ1 = Ψ1 |H|Ψ2 − b2よって

b2 = Ψ1 |H|Ψ2 = H12
これを繰り返して |Ψn+1  を|Ψn+1  = cn+1

H|Ψn  − an |Ψn  − bn |Ψn−1 − · · · − dn |Ψn−1−  − · · ·,  = 1, 2, · · ·と作る ( ≥ 1)。ここでもできた |Ψn+1  はこれまで作った |Ψj  (j = 1, · · · , n) に
すべて垂直すなわち |Ψj  方向の成分をすべてもたないようにし、cn+1 は規格化から定める。ただし帰納的に次の事実を仮定する。

Ψi |Ψj  = 0, i = j,

i, j ≤ nよって

Ψn |Ψn+1  = 0
より

an = Ψn |Ψn  = Hnn

物性論における数値計算について, 初貝、2001 年夏7またΨn−1 |Ψn+1 = 0= Ψn−1 |H|Ψn  − bn
これから

bn = Ψn−1|H|Ψn  = Hn−1,nさらにもう少し戻ってΨn−1− |Ψn+1  = 0より

Ψn−1− |H|Ψn  = Hn−1−,n = dn
一方で Ψn−1− |H は Ψn−−k |、 k = 0, 1, · · · の線形結合でかかれることから |Ψn との行列要素はなく
dn = 0となる。これはここで作った規格直交基底 {|Ψj } でのハミルトニアンの行列は３重対角行列となることを意味する。よってまとめると


Lanczos 法{H}ij ≡ Ψi |H|Ψj 
⎛
O
a1 b2
⎜ b∗ a2 b3
⎜ 2
H = ⎜
b∗3 a3 b4
⎝
.. ..
.
.
O

⎞
⎟
⎟
⎟
⎠

an = Ψn |H|Ψn , bn = Ψn−1 |H|Ψn 
ただし

|Ψ1 
は規格化された（基底状態と交わりのある）任意のベクトルで以下順次|Ψ2  = c2 H|Ψ1  − a1 |Ψ1 |Ψn+1  = cn+1 H|Ψn  − an |Ψn  − bn |Ψn−1  , n = 2, 3, · · · ,により定める。ただし cn は規格化の条件

Ψn |Ψn  = 1
より定める。

物性論における数値計算について, 初貝、2001 年夏

8

この H の絶対値最大、もしくはその近傍の状態をもって H の絶対値最大、もしくはその近傍の状態の「候補」とするのである。具体的にはこの行列を適当なところで切って（次元を有限として）絶対値最大の固有値が収束することを確認する。これをもってそれまでに構成したベクトルが張る空間に考える固有空間が含まれたとするわけである。
（さらに初期ベクトルをいろいろ変えてみて「候補」を実際的にはずす。）

2.3

固有ベクトルについて

前節の Lanczos 法により、基底状態のエネルギー Eg がわかった時、次に基底状態の状態ベクトル |G を求めることを考えよう。前節の Lanczos 法の手続きを繰り返すことにより、このベクトルはもとめられる。これについてここでは説明しよう。（より正確には逆反復法と共役勾配法等の組み合わせが使われる。)
まず |G を次のように展開する。
|G =cj |Ψj j

ここで Ψj  は Lanczos 法の構成過程で現れたベクトルである。よって

H|G = E|G
より

⎞
⎞
⎛
c1
c1
⎟
⎟
⎜
⎜
H ⎝ c2 ⎠ = E ⎝ c2 ⎠
..
..
.
.
⎛

これは係数ベクトルが H の絶対値最大の固有ベクトルであることを示す。よって
前節の固有値の収束性は数列
c1 , c2 , · · · が十分早く 0 に収束することを意味する。
この係数ベクトルは H が三角行列であることより容易に求められるがその具体的手続きをここでまとめておこう。
c1 = 1 として固有方程式を具体的に書き下して

a1 c1 + b2 c2 = Eg c1
b2 c1 + a2 c2 + b3 c3 = Eg c2
bj cj−1 + aj cj + bj+1 cj+1 = Eg cj

c1 = 1(Eg − a1 )c1→ c2 =b2(Eg − a2 )c2 − b2 c1→ c3 =b3···(Eg − aj )cj − bj cj−1
→ cj+1 =,bj+1

これと Lanczos 法からもとまる |Ψj  から |G が求まる。

j = 2, 3, · · ·物性論における数値計算について, 初貝、2001 年夏

39多粒子系におけるモンテカルロ法

一般に多変数のパラメーター σ ∈ S により配位空間 (S) の点が定まりその点ごとにある物理量の値が G(σ) と決まるような状況を考え、重みつきの集団平均 (アンサンブル平均)W (σ)G(σ)
Gσ = σσ W (σ)を計算する事を考えよう。

3.1

全積分の困難と確率的手法

実際の場合、配位空間の積分は高次の多重積分となることが多い。例えば 3 次
元の古典的な N 粒子系の場合
 = dx1 dy1 dz1 · · · dxN dyN dzNσと 3N 重の積分となる。特に熱力学的極限等を議論するためには N が十分大きな場合を考えなければならず、3N は非常に大きくなる。この積分例えば分配関数

Z = dx1 dy1 dz1 · · · dxN dyN dzN W (r1 , · · · , rN ), W = e−H(r1 ,··· ,rN )/kB Tを計算することを考えよう。ここで H(r1 , · · · , rN ) は粒子の位置を r1 , · · · , rN を決めれば決定する古典的エネルギーとする。系が一辺 L の立方体に入っているとして積分を長さ L の区間を M 等分して計算するとしよう。すると加えるべき項の数はM 3Nとなるが、これまた非常に大きな数となる。例えば N = 10, M = 20 として

log M 3N / log 10 = 39.0309と約 40 桁の数となり、とても計算できない。
そこで考えられるのが確率的手法である。よく使われる簡単な例として２次元の一様乱数を発生させて円の面積を計算する手法があげられる。これを単純に適用すれば区間 [ − L/2, L/2] の一様乱数を 3N × Nsample 個使ってZ≈

1
Nsample

NsampleWjj=1と計算することになる。しかし物理的な状況ではこの方法はほとんど役に立たない。11

物性論における数値計算について, 初貝、2001 年夏

10

そこで分配関数自体には適用できないが、最初にあたえたような、物理量の平均値の計算に限って考えよう。
ここでもしすべての配位で 0 ≤ W (σ) の場合（古典的な例では常にそうだが）Gσ =

P (σ)G(σ)σ

P (σ) =となり

0 ≤ P (σ) ≤ 1,W (σ)σ W (σ)P (σ) = 1σであるから P (σ) を配位 σ の実現確率と解釈すれば、各配位 σ を確率 P (σ) 発生させることができればその確率的期待値がここで必要とする積分値となる。

3.2

Importance サンプリング

一般に重要となるのは、配位空間 S が非常に大きくかつその分布関数 P (σ) が非常に偏った構造をしている場合がおおい。12 このような時に、ある程度の大きさを持つ有限系において物理量の集団平均を評価することを考えると σ の和をす
べてとることは現実的でなく、13 またランダムに配位空間の点を選んで単純平均
しても意味のある結果は得られない。14
そこで次のような Importance Sampling と呼ばれる効率的な方法を用いる。ポイントはそこでの時間平均と集団平均が等しくなるような確率過程を作ることにある。

1. 任意の初期状態 σi を選ぶ。
2. 状態 σi から状態 σj に確率 p(i → j) で状態を遷移する。(つまり遷移確率は一つ前の状態にしかよらない。)
この操作を繰り返すことにより確率過程を構成する。ただし遷移確率 {p(i → j)}
は次の２つの条件

1. エルゴード性
任意の２状態間を有限のステップで遷移できる。

2. 微細釣り合いの原理 (detailed balance)
P (σi )p(i → j) = P (σj )p(j → i)

物性論における数値計算について, 初貝、2001 年夏

11を満たすとする。このときこのマルコフ過程を十分なステップ数繰り返すと確率過程において配位 σi の現れる確率 P̃∞ (σ) は一定となり (確率過程は収束し) 15

P̃∞ (σ) = P( σ)1= GM C ( = lim
G(σj ))
T →∞ T
j=n
n+T

Gσ

となり、集団平均が仮想的時間ステップにかんする時間平均に等しくなる。
具体的に上記の条件を満たす遷移確率を構成する方法２つを紹介する。16

1. メトロポリス法(a) ランダムに状態 t を選ぶ(b) (1) P (σt ) > P (σi ) の場合無条件で σj = σt に遷移する。(2) P (σt ) ≤P (σi ) の場合、確率 r = P (σt )/P (σi ) で σj = σt に遷移する。
(c) 繰り返す。
2. 熱浴法(a) ランダムに状態 t を選ぶr(b) r = P (σt )/P (σi ) として確率 1+r
で σj = σt に遷移する。

(c) 繰り返す。

3.3

負符号の問題すべての配位で 0 ≤ W (σ) の場合前節の議論がつかえるがフェルミオン系の場
合 W (σ) = det OF であるから必ずしも正定値の条件が成立しない。そこで集団平均の式を次のように変形する。

Gσ =σ W (σ)G(σ)σ W (σ)|W (σ)|
sgnW (σ)G(σ)
σ |W (σ )||W (σ)|
sgnW (σ)σ
σ |W (σ )|σ P̄ (σ) sgnW (σ)G(σ)σ P̄ (σ) sgnW (σ)
sgnW (σ)G(σ)P̄sgnW (σ)P̄
|W (σ)|
σ |W (σ)|σ===
P̄ (σ) =物性論における数値計算について, 初貝、2001 年夏

12

この式で分母、分子をそれぞれ importance sampling することを考える。
P̄ (σ) = 1
0 ≤ P̄ (σ) ≤ 1,σであるから P (σ) を配位 σ の実現確率と解釈できるのである。つまり

Gσ =

sgnW (σ)G(σ)M C
sgnW (σ)M C

ただし符号の平均 sgnW (σ)M C が非常に小さくなる場合小さな数同士の比から有限の値を導くことになり誤差が非常に大きくなる。これがいわゆる負符号の問題である。

物性論における数値計算について, 初貝、2001 年夏4

密度行列繰り込み群

4.1(純粋状態に対する) 密度行列13

まず最初に基底状態に関する密度行列を復習しよう。
基底状態が縮退がなく唯一である場合を考えよう。この規格化された基底状態|ψ に対して密度行列 ρ を次のように定義する。

ψ=
|IψI
I

1 =ψ|ψ =

|ψI |2Iρ =|ψψ| =|IψI ψJ∗ J| =I,J|IρIJ J|I,J(ρ)IJ =ρIJ = ψI ψJ∗
ここで |I, I = 1, · · · は系の完全系である。この密度行列は次の関係式を満たす。
Tr ρ =
ψI ψI∗ = 1
I
2

ρ =ρ,

(ρ2 )IK =∗ψI ψJ∗ ψJ ψK
=

IJK



∗
ψI ψK
= (ρ)IK

IK

これより密度行列の固有値 ρ, λI は 1 か 0 であることが次のように示せる。

ρvI =λI vI ,ρ2vI = λ2I vI = λI vI ,
λI =Tr ρ = 1

λI (λI − 1) = 0

I

さらに固有値 1 の固有状態は縮退がなく他の固有値は 0 で
⎞
⎛
ψ1
⎟
⎜
vψ = ⎝ ... ⎠ψDψI ψJ∗ ψJ = ψI = (vψ )I(ρvψ )I =
J

となる。

4.2

部分的なトレースと (混合状態に対する) 密度行列

つぎに物理系が２つの部分系, in と out とにわかれている場合を考えよう。このとき系の完全系のラベル I は２つの組 I = (iin , jout ) で指定され|I =|iin |jout 物性論における数値計算について, 初貝、2001 年夏

14

全系の基底状態 |ψ は次のように展開されることとなる。
|ψ =ψI |I
I=

ψij |iin |jout 

iin jout

よって (全系の) 密度行列はρ=|iin |jout ψij ψi∗ j  i in |j  out |,ij,i j 

|ψij |2 = 1

ij

となる。
ここで out-系についての部分トレースを次のように定義し reduce された密度行
列 ρ̃ を次のように定義する。

ρ̃ =Tr out ρk out |==|iin (ii

≡

|iin |j out ψij ψi∗ j  i in |j  out ||k out iji j 

k






ψij ψi∗ j )i in |j|iin (ρ̃)ii i in |ii(ρ̃)ii =ψij ψi∗ jj

この 密度行列に ρ̃ 対しては次の関係式が成立する。
ρ̃ =|i(ρ̃)ii i |ii

Tr ρ̃ =1(=

|ψij |2 )ij

さらに任意の (in ) ベクトル |vin  = i |ii|vin  に対して
vin |ρ̃|vin  =vin |iψij ψi∗ j i |vin  =|vin |iψij |2 ≥ 0
ii j

i

i

よって

ρ̃ の固有値は全て非負であり、その和は 1 である。
ここで |i, i = 1, 2, · · · , M が in 空間での完全系を作るとすれば in 空間にのみ依存する物理量 Oin の期待値 Oin  は次のように in 空間のみで密度行列 ρ̃ を用

物性論における数値計算について, 初貝、2001 年夏

15

いて計算できることとなる。

Oin  =ψ|Oin |ψ ψij∗ ψi j  iin |jout |Oin |i in |j  out =
iin jout i in j  out=ψij∗ ψi j  iin |Oin |i in 

jout iin i in=Tr ρ̃Oin

4.3

近似関数

以上の準備の下で in 空間の基底数を M から減らして M  とし全系の波動関数|ψ を近似することを考えよう。

|ψ =

M 

i=1

ψij |i|j ≈ |ψ̄j

ここで全空間の一般の状態ベクトル |X を次のような M × N 行列として表示
すると (N は j 空間の次元)

|X =

M 
N


Xij |i|j

i=1 j=1

N

Xij =(X)ij , X =X

M|X, |Y  の内積が次のように書けること
Y |X =||X||2 =Xij∗ Yijij=Tr X † Y
および out 空間に対するトレースが次のように書けることに注意する。
Tr |XY | =|iXij Yi∗ j i |j=|i(XY † )ii i |ii

長方形行列 ψ に対する特異値分解と呼ばれる次のような標準型を用いると

物性論における数値計算について, 初貝、2001 年夏16ψ =U DV †M

U= M= u1 · · · uM : UnitaryX=(|u1 , · · · , |uα , · · · , |uM )αuα∗i ui = δααiD =diag (d1 , d2 , · · · , dM ), di ≥ dj ≥ 0, (i > j)
N

V

†



= M

V †



v1
..
.
vM

=

,



V † V  = IMviα∗ viα = δααi

全系の基底状態は次のように書ける。

|ψ =

M

ij

=

M


uαi dα vjα |i|j

α=1

dα |uα|v α 

α=1|u  =α

uαi |i, |v α  =iviα |jjαu |u  = δαα ,αv α |v α  = δαα

ここで α の和を 1 → M から 1 → M  に制限して近似的な波動関数 |ψ̄ を作ることとしよう。

|ψ̄ =Mij

uαi dα vjα |i|j

α=1

M

=



dα |uα|v α 

α=1ψ̄ =U D̄V †物性論における数値計算について, 初貝、2001 年夏

17

このときの波動関数の近似の程度は次の量から読み取れる。


r =|| |ψ − |ψ̄||2 = Tr (U δDV † )† (U δDV † ),=Tr δDδD=MδD = D − D̄†|dα |2α=M  +1

さらに reduce された密行列はρ̃ =ψψ † = U DV † V D† U †=U D2 U †となり、

ρ̃U =U D2と密度行列 ρ̃ を対角化するベクトルから新しい in 空間の基底 uαi が構成できることとなる。
この議論をくり返し用いることにより多体系の基底状態を近似的に表示する手法
を密度行列繰り込み群とよぶ。(S. White. PRL 69, 2863, (1992), PRB. 48, 10345(1993))

4.4

特異値分解

最後にここでの議論に用いた以下の行列 ψ に対する特異値分解についてまとめてみよう。ただし det ψψ † > 0 を仮定しよう。まず次の M × M 型であるエルミート行列を対角化するユニタリ行列を U とする。

ψψ † =ψ†ψψψ † =U D̃U † =

U

ψψ †

=

D̃

U†
†

ここで (ψψ † )† = ψψ † と仮定 ( det ψψ † > 0 ) よりこの対角行列 D̃ = D̃ は実正定値である。(対角要素全て正) これからV =ψ † U
V

=

ψ†

U

物性論における数値計算について, 初貝、2001 年夏

18と定義すればV † V =U † ψψ † U = D̃
V†

=

V

D̃

よって

V  =V D̃

−1/2

として
†

V  V  =I
まとめて

ψ =U V † = U D̃=U DV 
D =D̃

−1/2

†

−1/2

V

†

物性論における数値計算について, 初貝、2001 年夏

519

陥りがちな数値計算の非常識

ここでは陥りがちな数値計算の常識についていくつか列挙して注意を喚起したい。(2001 年時点) 他にも数値計算の本例えば (「数値計算の常識」伊理他 等) を参照されたい。

5.1

計算資源は有限である計算機のメモリー、計算可能な時間は明らかに有限である。

• 行列は大きい行列

A
をとしてその固有値問題を扱うことを考えよう。（実数倍精度として記憶容量を目安として計算してみよう。)
このときその次元が1000 × 1000, 15MB
はまあまあがんばれば扱えるが、10000 × 10000, 1.5GB
はかなり難しい、

10 万 × 10 万, 152GBは無理。メモリーも厳しいが計算時間のほうでも限界であろう。

• 多体問題のヒルベルト空間は大きい
粒子数を保存する系の単純に考えたときのヒルベルト空間の次元

DN M =

N
Mを計算してみよう。ここでは M = N/2 として実数倍精度で記憶容量を見積
もろう。

D8,4 : 0.5MB, D16,8 : 98MB, D24,12 : 20.6GB, D36,18 : 69T B
これはとても大きい。直接扱えるのは小さな系だけである。

物性論における数値計算について, 初貝、2001 年夏

20• 行列のかけ算は順々にA, B, C を N × N 次の行列としたとき
ABCを求めることを考えよう。このとき次の「BAD」なプログラムと「GOOD」
なプログラムを比べてみよう。
特に行列の次元 N が大きい時はこの差は顕著である。試してみると N = 200である計算機で 0.84 秒と 93 秒とでた。（原理的には 200 倍違うはずだ）また行列の書ける数が多いときにはますます違う。決してこれはやってはいけない。

物性論における数値計算について, 初貝、2001 年夏// BAD#include <iostream>using namespace std;
#define N 100
double a[N][N];double b[N][N];double c[N][N];double result[N][N];int main (){for (int i=1;i<N;i++) {for (int j=1;j<N;j++) {a[i][j]=1.0;b[i][j]=1.0;
c[i][j]=1.0;}}
for (int i=1;i<N;i++) {for (int j=1;j<N;j++) {double s=0.0;
for (int k=1;k<N;k++) {for (int m=1;m<N;m++) {s=s+a[i][k]*b[k][m]*c[m][j];}}result[i][j]=s;}}for (int i=1;i<N;i++) {for (int j=1;j<N;j++) {cout << result[i][j]<< "\n";}}return 0;}21物性論における数値計算について, 初貝、2001 年夏// GOOD#include <iostream>using namespace std;
#define N 100
double a[N][N];double b[N][N];double c[N][N];double result[N][N];double tmp[N][N];int main (){for (int i=1;i<N;i++) {for (int j=1;j<N;j++) {a[i][j]=1.0;b[i][j]=1.0;
c[i][j]=1.0;}}
for (int i=1;i<N;i++) {for (int j=1;j<N;j++) {double s=0.0;
for (int k=1;k<N;k++) {s=s+a[i][k]*b[k][j];}tmp[i][j]=s;}}for (int i=1;i<N;i++) {for (int j=1;j<N;j++) {double s=0.0;
for (int k=1;k<N;k++) {s=s+tmp[i][k]*c[k][j];}result[i][j]=s;}}for (int i=1;i<N;i++) {for (int j=1;j<N;j++) {cout << result[i][j]<<"\n";}}return 0;22物性論における数値計算について, 初貝、2001 年夏

23• 積分のたたみこみ前の例、「行列のかけ算は順々に」を理解しても次の例で同様の bad なプログラムをすることがある。まず、関数
a(k), b(k), c(k)が与えられているとき
g(k) = dk1 dk2 dk3 a(k − k1 )b(k1 − k2 )c(k2 − k3 )単純にそのまま計算すると積分区間の分割数 N として O(N 4 ) の手続きがいるが、これも次のように計算すると O(N 3 ) ですむ。
g(k) =
dk2 dk3
dk1 a(k − k1 )b(k1 − k2 ) c(k2 − k3 )=
dk2 dk3d(k, k2 )c(k2 − k3 )d(k, K) =dk1 a(k − k1 )b(k1 − K)

5.2

計算精度は有限である計算機が扱える数は特別な場合をのぞき、有限の精度しかない。

• 塵はつもらないこともある。
次の２つのプログラムを比較せよ。

物性論における数値計算について, 初貝、2001 年夏//bad#include <iostream>using namespace std;
#define N 20000000

int main (){double s=1.0e8;
double small=1.0e0/((double) N);
for (int i=0;i<N;i++) {s=s+small;}cout.precision(40);cout << s <<"\n";return 0;}24物性論における数値計算について, 初貝、2001 年夏//good#include <iostream>using namespace std;
#define N 20000000

int main (){double s=1.0e8;
double s0=0.0;
double small=1.0e0/((double) N);
for (int i=0;i<N;i++) {s0=s0+small;}s=s+s0;cout.precision(40);cout << s <<"\n";return 0;}25物性論における数値計算について, 初貝、2001 年夏

26• 覚えられる数には限界がある。
計算機のなかでは記憶できる実数には限界がある。例えば 1010000 は記憶できないし 10−10000 もだめである。このようなときは小分けして対数を記憶すればよい。例えば数列 {ai }, i = 1, 2, · · · の多数の積を計算する場合例えば
次の様にする。
17// bad#include <iostream>using namespace std;
#define N 200
int main (){double s=1.0;
double rat=0.5;
for (int i=0;i<N;i++) {s=s*rat;}cout.precision(40);cout << s <<"\n";cout << log(s)/log(2.0) <<"\n";return 0;}

物性論における数値計算について, 初貝、2001 年夏//good#include <iostream>using namespace std;
#define N 20000

int main (){double s=1.0;
double log_of_s=0.0;
double rat=0.5;
for (int i=0;i<N;i++) {s=s*rat;if((i+1)%100==0) {log_of_s=log_of_s+log(s);s=1;}}cout.precision(40);cout << log_of_s/log(2.0) <<"\n";return 0;}27物性論における数値計算について, 初貝、2001 年夏

6

ふれられなかった話題• 量子モンテカルロ法の詳細• 並列計算28

物性論における数値計算について, 初貝、2001 年夏

29

Notes
1

hatsugai@pothos.t.u-tokyo.ac.jp,

2

サイト表示の基底を構成せよ。更に、3 サイトの場合に分配関数をこの表示
で実際に計算せよ。
3

このトレースをとる空間の次元が 2N である理由を述べよ。

4

U のユニタリ性を示せ

5{dk=2π m, d†k=2π m N } = δkk = δmm これを示せ
N

6

１次元鎖でなく、一般の D 次元正方格子の場合一粒子エネルギーはどうなるか？バンド幅は次元とともにどう変化するか？
7

この式を導け8

3 サイトの場合に分配関数をこの方法で実際に計算せよ。

9

ここまでの議論はハミルトニアンがフェルミ演算子について２次形式であれ
ば最近接かつ周期的な飛び移り積分を持つ場合に限らずに成立する。これが今後の議論において重要となる。
10

V = 0 の時、波数表示でうまくいかない理由を示せ。

11

なぜか？

12

古典熱力学においては W (σ) = e−βE(σ) 、E(σ) は古典系のエネルギーとなる。
よって系がマクロな場合、エネルギーは示量的であるから系の大きさを増やすとともに、その分布は指数関数的にするどいピークを持つものになることが予想される。
13系の大きさとともに配位空間も指数関数的に増大することが多いので

14

分布のピークが非常に鋭いのでランダムに選んだ点はほとんど重要な重みを持つ点にあたらず、「はずれて」しまう。
15

配位 i に n ステップ目で系がある確率を pn (i) として行ベクトル pn を {pn }i =pn (i) と定義するとこの確率過程の遷移行列 W , Wij = p(i → j) ≥ 0 を用いて

物性論における数値計算について, 初貝、2001 年夏

30

pn+1 = pn W となる。なお確率の保存よりpn (i) = 1, ∀inWij = 1, ∀jiとなる。さらに最後の式は 1 をすべての要素が1 からなるベクトルとして

W1 = 1を意味し、W が固有値 1 を持つことを示す。さらにエルゴード性より、適当な整数 k について W̃ はすべての要素が正の行列となり、この行列に対しては収束性が言え、次の極限が存在する。(斉藤:基礎数学線形代数演習:東大出版会参照)

⎛

v1 v2
⎜ v1 v2n⎜lim W̃ = ⎜n→∞⎝ v1 v2
.. ..
. .

⎞
···
··· ⎟
⎟
⎟ , bi > 0
··· ⎠
..
.

W̃ に対してはペロン ・ フロベニウスの定理から W̃ の絶対値最大の固有値 1 は縮退がなく、その左固有ベクトル v は唯一で定符号で次の様に与えられる。
{v}i = vi
1617

これらの方法が微細釣り合いの原理を満たすことを確認せよ。

単純な応用だが小分けにすることがわかれば a100 = (((((a2 )2 )2 )2 )2 )2 ·((((a2 )2 )2 )2 )2 ·((a2 )2 )2 と計算するであろう。