docs/tutorial/ja/000-intro_optimization_and_Ising.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "1052b600",
"metadata": {},
"source": [
"# 組合せ最適化とイジング模型\n",
"このノートでは、組合せ最適化および、イジング模型を用いた最適化計算の背景について説明します。\n",
"すぐに、OpenJijを動かしてみたい人は、次のチュートリアルから進めても問題ありません。"
]
},
{
"cell_type": "markdown",
"id": "64953e94",
"metadata": {},
"source": [
"## 組合せ最適化とは\n",
"最適化問題とは、問題で与えられた制約条件を満たす解のうち、最も良い解を見つけろという問題です。\n",
"特に、組合せ最適化問題といった場合には、取りうる組み合わせの中から最も良い組合せを見つける問題です。\n",
"\n",
"組合せを見つける問題というものが、どのようなものかを理解するために、いくつか例を見てみましょう。"
]
},
{
"cell_type": "markdown",
"id": "99ae443b",
"metadata": {},
"source": [
"## 組合せ最適化の例\n",
"### ナップサック問題\n",
"ある探検家が、探検の末、$N$個の宝物を見つけたとします。\n",
"ですが、この探検家は一つのナップサックしか持っておらず、このナップサックには最大$W$までしか入れつことができません。\n",
"$N$個の宝物の価値がそれぞれ$c_i\\ (i \\in\\{0,1,\\dots,N-1\\})$で、その大きさが$w_i$であるとします。\n",
"この時、ナップサックに入れることが可能な宝物の組合せのうち、宝物の合計の価値が最大になる組合せを見つけろという問題が**ナップサック問題**です。\n",
"\n",
"ここでは、一つのナップサックだけが存在するときに、そこにどのようにうまく詰めるかを考えましたが、複数のナップサックがあるときに、$N$個の宝物を全て入れるのに必要な最小のナップサックの個数はいくつかという問題は、**ビンパッキング問題**と呼ばれています。\n",
"\n",
"また、類似の問題としては、トラックへの荷物の詰め込みのような、ランダムな大きさの長方形の箱をそれぞれが重ならないように、容器に詰めていく問題は、**長方形詰め込み問題**と呼ばれています。\n",
"こちらに似た問題として無限に大きい箱に球を詰めるときに、適切な球の詰め方を調べる問題は、有名なケプラー予想の問題になっています。\n",
"\n",
"### 巡回セールスマン問題\n",
"あるセールスマンが商品を売るために$N$都市をそれぞれ1回ずつ回って、最初に出発した都市に戻ってくるとします。\n",
"この時、移動にかかるコストを抑えるために、移動距離が最も短い経路を見つけろという問題が**巡回セールスマン問題**です。\n",
"\n",
"ちなみに、セールスマンの各都市における商談時間が決まっており、特定の時間に特定の都市にいなければいけないという制約条件が加わった問題は、**時間枠付き巡回セールスマン問題**と呼ばれています。\n",
"\n",
"また、巡回セールスマン問題では、セールスマンは一人しかいませんが、例えば、**配車計画問題**では、複数の車両が都市を回るという状況を考えます。\n",
"具体的には、物を配送するとき、倉庫から$N$箇所の配送先へ$M$台の車両で配送を行うというような問題です。\n",
"\n",
"### グラフ彩色問題\n",
"ここでは、グラフ彩色問題のうち最も単純な頂点彩色問題を考えます。\n",
"**頂点彩色問題**とは、あるグラフが与えられたときに、隣り合う頂点同士が同じ色にならないように、グラフに存在する全ての頂点の色を塗ることができるかという問題です。\n",
"与えられたグラフが、平面グラフであれば、この問題は「どんな地図も隣接する領域が違う色になるように塗り分けることができる」という有名な四色定理に関連した問題となります。\n",
"\n",
"また、有名なパズルの一つである数独は、$9\\times 9$のマス目の中に$9$この数字を縦、横、そして$3\\times 3$の領域の中に同じ数字が被らないように、適切に入れていくというゲームですが、これは、$9\\times 9$マス目を9色で塗り分ける問題だと思えば、数独もグラフ彩色問題の一種として解くことができます。\n",
"\n",
"以上の例から分かるように、組合せ最適化問題では、組合せや順序を調べる問題が多く、また、実社会に関連した問題も数多く存在します。\n",
"\n",
"このノートでは、イジング模型を用いた組合せ最適化計算について説明しますが、その他、様々な問題に対して、様々なアルゴリズムが組合せ最適化計算の文脈で研究されていますので、そちらに関して興味がある方は、ぜひ教科書等を参考にしてみてください[1]。"
]
},
{
"cell_type": "markdown",
"id": "c65d8b5c",
"metadata": {},
"source": [
"# イジング模型を用いた最適化計算"
]
},
{
"cell_type": "markdown",
"id": "a0b4546d",
"metadata": {},
"source": [
"## イジング模型\n",
"イジング模型は、元々磁性体の振る舞いを理解するための理論模型として導入された模型です。\n",
"イジング模型のハミルトニアン$H_J$は以下のように記述されます。\n",
"\n",
"$$H_J = \\sum_{i>j}J_{ij}\\sigma_i\\sigma_j + \\sum_i h_i \\sigma_i$$\n",
"\n",
"ここで、$\\sigma_i\\in \\{-1,1\\}$は、$i$番目のスピンを表すバイナリ変数です。\n",
"$J_{ij}$はスピン同士の結びつきの強さを示し、$h_i$は各スピンが1か-1のどちらかを取りたがる性質の強さを示しています。\n",
"イジング模型の物理的な意味については、教科書等を参照してください[2]。\n",
"このハミルトニアンは最適化問題におけるコスト関数に対応しています。\n",
"つまり、最適な組み合わせを求めるという組合せ最適化問題を、このイジング模型のハミルトニアンを最小化するスピン配置を求める問題に置き換えて解くのがイジング模型を用いた最適化計算です。"
]
},
{
"cell_type": "markdown",
"id": "d63684e0",
"metadata": {},
"source": [
"## Quadoratic Unconstraint Binary Optimization (QUBO) 定式化\n",
"イジング模型を用いて組合せ最適化問題を解くためには、組合せ最適化問題を$H_J$の形に落とし込む必要があります。\n",
"しかし、通常、数理最適化の文脈では、-1,1の値を取るスピン変数$\\sigma_i$よりも、0,1の値をとるバイナリ変数$q_i$を用いて定式化する方が、素直な場合が多いです。\n",
"このバイナリ変数$q_i$が、スピン変数$\\sigma_i$と\n",
"\n",
"$$q_i = \\frac{\\sigma_i + 1}{2}$$\n",
"\n",
"という関係を持つことを用いて、イジング模型を次の形に書き換えた\n",
"\n",
"$$H_Q = \\sum_{i>j}Q_{ij}q_iq_j + \\sum_i Q_{ii} q_i$$\n",
"\n",
"は、Quadoratic Unconstraint Binary Optimization (QUBO)定式化と呼ばれています。\n",
"スピン変数$\\sigma_i$とバイナリ変数$q_i$は相互変換が可能という意味で、互いに等価な模型です。\n",
"多くのイジング模型を用いた最適化計算の定式化では、このQUBO定式化が用いられています。\n",
"\n",
"### 巡回セールスマン問題のQUBO定式化\n",
"ここでは、イジング模型を用いた最適化の例として、巡回セールスマン問題をQUBOで定式化してみます。\n",
"時刻$t$に都市$i$にいるとき1となり、それ以外の場合では0となるバイナリ変数$q_{it} \\in \\{0,1\\}$を考えます。\n",
"この問題では、都市を回る経路の合計を最小化したいので、都市$i$と$j$の間の距離を$d_{ij}$とかくと、目的関数は\n",
"\n",
"$$\\sum_{ijt}d_{ij}q_{it}q_{j(t+1)} \\mod n$$\n",
"\n",
"と書くことができます。\n",
"ここで、この問題では一人のセールスマンが各都市を一度ずつ回る問題なので、制約条件として、ある時間ではセールスマンは一つの都市にしかいないという制約条件と、ある都市にはセールスマンは一度しか訪れないという、次の2つの制約条件ががつきます。\n",
"\n",
"$$\\sum_i q_{it} = 1,\\ \\forall t\\\\ \\sum_t q_{it} = 1,\\ \\forall i$$\n",
"\n",
"\n",
"ここまでは、ただ数理モデルをバイナリ変数で定式化しただけで、これを$H_Q$の形に落とし込む必要があります。\n",
"制約条件を一つのハミルトニアンにまとめるために、Penalty法と呼ばれる方法がよく用いられます。\n",
"Penalty法とは大雑把には、制約条件を2乗した式を目的関数に加えることで、制約条件の影響を目的関数の中に取り込む方法です。\n",
"なので、Penalty法を用いると、最終的な目的関数は、\n",
"\n",
"$$H = \\sum_{ijt}d_{ij}q_{it}q_{j(t+1)} + \\sum_t A_t(\\sum_i q_{it} - 1)^2 + \\sum_i B_i(\\sum_t q_{it} - 1)^2$$\n",
"\n",
"と書くことができます。この式変形していくことで、最終的に$H_Q$の形にすることができます。\n",
"ここで、$A_t,B_i$は各制約条件に対する重み係数で、これらの係数が大きいと、制約条件を満たしやすくなりますが、目的関数の最小化が無視される傾向があり、一方で、重みを小さくしすぎると、目的関数の最小化は実行されますが、制約条件を満たした解が得づらくなります。\n",
"このように、制約条件を含む問題では、これらの重み係数は適切に設定する必要があります。\n",
"\n",
"ここでは、巡回セールスマン問題について説明しましたが、他にも具体的な最適化問題に対する定式化に関しては、本チュートリアルでも扱うのでそちらも参考にしてください。\n",
"また、様々なNP完全問題をどのように定式化するかについては、[3]のような有名な論文も存在するので、そちらも参考に定式化を行うと良いです。"
]
},
{
"cell_type": "markdown",
"id": "24b62e5f",
"metadata": {},
"source": [
"## 焼きなまし法 (Simulated Annealing)\n",
"ここまでは、どのように組合せ最適化問題をイジング模型として定式化するかを考えました。\n",
"しかし、最も重要なのはどのようにコスト関数を最小化した状態を得るかということです。\n",
"次に、イジング最適化の基礎となる焼きなまし法について説明します。\n",
"焼きなまし法(Simulated Annealing)は、「焼きなまし」という言葉から想像できるように、徐々に物質の温度を下げていくという物理プロセスから着想を得て考案されたメタヒューリスティックアルゴリズムです[4]。\n",
"焼きなまし法には、十分にゆっくり温度を下げれば最適解に辿り着くことができるというGeman-Gemanの定理が存在[5]しますが、実用上は、この理論が要求するような遅さで温度を下げることはできないため、確率的に最適解を含め様々な解を得ることになります。\n",
"\n",
"続いて、具体的なアルゴリズムについて説明していきます。\n",
"あるスピン配置$\\vec{\\sigma}= \\{\\sigma_0,\\sigma_1,\\dots,\\sigma_n \\}$が得られており、そのときのコスト関数の値を$H(\\vec{\\sigma})$と書くことにします。\n",
"このときに、このスピン配置をわずかに変化させて新しいスピン配置$\\vec{\\sigma}_\\mathrm{new}$を得ます。\n",
"どのような方法で、新しいスピン配置を得ても良いのですが、よく用いられる方法としては、一つのスピンを選び、そのスピンを反転させるシングルスピンフリップと呼ばれる手法が用いられています。\n",
"それ以外にも、Swendsen and Wangアルゴリズムと呼ばれる手法も用いられています[6]。\n",
"新しいスピン配置のコスト関数の値$H(\\vec{\\sigma}_\\mathrm{new})$を計算して、少しでもコスト関数の値が下がっていれば、新しいスピン配置を受け入れ、コスト関数の値が上がっている場合には、$\\exp(-(H(\\vec{\\sigma}_\\mathrm{new}) - H(\\vec{\\sigma}))/T)$という確率で新しいスピン配置を受け入れるとします。ここで、$T$は温度に相当するパラメータで、この大きさに依存して、コスト関数の値が上昇するスピン配置が採択される確率が変化します。\n",
"焼きなまし法では、この$T$をゆっくりと下げながら、新しいスピン配置への更新を行なっていくアルゴリズムです。\n",
"\n",
"まとめると、焼きなまし法は次のような次の確率\n",
"\n",
"$$\n",
"P(\\vec{\\sigma}_\\mathrm{new}|\\vec{\\sigma}) = \\begin{cases}\n",
" 1 & (H(\\vec{\\sigma}_\\mathrm{new}) \\leq H(\\vec{\\sigma})) \\\\\n",
" \\exp(-(H(\\vec{\\sigma}_\\mathrm{new}) - H(\\vec{\\sigma}))/T) & (H(\\vec{\\sigma}_\\mathrm{new}) > H(\\vec{\\sigma}))\n",
" \\end{cases}$$\n",
" \n",
"でスピン配置の更新を行いながら、ゆっくりと温度$T$を下げていく手法です。\n",
"このコスト関数の確率的な振る舞いが、物理的には熱揺らぎの効果で、別の状態へと遷移していく振る舞いを表しています。\n",
"\n",
"焼きなまし法は、コスト関数を下げる状態だけでなく、コスト関数が上がる状態も確率的に受け入れることで、局所解にハマるのを防ぐ効果が生まれています。\n",
"コスト関数が上がる状態がどの程度受け入れられるかは、$T$に依存しています。\n",
"ですが、適切な温度範囲は解きたい問題に依存しているため、常に問題ごとに適切な温度範囲の設計が必要となります。\n",
"OpenJijでは、自動的にこの温度を調整する仕組みを導入している。"
]
},
{
"cell_type": "markdown",
"id": "2ff0a008",
"metadata": {},
"source": [
"## 量子アニーリング\n",
"焼きなまし法は、熱揺らぎの効果で局所解にハマるのを避けるアルゴリズムでした。\n",
"非常に小さい、ミクロな世界での物理を記述する量子力学によって記述される世界では、量子揺らぎという熱揺らぎとは別の確率的な効果が存在します。\n",
"熱揺らぎの代わりに、この量子揺らぎの効果を用いて、最適化計算を行おうというのが量子アニーリングです。\n",
"量子アニーリングに用いるハミルトニアンは、\n",
"\n",
"$$H(s) = \\left(\\sum_{i>j}J_{ij}\\sigma^z_i\\sigma^z_j + \\sum_i h_i \\sigma^z_i\\right)-\\Gamma(s)\\sum_i \\sigma_i^x$$\n",
"\n",
"です。\n",
"\n",
"ここで、$s$は一種の時間を表し、$\\Gamma(s)$によって、それぞれの項の影響度合いが調整されています。\n",
"具体的には、通常は初期時刻では、$\\Gamma(s)$が非常に大きくなるように設定されており、徐々に$\\Gamma(s)$を小さくしていき、二つの項の影響を入れ替えるように係数を調整していきます。\n",
"なぜこのような方法をとるかの詳細は、教科書等[7]に譲りますが、2項目は取りうるスピン配置全ての重ね合わせ状態を表しており、1項目は元々の最適化問題のハミルトニアンを表しています。\n",
"なので、量子効果によって生じている全てのスピン配置の重ね合わせ状態から、徐々に最適化したいハミルトニアンの影響を加えていき、最終的に最適化したいハミルトニアンにおける最適値にスピン状態が存在するように調整しているのです。"
]
},
{
"cell_type": "markdown",
"id": "16910aec",
"metadata": {},
"source": [
"## Simulated Quantum Annealing\n",
"量子アニーリングは純粋な量子効果を用いるため、大規模な計算を古典計算機で行うのは難しいです。\n",
"一方で、この量子効果を模した計算手法として、Simulated Quantum Annealingと呼ばれる手法が存在します。\n",
"導出は教科書等を参考にしてください[8]。\n",
"この手法では、$H(s)$に対する量子アニーリングを行う代わりに、\n",
"\n",
"$$\\tilde{H} = \\frac{1}{L}\\left(\\sum_{i>j}\\sum_{l=1}^{L}J_{ij}\\sigma_{il}^z\\sigma_{jl}^z + \\sum_{i}\\sum_{l=1}^{L}h_i\\sigma_{il}^z\\right)\\\\ \\quad - \\frac{T}{2}\\log \\coth\\left(\\frac{\\Gamma}{TL}\\sum_i\\sum_{l=1}^L\\sigma^z_{ik}\\sigma^z_{ik+1}\\right)$$\n",
"\n",
"というハミルトニアンに対して焼きなまし法を行う手法です。\n",
"ここで、$L$はトロッター数と呼ばれる量であり、理論的には$L$が十分に大きいときに、この手法は量子アニーリングに対応した結果を与えます。"
]
},
{
"cell_type": "markdown",
"id": "500dd6ca",
"metadata": {},
"source": [
"## 参考文献\n",
"[1] Korte, Bernhard H., et al. Combinatorial optimization. Vol. 1. Heidelberg: Springer, 2011.\n",
"\n",
"[2] Kardar, Mehran. Statistical physics of fields. Cambridge University Press, 2007.\n",
"\n",
"[3] [Lucas, Andrew. \"Ising formulations of many NP problems.\" Frontiers in physics (2014): 5.](https://www.frontiersin.org/articles/10.3389/fphy.2014.00005/full?ref=https://githubhelp.com)\n",
"\n",
"[4] [Kirkpatrick, Scott, C. Daniel Gelatt Jr, and Mario P. Vecchi. \"Optimization by simulated annealing.\" science 220.4598 (1983): 671-680.](https://www.science.org/doi/abs/10.1126/science.220.4598.671)\n",
"\n",
"[5] [Geman, Stuart, and Donald Geman. \"Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images.\" IEEE Transactions on pattern analysis and machine intelligence 6 (1984): 721-741.](https://ieeexplore.ieee.org/document/4767596)\n",
"\n",
"[6] [Swendsen, Robert H., and Jian-Sheng Wang. \"Nonuniversal critical dynamics in Monte Carlo simulations.\" Physical review letters 58.2 (1987): 86.](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.58.86)\n",
"\n",
"[7] [Kadowaki, Tadashi, and Hidetoshi Nishimori. \"Quantum annealing in the transverse Ising model.\" Physical Review E 58.5 (1998): 5355.](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.58.5355)\n",
"\n",
"\n",
"[8] Tanaka, Shu, Ryo Tamura, and Bikas K. Chakrabarti. Quantum spin glasses, annealing and computation. Cambridge University Press, 2017."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}