原論文和譯及び解說 Geman-Gemanの定理 (シュミレーテッドアニーリングの収束定理)

原論文:http://image.diku.dk/imagecanon/material/GemanPAMI84.pdf

まずは論文概要や導入の概略を見てみましょう。原論文そのものは、劣化画像の復元のための手段としてシュミレーテッドアニーリングを捉えていることを念頭に置いておくと理解しやすいでしょう。原論文に数式などの不備や書き間違いが多少あるので、そこは修正しています。また、シュミレーテッドアニーリングそのものの理解にあまり必要ないと思われる章は飛ばしています。補足は、私が書いたものです。

 

1)直近の「近傍」のピクセルと境界要素の現在の値に基づいて、画像に局所的な変更が加えられます。この変更はランダムであり、ローカル条件付き確率分布からサンプリングすることにより生成されます。

2)局所条件付き分布は、「温度」と呼ばれるグローバル制御パラメータ\(T\)に依存します。低温では、局所条件付き分布は目的関数を増加させる状態に集中しますが、高温では分布は本質的に均一である。極限を考えると、\(T=0\)及び\(T=\infty\)は、それぞれ貪欲アルゴリズム(勾配上昇など)及び無向(つまり「純粋にランダム」)な変更に対応します。(高温は、隣接するピクセル間の疎結合と画像の混沌とした外観を引き起こします。低温では、結合はより緊密になり、画像はより規則的に見えます。)


3)我々の画像復元法は、確率的変化の多くが起こる高温で開始することにより、局所的な最大値を回避し、実際に目的関数を減少させます。緩和が進むと温度は徐々に低下し、過程は繰り返し改善されていくように振る舞います。(「焼き鈍し」をシュミレートしたこの温度の段階的低下は、特定の化学システムを低エネルギーの非常に規則的な状態に遷移させる過程です。)


我々の「アニーリング定理」は、事後分布のグローバルな最大値への収束を保証する温度を下げるためのスケジュールを規定しています。
このスケジュールは実際に適用するには遅すぎ、温度と時間の依存関係の関数形式を選択する際のガイドとしてのみ使用します。


また、私自身は引用論文は読んでいませんが、シュミレーテッドアニーリングの考え方の源流がどこにあるかも見て取れます。


温度と焼き鈍しのシュミレートを導入する発想は、Cerny[8]及び Kirkpatrick et al.[40]に依る。彼らはどちらも最適化問題に用いていて、Kirkpatrickはコンピュータ設計にも適用している。

[8] V. Cerny, “A thermodynamical approach to the traveling salesman problem: an efficient simulation algorithm,” preprint, Inst.
Phys. & Biophys., Comenius Univ., Bratislava, 1982.

[40] S. Kirkpatrick, C. D. Gellatt, Jr., and M. P. Vecchi, “Optimization by simulated annealing,” IBM Thomas J. Watson Research
Center, Yorktown Heights, NY, 1982.


いよいよ、定理の主張を見ていきます。

\(T\)が減少すると、事後分布からのサンプルは最小エネルギー構成に向かって強制されていきます。これらは、分布のモードに対応します。定理Bでこれを厳密にし、私たちの知る限り、この性質の最初の理論的結果です。大まかに言えば、\(k\)番目のサイト置換の実行に使用される温度\(T(k)\)(つまり、反復スキームの\(k\)番目のイメージ)が次の条件を満たす場合、

\[T(k)\geq \frac{c}{\log(1+k)}\]

すべての\(k\)について、\(c\)は\(k\)に依存しない定数であり、確率が1に収束する(\(k\rightarrow \infty\)とした)場合、アルゴリズムによって生成される構成は最小エネルギーの構成になります。
別の言い方をすれば、アルゴリズムは、最小エネルギー構成で分布が均一な尺度に収束するマルコフ連鎖を生成します。(確率1の収束は、一般に不可能であることを強調する必要があります。)

また、実際の収束率に関する定数\(c\)の性質についても議論します。

基本的に我々は、対数レートが最適であると考えます。


では定理の証明に関する話題に入っていきましょう。まずは数学的準備です。

3.グラフと近傍

\(S=\{s_1,s_2,…,s_N\}\)をサイトの集合とし、\(G=\{G_s,s \in S\}\)を\(S\)の近傍、つまり\(s \notin G_s,かつs\in G_r \Leftrightarrow r \in G_s\)を満たすような\(S\)の部分集合のそれぞれの集まりと定義する。

明らかに、\(G_s\)は\(s\)の近傍であり、\(\{S,G\}\)は普通の意味でのグラフを成す。

部分集合\(c  \subseteq S\)は、\(c\)のそれぞれのサイトの全てのペアが隣接している場合、クリークであり、\(C\)をクリーク\(c\)全ての集合とする。

 

4.マルコフ確率場とギブス分布

 

3章と同じように、\(\{S,G\}\)は任意のグラフを表す。\(X=\{X_s,s \in S\}\)を\(S\)を添え字集合としたランダムな値の族とする。

簡単のために、全ての\(s\)に対して\(X_s\in A\)なる\(A\equiv \{0,1,2,…,L-1\}\)のような共通の状態空間を想定する。

(補足:普通のシュミレーテッドアニーリングにおける、各サイトの取りうる値\(\{0,1\}\)に対応します。(quboの場合。量子ビットでは\(\{-1,1\}\)))

 

\(\Omega\)を全ての可能な構成の集合とする。

\[\Omega =\{ \omega = (x_{s_1},..,x_{s_N}):x_{s_i} \in A ,1 \leq i \leq N\}\]

(補足:普通のシュミレーテッドアニーリングでは、サイト数\(N\)のとき\(\Omega=\{(0,0,..,0),…(1,1,..,1)\}で、|\Omega|=2^N\)になりますね)

 

\(X\)は以下を満たすとき、\(G\)に関するMRF(マルコフ確率場)である。

全ての\(\omega \in \Omega\)に関して、\(P(X=\omega)>0\)

全ての\(S^{‘},S^{”} \subset S(r \neq s)\)と\((x_{s_1},..,x_{s_N})\in \Omega\)に関して、グラフ\(\{S,G\}\)上で\(B\)が\(S^{‘},S^{”}\)を分断している、つまり、\(S^{‘}\)から\(^{”}\)に行く全ての経路が\(B\)を通るとすると、

\(X(B)\)が与えられた時、\(X(S^{‘})\)と\(X(S^{”})\)は条件付き独立になる。

(補足:大域マルコフ性。上のは書き直してありますが、論文を見るときは論文の説明がまずいので、ウィキペディアの方を参考にして下さい。)

 

\(\{S,G\}\)に関するギブス分布は、以下の表式を持つ\(\Omega\)上の確率測度である。

\[\pi(\omega)=\frac{1}{Z}e^{-U(\omega)/T}\]

\(T\)は定数であり、\(U\)はエネルギー関数である。\(Z\)は規格化定数であり、

\[Z=\sum_{\omega}e^{-U(\omega)/T}\]

 分配関数と呼ばれる。

 

11.ギブスサンプラー:概要

 

一般的な計算上の問題は、

A)分布\(\pi\)からサンプルする

B)\(\Omega\)上で\(U\)を最小にする

C)期待される値を計算する

である。

 

時刻\(t\)で、\(X_{s_t}\)は\(A\)のいずれかの値である。各更新では、一つのサイトのみが値を変えることができる。つまり、\(X(t-1)\)と\(X(t)\)では一か所のサイトしか変化していない。\(n_1,n_2,..\)を、どのサイトが変化するために「訪問」されたかを示す数列とする。つまり、\(n_t \in S\)であり、\(i \neq n_t\)の時\(X_{s_i}(t)=X_{s_i}(t-1)\)である。

 

12.ギブスサンプラー:数学的基礎

 

定理A(省略)

一様マルコフ連鎖の収束性について述べています。ペロン・フロベニウスの定理の証明をもう一回やっているだけです。別の文献がいくらでもあります。

 

定理Bは、「アニーリングスケジュール」または温度低下率であり、システムを強制的に最低エネルギー状態にする。

\(\pi\)の\(T\)依存性を\(\pi_T\)と書き、\(T(t)\)をステージ\(t\)での温度とする。

最小エネルギー状態\(\Omega_0=\{\omega \in \Omega:U(\omega)=\min _{\eta}U(\eta)\}\)とし、

\(\pi_0\)を\(\Omega_0\)上の一様分布とする。

(補足:最小エネルギー状態が\(n\)個あるとき、最小エネルギー状態の取りうる確率が\(1/n\),それ以外が0になるような確率測度)

最後に、

\(U^*=\max_{\omega}U(\omega)\)

\(U_*=\min_{\omega}U(\omega)\)

\(\Delta=U^*-U_*\)

と定義する。

 

定理B(アニーリング)

全ての\(t=0,1,2…\)で\(S \subset \{n_{t+1},n_{t+2},…,n_{t+\tau}\}\)となる\(\tau>N\)の存在を仮定する。

(補足:そういう風にマルコフ連鎖を構成すればいいだけなので、容易に前提条件は満たせます。例えばそれぞれのサイトを順番に回るなど。)

 

\(T(t)\)を、以下を満たすような温度の減少列とする。

a)\(t \rightarrow \infty\)で\(T(t)\rightarrow 0\)

b)\(t_0\geq 2\)なる適当な整数に関して、全ての\(t \geq t_0\)で\(T(t)\geq N\Delta/\log t\)

を満たすならば、

任意の初期分布\(\eta\in \Omega\)と全ての\(\omega \in \Omega\)に関して、

\(\lim_{t\rightarrow \infty }P(X(t)=\omega|X(0)=\eta)=\pi_0(\omega)\)

 

定理C(エルゴード性)

省略。論文を参照してください。

 

定理B-証明

 

まずは二つの補題を示す。

(補足:非一様マルコフ連鎖の収束性について述べています。シュミレーテッドアニーリングや経路積分モンテカルロ法による量子アニーリング、Green関数モンテカルロ法による量子アニーリングなど、古典的に量子アニーリングをシュミレートする、本質的に非一様マルコフ連鎖と関連するアルゴリズムの全ての収束性について関係があります。量子アニーリングの収束定理参照。)

 

補題2

全ての\(t_0=0,1,2,…,\)に対して、

\(\lim_{t\rightarrow \infty}\sup _{\omega,\eta^\prime,\eta^{\prime \prime}}|P(X(t)=\omega|X(t_0)=\eta^\prime )-P(X(t)=\omega|X(t_0)=\eta^{\prime \prime}))|=0\)

 (補足:補題2と補題3を証明すれば、定理Bが証明されます。補題2の必要条件として定理Bの条件bが得られます。補題3の必要条件として、条件aが使われます。)

 

補題3

||は\(L_1\)距離を表す。

\(\lim_{t_0\rightarrow \infty}\sup_{t\geq t_0}||P(t|t_0,\eta)-\pi_0||=0\)

 

