ゆるふわブログ

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

ボロノイ図を描く! ~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} のような行列のことです.

素因数分解の可視化

まずはこれをご覧いただきたい.

See the Pen Prime Factorization by Ysmr-Ry (@As2) on CodePen.

おわかりいただけただろうか? 素因数分解を可視化しているのである.
この記事ではこの可視化のルールと気づいたこと,またおまけとして素因数分解について考えたことを綴ろうと思う.

ちなみに元ネタはこれである.

matome.naver.jp

素因数分解ってなんだっけ (真顔)

素因数分解は,(たぶん) 中学校ぐらいで習うことに加え,下手したら小学生でも知っているので大丈夫だとは思うが一応おさらいしておく.
ある整数 n にはそれを割り切る数がある.例えば,6 であれば,2 や 3 は 6 を割り切る.
そのような数を約数という.
だが,例えば 7 のように 1 か自分自身の 7 でしか割り切ることのできない数もある.それを数の素となる数ということで素数と呼ぶ.
素数を小さい方から書き出すというのをやったことがあるのではないだろうか.
2, 3, 5, 7, 11, 13, 17, 23, 29, …
整数はなんでも,素数を掛け合わせた形で表すことができる.その形の中での素数 1 つ 1 つのことを素因数と呼び,そのような形に分解することを素因数分解と呼ぶ.
例えば,12 であれば 12 = 2×2×3 と表すことができ,これを素因数分解と呼ぶ.素因数は 2, 3 である.

素因数分解の可視化のルール

この可視化のルールは蓋を開けてみれば至極単純である.
まずは一番単純な場合,可視化する数が素数である時を見てみよう.
例えば 5 である.

f:id:Ysmr_Ry:20170621234858p:plain

素数の場合は,円の中心が正多角形になっている.この場合は正 5 角形に並んでいる.
次に素因数を一つ増やして,35 = 7×5 を図示してみよう.

f:id:Ysmr_Ry:20170621235143p:plain

5 の時の図形 (上の画像) が正 7 角形上に並んでいることがわかる.
次は 385 = 11×7×5 だ!

f:id:Ysmr_Ry:20170622000046p:plain

正 11 角形上に 35 の図形が並んでいるのがわかるだろう.
以下同様だ.5005 = 13×11×7×5 の場合も一応見てみよう.

f:id:Ysmr_Ry:20170622000400p:plain

385 の図形が正 13 角形上に並んでいるのがわかるだろう.
このようなマトリョーシカのように入れ子になっている構造を (難しいが) 再帰的構造という.

自分自身が入れ子なフラクタル

もしも,整数を構成する素因数が全て同じ場合 (例えば 16 = 2×2×2×2 など) の場合はどうなるだろうか.
先ほどのルールに従えば同じ形が入れ子になっていくはずである.
無限に自分自身が入れ子になっている図形をフラクタルという.
つまり.フラクタルっぽくなるはずである.
では,実際に見てみよう.
1024 = 2×2×2×2×2×2×2×2×2×2 である.

f:id:Ysmr_Ry:20170622002548p:plain

これはシェルピンスキーのカーペットという有名なフラクタルのようになっている.
また,729 = 3×3×3×3×3×3 の場合である.

f:id:Ysmr_Ry:20170622002828p:plain

これもまたシェルピンスキーのギャスケットという有名なフラクタルである.
名前は知らずとも,このような形を見たことがある人も多いのではないだろうか.

素因数分解ってなんの意味があるの (おまけだから飛ばしてもいいよ)

