ゆるふわブログ

東京大学理科 I 類 2 年の学生です.プログラミング,大学の勉強,日常生活で感じたことをゆるふわに書いていけたらなと思います.技術的に拙いところがあっても温かい目で見守っていただければ幸いです.

高速フーリエ変換 (FFT)

お久しぶりです,Ysmr-Ry です.
試験範囲なので高速フーリエ変換 (FFT) について僕の理解を記しておこうと思います.
高速フーリエ変換とは離散フーリエ変換を高速に計算するためのアルゴリズムで,コンピューターで (離散) フーリエ変換をするために重要な手法です.
フーリエ変換・離散フーリエ変換については,下記のようなサイトで調べればわかると思います.たぶん.僕読んでないけど.

www.ic.is.tohoku.ac.jp

成果物

見て楽しくないとあれかなぁと思ったので一応 Visualizer を作ってみました.
矩形波 (もどき) を Fourier 変換して表示してます.Enter を押すとマイクから音を拾ってきます.
もう一度 Enter を押すとその波形で固定して Fourier 変換してくれます.

CodePen - Fourier Visualizer

また,高速フーリエ変換の理解に欠かせないバタフライダイアグラムを自動で作ってくれるのも作ったので記事と合わせてみると良いかもしれません.↑↓ で数を増減できます.

CodePen - Butterfly Diagram Maker

4 つのポイント

以下の 4 点を押さえていればおそらく FFT を理解できるであろうと思います.

  • 分配法則
  • 約分
  • 二進数
  • 再帰的 (漸化式的) 定義

上の 3 つはすぐわかると思います.4 つ目が曲者なのでそのつもりで.
この 4 点を順に確認した後,離散フーリエ変換の高速化を考えます.

分配法則

ご存知だと思いますが一応.分配法則とは,
 2 \times (1+3) = 2\times1+2\times3 = 2+6 = 8 という風に,かっこの中の各項に掛け算を分配することを言います.

約分

これも大丈夫でしょう.
 \frac{2}{8} = \frac{1}{4} のようにすることです.

二進数

コンピューターと相性のいい数の表し方です.
我々が普段使っているのは十進数で,桁が上がるごとに桁の担っている数が 10 倍になっていきます.一の位,十の位,百の位,千の位…という風に.
二進数はこれが 2 倍になっていきます.一の位,ニの位,四の位,八の位…という風です.
例えば,11 を表したいなら,1+2+8 より,1011 となります.

再帰的 (漸化式的) 定義

ここまでは「何を当たり前のことを」という感じだったかもしれませんが,これはちょっと難しいです.
これは高校で習ったであろう漸化式を思い出してもらえるとわかりやすいです.
ご存知フィボナッチ数列が良さそうですね.
フィボナッチ数列とは,1, 1 から始めて直前の 2 つの数を足し合わせて次を作っていく数列です.
 1, 1, 2, 3, 5, 8, 13, \cdots
左から数えて n 番目のフィボナッチ数列の数を  f_n と書くと,「直前の 2 つの数を足し合わせて次を作っていく」は,
 f_n = f_{n-1} + f_{n-2}\quad(n\geq3)
と表せます.これが漸化式です.
フィボナッチ数列の作り方のように,最初 ( f_1, f_2) から順に作っていけば  f_n を得られますが,今は何もない状態からぽっと  f_n を聞かれた状況を考えてみましょう.最初から数列を作っていないので何もわかりません.
その時この漸化式を見ると,「n 番目を求めるには n-1 番目と n-2 番目を求めて足せばいいのさ!」と語りかけてきます.
そう,左から n 番目を求める問題を,n-1 番目と n-2 番目を求めると言い換えられるのです.
n-1 番目は n-2 番目と n-3 番目がわかればよく,n-2 番目は…と繰り返していけば,いつか 1 番目と 2 番目がわかればわかるというところまできます.1, 2 番目は 1 と決まっていますから,問題が解けたことになります.
このように,より小さな (この場合は左側) な問題の組み合わせに言い換えていくことができるのです.
よって,問題を解くのに必要なのは,より 小さな問題に言い換えるルール (漸化式)もっとも単純な問題の答え (1, 2 番目が 1) だけです.
これを覚えておいてください.

離散フーリエ変換 (DFT)

 W_N = e^{-i\frac{2\pi}{N}} とし,これを回転子と言います.(いきなり何)
難しそうですがなんのことはなく,名前の通り 360° を N 等分した分時計回りに回ることを表します.
例えば  W_4 によって円周を一周すると図のようになります.

f:id:Ysmr_Ry:20171108232907p:plain

はい.ちなみに円の右端を 0° として時計回りに回ってます.番号の順に回ってます.
 W_N^k W_N の回転を  k 回連続で行うことを表します.
ここでは簡単のため, W_N^k を円の中心から回転で行き着いた点を結ぶ矢印で表します.
例えば, W_4^2 ならば,円の右端から半周するので左端にいます.図の 3 の位置です.
よって,矢印 ← で表します.
すると,大きさ 8 の離散フーリエ変換は与えられた  x_i\,(i=1,2,\cdots,8) に対し,次の行列計算を行うことを指します.

 \begin{align*}
\begin{pmatrix}
X_0 \\ X_1 \\ X_2 \\ X_3 \\ X_4 \\ X_5 \\ X_6 \\ X_7
\end{pmatrix} = 
\begin{pmatrix}
\mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{↘} & \mbox{↓} & \mbox{↙} & \mbox{←} & \mbox{↖} & \mbox{↑} & \mbox{↗} \\
\mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} & \mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} \\
\mbox{→} & \mbox{↙} & \mbox{↑} & \mbox{↘} & \mbox{←} & \mbox{↗} & \mbox{↓} & \mbox{↖} \\
\mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} \\
\mbox{→} & \mbox{↖} & \mbox{↓} & \mbox{↗} & \mbox{←} & \mbox{↘} & \mbox{↑} & \mbox{↙} \\
\mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} & \mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} \\
\mbox{→} & \mbox{↗} & \mbox{↑} & \mbox{↖} & \mbox{←} & \mbox{↙} & \mbox{↓} & \mbox{↘} \\
\end{pmatrix}\begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 \\ x_7
\end{pmatrix}
\end{align*}

行列の 2 行目は 8 等分を 1 個飛ばしで,k 行目は 8 等分を k-1 個飛ばしで回っていることがわかります.
これは通常  O(n^2) の計算量になるのですが,隠れた対称性を見出すことでこの計算量を  O(n\log n) まで落とすのが高速フーリエ変換です.
バブルソートクイックソートぐらい違いますね.

高速フーリエ変換 (FFT)

高速フーリエ変換を一言でいうと,「 \frac{2}{8} って約分して  \frac{1}{4} にできんじゃん!」です.たぶん.
つまり, W_8^2 で 8 等分を 2 個飛ばしで回るのと, W_4^1 で 4 等分を順に回るのとでは区別がつきません.

f:id:Ysmr_Ry:20171109002725p:plain

f:id:Ysmr_Ry:20171108232907p:plain

これです.
これによって,行を半分に分けた時,回転していく様子が "対等" となり分配法則によってくくれるのです.
サイズの小さい場合から具体的に見てみましょう.

n = 2 の場合

離散フーリエ変換は,

 \begin{align*}
\begin{pmatrix} X_0 \\ X_1 \end{pmatrix} = \begin{pmatrix}
\mbox{→} & \mbox{→} \\ \mbox{→} & \mbox{←}
\end{pmatrix}\begin{pmatrix}
  x_0 \\ x_1
\end{pmatrix}
\end{align*}

となります.
ここで,行を半分 (1 個ずつ) に分けてみましょう.
1 行目はどちらも同じです.
2 行目は右側は左側を半周 (←) 回したものです.
よって,
 \begin{align*}