まず、これらの補題を如何に定理Bに適用するかを示す。\(\sup\)は\(\omega,\eta\)についてとる。

\(\limsup_{t\rightarrow \infty}||P(X(t)=\omega|X(0)=\eta)-\pi_0(\omega)||\)

\(=\limsup _{t_0 \rightarrow \infty }\limsup_{t\rightarrow \infty,t\geq t_0}\)

\(||\sum_{\eta^\prime} P(X(t)=\omega|X(t_0)=\eta) \cdot P(X(t_0)=\eta|X(0)=\omega)-\pi_0||\)

\(\leq \limsup _{t_0 \rightarrow \infty }\limsup_{t\rightarrow \infty,t\geq t_0}\)

\(||\sum_{\eta^\prime} P(X(t)=\omega|X(t_0)=\eta^\prime) \cdot P(X(t_0)=\eta^\prime|X(0)=\omega)-P(X(t)=\omega)|t_0,\pi_0)||\)

\(+||P(X(t)=\omega|X(t_0)=\pi_0)-\pi_0||\)

第二項は、補題3から0になる。第一項について、

 \(||\sum_{\eta^\prime} P(X(t)=\omega|X(t_0)=\eta) \cdot P(X(t_0)=\eta|X(0)=\omega))-P(X(t)=\omega)|t_0,\pi_0)||\)

\(=\sum_\omega \sup_{\eta^{\prime \prime}}|\sum_{\eta^\prime}(P(X(t)=\omega|X(t_0)=\eta^\prime) -P(X(t)=\omega|X(t_0)=\eta^{\prime\prime}) )\)

\(\times (P(X(t)=\eta^\prime|X(t_0)=\eta) -\pi_0(\eta^\prime))|\)

\(\leq 2\sum_\omega \sup_{\eta^\prime,\eta^{\prime\prime}}|P(X(t)=\omega|X(t_0)=\eta^\prime) -P(X(t)=\omega|X(t_0)=\eta^{\prime\prime}) |\)

よって、第一項も補題2から0になるから、補題2,3が成り立てば定理Bを示せることはわかった。

 

補題2の証明

\(\delta (t)=\inf_{1\leq i \leq N,(x_{s_1},…,x_{s_N})\in \Omega} \pi_{T(t)}(X_{s_i}=x_{s_i}|X_{s_j}=x_{s_j},i \neq j)\)

と定義する。

\(\delta (t)\geq \frac{e^{-U^*/T(t)}}{{e^{-U_*/T(t)}}}=e^{-\Delta/T(t)}\)

であることは容易に示せる。

(補足:ベイズの定理を用いる。論文では分母に\(L\)が来ていますが、必要ありません。どちらにせよ結果に影響はありません。)

 

今、\(k\)を固定して考え、\(m_i\)を以下のように定義する。

\(m_i=\sup \{t:t \leq T_k ,n_t=s_i\}\ ,\ 1\leq i \leq N\)

(補足:\(k\)は直観的に言えば、マルコフ連鎖で全てのサイトを今まで何周したかを表すものであることを思い出してください。つまり\(m_i\)は、\(k\)周目で”訪問”するサイト\(i\)の最後に”訪問”する時刻を表しています。)

 また、\(m_1>m_2>m_3>\cdots >m_N\)と仮定する。

(補足:そうなるように最初からマルコフ連鎖を構成すればよいだけです。逆に言えば、シュミレーテッドアニーリングのマルコフ連鎖を構成するときは、各サイトを順番に周って値を更新する必要があることを言っています。)

このとき、

\(P(X(T_k)=\omega|X(t_{k-1})=\omega^\prime)=P(X_{s_1}(m_1)=x_{s_1},\cdots X_{s_N}(m_N)=x_{s_N}|X(t_{k-1})=\omega^\prime)\)

\(=\prod _{j=1}^N P(X_{s_j}(m_j)=x_{s_j}|X_{s_{j+1}}(m_{j+1})=x_{s_{j+1}},\cdots ,X_{s_N}(m_N)=x_{s_N},X(T_{k-1})=\omega^\prime)\)(条件付き分布の定義から)

\(\geq \prod_{j=1}^N \delta(m_j)\)(\(\delta\) の定義から)

\(\geq  \prod_{j=1}^N e^{-\Delta/T(m_j)}\)

\(\geq \exp\{-\Delta N/(t_0+k\tau)\}\)(\(m_j \leq T_k=t_0+k\tau ,j=1,2,…,N\)であり、かつ\(Tは\)減少関数であるから)

\(\geq (t_0+k\tau)^{-1}\)(定理の仮定Bより。逆に言えば、以下の議論を追えば、このアニーリングスケジュールがシュミレーテッドアニーリングの収束の必要条件になることがわかります。)

 

\(t>t_0\)なる\(t\)に対して、\(K(t)\)を

\(K(t)=\sup\{k:T_k<t\}\)と定義する。

(補足:時刻\(t\)ですでに各サイトを何周回ったか)

 

\(\sup _{\omega,\eta^\prime,\eta^{\prime\prime}}|P(X(T_k)=\omega|X(t_0)=\eta^\prime)-P(X(T_k)=\omega|X(t_0)=\eta^{\prime\prime})|\)

\(\leq \prod_{k=1}^{K(t)} (1-(t_0+k\tau)^{-1})\)

 

\(t\rightarrow \infty \Rightarrow K(t) \rightarrow \infty\)であることと、

\(\prod_{k=1}^m (1-(t_0+k\tau)^{-1})\)

\(=\prod_{k=1}^m \exp \log (1-(t_0+k\tau)^{-1})\)

\(=\exp(\sum_{k=1}^m \log (1-(t_0+k\tau)^{-1}))\)

\(\leq \exp(\sum_{k=1}^m (t_0+k\tau)^{-1})\)(\(x>0\)で\(\log (1-x)<-x\)であることから)

\(\rightarrow 0(m \rightarrow \infty)\)

 

これより、補題2は示された。

 

補題3の証明

\(P(t,\omega|t_0,\pi_0)\)

\(=\sum_\eta P(t,\omega|X(t_0)=\eta)\pi_0 (\eta)\)

(補足:初期状態が\(\pi\)によって決定された後、時刻\(t\)で\(\omega\)である確率。証明を構成するためにこのような確率を考えますが、結局初期状態が何であっても収束することが示せればよいので、初期状態の選び方そのものは特殊な状況を考えて問題ないということを言っています。)

 

 \(P(t,\omega|t_0,\pi_0)\)を以下、\(P_{t_0,t}\)と略記する。

始めに、\(t>t_0\geq 0\)に対し、以下を示す。

\(||P_{t_0,t}-\pi_{T(t)}|| \leq ||P_{t_0,t-1}-\pi_{T(t)}||\)

 

簡単のために\(n_t=s_1\)として、

\(||P_{t_0,t}-\pi_{T(t)}||\)

\(=\sum _{(x_{s_1},…,x_{s_N})} |\pi_{T(t)}(X_{s_1}=x_{s_1}|X_s=x_s,s \neq s_1) \cdot P_{t_0,t-1}(X_s=x_s,s \neq s_1)-\pi_{T(t)} (X_s=x_s,s \in S)|\)

\(=\sum _{(x_{s_2},…,x_{s_N})}|P_{t_0,t-1}(X_s=x_s,s \neq s_1)-\pi_{T(t)} (X_s=x_s,s \in S)|\)

\(\leq \sum_{(x_{s_1},…,x_{s_N}) \in \Omega} |P_{t_0,t-1}(X_s=x_s,s \neq s_1)-\pi_{T(t)} (X_s=x_s,s \in S)|\)

\(=||P_{t_0,t-1}-\pi_{T(t)}||\)

 

次に、\(t \rightarrow \infty\)で\(||\pi_0-\pi_{T(t)}||\rightarrow 0\)を示す。

