はじめに
最近はD-Waveが量子アニーリング方式の量子コンピュータを商用利用し初めて、量子コンピューティング業界を賑わいつつありますね。実際のビジネス面への有用性はとりあえず横に伏せて、その理論的な側面について量子力学の基本的な考え方からできるだけ分かりやす説明したいと思います。
(※数学的な厳密さや内容の正確さは必ずしも保証しません。間違っていたらごめんなさい。やっつけ仕事感があるので時間に余裕があればさらに込み入った計算を追記しておきます。)
目的
多くの量子力学的な問題はハミルトニアンと呼ばれる系のエネルギーを決める量を設定し、そのハミルトニアンの最低エネルギー(基底状態の基底エネルギーとも呼びます)を求める問題として定義されます。
ハミルトニアンを設定→シュレディンガー方程式→固有値問題を解くといった流れになります。
シュレディンガー方程式は
$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
i\hbar \frac{d}{dt}\ket{\psi(t)} = \mathcal{H}(t)\ket{\psi(t)}
$$
と表されます。$\mathcal{H}(t)$をハミルトニアン、$\ket{\psi(t)}$を波動関数と呼びます。
ハミルトニアン$\mathcal{H}(t)$が時間に依存しない場合、波動関数を変数分離することで解くべき方程式は、
\begin{align}
\mathcal{H}\ket{\psi} &= E\ket{\psi} \\
\ket{\psi(t)} &= \exp[-i\frac{Et}{\hbar}]\ket{\psi}
\end{align}
となります。これは定常状態のシュレディンガー方程式と呼ばれ、線形代数の固有値問題に相当します。
$\mathcal{H}$に対して適当な基底を決めて、行列要素を求めれば後は対角化するだけで固有エネルギーを計算することができます。
しかし、この対角化が以下の理由により難しくなります。
- 行列の次元がシステムサイズに対して指数関数的に増大する
- 量子ゆらぎが存在する
これらの理由について順に説明して行きます。
行列の次元の増大について
電子を対象とする場合は電子の統計性(フェルミ粒子性)により同じ状態を占有できる電子は1種類に限られてしまいます。「同じ状態」とは同じエネルギー(準位)を取る状態であったり、同じ軌道(サイト)を占有する状態などを意味して、扱う問題(そこで取る基底)によって意味は多少変化します。
また電子にはスピンと呼ばれる内部自由度を持っており、アップとダウンの2種類の状態を取ることができます。
このスピンによって電子はさらに区別されるので、同じ状態を占有できる電子はアップとダウンの電子1つずつになります。
例えると、１つの建物に男女別々のシャワールームがあり、そのシャワールームを狭いので男女それぞれにつき1人までしか同時に使用することはできません。
１つの建物あたりで「同時」に起きる状況としては
- 誰もシャワールームを使用していない。
- 男性だけがシャワールームを利用している。
- 女性だけがシャワールームを利用している。
- 男女がシャワールームを利用している。
この4パターンが考えられます。
この例の状態を「軌道」で例えると、ある軌道$i$を占有している電子がとる状態としては
- １つもいない $\ket{0}$
- アップスピンの電子が1つ $\ket{\uparrow}$
- ダウンスピンの電子が1つ $\ket{\downarrow}$
- アップとダウンスピンの電子が1つずつ $\ket{\uparrow\downarrow}$
が考えられます。
ポイントは同じエネルギー軌道にアップが2つとかダウンが2つといった同じスピンの電子が占有できない点です。
この様に1つの軌道に対して4つの状態を取ることが可能なので、全体のヒルベルト空間は電子の軌道が$n$だと$4^n$にまで膨れ上がってしまいます。
このためシステムサイズを増大すると、ハミルトニアンの行列要素が指数関数的に増大してしまいます。
例えば2軌道の場合では$4^2=16$個の状態を考える必要があります。
状態の記述の仕方を$\ket{\psi}=\ket{\psi_1,\psi_2}$とすると、$\ket{\psi}$は
- $\ket{0,0} $
- $\ket{\uparrow,0} $
- $\ket{0,\uparrow} $
- $\ket{\downarrow,0} $
- $\ket{0,\downarrow} $
- $\ket{\uparrow,\uparrow} $
- $\ket{\uparrow,\downarrow} $
- $\ket{\downarrow,\uparrow} $
- $\ket{\downarrow,\downarrow} $
- $\ket{\uparrow \downarrow,0} $
- $\ket{0,\uparrow \downarrow} $
- $\ket{\uparrow \downarrow, \uparrow} $
- $\ket{\uparrow \downarrow, \downarrow} $
- $\ket{\uparrow,\uparrow \downarrow} $
- $\ket{\downarrow,\uparrow \downarrow} $
- $\ket{\uparrow \downarrow,\uparrow \downarrow} $
が考えられます。これらの16個の状態に対して16次元の単位ベクトルを上から順に対応する様にとればハミルトニアンの行列要素を導出することが可能です。
ハミルトニアンに対称性(粒子数保存やスピン保存)があれば上記の様に電子数ごとやスピンの大きさごとにあらかじめ基底の取り方を整えておくと、全体のハミルトニアンの行列がブロックに分かれて対角化するべき行列の次元を下げることも可能です。
量子ゆらぎについて
ただ単に次元が増えても行列要素の中身が解きやすい形であれば次元の大きさは怖くありません。しかし多くの量子力学的な模型には量子揺らぎという同時対角化を妨げる要素が入っているため対角化が難しくなります。
Hubbard Model
$$
\mathcal{H} = -t\sum_{i,j,\sigma}(c^{\dagger}_{i\sigma}c_{j\sigma} + c^{\dagger}_{j\sigma}c_{i\sigma}) + U\sum_i n_{i\uparrow}n_{i\downarrow}
$$
それぞれの演算子の意味などについてはこちらの記事を参考にしてください。
https://qiita.com/ot_yakan/items/0ed316e57e3673e43d0b
第1項目はそのままの軌道規定では対角的ではありませんが、フーリエ変換してエネルギー表示をして表した基底では対角化されています。一方で第2項は現在の軌道基底で既に対角化されています。エネルギー表示にするためにフーリエ変換してしまうと逆にこの対角性が失われています。この第1項と第2項の競合により「あちらを立てればこちらが立たず」的な状況に陥ってしまいます。
Quantum Spin Model
上では電子系を考えていましたが、次はスピン系について考えます。
$$
\mathcal{H} = J\sum_{i,j}\boldsymbol{S}_i \cdot \boldsymbol{S}_j = J\sum_{i,j}\biggl(\sigma^x_i\sigma^x_j + \sigma^y_i\sigma^y_j + \sigma^z_i\sigma^z_j \biggr)
$$
ここで、$\sigma$は
\sigma^x_i = \left( \begin{array}{cc} 0 & 1\\ 1 & 0 \end{array} \right), \;\sigma^y_i = \left( \begin{array}{cc} 0 & -i \\ i & 0 \end{array} \right),\;\sigma^z_i = \left( \begin{array}{cc} 1 & 0 \\ 0 & -1 \end{array} \right)
で定義されるパウリ行列になります。
量子アニーリングではこの量子スピン模型のハミルトニアンを対象とします。
アップスピン、ダウンスピンの2つの状態の基底を$\sigma^z$に対して対角的になるように
\ket{\uparrow}_i = \left( \begin{array}{c} 1\\ 0 \end{array} \right) \\
\ket{\downarrow}_i = \left( \begin{array}{c} 1\\ 0 \end{array} \right)
ととります。これらは$\sigma^z$に対して固有状態となるように定義したので
\begin{align}
\sigma^z_i \ket{\uparrow}_i &= \ket{\uparrow}_i \\
\sigma^z_i \ket{\downarrow}_i &= -\ket{\downarrow}_i
\end{align}
となります。一方で$\sigma^x$については
\begin{align}
\sigma^x_i \ket{\uparrow}_i &= \ket{\downarrow}_i \\
\sigma^x_i \ket{\downarrow}_i &= \ket{\uparrow}_i
\end{align}
となり$\sigma^x$が作用することでスピンのアップとダウンが反転してしまいます。
$\sigma^y$についても前にかかる係数は異なりますが同じ働きをします。
$\sigma^x, \sigma^y$の固有状態は
\begin{align}
\ket{x_{\pm}}_i &= \frac{\ket{\uparrow}_i \pm \ket{\downarrow}_i}{\sqrt{2}}\\
\sigma^x_i \ket{x_{\pm}}_i&= \pm \ket{x_{\pm}}_i
\end{align}
となります。基底を$\ket{x_{\pm}}_i$に変更すれば$\sigma^x,\sigma^y$を対角化できますが、$\sigma^z$が非対角になります。
このように$\sigma^x, \sigma^y$が作用するとスピンが反転してしまい、ハミルトニアンの行列要素は非対角成分を持つようになります。これが量子的な揺らぎと呼ばれて問題を難しくしていると原因になります。
$\sigma^z$だけで構成される模型を古典スピン模型またはイジング模型と呼びます。
組合せ最適化問題は古典スピン模型として定義できます。
$\sigma^z$の固有値は$\pm 1$だったのでこれを$0,1$に変更するには単純に
$$
x_i = \frac{\boldsymbol{1} + \sigma^z_i}{2}
$$
という変換を行うと$x_i$は$\ket{\uparrow}_i,\ket{\downarrow}_i$に対して固有値$0, 1$を取るようになります。
組み合わせ最適化問題の変数$x_i$を組み合わせとして選択する場合は1、選択しない場合は0を取ると読み替えればおkです。
この時のハミルトニアンの固有状態は
$$
\ket{\psi} = \prod^N_i \otimes \ket{x}_i \equiv \ket{x_1 x_2 \cdots x_N}
$$
と表すことができます。例えばかなり単純な例題ですが、
$$
\mathcal{H} = x_1 - 2x_2
$$
という問題を考えます。考えられる固有状態は
\begin{align}
\ket{\psi} = \begin{cases}
\ket{00} & \rightarrow E = 0 \\
\ket{10} & \rightarrow E = 1 \\
\ket{01} & \rightarrow E = -2 \\
\ket{11} & \rightarrow E = -1
\end{cases}
\end{align}
であり、基底状態は$\ket{01}$であり基底状態のエネルギーは$-2$であることが分かります。
今は2変数で$2^2=4$通りの固有状態のみを考えれば良いので全ての固有状態を列挙し、最低エネルギーをとる状態を調べることが可能ですが、システムサイズ$N$が増大すると固有状態の数は$2^N$のオーダーで増大するので全列挙が困難になります。
量子アニーリング
量子アニーリングの目的は古典的なイジング模型の基底エネルギーを求めることです。
$$
\mathcal{H}_0({ \sigma^z }) = \sum_i J_i \sigma^z_i + \sum_{i,j} J_{ij} \sigma^z_i\sigma^z_{j} ++ \sum_{i,j,k} J_{ijk} \sigma^z_i\sigma^z_{j}\sigma^z_{k} + \cdots
$$
このハミルトニアンは$\sigma^z$のみから構成されます。この基底状態を求めるのはNP困難な問題です。これに対して次の量子ゆらぎを$\mathcal{H}_0$に加えます。
$$
\mathcal{H}_1 = -\sum_i \sigma^x_i
$$
$\mathcal{H}_1$のみであれば基底状態は簡単に
$$
\ket{\psi}_x = \prod^N_i \otimes \frac{\ket{\uparrow}_i + \ket{\downarrow}_i}{\sqrt{2}}
$$
と簡単に求めることができます。そこで全体のハミルトニアンを
$$
\mathcal{H} = \frac{t}{T}\mathcal{H}_0 + (1-\frac{t}{T})\mathcal{H}_1
$$
と定義して、$t=0$では完全に基底状態が分かっている$\mathcal{H} = \mathcal{H}_1$に、$t=T$で実際に知りたいハミルトニアン$\mathcal{H} =\mathcal{H}_0$になる様に設定します。
これは時間に依存するハミルトニアンになったので、元のシュレディンガー方程式に戻って波動関数の時間発展を追いかける必要があります。
一般的に初期時刻の固有状態と最終時刻の固有状態は別物になりますが、以下の断熱定理を利用すると$T$を十分大きくとれば、$\mathcal{H}_1$の固有状態$\ket{\psi}_x$から出発して、$\ket{\psi}_x$が最終時刻$T$で$\mathcal{H}_0$の基底状態になる様に自動的に変化してくれるのです。
$t=0$では$\mathcal{H}(0)=\mathcal{H}_1$で$\ket{\psi}_x$が基底状態でした、$t=T$ではもともとエネルギーが知りたかった$\mathcal{H}(T)=\mathcal{H}_0$になり、その基底状態は$\ket{\psi}_x \rightarrow ?$になります。最終的な状態はハミルトニアンの具体的に形によってどうなるかは分かりませんが、出てきた状態は必ず基底状態になっていることが保証されているということです。
先ほどの例を量子アニーリングで計算すること考えると、まず次のハミルトニアン
$$
\mathcal{H}(t) = \frac{t}{T}(x_1 - 2x_2) - (1-\frac{t}{T})(\sigma^x_1 + \sigma^x_2)
$$
と第2項目の基底状態
$$
\ket{\psi}_x = \frac{\ket{\uparrow}_1 + \ket{\downarrow}_1}{\sqrt{2}} \otimes \frac{\ket{\uparrow}_2 + \ket{\downarrow}_2}{\sqrt{2}} = \frac{1}{2} \biggl(\ket{\uparrow_1 \uparrow_2} + \ket{\uparrow_1 \downarrow_2} + \ket{\downarrow_1 \uparrow_2} + \ket{\downarrow_1 \downarrow_2}\biggr)
$$
を用意して$t=0$から時間を$t=T$まで発展させると、最初の4つの線形和の中に含まれる状態のどれか1つに落ち着くことになります。
断熱定理
ハミルトニアンは時間に依存していますが、$T$が十分大きく系が断熱的にゆっくりと変化していると仮定し、時刻$t$でパラメトライズされた定常状態のシュレディンガー方程式
$$
\mathcal{H}(t)\ket{j(t)} = \epsilon_{j}\ket{j(t)}
$$
を考えます。ここで$j$はエネルギー固有値の順番を決めるインデックスであり
$$
\epsilon_0 \le \epsilon_1 \le \epsilon_2 \le \cdots
$$
の様に並んでいると仮定します。なので$j=0$が基底状態に対応します。
そうすると各時刻における系の状態$\ket{\psi}$は$\ket{0(t)}$に十分近くなる、というのが断熱定理になります。もう少し数学的に記述すると
|\braket{\psi(t)}{0(t)}|^2 = 1-\delta^2\;\;(\delta\ll 1)
になります。
$\braket{\psi(t)}{0(t)}$は$\ket{\psi(t)}$と$\ket{0(t)}$との内積を意味します。
$\delta$は次の関係から導かれる微小量になります。
\begin{align}
\frac{\max |\braket{1(t)}\partial_t \mathcal{H}(t){0(t)}|}{\min \Delta\epsilon^2_t} = \delta
\end{align}
ここで
$$
\Delta\epsilon = \epsilon_1(t) - \epsilon_0(t)
$$
であり基底状態の1つエネルギー準位の高い「励起状態」と基底状態のエネルギー差を意味します。$|\braket{1(t)}\partial_t \mathcal{H}(t){0(t)}|$を具体的に計算すると次の関係を導くことができます。
$$
T \propto \frac{1}{\delta \min \Delta\epsilon^2_t}
$$
$\Delta\epsilon^2_t$が小さい場合、つまり励起状態と基底状態のエネルギーギャップが縮まる様な系では計算時間が無限大になってしまいます。
より厳密な計算は量子アニーリングの提唱者である西森先生のレビュー論文を参考にしてください。
https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/189516/1/bussei_el_033203.pdf
(時間があればさらに詳しく上のレビュー論文を解説したいと思います。)
1次相転移付近ではこのエネルギーギャップが消失するケースが多く観測されています。
故に1次相転移をいかに回避してアニーリングするかが重要になってきています。
最近では多体の量子ゆらぎを導入して1次相転移を回避する手法を西森先生が提案されている様です。
これはまだD-Waveのマシンでは実現されていませんが。
QUBO
D-Waveで解くことのできるイジングモデルは$\sigma^z$に対して2次までの項で構成されるハミルトニアンになります。さらにこのハミルトニアンは変数間の制約条件を全て加味した形の定式化が必要となり、QUBO(Quadratic Unconstrained Binary Optimization)と呼ばれる形で定式化する必要があります。(エネルギーの値を決めるのにif文などの条件分岐が使えない。)
QUBOは次の様なハミルトニアンで定義されます。
\mathcal{H} = \sum^N_{i=1}c_ix_i + \sum^N_{i=1}\sum^N_{j=1}Q_{ij}x_i x_j\;\;\;(x_i \in \{ 0,1 \},\{c_i,Q_{ij} \} \in \boldsymbol{R})
まぁ$x^2_i=1$なので$c_i$は$Q_{ij}$の対角成分に含むことができるので
$$
\mathcal{H} = \sum^N_{i=1}\sum^N_{j=1}Q_{ij}x_i x_j = \boldsymbol{x}^T \boldsymbol{Q} \boldsymbol{x} = \sum^N_{i=1}\sum^N_{j=1}Q_{ij}x_{ij}
$$
と表すことが出来ます。
巡回セールスマン問題では$Q_{ij}$に都市$i$と$j$間の距離が格納されており、$x_{ij}$は訪問順序$i$番目に都市$j$を通る場合は1と通らない場合は0を取る変数になります。
巡回セールスマン問題では各訪問順序番号では1都市だけを周り、各都市は1回だけしか訪問されないという制限が付与されます。この制約条件をハミルトニアンに加える必要があります。
具体的には上記の制約を破ったときにエネルギーが増加する様に
$$
\sum^N_{i=1}(x_{ij} -1)^2 = 1, \sum^N_{j=1}(x_{ij} -1)^2 = 1
$$
この式に適当な大きさの係数をかけて全体のハミルトニアンに加えます。
ラグランジュの未定乗数法や、不等式の場合はカルーシュ・クーン・タッカー条件(KKT条件)を利用して制約をハミルトニアンに付与して、非制限な2次式を設定して行くという手順になります。
最後に
D-waveでは初期状態である重ね合わせ状態の実現や量子ゆらぎの導入に祭して超電導回路を利用しており、そのハードウェア的制約により1変数あたりに結合できる相互作用数が6~8程度に絞られるという制限がついてしまいます。(キメラグラフなどとも呼ばれています。)
一応、制限の結合数を超えた相互作用を利用したい場合は余分なビット介して実現することも可能です。
また3次以上からなる変数項も補助ビットを導入すれば実現することは可能ですが、それに伴い必要なビット数も増大してしまいます。
量子現象は利用していませんが、QUBOを解くためのアーキテクチャとして弊社もデジタルアニーラを提供しています。
デジタルアニーラの利点としては相互作用の結合数に制限がなく全結合に対応している点です。
欠点はやはりまだまだ対応可能なビット数が少ない点です(第1世代では1024ビット、第2世代では8192ビット)。
ビット数と相互作用の結合数はトレードオフの関係にあるそうですね。
現実的な組み合わせ最適化問題をわざわざイジング模型にマッピングして解くことの利点辞典もまだよくわからない状態ですね。
物理的には3次元イジング模型などがより効率的に解ける様になれば美しいと思いますが。
量子的なゆらぎを考慮したハミルトニアンをQUBOとして定式化しD-Waveのマシンで解く様な研究も量子化学の分野で行われています。
その話題についても後日解説したいと思います。
あとはQUBOの定式化についてなども時間あれば今後解説していきたいと思います。