素因数分解することで嬉しいこととはなんだろうか.
僕の中のイメージはこうだ.
素数とは の元 なのである.素因数分解とはいわば,数の化学式 なのである.
化学では,有限な種類の元素 (現在名前がついているのは 118 種) がどのように集まって化合物を作っているのかが重要である.
同じように,数がどのような素数でできているかが重要なのである.
(素数は無限にあるという違いはあるが…)
素因数分解をするといろいろなことがわかる.例えば, 360 = 2^3\times3^2\times5^1素因数分解すると,それぞれの素因数の肩に乗っている指数 (3,2,1) にそれぞれ 1 を足してかけると (3+1)(2+1)(1+1) = 24 であり,これが 360 の約数の個数である.また, 12 = 2^2\times3^1 など,それぞれの指数が 360 のものより小さい場合はその約数になる.(360 は 12 の約数) 約数の総和は (1+2+4+8)(1+3+9)(1+5) = 1170 である.
また,2 つの数の比較にも使える. 360 = 2^3\times3^2\times5^1 300 = 2^2\times3^1\times5^2 を考えよう.指数の列を書き出そう.
360 … (3,2,1)
300 … (2,1,2)
この中で,同じ列にあるものの中で小さい方をとって行くと,(2,1,1) である.これを数に直すと, 60 = 2^2\times3^1\times5^1 となり,360 と 300 の最大公約数である.
また,同じ列で大きい方をとって行くと (3,2,2) で, 1800 = 2^3\times3^2\times5^2 となり,360 と 300 の最小公倍数である.
このあたりは現在の 数学 A の教科書に載っていると思われるので,本棚から引っ張り出してみるのも一興だろう.
ちなみに,例えば 360 の約数の積は約数の個数 24 を用いて  360^{\frac{24}{2}} である.

ベクトルを知っているとこのようなことも言える.
最大公約数が 1 である場合,互いに素というが,a と b が互いに素の時,a ⊥ b と表す.
これは,指数を並べたものをベクトルとみると,内積が 0 となるからである.
例えば, 12 = 2^2\times3^1\times5^0\times7^0 35 = 2^0\times3^0\times5^1\times7^1 を考えよう.
指数を並べると,2, 3, 5, 7 に対応する指数を順に並べて,12 … (2,1,0,0), 35 … (0,0,1,1) である.
これらの同じ列にある数をかけて足し合わせると,2×0 + 1×0 + 0×1 + 0×1 = 0 であり,内積が 0 なので "直交" していると言える.よって,12 ⊥ 35 と表記する.
このようなベクトル的な考え方は「数学ガール」を参照されたい.

「互いに素」という概念

www.hyuki.com

素因数分解は暗号にも利用されているようである.

おわび