\(\pi_{T(t)} (\omega)=\frac{\exp(-U(\omega)/T(t))}{\sum_{\omega^\prime \in \Omega_0}\exp(-U(\omega)/T(t)+\sum_{\omega^\prime \in \Omega \backslash \Omega_0}\exp(-U(\omega)/T(t)}\)

\(=\frac{\exp(-(U(\omega)-U_*)/T(t))}{|\Omega_0|+\sum_{\omega^\prime \in \Omega \backslash \Omega_0}\exp(-(U(\omega)-U_*)/T(t)}\)

よって、\(t \rightarrow \infty\)の極限で、\(\pi_{T(t)} (\omega)\)は\(\Omega_0\)上の一様分布に近づく。

 

(原論文のこれ以降の説明は、ある一定より先は単調減少して0に近づく数列の和が有界であるなどと主張している(簡単な反例が挙げられる)ことや、\(\pi_T(t)\)の時間依存性についての議論がデタラメなので、修正しています。)

 

\(t>t_0\geq 0\)において、

\(||P_{t_0,t}-\pi_0||\leq||P_{t_0,t}-\pi_{T(t)}||+||\pi_{T(t)}-\pi_0|| \)

\(\leq ||P_{t_0,t_0}-\pi_T(t)||+||\pi_{T(t)}-\pi_0||\)(\(||P_{t_0,t}-\pi_{T(t)}|| \leq ||P_{t_0,t-1}-\pi_{T(t)}||\)であるから)

\(=||P_{t_0,t_0}-\pi_{T(t_0)}||+||\pi_{T(t_0)}-\pi_{T(t)}||+||\pi_{T(t)}-\pi_0||\)

 ここで、\(P_{t_0,t_0}=\pi_0\)であり、\(t\rightarrow \infty\)で\(||\pi_{T(t)}-\pi_0||\rightarrow 0\)であるから、上式の\(\lim_{t_0\rightarrow \infty}\sup_{t\geq t_0}\)をとれば良い。

ホップの定理(Hopf’s theorem)-行列

ページランクなどの、本質的にマルコフ連鎖に関連するもの(要は、べき乗法)の収束するかしないかの判定法として、ペロン・フロベニウスの定理(Perron-Frobenius theorem)があります。

しかし、この定理が言及するのは、あくまで\(n \rightarrow \infty\)の極限のみであり、収束の速さについては述べません。

べき乗法の収束の速さについて、行列の固有値の絶対値の大きさが本質的に関わって来ることはそこらへんにいくらでも書いてありますし、ちょっと考えればわかることなのでここには書きませんが、以下に示すHopfの定理を用いれば正行列に対し、べき乗法の収束の速さの下限を求めることができます。(べき乗法と固有値にについては、Wikipediaにも書いてあります。)

直観的に言えば、マルコフ過程の遷移行列が疎になる(0成分が多くなる)ほど収束が早いか遅いか判定しづらくなると言うことを述べたのが、以下、ホップの定理(Hopf’s theorem)です。

(\(\kappa \rightarrow \infty\)を考えてみてください。つまり残念ながら、定理の適用可能な行列のクラスは、厳密に正な正方行列に限られます。)

 

ホップの定理(Hopf’s theorem)

正方行列\(M\)の全ての要素が正定値の時、\(M\)の最大固有値\(\lambda_0\)と他の任意の固有値\(\lambda\)は、

\[|\lambda|\leq\frac{\kappa-1}{\kappa+1}\lambda_0 \]

ここで、\(\kappa\)は

\[\displaystyle \kappa=\max_{i,j,k}\frac{M_{ik}}{M_{jk}}\]

 

結構重要な定理だと思うので、日本語でまとめました。ペロン・フロベニウスの定理の証明について知っている方は、以下の定理1,定理4の部分だけ見るだけでもいいと思います。

ちなみに、ホップはもともと、一般に正の積分作用素についてこの定理を証明したようですが、ここでは正行列に注目して考えます。

量子アニーリングの数学的基礎 森田 悟史 西森秀稔 (Mathemathical Foundation of Quantum Annealing)という論文を参考にしています。(日本人が書いた英語の論文を、日本人が見ると言う…)


\(M\)を、厳密に正な(すべての成分が正の)\(m \times m\)行列とする。これを\(M>0\)と書く。

同様に、\(M \geq 0\)は、全ての成分が非負であることを意味する。

ベクトルにも同じ記法を用い、\(\boldsymbol{v} >0\)は\(\boldsymbol{v}\)の全成分が正であることを意味する。

\(M\)の厳密な正定値性は、以下と同値。

\(M\boldsymbol{v}>0 \ \ \ \ if\ \ \ \ \boldsymbol{v} \geq 0,\boldsymbol{v} \neq \boldsymbol{0} \tag{1}\)

 

あらゆる実ベクトル\(\boldsymbol{v}\)と正ベクトル\(\boldsymbol{p}>0\)は、以下を満たす。

\[\displaystyle \min_{i} \frac{v_i}{p_i} \leq \min_i \frac{(M \boldsymbol{v})_i}{(M \boldsymbol{p})_i} \leq \max_i \frac{(M \boldsymbol{v})_i}{(M \boldsymbol{p})_i} \leq \max_i \frac{v_i}{p_i}\]

ならば、

\[(M \boldsymbol{v})_i-( \min _i \frac{v_i}{p_i})(M \boldsymbol{p})_i  =\sum _{j=1}^m M_{ij} [v_j -(\min _i \frac{v_i}{p_i})p_j] \geq 0\]

\(\max\)の方も同様。

 

以下の記号を用いる。

\[\displaystyle \underset{i} {osc} \frac{v_i}{p_i}=\max _i \frac{v_i}{p_i}-\min_i \frac{v_i}{p_i}\]

これは振幅(oscillation)と呼ばれる。

複素ベクトルに対しては、以下の様に定義する。

\[\displaystyle \underset{i}{osc}\ v_i=\sup_{|\eta|=1}\underset{i}{osc}\ Re(\eta v_i)\]

あらゆる複素数\(c\)に対し、

\[\underset{i}{osc}\ v_i=|c|\underset{i}{osc}\ v_i\]

であることがわかる。 

 もし\({osc}_i\ v_i=0\)であるならば、\(v_i\)は\(i\)に依らないことも簡単にわかる。

 

行列\(M\)に対し、以下の仮定を置く。

\[\forall i,j,k\ \ ,\exists \kappa ,\ \frac{M_{ik}}{M_{jk}}\leq  \kappa\]

以下の様に書き換えることもできる。

\[\forall i,j, \boldsymbol{v}\ (\boldsymbol{v}\geq 0,\ \boldsymbol{v} \neq 0),\ \ \frac{(M\boldsymbol{v})_i}{(M\boldsymbol{v})_j}\leq \kappa,\ \  \tag{2} \]


定理1.

もし\(M\)が 条件(1),(2)を満たすならば、あらゆる\(\boldsymbol{p}>0\)と複素ベクトル\(\boldsymbol{v}\)に対し、以下が成り立つ。

 

\[\underset{i}{osc}\  \frac{(M\boldsymbol{v})_i}{(M\boldsymbol{p})_i} \leq \frac{\kappa-1}{\kappa+1} \underset{i}{osc}\  \frac{v_i}{p_i} \tag{3}\]


証)まず、実ベクトル\(\boldsymbol{v}\)について考える。\(i,j,\boldsymbol{p}>0\)を固定し、実数\(X_k\)を以下に定義する。

\[\displaystyle \frac{(M\boldsymbol{v})_i}{(M\boldsymbol{p})_i}-\frac{(M\boldsymbol{v})_j}{(M\boldsymbol{p})_j}=\sum_{k=1}^{m} X_k(v_k)\]

\(X_k=X_k (i,j,\boldsymbol{p})\)の正確な表式は知る必要はない。

\(\boldsymbol{v}=a\boldsymbol{p}\)のとき、上式の左辺は消えるから、

\[\displaystyle \frac{(M\boldsymbol{v})_i}{(M\boldsymbol{p})_i}-\frac{(M\boldsymbol{v})_j}{(M\boldsymbol{p})_j}=\sum_{k=1}^{m} X_k(v_k-ap_k) \tag{4}\]

ここで

\[\displaystyle a=\min_i\frac{v_i}{p_i},b=\max_i \frac{v_i}{p_i}\]

のように選ぶと、

\(v_k-ap_k=(b-a)p_k-(bp_k-v_k)\)であるから、\(v_k-ap_k\)は\(v_k=ap_k\)で最小値\(0,v_k=bp_k\)で最大値\((b-a)p_k\)をとる。

それゆえ(4)式は、

\[\boldsymbol{v}=a\boldsymbol{p}^- +b\boldsymbol{p}^+=a\boldsymbol{p}+(b-a)\boldsymbol{p}^+\]

で最大値をとる。

(論文中(222)式の\(a\boldsymbol{p}^- -b\boldsymbol{p}^+\)は間違い)

ここで、

\[p_i^- =p_i(X_i\leq0),0(X_i>0)\]

\[p_i^+ =0(X_i \leq 0),p_i(X_i>0)\]

のように定義した。

結果、

\[\displaystyle \frac{(M\boldsymbol{v})_i}{(M\boldsymbol{p})_i}-\frac{(M\boldsymbol{v})_j}{(M\boldsymbol{p})_j} \leq \left[\displaystyle \frac{(M\boldsymbol{p}^+)_i}{(M\boldsymbol{p})_i}-\frac{(M\boldsymbol{p}^+)_j}{(M\boldsymbol{p})_j}\right](b-a) \tag{5}\]

が得られる。

仮定より\(M>0\)かつ\(\boldsymbol{p}>0\)であるから、

\[M\boldsymbol{p}^- \geq 0,\ \ M\boldsymbol{p}^+ \geq 0,\ \ M\boldsymbol{p}^- +M\boldsymbol{p}^+ >0\]

さらに、\(\boldsymbol{p}^- =0\)或いは\(\boldsymbol{p}^+ =0\)であるとき、\(\boldsymbol{p}^+=\boldsymbol{p}\)または\(0\)となるから、右辺の括弧内が消えてしまう。

よって、\(M\boldsymbol{p}^->0\)かつ\(M\boldsymbol{p}^+ >0\)であるとして議論を進めてよい。(5)式は、以下の様に書き直せる。

\[\displaystyle \frac{(M\boldsymbol{v})_i}{(M\boldsymbol{p})_i}-\frac{(M\boldsymbol{v})_j}{(M\boldsymbol{p})_j} \leq \frac{1}{1+t}-\frac{1}{1+t’},t\equiv \frac{(M\boldsymbol{p}^-)_i}{(M\boldsymbol{p}^+)_i},t’\equiv \frac{(M\boldsymbol{p}^+)_i}{(M\boldsymbol{p}^-)_i}\]

式(2)より、\(t\)と\(t’\)は\(\kappa^{-1}\)から\(\kappa\)で抑えられるから、\(t’ \leq \kappa^2 t\)であることがわかる。

(論文中,\(t’ \leq \kappa t^2\)は間違い)

\[\frac{(M\boldsymbol{v})_i}{(M\boldsymbol{p})_i}-\frac{(M\boldsymbol{v})_j}{(M\boldsymbol{p})_j} \leq \frac{1}{1+t}-\frac{1}{1+\kappa^2 t}\]

\(t>0\)において、上式の右辺は\(t=\kappa^{-1}\)において最大値\(\frac{\kappa-1}{\kappa+1}\)をとる。

最終的に、以下を得る。

\[\frac{(M\boldsymbol{v})_i}{(M\boldsymbol{p})_i}-\frac{(M\boldsymbol{v})_j}{(M\boldsymbol{p})_j} \leq \frac{\kappa-1}{\kappa+1} \underset{i}{osc}\ \frac{v_i}{p_i}\]

複素ベクトル\(\boldsymbol{v}\)に対しては、\(M\ Re (\eta \boldsymbol{v})=Re (\eta M \boldsymbol{v})\)が成り立つから、\(v_i\)を\(Re(\eta v_i)\)と置き換え、\(|\eta|=1\)の条件の下、両辺のsupをとることで以下を得る。

 \[\underset{i}{osc}\ Re\left(\eta\frac{(M\boldsymbol{v})_i}{(M\boldsymbol{p})_i}\right) \leq \frac{\kappa-1}{\kappa+1} \underset{i}{osc}\ Re\left(\eta\frac{v_i}{p_i}\right)\]


この定理を固有値問題に適用する。

\[M\boldsymbol{v}=\lambda \boldsymbol{v} \tag{6}\]

ペロン・フロベニウスの定理により、非負な正行列\(M\geq 0\)に対し、他のあらゆる固有値\(\lambda\)に関して、\(|\lambda|\leq \lambda_0\)となる。この結果は厳密に正な行列\(M>0\)にも使え、以下の定理を示す。


定理2

式(1),(2)の仮定の下、固有値問題(6)は正の解\(\lambda=\lambda_0>0\),\(\boldsymbol{v}=\boldsymbol{q}>0\)を持ち、さらにどんな\(\boldsymbol{p}(\boldsymbol{p}\geq 0,\boldsymbol{p}\neq 0)\)に対しても、

\[\boldsymbol{q}_n=\frac{M^n\boldsymbol{p}}{(M^n\boldsymbol{p})_k}\]

は\(k\)が固定された下、\(\exists \boldsymbol{q}\)に収束する。


証)

非負で0ではないベクトル\(\boldsymbol{p},\boldsymbol{\bar{p}}\)を考えよう。また、以下の様に定義する。

\[\boldsymbol{p}_{n+1}=M\boldsymbol{p}_n,\ \ \boldsymbol{\bar{p}}_{n+1}=M\boldsymbol{\bar{p}}_n,\ \ \boldsymbol{p}_0=\boldsymbol{p},\ \ \boldsymbol{\bar{p}}_0=\boldsymbol{\bar{p}}\]

仮定(1)から、\(\boldsymbol{p},\boldsymbol{\bar{p}}_n\)はともに、正である。定理1を繰り返し適用することにより、\(n>1\)で以下を得る。

\[\underset{i}{osc}\ \frac{\bar{p}_{n,i}}{p_{n,i}}\leq \left(\frac{\kappa-1}{\kappa+1}\right)^{n-1}\underset{i}{osc}\ \frac{\bar{p}_{1,i}}{p_{1,i}}\ tag{7}\]

 ここで、\(p_{n,i}\equiv(\boldsymbol{p}_n)_i\)という記法を使った。結果として

\[\exists \lambda>0,\forall i,\lim _{n\rightarrow \infty}\frac{\bar{p}_{n,i}}{p_{n,i}}= \lambda \tag{8}\]

がわかる。

ベクトル\(\boldsymbol{p}_n,\boldsymbol{\bar{p}}_n\)を以下の様に規格化する。

\[\boldsymbol{q}_n=\frac{\boldsymbol{p}_n}{p_{n,k}},\ \ \boldsymbol{\bar{q}}_n=\frac{\boldsymbol{\bar{p}}_n}{p_{n,k}}\]

仮定(2)より、

 \[\kappa^{-1}\leq q_{n,i}\leq\kappa,\ \ \kappa^{-1}\leq \bar{q}_{n,i}\leq \kappa \tag{9}\]

が言える。

それゆえ、以下が言える。

\[|\bar{q}_{n,i}-q_{n,i}|=q_{n,i}\frac{p_{n,k}}{\bar{p}_{n,k}} \left|\frac{\bar{p}_{n,i}}{p_{n,i}}-\frac{\bar{p}_{n,k}}{p_{n,k}}\right| \leq \kappa \frac{p_{n,k}}{\bar{p}_{n,k}} \underset{i}{osc}\ \frac{\bar{p}_{n,i}}{p_{n,i}} \tag{10}\]

ここで特に、\(\boldsymbol{\bar{p}}=M\boldsymbol{p}=\boldsymbol{p}_1\)である場合、

\[\boldsymbol{\bar{p}}_n=M\boldsymbol{p}_n=\boldsymbol{p}_{n+1},\ \ \boldsymbol{\bar{q}}_n=\boldsymbol{q}_{n+1}\]

式(7),(8)を用い、(10)により\(\boldsymbol{q}_n\)が\(\boldsymbol{q}\)へと収束することが言える。(9)より、\(\boldsymbol{q}>0\)である。(8)から、

\[\frac{p_{n+1,i}}{p_{n,i}}=\frac{(M\boldsymbol{p}_n)_i}{(\boldsymbol{p}_n)_i}=\frac{(M\boldsymbol{q}_n)_i}{(\boldsymbol{q}_n)_i}\rightarrow \lambda_0\]

 よって、

\[M\boldsymbol{q}=\lambda_0 \boldsymbol{q}\]

となる。定理2は示された。


定理3

同様の仮定の下、固有値問題(6)は\(\lambda=\lambda_0,\boldsymbol{v}=c\boldsymbol{q}\)を除き、正の解をもたない。\(\lambda=\lambda_0\)のとき、(6)は\(\boldsymbol{v}=c\boldsymbol{q}\)を除く解はない。


証)

