Table of contents
  1. Quantum annealing for you 05
    1. ブラックボックス最適化
    2. スクリプト

Quantum annealing for you 05

ブラックボックス最適化

\(\mathbf{x} = (x_1, x_2, \dots, x_N), x_n = \{ 0, 1\}\)のバイナリ変数列があるとき、コスト関数\(E(\mathbf{x})\)を最小化(最大化)するようなバイナリ変数列を求めることが量子アニーリングでできるのでした。ただし、このコスト関数は

\[E(\mathbf{x}) = \mathbf{x}^T Q \mathbf{x} = \sum_{i=1}^N \sum_{j=1}^N Q_{ij} x_i x_j\]

のような形にする必要があり、現実の問題をこの定式化に落とし込む必要があります。
ではそもそもこの形で書けるかどうかわからない場合はどのように量子アニーリングで解けば良いでしょうか。それをなんとかする方法がブラックボックス最適化です。

今、\(E(\mathbf{x})\)は未知(知ることができない)としましょう。すなわち、バイナリ変数列\(\mathbf{x}\)をこの関数に入れたときに、どのような値が出力されるかその関数系はわからないものとします。

\[y = E(\mathbf{x})\]

我々にはQUBO行列を作ることしかできないので、この関数を真似するような\(E(\mathbf{x}) \simeq \sum_i \sum_j Q_{ij} x_i x_j\)を求めましょう。このとき何もヒントがなくては真似しようがないため「\(\mathbf{x}\)を入れると\(y\)は教えてくれる」、すなわちヒントとなるデータセット\((\mathbf{x}_d, y_d) \ d = 1\sim D\)が与えられ、それを手掛かりに\(Q_ij\)を推定します。ここでは

\[\min_{Q} \left\{ \sum_d \left( y_d - \sum_{i} \sum_j Q_{ij} x_i^d x_j^d \right)^2 + \lambda \sum_{i, j} Q_{ij}^2 \right\}\]

のような最小2乗法を考えます。最後の項、データ(手掛かり)が少ないのに対してパラメータが多い場合を考えるため、過学習を防ぐためのものL2ノルム(の2乗)正則化項です。このようにすると、これは\(Q_{ij}\)の2次関数であるため、この関数の微分が0となるところを探せばそこが最小値になることから、簡単な形であるからです。
例えば\(N=10\)だとすると、\(Q_{ij}\)の数は100個です。これに対して\(D=5\)だとしましょう。するとこの5個のデータに対して100個のパラメータ数は大き過ぎます。5個のデータを100個のパラメータを調整すると、それらのデータに合わせ過ぎてしまうことがあります。すると6個目のヒントデータが出てきた時に、その合わせ過ぎてしまった\(Q_{ij}\)では対処できないことがあります。よって不測の事態へのズレに対処するために、よそ見をさせる効果を持つものを正則化項と呼びます。
L1ノルムをつけても良いですが、その場合には近接勾配法を用いる必要があります。

\(\sum_i \sum_j Q_{ij} x_i x_j\)$をそのまま使うのは効率が悪いので書き換えを行いましょう。

