Table of contents

1. Quantum annealing for you 04
   1. ボルツマン機械学習
      1. データ整形(one-hot表示, one of k 表示)
      2. ボルツマン機械学習
      3. 制約
      4. 等式・不等式制約について
      5. D-Wave, OpenJijで解く
      6. QUBOの調整
   2. 量子アニーリングについて
   3. フォルクスワーゲンの論文再現
      1. 環境設定
      2. 論文のすごいところ
      3. 定式化
      4. QUBO行列の作成
      5. 表示方法

Quantum annealing for you 04

ボルツマン機械学習

データ整形(one-hot表示, one of k 表示)

\((3, 2, 2, 2)\)の選択肢を選んだ人のデータを\(0001, 001, 00100, 0010000\)のように表現します。発想としては2進数表現と同じです。
one-hot表示はすごく無駄が多いです。2進数表現の方がビットは節約できますが、最適化する部分でコスト関数が複雑になります。ただし、本当にその部分が問題になるようならば、その方法を疑う方が良いでしょう。

ボルツマン機械学習

246人の[0, 0, 1, … , 0, 1]が、実は量子アニーリングの出力結果ではないか、と見なすことができます。よってこの結果をD-Waveマシンに真似させることを考えます。そうなるようにQUBO行列を考えましょう。とあるデータセット結果を真似して、同じようなデータを生成できるものを生成モデルと呼びます。
この技術は似たような画像や音楽を生成させるようなプログラムを作成するのに用いることができます。
ボルツマンマシンは1970~80年代からと歴史は古くありましたが、2000年代に入ってGAN(敵対的生成)が盛り上がりを見せました。GANは生成されたデータを、全てを知る監督者がダメ出しをします。そのダメ出しを受けてよりよいデータを生成するという手法です。
ボルツマンマシンの構造をより複雑にしたディープボルツマンマシンなどの人工知能への応用も期待されています。
ボルツマンマシンは試しに今のQUBO行列で[0, 1]の配列を出すのに非常に時間がかかるという欠点があります。そのような欠点を補うためにD-Wave量子アニーリングマシンを用います。ただし、量子アニーリングはギブス・ボルツマン分布に従う様々な解を出します。D-Waveマシンの出力結果がギブス・ボルツマン分布に従うならば、それに従う結果を出力するようなボルツマンマシンに応用できるのではないか、というわけです。
OpenJijはギブス・ボルツマン分布に従うことを利用したシミュレータです。

制約

先ほどのone-hot表示から、例えば4つの選択肢のうちの1つは選ばなければならないという条件より\(x_1 + x_2 + x_3 + x_4 = 1\)と考えられます。次の3つの選択肢のうちの1つを選ぶ必要がある条件より\(x_5 + x_6 + x_7 = 1\)のように、それぞれの区間で足したら1になるというのが制約です。よってコスト関数は

\[E(\mathbf{x}) = \lambda (x_1 + x_2 + x_3 + x_4 -1)^2 + \lambda (x_5 + x_6 + x_7 -1)-2 + \lambda (x_8 + x_9 + x_{10} + x_{11} + x_{12} - 1)^2 + \lambda (x_{13} + x_{14} + x_{15} + x_{16} + x_{17} + x_{18} + x_{19} - 1)^2\]

しかしこれだけではランダムになってしまいます。よってここから「どの選択肢を選びそうか」の傾向を元のデータから真似る必要があります。それを機械学習で叶える必要があります。それには

\[Q_{ij} x_i x_j\]

を加えます。今の時点では\(Q_{ij}\)は未知なもので、良い答えが出せたときはそれを採用し、ダメな答えならそれをさらに調節・改良していきます。そうして求めた\(Q_{ij}\)を詳細に見ることで\(x_i, x_j\)の相関を見ることができます。

それでは罰金項を具体的に計算しましょう。

\[(x_1 + x_2 + x_3 + x_4 -1)^2 = x_1^2 + x_2^2 + x_3^2 + x_4^2 + 1 + 2 x_1 x_2 + 2 x_1 x_3 + 2 x_1 x_4 + 2 x_2 x_3 + 2 x_2 x_4 + 2 x_3 x_4 - 2 x_1 -2x_2 - 2 x_3 -2 x_4\]