\(\boldsymbol{v}\geq 0,\boldsymbol{v}\neq 0\)は固有値問題(6)の解とする。仮定(1)より、\(M\boldsymbol{v}>0\)であり、\(\lambda>0\)かつ\(\boldsymbol{v}>0\)なる解を持つ。これを定理2における初期ベクトル\(p\)として用い、定理の最後の部分を適用すると、

\[\frac{M^n\boldsymbol{v}}{(M^n\boldsymbol{v})_k}=\frac{\lambda^n\boldsymbol{v}}{\lambda^n\boldsymbol{v}_k}=\frac{\boldsymbol{v}}{v_k}\]

それゆえ、極限\(\boldsymbol{q}\)は\(\boldsymbol{v}/v_k\)と等しくなり、つまり\(\boldsymbol{v}=c\boldsymbol{q},\lambda=\lambda_0\)である。


定理4

同様の条件の下、固有値問題(6)のあらゆる固有値\(\lambda\neq\lambda_0\)は、以下を満たす。

\[|\lambda|\leq\frac{\kappa-1}{\kappa+1}\lambda_0\]

補足:\((\kappa-1)/(\kappa+1)\)は他の条件がない限り、最良の見積もりである。例えば、

\[M=\begin{pmatrix}\kappa & 1 \\ 1 &\kappa\end{pmatrix}\ \ \kappa>0\]

は固有値\(\lambda_0=\kappa+1,\lambda=\kappa-1\)を持つ。


 証)

次に、定理2の\(\lambda_0>0\)と\(\boldsymbol{q}>0\)なる解と\(M\boldsymbol{v}=\lambda\boldsymbol{v}\)の解を考える。\(\boldsymbol{q},\boldsymbol{v}\)に定理1を適用して、

\[\frac{|\lambda|}{\lambda_0}\underset{i}{osc}\ \frac{v_i}{q_i}=\underset{i}{osc}\ \frac{\lambda v_i}{\lambda_0q_i}=\underset{i}{osc}\ \frac{(M\boldsymbol{v})_i}{(M\boldsymbol{q})_i}\leq\frac{\kappa-1}{\kappa+1}\underset{i}{osc}\ \frac{v_i}{q_i}\]

\(\lambda\neq\lambda_0\)かつ\(\boldsymbol{v}\neq 0\)ならば\(\boldsymbol{v}\)は\(c\boldsymbol{q}\)と等しく成り得ないから、定理を得る。

【論考】グループ分けページランク

以下の内容は、転載を禁じます。

【導入】

ページランクとそれに類するアルゴリズムは一般に、ランク付けすべき対象(本論考では、ノードと呼ぶ)は同質な集団である。(例えば、ページランクの場合はウェブページ、EigenTrustではブロックチェーンに参加しているノード)

 

本論考では、ノードが本質的に違うようないくつかの集団に分けられるが、お互いに関係があって相対的な評価が可能である場合に、グループ分けしてランクを評価したい、あるいは、何らかの属性を持っているノードをクラスタリングし、グループ分けして異なるグループに属するノード間の相対的な評価を何らかの方法で定義した場合に、クラスタリングにより分けられたグループのそれぞれのノードの中でのみランクを評価したい場合などを考える。

例えば、店舗網があり、複数のアルバイトがいて、それぞれのアルバイトはいくつか複数の店舗で働いた経験をもっているような場合を考えると、店舗網同士、アルバイト同士での相対的評価値は存在しないが、各店舗と各アルバイトでの相対的評価値は何らかの方法であらかじめ計算しておけるような状況が考えられる。

また、会社などの組織中で、取締役、エンジニア、営業などの本質的に異なるグループがあり、取締役同士で直接評価値を出したくはないが、エンジニア、営業などからどれぐらい評価値があるのかを考えたうえで取締役同士の評価をしたい場合なども考えられる。

さらに、ある文章があり、それぞれの文ごとに区切り、これを一グループとして、それぞれの文中の単語のそれぞれにスコアを与えたい場合なども考えられるだろう。

 

以下、このような状況を数学的に抽象化、一般化、定式化して議論を行う。

 

【数学的議論】

 

グループの集合\(\textbf{G}=\{G_1,G_2,…,G_n\}\)

ノードの集合\(\textbf{n}=\{n_1,n_2,…,n_N \}(\forall n \in \textbf{n},\exists ! G \in \textbf{G},n \in G)\)

とし、ノード\(i\)から\(j\)への相対評価を\(r_{ij}\)と書く。

全てのグループの数\(n(=|\textbf{G}|)\)を長さとする任意の巡回置換\(\sigma\)を考え、全てのグループ\(G_i\)に関して、それに属する全てのノードからあるグループ\(G_{\sigma (i)}\)に属する全てのノードのみへの\(|G_i||G_{\sigma(i)}|\)通りの正の相対評価値\(\forall n \in G_i,m \in G_{\sigma(i)}\)に対して\(r_{nm}\)があらかじめ何らかの方法で計算されている場合を考える。

\(\forall i \in \{1,2,…,n\},\forall n \in G_i,m \in G_{\sigma(i)},\exists r_{nm} \in \mathbb{R}_{>0},\)

\(if \ n \in G_i,m \not\in G_ {\sigma(i)},r_{nm}=0\)

逆に、巡回置換でない場合、以下の議論を追えばわかるように、相対評価を行列で表現したものは既約でないため、べき乗法による収束は保証されない。

巡回置換である場合には、以下に示すように数学的、アルゴリズム的に面白い。さらに以下で示すように、複数の巡回置換を考えて行列に足しこむことで理論は一般化できるから、本論考ではこの場合のみに絞って考えていく。

 

ここで、以下のようにそれぞれのノードを評価する方法を考える。

 まず、相対評価を表す行列\(r_{ij}\)を行列で表したもの)\(A\)は、行、列の入れ替え(つまり、\(\sigma(i)=i+1\)になるようにグループ番号を変更し、行列の各行、列に早いグループ番号に属するノードから詰めていく。)によって以下の形の行列に変形できる。

\(A=\begin{pmatrix}0 & A_1 & 0 & \cdots &0\\ 0 & 0 & A_2 & 0 &0\\ \vdots & \vdots & \vdots  &  \ddots &\\ A_n & 0 & 0 & 0 &0\end{pmatrix}\)

\(A_i\)は、グループ\(i\)の属するノードからグループ\(i+1\)に属するノードへの評価を表した行列である。

 \(A^2\)は、

\(A^2=\begin{pmatrix}0 & 0 & A_1 A_2 & \cdots &0\\ 0 & 0 & 0 & \cdots &0\\ \vdots & \vdots & \vdots  &  \ddots &\\ 0 & A_n A_1 & 0 & 0 &0\end{pmatrix}\)

となる。

これを繰り返す(n-1回の掛け算を行う)と、以下の式になることは容易に想像できよう。

\(A^n=\begin{pmatrix}A_1 A_2 … A_n & 0 & 0 & \cdots &0\\ 0 & A_2 A_3 … A_1 & 0 & \cdots &0 \\ \vdots & \vdots & \vdots  &  \ddots \\ 0 & 0 & 0 & 0 &A_n A_1 … A_n-1\end{pmatrix}\)

つまり対角線上に、それぞれのグループの要素数ごとの大きさの正方行列があるような行列が作れる。

それぞれの正方行列はまた、\(i\)のグループから\(i+1,i+2,…,N,1,…,i-1\)のグループの相対評価を考慮した上での、それぞれのグループ内での相対評価値を求めたものとして解釈できる。

