《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

量子アニーリングの数理
東京工業大学 大学院理工学研究科 物性物理学専攻
西森 秀稔

物理学の考え方や手法が，本来の物理学の対象を超えて広範な問題に適用されている。生物物
理学は古くから知られている例であるが，情報科学にも統計力学や量子力学の方法が役に立つことが明らかになってきている。情報統計力学と呼ばれる分野はその典型例である。また，量子力学の原理を用いて情報処理の高速化を目指す量子情報も盛んに研究されている。本講義では，量子力学の手法を用いた最適化問題の汎用解法である量子アニーリングについてその基礎付けに重点を置いて解説する。量子アニーリングは量子情報科学の主流とはやや趣を異にする分野であるが，汎用性やノイズに対する強さなどの長所が注目されて，理論だけでなく実験的にも研究が盛んになってきている。量子アニーリングの最初の提案者のひとりとして，独自の視点からその特性の解析や問題点を丁寧に解説していきたい。

1

組み合わせ最適化問題
離散的な多変数の 1 価関数を最小化する問題を組み合わせ最適化問題という。わかりやすい例

としては巡回セールスマン問題がある。図 1 のように，地図上で与えられた N 個の都市をそれ
ぞれ 1 回ずつ巡って元の都市に戻るという制約条件の下に，経路長を最小にせよという問題である。最小化するべき関数を一般にコスト関数という。巡回セールスマン問題では経路長がコスト関数である。コスト関数にマイナスの符号を付けると最小化が最大化になるから，最小化のみを考察しても一般性は失われない。物理の例としては，イジング模型の基底状態を求める問題が典図 1: 巡回セールスマン問題の例。すべての都市を一度ずつ巡って元の場所に戻る最短経路を求める。

型的な例である (図 2)。ハミルトニアン（エネルギー）がコスト関数である。もちろん，相互作《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

図 2: イジング模型の基底状態。単純な強磁性相互作用の場合（左）とスピングラスの場合（右）。

用が強磁性的なら基底状態は自明（すべてのスピンが上向きか下向き）であるが，複雑な相互作用を持つスピングラスなどでは，基底状態を求めるのは一般に極めて難しい問題であり，非自明な組み合わせ最適化問題の一例と言える1 。
ある組み合わせ最適化問題が難しいかどうかの基準は，現在知られている最善のアルゴリズム
がその問題を解くのに要する時間が問題の大きさの指数関数に比例するか，あるいはべき関数に比例するかである。指数関数は引数の増加と共に急激に増大するから，コンピューター上で解を見つけるのはちょっと問題が大きくなるとすぐに出来なくなる。これに対して，計算時間がべきで増大するアルゴリズムが見つかっていれば，かなり大きな問題でも正しい解に到達することが出来る。問題が解けるかどうかではなくて，解くための時間が問題のサイズと共に急激に（指数関数的に）増大するかどうかが基準になるのである。指数関数的に時間が増大する問題を困難な問題という。巡回セールスマン問題や 3 次元以上のスピングラスの基底状態を求める問題は困難な問題として知られている。
以上は，与えられた問題の厳密に正しい解を求めるためのアルゴリズムの話であるが，実用上は近似解で十分な場合も多い。汎用性の高い近似解法として，シミュレーテッド・アニーリングがよく知られている。シミュレーテッド・アニーリングにおいては，コンピューター上で仮想的に温度変数を導入して温度を時間の変数として適切に制御することにより，最適解に出来るだけ近い状態を求める（図 3）。シミュレーションの基礎理論は付録 C に述べた。シミュレーテッド・
アニーリングは実用的に広く用いられており，統計力学が情報科学に有用な寄与をした顕著な例である。
量子アニーリングは，シミュレーテッド・アニーリングの量子力学版と言えよう [3] -[6]。温度の代わりに量子ゆらぎを導入し，古典的な確率過程の代わりに量子力学的な状態の重ね合わせをうまく利用して組み合わせ最適化問題の解を探索するのである (図 4)。この際，量子ゆらぎの種類や強さをどう選ぶかが成功の鍵となる。
本講義の目標のひとつは，熱ゆらぎと量子ゆらぎの最適化問題に対する有効性の比較である。

2

量子計算と量子アニーリング
量子計算理論の主要な目標は，組み合わせ最適化問題の高速解法アルゴリズムを量子力学の原理に基づいて構成することにある [7]。通常の量子計算の枠組みでは，量子ビットの組からなる入力に対してゲート（素子）においてユニタリ変換を逐次的に作用させていき，最終的に得られた量子ビットの組に観測をかけて出力を取り出す (図 5)。これに対して量子アニーリングでは，時
1 スピングラスについては [1, 2] を見よ。

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

図 3: シミュレーテッド・アニーリングでは。熱ゆらぎにより様々な状態の間を確率的に頻繁に遷移させる。ゆっくりと温度を下げていくと，各時刻において近似的に平衡状態が実現する。最後に温度を 0 にすると，絶対零度での平衡状態である基底状態（あるいはそれに近い状態）が高い確率で実現する。

図 4: 量子アニーリングでは量子ゆらぎによる量子トンネル効果を利用して様々な状態の重ね合
わせを作ることにより状態探索をする。最初は量子ゆらぎを大きく取って多くの状態を重ね合わせて同時に探索し，次第に量子ゆらぎの強さを押さえていくことにより最終的にコスト関数の最小状態に行き着くことを目指す。

間に依存した制御変数を持つハミルトニアンを用意し，シュレディンガー方程式に従って系の状態を連続的に変化させる。そして，最終的に得られる状態（波動関数）が高い確率でもとの最適化問題の解を表現しているよう初期状態や時間に依存したハミルトニアンを設計するのである。
通常の量子計算では，解を高い確率で得るのに必要なユニタリ変換の数（ゲートの数，ステップ数）が問題の大きさとともにどう増大するかが計算量の目安になる。これに対して量子アニー
リングでは，所期の最終状態に到達するまでの時間の長さが考察の対象となる。
量子アニーリングが最近注目を浴びているのにはいくつかの理由がある。通常の量子計算では，外界との相互作用による系の乱れ（デコヒーレンス）の影響をどれだけ抑えることが出来るかが，量子計算の成功の鍵となる。本質的にミクロな世界である量子系を特定の状態のまま乱れなしに保ち，さらに能動的な制御を施すことは実験的にきわめて難しい課題であり，多大な努力がなされている。量子誤り訂正符号などの技術が開発されているが，その実現に必要な技術的課題は多い。一方，量子アニーリングで利用する基底状態は，各種の擾乱に対して比較的安定である。例
えば温度効果（熱擾乱）については，基底状態と第一励起状態の間のエネルギーギャップと比較して温度が十分低い限り無視できる。

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

図 5: 量子計算のゲート模型では，量子状態に次々にユニタリ変換を施して状態を変換していく。
各横線はひとつの量子ビット，箱はひとつまたは 2 つの量子ビットに作用させるユニタリ変換を表す。左側が入力，右側が出力である。