\(\{ 0, 1\}\)のバイナリ変数では\(x = x^2\)です。またまた\(2 x_1 x_2 = x_1 x_2 + x_2 x_1\)です。これでQUBO行列を簡単に作ることができます。

もし1つではなく2つ選択するような複数回答制約を考える場合には、\(2 x_1 x_2 \rightarrow 4 x_1 x_2\)のようになるだけです。このようなクロスタームの係数を\(2k\)のようにしておくことで、\(k\)の値を変えるだけで応用することができます。

等式・不等式制約について

選択しなけらばならない数が1, 2の場合は1.5とすれば良いでしょう。何個でも自由なら罰金法は入りません。最適化問題で不等式制約はあまり聞かない場合があります。不等式で入れてもぴったりイコールが解であることが多いので、等式制約でも良いという場合があります。

D-Wave, OpenJijで解く

DWaveCliqueSamplerという全結合グラフ埋め込みで考えます。これはQという行列がどのように結合しているかわからない(全部繋がら雨可能性もある)ため、全結合をバランスよく埋め込む方法を用います。OpenJijのSQASamplerを用いる場合には埋め込みは必要ありません。
ボルツマンマシンの場合にはanswer_mode='raw'を選択すると良いでしょう。
もしOpenJijを用いる場合には

Qdict = {}
for i in range(M):
    for j in range(M):
        if Qmat[i][j] != 0.0:
            Qdict[i, j] = Qmat[i][j]

のようにして辞書を用意します。またOpenJijにはanswer_modeがないため、出てきた結果から「その解が出た頻度」を考慮する必要があります。

for k in range(len(sampleset.record))
    x[k, :] = sampleset.record[k][0]

QUBOの調整

それではD-Waveマシンの出す結果を用いて、QUBO行列を調整していきましょう。現在は

\[Q_{ij} x_i x_j + (\mathrm{penalty})\]

のように実装されています。目的は\(Q_{ij}\)を調整することです。仮に\(Q_{ij}\)を大きくすると、最小化するために\(x_i x_j\)を小さくしようとします(0になろうとします)。次に\(Q_{ij}\)を小さくすると\(x_i x_j\)が大きくなります(1になろうとします)。
元のデータを\(Z_{im}\)(N行M列), D-Waveマシンから出た結果を\(X_{im}\)(Nsample行M列)としましょう。これらを1つずつ比べるのではなくて平均的に比べていきます。\(m\)の添字はある項目についての結果です。元々\(Q_{ij}x_i x_j\)は\(i\)番目の項目と\(j\)番目の項目の関係を表すものだったので

\[\sum_i \frac{Z_{im} Z_{in}}{N}, \quad \sum_i \frac{X_{im} X_{in}}{\mathrm{Nsample}}\]

の2つを比べることにします。統計学では\(\sum_i \frac{X_{im} X_{in}}{\mathrm{Nsample}}\)を共分散行列と呼びます。この共分散行列の値が小さければ\(Q_{ij}\)を小さくすることで\(X_{im} X_{in}\)を大きくし、値が大きければ\(Q_{ij}\)を大きくすることで\(X_{im} X_{in}\)を小さくします。このようにして\(\sum_i \frac{Z_{im} Z_{in}}{N}\)と比べた結果から、\(Q_{ij}\)を自動調整するように機械学習させます。

もしOpenJijを使っている場合には

Xmat = np.zeros((M, M))
for k in range(Nsample):
    Xmat += sampleset.record[k][2] * np.dot(x[k, :].reshape(M, 1), x[k, :].reshape(1, M)) / Nsample

のようになります。よって次のステップのQUBO行列を

\[Q_{mn} - \eta \left( \sum_i \frac{Z_{im} Z_{in}}{N} - \frac{\sum_i X_{im} X_{in}}{\mathrm{Nsample}}\right)\]