それぞれの正方行列は一般に、全ての成分が正であるから、既約かつ非周期的である。よって、それぞれの正方行列でべき乗法を実行すれば、それぞれのグループ内での絶対的評価を計算することもできる。

つまり、対角線上のそれぞれの対角行列\(P_i=A_i A_{i+1}\cdots A_{i-1}\)と\(P_i\)の、最大固有値以外の固有値に属する固有ベクトルを除く、次元\(|G_i|\)の任意のベクトル\(x\)に対して、

\(\displaystyle \lim _{n \rightarrow \infty} (P_i)^n x/|(P_i)^n x|\)

を考えると、\(n \rightarrow \infty\)での収束が保証され、収束先のベクトルは、ページランクと同じようにそれぞれのノード固有の評価値として解釈できる。

(これは、\((A^n)^m\)を計算し、再度対角線上の正方行列を取り出すことと同じである。それぞれのグループの中でのみランク付けしたい場合、Aそのものの掛け算を計算する必要はないから、計算コストは削減できることに気付いてほしい。)

さらに、長さ\(n\)の巡回置換を複数考え、これを”道”のようなものと捉えることもできる。

道を複数考え、つまり大きさ\(n\)の巡回置換を複数考え、その巡回置換に対応するようなグループ間の評価値を行列\(A\)に足しこんだとしても、道に交わり(合流する部分、つまり2つの巡回置換\(\sigma_1,\sigma_2\)において、\(\sigma_1(i)=\sigma_2(i)\)となるような\(i\)が存在するとき)がない場合、行列の和の積は積の和に直せるから、計算コストを削減しつつ、有用な道をたくさん考えたうえで\(A\)に足しこむことも可能である。

交わりのある場合でも、繰り返しループの中でその交わりを考慮しなければならないとき以外で行列の和の積は積の和の形で書けるから、計算量は大幅に削減できる。

量子アニーリングの理論を学んでみる-その1最適化問題とイジング模型(数式)

「量子アニーリングの基礎 西森秀稔 大関真之」という本を読み始めたので、学んだことのアウトプットも兼ねつつ、記事を書いてみることにしました。量子力学についてある程度知っている方を対象にしていますが、数学的素養があれば理解できると思います。

本記事「量子アニーリングの理論を学んでみる-その1最適化問題とイジング模型(数式)」では、最適化問題をイジング模型に帰着できる、つまり、なぜ最適化問題を量子アニーリングの方法で探索できるのかを示します。

 

量子アニーリングとは、端的に言えば、\(\pm1\)の値をとるスピンを利用して、最適化問題を解こうと言う技術のことです。

電子の、スピンという、物質の磁性に関係する物理量は、量子力学的に\(\pm 1\)の値だけとりうることが知られています。そして、系(取り扱う物質全体のなすものを一塊と見たもの)のエネルギーに関係する、ハミルトニアン\(H\)という関数が最小になる固有状態(これを基底状態と言います。)を目指すように系を操作します。ハミルトニアンを、解きたい最小化問題の関数と同じになるように最初から設定しておけば、これは最小化問題が解けたことになります。

この”系”をD-waveという会社が実用化し、既に市販されているのはニュースでご存知かも知れません。

その理論的背景を、数式を用いながら説明していきます。

 

イジング模型

イジング模型とは、磁場の作用する環境下における、複数の電子のなす系を抽象化したモデルです。

おのおのの電子のスピン\(\sigma_i\)は\(\pm1\)の値をとり得ますが、これを\(\pm1\)の値をとり得るイジング変数\(\sigma_i=(\pm1)\)と考えます。他にもイジング変数を実現する方法はありますが、この記事では電子の系を考えます。

電子スピン相互作用を考えます。\(i,j\)間の相互作用の比例定数を\(J_{ij}\)とすると、ハミルトニアンへの寄与として、

\[-J_{ij}\sigma_i \sigma_j\]

と書けます。これは、スピンがそろうと系のエネルギーへの寄与は小さくなり、逆を向けば大きくなることを表しています。(\(J_{ij}>0\)の場合、これは強磁性的相互作用と呼ばれます。負の場合は今言ったことが逆になります。)

また、おのおのの電子における局所磁場に関するハミルトニアンへの寄与も考えます。局所磁場とは、外部磁場と系(自分と周りの電子たちの作るスピン)が作る磁場を合わせた磁場の、そのおのおのの電子のいる地点での磁場です。比例定数を\(h_i\)として、

\[-h_i \sigma_i\]

と書けます。これも、磁場にスピンの向きがそろえばエネルギーへの寄与は小さくなり、逆を向けば大きくなることを表しています。

よって、系全体のハミルトニアンは、

\[H(\mathbb{\sigma})=-\sum_{i<j}J_{ij}\sigma_i \sigma_j-\sum_{i=1}^N h_i \sigma_i\]

と書けます。

第一項の\(i<j\)は、同じ組を2回数えないためのものです。

 

フラストレーション

単純な系として\(J_{ij}=const>0\)の場合を考えれば、系の基底状態は自明、全てのスピンが同じ値をとれば最適化問題が解けたことになります。実際の最適化問題はそんなに単純ではありません。ここで、フラストレーションという概念を導入します。

まず、局所磁場のない(\(h_i=0\))、電子が2つ(\(N=2\))の簡単な場合を考えてみましょう。相互作用のある(\(J_{ij}\neq0\)の)電子対を、ボンドと呼びます。

\[H=-J\sigma_1\sigma_2\]この場合は簡単で、\(J_{ij}>0\)の場合はスピンがそろえば、\(J_{ij}<0\)の場合はスピンが逆になれば系全体の基底状態を実現できます。局所的な基底状態が系全体の基底状態を実現する、という言い方もできます。

電子が3つ(\(N=3)\)の場合はどうでしょうか。

\[H=-J(\sigma_1\sigma_2+\sigma_2\sigma_3+\sigma_3\sigma_1)\]

\(J_{ij}=const<0\)とします。うち2つは、スピンが逆になれば局所的には基底状態を実現できます。ところが、もう1つの電子を考えると、これはどっちを向こうが、どちらかもう一方の電子と同じ向きになってしまいます。これがフラストレーションです。フラストレーションとは、全体の基底状態が実現しているとき、必ずしも局所的な基底状態が実現されているとは限らないということです。

フラストレーションが発生すると、エネルギーの縮退(異なる状態が同じエネルギーを持ちうること)が起きます。例の電子3つの系では、\(\sigma_i=\pm1(i=1,2,3)\)ですから異なる状態の数は8つありますが、全てのスピンが同じ向きを持つ場合の2つ以外は、全て基底状態になります。

 

最適化問題とイジング模型

ところで、みなさんは電子のスピンと2進数が、相性がいいのに気がついたでしょうか。つまり、\(0\)をスピン\(-1\),\(1\)をスピン\(+1\)に対応させればよいのですから、2進変数を\(q=0,1\)とすれば、\(q=\frac{\sigma+1}{2}(\sigma=2q-1)\)と考えればよいだけです。

 

最適化問題(組み合わせ最小化問題)、つまり\(F(q_0,q_1,…q_N)\)の最小化(最大化の時は符号を逆転すればよい)について考えるために、テーラー展開

\[F(q_0,q_1,…q_N)=c_0+ \sum_{i=1}^N c_i^{(1)}q_i+\sum_{i,j}c_{ij}^{(2)}q_i q_j+\sum_{i,j,k} c_{ijk}^{(3)}q_i q_j q_k +\cdots\]

を、電子のスピンで表現することを考えます。上式に\(q_i=\frac{\sigma_i+1}{2}\)を代入すると、\(\sigma^2=1\)の関係から\(N+1\)次以上の項が消え、以下の式に変換できることがわかります。

\[F(\sigma_0,\sigma_1,…\sigma_N)=c^\prime_0+\sum_{i=1}^N c_i^{\prime(1)}\sigma_i+\sum_{i<j}c_{ij}^{\prime(2)}\sigma_i \sigma_j+\sum_{i<j<k} c_{ijk}^{\prime(3)}\sigma_i \sigma_j \sigma_k +\cdots+ c_{123\cdots N}^{(N)}\sigma_1 \sigma_2 \cdots \sigma_N\]

よって最適化問題は、多体相互作用を持つイジング模型の基底状態を探す問題に帰着できることが分かりました。

量子アニーリングの理論を学んでみる-その2量子アニーリングを語るための、数学的準備(数式)

に続きます。

量子アニーリングの理論を学んでみる-その2量子アニーリングを語るための数学的準備(数式)

量子アニーリングの理論を学んでみるシリーズ

量子アニーリングの理論を学んでみる-その1最適化問題とイジング模型(数式)

の続きです。

前回、最適化問題を量子アニーリングの手法で解ける、つまり最適化問題はイジング模型の基底状態を求めることに帰着できることを学びました。

次の記事では、量子アニーリングの手法でどのように計算するか、具体的にどのようにして基底状態を求めるのかを説明します。

まずこの記事で、そのために必要な程度に、量子力学を簡単に説明します。量子力学をある程度知っている人は、あまり読む必要はないかもしれません。下の()は興味のある方だけ読んでください。

 

(個人的に、論理を誤魔化して議論を進めるのが好きではないので、一応これを書いておきます。量子力学では、ヒルベルト空間\(\mathcal{H}\)における抽象的なベクトルの元として、系の状態を表現します。物理量は、ヒルベルト空間\(\mathcal{H}\)の元に対する作用素として定義されます。ただ、この記事では元が非可算無限個ある状態空間(例えば、粒子の運動)を考える必要がないので(可算無限でさえも)、この抽象的なベクトルは馴染み深い、普通のベクトル\(\mathbb{C}^N\)に対応付けて考えることができます。この場合、作用素は行列で表現できます。とくにスピンに関しては、状態空間の次元がたったの\(2\)なので、\(\mathbb{C}^2\)ととることができます。他の量子コンピュータの教科書を見てみても、量子力学の枠組みに関してはこの記事の程度の説明です。詳細は、「現代の量子力学 J・Jサクライ」などが詳しいです。以下は、このような枠組みのもとでの状態の一つの表現形式であり、多くの量子力学の教科書で採用されています。)

 

これから量子力学を語るために、まずはそのための言葉を準備しましょう。

\(z\)方向のスピンの\(\pm1\)に対応するベクトルを、

\[\mid \uparrow \rangle=\begin{pmatrix} 1\\0\end{pmatrix},\mid \downarrow \rangle =\begin{pmatrix}0\\1\end{pmatrix}\]

とおきます。

次に、パウリ行列と呼ばれる、スピンを表現するための行列を導入します。

\[\sigma_z=\begin{pmatrix}1 & 0 \\ 0 &-1\end{pmatrix}\]

これを、上で導入したベクトルたちに作用させると、

\[\sigma_z \mid \uparrow \rangle=\begin{pmatrix}1 & 0 \\ 0 &-1\end{pmatrix} \begin{pmatrix} 1\\0\end{pmatrix}=\begin{pmatrix} 1\\0\end{pmatrix}\]