\[\mathbf{x} = (x_1, x_2, \dots, x_N) \ \Longrightarrow \  \mathbf{x}' = (1, x_1, x_2, \dots, x_N, x_1 x_2, x_1 x_3, \dots, x_{N-1} x_N)\]

のように1とクロスタームも含めた形に書き直しましょう。また

\[\mathbf{a} = (Q_0, Q_{11}, Q_{22}, \dots, Q_{NN}, Q_{12}, Q_{13}, \dots Q_{N-1, N})\]

のように対角項と非対角項を並べたベクトルを作ります。ここで\(Q_0\)は\(E(\mathbf{x}) \simeq \sum Q_{ij} x_i x_j + Q_0\)のように、\(x_i\)の項で表現できない定数項です。
例えば\(Q_{11}\)は\(Q_{11} x_1 x_1 = Q_{11} x_1\)のように書けるため、上の\(\mathbf{x}\)の成分の並びと対応していることがわかります。このようなものを定義しておくと

\[\sum_i \sum_j Q_{ij} x_i x_j = \sum_k a_k x_k'\]

のようにベクトルの内積で表現されることがわかります。\(\mathbf{x}'\)の中に\(x_2 x_1\)がないのは、QUBO行列が対称行列であるため、考える必要がないからです。
このような手続きによって\(Q_{ij}\)を求める問題が\(a_k\)を求める問題に変換されました。この形にしておけば線形回帰(\(x\)の1次)として問題を解くことができます。
よって

\[\min_\mathbf{a} \left\{ \left( y_d - \sum_k a_k x_k'^d \right)^2 + \lambda \sum_k a_k^2 \right\} \tag{1}\]

のようになり、微分が0の場所もしくは平方完成することで最小値を求めることができます。

\(a_k\)で微分してみましょう。すると

\[-2 \sum_d (y_d - \sum_\ell a_\ell x_\ell'^d) x_k'^d + 2 \lambda a_k = 0 \ \Longrightarrow \  \sum_d \sum_\ell x_\ell'^d x_k'^d a_\ell + \lambda a_k = \sum_d y_d x_k'^d\]

ここで\(x_k'^d\)を行列だと考えて、\(X = (x_k'^d)\)のように書くと

\[X X^T \mathbf{a} + \lambda \mathbf{a} = X \mathbf{y} \ \Longrightarrow \ (X X^T + \lambda I_N) \mathbf{a} = X \mathbf{y}\]

逆行列を用いることで

\[\mathbf{a} = (X X^T + \lambda I_N)^{-1} X\mathbf{y} \tag{2}\]

のように、既知の量から\(\mathbf{a}\)を求めることができます。
\(\lambda\)はL2ノルム正則化のために導入したパラメータです。\(\lambda\)を大きくすると、(1)式から\((\cdots)^2\)部分をあまり信用しない(データを信用しない)結果を導くことができます。また\(\lambda\)を小さくするとデータを信用するような解き方になります。
まずデータセットから\(\mathbf{a}\)(QUBO行列)を計算します。それを元にD-Wave計算をすると最適解の候補が求められます。その最適解の候補をブラックボックスとなっている\(E(\mathbf{x})\)に入力すると、コスト関数の値が求まり、それが本当に最適解かどうかがわかります。これを繰り返すことで\(\mathbf{a}\)を洗練し、\(E(\mathbf{x})\)の関数を近似したものを求めます。
ブラックボックス最適化は\(x, y\)の両方が必要で、それらが良いか悪いかを判断する必要があります。ボルツマンマシンは\(x\)のみを用いて、\(x\)の傾向のみが必要です。
\(E(\mathbf{x})\)がQUBO行列で書かれるかどうかもわからないので、上述の議論で推定された結果が正しいかどうかは定かではありません。しかし限られたヒントから最適解の候補を見つけることができるというのが、ブラックボックス最適化の素晴らしい点です。ただしQUBOよりも複雑な形をしている場合には、近似がうまくいきません。

ちなみに平方完成では

\[\sum_d \left( y_d - \sum_k a_k x_k^d\right)^2 + \lambda \sum_k a_k^2 = \sum_d y_d^2 - 2 \sum_{k, d} a_k x_k^d y_d + \sum_{d, k, \ell} a_k x_k^d x_\ell^d a_\ell + \lambda \sum_k a_k^2 = \sum_d y_d^2 - 2 \sum_k a_k (\sum_d x_k^d y_d) + \sum_{k, \ell} a_k (XX^T + \lambda I_N)_{k\ell} a_\ell \tag{3}\]

のようになります。2次関数の頂点部分は今までのデータが一番良く説明できる箇所に当たります。ではその推定がどれくらい正しいと言えるでしょうか。それはこの下に凸の2次関数の鋭さからわかります。

  • 鋭いならば、その推定から少しずれただけでデータを説明できなくなることを意味する。よってこの場合はその推定は自信がある。
  • 鈍いと、その推定から少しずれてもデータをまあまあ説明することができる。よってその推定には自信がない。

(3)式から、その鋭さが\(XX^T + \lambda I_N\)であることがわかります。

スクリプト

P = 1 + N + N * (N-1) // 2

のように//2を用いると、Pがfloatではなくintになります。キャストする必要がないので便利です。

この計算のボトルネックは逆行列の計算です。またデータが多くなると、その通信量が大きくなります。

2次関数の鋭さと解の発見パターンには、関係があります。

var = 0.5 * XXinv

において係数を0にしてしまうと、極小解にハマりやすくなります。逆にこの係数を大きくすると、手がかりを信用せずにあてもなく探索をするため、良い解が求まりません。

最初に用意するデータはできるだけ似ていないものを選んだ方が、QUBOを推定しやすくなります。これは大局的な振る舞いを知るためです。ただし集められてきたデータがあまりにもバラバラだと、そこから得られるヒントもバラバラなので、上手く探索を深堀りできないことがあります。よってデータが集まってきたら、ある程度そのQUBO推定結果を信じる方が良いです(アニーリングのようにvarを小さくした方が良い?)。


Copyright © github-nakasho