とします。これを更新式と呼びます。\(\eta\)は程々の値(機械学習では0.1, これを学習率や勾配率と呼びます)にします。上式の()部分はボルツマンマシン学習の損失関数の勾配(微分)なので、ただ勾配法を用いているだけです。ボルツマンマシンラーニングの損失関数にはKL情報量を用いますが、ここでは説明を割愛します。
Chain-breakの結果はXmatにノイズが入ります。機械学習にはノイズを入れて学習する手法があるため、それに似たようなものと考えることもできます。ただし、制御のできないノイズであるため注意が必要ですが、このおかげで勝手に正則化のような効果がされることが知られています。
SMBCが量子アニーリングでオーバーサンプリングしたような結果があるが…機械学習の分野だとデータの水増し・生成などを行っているが多分同じこと。もしD-Waveで同じようなデータを作れるなら、取引データが100個しかなくても、もっとデータを作ることが可能という考え方です。大関さんは「結局、元の100個のデータ以上の意味はないはずだから、オーバーサンプリング自体にあまり強い意味はないのではないか?その外挿大丈夫?裏に物理モデルがあれば意味があるかも…?」との意見。データの構造が全くわかっていないときにボルツマンマシンを使うことは危険かもしれない?
結果更新が機械学習のSGDのようになっているため、Adamのような別の最適化式を用いることも可能です。

量子アニーリングについて

横磁場をだんだん弱めることで\(\uparrow\)と\(\downarrow\)の重ね合わさった状態から、\(\uparrow, \downarrow\)の状態どちらかに決まる。
しかしD-Waveはノイズの影響及び横磁場を切っているのに横磁場が残った状態の結果が出るfreezing現象により、理想的な量子アニーリングではない結果が出ることがあります。しかしこのおかげでギブス・ボルツマン分布に従うという特性があるため、ボルツマンマシンとして応用することができます。

フォルクスワーゲンの論文再現

環境設定

OSMnxの描画にはmatplotlibを用いているため、そのバージョンによっては描画が実行されないなどの不具合があります。

論文のすごいところ

タクシーを乱数で配置し、その出発点から到着点までの経路を最適化します。同じタクシーが同じ経路を通らないように避けることで、渋滞が起こらないようにします。

• 新しい問題設定を量子アニーリングで解けるように定式化したこと
• QUBOに実際のデータを用いたこと

の2点が特に新しい点です。また手法として

1. (古典コンピュータで)地図データを前処理
2. (古典コンピュータで)渋滞しそうなところを抜き取る
3. (古典コンピュータで)もしできるのであれば、それぞれの車に対して代わりになる道を探す
4. (古典コンピュータで)QUBOを作る
5. (量子・古典のハイブリッドで)全体の交通グラフから割り当てられた経路が混雑を回避するような最適化問題を解く
6. (古典コンピュータで)得られた結果を各車に再分配する
7. 2-6を繰り返す

量子コンピュータが得意な部分だけを量子コンピュータに任せるという戦略です。この論文の素晴らしい点はそれらも含めた手法全てを示したところにもあります。

量子アニーリングで解く問題として有名なのは巡回セールスマン問題です。ですが、実はこれには制約条件として罰金法が2つ加わるために、量子アニーリングで探索する良さが活かせず、性能が出ません。このことから「量子アニーリングは巡回セールスマン問題も解けない」というイメージが先行し、物流・交通などでその期待感が薄れました。

定式化

\(x_{ik}\)…\(i\): 車のindex, \(k\): ルートのindex

ここでルートというのは、「現在地から目的地までの長い道のり」という意味です。1量子ビットで広い範囲の道路を表現します。

まるで経路積分のように、\(k=1, 2, 3\)のルートのコスト関数を計算し、そのコストが最小になるようなコスト(作用)を選びます。この論文では1つのルートから1つを選ぶため

\[\sum_{k=1}^3 x_{ik} = 1 \ \Longrightarrow \  \sum_{i=1}^N \lambda \left( \sum_{k=1}^3 x_{ik} - 1\right)^2\]