\[\sigma_z \mid \downarrow \rangle =\begin{pmatrix}1 & 0 \\ 0 &-1\end{pmatrix} \begin{pmatrix}0\\1\end{pmatrix}=-\begin{pmatrix}0\\1\end{pmatrix}\]

となりますから、これらは固有値\(\pm1\)に対応する固有ベクトルになっていることがわかります。この固有値こそ、スピンの実際にとりうる値を表しています。

これら固有ベクトルの、規格化された(\(\alpha^2+\beta^2=1\)とした)線形結合

\[\mid \psi \rangle =\alpha \mid \uparrow \rangle+\beta \mid \downarrow \rangle \]

は、状態の重ね合わせを表現します。

規格化とは、全体の確率を\(1\)にするためのものです。今まで例に挙げてきた状態ベクトルも、全て規格化されています。実は、このベクトル\(\mid \psi \rangle\)に現れた\(\alpha^2\)は上向きスピン、\(\beta^2\)は下向きスピンを観測する確率になっています。それが量子力学の枠組みです。規格化とは、これらの和が\(1\)になるように設定することです。(正確には、一般に係数\(\alpha\)は複素数のため、係数の絶対値の二乗\(\mid \alpha \mid^2\)です。)

また、\(x\)方向のスピンを表す行列

\[\sigma_x=\begin{pmatrix}0 & 1 \\ 1 &0\end{pmatrix}  \]

を導入します。

この行列の\(\pm1\)に対応する固有ベクトルは、
\[\mid + \rangle=\begin{pmatrix}1\\1\end{pmatrix},\mid – \rangle=\begin{pmatrix}1\\-1\end{pmatrix}\]

です。確認してみてください。

\(\sigma_y\)も、もちろんありますがここでは使わないので省略します。

 

また量子力学の枠組みにおいて最も重要な仮定、時間に依存する状態ベクトルのシュレディンガー方程式は、

\[i\frac{\partial}{\partial t}\mid \psi(t)\rangle=\hat{H}(t)\mid\psi(t)\rangle\]

です。

時間に依存しないシュレディンガー方程式は、

\[\hat{H}\mid\psi\rangle=E\mid\psi\rangle\]

です。

詳しい解説は量子力学の専門書に譲りますが、これらの式にそこまで深く立ち入る必要もないので、これで十分でしょう。

量子アニーリングの理論を学んでみる-その3実際に、どのように量子計算を実行するのか(数式)

に続きます。

 

量子アニーリングの理論を学んでみる-その3実際に、どのように量子計算を実行するのか(数式)

前々回では、最適化問題を量子アニーリングの手法で解ける、つまり最適化問題はイジング模型の基底状態を求めることに帰着できることを学びました。

前回では、この記事の理解に必要な、多少の量子力学を解説しました。

この記事では、実際に、どのようにして量子アニーリングの計算を実行するのか、つまりイジング模型の基底状態をいかにしてもとめるかを解説します。

量子アニーリングの理論を学んでみる-その2量子アニーリングを語るための数学的準備(数式)

の続きです。

 

目標は、イジング模型

\[\hat{H}_0=-\sum_{i<j}\hat{\sigma}_i^z \hat{\sigma}_j^z\]

の基底状態を求めることでした。ちなみに量子力学では、表現したいものがただの数ではない(演算子や行列など)であることを強調するため、上に\(\hat{}\)をつける習慣があります。

そこで、時間とともに変化するハミルトニアン

\[\hat{H}(t)=-\sum_{i<j}\hat{\sigma}_i^z \hat{\sigma}_j^z-\Gamma(t)\sum_i^N \hat{\sigma}_i^x=\hat{H}_0-\Gamma(t)\sum_i^N \hat{\sigma}_i^x\]

を導入します。

添え字\(i\)は、それぞれの電子(サイト)のみにしか作用しないことを表します。もとの\(z\)軸方向のみのハミルトニアンに加え、電子に\(x\)方向に磁場がかかっている状況と同じですから、横磁場イジング模型のハミルトニアンと呼ばれます。

初期状態として、\(\Gamma(0)\)が第一項に比べて十分大きい場合を選ぶと、ハミルトニアンは

\[\hat{H}(0)=-\Gamma(0)\sum_i^N \hat{\sigma}_i^x\]

となり、この基底状態は自明、それぞれのサイトの基底状態の積(複数のサイトの状態を同時に表しすためのもの。\(\otimes\)を用いる。)であるから、\(\mid +\rangle=\frac{\mid \uparrow\rangle+\mid\downarrow\rangle}{\sqrt{2}}\)を用いて、

\[\mid \psi (0)\rangle=\mid++\cdots+\rangle=\mid+\rangle_1 \otimes \mid+\rangle_2 \otimes \cdots \otimes \mid+\rangle_N=\prod_{i=1}^N \mid + \rangle_i=\frac{1}{2^{\frac{N}{2}}}\prod_{i=1}^N(\mid \uparrow\rangle_i+\mid\downarrow\rangle_i)\]

となります。

これは、

\[\mid \psi (0)\rangle=\frac{1}{2^{\frac{N}{2}}}(\mid\uparrow\uparrow\cdots\uparrow\rangle+\cdots\mid\downarrow\downarrow\cdots\downarrow\rangle)\]

とも表せるように、サイト\(1\)から\(N\)までがそれぞれ上スピンか下スピンであるかの、全ての組み合わせ\(2^N\)個を等しく足し合わせたものです。もとめたいイジング模型のハミルトニアン\(\hat{H}_0\)を求めたいのだから、同じ重みの状態から出発するのは理に適います。

あとは、\(\hat{H}(t)=-\sum_{i<j}\hat{\sigma}_i^z \hat{\sigma}_j^z-\Gamma(t)\sum_i^N \hat{\sigma}_i^x\)における\(\Gamma(t)\)を徐々に\(0\)に近づければ、系は(状態ベクトルの時間に依存する)シュレーディンガー方程式

\[i\frac{\partial}{\partial t}\mid \psi(t)\rangle=\hat{H}(t)\mid\psi(t)\rangle\]

に従い、自然に、求めたいハミルトニアン\(\hat{H}_0\)の基底状態へと変化します。

各時刻において、\(\Gamma(t)\)の変化が十分小さい(これを、量子断熱と言います)ならば、各時刻における状態は、時間に依存しないシュレディンガー方程式

\[\hat{H}(t)\mid\psi(t)\rangle=E(t)\mid\psi(t)\rangle\]

と近似できる(ほぼ等しいと見なせる)から、基底状態の探索において、全ての時刻におけるハミルトニアンの基底状態をなめらかに移りながら変化していきます。

ちなみに、この近似ができる条件を、量子断熱条件と言います。

この条件の詳しい導出は、「量子アニーリングの基礎 西森秀稔 大関真之」を参照して下さい。

nemの楕円曲線暗号を基礎から学んでみる(数式)-その1 数学的準備

ブロックチェーンnemにも使われていると言うことで、僕も楕円曲線暗号を一から勉強してみることにしました。

ちなみに、nemのとか書いていますが楕円曲線暗号そのものは広く使われていて、基本的な考え方は全て同じです。

こんなpdfがあったので、大いに参考にしています。

https://written.4403.biz/source/ecc_rev30r.pdf

前提知識として、集合論とかはある程度知っておいた方が理解しやすいと思います。

 

まず、素数\(p\)について、集合

\[\mathbb{F}_p=\{ 0,1,..p-1\} \]

を考えます。

このちょっと変な集合の上では、要素\(n\)は、普通の整数の\( \{…,n-2p,n-p,n,n+p,n+2p,…\}\)と同じと考えます。

ここでは例として、\(p=5\)とします。

\[\mathbb{F}_5=\{ 0,1,2,3,4\} \]

つまり、この集合での足し算は\(3+4=2(=7-5)\)、

引き算は\(1-4=2(=-3+5)\)

掛け算は\(2\times 3=1(=6-5)\)

割り算は\(3\div 4=2(=(3+5)\div 4)\)

のようになります。

自分でいろいろ適当な数で計算して、遊んでみてもいいかもしれません。ここで不思議なのは、なぜ\(p\)を最初から素数にしたのかです。実際、\(p=9\)などの非素数にしてみると、割り算がうまくいきません。

\[1\div 3=(1+9)\div 3=(1+9\times 2)\div 3=… \]

ここまで来ればお気づきでしょう。9が3で割り切れてしまうからです。なぜ\(p\)が素数の時はすべてうまくいくのかは、ちょっと考えてみてください。

数学では、可算減算積算除算のできる集合を体(有理数体、実数体など)と言いますが、特に、ここで挙げた\(\mathbb{F}_p\)を、素体と言います。

 

次に、楕円曲線の定義をします。

\(y^2=x^3+ax+b (a,b\in \mathbb{F}_p,4a^3+27b^2\neq0)\)

なる曲線を、\(\mathbb{F}_p\)上の楕円曲線と言います。

\(4a^3+27b^2\neq0\)は、異なる\(x\)で同じ\(y\)になる(解が重複する)のを防ぐための条件で、二次方程式の場合の\(b^2-4ac\neq 0\)と同じものです。

例えば、\(a=2,b=3\)とした\(y^2=x^3+2x+3\)は、\(\mathbb{F}_5\)上の楕円曲線のひとつです。

そして楕円曲線上で、\((x,y)\)の両方が\(\mathbb{F}_p\)に属する点を、\(\mathbb{F}_p\)-有理点と言います。さらに無限遠点と言う特別な点\(\mathcal{O}\)も、この\(\mathbb{F}_p\)-有理点の仲間に入れてやります。

では、素体\(\mathbb{F}_5=\{0,1,2,3,4\}\)上の楕円曲線\(y^2=x^3+2x+3\)の有理点を列挙してみましょう。

\((1,1),(1,4),(2,0),(3,1),(3,4),(4,0),\mathcal{O}\)

確認してみてください。

(この場合は、\(x=0\)のときは有理点がありません。また、\(x\)座標が同じで\(y\)座標が\(\mathbb{F}_5\)の和に関する逆元同士(つまり足して5)になっている有理点の組がいくつかあります。この理由もちょっと考えてみて下さい。)

次に、この\(\mathbb{F}_p\)-有理点の集合での加法を考えてみます。以下は、このようにするとうまくいくよ~という計算法です。

1.2点\(P,Q\)を通る直線を引く

2.楕円曲線と直線の交点\(R^\prime \)とする。

3.\(R^\prime\)の\(x\)軸に対する対称点(ようは、\(y\)座標をひっくり返したモノ)\(R\)を和\(P+Q\)とする。

 

 

ホントにうまくいくの?と思いますが、例えば、

1.\((1,1),(2,0)\)を通る直線は、\(y=4x+2\)です。