\begin{pmatrix} X_0 \\ X_1 \end{pmatrix} &= \begin{pmatrix}
\mbox{→} x_0 + \mbox{→} x_1 \\
\mbox{→} x_0 + \mbox{←} x_1 \\
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→}(x_0 +x_1) \\
\mbox{→}(x_0 + \mbox{←} x_1) \\
\end{pmatrix} \\
&= \mbox{→}
\begin{pmatrix}
x_0 +x_1 \\ x_0 + \mbox{←} x_1
\end{pmatrix} \\
&=
\begin{pmatrix}
x_0 +x_1 \\ x_0 + \mbox{←} x_1
\end{pmatrix}
\end{align*}

 \mbox{→} = W_2^0 = 1 ですからこうなります.
これを図にしてみるとこうなります.

f:id:Ysmr_Ry:20171109004714p:plain

左から線をたどっていきます.合流したら足します.
三角があったら流れてきた数に傍に書いてある回転子を書きます.
この図は, X_0 = x_0 + W_2^0 x_1 = x_0 + x_1 X_1 = x_0 + W_2^1 x_1 = x_0 + \mbox{←} x_1 を表しています.

この図を,蝶のようなのでバタフライダイアグラムといいます.

n = 4 の場合

離散フーリエ変換は,

 \begin{align*}
\begin{pmatrix}
X_0 \\ X_1 \\ X_2 \\ X_3
\end{pmatrix} = \begin{pmatrix}
\mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} \\
\mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} \\
\mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} \\
\end{pmatrix} \begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3
\end{pmatrix}
\end{align*}

今度は偶数行目と奇数行目で考えてみましょう.偶数行目は行を半分に分けた時に同じものの繰り返しですが,奇数行目は半周分ずれています.
よって.

 \begin{align*}
\begin{pmatrix}
X_0 \\ X_2
\end{pmatrix} &= \begin{pmatrix}
\mbox{→} & \mbox{→} & \mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{←} & \mbox{→} & \mbox{←} \\
\end{pmatrix} \begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} x_0 + \mbox{→} x_1 + \mbox{→} x_2 + \mbox{→} x_3 \\
\mbox{→} x_0 + \mbox{←} x_1 + \mbox{→} x_2 + \mbox{←} x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} (x_0+x_2) + \mbox{→} (x_1+x_3) \\
\mbox{→} (x_0+x_2) + \mbox{←} (x_1+x_3) \\
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} & \mbox{→} \\
\mbox{→} & \mbox{←} \\
\end{pmatrix} \begin{pmatrix}
x_0+x_2 \\ x_1+x_3
\end{pmatrix} \\
\begin{pmatrix}
X_1 \\ X_3
\end{pmatrix} &= \begin{pmatrix}
\mbox{→} & \mbox{↓} & \mbox{←} & \mbox{↑} \\
\mbox{→} & \mbox{↑} & \mbox{←} & \mbox{↓} \\
\end{pmatrix} \begin{pmatrix}
x_0 \\ x_1 \\ x_2 \\ x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} x_0 + \mbox{↓} x_1 + \mbox{←} x_2 + \mbox{↑} x_3 \\
\mbox{→} x_0 + \mbox{↑} x_1 + \mbox{←} x_2 + \mbox{↓} x_3
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} (x_0+\mbox{←}x_2) + \mbox{↓} (x_1+\mbox{←}x_3) \\
\mbox{→} (x_0+\mbox{←}x_2) + \mbox{↑} (x_1+\mbox{←}x_3) \\
\end{pmatrix} \\
&= \begin{pmatrix}
\mbox{→} & \mbox{↓} \\
\mbox{→} & \mbox{↑} \\
\end{pmatrix} \begin{pmatrix}
x_0+\mbox{←}x_2 \\ x_1+\mbox{←}x_3
\end{pmatrix} \\
\end{align*}

ここまでが第 1 段階です.これは n = 2 の場合とやっていることが全く同じです.
n = 2 の場合の  x_0 +x_1 が偶数行目に  x_0 + \mbox{←} x_1 が奇数行目に対応しています.(n = 2 の場合はそれぞれ 1 個ずつだったのが,今は 2 個ずつ)
しかし,添字の番号が違っていますね. x_0 + x_1 に対応するのが  x_0+x_2, x_1+x_3 です.ここにどのような対応があるのでしょう?
答えは二進数にあります.添字を二進数で表してみると,(0,2) が (00,10), (1,3) が (01,11) です.二進数の一番上の桁に注目すると,どちらの組も (0,1) となっていることがわかります.これがルールです.
これが第 k 段階ならば上から k 桁だけに注目します.あとで詳しく考えます.
次に,
 \begin{align*}
\begin{pmatrix}
X_0 \\ X_2
\end{pmatrix} &= \begin{pmatrix}
x_0+x_2 + \mbox{→}(x_1+x_3) \\ x_0+x_2 + \mbox{←}(x_1+x_3)
\end{pmatrix} \\
\begin{pmatrix}
X_1 \\ X_3
\end{pmatrix} &= \begin{pmatrix}
x_0+\mbox{←}x_2 + \mbox{↓}(x_1+\mbox{←}x_3) \\ x_0+\mbox{←}x_2 + \mbox{↑}(x_1+\mbox{←}x_3)
\end{pmatrix}
\end{align*}

で完成です.(第 2 段階)
これをバタフライダイアグラムで表すと,

f:id:Ysmr_Ry:20171109093939p:plain

となります.左側の列の蝶の集まりが第 1 段階,右側が第 2 段階に対応しています.
第 1 段階は n = 2 の場合に対応すると言ったように,n = 2 の時のバタフライが 2 つ入っています.
サイズ n の FFT を実行することを, fft_n(x_i) と表すとすると, fft_4(x_0,x_2,x_1,x_3) = fft_2(x_0,x_2) + fft_2(x_1,x_3) + (\mbox{一番右のバタフライ}) です.これが漸化式です.

n = 8 の場合

これも同様です.バタフライダイアグラムだけ示すと,

f:id:Ysmr_Ry:20171109094547p:plain

となります.1 列目 (第 1 段階) には n = 2 の場合のバタフライが 4 つ,1, 2 列目を合わせると n = 4 の場合のバタフライが 2 つ隠れています.
この場合も, fft_8(x_0,x_4,x_2,x_6,x_1,x_5,x_3,x_7) = fft_4(x_0,x_4,x_2,x_6) + fft_4(x_1,x_5,x_3,x_7) + (\mbox{一番右のバタフライ}) が漸化式になります.

先ほど確認したように,問題を解くには漸化式の他にもっとも単純な問題の答えが必要です.
もっとも単純なのはこの場合 n = 1 の場合を考えると簡単です.n = 1 の場合は,フーリエ変換する前もなくそのままになります.
 fft_1(x_i) = x_i です.これで問題が解けます.

解けるんですが,あと問題なのは,バタフライダイアグラムの左側の  x_i の並びです.これにはある規則があります.

なぜビットリバースの順なのか

