量子アニーリングとスペクトラルクラスタリングの繋がり

# 量子アニーリングとスペクトラルクラスタリング
## 量子アニーリングとは
私は専門家ではないので詳しいことは[他所](https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/189516/1/bussei_el_033203.pdf)に丸投げします。他にも、[ここ](https://qiita.com/ab_t/items/8d52096ad0f578aa2224)とか[そこ](http://blog.brainpad.co.jp/entry/2017/04/20/160000)とか[あそこ](http://www-adsys.sys.i.kyoto-u.ac.jp/mohzeki/QA.pdf)とか良質な資料がたくさんあるのでそちらをご覧ください。

私が理解できる範囲でざっくりまとめると、量子アニーリングマシンとは次のようなコスト関数$H$を最小化する基底状態${\bf \sigma}$を探索する装置です。

```math
H = -\sum_{(i,j)}J_{ij}\sigma_i \sigma_j \tag{1}
```

ここで$\sigma_i, \sigma_j$はそれぞれ$i$番目の格子点と$j$番目の格子点におけるスピンを表す変数であり、状態${\bf \sigma}$は、その$i$番目の成分について$\forall i, \sigma_i \in \\{1, -1\\}$という制約を持ちます。また、$J_{ij}$は$i$番目の格子点と$j$番目の格子点の間における相互作用係数であり、任意の実数値を設定することができます。

量子アニーリングを「アルゴリズム」ではなく「量子アニーリングマシンとは〜探索する装置です」と表現したのは、量子アニーリングマシンはこのコスト関数を最小化するために数値計算的な計算手順を必要としないからです。実は$H$はハミルトニアンといって、系全体のエネルギーに対応しています。そこで外部から印加する磁場や量子場を調節して系の温度を徐々に下げてやると、水が低きに流れるがごとく、量子力学の法則にしたがって系全体がエネルギー最小の状態に落ち着きます。私たちは相互作用係数$J_{ij}$を解くべき問題に応じた値に固定し、温度が十分に冷えるまで待って、あとで${\bf \sigma}$を観測すれば$(1)$式を最小化する解を得られるという仕組みです。

最適化と言っていますが、量子力学的な性質から観測結果は確率的なものになるので、常に最適解が得られるとは限りません（実際に試すと準最適解が多く出てきます）。

## スペクトラルクラスタリングとは
スペクトラルクラスタリングはグラフカット最小化問題を利用して主に2値分類問題に適用されるクラスタリング手法です（[参考1](https://en.wikipedia.org/wiki/Spectral_clustering),[参考2](https://www.slideshare.net/pecorarista/ss-51761860),[参考3](http://d.hatena.ne.jp/mr_r_i_c_e/20121214/1355499195),[参考4](https://www.cs.cmu.edu/~aarti/Class/10701/readings/Luxburg06_TR.pdf),[参考5](https://www.slideshare.net/soyeon1771/spectral-clustering)）。

[参考5: SPECTRAL CLUSTERING](https://www.slideshare.net/soyeon1771/spectral-clustering)にめちゃくちゃ分かりやすいスライドがあったのでそのままお借りいたします。
![spectral-clustering-7-638.jpg](https://qiita-image-store.s3.amazonaws.com/0/47653/5aeb8703-8aeb-0027-5724-1aaae70deda0.jpeg)

まずは2種類に分類したいデータそれぞれに番号を振って$s_i$とし、$s_i$を一列に並べたベクトルを${\bf s}$と表記することにします。また、$s_i$と$s_j$の類似度$w(s_i, s_j)$を$ij$成分に持つ類似度行列$W$を定義します。

```math
W_{ij} = w(s_i, s_j)
```

類似度$w(s_i, s_j)$には、類似度を表す適当な関数（$s_i$と$s_j$が近いものであるほど大きい値を取り、遠いものであるほど小さい値を取るような関数）をユーザが定義します。たとえばガウス関数$w(s_i, s_j)=exp(-\\|s_i - s_j\\|^2)$などを用いることができます。

さて、今回は2値分類問題を解きたいので$s_i$に対応する$x_i \in \\{1, -1\\}$という分類を求めることを考えます。$x_i$を一列に並べたベクトル${\bf x}$と置いて分類を最適化したあとでこの${\bf x}$を参照し、$1$になっている$s_i$を集めるとグループ1、$-1$になっている$s_i$を集めるとグループ2になっているよ、といった具合です。

ではどうやって最適化するかという話ですが、ここでグラフカット最小化を行います。すべての$s_i$を、グループ$G_1,G_2$に分類し終わった状態を考えてください（上の図）。分類し終わっているということは$G_1,G_2$の中にはそれぞれ似通った$s_i$が含まれているということであり、$G_1$の中から適当に選んだ2点$s_i,s_j$の類似度は大きくなっているはずです（図では2とか3）。逆に、異なるグループ$G_1,G_2$からそれぞれ1点ずつ選んだ$s_i, s_j$の類似度は小さくなっているはずです（図では0.1）。

したがってすべての点を2つのグループに分割し終えたとき、異なるグループから2点を選ぶすべての組み合わせについて類似度を足してみて、それが最小になっていればうまく2つのグループに分割できたと言えそうです（異なるグループから選んだ2点は極力似ていないということになる）。グラフ理論の言葉では、2つのグループをまたいでいるエッジをカットエッジといい、カットエッジの総和をカットのサイズと言います。

私たちは適当に選んだ2点の間の類似度の情報を、類似度行列$W$という形で完全に把握しています。これを用いると実はカットのサイズを表現するうまい方法があります。

```math
S = \sum_{i, j} W_{ij}(x_i - x_j)^2 \tag{2}
```

$x_i$と$x_j$は分類に応じて$1,-1$のいずれかの値を取る変数だったことを思い出してください。もし$s_i$と$s_j$が同じグループであったならば$x_i = x_j$となり、$(x_i-x_j)^2 = 0$となって対応する重み$W_{ij}$は$S$に加算されません。一方で$s_i$と$s_j$が異なるグループであったならば$(x_i - x_j)^2 = 4$となって対応する$W_{ij}$の4倍が$S$に加算されます。

すなわち$S$はカットのサイズの4倍の値を持つ関数なので、この$S$を最小化すればよいということになります。余談ですがこの$S$は[グラフ信号処理](https://www.jstage.jst.go.jp/article/essfr/8/1/8_15/_pdf)の分野で「信号のなだらかさ」を表す量です。一般のグラフ構造に対する畳込みニューラルネットワークである[Graph Convolutional Network](http://tech-blog.abeja.asia/entry/2017/04/27/105613)にも登場します。

## 量子アニーリングとスペクトラルクラスタリングの繋がり
さて、ここまでで量子アニーリングとスペクトラルクラスタリングのそれぞれについては大まかに解説しましたが、この2つがどう結びつくかというのが今回の本題です。今までの話からして「なんだ大体同じじゃん」って感じがしたと思います。勘の鋭い人ならば$(1)$式と$(2)$式を見ればその時点でビビッと来ると思います。私は勘が鋭いのでビビッと来ちゃったんですね。天啓？っていうんですかこれ。

すみません図に乗りました。行列の2次形式について知っていればなんとなく気付くと思います。

まずは量子アニーリングのほうから取りかかりましょう。$(1)$式で相互作用係数を表していた$J_{ij}$を${ij}$成分に持つ行列$J$を考え、さらに$J$の対角成分だけを抜き出した行列を$T$と置きます。クロネッカーのデルタを用いて表現すれば

```math
T_{ij} = J_{ij} \delta_{ij}
```

です。こうして定義した$T$を用いてさらに$P = T - J$となる行列$P$を用いれば、$(1)$式は

```math
H = {\bf \sigma}^{\mathrm{T}} P {\bf \sigma} \tag{3}
```

と表すことができます。頑張って展開してみてください。

続いてスペクトラルクラスタリングのほうに取りかかりましょう。まずは重みつきグラフに拡張された次数を対角要素に持つ次数行列$D$を定義します。

```math
D_{ii} = \sum_{j=1}^n W_{ij}
```

ただし$n$はデータ${\bf s}$の個数であり、$D$の対角成分以外は$0$です。グラフのノードの次数とはそのノードに接続しているエッジの本数を表しますが、重みつきグラフでは接続しているエッジの重みの総和として次数を定義することがあります。ここで$L = D- W$となる行列$L$を用いると、$(2)$式は

```math
S = 2{\bf x}^{\mathrm{T}} L {\bf x} \tag{4}
```

と表すことができます。こちらも頑張って展開すれば$(2)$式との一致を確かめることができます。ここで$L$は[スペクトラルグラフ理論](https://www.slideshare.net/irrrrr/ss-25911553)([参考](http://www.math.ucsd.edu/~fan/cbms.pdf))においてグラフラプラシアンと呼ばれる量です。スペクトラルクラスタリングがスペクトラルクラスタリングと呼ばれる所以がここにあるのかどうかは知りません。

さて、$(3)$式と$(4)$式を見比べていただければ分かると思いますが、この2つは定数倍の違いこそあれ一致しています。これらが最小化問題のコスト関数であることを考えれば定数倍の違いは無視できるので、量子アニーリングとスペクトラルクラスタリングは本質的に同じ問題を解いていたということになります。

厳密に言えば類似度行列$W$に関しては、類似度の対称性から$W_{ij}=W_{ji}$となっており、$L$は実対称行列になっている一方、イジングモデルの相互作用係数$J$に関しては$J_{ij}$と$J_{ji}$が一致するのかどうか私にはちょっとよくわからないので$P$が実対称行列になるかどうかはわかりません（どなたかご教授くださるとありがたいです）。また、類似度を考える場合は普通は「どれだけ似ているか」を与える正の値だけを考えるので、$W$の各成分がすべて正であるケースが多いですが、$J$に関しては各成分に負の値が入っても構いません。

とまぁそんなこんなありつつも、量子アニーリングとスペクトラルクラスタリングが行列の2次形式によって繋がることが確認できました。逆に言えば、問題を行列の2次形式の最小化問題に落とし込むことさえできれば、この2つの強力な手法で解ける可能性があるということです。量子アニーリングの将来が期待されますね。

## おまけ　スペクトラルクラスタリングの解法
量子アニーリングとスペクトラルクラスタリングが一致することに偶然にも気がついてしまったみなさんが誰しも抱く疑問があると思います。

「じゃあ量子アニーリングいらなくね？」
「スペクトラルクラスタリング解けないんじゃね？」

ええ、どちらもごもっともな疑問です。古典コンピュータのスペクトラルクラスタリングで解けるんならなんでわざわざ量子アニーリングマシンなんか作ったんだ、って話になりますよね。実を言うとスペクトラルクラスタリングは数値計算によって厳密解を求める方法がまだ発見されていないので量子アニーリングマシンにも需要がある、って感じです。

$(4)$式では${\bf x}$の各成分は$1,-1$のいずれかの整数値をとるという制約がありますが、これは「整数計画問題」と呼ばれていて一般に解くことが困難です。そこで${\bf x}$の各成分の範囲をとりあえず実数で計算してよい、という緩和問題を解くことにします。また都合上、類似度は常に正の値を持つとします（たとえば類似度にガウス関数を用いていればこの仮定は成立します）。

類似度が常に正であるという仮定より$(3)$式から$S \geq 0$となります。ここで少し考えてみると、$(4)$式には$\forall i \forall j, x_i = x_j$によって$S=0$となる自明な解が存在します。要するに全部同じグループにまとめてしまえば、カットエッジが存在しないのでカットのサイズも$0$になるというわけです。これは無意味な解ですね。

また各成分の値が$0$に近づく（${\bf x} \to {\bf 0}$）ときも同じ理由で$S \to 0$となります。したがって${\bf x}$の大きさに制約を設けない限り$S$はいくらでも下界$0$に近づきうるので、${\bf x}^{\mathrm{T}}D{\bf x}=1$となる制約条件を付け加えて、ラグランジュの未定乗数法によって$(4)$式を解くことにします（心配であれば${\bf x}^{\mathrm{T}}D{\bf x}=r$のように制約を一般化しても構いませんが、結果には影響がありません）。するとラグランジュ関数が

```math
\mathscr{L}({\bf x}, \lambda) = {\bf x}^{\mathrm{T}}L{\bf x} - \lambda({\bf x}^{\mathrm{T}}D{\bf x}-1)
```

となるような最適化問題を解くことになります。これを${\bf x}$で微分して${\bf 0}$とおくと

```math
L {\bf x} = \lambda D {\bf x}
```

という一般化固有値問題を解くことに帰着されます。ただし先ほど説明したように固有値$\lambda=0$で${\bf x}\propto {\bf 1}$という自明な解が存在するので、$0$でない最小の固有値に対応する固有ベクトル${\bf x}$がこの問題の解になります。

さて、そもそも解きたいのは$\forall i, x_i \in \\{1, -1\\}$という整数の制約がついた問題だったので、実数に緩和して求めた${\bf x}$の各成分を何らかの方法で$\\{1,-1\\}$の空間に写像してやる必要があります。うまくいきそうなのは

```math
f(x) = \begin{cases}
    1 & (x \geq a) \\
    -1 & (otherwise)
  \end{cases}
```

という閾値関数を用いる方法です。閾値$a$の決め方にはいろいろな方法があり、たとえばグループのメンバー数がなるべく偏らないように設定するとか、クラス内の分散がなるべく小さくなるようにといった方法が提案されているようです。

スペクトラルクラスタリングは[ラプラシアン固有マップ法](http://d.hatena.ne.jp/yag_ays/20110621/1308584676)（[参考](http://nlp.dse.ibaraki.ac.jp/~shinnou/zemi2012/kernel/kernel-kuramochi-0511.pdf)）というカーネル主成分分析による次元圧縮手法とも関連があります。関連があるというかほぼ同じです。詳しくは『[カーネル多変量解析](https://www.amazon.co.jp/%E3%82%AB%E3%83%BC%E3%83%8D%E3%83%AB%E5%A4%9A%E5%A4%89%E9%87%8F%E8%A7%A3%E6%9E%90%E2%80%95%E9%9D%9E%E7%B7%9A%E5%BD%A2%E3%83%87%E3%83%BC%E3%82%BF%E8%A7%A3%E6%9E%90%E3%81%AE%E6%96%B0%E3%81%97%E3%81%84%E5%B1%95%E9%96%8B-%E3%82%B7%E3%83%AA%E3%83%BC%E3%82%BA%E7%A2%BA%E7%8E%87%E3%81%A8%E6%83%85%E5%A0%B1%E3%81%AE%E7%A7%91%E5%AD%A6-%E8%B5%A4%E7%A9%82-%E6%98%AD%E5%A4%AA%E9%83%8E/dp/4000069713)』を読むとよいでしょう。

以上、お疲れさまでした。