2.\(y=4x+2\)と楕円曲線\(y^2=x^3+2x+3\)の交点は、\((2,0)\)です。

3.対称点は、\((2,0)\)です。

なんと、ちゃんと和の結果が元の仲間に入っている(有理点の集合は和について閉じて)います。

 

一般に、楕円曲線の有理点の集合は加法群(加法、減法について閉じている)になっています。

この方法に、ちゃんとした定義を与えます。

和\(R=P+Q\)に関して、

\(P\ or\ Q=\mathcal{O}\)のとき、それぞれ\(R=P,Q\)とする。

\(y_P=-y_Q\)のとき、\(R=\mathcal{O}\)とする。(\(P\)は\(Q\)の逆元\(-P\))

\(y_P\neq -y_Q\)のとき、\(x_R=\lambda^2-x_P-x_Q,y_R=\lambda(x_P-x_R)-y_P\)

ここで、\(\lambda\)は2点を通る直線(\(P=Q\)のときは接線)の傾きを表す。

\(\lambda= \frac{y_P-y_Q}{x_P-x_Q}\ or\ \frac{3x_P^2+a}{2y_P}\)

適当に素体\(\mathbb{F}_p\)上の楕円曲線上の点を選びます。これをベースポイントと呼びます。適当な整数\(d\)に対し、

\[d\times P=P+P+…+P\]

を求めることを、スカラー倍算と呼びます。

では、素体\(\mathbb{F}_5=\{0,1,2,3,4\}\)上の楕円曲線\(y^2=x^3+2x+3\)の有理点の集合\(\{ (1,1),(1,4),(2,0),(3,1),(3,4),(4,0),\mathcal{O}\} \)において、\((1,1)\)のスカラー倍算を求めてみましょう。

\((1,1)\times 2=(1,1)+(1,1)\)

\(\lambda=\frac{3\cdot 1^2+2}{2\cdot 1}=0\)

よって、\((3,4)\)

 

次に、

\((1,1)\times 3=(3,4)+(1,1)\)

\(1+4=0\)でるから、\(\mathcal{O}\)

このような、楕円曲線上の有理点のなす加法群をMordel-Weil群と呼び、加法群の要素の個数を、群位数と呼びます。この場合は、群位数は7です。

また、楕円曲線上のスカラー倍算はじゅんぐり回っていって、かならず途中で計算結果が無限遠点\(\mathcal{O}\)に等しくなります。この、何倍して無限遠点に等しくなったかを、ベースポイントの点位数と呼ばれます。

例えばこの場合は、ベースポイント\((1,1)\)に対する点位数は、\(3\)です。

また、じゅんぐりになる群を、巡回群と呼びます。よって、楕円曲線上の有理点の集合からじゅんぐりになるものだけ抜き出し集合は、巡回群です。

この場合は、\(\{ \mathcal{O} ,(1,1),(3,4)\} \)は巡回群をなします。

さて、ここらあたりでブロックチェーンNEMのTecnical Reference(リンクは日本語版)をチラッと読んでみると、だいぶん理解が進んだことに気付くのではないでしょうか。

ちなみに、NEMのこの数字がなぜ分数なのか

\[\frac{121655}{121656}\]

も、除算の説明のところを復習すればわかるでしょう。

nemの楕円曲線暗号を基礎から学んでみる(数式)-その2 スカラー倍算の高速化

に続きます。

nemの楕円曲線暗号を基礎から学んでみる(数式)-その2 スカラー倍算の高速化

nemの楕円曲線暗号を基礎から学んでみる(数式)-その1 数学的準備続きです。

次は、スカラー倍算の高速化を考えます。

楕円曲線暗号を処理するプログラムでは、暗号を生成するにせよ、暗号を解読するにせよ、スカラー倍算の計算が重要です。そこで、方法を考えてみます。

素体\(\mathbb{F}_p\)上の四則演算では、加算・減算に比べ、乗算、さらには除算はもっと手間がかかります。自分で、適当な数字でやってみれば実感できます。しかし、楕円曲線の加法公式は除算の計算を必要とするため、その計算も長い時間を必要とします。そこで、この除算を避けるため、射影座標と言うものを導入します。

射影座標\((X,Y,Z)\)において、射影座標全体の集合を\(V\)とするとき、\((x,y,z),(X^\prime,Y^\prime,Z^\prime) \in V\)に関して、

\[(X^\prime,Y^\prime,Z^\prime)=(cx,cy,cz)\]

となる\(c\in \mathbb{F}_p\backslash \{0\}\)が存在するとき、この2点を同一視(同じと考える)します。つまり、

\[(X,Y,Z)=(2X,2Y,2Z)=\cdots =(X/Z,Y/Z,1)\]

となります。そこで、普通の座標(アフィン座標)上のある点\((x,y)\)と、射影座標における座標値\((x,y,1)\)を同一視します。このとき、射影座標における点\((X,Y,Z)\)は、アフィン座標上の\((X/Z,Y/Z)\)と同一視できます。

 

射影座標における楕円曲線の定義方程式は、\(y^2=x^3+ax+b\)に\(x=X/Z,y=Y/Z\)を代入して、

\[Y^2Z=X^3+aXZ^2+bZ^3\]

となります。

今までの様に、\(a=2,b=3\)として立体グラフを描いてみました。

目をボヤッと気を抜いた感じにして、重ねて見るヤツです。

\(z\)軸に平行な平面で輪切りにすると、前回のその1のグラフになることが想像できるでしょうか。

また射影座標では、無限遠点は\(X\)座標、\(Z\)座標がともに\(0\)の点と定義します。

射影座標における加法は、以下の様になります。

2点\(P(X_P,Y_P,Z_P),Q(X_Q,Y_Q,Z_Q)\)の和\(R=P+Q=(X_R,Y_R,Z_R)\)を以下の様に定義する。

\(P\neq Q\)のとき、

\(X_R=vA\)

\(Y_R=u(v^2X_P Z_Q-A)-v^3 Y_P Z_Q\)

\(Z_R=v^3 Z_P Z_Q\)

ただし、\(u=Y_Q Z_P-Y_P Z_Q ,v=X_Q Z_P-X_P Z_Q,A=u^2 Z_P Z_Q-v^3-2v^2 X_P Z_Q\)

\(P=Q\)のとき、

\(X_R=2hs\)

\(Y_R=w(4B-h)-8Y_P^2s ^2\)

\(Z_R=8s^3\)

ただし、\(w=aZ_P^2+3X_P^2,s=Y_P Z_P,B=sX_P Y_P,h=w^2-8B\)

では、この方法で\((1,1)+(1,1)\)を計算してみましょう。

射影座標では、\((1,1,1)+(1,1,1)\)

\(w=2\cdot 1^2+3\cdot 1^2=0\)

\(s=1\cdot 1=1\)

\(B=1\cdot 1\cdot 1=1\)

\(h=0^2-8\cdot 1=2\)

よって

\(X_R=4\)

\(Y_R=2\)

\(Z_R=3\)

ゆえに、\((x,y)=(\frac{4}{3},\frac{2}{3})=(3,4)\)

見事に、前回の計算と一致しました。

 

では、等しいはずの\((2,2,2)+(1,1,1)\)も、確認のために計算してみましょう。これは同じ点であることに注意して、

\(w=2\cdot 2^2+3\cdot 2^2=0\)

\(s=2\cdot 2=4\)

\(B=4\cdot 2\cdot 2=1\)

\(h=0-8\cdot 1=2\)

よって、

\(X_R=1\)

\(Y_R=3\)

\(Z_R=2\)

ゆえに、\((x,y)=(\frac{1}{2},\frac{3}{2})=(3,4)\)

やはり一致しました。

 

射影座標を用いてスカラー倍算を計算するには、アフィン座標\((x,y)\)を射影座標\((x,y,1)\)として計算し、最後にその結果\((X,Y,Z)\)をまたアフィン座標に戻せばよいので、除算は最後のみ必要になります。こうすることで、計算コストを削減できます。

さらに計算コストを削減するアルゴリズムもありますが、それらはここでは言及しません。

nemの楕円曲線暗号を基礎から学んでみる(数式)-その3 楕円曲線離散対数問題

に続きます。

nemの楕円曲線暗号を基礎から学んでみる(数式)-その3 楕円曲線離散対数問題

nemの楕円曲線暗号を基礎から学んでみるシリーズ

nemの楕円曲線暗号を基礎から学んでみる(数式)-その2 スカラー倍算の高速化の続きです。

次はいよいよ、楕円曲線暗号の解読が実質的に、ほとんど不可能な理由を考えていきましょう。

前回まで、素体\(\mathbb{F}_p\)上の楕円曲線の有理点の集合で、適当にベースポイント\(P\)を選びスカラー倍算を計算するのは、比較的簡単であることを学びました。

計算の分野では問題の難しさは、その問題が入力に対し、どれぐらいの時間で解けるかを目安に考えることが多いです。特に、ある入力の長さ\(n\)に対して平均でどれぐらい時間がかかるかとか、最悪どれぐらいかかるかなんかを考えたりします。ここでいう簡単とか難しいと言うのは、そういう意味です。

つまり、実質的に不可能とは、現在用いられているコンピュータの性能を以て問題に取り組んでも、入力の長さが大きい場合などは、現実的な時間(せいぜい数年)では解けないことを意味します。

 

ここでは今までとは逆に、点\(Q\)が点\(P\)の何倍になっているのかを求める問題を考えます。これは、楕円曲線離散対数問題と呼ばれます。

この問題に対する解法ですぐに思いつくのは、総当たり法です。

しかし、その1の\(p=5\)の場合など、点位数が小さい場合はいいですが、NEMのTechnical Referenceを見てみれば、

「曲線上の巡回群\(G\)に対応するベースポイントを\(B\)とします。この群は、\(q=2^{252}-27742317777372353535851937790883648493\)個の元を持ちます。」

とあります。CPUのクロック周波数は\(5GHz\)程度ですから、一秒に\(5\times10^{9}\)回命令を回せます。本当は命令と言うと語弊がありますが、ここでは無視しましょう。シングルコアを想定して、楕円曲線離散対数問題の正解が、\(d=2^{200}\)のところにあった場合を考えます。一年は\(3\times 10^7\)程度ですから、

\[{2^{200}}/{5\times 10^9}/{(3\times 10^7)}\simeq 10^{43}年\]

かかってしまいます。ベースポイントの位数を\(l\)とすれば、かかる時間は最悪の場合、最悪計算量\(O(l)\)(\(l\)に比例する程度の時間)になります。

(逆に今まで考えてきた掛け算の場合は、\(100\times P=(64+32+4)\times P\)のような計算ができますから、比較的簡単です。つまり、もとの問題の計算時間の\(log\)程度で抑えることができます。)

ちょっと現実的ではありません。もうちょっと簡単に暗号を打ち破る方法はないのでしょうか。

 

Baby-step Giant-step法