量子アニーリングにやや遅れて提案された量子断熱計算は，量子ゆらぎをうまく制御して組み合わせ最適化問題を解く方法という点で，量子アニーリングと本質的に同じ汎用アルゴリズムである [8]。ただし，量子断熱計算においては量子力学系の断熱的な（十分ゆっくりとした）時間発
展を追うために必要な計算量という側面が特に強調される。本講義では量子アニーリングと量子断熱計算を特に区別せずに扱う。

量子アニーリングの定式化

3

それでは，量子アニーリングは実際どのように定義されていて，どのようにして困難な最適化
問題の解法として役立つのかについてより詳しく見ていこう。最適化問題のコスト関数を統計力学のイジング模型で表現したのち，量子力学的な揺らぎを横方向の磁場として導入してその強さ
を時間の関数として上手に制御するのである2 。

3.1

最適化問題のハミルトニアン表現

まず，解くべき組み合わせ最適化問題のコスト関数をハミルトニアン H0 で表現する。基底状
態が最適解になるよう H0 を選ぶのである。例えば，N 個の項目の中から特定の 1 個を探すデー
タベース探索問題では，項目を |0i, |1i, · · · , |N − 1i で表し，H0 を

H0 = I − |wihw|(1)とすればよい。|wi (w は 0 から N − 1 のいずれか) が目的の項目であり，I は N 次元の単位行列である。解を表現する状態 |wi が H0 のエネルギー 0 の基底状態になっている。他の項目 |ki (k 6= w)に H0 を作用させると，励起状態のエネルギー 1 が得られる。
統計力学でよく知られたスピングラスの基底状態を求める問題では，ランダムな相互作用を持
つイジング模型のハミルトニアンを最小化する。

H0 = −∑Jij σiz σjz

(2)

hiji
2 人間にとって有用な目的のために何かを「制御する」という発想は通常の物理の考え方と異なっている。自然を理解するのではなく，能動的に働きかけて役立つものを引き出すのである。

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

右辺の和は適当な格子上の隣同士の格子点の組について取る。σiz は格子点 i に割り当てられたパ
ウリスピンの z 成分である。スピングラスのハミルトニアンの基底状態を求めるのは一般には困難な問題である。
一般の組み合わせ最適化問題も多くの場合複雑な相互作用を持つイジング模型で表すことが出来る。組み合わせ最適化問題は離散変数で表されるから，0 か 1 かのビット列の変数に焼き直すことが出来，同じ 2 値のイジングスピンに書き換えられる。ビット間に束縛条件（例えば巡回セー
ルスマン問題で 1 都市は 1 度しか巡れないという条件）があるときには，それを表現するためにやや込み入った考え方が必要になることもある。

3.2

量子ゆらぎの導入