のように罰金項を作ります。\(1 \leq k \leq 3\)なのは当時のD-Waveマシンのスペックのためです。
次にコスト関数を設定しましょう。

道路のセグメント\(e\)に乗っかる車の台数をコストとし、それを最小にするという考え方をします。\(C_{ek}\)は道路セグメント\(e\)をルート\(k\)が埋めるかどうかを表し、埋めるとき1, そうでないとき0となるような変数です。これを用いると、コストは

\[\sum_e \left( \sum_i \sum_k C_{ek} x_{ik} \right)^2 = \sum_e \sum_i \sum_j \sum_k \sum_{\ell \neq k} C_{e k} C_{e \ell} x_{ik} x_{j\ell}\]

のように書けます。ここで()部分は道路\(e\)を何台の車が通るか、ということを意味するものです。それを2乗するために、\(i\)番目の車と\(j\)番目の車のクロスタームが発生し、その2台の車が同じ道路\(e\)で占有することは避けるという効果(斥力)が生まれます。
もし2乗しなければ、「ある道路を埋める車の台数は減らす」「でもどこかは通らなければならない」ということから、道路をできるだけ埋めないような、最短経路を選択するようになります。\sum_e C_{ek}は\(k\)ルートの距離\(d_k\)を表すので、それを最小化してしまうというように考えることもできます。また2乗にすると、車の台数が増えるごとに2乗で大きくなることから、最小以外を選んだときのコストの増大を表現すること効果も入ります。
\(k \neq \ell\)の部分は、もし\(k = \ell\)ならば\(C_{ek} C_{e\ell} = C_{ek}\)となります。そして\(k = \ell\)が成り立つのは同じ車のときだけ(もし違う車ならば、そもそも出発地と目的地が異なるため、同じルートにはならない)なので、同時に\(i = j\)も成立します。よって\(k = \ell\)を入れると、\(c_{ek} x_{ik}\)の効果(最短距離を選ぶような効果)を入れていることになります。
もし時刻\(t\)まで考えるのであれば

\[\sum_t \sum_e \left( \sum_i \sum_k C_{etk} x_{ik}\right)^2\]

のように拡張すれば、\(C_{etk}\)は\(e\)を時刻\(t\)にルート\(k\)が占有するかどうかを考える\(\{ 0, 1\}\)変数とします。このようにすることで、遅く\(e\)に到達したものと早く\(e\)に到達したもので同じ渋滞に巻き込まれないような最適化問題にすることができます。
もし\(e\)を時間\(t\)にして、\(C_{tk}\)を時刻\(t\)にお店\(k\)を訪れるかどうかの変数とすれば、お店の混雑解消問題に応用することができます。
\(()^2\)を\(()^4\)とすれば4台の相互作用を考えることができます。そうすれば4台のうち3台通っても1台通らなければコストに寄与しないということをすることができます。
空間分割\(e\)と時間分割\(t\)の総数は量子ビットの添字とは無関係なので、古典コンピュータのスペックによります。

QUBO行列の作成

今、\(x_{ik}\)の添字を\(m = k + i * \mathrm{num_route}\)のようにしましょう。するとコストは

\[\sum_e \left( \sum_m C_{em} x_m\right)^2 = \sum_e \sum_m \sum_n C_{em} C_{en} x_m x_n\]

また制約から

\[\lambda \sum_i \left( \sum_k x_{ik} - 1\right) \left( \sum_\ell x_{i\ell} -1\right) = \lambda \sum_i \sum_k \sum_\ell x_{ik} x_{i\ell} - 2 \lambda \sum_i \sum_k x_{ik} + \mathrm{Const}\]

もし巨大な問題などを解き、メモリ消費量を気にするのであればQUBO行列をスパース化する方が良いでしょう(行列\(C\)は0の多いスパース行列なのでそこから工夫するのが良いかもしれません)。もしQUBO行列をそのまま使おうとすると、ハイブリッドでもとても時間がかかります。

表示方法

mnosxは内部でmatplotlibを用いています。内部をいじれば、経路ごとに色を変更することも可能です。