楕円曲線離散対数問題\(Q=d\times P\)に対し、ベースポイント\(P\)の位数\(l\)を用いて、\(m=[\sqrt{l}]\) ([]はガウス記号)とします。このとき、\(d=sm+t\)(\(s\)は\(d\)を\(m\)で割った商\(,t\)は\(d\)を\(m\)で割ったあまり)とします。

\(R=m\times P\)として、

\[Q=d\times P=(sm+t)\times P=s\times R+t\times P\]

となるから、\(s,t\)は以下を満たします。

\[Q-t\times P=s\times R\]

そこで、

\(\{Q-t\times P|t\in \mathbb{N},0\leq t\leq m-1\}\),

\(\{s\times R|s\in \mathbb{N},0\leq s\leq m-1\}\)

のデータベースをあらかじめ計算し、この中で互いに同じものがあれば、それは答えになります。

 

しかしこれでも、計算量は\(O(\sqrt{l})\)です。

ちなみに、\(y=log(x),y=x,y=\sqrt{x}\)のグラフはこんな感じです。

 

\(l\)がさらに大きいところだと、コイツらの差はアホみたいに大きくなります。それぐらい、計算量のオーダーは大事だと言うことです。

 

さらに、\(\rho\)法と言う、よりすぐれたアルゴリズムがあります。最近の、大きなサイズの離散対数問題を突破した手法に関しては、全てこの方法が用いられています。

楕円曲線離散対数問題\(Q=d\times P\)に対し、何らかの方法(実際にはランダム)で

\[s\times P+t\times Q=s^\prime \times P+t^\prime \times Q\]

となる、互いに異なる\(s,s^\prime,t,t^\prime\)が見つかったとします。

すると、

\[(t-t^\prime)\times Q=(s^\prime-s)\times P\]

よって、\(\frac{s^\prime-s}{t-t^\prime }\)(ベースポイント\(P\)の点位数\(l\)の\(mod(l)\)での除算)が答えになります。

\(i\)番目にランダムに選ばれた点を\(R_i=s_i\times P+t_i\times Q\)とすると、このような点同士をいかに早く見つけるかが問題です。この特別な点どうしを、コリジョンペアと言います。

ランダムな点を効率的に求めるためには、つぎの\(R\)として、以下の\(f(R)\)を用います。

\[f(R)=R+M_i(ここで、i=x(R)mod4)\]

ただし、\(M_0,M_1,M_2,M_3\)はあらかじめランダムにとった点であり、\(x(R)\)は\(R\)の\(x\)座標です。

 

この\(4\)という数字は、大規模な問題の場合はランダム性の向上のため\(20\)などの大きな数にします。

これは、線形合同法と言う基本的な乱数生成法で、C言語のrandは線形合同法を用いています。ちなみに、pythonはメルセンヌツイスタです。

 

この\(rho\)法を用いた攻撃は、誕生日攻撃と呼ばれます。

点の集合から2つの点をランダムに選んで試すのと、

クラスから2人同じ誕生日の人を選んでくるのが、同じ構造を持っているからです。

詳細は省きますが、誕生日攻撃において、\(50\%\)の確率で攻撃が成功する試行回数は\(O(\sqrt{N})\)であることが、それほど難しくない数理的解析から分かります。これはBaby-step Giant-step法と同じですが、ある方法を用いることで保存すべき情報が減らせるから、こちらの方が優れているとされています。

このように、暗号の生成(スカラー倍算)に比べ、楕円曲線離散対数問題を解くのにメチャクチャ時間がかかる(少なくとも、それよりいいアルゴリズムが見つかっていない)ため、楕円曲線暗号の安全性が保障されています。

ちなみに、量子計算機を用いた楕円曲線離散対数問題の効率的な解法はすでに知られているようで、これが最近ニュースでよく見る暗号の問題かと、やっと納得に至りました。

nemブロックチェーンに使われているEigenTrustの原論文の要約・和訳(数式)

nemブロックチェーンに使われている、EigenTrustというピアの評価、管理アルゴリズムについての原論文を要約・和訳しました。

 

原論文:https://nlp.stanford.edu/pubs/eigentrust.pdf

要約・和訳版:http://enomotokan.com/wp-content/uploads/2019/01/Eigrentrust.pdf

 

この記事では、要約・和訳版をさらにまとめ、論文の概略を示します。

 

1.導入

悪さをするピアがいるため、そいつらを排除したい。そのために、それぞれのピアに固有のグローバルな、全てのピアの経験を反映させた信頼指数を提案する。

 

本論文では、”ローカルな”ピアの信頼度と”グローバルな”ピアの信頼度を明確に区別しています。

この論文でいうローカルとは、一組のピア同士における、あるピアからあるピアに対する相対的な評価のことを指します。グローバルとは、全てのピアに割り当てられた、お互いの大きさの比較によってそのピアの信頼性を絶対的に評価できる指標という意味です。

 

2.設計思想

1.システムは中央権力ではなく、自己防衛する仕組みにしなければならない

論文では、一旦4.4「基本的なEigenTrust」、つまり中央集権的に信頼度を計算するアルゴリズムを説明しますが、4.6で修正を加え、分散環境で計算する方法を説明します。

 

2.システムは匿名でなければならない。

ピアの評価は、例えばピア自身のIPアドレスのような、外部的に関連付けられたものではなく、不透明な識別子に関連付けなければならない。

論文5.1「アルゴリズムの説明」では、あるピアが、他のピアの信頼度を計算するのがどのピアなのか知ることができないようになっていることがわかります。

 

3.システムは、新しい参加者に評価を割り当ててはならない。つまり、評価はいつくかのトランザクションにおける継続したいい行いによって得られるべきであり、低い評価の悪意のあるピアが彼らの不透明な識別子を変え、新しい参加者の地位を得るのに有利であってはならない。

新しい参加者に関しては、評価値が0になるようにしてあります。

 

4.システムは計算や構造、貯蓄、メッセージの複雑性において、最小限のコストでなければならない。

論文4.6「分散EigenTrust」では、ローカル信頼度がほとんど0であり、ほとんどのピア交流がないために、計算するためのメッセージを送る回数を減らせるから、アルゴリズムの実行コストが実際、低くなることが説明されています。

 

5.システムはお互いに知っていて、システムを転覆させようとする集団に対して堅牢でなければならない。

論文5「安全なEIGENTRUST」では、悪意のあるピアがいたとしてもその影響を受けないようにするため、”基本的なアルゴリズム”をさらに修正しています。

 

3.評価システム-概略

論文では、成功している評価システムの例としてオンラインオークションシステムであるeBayを例に挙げ、その評価法を最も基本的なローカル信頼度として採用します。

具体的には、ピア同士がトランザクションごとにお互いに評価します。おのおのの時刻において、ピア\(i\)が\(j\)からファイルをダウンロードするとき、正\((tr(i,j)=+1)\)または負\((tr(i,j)=-1)\)としてトランザクションを評価します。論文でいうローカル信頼度\(tr(i,j)\)とは、端的に言えば、各トランザクションにおけるピア\(i\)の\(j\)に対する信頼度評価をいいか悪いかのみで表したものです。

 

これらを過去から現在まですべて足し合わせたものを、\(s_{ij}\)とします。

\[s_{ij}=\sum{tr(i,j)}\]

 

次に、\(s_{ij}\)を正規化します。正規化とは、出力を確率的に扱いたいときによく使う手法です。

正規化されたローカル信頼度\(c_{ij}\)として、以下の様に定義します。

\[c_{ij}=\frac{max(s_{i,j},0)}{\sum_j{max(s_{ij},0})}\]

式を見て分かるように、\(0\leq c_{ij}\leq 1\)が保証されます。\(j\)に関する和も、もちろん\(1\)になります。

ここからが本論文のキモです。

\(t_{ik}\)を導入します。

\[t_{ik}=\sum_j{c_{ij}c_{jk}}\]

これは、\(i\)が\(j\)に、\(k\)の信頼度を「訊いた」ことに基づいた信頼度を表しています。ある人に関して、信用がある人に信用があると言われれば、そのある人の信用度は高く評価できるでしょう。

これは行列によっても表現できます。\(c_{ij}\)を\(C\)の\(i\)行\(j\)列目として考えます。\(t_{ik}\)も行列として捉え、\(i\)行目を抜き出し、縦ベクトルにしたものを\(t_i\)とします。\(c_i\)も同じようにとらえます。

すると、以下のように書けます。

\[t_i=C^T c_i\]

少し計算すればわかるように、\(t_{ij}\)も\(\sum_j{t_{ij}}=1\)を満たしていますから、正規化されています。

これを繰り返す、つまり

\[t_i=(C^T)^n c_i\]

とします。

この式は、\(i\)が他のピアに訊き、さらに他のピアに訊き…を繰り返しまくったすえ、目的のピアの信頼度を知ったという意味で解釈できます。

これを無限回繰り返すと、実は、\(t_i\)は\(i\)に依らず\(C^T\)の固有ベクトル(これを、\(C\)の左固有ベクトルと言い、論文ではこちらの用語を用いています。)に収束します。(論文7「実験」ではシュミレーションで実際、10回未満の繰り返しで十分に収束することが経験的に示されます。)

ちなみに、これはべき乗法と呼ばれる、計算機を用いた基本的な固有ベクトルの求め方です。この方法では、行列の最大固有値に対応する固有ベクトルが得られることが知られています。

\(t_{i}\)は数学的帰納法により、やはり\(\sum{t_{i}}=1(ベクトルの総和)\)を満たしていますから、正規化されていることがわかります。論文4.1「ローカル信頼度の正規化」には、この性質により、繰り返しループの中で再正規化を行う必要がないから、計算コストの削減が可能になっているとあります。

この\(t_i\)こそまさに、システムのグローバル信頼度です。つまり、これはもはやピア同士ではなく、あるピア固有の評価値の絶対的尺度として捉えることができます。

以上が基本的な考え方です。

別にこの考え方は、ブロックチェーンのみにとらわれません。例えば、会社の人同士で信用度を百点満点で評価してもらい、表にします。これは相対評価の表です。これを行列とみなして固有ベクトルさえ求めれば、これは絶対評価として解釈できる、この論文の前半は、そういうことを主張しています。

(基本的な考え方はGoogleの、ページランクアルゴリズムと本質的に同じです。そもそも、何もBrinとPageが思いついたわけでもなく、昔から考察されてきていた考え方のようです。この基本的な考え方の詳細はグラフ理論に関わってきますので、そちらを参照してください。)

以降の論文後半では、このアルゴリズムを分散環境で計算可能で、ピアの匿名性を保ちつつ、悪意のあるピアに対して堅牢であるように修正します。

 

原論文:https://nlp.stanford.edu/pubs/eigentrust.pdf

要約・和訳版:http://enomotokan.com/wp-content/uploads/2019/01/Eigrentrust.pdf