H0 の基底状態を求めるために，量子力学的なゆらぎを用いて状態探索を行うのが量子アニー
リングの基本戦略である。演算子の非可換性が量子性（量子ゆらぎ）の特徴だから，H0 と非可換
な H1 を導入して全ハミルトニアンを)(ttH(t) = H0 + 1 −
H1
T
T

(3)

とする。T は計算時間であり，初期状態 t = 0 で H(0) = H1 ，終状態 t = T で H(T ) = H0 になる。初期ハミルトニアン H1 の基底状態が自明に見つかるように H1 を選ぶことが重要である。

H1 の自明な基底状態を初期状態に選び，H(t) で記述される系を時間に依存したシュレディン
ガー方程式に従って時間発展させる。T が十分大きければ時間発展はゆっくりであり，量子力学の断熱定理が適用できる。すると，系の状態は H1 の基底状態から始まって各瞬間の基底状態を連続的に追って行って，最終的には H0 の基底状態に到達する。これは最適化問題が解けたことを意味する。自明な初期状態から出発して，自然な時間発展により非自明な状態に到達できるというわけである。話がうまいではないか。
式 (3) の代わりにH(t) = H0 + Γ(t)H1(4)と定式化することもよくある3 。初期状態においては制御変数 Γ を十分大きく選んで系が H1 だけで記述されるようにしておき (H(0) ≈ Γ(0)H1 )，Γ を時間とともにだんだん小さくしていって，最終的に 0 に持って行く (Γ(t) → 0, H(t) → H0 , t À 1) のである。t は 0 から始まって十分大きな値まで取る。この場合も，H1 の基底状態から出発して十分長時間の後には，H0 の基底状態である最適解に到達することを目指す。
例えば，(1) の H0 で表されるデータベース探索問題では，すべての状態間の量子力学的遷移を同じ重みで許すような演算子
N −1

H1 = I − |ψ0 ihψ0 |,1 ∑|ψ0 i = √
|ji
N j=0

(5)

がよく用いられる。解 |wi を知らなくても，H1 の基底状態 |ψ0 i は単にすべての状態を同じ重みで足せばよいので自明である。これが初期状態となる。また，(2) のスピングラスでは，横磁場3 量子アニーリングはもともと (4) の形で定式化された。これに対して，量子断熱計算は (3) を念頭に置いている。今日ではどちらも扱われるが，(3) を議論することが多いようである。

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）(x 方向にかけた磁場) を表す演算子
H1 = −

∑

σix

(6)

i

を選ぶことが多い。H1 の基底状態は，{σiz }i を対角化する表示ですべての状態を等確率で足し上げた自明な状態∏ |+ii + |−ii
√
2
i

(7)

である。ここで，|±ii は σiz |±ii = ±|±ii なる状態である。これらすべての状態の中から，目標とする H0 の基底状態が時間発展の後に自動的に選ばれるというのが基本的なアイディアである。
自動的にというところがポイントであり，そんなにうまくいくのかという疑問が当然わくだろう
が実際うまくいくのである。

3.3

断熱変化

ハミルトニアン H(t) は時間変数 t にあらわに依存する。そのため各時刻における系の状態 |ψ(t)i
は，定常状態のシュレディンガー方程式の解とは一般には異なったものになる。しかしながら，Tが十分大きくて系が断熱的に (ゆっくりと) 変化していると見なせる場合には，|ψ(t)i は t を固定したパラメータとみなしたときの定常状態のシュレディンガー方程式

H(t)|j(t)i = ²j (t)|j(t)i,(²0 (t) < ²1 (t) ≤ ²2 (t) ≤ · · · )(8)の基底状態 |0(t)i に十分近い (|hψ(t)|0(t)i| ≈ 1)。
各時刻における定常状態のエネルギー固有値 ²j (t) は，パラメータ t の関数と見なすことが出来る（図 6）。系のサイズが有限のときほとんどすべての場合において基底状態は縮退しておらEnergyt
図 6: 定常エネルギー固有値の t 依存性の模式図
ず，第一励起状態との間にエネルギーギャップ ∆t = ²1 (t) − ²0 (t) が存在する。断熱定理によれば，(3) の H(t) で系を時間発展させた後に，t = T で最適解 |0(T )i に十分近い状態に行き着く条
件 |hψ(T )|0(T )i|2 = 1 − δ 2 (δ ¿ 1) は

max |h1(t)|∂t H(t)|0(t)i|=δmin ∆2t(dH(t)∂t H(t) =
dt

)
(9)

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

である（付録 A 参照）。max と min は，0 ≤ t ≤ T の間での各引数の最大値あるいは最小値である。

(3) によると，上式の分子は 1/T に比例する。よって，
T ∝1
δ min ∆2t

(10)が満たされるように T を選べば，高い確率で正解（H0 の基底状態）に行き着くことになる。エ
ネルギーギャップ ∆t が小さいと状態変化の途中で基底状態から励起状態に遷移する確率が無視できなくなるので，T を大きく選んでゆっくりと系を制御しなければならない。また，δ を小さく取って最後の波動関数が正解である確率を高めようとすると，それに応じて T を大きく取らねばならない。
多くの難しい問題においては，最小のエネルギーギャップ min ∆t は系のサイズに対して指数関数的に小さくなる。一般に，1 次相転移点ではエネルギーギャップが指数関数的につぶれることが多い。(10) によると，エネルギーギャップが指数関数的に減少するとき，最適化問題の解を高い確率で得るための計算時間 T は系のサイズに対して指数関数的に増大し，問題は困難なものとなる。すなわち，量子アニーリングにおいては問題の困難さは 1 次相転移の有無によって左右されることが多い。

アニーリングの収束条件

4

量子アニーリングにおいてハミルトニアンを (4) のように選んだとき，断熱条件 (9) は制御変数 Γ(t) に対する微分方程式を与える。汎用最適化アルゴリズムとしての量子アニーリングの性能評価のために，一般の系に対してエネルギーギャップ ∆t や行列要素 h1(t)|∂t H(t)|0(t)i の振る舞いを評価して Γ(t) の微分方程式を解いてみる。そうすると，t が 0 から十分大きな値に変化する
につれて，Γ(t) が t のべきで減衰する (Γ(t) ∝ t−c ) と，ほぼ確実に正しい解に行き着くことが示される [9, 10, 13]。
一方，古典的な汎用最適化アルゴリズムであるシミュレーテッド・アニーリングでは，一般に，
温度 T を時間 t の関数として T (t) ∝ a/ log t のように減衰させると十分な長時間の後に最適解に到達することが証明されている [11] 4 。Γ(t) ∝ t−c と T (t) ∝ a/ log t を比べると，前者の方が制御変数をより高速に減少させている。この意味において，量子アニーリングがシミュレーテッド・
アニーリングより高速なアルゴリズムであることが言える。以上の論理の詳細を [12, 13] に従って見てみよう。

4.1

量子アニーリングの収束条件

最適化問題が，一般には多体相互作用を含むイジング模型の基底状態探索問題として定式化出来たとする。

H0 = −N∑i=1

Ji σiz −∑ij

Jij σiz σjz −∑Jijk σiz σjz σkz − · · ·ijk

4 温度と計算時間を同じ T で表してあるので，混同しないようにしていただきたい。

(11)

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

このハミルトニアンは全体として N に比例する示量性の条件を満たしているとする。これに量子ゆらぎを表す横磁場の項

H1 = −

N
∑

σix ,

(12)

i=1

を強さ Γ(t) で付け加えて，全ハミルトニアンをH(t) = H0 + Γ(t)H1(13)とする。Γ(t) は t = 0 で非常に大きな値に取り，t が増加するにつれて単調に減少していって

t → ∞ では 0 に収束する関数であるとする。状態ベクトルは，初期ハミルトニアン (12) の自明な基底状態から最終的な (11) の非自明な基底状態にシュレディンガー方程式に従って自然に発展すると期待される。Γ(t) の関数形をどのように選んで減少させていくと各瞬間において (13) の定常状態を近似的にたどっていくことが出来るかが問題になる。次の定理はその解答である。
定理 4.1. (13) で記述される系の断熱条件が満たされるための十分条件は次の式で与えらえる。
ある正の t0 が存在して，t > t0 なる t において

Γ(t) = a(δt + c)−1/(2N −1)(14)

ここで a と c は O(N 0 ) の量であり，δ は断熱条件の式 (9) に現れる微小量である。
証明には，次のホップによる定理 [15] を使う。
定理 4.2. 正方行列 M のすべての要素が正ならば (Mij > 0)，最大固有値 λ0 およびそれ以外のすべての固有値 λ は次の不等式を満たす。

|λ| ≤

κ−1
λ0
κ+1

(15)

ここで，κ は次のように定義される。

κ ≡ maxMiki,j,k Mjk

(16)

定理 4.1 の証明. 式 (14) が断熱条件 (9) を満たすことを示す。このためにはエネルギーギャップとハミルトニアンの時間微分の行列要素を評価しなければならない。後者は簡単である。H(t) の時間依存性は Γ(t) だけにしかないことと H1 は N 項の和であることから¯¯¯¯¯hj(t)| dH(t) |0(t)i¯ ≤ −N dΓ(t)
¯
¯
dt
dt

(17)

dΓ/dt は負であることに注意。
エネルギーギャップの評価をするために，定理 4.2 を行列 M ≡ (E+ − H(t))N に適用する。ここで，E+ は Γ0 ≡ Γ(t0 ) として E+ > Emax + Γ0 を満たす定数である。Emax は H0 の最大固有値である。E+ − H(t) は非負かつ既約（どの状態間もせいぜい N ステップで遷移可能）だから，

M のすべての行列要素は {σiz } を対角化する表現で正定値である。
Γ(t) < Γ0 を満たす t > t0 に対して，E+ − H(t) のすべての対角要素は，0 でないすべての非対角要素 Γ(t) より大きい。これより，M の最小要素は完全に反転した 2 つの状態間で与えられその具体形は N !Γ(t)N である。N ! はスピンを反転する順序の数に対応する。Γ(t)H1 を −N Γ0 で

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

置き換えると，M の最大行列要素は (E+ − Emin + N Γ0 )N であることが分かる。Emin は H0 の最小固有値である。こうして，次の式が成立することが示された。

κ≤(E+ − Emin + N Γ0 )N
N !Γ(t)N(18)

H(t) の固有値を εj (t) と書くと，(15) は次のように書き換えられる。
[E+ − εj (t)]

N

≤

κ−1
N
[E+ − ε0 (t)]
κ+1

(19)

(18) を上式に入れると∆j (t) = ²j (t) − ²0 (t) ≥2[E+ − ε0 (t)]N !
Γ(t)N ≡ AΓ(t)N
N (E+ − Emin + N Γ0 )N(20)

ここで， κ ≥ 1 かつ N ≥ 1 のとき 1 − ((κ − 1)/(κ + 1))1/N ≥ 2/N (κ + 1) であることを使った。
係数 A をスターリングの公式を使って評価する。
√)N(]
2 2πN [E+ − εmax
N
0
A≈
N eN
E+ − Emin + N Γ0

(21)

を maxt>t0 {ε0 (t)} と定義した。A が N の増加と共に指数関数的に減少することが分かる。
εmax
0(17) と (20) により，t > t0 での収束条件が

−

N
dΓ(t)
=δ¿1
A2 Γ(t)2N dt(22)であることが分かる。δ は任意の微小量。これを積分すると (14) になる。

4.2

シミュレーテッド・アニーリングの収束条件

今度は古典的なシミュレーテッド・アニーリングの収束条件を検討しよう。面白いことに，古典的な問題をわざわざ量子力学に焼き直すと顕著な結果が得られるのである [10]。
前節と同様に，イジング模型の基底状態探索として定式化された最適化問題を考える。

H=−

N
∑

Ji σiz −∑i=1

Jij σiz σjz −ij∑

Jijk σiz σjz σkz − · · ·

(23)

ijk

古典的な物理量 Q({σiz }) の統計力学的な期待値は

hQiT =1 ∑ −βHe
Q({σi })
Z(T )

(24)

{σ}

である。和は，パウリ行列の z 成分のすべての固有値 σiz = σi (±1) (∀i) について取ってある。{σi }は {σ1 , σ2 , · · · , σN } である。次の定理が中心的な役割を果たす。
定理 4.3. 統計力学的期待値 (24) は，同じ Q の次の波動関数による量子力学的な期待値と一致する。

|ψ(T )i = e−βH/2∑{σ}|{σi }i(25)

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

ここで，|{σi }i は各 σiz を対角化した表現における固有値 σi を持つ固有状態の集合である.
さらに，T > 0 として，波動関数 (25) は次の量子力学的ハミルトニアンの基底状態である。
∑∑
Hqj (T ) ≡ −χ(σxj − eβHj )(26)
Hq (T ) = −χjj

ここで，Hj は (23) において j 番目のサイトを含む項である。
∑∑
Jjkl σjz σkz σlz − · · ·Hj = −Jj σjz −Jjk σjz σkz −
k

(27)

kl

係数 χ は，p = maxj |Hj | として χ = e−βp で定義される。
証明. 前半は自明。

1 ∑ −βHhψ(T )|Q|ψ(T )i=eh{σi }|Q|{σi }i = hQiThψ(T )|ψ(T )i
Z(T )

(28)

{σ}

後半を示すためにまず次式に注目する。
∑∑
σxj|{σi }i =|{σi }i
{σ}

(29)

{σ}

演算子 σxj は単に和の順序を変えるだけだからこの式が成立する。次の関係も有用である。

σxj e−βH/2 = eβHj e−βH/2 σxj(30)

これを示すには，H および Hj が共に対角表示されていることと H − Hj が σjz を含まないことから [H − Hj , σjx ] = 0 に注意して，

σxj e−βH/2 σxj = e−β(H−Hj )/2 σxj e−βHj /2 σxj = e−β(H−Hj )/2 eβHj /2 = eβHj e−βH/2が成立することを使えばよい。以上より

Hqj (T )|ψ(T )i = (σxj − eβHj ) e−βH/2∑|{σi }i = 0

(31)

(32)

{σ}

が導かれる。したがって |ψ(T )i は Hq (T ) の固有値 0 の固有ベクトルである。今使っている行列表現では，−Hq (T ) の 0 でない非対角要素はすべて正であり |ψ(T )i の係数もすべて正だから，ペロン・フロベニウスの定理 [14] によると |ψ(T )i は Hq (T ) の基底状態である。
さて，高温極限では横磁場のみが残り

Hq (T → ∞) = −∑(σxj − 1)(33)jとなる。対応する基底状態波動関数 |ψ(T → ∞)i はすべての許される状態の単純な和になっているが，これは古典的な高温極限の描像とよく一致している。また，低温極限では∑
Hq (T ≈ 0) → χ
eβHj

(34)

j

であるが，この基底状態波動関数は (25) から明らかなように全体のハミルトニアン H の基底状態でもある。こうして古典的な熱ゆらぎを量子力学的な波動関数における量子ゆらぎに対応させることが出来た。
以上の古典・量子対応を使ってシミュレーテッド・アニーリングの収束条件を調べよう。温度を時間の関数 T (t) として単調減少させるとする。

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

定理 4.4. ハミルトニアン Hq (T ) に対する断熱条件は，N が大きい極限で T (t) についての以下の式に帰着される。

T (t) =pN
log(αt + 1)(35)

係数 α は N とともに指数関数的に減少する量である。
次の補題をまず議論しておこう。
補題 4.5. Hq (T ) の基底状態と第 1 励起状態間のエネルギーギャップ ∆(T ) は N が十分大きいとき次の下限を持つ。

√∆(T ) ≥ a N e−(βp+c)N(36)

a と c は N によらない正の数である。
証明. 前節の量子アニーリングの議論は，Γ(t) を χ = e−βp で置き換え ε0 (t) = 0 とすれば今の問題に適用できる。後者の条件は Hq (T )|ψ(T )i = 0 による。非対角要素 χ は対角要素に適当な正の数を足せば常に小さくできるから，条件 Γ(t) < Γ0 (t > t0 ) は不要である。(20) によれば

∆j (t) ≥ Ae−βpN(37)であり，(21) を使うと A が次の式で与えられることが分かる。

√
A ≈ b 2πN e−cN(38)

b と c は O(N 0 ) の定数である。
補題 4.6. Hq (T ) の温度微分は次の式を満たす。

hψ1 (T )|∂T Hq (T )|ψ(T )i = −∆(T )hψ1 (T )|H|ψ(T )i
2kB T 2(39)

ここで，ψ1 (T ) は Hq (T ) の第 1 励起状態である。
証明. すでに示した関係式

Hq (T )|ψ(T )i = 0

(40)

)
(
)
∂
1
∂
Hq (T ) |ψ(T )i = −Hq (T )|ψ(T )i = Hq (T ) −H |ψ(T )i∂T∂T
2kB T 2(41)を微分すると(Hq (T ) の基底エネルギーが 0 であることより Hq (T )|ψ1 (T )i = ∆(T )|ψ1 (T )i が成立することに注目すれば，上式で証明が完了することが分かる。
補題 4.7. H の行列要素は次式を満たす。

|hψ1 (T )|H|ψ(T )i| ≤ pN
証明. H =

∑

√
Z(T )

(42)
√
Z(T ) はj Hj は N の項の和で書かれており，各項の大きさは p で押さえられる。

|ψ(T )i の規格化から来る。

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

定理 4.4 の証明. 量子系 Hq (T ) に対する断熱条件は，δ を十分小さい正の数として
¯
¯
¯
1
dT ¯¯
¯
√
hψ1 (T )|∂T Hq (T )|ψ(T )i=δ
dt ¯∆(T )2 Z(T ) ¯(43)

行列要素に補題 4.6 の結果を使えば，左辺は¯ ¯|hψ1 (T )|H|ψ(T )i| ¯¯ dT ¯¯√
2kB T 2 ∆(T ) Z(T ) ¯ dt ¯

(44)

¯ ¯
¯ dT ¯pN¯ ¯ = δ̃ ¿ 1
2
2kB T ∆(T ) ¯ dt ¯(45)

分母を補題 4.7 の結果で置き換えてが断熱条件が満たされるための十分条件であることが分かる。補題 4.5 を使い，dT /dt < 0 に注意して上式を積分すれば定理 4.4 に至る。

実際の例

5

一般論はこのくらいにして，具体的な問題に適用した例を紹介しよう。

5.1

データベース探索問題(1) の H0 と (5) の H1 から構成されるハミルトニアン (3) を使って，データベース探索問題を解くことにする [16]。このハミルトニアンは単純な形をしており，すべてのエネルギー固有値を求めることが出来る。計算の詳細は付録 B で述べることにして，結果は N − 2 重に縮退した固有値 1 と，縮退のない 1 より小さい固有値 2 つ（基底状態と第 1 励起状態）があることが分かる。
基底状態と第一励起状態との間のエネルギーギャップは
√
(
) (
)
1
t
t
∆t = 1 − 4 1 −
1−
N T
T

(46)

√
である (図 7)。このエネルギーギャップは t = T /2 で最小値 min ∆T = 1/ N を取る。断熱条件
D
0.4

0.2

0.2

0.4

0.6

0.8

1.0

s

-0.2

-0.4

図 7: N = 64 のときの基底エネルギーと第 1 励起エネルギー。横軸は s = t/T 。

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）(9) の分子の行列要素は

¯
¯
¯
¯
dH
max ¯¯h1(t)||0(t)i¯¯ ≤ T(47)
dtと押さえられるから，(10) より計算時間は N/δ に比例する。これは，ひとつずつ順番に探していく古典的な探索と同じ計算量である。
ところで，エネルギーギャップが小さいときに励起状態への遷移が起きやすいことを考えると，

t = 0 から t = T までずっと同じ速さでハミルトニアンを変化させていくのではなく，ギャップが小さいところだけでゆっくり動かせば効率が上がるものと思われる。そこで，一様な変化を表す(3) の代わりに，初期条件 f (0) = 0 および終端条件 f (T ) = 1 を満たす単調増加関数 f (t) を導入
してH(t) = f (t)H0 + (1 − f (t))H1(48)としてみる。そして，0 ≤ t ≤ T 全体での最悪評価 (9) の代わりに，各時刻で断熱条件を課して|h1(t)|∂t H(t)|0(t)i|
=δ
∆2t

(49)

とする。この式は f (t) に対する微分方程式を与えるから，それを解くことにより各時刻での最適な制御法が決まる。断熱条件と行列要素の上限 (47)，およびギャップの式 (46) を組み合わせると，

f (t) が次の式で決まる関数よりもゆっくり変化すればよいことが分かる。
(
( )(
( )))
df
N −1
t
t
2
= 2δ∆t = 2δ 1 − 4
f
1−f
dt
N
T
T
これを解くと

t=

{
}
√
√
1
N
√
tan−1 ( N − 1(2f (t) − 1)) + tan−1 N − 1
2δ N − 1

(50)

(51)

図 8 のように，ギャップが小さいところではゆっくり動かすというもっともな結論になっている。

f (t = T ) = 1 とおくと最終時刻 T の N 依存性が出る。
√
1
π√
N
T = √
tan−1 N − 1 ≈N(52)
δ N −12δ
√となるから計算時間 T は N に比例する。古典的な計算量 N に比べて小さくなるのである。グローバーのアルゴリズムとして知られる，量子計算の伝統的な枠組みでのデータベース探索アル
ゴリズムの計算量と同等になったのである。

5.2

数値計算や実験の例

具体的な最適化問題におけるシミュレーテッド・アニーリングの比較が多数報告されており，ほとんどすべての場合にシミュレーテッド・アニーリングに比べて高速になることが知られている。
本講義ノートには具体例は示さないが，講義では多数紹介する [17]–[31]。なぜこれだけうまくい
くのか直観的にどうもよく分からないが，ほとんどすべての例において熱ゆらぎを用いるより量子ゆらぎを使う方が高速（あるいは確実）な状態探索が出来るのである。
量子アニーリングの実験的な検証も試みられている。ランダムさを持つ磁性体で低温の状態に達するのに，量子アニーリングに相当する磁場・温度制御のほうがシミュレーテッド・アニーリ
ングに相当する方法より確実に平衡状態に行き着けることが示されている [32, 33]。また，微小な SQUID で量子ビットを表現し，それらを適切に結合することによって多数の量子ビットの結合系を実現して量子アニーリングを実行したという報告もされている [34, 35]。

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

1.0

0.8

0.6

0.4

0.2

0.0

0

10

20

30

40

50

60

70

図 8: 時間変化の関数 (51) の振る舞い。N = 30, T = 78 の場合。横軸は t，縦軸は f (t)。

6

おわりに量子力学的なゆらぎを用いた状態探索により最適化問題の解を高速に見つけるための汎用アルゴリズムである量子アニーリングの初歩的な解説をした。量子アニーリングは今や大きな分野に育っており，全貌を紹介することはなかなか出来ない。特に，量子断熱計算としての観点からの計算量の問題には全く触れることが出来なかった。量子系の基底状態で系のパラメータを変化させたときに生じる相転移（量子相転移）と計算量の問題が盛んに研究されており，相転移が 1 次のときには相転移点を安定的に通過する際に要する計算量が系の大きさに対して指数関数的に増大することが多い。つまり解決が困難な問題になる。このため，基底状態をたどるというやり方
をする量子アニーリングの過程において 1 次転移が起きるかどうかが重要なポイントとなる。この点についての研究も余裕があれば講義で紹介したい [36, 37, 38]。
量子アニーリングの数理的な研究が量子計算機の構築や量子力学の原理のさらなる解明にどのような寄与をしていくのか，今のところ明らかとは言えない。この分野の研究は未だ海図のない航海であり，それだけに未知の大陸がどこかで待っている期待は大きい。若い人たちの参加を期待したい。

付録 A. 断熱条件
ハミルトニアンの時間変化が s = t/T (t : 0 → T ) の関数で表される (3) のような場合を念頭に置いて断熱条件 (9) を導こう。[13] にしたがって話を進める。各瞬間における定常状態の固有値，
固有ベクトルを ²k (s), |k(s)i と書くことにする。

H(s)|k(s)i = ²k (s)|k(s)i (²0 (s) < ²1 (s) ≤ ²2 (s) ≤ · · · )(53)

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

この式の両辺を s について微分すると，j 6= k のとき

hj(s)|

1dH(s)d|k(s)i =hj(s)||k(s)ids²k (s) − ²j (s)
ds

(54)

が導かれる。j = k のときには

hk(s)|

d|k(s)i = 0
ds

(55)

と取ることが出来る。というのは，|k̃(s)i = eiθ(s) |k(s)i と定義して微分すると

hk̃(s)|dθdd|k̃(s)i = i + hk(s)| |k(s)i
ds
ds
ds

だが，右辺第 2 項は純虚数である。
(
)∗
d
d
d
hk(s)| |k(s)i + hk(s)| |k(s)i =hk(s)|k(s)i = 0
ds
ds
ds

(56)

(57)

よって (56) より θ を適切に選べば (55) が成立する。
さて，時間依存のシュレディンガー方程式の解を |j(s)i で展開しよう。
(
)
∫ s
∑
cj (s)e−iT φj (s) |j(s)iφj (s) =du²j (u)|Ψ(s)i =

(58)

0

j

これをシュレディンガー方程式に入れて (54) と (55) を使うと cj (s) についての微分方程式が得ら
れる。

∑
dcj
eiT (φj (s)−φk (s))dH(s)=
cj (s)
hj(s)||k(s)ids²j (s) − ²k (s)
ds

(59)

k6=j

積分して

cj (s) = cj (0) +

∑∫ s
k6=j

duck (u)0eiT (φj (u)−φk (u))dH(u)hj(u)||k(u)i²j (u) − ²k (u)
du

(60)

初期条件として c0 (0) = 1, cj6=0 (0) = 0 とする。最初は完全に基底状態にあるとするのである。T À 1 で上式の右辺第 2 項は O(T −1 ) になる（下の (64) のあたりの議論を見よ）。よってcj (s) (j 6= 0) は O(T −1 ) であり，上式右辺第 2 項の和で k = 0 の項のみが T À 1 で主要項として残る。したがって j 6= 0 として

c0 (s) ≈ 1∫ seiT (φj (u)−φ0 (u))dH(u)
ducj (s) ≈hj(u)||0(u)i
∆
(u)
du
j
0
)
i (
≈
Aj (0) − eiT (φj (s)−φ0 (s)) Aj (s)
T

(61)
(62)
(63)

ここで

Aj (s) =1dH(s)hj(s)||0(s)i∆j (s)2ds(64)(62) と (63) が T À 1 で等しいことを検証するには，後者を微分して前者と一致することをみるとよい。|cj (s)| ¿ 1 が断熱条件である。s = t/T に注意し，さらに 0 ≤ t ≤ T で (64) の最大値を取ると (9) が得られる。

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

付録 B. データベース探索問題の固有値
ハミルトニアン

(
)
t
s=
T

H(t) = sH0 + (1 − s)H1(65)の定常状態固有値を求める。H0 と H1 はそれぞれ以下で定義されている。
N −1

H0 = I − |0ih0|,

H1 = I − |Ψ0 ihΨ0 |,1 ∑|ji|Ψ0 i = √
N j=0

(66)

一般性を失わずに，探し出すべき項目の番号を 0 とした。
まず，固有値 1 の固有ベクトルが次の式で与えられる。
N −1

1 ∑ i2πjk/(N −1)|Ψk i = √
e
|ji
N j=1

(k = 1, 2, · · · , N − 2)(67)

これを確かめるために，H0 と H1 に作用させると hj|0i = 0 (j 6= 0) および hΨ0 |Ψk i = 0 より5

H0 |Ψk i = |Ψk i,

H1 |Ψk i = |Ψk i

(68)

これより

H(t)|Ψk i = |Ψk i(69)

固有値 1 は N − 2 重に縮退している。次に|Ψi = a|Ψ0 i + b|0i(70)が固有ベクトルであるための条件を求めよう。H0 と H1 に作用させるとb
H0 |Ψi = a|Ψ0 i − √ |0i,NよりH(t)|Ψi =b
H1 |Ψi = b|0i − √ |Ψ0 i
N

(
)
(
)
b(1 − s)assa − √|Ψ0 i + (1 − s)b − √
|0i
N
N

(71)

(72)

これが E|Ψi に等しいという条件

Ea = sa −b(1 − s)√,N

を E について解くと

√
1 1
E= ±
2 2

as
Eb = (1 − s)b − √N()
1
1−4 1−s(1 − s)
N

(73)

(74)

が得られる。従ってエネルギーギャップは√∆=で与えられる。
5 後者は

PN −1
j=1

ei2πjk/(N −1) = 0 による。

()
1
1−4 1−s(1 − s)N(75)

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

付録 C. 古典スピン系のシミュレーションの基礎
シミュレーテッド・アニーリングは古典系のモンテカルロ・シミュレーションに基づいている。
この分野の基礎知識を持たない人のために，文献 [39] に基づいて初歩的な解説を書いておく。

マスター方程式
モンテカルロ・シミュレーションは，マスター方程式で記述される確率的ダイナミクスを数値計算で実現する方法である。イジング模型の言葉を使ってその考え方を解説しよう。
イジング模型は古典的な離散変数であるイジングスピンで表現されており，その状態の時間発
展はニュートン方程式やシュレディンガー方程式によって直接的に記述されるわけではない。系
が熱浴との相互作用により状態を変化させていると考えて，各時刻における状態の確率を記述する関数の時間変化を導入することにより，自然なダイナミクスを導入することが出来る。
スピンの状態を a = {1, −1, −1, · · · , 1} のようにアルファベットで表すことにしよう。時刻 t において系が状態 a にある確率を P (a, t) と書くことにする。N 個のイジングスピンからなる系の状態総数は 2N だから，a は 2N 通りある。この個数の状態について P (a, t) の値が分かれば，時
刻 t における系の確率的な記述が完成する。
短い時間 ∆t の間に状態が a から b に変化する遷移確率が w(a → b)∆t (≥ 0) で与えられるとしよう。このとき，系が状態 a にある確率は w(a → b)∆t · P (a, t) だけ減少する。同様に，状態

b にあったとすると状態 a の確率は w(b → a)∆t · P (b, t) だけ増加する。したがって，状態 a にある確率の変化は次式で表される。
P (a, t + ∆t) − P (a, t) = −∑w(a → b)P (a, t)∆t +b(6=a)∑

w(b → a)P (b, t)∆t

(76)b(6=a)

これがマスター方程式（の離散時間版）である。連続時間版は ∆t → として差分を微分にすれば導かれるが，ここではモンテカルロ・シミュレーションを念頭に置いて，離散時間版で話を進める。
マスター方程式を次の形に書き換えると便利である。

P (a, t + ∆t) =
≡

∑
∑
)
(
w(b → a)∆t · P (b, t)w(a → b)∆t P (a, t) +1−∑b(6=a)b(6=a)Lab P (b, t)

(77)

b

ここで

Laa = 1 −∑w(a → b)∆t, Lab = w(b → a)∆t (b 6= a)(78)b(6=a)と定義した。行列とベクトルで書くと

P (t + ∆t) = LP (t).

(79)

P (t) の a 番目の成分を P (a, t) と考えるのである。L を確率行列と呼ぶことにする。確率行列の右固有ベクトルと固有値を次のように書こう。
Leα = λα eα(80)

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

固有ベクトルが完全系をなすとして P (t) を展開する。

P (t) =∑cα (t)eα(81)α

このとき，時間発展の式 (79) はP (t + ∆t) =∑cα (t)Leα =α∑cα (t)λα eα

(82)

α

t = 0 における初期状態から n ステップ経過した後の確率はP (n∆t) =∑cα (0)λnα eα

(83)

α

と表される。こうして，確率の時間変化は確率行列の固有値と固有ベクトルで記述される。
∑
さて，確率の保存の条件 a P (a, t) = 1 は固有値に強い制約を課す。(83) の左辺のすべての要素の和は 1 だから，右辺も同様に∑

cα (0)λnαα∑eα (a) = 1

(84)

a

を任意の n に対して満たす。ここで eα (a) は eα の a 成分である。(84) の右辺は n によらないか∑
ら，左辺も n によらない。このためには最大固有値 λ0 が 1 で a e0 (a) = 1 を満たし，そのほか∑の固有値は a eα (a) = 0 (α 6= 0) を満たさなければならない。さらに，P (a, n∆t) が n とともに
限りなく増大しないという条件から，α 6= 0 に対して c0 (0) = 1 と |λα | < 1 が必要である。

L の行列要素はすべて非負で6 確率行列 L が既約（どんな状態も任意の別の状態から到達できる）ゆえ，ペロン・フロベニウスの定理により最大固有値 λ0 は縮退してない。以上の議論より，n → ∞ で系は平衡状態 e0 に指数関数的に近づく。
P (n∆t) = e0 + c1 (0)e−n| log λ1 | e1 + · · ·(85)

第 2 固有値を λ1 (< 1) として，緩和時間は 1/| log λ1 | である。こうして，マスター方程式によって記述されるダイナミクスが系を平衡状態 e0 に緩和させることが保証された。
以上の議論は極めて一般的であり，平衡状態がボルツマン分布 e−βH /Z であることは保証されてない。この後者の条件ついて次に検討しよう。

モンテカルロ・シミュレーション
マスター方程式に従って系が時間発展していき，平衡状態が実現して P (a, t) = Peq (a) になったとしよう。(76) の左辺は 0 ゆえ，∑b(6=a)

w(a → b)Peq (a) =∑w(b → a)Peq (b)(86)

b(6=a)

これが成立するための十分条件は各項が等しくなることである。

w(a → b)Peq (a) = w(b → a)Peq (b)6 (78) に現れる対角要素 L
aa は状態 a に留まる確率だから非負。

(87)

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

これが詳細つり合いの条件である。ボルツマン分布 Peq (a) を平衡分布に使うとw(a → b)= e−β(H(b)−H(a))
w(b → a)(88)

これを満たすように遷移確率を選ぶ。よく使われるのは熱浴法とメトロポリス法の 2 種である。
熱浴法ではw(a → b) =e−βH(b)e−βH(a) + e−βH(b)(89)

w(a → b) = e−β(H(b)−H(a))+(90)とし，メトロポリス法ではと選ぶ。f ≥ 0 なら (f )+ = f でありそれ以外では (f )+ = 0 である。
問：熱浴法とメトロポリス法が詳細つり合いの条件を満たすことを確かめよ。
次に，w(a → b) > 0 となる a と b の組み合わせを指定する。一度にひとつのスピンだけを反転させるシングルスピンフリップ法が一番よく使われる。サイト i を選んで確率 w(a → b) によりSi を −Si とするのである。ただしa = {S1 , S2 , · · · , Si , · · · , SN },b = {S1 , S2 , · · · , −Si , · · · , SN }(91)

他のすべての a と b の組み合わせは 0 とする。シングルスピンフリップ法により確率行列が既約になることは明らかである。シングルスピンフリップ法によると，熱浴法とメトロポリス法の遷移確率は次のようになる。

w(Si → −Si ) =e−βH(−Si )e−βH(Si ) + e−βH(−Si )=e−β∆Ei1 + e−β∆Ei

(92)

および

w(Si → −Si ) = e−β(∆Ei )+(93)

ここで ∆Ei = H(−Si ) − H(Si ) であり，H の Si 以外の変数への依存性は省略した。∆Ei がサイ
ト i に関わる局所的な変数だけで書かれるとき，(92) と (93) は実用的に有用である。ハミルトニ
アンが次のように書かれるとしよう。

H = Hi + H 0 = −∑Jij Si Sj + H 0

(94)

j

ここで H 0 は Si を含まない部分である。このとき，シングルスピンフリップによるエネルギーの変化は∆Ei = H(−Si ) − H(Si ) = 2∑

Jij Si Sj(95)j

このエネルギー変化量は i の周りの局所変数だけで書かれており，数値計算において容易に計算できる。これを (92) や (93) に代入して簡単に遷移確率を求められる。
実際の手順は次の通り。∆t = 1 とする。

1. スピン配位 {Si } を初期化する。完全にランダムにしたり，完全にそろった強磁性状態にすることが多い。

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

2. シングルスピンフリップ試行をさせるサイト i を選ぶ。
3. エネルギー変化量 ∆Ei を計算する。
4. 0 と 1 の間の一様乱数 r を発生させ，w(Si → −Si ) と比較する。r < w(Si → −Si ) ならば
Si → −Si を実現し r > w(Si → −Si ) なら何もしない。これにより確率 w(Si → −Si ) でフリップが実現される。
5. スピン配位 {Si } を使って物理量を計算する。
6. 統計精度が上がるまで 2 から 5 までを繰り返す。
いくつか注意点を述べておこう。緩和時間 τ = 1/| log λ1 | が経つまでは平衡状態に達しきってないので，初期条件の影響を除くために最初の適当な数のステップ数については平衡物理量の測定は避けなければならない。また，平衡状態に達した後もボルツマン分布の統計的なサンプリングをしているのだから統計誤差の（場合によっては系統誤差も）影響を十分注意して評価しなければならない。有限サイズ効果についても同様である。
シミュレーテッド・アニーリングにおいては温度をステップ数の関数として下げていく。完全に収束するためには，本文で議論したように十分長い時間かけて十分ゆっくり下げて行く必要があるが，実際には結構さっと下げても実用上は差し支えない解に到達することも多い。

プログラム例‘genrand()’ は 0 と 1 の間の一様乱数を発生させる関数。自分で用意する。
/* ---------------------------------------------------------------強磁性的相互作用 J=1 を持つイジング模型のメトロポリス法によるシミュレーション------------------------------------------------------------------*/#include <stdio.h>#include <stdlib.h>#include <math.h>#include <time.h>#define ls 40// linear size of the system

int main(void){
int mcs=100000;
int discard=500;
int measure=10;// total MC steps (per spin)// # of steps to remove initial effects// measurement interval

int i,j;int i1,i2;
int mcprocess;// variables to control loops// site index for flip trial// current Monte Carlo step

int spin[ls][ls];int ip[ls],im[ls];// spin configuration// table of right and left neighbours

double t=2.0;// temperature《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

double delta_E,energy;// energy change and energy

double ecurrent;
double field;// local energy// local field

double denominator;// measurement normalizer
double genrand(void);// random number generator//--------------------Initialization----------------------// Table of nearest neighbour sites for periodic boundaries
for (i=0;i<ls;i++){ip[i]=i+1;
im[i]=i-1;// right (upper) neighbour// left (lower) neighbour}
ip[ls-1]=0;
im[0]=ls-1;// right (up) of ls-1 is 0.
// left (bottom) of 0 is ls-1.

// Initial configuration.
for (i=0;i<ls;i++){for (j=0;j<ls;j++){All up.

spin[i][j]=1;}}
srand((unsigned)time(NULL));// random number initialization//--------------------- Main loop ------------------------energy=0.0;
for (mcprocess=1;mcprocess<=mcs*ls*ls;mcprocess++){i1=(int)(genrand()*ls);// Choose a site for flip trial. x coordinate.
i2=(int)(genrand()*ls);// Same for y coordinate.
field=spin[ip[i1]][i2]+spin[i1][ip[i2]]+spin[im[i1]][i2]+spin[i1][im[i2]];// Sum of spin states around spin[i1][i2]delta_E=field*spin[i1][i2]*2.0;// Energy change caused by the flip of spin[i1][i2].
// Execution of the Metropolis method----------------------if(delta_E<0){spin[i1][i2]*=-1; // Flip if energy decreases.
}
else{
if(exp(-delta_E/t)>genrand()) spin[i1][i2]*=-1;} // Probabilistic flip for energy increase.
//---------------------Measurement -----------------------if(mcprocess%(measure*ls*ls)==0){
if(mcprocess>(discard*ls*ls))//Skip until the system equilibrates.
{《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）

ecurrent=0.0;
for (i=0;i<ls;i++){for (j=0;j<ls;j++){field=spin[ip[i]][j]+spin[i][ip[j]]+spin[im[i]][j]+spin[i][im[j]];ecurrent+=spin[i][j]*field;}}
ecurrent/=-2.0;//Divide by two since each bond is counted twice.
energy+=ecurrent;}}//-----------------End of the Monte Carlo loop -----------}//------------------------Avereage------------------------denominator=(mcs-discard)/(double)measure*ls*ls;//Number of data points divided by the system sizeenergy/=denominator;// Simple average approximates the canonical average.
//-------------------------Output-------------------------printf("Average energy per spin :%f\n",energy);
return 0;}

参考文献[1] 西森秀稔: 相転移・臨界現象の統計物理学 (培風館, 2005)[2] 西森秀稔: スピングラス理論と情報統計力学（岩波, 1999）[3] T. Kadowaki and H. Nishimori: Phys. Rev. E 58, 5355 (1998).
[4] A. Das and B. K. Chakrabarti: ‘Quantum Annealing and Related Optimization Methods’(Springer) (2005).
[5] G. E. Santoro and E. Tosatti: J. Phys. A 39, R393 (2006).
[6] A. Das and B. K. Chakrabarti: Rev. Mod. Phys. 80, 1061 (2008).
[7] M.A. Nielsen and I. L. Chuang: Quantum Computation and Quantum Information (Cambridge, 2000).
[8] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren and D. Preda: Science 292,
5516 (2001).
[9] S. Morita and H. Nishimori: J. Phys. Soc. Jpn. 76, 064002 (2007).

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）[10] R. D. Somma, C. D. Batista, and G. Ortiz: Phys. Rev. Lett. 99, 030603 (2007).
[11] S. Geman and D. Geman: IEEE Trans. Pattern Anal. Mech. Intell. PAMI-6, 721 (1984).
[12] S. Morita and H. Nishimori: J. Math. Phys. 49, 125210 (2008).
[13] 森田悟史: 東京工業大学 学位論文 (2008).
[14] E. Seneta: Non-negative Matrices and Markov Chains (Springer, New York, 2006)[15] E. Hopf: J. Math. Mech. 12, 683 (1963).
[16] J. Roland and N. J. Cerf: Phys. Rev. A 65, 042308 (2002).
[17] G. E. Santoro, R. Mortonak, E. Tosatti and R. Car: Science 295, 2427 (2002).
[18] R. Mortonak, G. E. Santoro and E. Tosatti: Phys. Rev. B 66, 094203 (2002).
[19] R. Mortonak, G. E. Santoro and E. Tosatti: Phys. Rev. B 70, 057701 (2004).
[20] D. A. Battaglia, G. E. Santoro and E. Tosatti: Phys. Rev. E 71, 066707 (2005)[21] L. Stella, G. E. Santoro and E. Tosatti: Phys. Rev. B 72, 014303 (2005).
[22] L. Stella and G. E. Santoro: Phys. Rev. E 75, 036703 (2007)[23] L. Stella, G. E. Santoro and E. Tosatti: Phys. Rev. B 73, 144302 (2006).
[24] T. Cavena, F. Fazio and G. E. Santoro: Phys. Rev. B 76, 144427 (2007).
[25] M. Sarjala, V. Petaja and M. Alava: J. Stat. Mech, P01008 (2006).
[26] A. Das, B. K. Chakrabarti and R. B. Stinchcombe: Phys. Rev. E 72, 026701 (2005).
[27] S. Suzuki and M. Okada: J. Phys. Soc. Jpn. 74, 1649 (2005).
[28] Y.-H. Lee and B.J. Berne: J. Phys. Chem. 104, 86 (2000).
[29] Y. Matsuda, H. Nishimori and H. G. Katzgraber: New J. Phys. 11, 073021 (2009).
[30] T. Neuhaus, M. Peschina, K. Michielsen, and H. De Raedt: Phys. Rev. A 83, 012309(2011).
[31] J. Inoue: Phys. Rev. E 63, 046114 (2001).
[32] J. Brooke, D. Bitko, T. F. Rosenbaum and G. Aeppli: Science 284, 779 (1999).
[33] H. M. Ronnow1, R. Parthasarathy, J. Jensen, G. Aeppli, T. F. Rosenbaum and D. F.
McMorrow: Science 308, 389 (2005).
[34] Z. Bian et al: arXiv:1201.1842.
[35] M. W. Johnson et al.: Nature 473, 194 (2011).

《講義ノート》

物性研究・電子版 Vol. 3, No. 3, 033203（2014年8月号）[36] T. Jorg, F. Krzakala, J. Kurchan and A. C. Maggs: Phys. Rev. Lett. 101, 147204 (2008).
[37] A. P. Young, S. Knysh and V. N. Smelyanskiy: Phys. Rev. Lett. 104, 020502 (2010).
[38] Y. Seki and H. Nishimori: arXiv:1203.2418. Phys. Rev. E に出版予定.
[39] H. Nishimori and G. Ortiz: Elements of Phase Transitions and Critical Phenomena (Oxford
University Press, 2011) (いわゆる「太陽系最強の相転移理論入門書」)