深夜なのでテンションがおかしくてすみません.短くあっさりと仕上げるつもりだったのですが,主に最後の方が冗長になってしまいましたすみませんなんでもし(ry
素因数分解は中学生でも知ってますが,それをこのように美しく Visualize できること,また,ベクトル的な見方など面白い捉え方もあるのだということを知ってもらい,素因数分解に対して楽しいイメージを持っていただければ幸いです.
何か不明な点があれば,どんな些細なことでも結構ですので,コメントで聞いてください.
ブクマやスターなどつけていただくと泣いて喜びます.

凸ってなに?(数学)

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

ふと以前に見たキラキラネームについての動画を思い出しました.
凸って書いてテトリスって読むそうです.
これは特にひどい例ですが,名前つけるならちゃんと考えなきゃいけないよなぁと.
考えてみればみんなうまいこと子供の名前つけてるなぁって.昔小学校の授業で自分の名前の意味を調べる的なのやった気がします.
なんかアイデア集的なのがあるんですかね.

kazuyahkd.com

あと凸って言葉で印象に残ってることといえば… あれは中学 2 年生ぐらいだったか… ちょうど教育実習の先生がきていて,数学の単元が 2 次関数だったのです.
そこで先生は, y=ax^2 a>0 ならば「下に凸」, a<0 ならば「上に凸」というと強調していたのを覚えています.

実はこの,上に凸・下に凸,さらにそれが初めて出てくるのが 2 次関数であるということ,は大学に入ってからの微積分や熱力学で実は出てきていて,大事なことだなと思ったので,数学でいう凸性について僕が知っている範囲で面白いことを書き連ねていこうと思います.

数学でいう凸って?

凸というと,オウトツのトツ,デコボコのデコのことですが,数学的にはこうとらえます.
グラフ上の二点を結んだ線分が常にグラフの上側にあるような関数が下に凸な関数,逆に常にグラフの下側にあるような関数が上に凸な関数です.
図にするとこんな感じです.

f:id:Ysmr_Ry:20170609222327p:plain

f:id:Ysmr_Ry:20170609222347p:plain

数式にすると,線分の端点の  x 座標を  x_1,\ x_2 とおくと,線分上の点は  0 \leq t \leq 1 なる実数  t を用いて  t:(1-t) に内分した点と捉えられるので,次のようになります.
下に凸

f:id:Ysmr_Ry:20170612165844p:plain

上に凸

f:id:Ysmr_Ry:20170612165856p:plain

凸から相加・相乗平均の不等式がわかる?

 f(x) = \log x の場合を考えましょう. \log x は上に凸な関数なので*1,先ほどの 2 つ目の不等式を使って,特に  t = \frac{1}{2} のとき,

f:id:Ysmr_Ry:20170612170855p:plain

となります.第 1 式から第 2 式へは  \log の性質  \log a+\log b = \log ab,\ k\log a = \log a^k を,第 2 式から第 3 式へは  \log x (\log x)' = \frac{1}{x} > 0\,\,(x>0) から単調増加関数であることを利用し, \log をかけた不等式が成り立っていれば  \log を外しても成り立つことを使いました.例えば,2 人のテストの点数がそれぞれ  x_1,\ x_2 点だったときに,その平均にはいろいろあって,普通使う平均である  \frac{x_1+x_2}{2} を相加平均,ちょっとマイナーな  \sqrt{x_1x_2} という平均を相乗平均というのでした.上の不等式はそれらの関係,すなわち,(相乗平均)  \leq (相加平均) という関係を示しています.人数が  n 人になっても同様にできます.相加・相乗平均の不等式は  \log の凸性から出てくるのです.

凸で極大・極小値か見極める! (1 変数)

 f(x)極値を求めるとき, x_0極値をとるならば  f'(x_0)=0 ですが,その逆は成り立ちません.例えば  y=x^3 x=0 のように変曲点になっている場合があるからです.

増減表をかけば極値か変曲点かはっきりするのですが,他の方法として, f(x) の 2 階微分を用いて, f''(x) > 0 ならば極小値, f''(x) < 0 ならば極大値, f''(x) = 0 ならば変曲点であると判定できると数 II の教科書に確か書いてあったと思います.実はこの判定法がやっていることが「凸で極大・極小値か見極める」ことなのです.
この判定法は,極値の候補である  f'(x)=0 を満たす点の付近にぴったりとフィットする放物線を考えています. f(x) x=x_0 付近にフィットする放物線は  y=f(x_0)+f'(x_0)\,x+\frac{f''(x_0)}{2}x^2 と求められます.特に今回は  f'(x_0)=0 が成り立っているので,1 次の項は消えます.
極小値付近は形的に下に凸なはずですから,その付近にフィットする放物線の 2 次の係数は正なはずです.逆に極大値は負です.変曲点においては 2 次の係数も 0 なので単なる  y=f(x_0) になってしまいます.

 y=x^4-2x^2 のグラフです. x=0 で極大値, x=\pm1 で極小値を取ります.紫がフィットする放物線です.

 y=x^3 x=0 付近ではフィットする 2 次関数は単なる  y=0 になってしまいます.

つまり, y=ax^2\,\,(a>0) ( x=0 で極小) と  y=ax^2\,\,(a<0) ( x=0 で極大) と  y=a (変曲点) が基本形になっているのです.
「凸を知りたければ 2 階微分を見よ!」です.

グラフが位置エネルギーを表しているとすると,極小付近でフィットする放物線を求めることは対応する平衡点付近の運動を単振動に近似することに相当します.
下に凸なのがポイントで,これを安定な平衡点といいます.

もっと知りたい人のために

イェンゼンの不等式

 n 人の場合も相加・相乗平均の不等式を求められるといいましたが,どうやるのか詳しくみましょう.
3 人の場合を考えます.先ほどまでの凸の定義では,関数  f(x) 上の 2 点を結んだ線分を考えていましたが,今度は 3 点を結んだ三角形を考えます.

f:id:Ysmr_Ry:20170612172757p:plain

するとその三角形内の全ての点は同じ  x 座標の  f(x) 上の点より下にあることがわかります.三角形内部の点は重心座標*2を用いて  t_1x_1+t_2x_2+t_3x_3\,\,(t_1+t_2+t_3=1) と表せるので,不等式は,

f:id:Ysmr_Ry:20170612173347p:plain

となります.特に  f(x)=\log x,\ t_1=t_2=t_3=\frac{1}{3} のとき,

f:id:Ysmr_Ry:20170612173649p:plain

と導出することができました. n 人の場合も同様に, n 個の点を結んだ  n 角形内の点を和が 1 となるような  n 個の重心座標  t_i\,\,(i=1,2,\cdots,n) で表せばよいことがわかります.これをイェンゼンの不等式といいます.

凸で極大・極小値か見極める! (2 変数)

先ほどは変数が  x と 1 つだけでしたが,ここでは 2 変数にして  f(x,y) を考えます.次元の数を 1 つ増やして 3 次元とし, z=f(x,y) を満たす点を集めればこれのグラフを書くことができます.極値の候補になる点は,1 変数の場合と同様に  f_x(x,y) = f_y(x,y) = 0 と求められます*3.この候補の周囲の状況は先ほどの 1 変数の場合と同じように,極大か極小か鞍点かその他かに大別できます.そのうち極大,極小,鞍点にはやはり同様に基本形となる曲面があります.それらは 2 次形式という 2 変数における 2 次関数に対応するもの*4で,次のような式と形をしています.

極小:  z=x^2+y^2
f:id:Ysmr_Ry:20170612183511p:plain

極大:  z=-x^2-y^2
f:id:Ysmr_Ry:20170612183531p:plain

鞍点:  z=x^2-y^2
f:id:Ysmr_Ry:20170612183545p:plain

 f(x,y) 上の点  (x_0,y_0) 付近にフィットする 2 次形式は,行列を用いて次のように表されます.

f:id:Ysmr_Ry:20170612210422p:plain

今, f_x(x,y) = f_y(x,y) = 0 なのでやはり 1 次の項は消え, f(x_0,y_0) は所詮  z 方向の平行移動なので,重要なのは 2 次の項です.
この 2 次の項のあらわす 2 次曲面を  z 軸に平行な回転軸で回転して,上にあげた基本形のどれと同じかを見極めれば判定ができます.
詳細は省略しますが,行列の部分を対角化した固有値 \lambda_1,\ \lambda_2 とすると,回転後の式は  z=\lambda_1x^2+\lambda_2y^2 になるので,  \lambda_1 \lambda_2 > 0,\ \lambda_1 > 0,\ \lambda_2 > 0 ならば極小 ( z=x^2+y^2) 型,  \lambda_1 \lambda_2 > 0,\ \lambda_1 < 0,\ \lambda_2 < 0 ならば極大 ( z=-x^2-y^2) 型, \lambda_1\lambda_2 < 0 ならば鞍点 ( z=x^2-y^2) 型とわかります.

これが大学の微積分でやる Hessian による極値判定の本質です.1 変数の場合の綺麗な拡張になっています.

Legendre 変換と凸 (おまけ)

Legendre 変換とは, y=f(x) を通常  (x,f(x)) という点の集合と捉えるのに対し, x における接線の傾き  f'(x) とその  y 切片  f(x)-f'(x)\,x の組  (f'(x), f(x)-f'(x)\,x) で表わそうという変換です.(いわば "接線座標") このとき,グラフに凹凸があると同じ接線の傾きに複数の点が対応してしまい面倒です.それに対し,凸であると,1 対 1 の対応になるので都合が良いのです.熱力学では状態曲面  U(S,V) が凸な曲面であることを利用し,Legendre 変換の関係で  G,\ F,\ H,\ U を結んでいます.

*1:グラフを描いても確かめられますが, (\log x)''=-\frac{1}{x^2} < 0\,\,(x > 0) からもわかります. y=ax^2\,\,(a<0) と同じ状況なので上に凸です.

参考:  y=\log x のグラフ

*2:くわしくは
ysmr-ry.hatenablog.com
をごらんください .

*3: f_x(x,y) f(x,y) y を定数と見て  x微分したもの.(偏微分)

*4:これは統計とかこれとかを見るとわかるかもしれません.
ysmr-ry.hatenablog.com