バタフライダイアグラムの左側の  x_i の右に書いてあるかっこの中の数は添字を二進数で表したものです.
これをそれぞれ逆から読んでみると,000, 001, 010, 011, 100, … で,0, 1, 2, 3, 4, … と順に増えていることがわかります.これをビットリバースの順と言います.
なぜこうなるのか考えてみましょう.漸化式に沿って考えてみると,右から 2 番目のバタフライが 2 つある段階 (第  \log_2 n-1 段階) において,添字を二進数で表した時の上から  \log_2 n-1 桁だけ見ています.言い換えると,一番右の桁だけ無視しています.なので,( \log_2 n-1 桁)0 と ( \log_2 n-1 桁)1 の 2 つのセットがあります.これが 2 つのバタフライに対応します.
漸化式から一番右のバタフライと 2 つの半分のサイズの  fft に言い換えられるため,半分のサイズのバタフライダイアグラムに今と同じことを順に適用していくと,n = 2 の場合,一番上の桁が 0, 1 なので,ビットリバースの順に増えていくことに対応します.
例えば,n = 4 の場合,第 1 段階は 0, 1 ですが,第 2 段階まで考えると,その後ろに 0 をつけるか,1 をつけるかの 2 セットあります.これを順に書くと,00, 10, 01, 11 でビットリバースの順になっています.
n = 8 ならば,00, 10, 01, 11 の後ろに 0 をつけるか 1 をつけるかで,000, 100, 010, 110, 001, 101, 011, 111 です.
ビットリバースの順に増えていくのは,そこまでの列の後ろに 0 を加えたものと 1 を加えたものをくっつけて拡張するということを繰り返していくと作れます.
ここにも再帰的定義が隠れているのです.
(伝わる?)

まとめ

個人的にネットに転がってる FFT の説明がわかりにくかったので自分なりにまとめてみました.
徐々に難しくしていったつもりなんですが,ビットリバースのところがうまく伝わったか心配です…
不明な点があったらコメントください.間違っていた場合もお願いします.

行列と複素数 (小ネタ)

とりあえずなんでもいいから書かなきゃなぁと思ったので簡単な小ネタ.


複素数って回転拡大を表すじゃないですか.
 r(\cos\theta+i\sin\theta) って  r 倍拡大, \theta 回転を表しますね*1
ちょっと書き方を変えると  a+bi = \sqrt{a^2+b^2}\left(\frac{a}{\sqrt{a^2+b^2}}+\frac{b}{\sqrt{a^2+b^2}}i\right) という複素数 \sqrt{a^2+b^2} 倍拡大, \phi = \arctan \frac{b}{a} 回転を表します.


回転で行列といえば回転行列というのもありますね.
 \theta 回転ならば  \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} です.
拡大はこの行列の前に拡大率をかければよいので, a+bi の表す回転拡大に対応する行列は  \sqrt{a^2+b^2}\begin{pmatrix} \cos\phi & -\sin\phi \\ \sin\phi & \cos\phi \end{pmatrix} = \begin{pmatrix} a & -b \\ b & a \end{pmatrix} となります.(というか基底を直交する  (a,b) (-b,a) に張り替えると考えた方が早い.)


逆にこの行列を  a+bi とみなせば見通しがよくなります.

例えば,線形代数で †随伴行列† ってやるじゃないですか.中二病みたいなマーク使うやつ.
複素行列  A に対して  A^† A を転置して各成分の複素共軛をとります.
これは実数行列の転置に当たる操作なんですけど,複素数バージョンの場合はなんで共軛をとるんだろうって不思議じゃないですか.(初めて習ったときそう思った)
これが各成分の複素数を対応する行列に書き換えるとわかるんです.
例えば,
 \begin{align*}
\begin{pmatrix}
  1 & i \\ i & 1
\end{pmatrix}
\end{align*}
は,対応する行列が  1 \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} i \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} なので,
 \begin{align*}
\begin{pmatrix}
  1 & i \\ i & 1
\end{pmatrix} &= \begin{pmatrix}
  \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} \\
  \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} & \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}
\end{pmatrix} \\
&= \begin{pmatrix}
 1 & 0 & 0 & -1 \\
 0 & 1 & 1 & 0 \\
 0 & -1 & 1 & 0 \\
 1 & 0 & 0 & 1
\end{pmatrix}
\end{align*}

こいつの転置をとると,
 \begin{align*}
\begin{pmatrix}
 1 & 0 & 0 & 1 \\
 0 & 1 & -1 & 0 \\
 0 & 1 & 1 & 0 \\
 -1 & 0 & 0 & 1
\end{pmatrix} &=  \begin{pmatrix}
  \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \\
  \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} & \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}
\end{pmatrix} \\
&= \begin{pmatrix}
1 & -i \\ -i & 1
\end{pmatrix} 
\end{align*}

 \begin{pmatrix}
  1 & i \\ i & 1
\end{pmatrix}^† になってます.複素共軛は対応する行列の転置に当たるので,複素数という行列のブロックを転置して,各々の共軛を取れば全てを展開した実行列の単なる転置に当たるのです.†納得†.

あとは指数行列ですね.指数関数に行列をぶち込むっていう.
 e^x = 1+x+\frac{x^2}{2}+\cdots から類推して,
 e^A = E+A+\frac{A^2}{2}+\cdots と定義します.

例えば  e^{\begin{pmatrix} a & 0 \\ 0 & a \end{pmatrix}} = \begin{pmatrix} e^a & 0 \\ 0 & e^a \end{pmatrix} ってなります.
で, e^{\begin{pmatrix} 0 & -\theta \\ \theta & 0 \end{pmatrix}} = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} になるんですね.定義に代入して確かめればいいんですけど,直感的に収束先を捉えたくて.そこでこの小ネタの出番です.


 1 \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} i \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} だったじゃないですか.それぞれ  E,\ I と表します.すると,左辺は  e^{I\theta},右辺は  E\cos\theta + I\sin\theta じゃないですか.どっかで見た事ありますよね?*2 そういう事です.

 e^{\begin{pmatrix} a & -b \\ b & a \end{pmatrix}} = e^a\begin{pmatrix} \cos b & -\sin b \\ \sin b & \cos b \end{pmatrix} もわかると思います.


ちなみに  \begin{pmatrix} 0 & \theta \\ \theta & 0 \end{pmatrix} みたいなのは複素数を拡張した分解型複素数なる "双曲線の" 複素数を考えれば納得できそうです.いつもの複素数は "円の" です. j^2=1 なる虚数単位チックな  j を導入して  a+bj とかするみたいです. j J=\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} に対応付ければ 2 乗と 1 ( E) なのでつじつまが合います.

 e^{j\theta}=\cosh\theta+\sinh\theta なる関係が成り立つので同じように  e^{\begin{pmatrix} 0 & \theta \\ \theta & 0 \end{pmatrix}} = e^{J\theta} = E\cosh\theta+J\sinh\theta = \begin{pmatrix} \cosh\theta & \sinh\theta \\ \sinh\theta & \cosh\theta \end{pmatrix} です.

Lorentz 変換とかに関係ありそう.


あと全然関係ないんですが,なんで共軛とるかわからないといえばあれです.複素ベクトルの内積をとるときって片方の共軛をとるじゃないですか. \vec{a}=\begin{pmatrix} a_1 \\ a_2 \\ a_3 \end{pmatrix}, \vec{b}=\begin{pmatrix} b_1 \\ b_2 \\ b_3 \end{pmatrix} なら, \vec{a}\cdot\vec{b} = \bar{a_1}b_1+\bar{a_2}b_2+\bar{a_3}b_3 みたいに.あれがなんでかって話ですね.まぁ  \bar{z}z=|z|^2 なので, \vec{a}\cdot\vec{a} = \lvert a_1\rvert^2+\lvert a_2\rvert^2+\lvert a_3\rvert^2 = \lvert\vec{a}\rvert^2 となるようになってるといえばおしまいです.1 次元複素ベクトル (複素数) 
 \alpha = \lvert\alpha\rvert e^{\arg\alpha}, \beta = \lvert\beta\rvert e^{\arg\beta}内積を考えると, \bar{\alpha}\beta = \lvert\alpha\rvert\lvert\beta\rvert e^{\arg(\beta-\alpha)} となります.これの実部をとると 大きさ x 大きさ x cos (なす角) なので内積っぽいです.

計算幾何とかで複素数クラスをベクトルに流用する場合,このようにして内積外積をスマートに記述したりします.

以上.行列を複素数に対応付ければ見通しが良い例,無数にありそうですが僕が観測した範囲ではこんな感じです.他にあったら教えてください.

もっと本格的な記事も気が向けば書きます (書くとは言ってない)

*1:知らない人は複素平面で調べればわかるかも

*2: e^{i\theta}=\cos\theta+i\sin\theta,Euler の公式で調べてみてください

ボロノイ図を描く! ~Fortune's Algorithm~

こんにちは! Ysmr-Ry です.
Twitter の知り合いにボロノイ図の論文を見せてもらって興味が湧いたのでボロノイ図の描き方を勉強してみました!
ボロノイ図とはなんなのか,そしてその衝撃の描き方を見てみましょう.
ビーチに波が押し寄せるようにボロノイ図ができてゆきます.(面白そう)
長くなってしまいましたが,なぜ放物線なのかがわかれば十分な気がします (ごめんなさい)

ボロノイ (Voronoi) 図とは

ボロノイ図とは,与えられたいくつかの母点のうちどれに一番近いかで領域を分けた図のことです.
実際にボロノイ図を見てみましょう.

f:id:Ysmr_Ry:20170906172456p:plain

これがボロノイ図です.黒い点が母点 (ボロノイ図を描くために与えられた点) で色分けされてるのがそれぞれのボロノイ領域です.
よく見てると気づくとは思いますが,隣接するボロノイ領域の境界線は対応する母点を結んだ線分の垂直二等分線です.

球面上で作ったり,いろいろ拡張はあるようです.
次のリンク先は地球上の空港を母点として球面上にボロノイ図を描いています.

World Airports Voronoi

Fortune's Algorithm

ボロノイ図の描き方は当然いろいろあるとは思いますが,今回紹介するのは Fortune's Algorithm という方法です.
Fortune's Algorithm は平面走査法の一種で,円錐曲線としての放物線に着想を得ています.
簡単言うと,ボロノイ図を描きたい平面上を横線が上から下に動いていってスキャンし,それに放物線を重ね合わせた "波" が続いて,波の継ぎ目がボロノイ領域の境界をなぞっていきます.
実際見た方が早いですね.こんな感じです.( 点の配置によってはバグるのでそういう時はうまくいくまでページを更新してください((((( )

Voronoi Visualizer

水色の線は放物線の組み合わせでできている "波" で,それよりも上の領域が "海" と考えると良いです.
黒い横線が平面上を走査して行って,緑色のボロノイ図が描かれていくのがわかります.
赤い丸はなんなのか,右の木 (こういう形のデータ構造を木と言います) は何を表しているのか,はこの後説明します.

前提知識 1 (円錐曲線としての放物線)

何と言ってもこのアルゴリズムの面白いところは放物線が生きてくるということです.
放物線,2 次関数自体は中学生でもやると思いますが,ここで使う放物線の性質は現在の課程では数 III の 2 次曲線のところで出てきます.
それは,「ある定直線  d とその上にない定点  F があるとき, d F から等距離にある点の軌跡は放物線である」です.
定直線  d を準線 (directrix),定点  F を焦点 (focus) と言います.
アニメーションで見るとこんな感じです.

紫の定点と定直線がそれぞれ焦点と準線,赤いのがそれによってできる放物線,オレンジ色の円は焦点を通り定直線に常に接しているような円です.
ということはこの円の中心は焦点と準線から同じ距離 (半径) にあるので,放物線の上を動くというわけです.
ちなみにこれを一般化して焦点と準線からの距離の比が  \varepsilon となるような点の軌跡を離心率  \varepsilon の二次曲線といい,放物線・楕円・双曲線になります.

前提知識 2 (Binary Search Tree)

先ほどの Voronoi Visualizer の右側に表示されていた木を二分探索木 (Binary Search Tree; BST) といいます.
データ構造とか言ってますがそんなに難しい話じゃないので見ていきましょう.
木というのは先ほど見ていただいた通りの形をしたデータ構造なのですが,中でも二分木はそれぞれの頂点 (丸のこと) から伸びる枝が高々二股に分かれているものを言います.
それぞれの頂点には例えば数字が入っています.二分探索木は,「ある頂点に注目した時にその中に入っている値に対して,そこから左に生えている枝の先の部分木に入っている値は全てそれよりも小さく,右は全て大きいという性質を保つもの」です.
図で見ていきましょう.次のような木は二分探索木と言えます.

f:id:Ysmr_Ry:20170906184935p:plain

一番上の 5 の入った頂点 (根と言います) に注目します.すると,左の部分木は下の四角で囲った部分です.

f:id:Ysmr_Ry:20170906185322p:plain

見てみるとこの中のどの値も 5 より小さくなっていることがわかります.(1, 2, 3 < 5)
また右の部分技についてはどの値も 5 より大きくなっていることがわかります.(5 < 6, 8, 11)
今回は根に注目しましたがこれの性質はどの頂点についても成り立っています.(例えば 2 の頂点なら 1 < 2 で 2 < 3 です)
こういう風に数値のデータを持っていて何がうれしいかというと,探索木ですから,この中からある値を探す時に効率がいいのです.
例えばこの二分探索木から 6 を探したいとしましょう.
根から順にたどっていきます.根に入っている値は 5 です.5 < 6 なので,先ほど見た性質から,6 は少なくとも右の部分木にあることがわかるので,右の枝を辿ります.
これを繰り返せば,最悪の場合でも木の深さ分だけ比較をすれば見つかることがわかります.

f:id:Ysmr_Ry:20170906190834p:plain

難しいことを言うと木の深さ  d に対してその深さにある頂点はだいたい  2^d 個あるので計算量は (平衡状態であれば)  O(\log n) となります.
Fortune's Algorithm では二分探索木を放物線の波の最前線を管理するのに使います.

Fortune's Algorithm (なぜ放物線を考えるのか)

では本題の Fortune's Algorithm について考えていきましょう.
Fortune's Algorithm は平面走査法なので,走査線 (Sweep Line) が平面上を上から下に動きます.線が掃いていくから Sweep です.
ただこのアルゴリズムの場合は特殊で,もう一つ放物線の組み合わせの "波" である汀線 (Beach Line) を考えます.
これがなぜかということです.
今回はボロノイ図を構成するためのアルゴリズムなので,走査線よりも上の掃き終わった領域においてはボロノイ図の一部が組みあがっていて欲しいのです.
しかしながら,走査線よりも上ならば下にある母点たちの影響を全く受けないとは言えません.
そこで,それよりも上ならば下の母点たちから一切影響を受けないというもう一つの境界線,汀線を考えるのです.
先ほど言ったように汀線の上側の領域を "海" と呼称すると,ボロノイ図の出来上がった領域である "海" がじわじわと平面全体を覆っていくということです.
ではなぜその汀線が放物線となるのか.それは前提知識 1 で見た放物線の性質によるものです.
汀線は走査線が掃き終わった,上にある全ての母点をそれぞれ焦点とし,現在の走査線を準線とした放物線よりも上の領域のうち,もっとも下にある境界線をつぎはぎしてできた最前線として決めます.(日本語が不自由)
例えば,先ほどの Visualizer で放物線を全て表示すると

f:id:Ysmr_Ry:20170906193730p:plain

となるのですが,実際に考えるのは一番下にある,最前線の放物線の組み合わせです.
そうするとそれらの放物線は「走査線と対応する母点から等距離にある点の軌跡」なので,放物線の上の領域は「走査線よりも母点の方が近い点の集まり」になります.走査線よりも下の点は当然走査線よりも遠いですから,"海" は走査線よりも下の点から影響を受けないことがわかります.
また,隣接する放物線の交点を Breakpoint と呼ぶのですが,その点は準線と二つの母点からの距離が等しいのでボロノイ図の辺の上に乗ります.

Fortune's Algorithm (2 つのイベント)

平面走査法はアニメーションとしては連続的に走査線が動いていくのですが,実際のプログラムでは特別なイベントが起こる時点をとびとびにたどっていきます.
今回の場合そのイベントは 2 つあります.

Circle Event

f:id:Ysmr_Ry:20170906201010p:plain

図のように 3 つの母点が同一円周上にあり,その円の下側に走査線が接する時点を言います.(したがって上の図は Circle Event が起こる少し前)
このまま走査線が下に行くと,赤い円の中心に集まってきている 3 つの放物線の一部のうち,真ん中が完全に汀線からなくなります.
また,円の中心は Circle Event のまさにその時点では 3 つの放物線に同時に乗っていると言うことができ,それに対応する 3 つの母点 (と走査線) から等距離にあるため,ボロノイ図における点 (ボロノイ点) になります.
したがってこのイベントが起きた時は真ん中の放物線を汀線の放物線のリストの中から消去し,また,ボロノイ点で終わっている辺の終わりにボロノイ点をセットし,ボロノイ点から始まる辺の始点をセットします.

Site Event (NewPoint Event)

f:id:Ysmr_Ry:20170906210749p:plain

Site Event は母点の一つに走査線が触れた時点で起きます.上の図は走査線が 3 の母点を掃いた直後です.
母点 3 を焦点とする放物線が生じるため,それを汀線に組み込まねばなりません.

Fortune's Algorithm (実装の詳細)

先ほどのイベントは上から順番に並べて順に処理します.この時 y 座標が最小のものを効率よく取り出すデータ構造が必要で,Priority Queue というのを使うみたいです.ちなみに最初は母点から Site Event だけを作って並べて置いて,それを処理ながら Circle Event を発生させてイベントリストに追加します.
このアルゴリズムで使うもう一つのデータ構造は先ほども紹介したように二分探索木です.
頂点の持つ値として何を入れているかというと,Breakpoint の x 座標です.
Visualizer で (1,2) と頂点に表示されていたら,母点 1 と母点 2 を焦点とする放物線が左からこの順で交わっている時の交点の x 座標が入っています.(実質的に)
Site Event のところの図がわかりやすいと思います.(再掲)
f:id:Ysmr_Ry:20170906210749p:plain

この場合は,(0,2) (2,3) (3,1) (1,0) の順で左から Breakpoint が並んでいるのがわかります.(このうち画面内に見えているのは (2,3) (3,1) だけですが…)
ちなみにこの順序は二分探索木を深さ優先探索 (Depth-First Search) すると得られます.

Circle Event においては真ん中の放物線を消し (つまり真ん中と左,真ん中と右の作る Breakpoint を消す),新たな Breakpoint として円の中心を追加します.
Site Event では新たに放物線ができるので,既存の汀線とどこでぶつかるかを調べて,新たに Breakpoint を追加しなければなりません.
関係する汀線の部分を探索するのに効率が良いのです.

また実装の際には母点の座標と準線の位置さえあれば放物線も Breakpoint も計算できるのでそれさえ持っておけば OK です.
辺は始点だろうが終点だろうがどのみち Circle Event の時の円の中心になるのでその時決定すれば良いです.

詳しくは Haskell で実装してる人がいるのでそちらでどうぞ.

github.com

また,他の言語での実装もあります.(英語だけど)

Fortune’s algorithm and implementation

おわり

日本語の資料少ないので苦労しました…

http://www.cs.sfu.ca/~binay/813.2011/Fortune.pdf
http://cs.nyu.edu/yap/bks/egc/09/7Vor.pdf

ちなみに僕の Visualizer のソースコードは以下にあります.Haskell の Haste というライブラリを使って JavaScript のコードを生成しています.
したがって描画は HTML5Canvas です.
僕が確認しただけでも以下のようなバグがあるので,Pull-Request をいただけると助かります…(((((

・最初の母点 0, 1 だけの放物線の場合,片方が表示されずかくかくする (治った?)
・点の配置による原因不明のフリーズ
・放物線の継ぎ目がはみだしてる (CircleEvent が発生していない) ← CircleEvent がおきその先は画面によって切り取られているときにおかしくなる

これアルゴリズムの性質上 Circle Event が起きないと辺が終わらないので画面に切り取られているような辺は自分で計算しなきゃいけないんですよね…
半直線的な.それ僕が改造したんですが,なんかボロノイ辺がありえない感じになる場合があるみたいで… (僕のせい)

github.com

だいぶ下調べはしたはずなんですけど,間が空いちゃって今は Haskell のソースちらちら見ながら書いてるだけなので,アレなところがあったらご指摘お願いします(((((

プラネタリウム作ってみた

こんにちは!Ysmr-Ry です!

まあ知っての通り「放課後のプレアデス*1を愛してやまないので星に関わることをしたくって.

この度ヒッパルコス星表のデータを使い,JavaScriptcanvasプラネタリウムというか動く天球を作ってみました!
成果物で遊ぶだけでも楽しいので見ていってください〜

成果物

これです.

https://ysmr-ry.github.io/Hipplanetarium/

f:id:Ysmr_Ry:20170806104120p:plain

操作方法は左の方に書いてある通りなんですが,一応.

キー 操作内容
A 星座強調 On / Off
S 星座線 On / Off
X 星の名前表示 On / Off
C 星座名 非表示 / 和名 / 英名
Q モード切り替え 天球 / 天球距離反映 / 宇宙空間
Z 拡大
星をクリック 星にフォーカス / 詳細を表示
Space 星座ツアー開始 / 次の星座へ
Shift + Space 前の星座へ
矢印 (←↑↓→) キー 天球を回転

ちなみに星をクリックしてフォーカスしている時に,その星がどこかの星座に属していたならばそのまま Space を押すとその次の星座から星座ツアーが始まります.
クリックしていた星が属している星座にフォーカスしたい時は,Space → Shift + Space とするとできます(
星座線表示しながらモードを SphereDist にすると,すばるの「よりそうように輝く星も、本当は、ひとつひとつが、何光年も遠く遠く離れています。」というセリフも理解できるんじゃないでしょうか.
あ,完全に PC でやることを前提としてます.スマホ勢のみなさんごめんなさい(((

ヒッパルコス星表

ヒッパルコス星表とは,1988 年に打ち上げられた位置天文衛星ヒッパルコスが 4 年間かけて収集した恒星のデータです.
ちなみにヒッパルコスとは古代ギリシア天文学者の名前なのだそう.(英語の綴り違うみたいだけど)
今回の成果物はそのヒッパルコス星表の基本的な Visualizer です.
ちなみに Hipplanetarium という名前は Hipparcos + Planetarium ということです.
(ネーミングセンスが絶望的にない)

ちょっと舞台裏

ちょっと中身の話についてすると,canvas で描画してるので,遠近法の計算は自前でやっていて,3 次元回転には四元数 (Quaternion; クォータニオン) を使ってます.
複素平面とか知ってる人はその 3 次元 Ver. と思っていただければいいと思います.
まあその程度しか知らないんですが.

あとは基本データをうまく処理して描画するだけだったんですが,ちょっと苦労したのが星の色ですね.
ヒッパルコス星表には安直に星の色が RGB でいくら,みたいなことは載ってないので,そこはあるデータから自分で計算せねばならんのです.
今回は bv Color という値から下記の式で表面温度を推定し,そこから Plank の黒体輻射の分布式を XYZ 等色関数で数値積分して XYZ 表色系での値を求め,RGB に線形変換してます.

f:id:Ysmr_Ry:20170806102538p:plain

このサイトを参考にしました.

www.geocities.jp

ただ,これで計算すると RGB が負になったりしてよくわかんなかったところはあるんですね.
ちょっとだけ色のお勉強をして見た*2んですが,どうも RGB ですべての色が表せるという単純な話ではないようで,黒体の色温度がディスプレイで再現できる色の範囲を出てしまってるってことなんだと思いますたぶん.こんな風に黒体軌跡が色域の三角形を出てる感じ?

f:id:Ysmr_Ry:20170806103148p:plain

あとは色相を保って明るさ最大にしてるので,色の違いが目立っていいのかなぁと.
星の色ちゃんとあってますかね… さそり座のアンタレスが赤いから大丈夫だろうという謎の基準()
おかしかったら星に詳しい人教えてください〜

あとの詳細が知りたい人は僕にコメントとかで聞いてもらうか,Github にソースが上がってるので解読して貰えばいいと思います.(拙くて申し訳ないですが)

github.com

これの docs ってフォルダーの中です.

ちなみに

やっぱり同じようなこと考えた人はいて,もっとやばいことやってるのがありました.

HippLiner -- A 3D interstellar spaceship simulator with constellation writing function

でも星の知識とか絶望的に足りないし今回はとりあえずこんな感じで… orz

おわり

まあ,僕大昔に同じようなの作ったことあるんですが,その時は星の色とかかなり適当だったので作り直してみようかと思いまして.
それに canvas 使えばまあ Web ブラウザが動けばどこでも見れますし.
あと夏休みで暇だったので()

*1:いい加減しつこいのですがいいものはいいのです笑
なんか昨日はいつきちゃんの誕生日だったらしくて上映会をやってたみたい (?) ですね.
よかったら見てみて〜
sbr-gx.jp

*2:このサイトがわかりやすかったです.
qiita.com
qiita.com
あと,このサイトが参考にしてる本も特に 3 章の表色系の話が数学と絡んできて面白かったです.
( というかそれしか読んでない() )

色彩工学入門 定量的な色の理解と活用 | 森北出版株式会社

相対性理論の本を読んだ

こんにちは,Ysmr-Ry です.

大学が夏休みに突入したので春休みに買って挫折していた相対論の本に再チャレンジしました.
これです.

f:id:Ysmr_Ry:20170730211259j:plain

www.amazon.co.jp

これ,発売日が 2017 年 3 月 27 日だったんですね.
大学の書籍部に電話して予約してたのを思い出します…
1 年の時周りに相対論の授業をとってる人がいて,自分も勉強してみたいな,でも授業受けるだけじゃどうせわからないだろうしな,自分で歩み寄らなきゃなぁと思っていた僕にとっては春休みにこの本が発売されるのは天啓のように感じられましたw
「知る限りこの分野で一番分かり易い本」という評判がついてるくらいですしw
ちなみに再チャレンジしようと思ったきっかけは,「放課後のプレアデス」を見たからなんですね.なんか劇中で,「ブラックホールのシュワルツシルト半径がこんなにくっきり見えるなんて」「宇宙項 Λ が…」とか言ってて,似非物理屋としては「やらねば!」と思ったんですね.
物理屋は食い尽くし,もちろん内容も素晴らしいので,「放課後のプレアデス」一度見てみてはいかがでしょう.

sbr-gx.jp

しかも書いているのは石井俊全さん.先生の線形代数微積分の本には本当に助けられました.僕が大学に入ってからこれらの数学を楽しんで勉強できたのも先生のおかげです!

www.amazon.co.jp

www.amazon.co.jp

また,先生の統計学の本のおかげで統計検定 2 級を S ランクで合格できました…
 \chi^2 分布や  t 分布, F 分布の定義は何か,なぜこの統計検定量がこれらの分布に従うのかが式でわかってよかったです.

www.amazon.co.jp

www.amazon.co.jp

大学入ってから石井さんの本にばっかり助けられてる…

まあ,宣伝はともかく,中身の話を簡単にしましょう.
どこまで理解できてるか怪しいですが… (-_-;)
おかしかったらっょぃ人教えてください!

まず最初の難関はテンソルが理解できるかどうかだと思うんです.
僕の理解はこうです.
僕らの視界は首を左にひねれば右に動きますよね?
つまり,"反"対方向に"変"わるです.また,当たり前ですが視点はひねった方向と"共"に"変"わるわけです.
でも,僕らが見ている世界は実体としては 1 つですよね?
これを逆に考えて,首を左にひねったら視点は同じ方向に,視界は逆の方向に動く.このようにつじつまのあった動き方をする時に,それで捉えているもの (この場合世界) は普遍の実体であることがわかるのです.それを書きあらわす数学的な道具がテンソル
相対性理論では首をひねるのが座標系を変えること,普遍の実体が物理法則です.
普遍な実体を指しているので,テンソルは物理法則を表す数式に用いるのにもってこいというわけです.
この本のテンソルの説明はとてもわかりやすいです.(と言っても僕はこれで初めてテンソルというものを学んだわけですが)
通常物理の本では成分の書き換えの規則で定義するようなのですが,この本では数学寄りらしく,基底を出して説明しています.わかりやすい.

また,高校時代にやった曲率半径の話が再登場するのは感動しましたね.
高校数学はやっぱりすごい.
それを拡張したガウス曲率,それにつながるリーマンの曲率テンソル,接続係数,計量テンソルが時空における我々の空間の曲がり具合の情報を持っていて,それが重力の式や Einstein の重力場方程式に現れるということは,重力は空間の曲がり具合であると言って差し支えないということでしょうか.
曲面人の例えもすごくわかりやすかったです.
3 次元空間内の曲面に住んでいる人は,自分が 3 次元空間にはめ込まれているなんて思わないけど,入手可能な計量テンソルからガウス曲率を計算してみると 0 ではなく,自らの世界が曲がっているのを自覚できると.ガウスはこれを指して "驚きの定理" と言ったそうです.
我々の 4 次元時空もより高次元の空間では曲がっていて,その曲がり具合が重力を生み出すってことですかね.
高次元から見たら低次元の見えないものが見える,っていうのは「正解するカド」の40 (41?) 次元の異方存在を思い出しました.
最後急展開で賛否が分かれますが見てみてはどうでしょう.

seikaisuru-kado.com


また,この本ではあまり触れてませんでしたが,Lorentz 変換で普遍な Minkowski 距離はいわば双曲線的な距離の測り方で,それは Lorentz 変換とその逆変換の対等性から出てくるという話は実は挫折する前に書いていました.よかったらご覧ください.

ysmr-ry.hatenablog.com

この中で登場する「小人さんの妄想」ってブログはものすごくわかりやすいです.
こんなのもあった.

d.hatena.ne.jp

空間の曲がり具合につながる計量テンソルが場所によって違うということは重力が場所によって違うってことですかね.
でも 1 点においては計量テンソルが Minkowski 計量な直線座標である局所ローレンツ系を取れる?
テンソルが座標系の境界,変換表ってどういうことなんですかね.っょぃ人教えて〜

生存数の分布とカノニカル分布の類似性

このような記事を拝見しました*1

www.procrasist.com

この「何ヶ月ブログが続いているか」のヒストグラムが,最近統計力学の試験勉強をしていたこともあり,カノニカル分布*2 (要するに指数分布) に似ていることに気づいたので,考えたことを綴ろうと思います.

統計力学の理解ゆるゆるなので,おかしかったら教えてください.

カノニカル分布とは

統計力学ではミクロな世界の仕組みでマクロな世界の熱力学的物理量を説明するわけですが,そこに確率モデルを用います.
例えば気体分子の数はアボガドロ数 10^{23} オーダーですから,ここに運動方程式を立式して解くというのは現実的ではありません.数が莫大すぎるのです.
そこで,その数の莫大さを逆手に取り,莫大がゆえに確定的な議論ができるように確率モデルを導入するのです.
個々の量子状態に対応する確率を重み付けし,その総和を取れば期待値 (平均) が求められます.ですが,確率を用いている以上,それには不確かさ (ゆらぎ) が存在します.しかしながら,数が莫大がゆえにそのゆらぎは限りなく小さくなり,確率といういわば "不確かな" ものを用いて計算しているにもかかわらず,実質的には確定的な物理量を議論できるのです.そのような期待値として,ミクロの理論からマクロな物理量をひねり出します.
さて,ここで,各量子状態に対して,どのように確率を重み付けするかということが問題になります.
熱平衡状態として実現するマクロな状態を司るミクロな量子状態は,全ての量子状態の中で圧倒的大多数を占めており,また,その中でミクロな量子状態が異なってもマクロには区別がつきません.
そこで一番愚直には,全ての量子状態に一様に同じ確率を重み付けするというのが考えられます.すなわち,全ての状態数を  W とすると, \frac{1}{W} で確率が一定として重み付けするのです.
これをミクロカノニカル分布といいます.
もう少し凝ると,熱平衡状態を注目している系に対して圧倒的に巨大な熱浴にくっつけて実現していると考えて,頑張ると,量子状態  i のエネルギーが  E_i として,次のような確率モデルが考えられます.

f:id:Ysmr_Ry:20170720234055p:plain

 Z(\beta) は分配関数といい,単に確率の総和を 1 にするためのつじつま合わせの定数です.
要するに,エネルギー  E を持つ確率は  e^{-\beta E} (指数関数) に比例するということです.
これをカノニカル分布といいます.

カノニカル分布と貧富の格差

このように,とある分布が非自明に指数分布に従うと言ったことは,他の場面でも起きていて,(ブログのヒストグラムもその一種ではないかと考えているのですが) そのなかで最たるものが "貧富の格差の分布" です.

www.procrasist.com

d.hatena.ne.jp

カノニカル分布は,大雑把には,全エネルギーが一定の条件下の元でランダムに分配されると,エネルギー  E を持っている確率が  e^{-\beta E} に比例するという分布 (指数分布) に落ち着くということです.
エネルギーをお金だと思うと,世の中に流通するお金の総和は一定として,ランダムに分配がなされると,その分布は指数分布となり,貧乏 ~ 平凡がたくさんいる一方,持っているお金が上がれば上がるほど,それこそ指数関数的に急激に数が減っていくという,この世の理不尽がわかると思います.
つまり,お金の例とカノニカル分布の裏には共通の理屈が走っていることになります.
これがカノニカル分布を "直感的" に捉える上で重要だと個人的には思っています.
では,その理屈とはなんなのか.
それは,「総和一定という条件付きで,一様分布にもっとも "近い" 分布が指数分布である」ということです.

カノニカル分布の情報理論的導出

まあ,全てこの資料を見ていただければ解決なのですが,

http://genkuroki.web.fc2.com/20160616KullbackLeibler.pdf

簡単に説明します.
母集団分布  q に従っている母集団から,標本を抽出し,その経験分布を  p とします.
大数の法則から,抽出する標本の数を無限に飛ばせば  q p は一致します.
しかしながら,今考えたいのは,条件付きで抽出する標本の数を無限に飛ばせば母集団分布  q に対して経験分布  p はどうなるのかということです.
そこで,2 つの分布  p,\ q の "距離" を測る,カルバック・ライブラー (Kullback-Leibler) 情報量  D(p||q) を導入します.

f:id:Ysmr_Ry:20170721002235p:plain

これがなぜ分布の距離を測るのかは,上記資料の 1.1 ~ 1.3 節を読んでいただければ,多項分布と Stirling の公式から鮮やかに導出される様子が見れると思うのですが,
要するに,経験分布がほぼ  p になる確率  r というのが, \log r = -nD(p||q) という関係で結ばれていて, D(p||q) が最小になれば  r が最大となり,その  p が経験分布として実現することを意味することになるのです.
直感的には  D(p||q) が分布  p q の距離を表しており,実際に実現する経験分布  p はこの量を最小化する (もっとも近い) 分布に落ち着くということです.
また,情報理論的な相対エントロピー S(p||q) = -D(p||q) と定義でき,これは最大化します.このエントロピーは熱力学・統計力学に登場するエントロピーと定数倍の差しかなく,エントロピー増大の法則がわかると思います.

今,考えたいのは,条件付き大数の法則ですから,つまり,条件付きで  D(p||q) を最小化するということです.
まあ,お馴染みの Lagrange の未定乗数法を使うわけです.
 p の総和が 1 であること, p によるエネルギーの期待値が内部エネルギー  U であること (期待値 (平均) ではなく,総和でも OK) を等号条件として課し, D(p||q) を最小化すると,

f:id:Ysmr_Ry:20170721004457p:plain

と,元の分布  q e^{-\beta E} をつけただけになります*3
このあたりの詳細が知りたい人は,資料の 2.1 ~ 2.2 節らへんを見ると幸せになれます.
つまり,統計力学で考えているカノニカル分布は,母集団分布  q が全て一定の一様分布になっていることがわかります.
ミクロカノニカル分布が一様分布ですから,ちょうどエネルギーの期待値が内部エネルギーと一致するようにした上で,ミクロカノニカルにもっとも近い分布を自らの拡張として選んだ,ということでしょうか.
よく統計力学の本とか授業とかで Stirling の公式と Lagrange の未定乗数法を組み合わせて導出しますが,その本質は,"分布の距離の最小化" なのです.

つまり何が言いたいの

つまり,ブログのヒストグラムは,ブログの生存数の総和がほぼ一定なのではないかということです.その上で完全に一様に分配されるのなら,上記で語ってきたメカニズムで持って  m ヶ月生き残ったブログがいる確率は  e^{-m} に比例しているのではないかと.
ブログの生存数は,気体分子の中の速度分布や貧富の格差の分布と,本質的に同じ形をしているかもしれないという話です.

*1:リンク貼らせてもらいましたが,問題があればお申し付けください (他のリンクも同様)

*2:統計力学でもっとも有用と言っても良い確率モデル

*3:これから,カノニカル分布にとって当重率の原理 (母集団分布を一様分布にするということ) は本質的ではないことがわかります.

転送行列の意味 (1 次元 Ising モデルの厳密解)

こんにちは!Ysmr-Ry です!

最近グラフに関する記事を読んだので,それに関連することを書こうと思います.
最近僕は統計力学という物理の 1 分野を勉強しているのですが,その中でグラフ的な考え方が登場したのでそれを紹介しようと思います.

それが転送行列というやつです.

僕の読んでいる本では,これを使った方法がグラフ的な説明なしにさも "当然" というかのように書かれていたので,そのあたりを補完できたらと思います.
(個人的には) このような行列の見方はあまりしないので,僕と同じように戸惑う人もいるのかなと思いまして.

今回は割と難しい (?) ので,わかってくれる人にわかってもらえれば満足という感じです ( ごめんなさいなんでも(ry )

物体の磁性と Ising モデル

統計力学の応用の中で魅力的なものに,物体の磁性があります.物体の磁性とは物体が磁石にくっつくかということです.
物体の磁性の起源は電子のスピンと言われるものです.
スピンは簡単に言えば電子が上向きか下向きかということで,スピン変数  \sigma = \pm1 で表し, \sigma=+1 ならば上向き, \sigma=-1 ならば下向きということにします.
簡単のため,ここでは 1 次元,つまり  L 個の電子が一列に並んでいる状況を考えます.
例えばこんな感じです.

f:id:Ysmr_Ry:20170705204451p:plain

上向き・下向きの矢印が電子を表しています.この場合は  L=5 です.
以下周期的境界条件を考え,この場合ならば  L=5 の列が無限に繰り返されているという状況を考えます.
すると,左から  i 番目のスピン変数を  \sigma_i と表すと, \sigma_{L+1} = \sigma_1 となります.
さて,Ising モデルは,これらのスピンのうち,隣り合っているもの同士が相互作用を持つ時を考える簡単なモデルです.
スピン配位  (\sigma_1,\cdots,\sigma_L) に対し,そのエネルギー  E_{(\sigma_1,\cdots,\sigma_L)} は以下のように定義されます.

f:id:Ysmr_Ry:20170705205232p:plain

 J は交換相互作用定数という物質固有の正の定数です.
第 2 項は磁場  \vec{H} = (0,0,H) が物質にかかった時の各スピンに対応するエネルギーです.
第 1 項が隣り合うスピン同士の相互作用に対応しています.
今,スピン変数は  \pm1 の値しか取らないため,その積  \sigma_1\sigma_2 \sigma_1, \sigma_2 が同じ値ならば 1,違うなら -1 になります.よって, -J\sigma_1\sigma_2 は同じ値ならば  -J,違うならば  +J になります.
エネルギーは低いほど安定*1ですから,隣り合うスピンが同じ方が安定になるようになっています.このように相互作用を組み込むのが Ising モデルです.

転送行列が出てくるまで

先のエネルギーの式を考えやすいように少し変形しましょう.周期的境界条件が成り立っていますから,

f:id:Ysmr_Ry:20170705210600p:plain

です.
ここで,カノニカル分布を考えます.
カノニカル分布では,エネルギー  E_i に対する確率を  \frac{e^{-\beta E_i}}{\sum_i e^{-\beta E_i}} とします.
 Z(\beta) = \sum_i e^{-\beta E_i} は確率の総和を 1 にするためのつじつま合わせの定数 (規格化定数) で,分配関数と呼ばれています.
カノニカル分布は, U = -\frac{\partial}{\partial\beta} \log Z(\beta),\ F = -\frac{1}{\beta} \log Z(\beta) のように分配関数  Z(\beta) からマクロな熱力学関数をひねり出すので,分配関数が解析的に閉じた形で表せれば勝ちなのです.
この場合の分配関数は,

f:id:Ysmr_Ry:20170705211844p:plain

と変形できます.あとはこの積を解析的に表せば良いのです.
(ちなみに  \prod \sum の積バージョンです.)
ここで, \sigma_i, \sigma_{i+1} = \pm1 しかとりませんから, f(\sigma_i,\sigma_{i+1}) のとり得る値は  f(\pm1,\pm1) の 4 通りです.
これを 2 次の正方行列で表すと,

f:id:Ysmr_Ry:20170705212503p:plain

となります.この行列  T を転送行列というのです.

転送行列の意味

さて,この行列を使って先ほどの積を計算するのですが,まず, \prod 以降の式に注目してみましょう.

f:id:Ysmr_Ry:20170705213013p:plain

となっています. f の引数に着目すると, (\sigma_1,\sigma_2) \to (\sigma_2,\sigma_3) \to \cdots \to (\sigma_{L-1},\sigma_L) \to (\sigma_L,\sigma_1) と, \sigma_1 \to \sigma_2 \to \cdots \to \sigma_L \to \sigma_1 のような値の列,道順のようなものが決まれば  \prod で積をとっている各々の項は決まり,全ての道順, 2^L 通り全てについて総和を取れば良いことがわかります.
例えば, L=5 の時, +1\to-1\to+1\to+1\to-1\to+1 という道順がありえます.
(上の図に対応しています)
この時,対応する項は  f(+1,-1)\,f(-1,+1)\,f(+1,+1)\,f(+1,-1)\,f(-1,+1)\,f(+1,+1) です.
また,道順の総数は, \sigma_i\,(i=1,2,\cdots,5) の各々が  \pm1 の 2 通りの値を取りうるので,重複順列より, 2^5=32 通りです.
それら各々に対応する項をこの要領で求め,総和を求めるということです.
これをグラフに帰着します.
下図のように +1, -1 の 2 つの頂点を持つグラフを考えます.

f:id:Ysmr_Ry:20170705215510p:plain

有限オートマトン的な.
辺を通ったら,辺の上に書いてある式を掛けると約束します.
すると, +1\to-1\to+1\to+1\to-1\to+1 とこのグラフ上を移動することが, f(+1,-1)\,f(-1,+1)\,f(+1,+1)\,f(+1,-1)\,f(-1,+1)\,f(+1,+1) の項を作ることに対応します.
さらに,転送行列  T がこのグラフの隣接行列になります.
グラフに対する隣接行列  A=(a_{ij})とは,上のグラフで言えば,頂点 1 を 1 番,頂点 -1 を 2 とすると, a_{ij} = (\mbox{頂点 i から 頂点 j に行くコスト (ここでは通ったら掛ける式)}) と定義される行列のことです.
隣接行列には便利な性質があって, A^n を考えると,その  (i,j) 成分は頂点 i から頂点 j まで,合計辺を  n 本通って行く全ての道順に対応する項の和となります*2
よって, \sigma_1 \to \sigma_2 \to \cdots \to \sigma_L \to \sigma_1 の道順は出発した頂点と同じ頂点に戻る長さ  L の道順 (これを閉路といいます) なので, T^L (1,1) 成分と  (2,2) 成分を足したもの,つまり, \mathrm{Tr}(T^L) を求めればよいのです.
ここまでが転送行列の意味 (この記事のメインテーマ) です.あとは,線形代数の知識をもとに計算をするだけです.

トレースを対角化によって求める

対角化は大学初年時で習う線形代数の中で 1 番大事と言っても過言ではない操作です.
行列  A に対する対角化とは, D=P^{-1}AP が対角行列*3となるような正則な行列  P を見つけることです.
逆に言えば, A=PDP^{-1} と表せるということです.今,転送行列  T は対称行列なので,直交行列  P を用いて  T=PDP^{-1} と対角化できます.
また, \mathrm{Tr}(AB) = \mathrm{Tr}(BA) が成り立つので,

f:id:Ysmr_Ry:20170705221848p:plain

となります.ただし,対角行列の表す線形変換は拡大変換であるため, D^L は単に  D の対角成分 (固有値といいます) を  L 乗するだけで OK です.
これが分配関数になるので,あとは具体的に転送行列の固有値を求め,分配関数から熱力学関数をひねり出せばよろしいのです.
やってみると,

f:id:Ysmr_Ry:20170705224550p:plain

と,磁化率などが不連続点を持たないため,1 次元 Ising モデルでは相転移は起こらないことがわかります.

おわり

どうだったでしょうか…
これ,決して簡単な話ではないと思うのですよ.現に Ising の学位論文ではこんなスマートな方法は使ってなくて,分配関数の求め方がかなり複雑になっているそうですし.
どいつもこいつも「え?こんなの当然わかるよね?」みたいに書いてあるのがなんかなぁと思ったのでまとめてみました.
誰かしらの役に立てば幸いです…
(統計力学勉強中なので間違えてるところがあったら教えてください)

例によって,コメント・ブクマ・スターとくれると泣きながら飛び上がって喜ぶのでよろしければレスポンスをください.

*1:エネルギーは変化する能力なので,それが低ければ変化しにくく安定.

*2:行列の積の定義を考えるとわかります. A^2 (i,j) 成分は  \sum_{k=1}^N a_{ik}a_{kj} なので,これで総和をとっている各項は頂点 i → 頂点 k → 頂点 j と行く通り数であり,k を全ての頂点について動かしているので,頂点 i → 頂点 j へと行く長さ 2 の全ての道順を試していることになります.以降帰納的に示せます.

*3:対角成分のみがあるような行列, \begin{pmatrix}\lambda_1 & 0 \\ 0 & \lambda_2\end{pmatrix} のような行列のことです.