2018-07-21

幾何学的モーメントの不変式とshape analysis

3DCGでよくやるように、何か(向き付けられた)閉曲面があって、その曲面は、三角形に分割(単体分割)されているとする。この時、曲面内部の領域の体積を計算するには、適当な一点Oを取って(特に理由がなければ原点でいい)、Oと各三角形がなす四面体の符号付き体積の和を計算すればいい(三角形の頂点の順序が、面の向きと合っていなければならない)。3DCGで使われるモデルでは、境界が連結でなく、複数の閉曲面からなることもあるけど、その場合でも、この体積計算は有効。

同様の計算は、三次元空間上で定義された関数を曲面内部の領域で積分するのに使える。関数として、最も単純な単項式を選ぶと、これは、曲面内部の領域の幾何学的モーメントで、体積は特に0次のモーメント。この積分の計算自体は、大学一年生の算数で、ちょろいけど、計算結果が書いてあるものを見つけられなかったので、以下に書いた

単体の幾何学的モーメントの計算

https://vertexoperator.github.io/2018/04/30/simplex_moment.html

実際に計算する時、世の中で出回ってる3Dメッシュデータは、non-manifold edge(本来、全てのedgeは丁度2つの異なる面にのみ含まれるべき)を持つことがよくあるので、気をつける必要がある



幾何学的モーメントは、曲面を合同変換した時に、どのような変換を受けるか分かるので、幾何学的モーメントを組み合わせて、形状不変量を作ることができる。これは、Huモーメントの場合と同じ考え方。

Image moment

https://en.wikipedia.org/wiki/Image_moment

Huモーメント不変量は、通常、二次元で定義されていたけど、三次元でも同様に定義できる。世の中には、形状データベースみたいなものを作りたいという需要もあるらしいから、そういうので、使えそうな気もする。何に使うのか、よく分からないけど。物体認識とかで使えそうな気もするけど、一枚の2D画像から3Dデータを再構成するのはしんどいし、自然界には、同一形状の物体というのはあまり存在しない(岩とか木とか)

低次の場合を考える。0次のモーメントは体積でこれは自明に合同変換で不変な量。一次のモーメントは、重心座標(と体積の積)を表すので、これから合同変換で不変な量は作れない。二次のモーメントは、重心が原点と一致するように動かしておくと、残りの自由度は回転のみで、モーメントを適当に並べると、2階の対称テンソルをなす。これを回転で動かすと、対角化できる。固有関数は互いに直交する三軸で、固有値は、各軸方向の"広がり"を表す。つまり、二次のモーメントは、物体を楕円体で近似した時の形状を表すと思える。回転不変な量は、対称テンソルの固有値の対称関数で、トレースや行列式などを含む。

注)ここの対称テンソルは、正確には、群SO(3)の自然表現の対称テンソル積表現空間の元であることを意味する。SO(3)の表現空間には、SO(3)作用で不変な内積が定数倍を除いて一意に存在(コンパクト群の有限次元表現がユニタリ表現になるのと同様の論法)して、特に、双対空間との同一視ができるので、対称テンソルと対称行列が同一視できる。モーメントが対称テンソルになるというのは、物理では、多重極モーメントとかで使う考え方

従って、2次までの幾何学的モーメントは、0次のモーメントで、大体の大きさ、一次のモーメントで、おおよその位置、二次のモーメントで大雑把な形状(平べったいとか丸っこいとか、細長いという程度の)を表現しており、人間の直感的な捉え方に近い感じがする。より詳細な構造の情報は、もっと高次のモーメントに含まれる。高次のモーメントからも、合同変換の不変量が作れて、これらは"実質的に任意の形状"を分類するのに十分な不変量を与える。尤も、これらの不変量を計算する不変式を一般的に求めるのは、多分難しい


どれくらい難しいか。高次のモーメントから幾何学的Huモーメント不変量を作る問題は、並進自由度は簡単に除けるので、残る回転自由度に関する問題となり、数学的には、

3次元のHuモーメント不変量の計算:「SO(3)の自然表現の高階対称テンソル積表現から、表現空間上の多項式関数で、SO(3)作用で不変なものを見つける」

という形に定式化される。一方、数学では、19世紀に不変式論が研究され、そこでの主要な問題は、現代の言葉では、

19世紀不変式論の基本問題:「SL(2,C)の自然表現の対称テンソル積表現(既約表現になる)の表現空間上の多項式関数でSL(2,C)作用で不変なものを見つける」

というものだった(当時は、既約表現という概念もテンソル積という概念も定式化されてないので、2元n形式へのSL(2,C)作用という形で理解されていた)。SO(3,C)とSL(2,C)はLie環を取れば同型であり、こういう風に定式化すると、3次元のHuモーメント不変量を決定する問題と、19世紀の不変式論で扱われてた問題が非常に似た種類の問題だと分かる。Huモーメント不変量の決定のほうが、次元が大きい分、難易度が高そうに思える。ところで、後者の問題は、多分、殆どの数学者が特に重要な問題じゃないと考えるようになって久しく、現在でも、一般的な答えが分かっているわけではない(確か12次くらいまでは、不変式の生成元が決定されていた気がする)


(離散)曲面上のラプラシアンは、1993年のPinkallとPolthierの論文以来、色々な問題で、よく使われるようになった。cotangent formulaで検索すれば、沢山解説が出てくる

Computing discrete minimal surfaces and their conjugates

https://projecteuclid.org/euclid.em/1062620735

曲面の形状の不変量を得るのに、ラプラシアンの固有値と固有関数を見るのは、自然に思える。このような方法は、spectral shape analysisという名前が付いている程度には、ポピュラーらしい。けど、例えば、以下の論文のFig4とか見ると、幾何学的Huモーメント不変量ほど、直感的な情報を与えてくれそうな感じはしない。

Spectral Mesh Processing

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.229.4191


うまいことやれば、2つの曲面の剛体位置合わせも多分できる。2つの曲面が完全に合同であれば、適当な合同変換で、幾何学的モーメントが一致するようにできる。簡単のため、二次元で考えると、座標のアフィン変換

x’ = ax + by + p

y’ = cx + dy + q

は6つの量a,b,c,d,p,qで定まるので、最低でも、二次までのモーメントを見る必要がある。アフィン変換に対して、二次までの幾何学的モーメントは

M_{00}’ = (ad - bc)M_{00}

M_{10}’ = (ad - bc)(aM_{10} + b M_{01} + p M_{00})

M_{01}’ = (ad - bc)(cM_{10} + d M_{01} + q M_{00})

M_{20}’ = (ad - bc)(a^2 M_{20} + 2ab M_{11} + b^2 M_{02} + 2ap M_{10} + 2bp M_{01} + p^2 M_{00})

M_{11}’ = (ad - bc)(ac M_{20} + (ad+bc) M_{11} + bd M_{02} + (cp + aq)M_{10} + (dp + bq)M_{01} + pq M_{00})

M_{02}’ = (ad - bc)(c^2 M_{20} + 2cd M_{11} + d^2 M_{02} + 2cq M_{10} + 2dq M_{01} + q^2 M_{00})

のように変換する。2つの曲面が合同であれば、0次のモーメントは等しいはずだけど、多分殆どの場合は、浮動小数点誤差のために、完全には一致しない。拘束条件として、ad-bc=1を課すか、ad-bcのスケールも不定とするかは問題に依存する選択だと思う。


合同変換であれば、ad-bc=1かつa=dかつb=-cである。この条件を課した上で、互いの幾何学的モーメントがなるべく一致するように、パラメータa,b,c,d,p,qを決定すれば、剛体位置合わせができる。一般に、解は一つとは限らない(例えば、球とかの場合)。何らかの評価関数を決めて最小化するというのが一番オーソドックスに思いつく。


というようなことを考えたけど、差し当たって、何かに使おうと思ってたわけではないので、本当に有用かどうかは知らない

IgnatowskiによるLorentz変換の導出のvariation

Lorentz変換の代数的導出

https://vertexoperator.github.io/2018/07/05/ignatowski.html

というのを書いた。

正しいLorentz変換の式を書いたのは、Larmorが最初っぽい(1897年)けど、電磁場と電荷・電流の変換まで考えたのはLorentzらしい。特殊相対論では、光速度不変の原理と特殊相対性原理から、Lorentz変換を説明する。1910年頃、Ignatowskiという人が、光速度不変の原理を仮定せずに、Lorentz変換の導出を与えた。と言っても、光速度が決まるわけないので、任意パラメータが一個残る(上のリンク先の定数γ)。このパラメータγは、大きく分けると、正か0か負になり、それぞれ、対応する変換群は、SO(3,1),SE(3),SO(4)となる。場の理論としては、それぞれ、相対論的場の理論、ガリレイ不変な場の理論(※)、Euclidean field theoryが対応する

(※)日本語的には、非相対論的場の理論だと、対称性がISO(n,1)じゃないやつは全部当てはまりそうな気がする

Ignatowskiの論文は、不完全だったとかいう話もあるけど、問題はすぐにfixされて、それ以後、この方法の色んなvariationが出ている。この方法に対するPauliの評価は"from the group theoretical assumption it is only possible to derive the general form of the transformation formulae, but not their physical content"という感じで、特に高くはない。個人的には、光速度不変の原理は、原理と呼ぶには、dirtyすぎる気もするけど、Einsteinの特殊相対論への貢献が、なくなってしまう


で、たまたま導出を見たところ、見通しが悪いと感じたので、自分で考えなおしたのが、上の話。見通しはよくなったと思うけど、あんまり初等的でもなくなった。上の説明だと、一次元formal group law(FGL)の可換性は、本質的に効いていて、回避する術はなさそうに思える(適当な係数環の下での、一次元formal group lawの可換性の証明は、難しいわけではないけど、何も知らずに自分で思いつくのは割と難易度高めな気がする)。Ignatowskiの時代には、formal group lawはなかったし、他の証明見ても、formal group lawは出てこないので、通常は、何か無意識に仮定している物理的条件が存在してるんじゃないかと思ったけど、よく分からなかった。まぁ、どうでもいいけど


Fizeauが1851年に実験によって確認した速度合成則は、今から見ると

f(u,v) = u + v - ¥frac{u^2 v}{c^2}

だったわけで、可換でないだけでなく、結合的でもない。可換にする最も安直な方法は

f(u,v) = u + v - ¥frac{u^2 v + v^2 u}{c^2}

とすることで、実験結果とは矛盾しないけど、結合則は満たさない。とか考えていくと、発見的に、正しい速度合成則にたどり着けても、よさそうなもんだけど、そうはならなかった。1851年には、FGLの概念とかなかったので仕方ない

フィゾーの実験

https://en.wikipedia.org/wiki/Fizeau_experiment#Fresnel_drag_coefficient


ところで、Lorentz変換というと、SO(3,1)と思うけど、フェルミオン場との相互作用の記述にスピン接続を使うことを考えると、Spin(3,1)が"本体"じゃないかという気もする。普段、Lie環でばっか考えるので、気にしたことなかったけど。

CUDAの共役勾配法サンプルが遅かった話

https://gist.github.com/vertexoperator/e2550b012978e87edc2454a77f6c70bd

鉄剣は必ずブロンズソードより強いのか

RPGだったら、ブロンズソードより、鉄剣の攻撃力が高くて高性能となるとこだけど、現実の人間は、HPと防御力が低すぎるので、石斧だろうが、ブロンズソードだろうが、当たりどころ(切られ所)が悪いと、簡単に死ぬ。

もう少し真面目な話としては、春秋戦国時代の中国では、岩鉄鉱石か沼鉄鉱石か知らないけど、鉄鉱石を一度溶融して(一酸化炭素で)還元し、出来た銑鉄を、そのまま鋳造していたらしい。こうして作られた鉄製品は、炭素含有量が高く、硬いけど脆くなると一般に書かれている(※)。それで、こうして作った鉄製の武器は脆いので、「悪金」と呼ばれて、農具などに利用されていたらしい。中国の製鉄技術の起源は、外来のものっぽく、初期の頃から、他の地域と同じように、半溶鉄を低温還元した後、鍛造して作った(鉄を溶かせない場合は、鋳造はできない)ものもなかったわけではないらしい。何故、鋳造に拘ったのかは謎。鋳造する方が生産効率は圧倒的に高いんだろうけど、それが理由かは分からない

(※)鋳鉄という用語も、銑鉄とほぼ同じ意味で使う。多分、鋳鉄の元来の意味は、鋳造用の鉄だったのだと思うけど、現代の工業製品としての鋳鉄は、大抵、1%以上のSiを含んでいるのに対して、紀元前の中国の鉄は、Siを1%未満しか含んでいない。Siの有無が物性にどう影響するのかは知らない。また、普通の鋳鉄も、組織構造によって、ねずみ鋳鉄/ねずみ銑鉄と白鋳鉄/白銑鉄などに分類されている。両者は、組織的には結構違うので、区別すべきなのだと思う。鉄関係の用語は、色々紛らわしい

春秋戦国時代の武器の鉄製と青銅製の割合がどれくらいだったかは分からない。中国では、紀元前500年頃作られたらしい青銅剣である越王勾践剣というのが発掘されていたりする(勾践は紀元前473年に上海のあたりにあった呉を滅ぼした越の王。越は呉の南にあって、この二国は呉越同舟の故事の起源)し、秦の兵馬俑からも、青銅製の剣が発掘されているらしいから、青銅製の武器も広く使われていたと思われている。始皇帝は、鉄官という、何をするのか知らないけど、鉄の生産に関連する役職を作ったとされる。キングダムなどの漫画で出てくる干将・莫耶は伝説では鉄剣ということになっている

歴史的なことは古い話なので、よく分からないけど、鉄剣は、ブロンズソードより優れていると無条件には言えない。といって、硬いとか脆いとか感覚的なことを言っても水掛け論なので、剣の性能を測る材料的なパラメータが知りたい、と思った。当然、剣の性能を決めるのは材料以外の要素もあって、大きく分厚く重く作れば頑丈になるだろうけど、そのへんは無視する。古代ローマの鉄剣とか日本刀は、複数の異なる性質の鉄を鍛接して作っていたらしいから、単一の材料パラメータだけ見ればいいというものでもないかもしれないけど、そのへんも面倒なので無視する


日本刀は、高硬度で高靭性とか書いてあったりする。Wikipediaを見ると、"靱性とは、物質の脆性破壊に対する抵抗の程度、あるいはき裂による強度低下に対する抵抗の程度のことで、端的には破壊に対する感受性や抵抗を意味する。材料の粘り強さとも言い換えられる"と書いてある。何を言ってるんだって感じだけど、脆さ(脆性)の反対の性質らしいことが分かる。英語だと、靭性はtoughnessらしいので、分かりやすい。

金属の応力歪み曲線は、変形量が小さい時は、弾性的に振るまう領域があり、ある程度まではHookeの法則が成立する。その後、変形量を大きくしていくと、弾性的に振る舞いつつもHookeの法則の当てはまりは悪くなり、やがて塑性変形する領域に移行する。更に変形していくと、やがて破断を迎える。この応力歪み曲線を、歪み0から破断点まで積分した量は、圧力=単位面積にかかる力=単位体積あたりのエネルギーと同じ次元を持ち、材質を破壊するのに、どれくらいのエネルギーが必要かという指標になる。これが大きければ折れたり刃こぼれしにくいということになるので、武器に使う材質としては、このエネルギーが大きい方がいい。

弾性変形する領域がなくて、塑性変形する領域が大きいと、ちょっと強い力を加えただけで、変形したまま戻らないということになって、それはそれで困るけど、このへんの問題は、あとで考える。また、材質によっては、塑性変形する領域が殆どなくて、弾性変形できる限界を迎えると、いきなり壊れたりする。これは脆性破壊と呼ぶ。金属ではないけど、ガラスとかは脆性破壊を起こす材料。これは、一般的な「脆い」というイメージに対応する。十円玉(銅・亜鉛・錫の合金)とかは(曲げたことないけど)多少曲げても、いきなり折れたりはしない。何回も曲げたり戻したりしてると、金属疲労を起こして折れやすくはなるけど、それは今は考慮の対象としない

この応力歪み曲線を積分して得られるエネルギーは破壊エネルギーとか破断エネルギーとか呼ぶらしい。破断エネルギーを大きくするには、破断する歪み量は大きく、かつ応力歪み曲線の最大応力値が大きいことが望ましい。前者は、伸び率、後者は引張り強度で定量化されることが多い。ゴムなんかだと、伸び率は、数百%に達しうるけど、鉄系合金だと、数十%がいいとこのよう。現実に剣を引っ張ったりすることは殆どないと思われるので、曲げ強度という量を使うほうが、より適切かもしれない


硬さ。材料界隈では、「硬い」と「固い」は、違う用語なのかもしれない。一般的な硬いというイメージは、変形しにくいということで、変形のしにくさを測る指標としては、Young率(縦弾性係数)がある。ゴムのYoung率は、1?10MPa程度とされ、銅が125GPa、純鉄では205GPa、ダイヤモンドでは1000GPaとなるので、常識的な「硬い」という認識に対応しているように見える。応力歪み曲線では、Young率は、歪みが0付近での曲線の傾きの大きさとなる。まぁ、ゴムのように簡単に変形するようでは剣の素材として使えない(武器としては、鞭とかあるけど)けども、実際のところ、殆どの金属では、Young率が50〜200GPaくらいの範囲なので、金属間での差異が問題になる要素とは言い難い(銅と鉄で2倍も違うと言えなくもないけど)。この意味での「かたさ」は、英語では、stiffnessと呼ばれているけど、紛らわしいことに、いわゆる、工学で使われる硬度とは別の概念となっている。

通常、ダイヤモンドが世界一硬いという時には、モース硬度のことを指している。モース硬度は、「ひっかいたときの傷のつきにくさ」と定義されているけど、定量的な尺度とは言い難い。この硬さは、英語では、hardnessと呼ばれて、stiffnessとは区別されている。硬度には、モース硬度以外に、色んな指標があり、材料界隈では、ビッカース硬度やブリネル硬度などがよく使われるっぽい。これらの硬度は定量的でもある。ビッカース硬度は、「ダイヤモンドでできた剛体(圧子)を被試験物に対して押込み、そのときにできるくぼみ(圧痕)の面積の大小で硬いか柔らかいかを判断する」らしい。ブリネル硬度も、なんか似たようなもので「直径D の球形の金属球を圧子として、圧子を試験面にP の力で一定時間押し当てた後、荷重を除いたあとに残った永久くぼみの面積を測定する」らしい。どちらも、押し込み硬さというものに分類されるよう。

それらしい説明とか見ると、硬さは、工業量(複数の物理的性質が関与する量で、測定方法に依存する)とか書いてあったりする。ビッカース硬度やブリネル硬度の測定の仕方を見ると、材質に塑性変形するまで圧力をかけて、塑性変形が開始する応力を決定しているように見える。塑性変形の開始点は、弾性限界として知られるけど、これを特定するのは難しい。そもそも、この名前は、弾性変形する領域と塑性変形する領域が明確に分かれる印象を与えるけど、多分、そういう理解は正しくない。

材質によって(軟鋼など)は、歪みが大きくなるのに引っ張り応力は下降する降伏現象が見られることがあり、その場合、降伏中の最大応力を、降伏応力や降伏強度(yield strength)と呼ぶ。明確な降伏点を持たない材料の場合は、耐力というものが代用として使われる。正確には、耐力は、永久に残る歪みの大きさを指定しないと定まらないので、0.2%耐力や0.5%耐力などという形で決められる。耐力も、降伏強度(yield strength)と呼ぶらしい。降伏強度は、弾性限界より特定しやすいので、よく使われる。硬度は測定の仕方は明白だとしても、何を測定しているのかよく分からないので、降伏強度の方が、何かを考える時には、便利な量だと思う(慣れている人にとっては、硬度が加工のしやすさの目安になったりするらしいけど)

ビッカース硬度やブリネル硬度の単位は、N/mm^2でMPaと同じ。ビッカース硬度とブリネル硬度の換算表というものがあって、それを見ると前者のほうが後者より値が多少大きくなるらしい。おそらく、圧子を押しこむ時、圧子の方も変形するだろうから、圧子には、なるべく変形しにくい物質を使う方が、望ましいはず。なので、圧子としてダイヤモンドを使うビッカース硬度の方が降伏強度に近いかもしれない。とはいえ、降伏強度は、引張強度より必ず小さいはずだけど、引張強度より硬度が大きくなっている例もある。hardnessは塑性変形の起きにくさの指標であり、上の方で「弾性変形する領域がなくて、塑性変形する領域が大きいと、ちょっと強い力を加えただけで、変形したまま戻らない」とか書いた話と関連する量でもある。

ビッカース硬度は、ダイヤモンドが10000前後(物や測定によって、ばらつきが大きいようだけど)なのに対して、純鉄(工業用純鉄は、炭素含有量が、フェライト相に固溶する最大量0.0218%以下の鉄・炭素二元合金という定義)では100程度らしいので、この基準でも、ダイヤモンドは硬い。純鉄の降伏点は98MPaらしいので、ビッカース硬度とオーダーは合っている。鉄炭素化合物Fe3Cはセメンタイトと呼ばれ、そのビッカース硬度は約1100らしいので、かなり硬い。セメンタイトの炭素量は、6.7重量%(鉄3原子に対して、炭素1原子の割合なので、鉄の原子量55.845と炭素の原子量12から12.0/(55.845*3+12.0)=0.0668)となる

まとめると、破断エネルギーや引張強度、伸び率は、壊れにくさ(折れにくさ、刃こぼれしにくさ)の指標となり、硬度や降伏強度は曲がりにくさ(塑性変形のしにくさ)を表す。ついでに、Young率は、弾性変形のしにくさの指標となる。



というわけで、調べれば良い量と物理的意味が何となく分かった。目安として、木材の引張強度は、数十〜150MPaくらいまで、ばらつきがあり、スギは90MPa、ケヤキは130MPaくらいらしい。伸び率は30%くらいにはなるよう。木材は、普通あまり塑性変形しないので、これも脆性破壊する物質。ガラスの引張強度も、30-90MPa程度らしく、木材と大差ない。ガラスは圧縮強度が900MPaくらいあるらしいけど。石材も傾向としては、ガラスと同じで引張強度は数十MPa、圧縮強度は1000MPa前後のことが多い。ヤング率に関しては、木材は、5-15GPa程度で、ガラスは70GPa程度、鉄系合金は大体200GPaとなる。引張強度は同じような値でも、ガラスや石の伸び率は、多分非常に小さい(1%以下だろう)と思うので、破断エネルギーで見れば、木材は、ガラスや石材より強靭である可能性が高い。

現在、工業的に使われるステンレス鋼も色々種類があるけど、包丁などに使われることのあるSUS440A(Amazonで適当に検索して出てきたステンレス包丁の素材だった)は、焼きなまし状態で、引張強度590MPa、0.2%耐力245MPa、伸び率15%、ビッカース硬度270らしい。


具体的に、青銅の物性を調べようとしたけど、あまり良いデータが見つからなかった。

Bronze#Transition_to_iron

https://en.wikipedia.org/wiki/Bronze#Transition_to_iron

には、"Though bronze is generally harder than wrought iron, with Vickers hardness of 60-258 vs. 30-80,[12] the Bronze Age gave way to the Iron Age after a serious disruption of the tin trade: the population migrations of around 1200-1100 BC reduced the shipping of tin around the Mediterranean and from Britain, limiting supplies and raising prices."と書いてある(適当訳:一般に、青銅は錬鉄より硬く、ビッカース硬度は60-258と30-80であるけども、錫の交易が著しく減退した後、青銅器時代は鉄器時代に移行した。紀元前1200-1100年頃の人口移動が、地中海及びイギリスからの錫の運送を滞らせ、供給の減少と価格の高騰を引き起こした)。とりあえず、青銅のビッカース硬度は、60〜258らしい。錬鉄は、炭素含有量の低い鉄と一般的に書かれている。明確な定義はないけど、炭素含有量0.05%以下とか0.1%以下となっていることが多く、厳密には鋼の一種(鋼は炭素含有量が0.0218%から0.214%のFe-C二元合金と定義されている)であるか純鉄であるかということになる。鉄鉱石を直接還元して出来る海綿鉄の鍛打で不純物を除去して作る方法と、銑鉄を脱炭して作る方法がある。Googleが教えてくれた所によると、錬鉄の引張強度は234-372MPa、降伏強度は159-221MPaらしい


青銅であれば錫含有量、鉄/鋼であれば炭素含有量は、性能に影響する。剣に使用された青銅や鉄の組成については、考古学的なデータから得るしかない

EPMA法による殷墟青銅器の分析と古代中国青銅器鋳造法の解明

https://www.jeol.co.jp/applications/detail/1054.html

に、古代中国の青銅器の銅・錫比の調査結果があり、錫含有量は、10〜20%程度らしい(Fig8)。Fig9に、組成ごとの物性値がある。出典の一つとして、"Metallography and microstructure of ancient and historic metals"という本が引かれていたので、Google booksで眺めてみたところ、121ページのAppendix G,Fig197に同じような図が見つかった。本の方は、グラフデータの内、一部しか数値の大きさが分からないけど、上の報告にあるグラフとは、一部傾向が異なる。本の方では、錫含有率が10%の時のブリネル硬度は100を超えているが、上の報告の方では、6〜70くらいのようになっている。伸び率は、本の方では、錫含有率が7〜8%くらいの時がピークとなっているけど、上の報告では、3%の時がピークとなっている。試験片によって、これくらいの違いが出るものなのか、どちらかが正しくないのか私には分からない

国際的な金属材料の規格の一つにUNSというのがあるらしく、それのC90200というのは、錫青銅らしい。錫と銅以外は、微量の元素しか含まないので、定義上の青銅に近い。錫含有率は、6〜8%くらい。引張強度262MPa、0.5%耐力110MPa、伸び率は30%、ブリネル硬度は70と書いてある

https://alloys.copper.org/alloy/C90200

C90700は、錫含有率が10-12%で、引張強度250-350MPa、0.5%耐力120-200MPa、伸び率10-20%、ブリネル硬度60-100。

https://alloys.copper.org/alloy/C90700

他にも、Tin Bronzeと分類されているものはあるけど、亜鉛含有率が5%くらいあったりするので、回避。

上の報告にあるグラフだと、6〜8%の錫含有率では、目視で、引張強度280MPa、ブリネル硬度60、伸び率は10%程度。伸び率の値が3倍くらい異なっている。上の報告のグラフで、錫含有率が、古代中国で武器に使われた15%くらいの場合、目視で、引張強度340MPa、ブリネル硬度90、伸び率2.5%。古代中国以外の地域では、青銅の錫含有率は10%前後のことが多かったようなので、C90200やC90700に近い数値だったと推測される。古代中国で、青銅の錫含有率が高い理由は謎

錫青銅は、現在は、工業的には殆ど使われてないらしく、金属材料のJIS規格には存在しない。代わりに、微量のリンを含むリン青銅というものの物性を見てみる。今の所、古代の青銅器の分析でリンが含まれていたという話は見ないけど。C5210は、Snが7-9%、Pが0.03-0.035%含まれるらしい。C5210には、何を表しているのか分からないけど、質というのがある(熱処理の仕方が違うんだろうけど)。質がHだと、引張強度は590-705Mpa,伸び率は20%以上、ビッカース硬度185-235、弾性限界390MPaらしい。何だかよく分からないけど、性能が随分強化されている。UNS規格では、錫含有量5%程度のリン青銅C51000や10%程度のC52400という材料がある。

https://alloys.copper.org/alloy/C52400

1/2Hardだと、引張強度572MPa、0.5%耐力は記載なし、伸び率32%、ロックウェルB硬度92(ビッカース硬度だと200くらい)という数値が書かれている。錫含有量5%のC51000の1/2Hardで、引張強度450MPa、0.2%耐力372MPaなので、C52400の0.2%耐力も400MPa以上はあるだろうと思われる。性能的に、C5210と大体同じ数値。リンは、青銅内部の酸化物を除去し、機械的性質を向上させるとか書かれてるので、むしろ、青銅本来の性質は、こっちなのかもしれない。


紀元前6世紀頃の中国では、白銑鉄が作られてたらしいけど、ねずみ鋳鉄があったのかどうかは分からない(多分、なかったことはないと思うけど)。炭素含有量の高い高温の固体鉄を、冷却する際に、除冷すると、セメンタイトがFeとCに分解して、組織中に遊離グラファイト(鉄の方は、低温ではフェライト相/組織になる)を生じるのに対して、急冷すると、セメンタイトFe3Cが、そのまま残る。前者が、ねずみ鋳鉄で、後者が白銑鉄らしい。定義は結構曖昧で、両者が混じったまだら鋳鉄というのもあるらしい。

ねずみ鋳鉄の方がデータが多いので、ねずみ鋳鉄を調べると、以下の1967年の古い論文によれば、歪み1%以下で脆性破壊するっぽい。引張強度は200MPa前後になるようなので、これは青銅より低いと思われる。また、青銅の伸び率は、文献値のばらつきが大きいとはいえ、錫含有量を相当増やさない限り、歪み1%で破断とはならなそうなので、多分、破断エネルギーも、ねずみ鋳鉄剣よりは、ブロンズソードの方が、ずっと大きいだろうと推測される

ネズミ鋳鉄における応力-ひずみ曲線と組織との関連について

https://www.jstage.jst.go.jp/article/imono/39/6/39_480/_article/-char/ja/

JIS規格では、ねずみ鋳鉄が、FC100〜FC350まである。3桁の数字は引張強度を表すようで、ブリネル硬度は201〜277とかのよう。以下のWikipedia情報によれば、ネズミ鋳鉄の引張強度は350MPa、伸び率0.5%、ブリネル硬度260。これは前述の情報に近い。やはり青銅と比べると、伸び率の小ささが目立つ

白銑鉄は引張強度170MPa、伸び率は0.1%以下(?)、ブリネル硬度450とあるので、武器として使える感じではなさそう。春秋時代の中国では、銑鉄の生産が開始してすぐに、可鍛鋳鉄(Malleable iron)もあったと言われていて、これは、Wikipedia情報によれば、引張強度360MPa、伸び率12%、降伏強度230MPa、ブリネル硬度130とある

Table of comparative qualities of cast irons

https://en.wikipedia.org/wiki/Cast_iron#Table_of_comparative_qualities_of_cast_irons

古代中国で可鍛鋳鉄があったと言っているのは、黒心可鍛鋳鉄じゃないかと思うけど、JIS規格では、FCMBで始まる一連の金属材料が相当する。何種類かあるけど、FCMB35では、引張強度350MPa、0.2%耐力200MPa、伸び率10%、ブリネル硬度100程度であるらしい。概ね、Wikipediaの数値と合っていて、錫含有率10~12%の青銅C90700(引張強度250-350MPa、0.5%耐力120-200MPa、伸び率10-20%、ブリネル硬度60-100)と性能的に大差なさそうな値となっている。春秋戦国時代に、可鍛鋳鉄で作った剣は出土してないらしいけど


紀元前に銑鉄を使っているとこは、中国以外の地域では、今の所、知られていないと思う。金属材料のJIS規格では、炭素含有量が0.3%前後の炭素鋼S30Cというのがある。物性値は熱処理によって変わるらしいけど、焼ならしを行った場合、引張強度470MPa、降伏強度285MPa、伸び率25%、ブリネル硬度150前後らしい。低炭素鋼として、SS400という材料もあり、引張強度は400MPa(以上)、降伏強度200-250MPa、伸び率35%前後で、炭素含有量の規定はないけど、0.15-0.2%くらいと言われる。これらの数値は、青銅よりは、優れていると言っていい程度のものと思う

最古の鋼片の検出とその意味

http://www2.pref.iwate.jp/~hp0910/tayori/106p2.pdf

というのを見ると、ヒッタイトの遺跡から出土した鉄片の一つは、推定炭素含有量が0.1〜0.3%だったろうとか書いてあるので、低炭素鋼の部類で、錬鉄を、そのまま使ってたわけではなさそう(?)。S30CやSS400に近い物性値を持っていたかもしれない。この鉄片が武器用を想定されたものかどうかも分からないけど。

Metallurgical Investigations on Two Sword Blades of 7th and 3rd Century B.C. Found in Central Italy

https://www.jstage.jst.go.jp/article/isijinternational/45/9/45_9_1358/_article/-char/ja/

には、紀元前7世紀の古代ローマ(王政ローマがあったとかいう伝説的な時期)の鉄剣に関する記述があり、5つの炭素含有量の異なる鉄から出来ているとある(ローマのpattern-welded swordは、こういうものらしい)。一番炭素含有量が高いところで、平均0.15-0.25%、低い方では、0.05-0.07%とある。多分、現在の炭素含有量を測っているようなのだけど、数千年経って炭素含有量が変わらないものなのか不明。炭素含有量が減ったりしてないのであれば、古代ローマの鉄も、炭素含有量は、0.3%を大きく超えない程度だったのかもしれない


ダマスカス鋼になると、降伏強度が1GPaを超えるものがあったよう(要出典。逸話を信用すれば、ダマスカスブレードは、しなやかで曲げても折れないと言われる)なので、ダマスカスブレードはブロンズソードよりも圧倒的に高性能だったかもしれない。ダマスカス鋼も、アレキサンダー大王に献上されたとかいう話や、19世紀にファラデーが研究したという話もあることから、多分、生産の歴史は2000年以上に及ぶので、ずっと同じ方法で作られてたかどうかは分からないけど。

隕鉄の場合。情報は、あまり多くないけど、以下の論文のintroductionで、ギベオン隕石から鍛造したロッドが、引張強度402MPa、伸び率5.6%と報告されているとある。当然、使用する隕石や鍛造の仕方によって値は変わるだろうけど、青銅と性能的に大差ないかもしれない。青銅の伸び率の文献値が幅があるので、何とも言えないけど、破断エネルギーで見れば、隕鉄は青銅より脆い可能性もある

The Yield Strength of Meteoritic Iron

http://adsabs.harvard.edu/full/1970Metic...5...63K


青銅より前は、銅の時代があったと言われている。銅斧は出土しているらしいけど、「どうのつるぎ」が実際にどれくらい作られたか分からない。純銅は、C1020(JIS規格)というのが、銅99.96%以上の無酸素銅というものらしい。質Oだと、引張強度250MPa、降伏強度200MPa、伸び率50%、ビッカース硬度50。質Hだと、引張強度350MPa、0.2%耐力350MPa、伸び率4〜5%、ビッカース硬度115らしいので、機械的性質にも相当に差がある。最初に使われた銅は、自然銅だと一般に考えられていて、純度は98%以上はあり、次いで多い成分が酸素というのが一般的なようなので、純銅というよりは、粗銅という感じではある。自然銅は、焼きなましをしてない素の状態では、展延性には乏しいらしいけど、物性値などは分からない。


あと、リン青銅剣を作れば、並の鉄剣と同等以上の性能になる可能性がある。Wikipediaの「りん青銅」の項目には、19世紀頃に、鉄で大規模な鋳造が難しかったので、大砲の鋳造に用いられていたとある。検索すると、1874年に、Phospher-Bronzeというタイトルの論文が出ていた。

Phosphor-bronze

https://www.sciencedirect.com/science/article/pii/0016003274903421

論文の著者のC.J.A.Dickの兄弟であるGeorge Alexander Dickという人物は、1874年に、イギリスで、Phosphor bronze companyという会社を立ち上げたらしい。Dick兄弟がリン青銅の発見者というわけではないようで、論文冒頭には、"The invention is due to the owners of the Belgian Nickel Works of Val-Benoit," とあり、個人名などはない。青銅砲の開発のために、色々な配合の青銅を試していた的なことが書いてある。

論文のTABLE IやIIには、pulling stress,elastic stress,ultimate stressという用語が出ているけど、それぞれ、どの概念に対応するか分からん(直訳では、引っ張り応力、弾性応力、破壊応力だけど)。単位は、lbf/in^2かと思う(多分)。TABLE Iは、"Results of ezperiments to ascertain the Tensile Strength and resitance to Torsion of various Wires."あるけど、lbf/in^2としてMPaに換算すると、焼き鈍した状態で、pulling stressが、

銅: 255MPa
真鍮/黄銅: 355MPa
charcoal iron(木炭で精錬した鉄っぽいけど詳細不明): 318MPa
coke iron(コークスで精錬した鉄っぽいけどetc): 422MPa
鋼: 514MPa
リン青銅No.1: 405MPa
リン青銅No.2: 445MPa
...

となる。リン青銅の錫含有量の記載は見当たらない。多分、文章的にも数値的にも引張強度を測定したんじゃないかと思うけど嘘かもしれない。

同様に、TABLE Iによると、焼き鈍した状態でのultimate extension in per ct.(破断する限界の歪み?)は、charcoal iron,coke iron,steelが28,17,10.9%なのに対して、リン青銅は、42.4%〜46.6%となっている。これは、銅の34.1%より高い。リン青銅も、錫青銅と同じように、錫含有率を上げると、伸び率を犠牲にして、引張強度と降伏強度を高めることができるはずなので、リン青銅で、上記の鋼と同等以上の性能を実現できた可能性は高い

ついでに、19世紀のヨーロッパで作られていた鉄の性能が、こんな程度のものだったと分かる。多分、coke ironとcharcoal ironは共に錬鉄だろうと思う。この時代には、橋やレールの鋼材として、錬鉄が一般的に使われていたらしい。coke ironは、パドル法で作った錬鉄だろうけど、charcoal ironは、コークス高炉による銑鉄製造+木炭精錬炉による脱炭なのか、木炭高炉による銑鉄製造+木炭精錬炉による脱炭なのか分からない。

人類がブロンズソードの真の力を引き出すことなく、鉄器時代に移行してしまったのは、残念

2017-11-05

Peter Woitという弦理論の批判で知られた人が、数年前に

Use the Moment Map, not Noether's Theorem

http://www.math.columbia.edu/~woit/wordpress/?p=7146

という記事を書いていた(私には要点がイマイチ分からなかったが)


ネーターの定理は、「保存量と対称性が対応する」という標語で知られる有名な定理で、普通、Lagrangianを使って定式化される。moment mapは、ネーターの定理のHamilton形式版というかsymplectic版といえる。数学では、Lagrangianそのものが使われない傾向にあるので、ネーターの定理を見かけることは、少ない。一方、物理学者が、moment mapを使っているのを見たことは、記憶にある限りではない。

物理学者は(数学者とは逆に)Lagrange形式を好んで使うこと、また、彼らが考えたい系の多くは可積分でなく、エネルギー、運動量、角運動量以外の"自明でない"保存量は、"some kind of charge"以外出てこないので、ネーターの定理で得られる保存量を重視する必要性が低い(そのような自明でない保存量は、Runge-Lenz vectorやKowalevski topで出てくるものがあるけど、教科書に載ってないことも多い)ことなどが理由として考えられる


moment mapを使うべき理由として、私が思いつくのは、以下のようなもの

(1)余接空間以外の相空間(≒一般のsymplectic多様体)でも適用可能(このような力学系は、特に数学では、しばしば扱われ、Lagrangianが不明である)

(2)個別の保存量や"対称性"ではなく、保存量や対称性の集合に存在する構造(つまり、Lie群やLie環)について記述する

(3)運動量写像は、"力学には依存しない"ことが定義から分かる

(4)運動量写像の量子化も、数学的に厳密に定義できる


#moment mapでは、時間並進対称性に対応する物理量がエネルギーだという事実が出ないように見えるけど、これは、moment mapの問題ではなく、相空間の取り方の問題だと思う。空間並進対称性に付随する物理量が運動量であることと、時間並進対称性に付随する物理量がエネルギーであることは、明らかに同じ構造を持っているので、通常の扱いでは、相空間が時間とエネルギーに対応する座標を持ってないことがトラブルの原因と思われる。

#時間・エネルギー座標を持つ相空間として、extended phase spaceというものが使われることがある。extended phase spaceは、ガリレイ群のmoment mapを定義しようと思えば必要だし、相対論を考えても不自然ではない。時間/エネルギー変数を含む正準変換を行いたいケースもある(以下の文献などで使われている)

Stackel系の全ての保存量を保つ離散化 (可積分系研究の新展開 : 連続・離散・超離散)

https://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/42747


(3)について。ネーターの定理は、(少なくとも、見かけ上は)Lagrangianに依存するが、Lagrangianは、力学系ごとに異なる。一方、運動量写像は、symplectic多様体と、多様体への(symplectic形式を保つ)群作用のみで決まり、Hamiltonianには依存していない。従って、Hamiltonianを取り替えても、相空間と群作用が同じであれば、自明に同じ運動量写像を得る。


群作用が、Hamiltonianを保つなら、HamiltonianとPoisson可換な保存量を得ることができるけど、運動量写像の"型"は

¥mu:M ¥to ¥mathfrak{g}^{*}

で(Mはsymplectic多様体でetc.)、保存量はM上の関数でないと困る。"型"だけから類推してもpullbackを考えればよいというのは想像に難くなく、これはcomeoment mapと呼ばれることがある。comoment mapの"型"は

¥mu^{*} : S(¥mathfrak{g}) ¥to C^{¥infty}(M)

となる。S(¥mathfrak{g})は、Lie環の対称代数であるけど、Poisson代数になり、また、S(¥mathfrak{g}) ¥subset C^{¥infty}(¥mathfrak{g}^{*})に注意する。comoment mapは、群作用と可換で、Poisson括弧を保つ環準同型。環準同型なので、一次の元¥mathfrak{g} ¥subset S(¥mathfrak{g})の行き先が決まれば、残りも全部決まる。例えば、symplectic形式を保つ回転対称性SO(3)があれば、M上の関数L1,L2,L3で、{L1,L2}=L3,{L2,L3}=L1,{L3,L1}=L2を満たすものが得られ、これは、so(3)の基底のcomoment mapによる像である


#moment mapの定義によっては、保存量に定数を足す不定性が残ることがある。角運動量の例だと、L1,L2,L3をL1+c1,L2+c2,L3+c3に置き換えても(c1,c2,c3は異なる定数)、それぞれから定まる時間発展は変化しない。このような場合、moment mapのpullbackはPoisson括弧を保たない。moment mapの定義で、comoment mapがPoisson準同型になることを要請する場合もある。このへんは、上の(2)と関係する話


(4)について。いわゆるquantum moment mapも、comoment mapと対比すると分かりやすい。quantum moment mapの概念や名前を誰が最初に明文化したのかは分からない(概念自体は、誰でも思いつくようなものなので、多くの人が同時期に知っていたに違いない)。特に難しい概念ではないけど、日本語の教科書で見かけたことはないので、一応文献をあげておく

Lectures on Calogero-Moser systems

https://arxiv.org/abs/math/0606233

の25ページDefinition4.1を参照

#私は、昔誰かが"量子化運動量写像"と訳していたのを見て、その呼称を使っていたけど、今検索したら、このサイトしか出なかった。そもそも、普通に訳すと、"量子運動量写像"の気がする(これは、一件あった)。

quantum moment mapの"型"は、comoment mapの型に於いて、対称代数を普遍展開環に変え、関数環を非可換環に変えたものとなる(quantum comoment mapと呼んだほうが適切でないかとも思うが、非可換の場合の空間概念がないので、量子の場合、comomentの方しか存在しない)

¥kappa : U(¥mathfrak{g}) ¥to A

これは環準同型で、群作用と可換になるという条件が付く。この場合も、一次の元(つまり、Lie環の元)の行き先だけ決まれば、残りは全部決まる。上のEtingofによるLecture Noteでは、大域的な変換ではなく、無限小変換を考えている。この定義だと、ネーターの定理は、ほぼ自明である。

#P,XをAの元として、ad(P)(x) = [P,x]と置くと、ad(P)は線形かつ非可換ライプニッツ則を満たす(ad(P)(QR) = Q(ad(P)(R)) + (ad(P)(Q))Rとなる)。 この条件を満たす元XをDer(A)としているが、X=ad(P)となるPがいつでも存在するとは限らないけど、存在すれば、正に保存量となる(逆に、保存量Pを与えれば、ad(P)によって、Der(A)の元が決まるので、非可換環Aへの無限小変換も自動的に定まることになる。従って、普遍展開環からの環準同型は、いつでも何らかの無限小変換に付随するquantum moment mapとなっている)。一方、群作用を微分して直接得られるのはDer(A)の元である


ここでの量子化の意味は、"filtered quantization"と呼ばれるものを考えるのがいい。一般に、フィルター環A(filtered quntizationではincreasing filterを考える)に対して、associated graded algebra(日本語訳は次数化環?)gr(A)を定義することができて、gr(A)は可換とは限らないけど、可換代数である時には、自然なPoisson構造が入る。このような時、Aはgr(A)のfiltered quantizationと呼ばれる(gr(A)のPoisson構造まで込みで一致する必要がある)。filtered quantizationは、多分、最近使われるようになった用語で、あまり浸透していないっぽい(検索すると、arXiv:1212.0914やarXiv:1704.05144で見つかる)。


名前が付いてなかったものの、かなり昔から、界隈の数学者は、この意味で(暗黙のうちに?漠然と?)量子化を捉えていたっぽい(上記EtingofのLecture Noteでも、filtered quantizationという用語は使われていないが、正にこの意味での量子化になっていることが分かる)。この用語に従えば、Lie環の普遍展開環は対称代数のfiltered quantizationであり、Weyl代数は、T^{*}¥mathbf{C}^n上の多項式関数環のfiltered quantizationである

filtered quantizationで得たフィルター環から完備なRees代数を構成すれば、変形量子化になる。1990年以降に見つかってきた新しい例も多く、例えば「有限W代数はSlodowy sliceの量子化である」という時も、filtered quantizationの意味での量子化と理解できる(arXiv:math/0105225)。また別の例としては、(ある種の?)3-Calabi-Yau代数などがあるらしい(よく知らない)

Noncommutative del Pezzo surfaces and Calabi-Yau algebras

https://arxiv.org/abs/0709.3593

Calabi-Yau algebras viewed as deformations of Poisson algebras

https://arxiv.org/abs/1107.4472

associated graded algebraは、可換にならない場合でも、例えばClifford代数からは外積代数が得られる。外積代数には、Poisson superalgebraの構造が入り、『Clifford代数は外積代数の量子化である』と言われるのも、filtered quantizationの一種として理解できる。


#量子化すると、非自明な中心拡大が生じることがある(ガリレイ代数の質量項や、Virasoro代数のcentral chargeなど)。これらは、filtered quantizationとは別の枠組みで理解されるべきようにも見えるけど、よく分からない


(3)に戻ると、運動量写像にHamiltonianは必要ないので、例えば、相空間がsymplectic等質空間であれば、変換がHamiltonianを不変に保つか否かに限らず、運動量写像を考えることができる。そして、量子論では、このような運動量写像の量子化によって、spectrum generating algebraが作られ、波動関数に作用する。このような観察から得るべき教訓は、Hamiltonianを不変にするような"力学的対称性"だけでなく、相空間の幾何構造を不変にする幾何学的対称性に付随する物理量にも目を向けたほうがいいということだと思う


私は、普通のmoment mapで考えると混乱するので、comoment mapしか使わない(computer algebraを利用しやすいように、物事を代数的に考えておきたいというのもある)けど、comoment mapは、保存量が直接得られる点とquantum moment mapとの自明な類似性の二点に於いて、moment mapより優れていると思う。なので、「ネーターの定理より、moment mapより、comoment mapを使おう」といえる(数学では、moment mapの像を見たいこともあるけど)

雑survey: 可積分性の必要条件と十分条件

Lax方程式は、可積分系で頻出する道具の一つで、現在知られている多くの(古典)可積分系で、等価なLax方程式が得られている。

Lax equation

https://ncatlab.org/nlab/show/Lax+equation

には"Lax equation is used in integrable systems; namely some systems are equivalent to the Lax equation."などと書いてある。

#非線形可積分系の差分化とその現状

http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0889-07.pdf

は、1994年のものではあるけど、可積分性の性質として、Lax pairがあげられている(著者に、広田良吾が入っているので、"権威付け"として)。


原理的には、(有限自由度の古典可積分系では)作用・角変数から、Lax pairを作ることができる

L = ¥left( ¥begin{matrix} I & 2I ¥theta ¥omega ¥¥ 0 & -I ¥end{matrix} ¥right)

M = ¥begin{pmatrix} 0 & ¥omega ¥¥ 0 & 0 ¥end{pmatrix}

に対して(Iとθが時間依存する変数で、ωは定数)

¥dot{L} = ¥[L,M¥]

¥dot{I} = 0 , ¥dot{¥theta} = ¥omega

は同値。一般には、作用・角変数に対して、上のLとMを対角的に並べたものを、Lax pairとすればいい。勿論、作用・角変数が具体的に分かっているなら、もう系は解けているといってよく、Lax方程式は必要ないので、こんな話は、実用上は何の意味もない(つまり、このLax方程式は役に立たないものである)けど、大雑把には「Liouville可積分系はLax方程式で書ける」と言ってもよさそうである。

この議論は

Hamiltonian structures and Lax equations

http://www.sciencedirect.com/science/article/pii/037026939091198K

にあるものと同じ("any hamiltonian system, which is integrable in the sense of Liouville, admits a Lax representation, at least locally at generic points in phase space")。


可積分でなくても、Lax方程式で書ける系は、無数に存在する。というのも、非退化な不変内積を持つLie環¥mathfrak{g}の双対空間(標準的なPoisson構造がある)上で定義されるハミルトン力学系は、常にLax形式で書くことができるけど、当然、その殆どは、可積分でないだろうと思われる。Lax方程式の偉いところは、機械的に、沢山の第一積分が得られる点にあるけど、この場合、得られる保存量は、Casimir関数であり、例えば、Lie環がgl(n)であれば、次元は$n^2$に対して、独立なCasimir関数はn個しかないので、一般に全く足りない。というわけで、これも役に立たないLax方程式である(※)

※)しかし、この見方は、役に立つLax方程式への第一歩でもある。例えば、非周期有限戸田格子を"標準的な"Lax形式で書いた時、Liouville可積分であるための第一積分は、Casimir関数で与えられる(実際、tr(L^n/n)の形で書け、n=2の時が通常のHamiltonian)。Casimir関数は、classical r-matrixを通じて定義されるPoisson括弧に関しても包合的であるが、他の関数とPoisson可換とは限らなくなる。そして、Casimir関数を行列Mにmapする関数も、classical r-matrixで書ける。このような構成は、AKS(Adler-Kostant-Symes)の定理という名前で知られる。非周期有限戸田格子の場合は、これで可積分になるために十分な数の第一積分を得られる。classical r-matrixとLax方程式の一般論について、以下に書いた

classical r-matrixとLax方程式

https://vertexoperator.github.io/2017/11/04/rmatrix.html



物理的に意味のある例として、2D/3D Euler(流体)方程式はLax方程式で書け、かつ非可積分である。

Lax pair formulation for the Euler equation

http://www.sciencedirect.com/science/article/pii/0375960190908093

A Lax pair for the 2D Euler equation,

https://arxiv.org/abs/math/0101214

Lax Pairs and Darboux Transformations for Euler Equations

https://arxiv.org/abs/math/0101214


Euler方程式は無限自由度であるけど、Zeitlinによって、2D Euler方程式を有限自由度に"truncate"した系が知られていて(よく使われる"正式名称"みたいなのはないようなので検索しづらい)、これは、やはりLax方程式で書けるが、十分な数の第一積分が出てこない例となっている(多分、非可積分だと思うけど、証明があるかは知らない)

Finite-mode analogs of 2D ideal hydrodynamics: Coadjoint orbits and local canonical structure

http://www.sciencedirect.com/science/article/pii/016727899190152Y


Euler-Arnold方程式(有限次元or無限次元"Lie群"上の左不変or右不変計量に関する測地流)でLax形式で書けるものは沢山あるはず(Lie群の余接空間のsymplectic reductionは余随伴軌道なので)で、Euler流体の方程式は、その例となっている。



[補足1]『非線形可積分系の差分化とその現状』には、他に可積分性の条件として"N-ソリトン解を持つこと"が挙げられている。これについては、(そもそも、ソリトン解の数学的な定義は何かということを脇に置いても)微妙なとこではある。少なくとも、1−ソリトン解を持つ非可積分系は存在する(e.g. double sine-Gordon方程式)。「任意のNについて、N-ソリトン解が存在するなら、可積分か」という問いについては、私は答えを知らない

Exact, multiple soliton solutions of the double sine Gordon equation

http://rspa.royalsocietypublishing.org/content/359/1699/479

では、高次元double sine-Gordon方程式を考えて、時空の次元qに対して、(2q-1)-ソリトン解まで存在すると書いてある


まぁ、"少なくとも1−ソリトン解を持つ"という意味では、"non-integrable soliton equations"という言葉を使えるので、「ソリトン方程式は可積分である」という言明は、慎重に使う必要があると思う。以下のようなレビューもあるし

Dynamics of classical solitons (in non-integrable systems)

http://www.sciencedirect.com/science/article/pii/0370157378900741



[補足2]『非線形可積分系の差分化とその現状』に挙げられていないが、多くの古典可積分系はbiHamilton系であることが知られている("2つ"あるのは、HamiltonianではなくPoisson括弧なので、biPoissonと呼ぶほうが適切な気がするけど、biHamiltonianという名前が広く使われている)。biHamilton系に於いても、Lax方程式と同様、十分多くの第一積分を得るための処方箋がある(recursion operatorによって、"低次"のHamiltonianから高次のHamiltonianを作れる。hamiltonianとして、recursion operatorのべき乗のトレースがしばしば選ばれる)。けど、biHamilton系は、いつでも完全可積分かという問いは、結論としては、正しくない

MathOverflowの以下の質問に対する回答が参考になる

Connection between bi-Hamiltonian systems and complete integrability

https://mathoverflow.net/questions/14740/connection-between-bi-hamiltonian-systems-and-complete-integrability


リンクが切れている論文があり、タイトルも不明だったので、探したものを貼っておく。

Lax方程式とbihamiltonian系の関係に関する、F. MagriとY. Kosmann-Schwarzbachの論文というのは、多分これ(Journal of Mathematical Physics 37, 6173 (1996))

Lax–Nijenhuis operators for integrable systems

http://aip.scitation.org/doi/abs/10.1063/1.531771

Magri, Morosi, Gelfand and Dorfmanをsummarizeしている、R.G.Smirnovの論文というのはこれ

Magri–Morosi–Gelfand–Dorfman's bi-Hamiltonian constructions in the action-angle variables

http://aip.scitation.org/doi/abs/10.1063/1.532221


あと、回答にあげられていないけど、以下の論文も参考になると思う

Canonical Forms for BiHamiltonian Systems

https://link.springer.com/chapter/10.1007/978-1-4612-0315-5_12


また、全ての完全可積分系がbiHamiltonianではないという話も、以下にある

Completely integrable bi-Hamiltonian systems

https://link.springer.com/article/10.1007/BF02219188

論文で挙げられている"可積分だけどbiHamiltonian構造を持たない"例は、MIC-Kepler問題として知られる系(普通のKepler系はbiHamiltonianだけど、摂動を入れると、biHamiltonianでなくなる)。


[補足3]上記のMathOverflowの記事を見ると、完全可積分性に必要な運動の積分を構成する3つのアプローチとして、Lax形式、biHamiltonian構造、そして、変数分離(separation of variables/SoV)が挙げられている。変数分離は、解析力学の教科書に書いてある古い方法で、19世紀に知られていた多くの可積分系は、この方法で解かれている(と思う)けど、職人芸的なものである。最近(20世紀後半)になって、変数分離の理解にも進展があった(私は、あまり理解してないけど)。


Hamilton力学系が、スペクトルパラメータ付きのLax方程式と等価な場合、Lax行列の固有ベクトルである"properly normalized" Baker-Akhiezer関数の極と対応する固有値が、分離座標を与えるというのが、Sklyaninのmagic recipe(とSklyanin自身が書いている)というもの。Sklyanin自身は、むしろ量子可積分系への応用を念頭に置いていたらしい

Separation of Variables. New Trends.

https://arxiv.org/abs/solv-int/9504001

Sklyaninは以下のように書いているので、"recipe"として完成していると言えるのかは不明;

Amazingly, it turns out to be true for a fairly large class of integrable models, though the fundamental reasons responsible for such effectiveness of the magic recipe:“Take the poles of the properly normalized Baker-Akhiezer function and the corresponding eigenvalues of the Lax operator and you obtain a SoV”, are still unclear. The key words in the above recipe are “the properly normalized”. The choice of the proper normalization ~α(u) of Ω(u) can be quite nontrivial (see below the discussion of the XYZ magnet) and for some integrable models the problem remains unsolved

この場合、スペクトル曲線は代数曲線なので、系を調べるのに代数幾何が役に立つことになる。古典可積分系の解が、しばしば楕円関数や超楕円関数で書ける理由の説明にもなる



変数分離にはbiHamilton系からのアプローチもある(bihamiltonian theory of SoV)。鍵になるのは、よいbiHamilton系では、Nijenhuisテンソル(※)が存在して、その固有値が、Darboux–Nijenhuis座標というものの半分を与えるという感じらしい。こっちの方は、微分幾何学的で、理論的な理解はより進んでいる模様。一方で、量子可積分系への"移植"は不明

※)recursion operatorと呼ばれることもある。このようなoperatorの存在は、KdV階層に於いて、Lenardという人によって認識されたのが最初らしい

Separation of variables for bi-Hamiltonian systems

https://arxiv.org/abs/nlin/0204029

About the separability of completely integrable quasi-bi-Hamiltonian systems with compact levels

http://www.sciencedirect.com/science/article/pii/S0926224508000296

Generalized Lenard chains and Separation of Variables

https://arxiv.org/abs/1205.6937


現状、与えられたHamilton力学系を"解く"のに役立つLax形式やbiHamiltonian構造を見つけるのは、今の所、職人芸に頼るしかないと思う(し、存在するかどうかも明らかではない)けど、Lax形式やbiHamiltonian構造が分かれば、Liouville可積分性が要求する第一積分を全て発見することは、しばしば可能となる。けど、十分な数の独立な第一積分があれば、Liouville-Arnoldの定理によって求積可能であると言っても、これは原理的なものでしかなくて、実際に解く上では役に立たない(例えば、第一積分の等位集合の位相構造を調べるのは、一般には相当難しい)。変数分離は、実際に系を解くという目的に、より近く、実際的な方法であると言える(解を既知の関数で具体的に書けることと、物理的に意味のある情報を得ることは、あまり関係ないことが多いけど)

結局、Lax方程式やbihamiltonian構造を見つけるところが職人芸なので、可積分系の求積は職人芸のままだと思うけど、Lax方程式で書ける可積分系とか、biHamiltonian構造を持った可積分系を、ある程度系統的に作る方法が発見されており(AKSの定理だとか、Frobenius多様体からbi-Hamiltonian hierarchyを作れるとか。そうして作られる系には面白い例がある)、そのへんが重要な進展と言えるのでないかと思う


[補足4]

Canonicity of Baecklund transformation: r-matrix approach. I

https://arxiv.org/abs/solv-int/9903016

Canonicity of Baecklund transformation: r-matrix approach. II

https://arxiv.org/abs/solv-int/9903017

線形代数と解析力学

有限次元ベクトル空間Vに対して、対称代数S(V ¥oplus V{*})と対称代数S(V ¥otimes V^{*})には、それぞれ自然なPoisson構造が入り、前者は解析力学で基本的なPoisson代数として現れる。後者は、行列上で定義された多項式関数の集合と同一視できるので、線形代数に於ける基本的な考察の対象とみなされる(例えば、行列式は、行列上定義された多項式関数である)。また、前者の量子化はWeyl代数、後者の量子化はgl(V)の普遍展開環となり、それぞれ解析学と(Lie環の)表現論で基本的な対象といえる(gl(n)は半単純でないので、Lie環の表現論では、やや人気がないけど)。こうして、基本的な4つの代数を得る


S(gl(n))は、gl(n)の双対空間上の多項式関数と同一視できる。S(gl(n))は、n^2個の変数x_{ij}で生成される複素係数の多項式環に、以下で定義されるPoisson括弧を定義したものと同型

¥{ x_{ij} , x_{kl} ¥} = ¥delta_{jk} x_{il} - ¥delta_{il} x_{kj}

当然、この定義は、行列単位の交換関係

 ¥[E_{ij} , E_{kl}¥] = ¥delta_{jk} E_{il} - ¥delta_{il} E_{kj}

から来ている。


簡単な計算によって

 ¥{x_{ij} , ¥sum_{k=1}^n x_{kk} ¥} = ¥sum_{k=1}^n( ¥delta_{jk} x_{ik} - ¥delta_{ki} x_{kj} ) = x_{ij} - x_{ij} = 0

が分かる。明らかに

tr(¥sum_{i,j} x_{ij} E_{ij}) = ¥sum_{k=1}^n x_{kk}

である。同様にして、多項式

det(¥sum_{i,j} x_{ij} E_{ij})

を考えることができるが、これは、

 ¥{x_{ij} , ¥sum_{a,b} x_{ab} E_{ab}¥} = 0

を満たす。


Poisson代数Aに対して、

C(A) = { x ∈ A | ∀y ∈ A,{x , y} = 0 }

で定義される集合C(A)を、Poisson centerと呼ぶ。Poisson centerは、明らかにPoisson可換なPoisson代数をなす(Poisson可換とは、任意の元同士のPoisson括弧が0になること)。なので、Poisson centerの生成元を全部決定するという問題を考えることができる。

まぁ、detとtrが含まれることから想像される通り、直接計算で、

tr( (¥sum x_{ab} E_{ab})^p )

は、S(gl(n))のPoisson centerに含まれることが分かる。これで生成されることを言うのは、

ad(f)(g) = {f,g}

で定義されるPoisson随伴作用と、GL(n)の随伴作用を微分したものを比較すると、Poisson centerとGL(n)不変式の全体が一致することが分かるので、上の形の元は、S(gl(n))のPoisson centerを生成することが分かる


Chevalleyの制限定理は、簡約Lie環に対して

S(¥mathfrak{g})^G ¥simeq S(¥mathfrak{h})^W

で、Harish-Chandra同型というのは、簡約Lie環に対して

Z(U(¥mathfrak{g})) ¥simeq S(¥mathfrak{h})^W

となる。従って明らかに

S(¥mathfrak{g})^G ¥simeq Z(U(¥mathfrak{g}))

が言える(この同型もHarish-Chandra同型と呼ばれていることがある。Harish-Chandraが最初にどれを主張したのかは知らない)。この同型は、Poisson centerと量子化した代数の中心とが環同型と解釈できる。


Harish-Chandra同型やChevalleyの制限定理の場合と違って、

S(¥mathfrak{g})^G ¥simeq Z(U(¥mathfrak{g}))

は、任意の有限次元Lie環に対して、正しいことが証明されていて、Duflo isomorphismという名前が付いている。



S(gl(n))のPoisson centerは不変式論的に見ても同じものを得ることができるので、S(gl(n))のPoisson構造を考える必要は必ずしもない。Liouville可積分性を思い出せば、Poisson可換な代数は重要であるので、中心に限らず、もっと大きなPoisson可換な部分Poisson代数を考えるのは自然に思える。S(gl(n))の場合、Poisson可換な極大部分Poisson代数として、Mishchenko-Fomenko subalgebraというものと、Gelfand-Tsetlin subalgebraのclassical versionとが知られている

#Tsetlinさんはロシア人っぽいけど、Cetlin,Tsetlin,Zetlin,Zeitlinなど、名前の表記ゆれが激しい。ここでは、Tsetlinを採用


Gelfand-Tsetlin subalgebraは、歴史的に、quantum version(つまり、U(gl(n))の可換部分環)が先に作られたっぽい。誰が最初に考えたのか分からないけど、arXiv:math/0503140には、Vershik & Kerov(1985)とStratila & Voiculescu(1975)が独立に考えたと書いてある(私は全く知らないけど、Voiculescuは自由確率論の開拓者で、Kerov-Vershikは漸近表現論を開拓し、自由確率論の応用先となっているそうな。多分、n→∞の時が興味の中心なのでないかと思うけど)。ロシア語論文と古い本なので、確認してみる気にはならなかった

Gelfand-Tsetlin代数の定義は簡単で、gl(n)⊃gl(n-1)⊃...⊃gl(1)なので、Z(U(gl(n)),Z(U(gl(n-1)),...で生成されるU(gl(n))の部分環というのが、定義。これが可換環になることは、すぐ分かる。so(n)⊃so(n-1)⊃...⊃so(1)でも同じことが出来るし、量子普遍展開環やYangianでも、同様の定義ができる。classical versionは、中心の代わりに、S(gl(n))のPoisson centerを取っていけばいい。

"Gelfand-Tsetlin"という名前は、Gelfand-Tsetlin基底というものに由来する(1950年代の仕事らしい)。gl(n)有限次元既約表現空間に、U(gl(n))のGelfand-Tsetlin algebraの作用が自然に定まり、同時固有関数を考えると、それがGelfand-Tsetlin基底(今の文脈で言うと、Gelfand-Tsetlin algebraがsimple spectrumを持つというのが、Gelfand-Tsetlinの示したこと?)。物理の人とかは、so(3)のスピンl表現の基底として、|l,m>と書くものを取るけど、あれも(so(3)⊃so(2)なので)Gelfand-Tsetlin基底。


#シンプレクティック代数sp(n)の場合、sp(n-1)に制限していっても、Gelfand-Tsetlin基底は得られない。これは、sp(n)の有限次元既約表現を制限した時に、multiplicity-freeに分解しないため。Wikipedia見たら、何かYangian使って、Molevがsp(n)の場合にGelfand-Tsetlin基底を拡張したと書いてある

https://en.wikipedia.org/wiki/Restricted_representation#Gelfand.E2.80.93Tsetlin_basis


#Gelfand-Tsetlin基底は、無限次元表現でも得られることはある(例えば、so(4,2)の極小表現や、so(3,1)の既約ユニタリ表現、3次元Euclid代数の既約ユニタリ表現など)。一般には、単一のHamiltonianが与えられた時、それと交換する保存量の全体は可換になるとは限らない(例えば、角運動量などは互いに可換でない)けど、Liouville可積分性では、その中の可換な部分代数のみが問題となる。このような代数で極大なものが複数ある場合もあるかと思う。第一積分の包合性以外に、独立性が問題になるけど、今は代数的に考えているので、理論的には、超越次数やKrull次元などで、第一積分の数が測られる


Gelfand-Tsetlin代数のclassical versionについて最初に考えたのは、Vinbergのよう。

On certain commutative subalgebras of a universal enveloping algebra

http://iopscience.iop.org/article/10.1070/IM1991v036n01ABEH001925

Mishchenko-Fomenko algebraも、ここで初めて定義されたっぽい。


Poisson可換なPoisson部分代数を考える動機として、Liouville可積分性をあげたので、可積分性の話。そのまま、Gelfand-Tsetlin可積分系と呼ばれる古典可積分系がある。元々は、GuilleminとSternbergが、GT代数の実Lie環版(つまり、u(n) ⊃ u(n-1) ⊃ ....で考える)に付随する可積分系を考えて、その後、2004年頃と割と最近になって、Kostant-Wallachによって、gl(n)版が考えられた(ので、後者は、Kostant-Wallach theoryとか呼ばれていることもある)

Gelfand-Zeitlin theory from the perspective of classical mechanics. I

https://arxiv.org/abs/math/0408342

Gelfand-Zeitlin theory from the perspective of classical mechanics II

https://arxiv.org/abs/math/0501387

まぁ、通常の可積分系とは、大分趣は違う。so(n)版のGelfand-Tsetlin系も作れそうだけど、知らないなと思ったら

The Gelfand-Zeitlin integrable system and its action on generic elements of gl(n) and so(n)

https://arxiv.org/abs/0811.0835

あたりでやられていた。

あと、

Linear algebra meets Lie algebra: the Kostant-Wallach theory

https://arxiv.org/abs/0809.1204

はsurveyじゃないと書いてあるけど、"線形代数のspecialist"向けに書かれているようで分かりやすい。GT代数を積分して得られる群作用は、Ritz値を保つ(疎行列の数値計算界隈の技法であるKrylov部分空間法で出てくる用語)ということで、線形代数と数値計算にも応用があるかもしれない(実際に役立つ何かがあるのかは不明。可積分系に興味がないという人への導入にはなると思う)


Mishchenko-Fomenko algebraの方。この代数の定義は、上のVinbergの論文でされたようだけど、その起源は、1970年代のManakov topという古典可積分系の研究に遡る。Manakov topは、Euler topの一般化で、Euler topが、S(so(3))で定義されるので、S(so(n))への一般化を考えようというもの。Poisson多様体のLiouville可積分性は謎いので、余随伴軌道に制限すると、Poisson centerは定数関数となって、第一積分として役に立たない。Hamiltonian自身が余随伴軌道上で定数となる場合は無視することにする

#Manakov topは、SO(n)上の測地流として定式化できる(多分、定式化自体はArnoldによるものじゃないかと思うので、Manakov topという名前が適切かどうかは不明)。計量は慣性モーメントによって決まると見ることができる。相空間は、SO(n)の余接空間で、自然なSO(n)作用に関する運動量写像を関数環の方で見ると、S(so(n))->Fun(T^{*}SO(n))という型になる。関数のクラスFunはPoisson代数の構造が入れば何でもいい。この写像は、Poisson括弧を保つ環準同型になっている。Hamiltonianとかは、S(so(n))の像に入っていて、S(so(n))を"相空間"と思うこともできるし、よく知られる通り、T^{*}SO(n)のsymplectic reductionとしても、余随伴軌道が出る。


最も簡単なso(3)の場合、genericな余随伴軌道の次元は2なので、あと一つ保存量があれば、Liouville可積分になる。これはHamiltonian自身があるのでOK。一般のso(n)の場合、Hamiltonian以外の保存量が必要となる。Manakovは、スペクトルパラメータ付きのLax方程式で、Manakov topを表し、必要な保存量を得たらしい。Mishchenko-Fomenkoは、その構成を、argument shift/shift of argument methodと現在、呼ばれている方法で捉え直した。これは、S(g)の元を、Lie環の双対空間上の(多項式)関数と思った時、Poisson centerに属する関数fに対して、固定された元Zと形式パラメータtを与えて、f(X + t Z)をtの冪で展開した係数のなす関数という定義(つまり、Z方向にテイラー展開した微係数)。Zごとに異なる代数が得られるけど、Poisson可換となる

EULER EQUATIONS ON FINITE-DIMENSIONAL LIE GROUPS

http://iopscience.iop.org/article/10.1070/IM1978v012n02ABEH001859/meta


Mishchenko-Fomenko代数の量子化の問題を提起したのは、上のVinbergの論文(ここでの量子化は、filtered qunatizationの意味。Vinbergは、そのような言葉は使ってないけど。量子化した方は、shift of argument subalgebraとかいう方が、通りが良いかもしれない)。これに対して、"generic"な部分を解決したのが、以下の論文と思う。Yangianあるいはtwsited YangianのBethe subalgebraというものからやってくるらしい

Bethe Subalgebras in Twisted Yangians

https://arxiv.org/abs/q-alg/9507003

上の論文で除外されていた特殊ケースについても、ある種の極限を取ることによって、極大な可換部分代数を対応付けることができるらしい。Gelfand-Tsetlin代数も、このようにしてshift of argument subalgebraの退化したケースとして出る。一般に、簡約Lie環に対して、shift of argument subalgebraのモジュライ空間を考えることができ、gl(n)の場合は、種数0でn+2個の点付きリーマン面のモジュライ空間(Deligne-Mumfordの安定曲線のモジュライ空間)に一致する。

Degeneration of Bethe subalgebras in the Yangian of gln

https://arxiv.org/abs/1703.04147

などを参照。gl(n)以外の場合も同様に、shift of argument subalgebraのモジュライ空間を考えることができるらしいのだけど、知る限り、文献はない


gl(n)のGelfand-Tsetlin代数は、gl(n)の有限次元既約表現上でsimple spectrumを持ち、同時固有関数はGelfand-Tsetlin基底となった。Gelfand-Tsetlin代数は、shift of argument代数の退化したケースなので、shifto of argument代数が同様にgl(n)の有限次元既約表現上でsimple spectrumを持つかどうかという問題を考えられる(答えは知らない)。yesであれば、Gelfand-Tsetlin基底の変形が得られる。


#Manakov topは、Reyman&Semenov-Tian-Shanskyでも系統的に得ることはできる。この方法では、Lax形式が分かる(Poisson可換な代数を与えるのは、より単純ではあるけども、それだけでは、対応するLax pairを得ることはできない)。詳述してる文献がないけど、

A new integrable case of the motion of the 4-dimensional rigid body

https://projecteuclid.org/euclid.cmp/1104115435

など



線形代数は、ベクトル空間と線形写像の理論ではあるけど、線形性に全てを押し付けて理解できない側面がある。例えば、行列式は線形な関数ではないし、相似変換の軌道は、複雑な幾何構造を持つ。対角化などは、相似変換による軌道の代表元を取る操作と理解されるので、あんまり線形な過程ではない。線形代数の線形性を超えた情報の多くは、S(gl(n))に含まれていると考えられる。S(gl(n))と関係の深い力学系は、非周期有限戸田格子やGelfand-Tsetlin系、gl(n)-Manakov topなどがある。これらは、物理的な重要性は、今の所、あまりないように思う(戸田格子とかは、一応物理的な動機があって調べられたもののようだし、他にもなんかあるのかもしれないけど)。そういう系であっても、全然別の観点からの有用性があるかもしれない。

#まあ、一応、以下のような話を念頭に置いている

The QR algorithm and scattering for the finite nonperiodic Toda lattice

http://www.sciencedirect.com/science/article/pii/0167278982900690


というわけで、線形代数と解析力学は、それぞれ、2つの基本的なPoisson代数を調べる分野だという見方もできるようになる

2017-08-12

最近の老化研究

2017年3月の『Cell』の論文

Targeted Apoptosis of Senescent Cells Restores Tissue Homeostasis in Response to Chemotoxicity and Aging

http://www.cell.com/fulltext/S0092-8674(17)30246-5

大雑把には、(マウスの)"老化細胞"(Senescent cells)でアポトーシスを誘導してやって除去すると、若返るよという話(例えば、Fig7(D)では、毛が明らかに増量している)。ここ数年、老化細胞の除去による"若返り"の報告はいくつかあるが、いずれも老化細胞の除去に遺伝子改変を必要とするものばかりだったのに対して、この論文の方法では、遺伝子改変を必要としない(マウスの実験であるから、ヒトなどでは、別途検討が必要ではあるが)

cf)細胞老化が生体機能に及ぼす影響とその機構に関する研究

http://www.ncgg.go.jp/ncgg-kenkyu/documents/25-18.pdf

http://www.ncgg.go.jp/research/news/20160810.html

肺組織で特異的に老化細胞を除去できるように、遺伝子改変を行ったマウスでの類似実験を行っている

cf)Naturally occurring p16Ink4a-positive cells shorten healthy lifespan

https://www.nature.com/nature/journal/v530/n7589/full/nature16932.html


背景。正常な細胞は、老化刺激(テロメア短小化、DNA損傷、酸化ストレスetc)を受けると、恒久的に細胞周期を停止した細胞老化と呼ばれる状態になることが知られている。こうしてできた老化細胞は、増殖を再開することはないと考えられているが、ただ機能を失って存在するだけでなく、炎症性サイトカイン、増殖因子、細胞外マトリックス分解酵素など、様々な分泌因子を分泌するSASPと呼ばれる現象を起こし、周辺組織に炎症状態を誘発する。SASPは、周囲の細胞で細胞老化や癌を誘導するなど生体にとって有害な作用をもたらすと現在考えられている。

cf)老化した細胞ががん化を促進する仕組みをハエで解明(プレスリリース)

http://research.kyoto-u.ac.jp/research/141011_1/


論文では、老化細胞でFOXO4の発現が亢進することを確認し(Fig1(D))、また、FOXO4を阻害すると、老化細胞の生存率が下がる(Fig(1(H))ことから、FOXO4が老化細胞のアポトーシスを抑制し、老化細胞の生存・維持に寄与していると推測し、FOXO4のp53結合部位を含む部分ペプチドのD-retro inverso (DRI)-isoform(FOXO4-DRI)を設計し、マウスに腹腔投与(i.p. injection)したという感じ(5mg/kgを一日一回x3回)。

#D-retro inverso-isoformというのは初耳だったのだけど、配列とアミノ酸のキラリティ(つまりD型アミノ酸から構成される)が逆になってるような異性体のことらしい。なんだかよく分からないけど、DRI-isoformはしばしば、元のペプチドと同じ活性を保持し、かつ、安定であるらしい(多分プロテアーゼなどで分解されにくい?)。今の場合は、p53とFOXO4の結合を競合阻害するのが目的だろうと思う(部分ペプチドなので、p53を抑制する機能はない)


投与を行った"老いたマウス"の年齢は104〜110週齢程度のよう(Fig7 captionなど)。20ヶ月齢のマウスが人間でいうと60歳相当程度らしいので、人間換算で30歳弱と若干若い感じがする。マウスを直接扱ったことがないので、何とも言えないけど、老化研究では、もう少し高齢のマウスを対象にすることが多いように思う。一応、Fig6(B)の腎臓組織のSA-β-GAL染色では、130week miceでも、結構、老化細胞が増えているらしきことが見て取れるし、最初の方のcfであげた

細胞老化が生体機能に及ぼす影響とその機構に関する研究

http://www.ncgg.go.jp/ncgg-kenkyu/documents/25-18.pdf

でも、12ヶ月齢のマウスが使われている(これは、肺組織に於ける老化細胞の除去実験)。


毒性について。FOXO4は、FOXO1とFOXO3に比べると、注目されることが少なく、非老化細胞ではあまり重要な役割を持っていないらしい。実際にFOXO4をノックアウトしたマウスは、目立った表現型を示さない("FOXO4 knockout mice do not show a striking phenotype")らしいので、FOXO4の阻害は、非老化細胞に対する影響は少ないと予想される。


[疑問]FOXO4の阻害が老化細胞のアポトーシスを引き起こすなら、FOXO4ノックアウトマウスでは、老化細胞の数が少なく(論文の結果が正しければ)見た目の老化も抑制されそうに思えるが、目立った表現型は示さないという報告とは矛盾する気もする。老化が遅いなどの表現型は注意して観察しなければ気付かない微妙なもの(致死とかに比べればstriking phenotypeでない)で単に見落とされているという可能性もあるけど


[補足]癌抑制遺伝子として有名な転写因子p53は、細胞周期を停止し細胞増殖の抑制を行うだけでなく、細胞老化を引き起こし、アポトーシスを誘導することも知られている。これらの機能は、独立しているわけでなく、細胞やDNAがダメージを受けた時に細胞周期を停止し修復を試みると同時に、ダメージが重度で(癌化しそう)であればアポトーシスを起こし、中度であれば細胞老化を引き起こすのだろうと思われる(回復可能な程度に軽度であれば、細胞周期を再開する)。というわけで細胞老化を誘導する薬剤を高濃度で与えると、アポトーシスを誘導するっぽい。培養細胞に於いて、p53の発現を亢進すれば、癌は減るが老化は進行し(寿命も縮むらしい)、一方、抑制すれば、老化は抑えられるが癌は増えるらしい(マウス個体でp53をノックアウトすると早期に死亡するので解析は難しいが、培養細胞での実験と矛盾しない結果が知られている)。このような癌化と細胞老化のトレードオフは、テロメラーゼの高発現/ノックアウトでも観察することができる(マウスはテロメアが長く、テロメラーゼノックアウトマウスで個体老化が観察できるようになるまで数世代かかる)らしい

cf)癌抑制遺伝子p53 による老化の制御

http://www.pharmacol.or.jp/fpj/topic/topic121-2-130.htm


#細胞周期の一時停止、細胞老化の誘導、アポトーシスの選択がどのようにされているのか詳細は不明。p53には、リン酸化部位が沢山あり、中でもSer46がリン酸化される(このリン酸化を行う酵素はHIPK2,ATK kinase,DYRK2など複数ある?)と、アポトーシス関連遺伝子(AIP1という遺伝子が知られる)の発現が誘導されるらしい。一方、Ser15のリン酸化(AMPKなどによる?)は細胞老化関連遺伝子の発現と関連するらしい。DNA損傷が起きた時、Ser46のリン酸化はSer15のリン酸化に比べると、少し後になってから起こるらしい

The p53 response to DNA damage.

https://www.ncbi.nlm.nih.gov/pubmed/15279792

Ser46 phosphorylation of p53 is not always sufficient to induce apoptosis: multiple mechanisms of regulation of p53-dependent apoptosis.

https://www.ncbi.nlm.nih.gov/pubmed/17584297

Palmdelphin, a novel target of p53 with Ser46 phosphorylation, controls cell death in response to DNA damage

http://www.nature.com/cddis/journal/v5/n5/full/cddis2014176a.html?foxtrotcallback=true


素朴な疑問として、何のために細胞老化などという機構があるのか。癌抑制のためだけであれば、とっととアポトーシスを起こすなりして除去されてもよさそうなものだけど、わざわざアポトーシスを抑制してまで生存を図る意味が分からない。最近では、創傷治癒に於いて細胞老化の果たす役割が報告されたりしていて、他にもまだ知られてない生物学的機能が存在する可能性はある

cf)Cellular senescence controls fibrosis in wound healing

https://www.ncbi.nlm.nih.gov/pubmed/20930261



それらの話とはまた別に組織幹細胞の減少と機能低下も個体老化の一因として挙げられることがある。細胞老化と幹細胞は、共に創傷治癒に寄与する、という共通点があるっぽい。7月のNatureに

Hypothalamic stem cells control ageing speed partly through exosomal miRNAs

http://www.nature.com/nature/journal/v548/n7665/full/nature23282.html

という論文が出ていた。これは、視床下部神経幹細胞(htNSC)の減少が老化のトリガーとなっており、特に、exosomal miRNAの放出量が減ってしまうことが主要因なのでないかとのこと(生物種はマウス)。

ちなみに、論文中では"mid-aged mice (11–16 months old)" "aged mice (22 months and older)"と書かれていて、やっぱCell論文の110wk miceとかは中年手前なんだな、という感じ。こっちの論文のFig1(a)などを見ると、16 months-old miceでも、視床下部幹細胞の減少が見て取れるし、高齢マウスでは殆ど消失しかかっている。Cell論文に比べると、視床下部幹細胞の注入("炎症"のせいか、単純に移植しても、すぐに死んでしまうので、ドミナント・ネガティブIκB-αを発現する細胞を作ったと書いてある)などにより、老化の進行が抑制されて寿命が延長してる(16ヶ月齢雄マウスに移植して4ヶ月後の生理機能(Fig3(c))や18ヶ月齢雄マウスに移植して寿命計測(Fig3(d))


この論文では、視床下部幹細胞に注目しているが、

Muscle-derived stem/progenitor cell dysfunction limits healthspan and lifespan in a murine progeria model

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3272577/

のような論文も以前に出ている(タイトルを直訳すると、『筋由来幹/前駆細胞の機能不全が、早老症モデルマウスの健康寿命と寿命を制限する』)。若いマウスのmuscle-derived stem/progenitor cellを移植することにより、早老症モデルマウスの寿命を大幅に延長できたようである(Fig4(a))。通常の老齢マウスに対して、同様の実験を行って寿命延長などが起こるかは書いていない。まぁ、視床下部幹細胞が特別というわけでもなさそうではある。また、移植後、比較的短期間で効果が出ていることから、移植した幹細胞そのものが分化しているわけではなく、何らかの分泌因子が重要らしい(が、分泌因子の正体は不明)というところで、終わっている。視床下部幹細胞の方も事情は似ていて、移植した細胞が直接分化するわけではなく、こちらは、何らかのmiRNAが重要らしいということになっている

物凄く単純化していうと、老化細胞除去は、"悪いもの"を減らすという考え方で、一方、若い幹細胞を移植するというのは、"良い物"を増やすという考え方なので、割と真逆ではある。もし後者が重要であるならば、単純に老化細胞を除去しても、機能低下した幹細胞が回復するとは限らず、高齢になってからの老化細胞除去は、効果が薄いか、手遅れである可能性もある。老化細胞除去が比較的若いマウスに対して行われている理由かもしれないけど、どうだろう。

[] 体積粘性係数の分子動力学計算

等方的な流体の場合、粘性係数には、shear viscosityとbulk viscosity(体積粘性係数)の2種類がある。水のような液体の場合、非圧縮性近似をするのが一般的なので、Navier-Stokes方程式で体積粘性係数がかかる項は0になる。気体の場合も、流速が音速よりずっと遅い場合には非圧縮近似はよい近似になっているとされている。固体の場合は、弾性が支配的なので、通常、粘性そのものが考慮されない。


Wikipedia情報によれば、体積粘性係数の値は多くの場合、あんまり知られてないらしい

Volume viscosity/Measurement

https://en.wikipedia.org/wiki/Volume_viscosity#Measurement

The volume viscosity of many fluids is inaccurately known, despite its fundamental role for fluid dynamics at high frequencies. The only values for the volume viscosity of simple Newtonian liquids known to us come from the old Litovitz and Davis review, see References. They report the volume viscosity of water at 15 °C is 3.09 centipoise.


そもそも、通常の流体運動の範疇では体積粘性係数項は消えてしまって観測量と関わってこない(というか、NS方程式から消えてるだけだけど)から、測定が難しいのは仕方ない。多分、体積粘性係数を考慮に入れるべき身近な現象は音の減衰だと思う。なので、体積粘性係数を測定する方法として、音を利用するのは自然に思える。

Bulk viscosity and compressibility measurement using acoustic spectroscopy

http://aip.scitation.org/doi/abs/10.1063/1.3095471?journalCode=jcp

という論文には"To date, measurement of ultrasound attenuation is the only known way of experimentally studying bulk viscosity."と書かれている(Section Bの冒頭)。この論文によれば、25度に於ける水の体積粘性係数は、2.4cPであるらしい(TABLE III)

医療分野では、流体ではない生体組織の弾性や粘性(の分布)を測定するエラストグラフィと総称される手法がいくつか存在し、知る限り、いずれの方法でも超音波を何らかの形で利用している。

Elastography

https://en.wikipedia.org/wiki/Elastography

elastographyという言葉は、名前の通り、弾性の計測を主要な目的としている。癌などの腫瘍組織は一般に周囲の組織より硬くなる(乳がんが、しこりとして検出できる理由でもある)ことから、病変検出に有用と考えられている。一方、音の伝播の方程式に粘性項を入れるのは、固体であっても流体と変わるところはないので、同じ方法で粘性係数を決めることも理屈上はできる。実際に、そのような研究は存在するものの、弾性に比べると粘性を決定する医学的意義はよく分からないので、広く行われているとは言いがたいっぽい。しかし、何はともあれ、固体・流体を問わず、実験的に体積粘性係数を知る方法は、超音波の伝播や減衰を調べる以外にないらしい



現在は21世紀なので、簡単な分子の場合は、実験によらず、分子動力学で計算するという道がありえる。原理的には、shear viscosityも体積粘性係数も、Green-Kubo公式で計算できるはずなので。shear viscosityの計算は色々やられていて、手法もいくつか存在する

[1]Determining the Bulk Viscosity of Rigid Water Models

http://pubs.acs.org/doi/pdf/10.1021/jp211952y

は、分子動力学で、shear viscosityとbulk viscosityを決める論文。水のモデルとしては、TIP4P/2005がよいらしい。これは、他の論文でも、精度高いと書いてある


というわけで、Gromacsを使って、上の論文と同じことができないか試した。Gromacsのバージョンは5.1.2(aptで入れたら、これになった)。マニュアルには、gmx viscosityコマンドで、上の論文の手法が実装されてるよって書いてあるけど、そんなコマンドは存在しない。というわけで、手動で頑張ってデータを解析するしかない

また、水のモデルは、デフォルトでは、TIP4P/2005が入ってない

GROMACS files for the TIP4P/2005 model

http://www.sklogwiki.org/SklogWiki/index.php/GROMACS_files_for_the_TIP4P/2005_model

にあるのが、TIP4P/2005モデルらしいけど、面倒なので、普通のTIP4Pを使うことにした。



論文[1]によれば、

At the beginning, a MD simulation was performed at constant temperature and pressure conditions (NPT ensemble) for a period of 2 ns, in order to determine the density of liquid water as it is predicted by each water model.

The average energy ⟨E⟩ of the system at that density was computed from a second constant volume and temperature (NVT ensemble) simulation of 2 ns.

For the calculation of the bulk and shear viscosity at each temperature, we have performed a set of 100 independent constant volume and energy (NVE ensemble) simulations. For each trajectory, initial atomic momenta were sampled from a Maxwell distribution at temperature T.

らしいので、NPTアンサンブルで2ns→NVTアンサンブルで2ns→独立なNVE計算を100回(10psの平衡化の後、300psの計算をやると書いてある)という順番で計算するらしい。NVE計算では、長い時間やると、enerygy driftによって、エネルギー保存が成立しなくなるので短時間の計算を沢山やるということのよう。平衡化する場合、NVT→NPTの順が一般的な気がするけど、今回は逆


論文では、水256分子で計算しているので、適当に250分子を用意した。まあ、大体、以下のような感じのコマンドを叩く。エネルギー最小化は論文には書いてないので、特にやらなくてもいいかもしれない。NVE計算は、面倒なので100回やらずに30回しかやってない。

#-- We will get 250 molecules
gmx solvate -cs tip4p -o water.gro -box 2.1 2.0 2.0 -p conf/topol.top

#-- Energy minimization
gmx grompp -f conf/em1.mdp -o em1 -c water.gro -p conf/topol.top
gmx mdrun -s em1.tpr -deffnm em1 -v 
gmx grompp -f conf/em2.mdp -o em2 -pp em2 -po em2 -c em1 -t em1 -p conf/topol.top
gmx mdrun -s em2.tpr -deffnm em2

##gmx energy -f em1.edr -o min-energy.xvg
##gmx energy -f em2.edr -o min2-energy.xvg 

#-- NPT ensemble (to determine the density of liquid water)
gmx grompp -f conf/npt.mdp -o nptin -c em2.gro -p conf/topol.top -maxwarn 1
gmx mdrun -s nptin.tpr -deffnm nptout

#-- check density ,pressure, energy, temparature
echo -e "7\n8\n10\n14\n15\n" | gmx energy -f nptout.edr -o nptout.xvg

#-- NVT ensemble (to determine the average energy)
gmx grompp -f conf/nvt.mdp -o nvtin -c nptout -t nptout -p conf/topol.top
gmx mdrun -s nvtin.tpr -deffnm nvtout

#-- check energy
echo -e "7\n9\n11\n" | gmx energy -f nvtout.edr -o nvtout.xvg

#-- NVE ensemble
for i in `seq 1 30`;do
  gmx grompp -f conf/nve.mdp -o nvein -c nvtout -t nvtout -p conf/topol.top
  gmx mdrun -s nvein.tpr -deffnm nveout_$i
  echo -e "7\n8\n10\n" | gmx energy -f nveout_$i.edr -o nveout_$i.xvg    #--energy,temperature,pressure
done

使用したmdpファイルとtopol.topは、以下に置いた

https://gist.github.com/vertexoperator/aebb9676509e5eba608dfff6ff74bfc5




NPT計算の結果

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Total Energy               -9772.19        4.7    154.428    -16.766  (kJ/mol)
Temperature                 298.147      0.062    11.6618   -0.23701  (K)
Pressure                   -4.63437          2    868.369   -8.01756  (bar)
Volume                      7.51647     0.0023   0.124593 -0.00247449  (nm^3)
Density                     995.296        0.3    16.3923   0.194343  (kg/m^3)

とかなった。密度が大体正しそうなので多分OK。温度は298.15K=25度を指定しているので、そのへんで動く。論文[1]のTable Iでは、25度でTIP4Pを使用した時、0.9953(g/cm^3)という値を得たと書いてあるので、妥当な数値っぽい


NVT計算の結果は

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Total Energy               -9755.22        3.2    134.505    -16.928  (kJ/mol)
Temperature                 298.143       0.11    11.1177   -0.10815  (K)
Pressure                   -386.763        9.7    777.601    18.4497  (bar)

まぁ、よさげ?

Gromacsには、tcafコマンドというのがあって、私が理解してない謎の手法で、shear viscosityを計算できるらしい。

gmx tcaf

http://manual.gromacs.org/programs/gmx-tcaf.html


それで、NVT計算の結果とかに対して、

gmx tcaf -f nvtout.trr -s nvtin.tpr

とかやると、visc_k.xvgに(一部ヘッダ省略)

@    title "Fits"
@    xaxis  label "k (nm\S-1\N)"
@    yaxis  label "\xh\f{} (10\S-3\N kg m\S-1\N s\S-1\N)"
@TYPE xy
@    s0 symbol 2
@    s0 symbol color 1
@    s0 linestyle 0
 3.113 0.51114
 3.269 0.493288
 3.269 0.48944
 4.514 0.386852
 4.514 0.387317
 4.514 0.398109
 4.514 0.389984
 4.623 0.392635
 4.623 0.385249
 5.573 0.329574
 5.573 0.330219
 5.573 0.326714
 5.573 0.325943
 6.226 0.296062
 6.538 0.28417
 6.538 0.284974

などと出力される(デフォルトでは、visc_k.xvgにkとetaのみの出力がされる)。マニュアルによると、これを適当にfittingすると、shear viscosityが計算できるらしい

以下のようなコードで、実際にfittingしてみた

import numpy as np
from scipy import stats

k,y=[],[]
for line in open("visc_k.xvg"):
    if line[0]=="@":continue
    if line[0]=="#":continue
    ls = line.strip().split()
    k.append(float(ls[0]))
    y.append(float(ls[1]))


x=np.array(k)**2
y=np.array(y)
r=stats.linregress(x,y)
print(r.intercept)

計算すると、例えば0.540536022685という値を得たので、0.54e-3(kg/(m s))という感じであるらしい。水の粘性係数は、25度で0.000894(Pa s)らしいので、2/3くらいの値になっているが、論文[1]の方法では、TIP4Pで計算した場合、25度で、(4.83±0.09)e-4 (Pa s)という値になったらしいし、まぁTIP4Pの粘性に関する精度はこんなものであるらしい。NPT,NVEの結果から、tcafで計算しても、同じような結果を得た。TIP4P/2005なら、もっとよい結果が得られるらしい


Green-Kubo公式で、体積粘性係数を計算するために

B_{ACF}(t) = ¥langle ¥delta P(t) ¥delta P(0) ¥rangle

を計算する。

After an equilibration period of 10 ps, a configuration with total energy equal to the ⟨E⟩ decided before was selected. For that configuration and for the next 300 ps, the instantaneous stress tensor elements were stored in the disk every 20 fs for further analysis.

とか書いてある(<E>はNVT計算の結果を使うらしい)けど、<E>に到達するのに、ものすごく時間がかかるっぽかったので無視した(なんか勘違いしてる?)


よくやるように、

B_{ACF}{(t) = ¥langle ¥lim_{T ¥to ¥infty} ¥frac{1}{T} ¥int_{0}^{T} dt’ ¥delta P(t+t’) ¥delta P(t’) ¥rangle

で計算する。これが得られたら、論文の式(4)でfittingを行う。次のpythonコードが実装例

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt


def func(t , C , tau_f , tau_s , omega , beta_f , beta_s):
   fast_term = (1-C)*np.exp(-np.power((t/tau_f) , beta_f))*np.cos(omega*t)
   slow_term = C*np.exp(-np.power(t/tau_s , beta_s))
   return (fast_term + slow_term)


BACF = []
for i in range(1,31):
   pressures = []
   for line in open("nveout_{}.xvg".format(i)):
        if line[0]=="@" or line[0]=="#":continue
        ls  = line.strip().split()
        assert(len(ls)==4)
        t,E,_,P = ls
        pressures.append( float(P) )
   Pbar = sum(pressures)/len(pressures)
   deltaP = np.array([P-Pbar for P in pressures])
   ACF = []  #--auto-correlation function
   for t in range(len(deltaP)):
      if t==0:
          sv = np.sum(deltaP*deltaP)/(len(deltaP))
      else:
          sv = np.sum(deltaP[t:]*deltaP[0:-t])/(len(deltaP) - t)
      ACF.append( sv )
   if len(BACF)==0:
      BACF = np.copy(ACF)
   else:
      BACF = BACF+ACF


dt = 0.02  #-- pico seconds
xdata = np.arange(len(BACF))*dt
ydata = BACF/30

popt , pcov = curve_fit(func , xdata , ydata/ydata[0])
plt.plot(xdata[:1000], func(xdata[:1000], *popt))
plt.show()

#--積分値(total relaxation time)
print( np.sum(func(xdata , *popt)) * dt)

あと、単位に注意する(粘性係数の単位をPa sとなるように合わせる)。上のdtはpsなので、10^{-12}(s)で、Gromacsが出力する圧力の単位は、bar(1bar = 10^5 Ps)。体積はNVT計算の平均値、温度も同様にして、計算する。上のコードで得られた諸々の数値をもとに、

(7.51647e-27/(1.38e-23 * 298.15))*(ydata[0]*1.0e10)*np.sum(func(xdata , *popt)) * dt*1.0e-12

で、求めたい体積粘性係数の値が出る。私の場合は、0.0018268112792585417という値を得た。論文のTABLE2によれば、298KでTIP4Pを使った場合、0.001173という値を得ている。オーダーは合っている。ずれは、アンサンブル数の少なさ、方法の若干の違いによるものと思われる。最初のほうで見た2.4cPという実測値と比較しても、悪くない。論文を信じれば、TIP4P/2005の方が精度はよいらしい。おしまい

Boltzmannの壺とEhrenfestの壺

最近読んだ論文の読書感想文コーナー

[1]Statistical mechanics of money

https://arxiv.org/abs/cond-mat/0001432


N個の壺とE個のボールがあって、ボールはどれかの壺に入っているとする(論文には、ボールではなく、moneyと書いてあるし、壺とは書いてないけど、確率論に於いて、壺とボールはポピュラーなモデルなので壺にする)。(少なくとも一個以上のボールが入っている)壺をランダムに一つ選んで、やはりランダムに選んだ別の壺に移すということを繰り返すと、壺に入っているボールの平均数は当然E/Nだけど、Nがある程度大きい時、壺に入っているボールの数xの従う確率分布は、C exp(-x/T)で近似できる(Tは平均がE/Nになるように決める。Cは規格化定数)。ボール数をエネルギーと思い、壺と壺の間のボールの交換は、粒子と粒子の間のエネルギー交換の単純化と思うと、C exp(-x/T)はBoltzmann-Gibbs分布と見なせるというようなことが論文に書いてある。


統計力学のボルツマンの原理の類似を信じると、ボールがi個入っている壺の個数をn_iとした時

N=¥sum_{i=1}^B n_i

E=¥sum_{j=1}^B j n_j

という条件下で、

W = ¥frac{N!}{n_1! ¥cdots n_c!}

を最大化するような[n_1,...n_E]が、平衡状態では実現される。この条件から、Nが大きい時は、

 n_j/N ¥propto e^{-¥beta j}

であることが示せるというのは統計力学の教科書に書いてある通りで、dynamicsの話がなければ、まぁそういうもんかという感じであるものの、[1]の論文のように時間発展の仕方が定まれば、定常分布で巨視的状態[n_1,...n_E]が実現される確率は(原理的には)厳密に計算できる量となる(NとEが有限であれば、ただの有限Markov連鎖であるので)。その場合でも、(熱力学極限では)平衡状態が教科書通りに実現されるのかは自明でない気がする。で、[1]の論文にはこれが正しいことの説明があるが、数学の論文でなく、厳密な証明にはなってないと思うし、厳密な証明が知られているのかどうかは不明


このモデルを最初に考えたのは誰かと思ったら

非結合的代数によるランダムな衝突モデルの表現

https://ismrepo.ism.ac.jp/?action=pages_view_main&active_action=repository_view_main_item_detail&item_id=32615&item_no=1&page_id=13&block_id=21

によると、Boltzmann自身が、ほぼ同じモデルを考えていたっぽい。ただ、Boltzmannの結論だと¥sqrt{x} exp(-x/T)が出るらしい。ミクロカノニカル分布を仮定しても、この形は出ない。これは、本文中の条件(1.2)を仮定したためだと思うのだけど、ボールを一個ずつ交換する単純なモデルでは、これが成立する根拠はない。代わりに

A_{¥kappa ¥lambda}^{kl} = A_{kl}^{¥kappa ¥lambda}

を仮定して、Boltzmannと同様の議論を行えば、C exp(-x/T)は出る。これは良さそうな気もするけど、そもそも、このような計算が厳密には正しくないので微妙。例えば、N=4,E=10として、4つの壺に入っているボールの数(微視的状態)を(n0,n1,n2,n3)で表すとして、(1,2,7,0)->(0,3,7,0)と(1,2,3,4)->(0,3,3,4)の2種類の遷移を考える。1つめと2つ目の壺では、いずれも(1,2)->(0,3)という遷移が起きているけども、(1,2,7,0)から(0,3,7,0)へ遷移する確率と(1,2,3,4)から(0,3,3,4)へ遷移する確率は等しくない。従って、上のAがきちんと決まらない。


#元の論文では、ボールの数をmoneyと解釈しているが、現実の取引のモデルとしての妥当性は疑問(例えば、このモデルでは一旦金持ちとなっても、比較的容易に没落する)



一応、簡単な数値実験。

import numpy as np
from scipy.special import gammaln


N = 4  #--number of urns
urns = np.ones(N , dtype=np.int)*3
E = sum(urns)
p = {}
Q = {}          #-- N matrix
Nstep = 1000000

for i in range(Nstep):
    n0 = np.random.randint(0 , len(urns))
    if urns[n0]==0:continue
    n1 = np.random.randint(0 , len(urns))
    if n0==n1:continue
    bef = (urns[n0],urns[n1])
    urns[n0]-=1
    urns[n1]+=1
    aft = (urns[n0],urns[n1])
    #-- 本当は十分定常になってからcountすべき
    Q[(bef,aft)] = Q.get((bef,aft),0) + 1
    macro_state = tuple([sum(urns==b) for b in range(E+1)])
    p[macro_state] = p.get(macro_state,0) + 1



vtot = sum(p.values())
ctot = np.exp(np.sum(np.log(np.arange(N,E+N))) - np.sum(np.log(np.arange(1,E+1))))
w = np.zeros(E+1)
for k,v in p.items():
   t = np.array(k)
   c0 = np.exp(np.sum(np.log(np.arange(1,N+1))) - np.sum(np.where(t>0 , gammaln(t+1) , 0)))
   #--macrostate , observed probablitity , expected probability by microcanonical ensemble
   print( k , v/vtot , c0/ctot)
   w += t*(v/vtot)


for (s0,s1),c in Q.items():
    print (s0,s1,Q.get(((s1[1],s1[0]),(s0[1],s0[0])),0),c)


#AlderとWainwrightが1957年に始めた分子動力学計算(hard-sphere MD)では、衝突が起きるまで直進運動し、衝突時に保存則を用いて衝突後の速度を計算するというものであったらしい。MaxwellやBoltzmannが仮定していた分子運動論のイメージに近いのだと思う。この系は、決定論的でエネルギーも連続的であるので、Boltzmannの離散的なモデルよりは、現実寄りのモデルではあると思うけど、衝突によるエネルギー交換という単純なルールは保たれている



『Boltzmannの壺』モデルで、ある特定の壺に入っているボールの数の時間変化を追うと、時間平均はE/Nで、沢山サンプルを集めると、分布は、やはりC exp(-x/T)で近似できるようになるっぽい(数値実験の結果)。勿論、ボールの数はEが最大なので、NやEが有限の時は、C exp(-x/T)という分布は近似にしかならない。ところで、最初の分布は、集団に関するもので、後者は時間に関するものなので、この2つが同分布で近似できるというのは、"エルゴード仮説"っぽい(※)。これもやはり、E/Nを一定にして、NとEを無限大に飛ばした極限を考えると、この2つの分布のずれは0に収束していくと予想される。これも自明ではないと思うけど、厳密な証明は考えてないし、誰かが考えているかどうかも知らない


※)普通のエルゴード仮説は、任意の物理量の相空間上の空間平均と、相空間上の時間発展を追って時間平均を取ると等しいというものだけど、分子動力学屋は、よく、一粒子の物理量の時間方向の分布と、集団に対する分布が一致するという意味で、エルゴード仮説という用語を使っている気がする。MD屋のエルゴード仮説は、通常のエルゴード仮説よりも、強い主張になってる気がするけど、数値計算では成立することが多々あるらしい。Boltzmannのモデルでも、MD屋のエルゴード仮説が(近似的に)成立しているっぽい


#特定の壺のボール数の定常分布は、理屈上は、ボール数がk->k+1,k->k-1に変化する確率を求めれば厳密に計算できる。が、この計算は意外と面倒そうである(注目している壺にk個のボールが入っている時、残りの壺のいくつかは空であるかもしれないというのが面倒の原因)。別の直感的な考え方としては、Nが十分大きければ、2つの壺のボール数にほぼ相関はないと考えられる(例えば、N=2の時は、2つの壺のボール数には強い逆相関がある。Nが増えるにつれて、この相関が弱くなるのは、それほど非直感的ではないと思う)。一方、壺は、どれも特別でないので、各壺のボール数の定常分布は、どれも同じ分布に従うと考えられる。なので、集団としての分布は、独立同分布からサンプリングしたものとほぼ同じになるはず。定常分布に於いて、各壺にボールが(x_1,x_2,...,x_c)個入ってる確率をp(x_1,x_2,...,x_c)とすると、この確率分布は、exchangeable(x_1,...x_cの任意の置換によって値が不変)だから、De Finettiの定理によって、Nを無限大に飛ばした時、独立同分布に従うかのように振る舞い、この直感は正当化できそうな気がする(うさんくさい)


#等重率仮定を使った、よくあるタイプの計算「ある壺に入っているボールの数がxである確率は、残りの壺にE-x個のボールの数を分配する仕方の数W(E-x)に比例する」。例として、N=3,E=6の時を考えると、W(n)=n+1なので、p(x)=(7-x)/28が等重率原理から導かれる結果(平均はE/N=2になっている)。一方、定常分布を、厳密な遷移確率から数値的に計算すると、[p(0),p(1),p(2),p(3),p(4),p(5),p(6)] = [ 0.19047619, 0.25396825, 0.20634921, 0.15873016, 0.11111111, 0.06349206, 0.01587302]となって(検算として平均が2.0であることを確認)、当たり前だけど等重率仮定は全然成立していない。まあ、この例は、NやEが小さすぎて、そもそもBoltzmann-Gibbs分布で全然近似されてないので、問題とは言えないけど(ちなみに数値的に実験した限りでは、あるNまでは、p(0)<p(1)になっていて、Nが大きくなると、p(0)>p(1)と逆転するようである)


論文[1]には、エネルギー交換ルールの時間反転対称性を破ると、Boltzmann-Gibbsでない分布が出ると書いてある。私が最初同じモデルと誤解した、類似のモデルとして、multi-urn Ehrenfest modelがある。こっちは、毎回、ボールをランダムに一つ選んで取り出し、ランダムに選んだ壺に入れるというモデル。Boltzmannのモデルと違って、ボールをランダムに選ぶので、ボールが沢山入っている壺からはボールが取り出される確率が高くなり、少ししかボールのない壺からは、ボールが取り出されにくくなる。結果として、Boltzmannのモデルと違って、極端にボールが増えることもなければ減ることもない。十分時間が経った後は、特定の壺に入っているボールの数の時間分布は、B/Nの周りで正規分布(?)するっぽい。状態の集合は、Boltzmannのモデルと同じで、状態遷移の確率がBoltzmannのモデルと異なる。


multi-urn Ehrenfest modelは、元々は、N=2の時に、粒子の拡散の単純化したモデルとして提案されたらしい(N=2の時は単にEhrenfest modelという)。適当な容器内に沢山粒子があって、容器を2つの領域に分割し、各粒子がどっちの領域にいるかという問題を考えると、Ehrenfest modelとの類似性はイメージできる。


『Boltzmannの壺』とmulti-urn Ehrenfest modelの2つのモデルの違いはコードで書くと割と明確。

import numpy as np
import matplotlib.pyplot as plt

model_type = 0
N = 60  #--number of urns
urns = np.ones(N , dtype=np.int)*20
dist = []

for i in range(50000000):
    if model_type==0:
       n0 = np.random.randint(0 , len(urns))
    else:   #--multi-urn Ehrenfest model
       n0 = np.random.choice(np.arange(0,len(urns)) , p=urns/sum(urns))
    if urns[n0]==0:continue
    n1 = np.random.randint(0 , len(urns))
    if n0==n1:continue
    urns[n0]-=1
    urns[n1]+=1
    dist.append( urns[N//2] )  #--適当な壺の時間発展を追う


cnts = np.bincount(dist)
print( np.sum(np.arange(0,len(cnts))*cnts)/np.sum(cnts) )   #--print average

plt.plot(np.arange(0 , len(cnts)) , cnts)
plt.show()

Ehrenfest modelは、カットオフ現象というものが起きることが証明されてる。カットオフ現象は、多分1980年代に、Diaconisという数学者とその周辺の人々が発見した話らしい。multi-urn Ehrenfest modelにしても、最初のモデルにしても、有限Markov連鎖として定式化できる。有限Markov連鎖では、定常分布に到達する(定常分布に収束するとしてだが)のにかかる時間(混合時間)というものを考えることができるけど、定常分布に一様に近づいていくのかというと、必ずしもそうではなく、ある時間で急に定常分布に近付くということが起きるというのが、Diaconisらの発見らしく、カットオフ現象と呼ばれている(分布間の近さは全変動距離というものを使って測る)。


尤も、日常的な感覚からすると、ボールは区別することなく、また直接目にするのは、それぞれの壺に入っているボール数なので、それは、予想される通り、概ね単調に平均に収束していく。分布そのものの近さを測るということに意味があるのかどうかは、よく分からない


カットオフ現象を示すことが分かっているモデルはいくつかあって、最初に示されたのは多分、リフルシャッフルの数学的モデルと思う

Shuflling cards and stopping times

http://statweb.stanford.edu/~cgates/PERSI/year.html#86

マルコフ連鎖と混合時間 ー カード・シャッフルの数理 ー

http://www.kurims.kyoto-u.ac.jp/~kenkyubu/kokai-koza/H23-kumagai.pdf

2つ目の文献には、「一般のマルコフ連鎖について、カット・オフ現象が起こるための満足のいく必要十分条件は知られていないと言ってよい」と書いてある(数年前の話だけど、多分状況は変わってなさそう?)

Ehrenfest modelでカットオフ現象を示したのも、Diaconisらしい。ちなみに、どっちの解析でも、有限群の表現論が使われている。multi-urnだと、どうなるのか知らない

Asymptotic analysis of a random walk on a hypercube with many dimensions

http://onlinelibrary.wiley.com/doi/10.1002/rsa.3240010105/abstract

Ehrenfestモデルでは、"エントロピー"(正確には、定常分布との相対エントロピー)が単調増加することが厳密に証明されている。これ自体は、もっと広いクラスの有限Markov連鎖で成立するので、特筆すべき性質でもないけど、これはこれで全変動距離とは違う尺度で、定常分布への近さを測っていると見ることもできる。尺度を変えると、カットオフ現象とかがどう見えるのか(あるいは見えないのか)は知らない。ランダムに生成した連結なグラフ上のsimple random walkで実験した限り、KL divergence(負の相対エントロピー)の対数は、かなりキレイに直線的に減少していく(つまり、KL divergenceそのものは指数関数的に減衰する)(数値実験の詳細は下記コード参照)。こういう振る舞いをする理由が十分説明されているのかは不明



import networkx as nx
import numpy as np
import matplotlib.pyplot as plt


G=nx.fast_gnp_random_graph(950,0.02)
P=sorted(nx.connected_components(G), key = len, reverse=True)[0]

#--最大連結成分だけ残す
for c  in G.nodes():
   if not c in P:
       G.remove_node(c)


#-- transition matrix for simple random walk on G
W=np.zeros( (len(G.nodes()) , len(G.nodes())) )
for n , src in enumerate(G.nodes()):
   adjs = G.adj[src].keys()
   cnt = len(adjs)
   for tgt in adjs:
     m = G.nodes().index(tgt)
     W[n,m] = 1.0/cnt


#-- 定常分布の計算
eval,evec = np.linalg.eig(W.T)
idx = np.argmin(np.abs(eval-np.ones(len(G.nodes()))))  #--1.0に最も近い固有値のindex
pi = evec.T[idx]
pi /= np.sum(pi)

#--初期状態(任意)
p = np.zeros( len(G.nodes()) )
p[0] = np.random.random()
p[-1] += 1.0 - p[0]

#--適当なステップ数追跡
maxstep = 50
entropies=[]
for istep in range(maxstep):
    entropy = np.sum(p*np.where(p>0, np.log(p/pi), 0)) #-- KL divergence
    entropies.append(entropy)
    p = np.dot(W.T , p)



#--片対数グラフ
plt.plot(np.arange(0 , len(entropies)) , np.log(entropies) )
plt.show()


2017-06-11

球面上の測地流は、古典力学系として見ると、球面の余接空間を相空間としている。その関数環はPoisson代数であり、Poisson括弧の具体的な計算も、Diracの方法によってDirac括弧として計算できる。このPoisson代数を眺めると、Euclidean Lie algebraの対称代数からの全射準同型があることに気づく。というわけで、量子化した代数は、Euclidean Lie algebraの普遍展開環の何らかの商であると考えるのが自然。このことに最初に気付いた人が誰なのかは知らないけど、

Fundamental algebra for quantum mechanics on S^D and gauge potentials

http://aip.scitation.org/doi/abs/10.1063/1.530099

には、そういう風に書いてある。


Euclidean Lie algebra $e(N+1)$の既約ユニタリ表現は、Wignerの方法で群の表現から作ることができて、¥mathbf{R}^{N+1}への$Spin(N+1)$作用への軌道は、一般に球面で、little groupは$Spin(N)$となる(原点の軌道は一点からなりlittele groupはSpin(N+1)となるので、この場合は誘導表現はSpin(N+1)の既約表現である)。なので、$e(N+1)$の無限次元既約ユニタリ表現は、球面の半径$r$と$Spin(N)$の既約ユニタリ表現(WeylのユニタリトリックによってLie環$so(N,C)$の有限次元既約表現と言っても同じ)でラベルされる。Spin(N)の自明表現に対応するのは、標準的なL^2(S^N)への作用である。


$e(3)$の場合は、自明でないcoadjoint orbitは全てT^{*}S^2と微分同相になる(O(a,b)=¥{(x,p)¥in ¥mathbf{R}^3 ¥times ¥mathbf{R}^3 | x^2=a , x¥cdot p=b¥}とすると、pをp-(b/a)xに移すことで、O(a,0)とO(a,b)の微分同相写像が得られる)。$e(3)^{*}$やT^{*}S^2上で定義される力学系は色々あると思うけど、(Liouvilleの意味での)可積分系としては、球面上の測地流、spherical pendulum、Lagrange top、Kowalesvki topなどがある。後者2つのintegrable topは$e(3)^{*}$上で定義されているが、coadjoint orbitに制限することで、T^{*}S^2上で定義されていると思うこともできる($e(3)^{*}$上のHamilton力学系は非正準ハミルトン力学系であり、coadjoint orbitはいわゆるCasimir不変量で特徴付けられるsymplectic leavesとなっている)。


参考)3次元Euclid代数の無限次元既約ユニタリ表現については、以下に要点をまとめた

http://formalgroup.tumblr.com/post/160390024555/3%E6%AC%A1%E5%85%83euclid%E4%BB%A3%E6%95%B0%E3%81%AE%E6%97%A2%E7%B4%84%E3%83%A6%E3%83%8B%E3%82%BF%E3%83%AA%E8%A1%A8%E7%8F%BE


$e(3)^{*}$上で定義されたHamilton力学系を量子化することを考えると、自然な波動関数の空間/状態空間は、$e(3)$の既約ユニタリー表現である。これは、上にも書いた通り、一つではなく、無数にある。物理の実際の問題では、波動関数の空間はどれか一つに決める必要があるが、今の場合、$e(3)$のユニタリ表現であり、それを既約分解すると(※)、既約ユニタリ表現を得る。どのような既約ユニタリ表現が出るかは実際の物理の問題設定に依存する。

(※)I型の局所コンパクト群の場合、任意のユニタリ表現は、直積分の形に一意に分解できることが知られている(らしい)

#Casimir作用素は、U(e(3))の中心の元を指す。一方、$e(3)^{*}$上の関数環は対称代数S(e(3))=gr(U(e(3)))であり、grによって、Casimir作用素をS(e(3))に持って行った像がCasimir不変量となる。一般には、Poisson代数の元の内、任意の要素とPoisson可換なものをCasimir不変量と呼ぶのだけども、Poisson代数がLie環の対称代数の場合は、Casimir不変量は、このようにして得られる。



数学としては、各既約ユニタリ表現ごとにHamiltonianのスペクトルを計算できれば、問題は解けたことになる。Kepler系や球面上の測地流の時は分岐則から答えがわかるのでほぼ自明な感じだったが、Kowalevski topやLagrange topはそこまで簡単ではない。こうした量子可積分系の解法が理解されてきたのは1980年代以降のことだと思う。例えば、Kowalevski topのLax行列は1980年代後半に複数の人によって発見されている。

Lax representation with a spectral parameter for the Kowalewski top and its generalizations

http://link.springer.com/article/10.1007/BF00403470

The Kowalewski Top 99 Years Later: A Lax Pair, Generalizations and Explicit Solutions

https://projecteuclid.org/euclid.cmp/1104178400

不思議なことに、Kowalevski topのLax行列は、so(3,2)(のループ代数)に値を取る(Lagrange topの場合は、so(3,1)で、一般化して、so(p,q)-topというものが得られる)。この2つのintegrable topはGASDYM(generalized anti self-dual Yang-Mills)方程式の次元簡約によっても得られるが、その時ゲージ群は、Kowalevski topの場合はSO(3,2)、Lagrange topの場合はSO(3,1)とする。このへんの理由は、私にはよく分からない。



やや風変わりなintegrable topとして、Goryachev-Chaplygin(GC) topがある。これは、$e(3)^{*}$全体では可積分系ではないが、これが可積分になるようなcoadjoint orbitが存在する(大体文献見ると、$r=1$に固定しているのだけど、$r$の値は0でなければ何でもよくて"スピン"が0になることのみが可積分性の条件の気がする。半径の値はスケーリングで規格化してしまえば、何でも同じということかもしれないけど、一応同値でない既約表現が出るので、気にはなる)。e(3)のスピン0の表現空間(完備化すると、L^2(S^2)になる)には、so(3,2)の極小表現あるいは、同じことであるがsp(4,R)の振動子表現(の片割れ)が実現できる。この事実は、quantum GC topを"解く"のに利用された(最終的に、2つの無限次三重対角行列の固有値を計算する問題に帰着する)

Exact solution of the Goryachev-Chaplygin problem in quantum mechanics

http://iopscience.iop.org/article/10.1088/0305-4470/15/6/015/meta

など。ただ、この論文の計算は、Lie環の元の平方根を取るなど、やや危うい操作がある

Goryachev-Chaplygin top and the inverse scattering method

http://link.springer.com/article/10.1007/BF02107243

では、この問題は回避されていて、ついでにL-operatorを与えるなど、系統的な感じにもなっているが、最終的に出る式は同じ。ところで、後者の論文のR-matrixの形とRLL関係式を見ると、Yangian Y(sl(2))との関係を想像させる(当時はまだYangianは知られてなかったはず)。detailをcheckしてないので正しいのかどうか保証できないけど、以下の論文は、まさにその問題を考えているように思う

Transfer Matrix and Yangian Related to Chaplygin-Goryachev Top

http://iopscience.iop.org/article/10.1088/0256-307X/13/9/001

Yangian, Truncated Yangian and Quantum Integrable Systems

https://arxiv.org/abs/cond-mat/9509081

上の論文では引用されてないけど、truncated Yangianは、1988年にCherednikが考えたものらしい(当時は"truncated"という言い方はされなかった模様)。T(u)を有限次の展開で打ち切るので"truncated"と言われているが、打ち切る次数の大きさをレベルと呼ぶ。その言い方に従えば、上で考えているのは、レベル3のtruncated Yangianということになると思う


また、2000年頃になって、ある種の有限W代数が、truncated Yangianと同型になるということが発見された(数学的な有限W代数の定義では、gl(Np)のサイズpのJordan blockがN個並んでるような冪零元eを取ってきて、W(gl(Np),e)を考えると、これがY(gl(N))のレベルpのtruncated Yangianと同型になる、ということらしい(多分...))

Yangian realisations from finite W algebras

https://arxiv.org/abs/hep-th/9803243

RTT presentation of finite W-algebras

https://arxiv.org/abs/math/0005111

Yangians and W-algebras

https://arxiv.org/abs/1305.4100


#Higgs algebraは、(gl(N)ではなく)sl(N)の場合であるが、上のタイプの有限W代数であり、Y(sl(2))の(レベル2の)truncated Yangianとして実現されるのだと思う(要確認)

#現在、様々なタイプの"量子群"が知られており、それらの関係を理解するのに、RLL/RTT関係式を使うのは見やすいように思う。"sl(2)シリーズ"の場合のみではあるけど、以下の論文は各種"量子群"の関係が書いてあるので、沢山あるなぁと思って、ゲンナリするのに丁度いい

Cladistics of Double Yangians and Elliptic Algebras

https://arxiv.org/abs/math/9906189

#上の論文に出てくるR-matrixは、全て差法性R(u,v)=R(u-v)を持つが、応用上は、そうでないR-matrixが出てくることもある(ShastryのR-matrixなど)。こうしたR-matrixでも、RLL関係式を通して、何らかの代数構造が得られるらしいけど、謎が多い


その後、一般化として、gl(N)の場合、一般の有限W代数がtruncted shifted Yangianというもので実現されるという話があり、

Shifted Yangians and finite W-algebras

https://arxiv.org/abs/math/0407012

Representations of shifted Yangians and finite W-algebras

https://arxiv.org/abs/math/0508003

有限W代数やtruncated Yangianの表現論は、まだわからないことも多い(私が理解してないことが多いというわけだけでもなく)ので、現状、GC topを解くのに直接役立つわけではないし、今後も役立つのかどうかわからない。まぁでも、動機として、可積分系のLaxペアは、ヒマな人が試行錯誤して勘で見つけてるようにしか私には見えないものがあるので理解できる形にしたいというのもあるし、これらの代数から新しい可積分系を得ることもあるかもしれないし


Kowalevski topやLagrange topの場合は、classical r-matrixの形は分かっている

Classical R-matrices for generalized so(p,q) tops

https://arxiv.org/abs/nlin/0401025

これらの系に於いて、古典論のレベルでは、通常のループ代数ではなく、twisted loop algebraが出てくるので、量子化すると、多分twsited Yangian(のquantum double?)になるのでないかと思っているけど、定かではない(つまり、今までの話を素朴に総合するならtruncated twisted Yangianのようなものを考えるべきなのでないかと思うけど....)。この手の対称対から得られる力学系としては、他にEuler-Manakov topなどがある


#Yangianは、(1)XXX模型や(一次元)Hubbard模型のような格子模型、(2)(massiveな)integrable field theoryなどの対称性として出てきて、色んな定義があるけど、有限自由度の古典可積分系との関係では、half-loop algebraの量子化として導入される(Lie bialgebraの量子化に関する一般論から出るらしい)。古典可積分系で、ループ代数が出るというのは、多分1970年代から知られていた(ループ変数は、Laxペアのスペクトルパラメータ)。有限自由度の可積分系でYangianのような巨大な代数がそのまま対称性として働くわけではなく、truncationが効いてくるのだろうと私は理解している

・half-loop algebra <-> Yangian

・loop-algebra <-> quantum double of Yangian

・twisted half-loop algebra <-> twisted Yangian

・twisted loop algebra <-> ???


#twisted YangianはOlshanskiiが考えたのが始まりで、その後、Mackayらによって一般のsymmetric pairに対するtwisted Yangianが定義された模様

Twisted Yangians and symmetric pairs

https://arxiv.org/abs/math/0305285

Twisted Yangians for symmetric pairs of types B, C, D

https://arxiv.org/abs/1407.5247

Representations of twisted Yangians of types B, C, D: I

https://arxiv.org/abs/1605.06733

俺たちは何回コインを投げるべきなのか

表と裏が同じ確率で出ることが期待されているコインがあり、コインにイカサマがないことを確認したいとする。なんやかんやあって、統計的手法で、不正の有無(つまり、コインの表が出る確率が50%か否か)に決着をつけることになった場合、コインを何回投げれば十分といえるか


統計的手法にも色々あるかもしれないけど、不正の有無を統計的にチェックしろと言われた場合、オーソドックスに使われるのは仮説検定と思う。表が出る回数は二項分布で表されるが、コインを投げる回数Nが十分大きい時、これは正規分布で近似できるようになって、実際に表が出た回数をx、有意水準は、例えば95%とすると

¥frac{|x - Np|}{¥sqrt{Np(1-p)}} > 1.96

であれば、コインにイカサマがないという仮説は棄却される(但し、p=1/2で、不正がない時、表が出る確率)。二項分布を正規分布で近似した時の分散σは

¥sigma^2 = Np(1-p)

で、おおまかに、イカサマがない時の期待値Npから2σ以上外れていればアウトという判定基準になっている(偏差値でいうと、70以上か30以下)。イカサマを疑われている側からすると、イカサマしてない場合でも、冤罪をかけられる確率が5%あるということなので、もっと厳しくしたいと思うかもしれないが、そこらへんは話し合いで同意が取れたとする。もし、表が出る確率が本当に寸分の狂いもなく1/2だったら、(Nが十分大きい時には)コインを投げる回数が100回だろうが5000兆回だろうが、冤罪発生率が5%であることは変わらない。どんなに公正なコインでも、たまたま偏った出目になる可能性は0でないので、残念ながら、これを0%にすることはできない


逆に、本当にイカサマがあって、表が出る確率が60%くらいあった場合、コインをN=5回しか投げないで、有意差がないという結論になっても、明らかに、Nが小さすぎるので、もっとNを増やす必要がある。こういうのは検出力(statistical power)が足りないと言われる。また、極端な話、表が出る確率が90%であれば、少し試してみれば不正があることを見抜けそうであるけど、表が出る確率が60%などの場合は、それより多くのNを必要とすると予想される。確率的なものだから、表が出る確率が0%でない限り、たまたま全部表が出て、イカサマがあるという結論になってしまう可能性はある。けど、例えば、N=10として、表が出る確率が90%あるコインを投げた時、全部表になる確率は35%くらいあるのに対して、表が出る確率が60%のコインを投げて全部表になる確率は0.6%くらいしかない。


数学ばかり見てないで、現実を見ると、コインは表裏対称でないことが多いし、そもそも、人間が肉眼で区別できないほど表裏対称だったら、コインの表が出たのか裏が出たのか判定できない。表面にマーキングしたり汚れがある場合など、人間にとっては微細に思える差であっても、表が出る確率が、1/2から僅かにずれるには十分と思われる。何にせよ、コインが表が出る確率が完全に1/2であることは多分期待できない。仮説検定は、差があるかないかしか調べないので、表が出る確率が50.0001%とかで、わずかに裏より表が出やすいだけという場合でも、Nが大きければ、有意差ありという結果になる可能性は非常に高くなりうる。実際、最初に書いた仮説検定の式を少し変形すると

 p - 1.96 ¥frac{¥sqrt{p(1-p)}}{¥sqrt{N}} < ¥frac{x}{N} < p + 1.96 ¥frac{¥sqrt{p(1-p)}}{¥sqrt{N}}

のようになって、表が出た割合x/Nが許される範囲は、Nが大きくなるに従って狭くなっていく。一方、Nが大きくなれば、x/Nは、大数の法則によって、表が出る真の確率に収束していく。結果として、表が出る確率が寸分の狂いもなく1/2でない限り、有意差ありという結論になる可能性は、Nが大きくなると、どんどん高まっていくことになる。しかし、表が出る確率が、実は50.0001%で、統計的に有意な差があると言われても、常識的な感覚からすると、意味のある差とは思えない


ということを考えると、不正を検出するために必要なNの大きさは、どのくらいの不公平まで許容範囲と考えるかに依存する。つまり、表が出る確率が50±ε(%)の範囲外であればよくないコインであるが、範囲内であれば問題ない水準と考えられるεの大きさを事前に決めておく必要がある。このεの大きさは、有意水準αと同じく客観的な決め方はない。


更に冤罪発生率を0%にできないのと同様、コインに不正がある場合でも、たまたま、表と裏が同じ回数ずつ出てしまう可能性は0にはならない。つまり、不正がある時に不正を100%検出することはできない。とはいえ、Nを大きく取れば、表が出る確率が真に1/2でない限り、このようなことが起きる確率は減っていくと考えられる。というわけで、不正がある時に、正しく不正が検出できる確率を、有意水準と同様に設定する必要がある。これが検出力の定義で、1-βで表されることが多い。というわけで、α、β、εの3つのパラメータを設定すると、不正検出のためにコインを投げる回数Nを決められそうな気がする。


もう少し定量的に評価することを考える。さっきの式に戻ると、有意水準α=0.05と設定して、コインをN回投げた時、表が出た回数xが

 p - 1.96 ¥frac{¥sqrt{p(1-p)}}{¥sqrt{N}} < ¥frac{x}{N} < p + 1.96 ¥frac{¥sqrt{p(1-p)}}{¥sqrt{N}}

を満たせば、コインは公平であると帰結される(前回同様、p=1/2)。

コインに不正があって、表が出る確率がp+εであったとする。このとき、コインが表が出る回数は、やはり二項分布で記述され、上の区間に入る確率を計算することができる。この確率は、表が出る確率がp+εであるにも関わらず、不正はないと結論づけてしまう確率とみなせる。この二項分布は、Nが十分大きい時、平均N(p+ε)、分散N(p+ε)(1-p-ε)の正規分布で記述されることに注意すれば、不正検出ミスの確率は簡単に計算できる。表が出る確率が、例えばp+2εだとか、p+2.4εだというケースもありえるけど、検出力が最も低くなるのは、本当の表が出る確率がp+εか、p-εの場合なので、この2点のみチェックすればいい。更に、p+εの場合と、p-εの場合で、結果は同じことになるので、片方だけ見ればいい。この理由から、以下、εは正としておく


自然言語で書くと状況は複雑だけど、数式で書くと、表が出る真の確率がp+εの時に

¥int_{p-u}^{p+u} ¥frac{1}{¥sqrt{2¥pi¥sigma^2}} e^{-(t-¥mu)^2/2¥sigma^2} dt < ¥beta

という条件を満たすということになる。但し、p=1/2で

u = 1.96 ¥frac{¥sqrt{p(1-p)}}{¥sqrt{N}}

¥mu = p+¥epsilon

¥sigma^2 = (p+¥epsilon)(1-p-¥epsilon)/N

とする(積分変数のtはコインを投げて表が出た"割合"を表すので、平均と分散を、それぞれN、N^2で割っていることに注意)。で、積分の変数変換をすると

¥frac{1}{¥sqrt{¥pi}} ¥int_{(-u-¥epsilon)/¥sqrt{2¥sigma^2}}^{(u-¥epsilon)/¥sqrt{2¥sigma^2}}e^{-x^2} dx

のようになる。erf関数

erf(x) = ¥frac{2}{¥sqrt{¥pi}} ¥int_{0}^{x} e^{-t^2} dt

を使うと、最終的に

¥frac{1}{2} ¥left( erf((u-¥epsilon)/¥sqrt{2¥sigma^2}) - erf((-u-¥epsilon)/¥sqrt{2¥sigma^2}) ¥right) < ¥beta

となる。とりあえず左辺を評価するためにεを決める必要がある。まぁ、ε=0.01としておく。


数値的に、いくつか計算してみる。pythonで以下のようなコードを実行すると

import math

for N in [10000,15000,20000,25000,30000,35000,40000,50000]:
   p,eps = 0.5,0.01
   u = 1.96/math.sqrt(N/p/(1-p))
   sigma2 = (p+eps)*(1-p-eps)/N
   print(N , 1.0-(math.erf((u-eps)/(math.sqrt(2*sigma2))) - math.erf((-u-eps)/(math.sqrt(2*sigma2))))*0.5)

以下のような結果を得た。

10000 0.5159939775494847
15000 0.6877923080943419
20000 0.8074680951250419
25000 0.8854187382561318
30000 0.9337611517760266
35000 0.9626265172738802
40000 0.9793451544602965
50000 0.9940083977523118

大体

N=20000回で、検出力0.8ちょっと

N=25000回で、検出力88.5くらい

N=50000回で、検出力0.994...

などとなる。つまり、20000回くらいやれば、コインが表が出る確率が51%の時、80%以上の確率で有意差がある(有意水準は5%)という検定結果を得られる計算。まぁ、数万回くらいはコインを投げることになりそう。一日が86400秒なことを考えると、実行するには、くそだるい数字である。


検出力は、0.8とか0.9にすることが多いらしいが、αやε同様、客観的に決める基準はないので慣習的なものに過ぎない。また、有意水準を厳しくすると、同じ検出力を得るには、より多くの試行回数を必要とする。例えば、3σを要求すると

10000 0.15860713527390746
15000 0.29094698846888134
20000 0.43187317605341713
25000 0.5644691796250857
30000 0.6787457862331097
35000 0.7708974857449786
40000 0.841393149895479
50000 0.9295476650393496

のようになった。εが厳しすぎる気もするので、5%くらいまで緩めると、数千回程度でも十分な検出力を期待できるようになる。大雑把に言って、有意区間の幅は、Nの平方根に反比例するので、(αやβを変えずに)εを1/10にしたければ、Nを100倍にする必要があり、逆にεを10倍にしてよければ、Nは1/100程度まで小さくできる。このおかげで、今の場合、(例えば、濡れ衣を着せてやろうと思って)Nを大きくして小さな差でも有意になるようにしようとしても、実際に実行するのはかなり難しくなる。例えば、εを0.001%にして(上で計算したεが1%の時と)同じ有意水準、検出力で検定しようと思えば、100万倍ものデータ量が必要となる。



元ネタ(検定力と検出力は同じもの):

【翻訳】ダメな統計学 (3) 検定力と検定力の足りない統計

http://id.fnshr.info/2014/12/17/stats-done-wrong-03/

これは本にもなってる。Webページには、コインの例が載っているのだけど、上に書いたような計算は載っていなかった(数式を載せるような本ではないせいだろうと思う)ので、計算してみた。適切な統計学の教科書には例題として載っていそうな気もするけど、私が大学で受けた授業では、検出力/検定力とか聞かなかった気がする。



演習問題:各面が同じ確率で出ることが期待されている(6面)サイコロに不正がないかを統計的に調べたいとすると、俺たちは一体、何回サイコロを投げるべきか?

方針:教科書的には、カイ二乗(適合度)検定を行うことが第一歩となるだろうと思う(他の検定手法もあるかもしれないが)。

次に、コインの場合のεに相当する量を何らかの形で定義する必要がある。これに正しい方法というのがあるわけではないけど、サイコロの各面の期待される出現確率をp_i(i=1,...6)として、実際の出現確率をp_i'とする時、次のような量を定義すると便利

w = ¥sqrt{ ¥sum_{i=1}^6 ¥frac{(p_i - p_i’)^2}{p_i} }

各面が同じ確率出ることを期待されているという条件は、勿論

p_i = ¥frac{1}{6}

を意味する。wは効果量(effect size)と呼ばれる。例えば、1〜3が出る確率が1/6+1/100で、4〜6が出る確率は1/6-1/100である(各面1%くらいの誤差)とすると、wは0.06となる。実際の研究では、wは、0.1,0.3,0.5などの値が選ばれることが多いらしい


効果量の定義の仕方は、検定の手法ごとに異なりうるし、検定の手法を決めても、必ずしも標準的なものが一つあるというわけでもないらしいので、Cohen's wという名前で呼ばれることもあるらしい

Cohen's w

https://en.wikipedia.org/wiki/Effect_size#Cohen.27s_w


それで、サイコロを投げる回数をNとして、wを使うと、カイ二乗値は、"平均的には"

¥chi^2 = ¥sum_{i=1}^6 ¥frac{(O_i - E_i)^2}{E_i} = N w^2

となり、Nに比例して大きくなることが分かる(w=0でも、期待値は自由度kに等しくなるので、本当の意味での平均ではない)。カイ二乗値は、帰無仮説が正しい時(つまり、p_i=p_i'が成立する時)には、近似的にカイ二乗分布に従う(今の場合、自由度は5)というのが、カイ二乗検定の要旨だったけど、p_iとp_i'がずれている時には、非心カイ二乗分布で近似されることは直感的に分かる。検出力は、この非心カイ二乗分布を積分して計算できるけど、非心度λを知る必要がある。ところで、非心カイ二乗分布の期待値が自由度kとλの和で書けることに注意すれば、λと、上で計算したカイ二乗値の値Nw^2とが関係することは意外ではない。etc. 具体的な計算は面倒なので省略

ノンパラメトリックベイズ(1)infinite GMM

元論文

[1]The Infinite Gaussian Mixture Model

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.9111

論文中に"The likelihood for α may be derived from eq. (12)"とあるけど、何度読んでも、式(15)が式(12)から出るというのは無理がある気がするし、そもそも式(15)の1つ目の式は多分正確ではない。


おそらく、式(16)のChinese Restaurant Process(CRP)から出すのが正しい

具体的には

p(c_1,¥cdots,c_n) = p(c_1)p(c_2|c_1)¥cdots p(c_n|c_1,¥cdots,c_{n-1})

で、j番目のクラスに割り当てられた要素の数を$n_j$として、そのindexを昇順にI(1,j),...,I(n_j,j)とする。

P_i = p(c_i|c_1,¥cdots,c_{i-1})

と略記すると、クラスの数を$k$として

p(c_1,¥cdots,c_n) = ¥prod_{j=1}^{k} ¥prod_{l=1}^{n_j} P_{I_{l,j}}

と積を並び替えられる。そして、CRPの定義(論文中の式(16))から

¥prod_{l=1}^{n_j} P_{I_{l,j}} = ¥frac{¥alpha (n_j-1)!}{(I_{1,j} - 1 + ¥alpha)(I_{2,j} - 1 + ¥alpha)¥cdots(I_{n_j,j} - 1 + ¥alpha)}

と計算できる。最終的に

p(c_1,¥cdots,c_n) = ¥frac{¥Gamma(¥alpha)}{¥Gamma(n+¥alpha)} ¥alpha^k ¥prod_{j=1}^k (n_j-1)!

という分布を得る。これは、Ewens分布と呼ばれるものになっている。

で、この分布に従って、n個の要素をk個のクラスに分割する確率は

p(k|¥alpha,n) = ¥frac{¥Gamma(¥alpha)}{¥Gamma(n+¥alpha)}S(n,k) ¥alpha^k

となる。但し、S(n,k)は第一種Stirling数で、Γ(n+x)/Γ(x)をxのn次多項式と見た時のk次の係数に等しい。これが、式(15)の1つ目の式の正しい形と思う。第一種Stirling数の定義から、

¥sum_{k=1}^n p(k|¥alpha,n) = 1

が容易に分かる。nとkを固定してしまえば、式(15)の2つ目の式は、正しいので、論文通り、実装しても問題は生じない


逆に、Ewens分布からCRPを導出するのも、条件付き確率分布の定義通り計算すればいいので容易。


Ewens分布は、exchangeableなランダム分割のモデルの一つで、他の例として、Ewens分布の1パラメータ変形であるEwens-Pitaman分布は、Pitman-Yor過程とかで使われるらしい。exchangeableという性質は、上のEwens分布で、c_1,...,c_nの順番を任意に置換しても結合確率が変わらないことを指している。exchangebilityは、今はGibbsサンプリングが使えるために要求される性質と思っておけばいいっぽい。

cf)Exchangeable random variables

https://en.wikipedia.org/wiki/Exchangeable_random_variables

原理的には、exchangebleなランダム分割のモデルに対してinfinite GMMができるのだと思う(というわけで、Pitman-Yor mixture modelというものもあるらしい)。infinite GMMに於いて、Ewens分布を選ぶ理由は特になさそうに見えるけど

Exchangeable Gibbs partitions and Stirling triangles

https://arxiv.org/abs/math/0412494

によれば、Ewens-Pitman分布は、exchangeableなランダム分割としては、かなり一般的なモデルとなっているらしい


#論文[1]と同じようにして、Pitman-Yor mixture modelを実装してみようかと思ったけど、hyperparameterの更新式が複雑になった。これをどう乗り越えるべきか確信がなかった(力技で計算する以外の方法が思いつかなかった)ので、今回はpending


[余談]ちゃんと読んだわけではないけど、以下の論文2.4.2節によれば、Ewens分布は、Jack測度というものと同一視できるらしい(Jack測度は、自然数nの分割上で定められた測度)。論文で扱われている測度は、Macdonald測度とでも呼ぶべきだろうか(一般的に通りのいい名前は付いていないっぽい?)。Schur多項式->Jack多項式->Macdonald多項式という(1パラメータ、2パラメータ)変形に対応して、Plancherel測度->Jack測度->"Macdonald測度"という変形があるっぽい。

A probabilistic interpretation of the Macdonald polynomials

https://arxiv.org/abs/1007.4779

[余談2]arXiv:hep-th/0306238によれば"In many aspects, the set of partitions equipped with the Plancherel measure is the proper discretization of the Gaussian Unitary Ensemble (GUE) of random matrices"だそうである。Jack多項式のαは任意パラメータであるが、α=1,2,1/2の時は、ある対称空間上の球関数となる。一方、ランダム行列側では、ガウス型アンサンブルとして、GUE,GOE,GSEが知られている。α=1のJack測度がPlancherel測度で、その適当な極限は、GUEランダム行列の固有値の極限分布に一致するらしい(arXiv:math/9903176,arXiv:1402.4615)。というわけで、Ewens分布は数学的に見ても、非常に素性のよさそうな分布である


#ランダム行列について何も知らないけど、"hermitian random matrix ensembles"については、ある種の対称空間のCartan classとの対応があって、GUE/GOE/GSEは、その一例であるらしい(arXiv:cond-mat/0304363,arXiv:0707.0418)。今やランダム行列ユーザーの最大勢力は物理屋だと思うけど、これらのrandom matrix ensemblesは、class Bとclass DIII-oddを除けば、トポロジカル絶縁体・超伝導体の分類にも現れる(arXiv:0905.2029)。class Bとclass DIIIが物理で現れないというわけでもないっぽい(arXiv:cond-mat/0103089,arXiv:cond-mat/0103137)けど、ハミられている理由はよく分からない。最近は、"non-hermitian random matrices"も、応用があるそうな。


もうひとつ、論文[1]の式(17)の後半で、新しいクラスを割り当てる確率の計算で、積分が必要になるけど、

[2]Markov chain sampling methods for Dirichlet process mixture models

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.115.3668

のAlgorithm8を使って積分を回避する。この論文には"However, the equilibirum distribution of the Markov chain defined by Algorithm 8 is exactly correct for any value of m"とあるのだけど、何故か分からない(mが大きければ、積分のモンテカルロ評価になるというのは論文に書いてある通りで明らかだけど)。このあとのステップで、同じ基底分布G0からφ_cをサンプリングすることになるので、そうなりそうな気もするけど...


全体の流れとしては、初期化が終わったら、以下の処理を繰り返す感じになる

(1)クラス割当の更新。論文[2]のalgorithm8を使う

(2)各クラスのcomponent parameters(今の場合、個々のクラスは、単一のガウス分布に対応しており、その平均と分散)をサンプリング。論文[1]の式(4)と(8)を使う

(3)hyperparametersの更新。論文[1]の適当な分布に従ってサンプリング。更新するhyperparameterは大きく分けると、基底分布G0のパラメータ(4つある)とDirichlet過程の集中度αの二種類となる

これを十分に繰り返せば、クラス割当/component parameters/hyperparametersの分布が得られるが、一番尤もらしい割当とパラメータが欲しいだけなので、それぞれの尤度を計算して一番いいものを選ぶ。十分多く繰り返せば、MAP推定値に近づくはず。


尤度は、DPGMMで対象データが生成される確率を計算する。ハイパーパラメータを決めれば、DPGMMに従うデータ点を生成していけるが、そのようにして実際に生成される点は、元の学習データとは、あまり似てないものになる公算が高そうな気がするけど、大丈夫なんだろうか。各クラスターの平均(と分散)は、基底分布G0に従って生成するので、特にクラスター数が少ない時には、元の学習データにあるクラスターとずれたところにクラスターを生成しても不思議ではなiい。そのへんは、通常の混合ガウス分布の方が、ずっと有能そうではある


#infinite GMMやDPGMMはノンパラメトリックベイズ的手法の一つと言うけど、ノンパラメトリックという言葉は全く正しくないように思うし、誤解と混乱を招く言い方の気がする。infinite GMMで使われているDPGMMのパラメータ(上でhyperparameterと書いてるもの)はαと基底分布のパラメータからなる。データの次元が一次元の場合は、基底分布のパラメータは4つ(正規分布の分が2つとガンマ分布の分が2つ)なので、合計5つのパラメータを与えれば、モデルを指定できる。一方、(一次元の)混合ガウス分布の場合は混合数をkとすると、平均と分散がk組、混合比がk-1個のパラメータで指定されるので、3k-1個のパラメータで定まる。一次元では、DPGMMは、混合数2の混合ガウス分布と同じ数のパラメータを持つ。


現在では、MCMC以外の、より高速な推定手法も提案されてはいるようだけど、特に調べてはいない。


問題として、βの分布が、k=1の時、log-concaveではない気がする(ので、ARSが失敗する場合がある)。以下のコードで、対数確率密度関数のplotをしてみたけど、k=1の時は、xが2か3あたりを超えると、微妙に凹でない

import matplotlib.pyplot as plt
from scipy.special import gammaln

s,w=np.array([ 0.69756949 ]) , 1.159287098021641
k = len(s)
x = np.arange(1 , 100)*0.1
y = -k*gammaln(x/2) - 1/(2*x) +((k*x-3)/2)*np.log(x/2) + (x/2)*np.sum(np.log(s*w)) - (x*w/2)*np.sum(s)
plt.plot(x, y)
plt.show()

f:id:m-a-o:20170611142103p:image

αも同様に、log-concaveでないと思う(プロットして見ても、かなり分かりにくいが、kが小さい時の方が、"凹度"は高い気がする)

import matplotlib.pyplot as plt
from scipy.special import gammaln

k,n=1,500
x = np.arange(1,101)*0.5
y = (k-3/2)*np.log(x) - 1/(2*x) + gammaln(x) - gammaln(n+x)
plt.plot(x, y)
plt.show()

(私が間違ってるのでないとすると)結構イイカゲンな論文だな、という印象なんだけど、どうなんだろう


以下に実装したコードを置いておく。面倒なので、データが一次元空間上に分布している場合のみを考えている。αとβのサンプリングは、論文に従ってARSで行っているものの、上で書いた理由で、正しくはない(し、エラーを起こす場合がある)。今回の目的はサンプリングではなく、MAP推定なので、αとβのサンプリングが多少悪くても、よかろうと思って、そのままにしてある


実際に、DPGMMでデータを生成して、infinite GMMでfittingして見ると、データが一次元のせいもあって、クラス間の分離が悪い場合は、本来存在するより少ないクラスしか出てこない場合があるように思える。対数尤度を計算してみると、infinite GMMで計算した方が高くなっているので、実装が間違ってるわけでもないと思う


まぁ、さしあたって、infinite GMMは、ノンパラメトリックベイズの一番簡単な例題として勉強しただけで、本来の目的ではないので、とりあえず、諸々の問題点・課題(α、βのサンプリング、高次元データ対応etc.)は、そのままにしておいて、次に進む


"""
This is an implementation of the following paper written in python 3.5
(The only univariate case is supported)

The In&#64257;nite Gaussian Mixture Model
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.9111

"""
import numpy as np
import math
from scipy.special import digamma


def rnd_from_piecewise_exponential(coeffs , intercepts , points):
    """
       density function h(x) in [points[i] , points[i+1]]
       h(x) = exp(coeffs[i]*x + intercepts[i])
    """
    assert(len(coeffs)==len(intercepts)),(coeffs,intercepts)
    assert(len(coeffs)+1==len(points)),(coeffs,points)
    u = np.random.random()
    CDF = [0.0]  #--cumulative density value at each points
    for i in range(len(points)-1):
       z_cur = points[i]
       z_next = points[i+1]
       assert(z_cur<z_next),(points)
       a = coeffs[i]
       b = intercepts[i]
       if a!=0.0:
           r = np.exp(b)*(np.exp(a*z_next) - np.exp(a*z_cur))/a
       else:
           r = (znext-z)*np.exp(b)
       assert(not np.isinf(r)) , (coeffs , intercepts , points)
       r0 = CDF[-1]
       CDF.append( r + r0 )
    u = CDF[-1]*u    #--CDFが正規化されてないので
    for i in range(len(points)-1):
        if CDF[i] < u and u <=CDF[i+1]:
             a = coeffs[i]
             b = intercepts[i]
             z = points[i]
             """
                compute t such that
                ( exp(a*t+b) - exp(a*z+b) )/a == u - CDF[i]
             """
             if a!=0.0:
                 t0 = a*(u - CDF[i]) + np.exp(a*z+b)
                 assert(t0 > 0.0),t0
                 t = ( np.log(t0) - b )/a
             else:
                 t = (u - CDF[i])/np.exp(b) + z
             assert(t >= z and t <= points[i+1]),(t,z,points[i+1])
             return t
    assert(False),("rnd_from_piecewise_exponential",u,CDF,points,coeffs,intercepts)


def ars(h , hprime , points , support=(0 , np.inf)):
     xs = [x for x in points]
     xs.sort()
     assert(len(xs)>0)
     assert(xs[-1]>0)
     while hprime(xs[-1])>=0.0:
         xs.append( 2*xs[-1] )
     while 1:
         xs.sort()
         x0 = xs[0]
         xN = xs[-1]
         zs = [support[0]]
         coeffs = []
         intercepts = []
         for i in range(len(xs)-1):
             x = xs[i]
             xnext = xs[i+1]
             assert(xnext>x)
             zi = (h(xnext) - h(x) - xnext*hprime(xnext)+x*hprime(x))/(hprime(x) - hprime(xnext))
             assert(zi>zs[-1]) , (zi,zs,xs , xnext,h(x) , h(xnext) , hprime(x) , hprime(xnext))
             zs.append( zi )
             coeffs.append( hprime(x) )
             intercepts.append( h(x) - hprime(x)*x )
         zs.append( support[1] )
         coeffs.append( hprime(xN) )
         intercepts.append( h(xN) - hprime(xN)*xN )
         y = rnd_from_piecewise_exponential(coeffs , intercepts , zs)
         w = np.random.random()
         ukval = 0.0
         assert(not np.isinf(y)),(coeffs,intercepts,zs,y)
         for i in range(len(zs)-1):
             if zs[i]<=y and y<=zs[i+1]:
                 ukval = h(xs[i]) + (y-xs[i])*hprime(xs[i])
                 break
         if w<=np.exp(h(y) - ukval):
             return y
         else:
             xs.append( y )



def gen_alpha(k , n):
   def pdfln(x):
      C = sum([np.log(i) for i in range(1,n+1)])
      return C + (k-3/2)*np.log(x) - 1.0/(2*x) + math.lgamma(x) - math.lgamma(n+x)
   def pdfln_prime(x):
      return (k-3/2)/x + 1.0/(2*x*x) + digamma(x) - digamma(n+x)
   while 1:
      try:
         return ars(pdfln , pdfln_prime , [0.1 , 1.0 , 3.0])
      except:
         continue



def gen_beta(s , w):
   def pdfln(x):
       k = len(s)
       return -3*math.log(x)/2 - 1.0/(2*x) - k*math.lgamma(x/2)+(k*x/2)*math.log(x*w/2) + np.sum((x/2-1)*np.log(s) - (x*w/2)*s) 
   def pdfln_prime(x):
      k = len(s)
      return -(k/2)*digamma(x/2) + 1.0/(2*x*x) + (k/2)*math.log(x/2) + (k*x-3)/(2*x) + np.sum(np.log(s) - w*s)/2 + k*math.log(w)/2
   while 1:
     try:
        return ars(pdfln , pdfln_prime , [0.1 , 1.0 , 3.0])
     except:
        continue


def rnd_gamma(alpha,theta):
   return np.random.gamma(alpha/2 , 2*theta/alpha)


def pdf_normal(x , mean , precision):
    return math.sqrt(precision/(2*np.pi))*np.exp(-precision*(x-mean)**2/2)


def pdf_gamma(x , alpha,theta):
    lprob = (alpha/2-1)*np.log(x) - alpha*x/(2*theta) - math.lgamma(alpha/2) - (alpha/2)*(math.log(2*theta/alpha))
    return math.exp(lprob)


#-- normal gamma distribution pdf
#-- for multivariate case, this should be normal wishart distribution pdf 
def pdf_G0(mu , s , lam , r , beta , w):
    return pdf_normal(mu ,lam , r)*pdf_gamma(s,beta,1.0/w)


def pdf_lewens(k,n,alpha,ns):
    return (k*np.log(alpha)) + sum([math.lgamma(n) for n in ns]) - (math.lgamma(n+alpha) - math.lgamma(alpha))



class DPGMM:
   def __init__(self , alpha , lam , r , beta, w):
       self.alpha = alpha
       self.lam = lam
       self.r = r  
       self.beta = beta
       self.w = w
       self.cnts = []     
       self.params = []   #-- component parameters
       self.indicators = []  #-- for debug
   def __iter__(self):
       return self
   def __next__(self):
       N = sum(self.cnts)
       k = len(self.cnts)
       probs = [n/(N+self.alpha) for n in self.cnts]
       probs.append( (self.alpha)/(N+self.alpha) )
       c = np.random.choice( list(range(k+1)) , p=probs)
       self.indicators.append(c)
       if c<k:
          mu,s = self.params[c]
          self.cnts[c] += 1
       else:
          mu = np.random.normal(self.lam , math.sqrt(1.0/self.r))
          s = rnd_gamma(self.beta , 1.0/self.w)
          self.params.append( (mu,s) )
          self.cnts.append( 1 )
       return np.random.normal(mu , math.sqrt(1.0/s) )
   def likelihood(self,xs):
        k_rep = len(self.cnts)
        N = len(self.indicators)
        Ptmp = pdf_lewens(k_rep,N, self.alpha, np.array(self.cnts))
        for j,(m,s) in enumerate(self.params):
            mask = (np.array(self.indicators)==j).astype(int)
            Ptmp += np.sum( mask*np.log(pdf_normal(xs,m,s)) )
            Ptmp += np.log( pdf_G0(m,s,self.lam , self.r, self.beta, self.w) )
        return Ptmp


def iGMM(xs , iter=200):
    assert(type(xs)==np.ndarray)
    N = len(xs)
    mean_all = np.sum(xs)/N
    var_all = np.sum(xs**2)/N - mean_all*mean_all
    stdev_all = math.sqrt(var_all)
    best_indicators = np.zeros(len(xs) , dtype=np.int)
    best_alpha = 1.0/rnd_gamma(1.0 , 1.0)
    best_beta = 1.0/rnd_gamma(1.0 , 1.0)
    best_lambda = np.random.normal(mean_all , stdev_all)
    best_r = rnd_gamma(1 , 1.0/var_all)
    best_w = rnd_gamma(1 , var_all)
    best_mparams = [mean_all]
    best_sparams = [1.0/var_all]
    k_rep = 1
    Ptmp = 0.0
    for j,(m,s) in enumerate(zip(best_mparams,best_sparams)):
        mask = (best_indicators==j).astype(int)
        Ptmp += np.sum( mask*np.log(pdf_normal(xs,m,s)) )
        Ptmp += np.log( pdf_G0(m,s,best_lambda,best_r,best_beta,best_w) )
    Pmax = Ptmp
    #print(Pmax)
    for it in range(iter):
        #-- step(1) update class indicators
        h = N
        k_rep = np.max( best_indicators ) + 1
        assert(h >k_rep)
        new_indicators = np.copy( best_indicators )
        top_mparams = np.zeros(h)
        top_sparams = np.zeros(h)
        top_mparams[:k_rep] = np.copy(best_mparams)
        top_sparams[:k_rep] = np.copy(best_sparams)
        top_mparams[k_rep:] = [np.random.normal(best_lambda , math.sqrt(1.0/best_r)) for _ in range(0,h-k_rep)]
        top_sparams[k_rep:] = [rnd_gamma(best_beta , 1.0/best_w) for _ in range(0,h-k_rep)]
        for (n,x) in enumerate(xs):
             top_nums = np.bincount( new_indicators )
             k_rep = len(top_nums)
             c = new_indicators[n]
             top_nums[c] -= 1
             Fvals = np.array([pdf_normal(x , mu , s) for (mu,s) in zip(top_mparams,top_sparams)][:h])
             if top_nums[c]==0:
                 k_rep -= 1
                 weights = np.ones(h)*(best_alpha/(h-k_rep))
                 weights[:k_rep+1] = top_nums
                 weights[c] = (best_alpha/(h-k_rep))
             else:
                 weights = np.ones(h)*(best_alpha/(h-k_rep))
                 weights[:k_rep] = top_nums
             Fvals *= weights
             Ftot = np.sum(Fvals)
             new_c = np.random.choice(list(range(h)) , p=Fvals/Ftot)
             if top_nums[c]==0:
                 if (new_c==c or new_c>=k_rep):
                     new_indicators[n] = c
                 else:
                     new_indicators[n] = new_c
                     mask = (new_indicators>=c).astype(int)
                     new_indicators -= mask
                     top_mparams[c:h-1] = top_mparams[c+1:h]
                     top_sparams[c:h-1] = top_sparams[c+1:h]
                     top_mparams[h-1] = 0.0
                     top_sparams[h-1] = 0.0
                     h -= 1
             elif new_c>=k_rep:
                 #-- swap new_c and k_rep
                 new_indicators[n] = k_rep
                 top_mparams[k_rep],top_mparams[new_c] = top_mparams[new_c] , top_mparams[k_rep]
                 top_sparams[k_rep],top_sparams[new_c] = top_sparams[new_c] , top_sparams[k_rep]
             else:
                 new_indicators[n] = new_c
        #-- step(2) update component parameters  
        k_rep = np.max( new_indicators ) + 1
        new_mparams = []
        new_sparams = []
        for j in range(k_rep):
            mask = (new_indicators==j).astype(int)
            nj = np.sum(mask)
            ybar_j = np.sum(mask*xs)/nj
            sj = top_sparams[j]
            mu_j = np.random.normal((ybar_j*nj*sj*best_lambda*best_r)/(nj*sj+best_r) , math.sqrt(1.0/(nj*sj+best_r)))
            var_j = np.sum( mask*(xs-np.ones(N)*mu_j)**2 )
            sj = rnd_gamma(best_beta+nj , (best_beta+nj)/(best_w*best_beta+var_j))
            new_mparams.append( mu_j )
            new_sparams.append( sj )
        #-- step(3) update hyperparameters
        new_lambda = np.random.normal( (mean_all/var_all+best_r*sum(new_mparams))/(1.0/var_all+k_rep*best_r) , 
                                       math.sqrt(1.0/(var_all+k_rep*best_r)) )
        new_r = rnd_gamma(k_rep+1 , (k_rep+1)/(var_all + sum([(x-new_lambda)**2 for x in new_mparams])))
        new_w = rnd_gamma(k_rep*best_beta+1 , (k_rep*best_beta)/(1.0/var_all + best_beta*sum(new_sparams)))
        new_beta = gen_beta(np.array(new_sparams) , new_w)
        new_alpha = gen_alpha(k_rep , N)
        #-- compute P
        Ptmp = pdf_lewens(k_rep,N,new_alpha,np.bincount(new_indicators))
        for j,(m,s) in enumerate(zip(new_mparams,new_sparams)):
            mask = (new_indicators==j).astype(int)
            Ptmp += np.sum( mask*np.log(pdf_normal(xs,m,s)) )
            Ptmp += np.log( pdf_G0(m,s,new_lambda,new_r,new_beta,new_w) )
        P = Ptmp
        #print(P , k_rep)
        if P>Pmax:
            best_indicators = new_indicators
            best_mparams = new_mparams
            best_sparams = new_sparams
            best_alpha = new_alpha
            best_beta = new_beta
            best_r = new_r
            best_w = new_w
            best_lambda = new_lambda
            Pmax = P
        else:
            pass
    return (best_indicators , best_mparams , best_sparams , best_alpha,Pmax)

#-- simple test
if __name__=="__main__":
   m = DPGMM(1.1,0.0 , 15.0 , 1.1,0.5)
   dat=[]
   for _ in range(200):
      dat.append(next(m))
   dat=np.array(dat)
   print( iGMM(dat) )




今後の目標としては、

・infinite HMMを理解する。階層Dirichlet過程というものを使うらしい

・infinite PCFGを理解する

The Infinite PCFG using Hierarchical Dirichlet Processes

https://cs.stanford.edu/~pliang/papers/hdppcfg-emnlp2007.pdf

あたり。Dirichlet過程を使っているものは、須く、Pitman-Yor過程に一般化できるのだと思うので、必要に応じて理解したい。

やりたいこととしては、CCG(Combinatory Categorial Grammar)の確率版であるPCCGを作れるので、"infinite PCCG"ができないだろうかというあたり。CCGはCFGより強いが、どっちもCYK parsingでparseできるし、infinite PCFGができるなら、infinite PCCGができない道理はなさそうな気がする。まぁ、未知の言語の文法構造を、人出によるタグ付けとかなしで、推定できるようになれば嬉しいのだけど、infinite PCFGが流行ってるようにも見えないので、どうなのか(そもそも、文法推論業界で得られてる知見によれば、正規表現程度でも、"ちゃんと"学習するには、正例のみだけではダメで、負例が必要ということなので、もっと強いCFGやCCGでは、尚更、負例の提示が必要そうだけど、今のinfinite PCFGとかは多分正例しか食べ(れ)ないので、そのへんどうなっているか理解したい)

2017-04-09

大雑把ではあるが

お住まい Gr(2,N) $M_{0,N}$ P^1(C) 
線形 wave equation (一般化)超幾何関数 Gauss超幾何関数
非線形 GASDYM方程式   Schlesinger   Painleve VI
量子論(非線形) KZ equation 量子Painleve VI

というような対応表がありそう。


非線形と書いているのは、ゲージ群がSL(2,C)の場合。Gr(2,N)はGrassmann多様体で、Wolf spaceというものである("Representation theory and ADHM-construction on quaternion symmetric spaces"という論文で、ADHM構成が一般化されている)と同時に、(一般化)超幾何関数も定義できる。また、$M_{0,N}$は種数$g=0$のN点付きRiemann面のモジュライ空間で、Gr(2,N)のChow quotientあるいはHilbert quotientというものとして実現できる(arXiv:alg-geom/9210002)。特に、N=4で特殊化したのが、最後の列。


#今の所、$M_{0,N}$の大域的な情報が必要になる場面はない気がするけど。

Equations for Chow and Hilbert quotients

https://arxiv.org/abs/0707.1801


それぞれの方程式がどこで定義されていると考えるのが一番自然であるかについて、コンセンサスがあるのかは分からない。比較的最近の論文を見ても、しばしば適当な局所アフィン座標上でのみ考えられていて、大域的な空間が何であるのかは書かれていないことも多い。表の意味としては、Gr(2,N)上の方程式を適当に簡約すると、$M_{0,N}$上の微分方程式が出るのだけど、希望的観測としては、Gr(2,N)上の方程式の解の"基底"をなす(古典非線形の場合は、線形性がないので、このようには言えないけど)のじゃないかと期待している。N=4の時、Gr(2,N)は複素Minkowski時空の共形コンパクト化であり、共形対称性も複素化して、sl(4,C)を考えるべきだと思うが、Gauss超幾何関数に於いて、この時空の対称性は、近接関係式として現れる。他のケースでは、時空の対称性は、どこに行ってしまったのか(私の知る限り)定かでないのだけど、超幾何関数同様、全部のパラメータをまとめて考えると、時空の対称性が見えてくるのでないかと思う(特に根拠はないけど)。CFTから時空の対称性を抽出する方法が分かれば僥倖という感じ


GASDYM(Generalized anti-self-dual Yang-Mills)方程式については、4次元以外では、2形式のHodge双対が2形式でないので、高次元化は自明でない。

Twistor Theory and the Schlesinger Equations

http://link.springer.com/chapter/10.1007%2F978-94-011-2082-1_3

では、twistor spaceの方で考えている。GASDYM方程式という名前が出てくるのは同著者らによる"Integrability, Self-duality, and Twistor Theory"という本のようである。4次元以外では、自己双対条件と反自己双対条件は、(少なくとも見た目は)互いに"双対"って感じではない


"量子論"の行は、古典論で対応するものは、量子化した後でも対応するだろうというだけのもので、そもそも、GASDYM方程式の量子化がよく分からないので、根拠は乏しい。まぁ、どちらかというと、この対応を手掛かりとして、量子GASDYM方程式を理解したいという気持ちがある

ノンパラメトリックベイズ(0)適応的棄却サンプリング

もうブームは去った気がするけど、ノンパラメトリックベイズの勉強をすることにした。統計学とかで、パラメトリックモデル/ノンパラメトリックモデルという概念があるけど、それとは関係ないように見える。infinite GMMを実装しようとしたら、適応的棄却サンプリング(ARS,Adaptive rejection sampling)を使うとか書いてあって、よく分からなかったので、そこから。


論文

adaptive rejection sampling for gibbs sampling(PDF)

http://www.math.chalmers.se/Stat/Grundutb/CTH/mve186/1415/adaptive.sampling.pdf

によると、元々はDevroyeという人の本に書いてある方法を改良したものらしい(?)。adaptive rejection samplingという名前は、この論文で命名されたのじゃないかと思うけど。


確率密度関数p(x)が解析的な式で与えられている状況で、普通の棄却サンプリングだと、

Rejection sampling

https://en.wikipedia.org/wiki/Rejection_sampling#Theory

提案分布(proposal distribution)は人間がいい感じに設定することになっているが、そんないい分布わかんねーよとかいうケースで、p(x)がlog-concave(文字通り、p(x)の対数が凹関数であるという条件)を満たす時に、p(x)の定義域上で、複数のxを選んで、各xにおいて、log p(x)の接線を引く。それらをつなぐと、log p(x)を上から近似する区分的線形関数が得られる(log-concaveである必要がここで生じる)ので、これのexponentialを取って正規化したものを提案分布とする、ということらしい(論文では、正規化してない関数を$u_k(x)$、正規化した確率密度関数を$s_k(x)$と書いている)。また、接線を計算するため、p(x)は一回以上微分可能である必要がある。逆関数法では、p(x)を積分する必要があるが、積分はできないけど、微分はできるというケースはよくある


区分的線形関数のexponentialを確率密度関数とするような乱数は、(少し面倒ではあるが)普通に逆関数法で生成できる。

逆関数法

https://ja.wikipedia.org/wiki/%E9%80%86%E9%96%A2%E6%95%B0%E6%B3%95

あと、棄却サンプリングで、提案分布に従う乱数xを発生させて、最終的に棄却された時、提案分布を作る時に選んだ"複数のx"に、これも追加して、提案分布を改善する(一番最初の点は人間が選ぶ必要がある。最低二点必要)。提案分布とp(x)の乖離が大きいxでは、棄却率が高く、これによって乖離の大きい部分が改善される可能性が高い。繰り返し行うと、提案分布の近似は、どんどんよくなっていくので"適応的"という名前が付いている。infinite GMMでは、毎回(Gibbsサンプリングの一ループ毎に)パラメータが変わって、提案分布も作りなおしになるけど

論文では、上から抑える関数とは別に下から抑える関数(論文でsqueezing functionと呼んでているもの)も作っていて、多分、元の確率密度関数を評価するよりも、区分線形関数のexponentialを評価する方が、一般的に計算量が少ないから、ということで、そうしているのだと思う。色々な改良やバリエーションがあるようだけど、区分線系関数のexponentialを使って、上(や下)から確率密度関数を抑える関数を作るという点は共通しているっぽい

※)論文2.2.2において、

¥exp(l_k(x^{*}) - u_k(x^{*})) ¥leq ¥exp(h(x^{*}) - u_k(x^{*}))

が成立するので、squeezing testはスキップしても結果は変わらない


最近、頭が悪くなった気がするので、これくらいだと理解できて良い。


以下に、単純な実装を置いておく。棄却されるたびに、h(x)とh'(x)を計算し直すのは、無駄の極みだけど、わかりやすさ重視。例として、ガンマ分布(np.random.gammaで生成できるが)

"""
This is an implementation of the following paper written in python 3.5

Adaptive Rejection Sampling for Gibbs Sampling
https://www.jstor.org/stable/2347565?seq=1#page_scan_tab_contents
"""

import numpy as np

def rnd_from_piecewise_exponential(coeffs , intercepts , points):
    """
       density function h(x) in [points[i] , points[i+1]]
       h(x) = exp(coeffs[i]*x + intercepts[i])
    """
    assert(len(coeffs)==len(intercepts)),(coeffs,intercepts)
    assert(len(coeffs)+1==len(points)),(coeffs,points)
    u = np.random.random()
    CDF = [0.0]  #--cumulative density value at each points
    for i in range(len(points)-1):
       z_cur = points[i]
       z_next = points[i+1]
       a = coeffs[i]
       b = intercepts[i]
       if a!=0.0:
           r = (np.exp(a*z_next+b) - np.exp(a*z_cur+b))/a
       else:
           r = (znext-z)*np.exp(b)
       r0 = CDF[-1]
       CDF.append( r + r0 )
    u = CDF[-1]*u    #--CDFが正規化されてないので
    for i in range(len(points)-1):
        if CDF[i] < u and u <=CDF[i+1]:
             a = coeffs[i]
             b = intercepts[i]
             z = points[i]
             """
                compute t such that
                ( exp(a*t+b) - exp(a*z+b) )/a == u - CDF[i]
             """
             if a!=0.0:
                 t0 = a*(u - CDF[i]) + np.exp(a*z+b)
                 assert(t0 > 0.0),t0
                 t = ( np.log(t0) - b )/a
             else:
                 t = (u - CDF[i])/np.exp(b) + z
             assert(t >= z and t <= points[i+1]),(t,z,points[i+1])
             return t



def ars(h , hprime , points , support=(0 , np.inf)):
     xs = [x for x in points]
     while 1:
         xs.sort()
         x0 = xs[0]
         xN = xs[-1]
         zs = [support[0]]
         coeffs = []
         intercepts = []
         for i in range(len(xs)-1):
             x = xs[i]
             xnext = xs[i+1]
             zi = (h(xnext) - h(x) - xnext*hprime(xnext)+x*hprime(x))/(hprime(x) - hprime(xnext))
             zs.append( zi )
             coeffs.append( hprime(x) )
             intercepts.append( h(x) - hprime(x)*x )
         zs.append( support[1] )
         coeffs.append( hprime(xN) )
         intercepts.append( h(xN) - hprime(xN)*xN )
         y = rnd_from_piecewise_exponential(coeffs , intercepts , zs)
         w = np.random.random()
#--      omit squeezing test
#         lkval = 0.0
#         for i in range(len(xs)-1):
#             if xs[i]<=y and y<=xs[i+1]:
#                 xcur = xs[i]
#                 xnext = xs[i+1]
#                 lkval = ( (xnext - y)*h(xcur) + (y-xcur)*h(xnext) )/(xnext-xcur)
#                 break
#         if w <= np.exp(lkval - ukval):
#             return y
         ukval = 0.0
         for i in range(len(zs)-1):
             if zs[i]<=y and y<=zs[i+1]:
                 ukval = h(xs[i]) + (y-xs[i])*hprime(xs[i])
                 break
         if w<=np.exp(h(y) - ukval):
             return y
         else:
             xs.append( y )


if __name__=="__main__":
   #-- gamma distribution
   alpha = 2.0       #-- shape parameter
   beta = 2.0        #-- rate parameter
   C = np.log(4.0)   #--  log( beta^alpha/Gamma(alpha) )
   def h(x):
      return ((alpha-1)*np.log(x) - beta*x + C)
   def hprime(x):
      return (alpha-1)/x - beta
   print( ars(h , hprime , [0.1 , 1.0 , 5.0]) )


ノンパラメトリックベイズ(0.5)Dirichlet Process Mixture ModelとInfinite GMM

Dirichlet process mixture model(DPMM)とinfinite Gaussian mixture model。まだ実装していないので(0.5)。infinite GMMを実装したら、次はinfinite HMMを理解したい



普通、サンプルデータを混合ガウス分布でfittingしようとする時、混合数は、人間が手動で与えるけど、混合数はいくらでも大きくなり得る(つまり∞)ようにしておいて、うまく推定する方法(母集団の分布として、無限混合モデルを仮定するようなものなので、パラメータが原理的には無限個あるけど、パラメトリックモデルではあると思う。混合数を決めないので、ノンパラメトリックらしいけど、しっくりこないネーミングではある)。Dirichlet process mixture model(DPMM)を最初に考えたのが誰かはよく分からない(多分、1974年のAntoniakの論文?)けど、Rasmussenという人が、2000年にハイパーパラメータも込みで推定する方法を考えて、Infinite Gaussian Mixture Modelという名前が付いた。これ以後、色んなInfinite XXXが生まれることになったっぽい。

The Infinite Gaussian Mixture Model

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.9111


多分、DPMMだと、正規分布の混合モデルである必要もないし、理屈上は、Infinite GMMでも、正規分布以外の混合モデルでも使うことはできるのだろうと思う。が、実用上は正規分布以外を考えることは、あんまりない気もする。Rasmussenの論文だと、MCMCでパラメータ・ハイパーパラメータの分布をサンプリングしてるだけのように見えるけど、多くのケースでは、"一番いい"パラメータ・ハイパーパラメータが欲しい気がする(つまり、ベイズ推定じゃなく、単にMAP推定したい)。これについては、MCMCで十分な数のサンプリングしてやって、その中で確率密度が最大になってるものを選べば、近似的に、真の最大確率密度を与える点を得られるというのが常套手法のよう



Rasmussenの論文では、ハイパーパラメータ(事前分布のパラメータをハイパーパラメータと呼ぶらしい)αの事前分布の導出が省略されすぎてて理解不能なのだけど

Hyperparameter estimation in Dirichlet process mixture models

http://citeseerx.ist.psu.edu/viewdoc/citations;jsessionid=A8264D1702C8FD7E1DDAEF8B1F6BD6B9?doi=10.1.1.51.6826

を読むといい気がする。この論文の一番最初に、Antoniakが与えた式として

P(k|¥alpha,n) = c_n(k) n! ¥alpha^k ¥frac{¥Gamma(¥alpha)}{¥Gamma(¥alpha+n)}

が出てくる。Antoniakの論文は数学の論文なので全体的に読みたくないが、この式は、Ewensの抽出公式というのと、第一種Stiring数の組み合わせ論的解釈からでる。

第一種Stirling数の組み合わせ数学における意味

https://ja.wikipedia.org/wiki/%E3%82%B9%E3%82%BF%E3%83%BC%E3%83%AA%E3%83%B3%E3%82%B0%E6%95%B0#.E7.B5.84.E3.81.BF.E5.90.88.E3.82.8F.E3.81.9B.E6.95.B0.E5.AD.A6.E3.81.AB.E3.81.8A.E3.81.91.E3.82.8B.E6.84.8F.E5.91.B3

Ewens's sampling formula

https://en.wikipedia.org/wiki/Ewens%27s_sampling_formula


#Ewensの抽出公式は、1972年に、集団遺伝学の研究で導かれたものらしい。

これで、例えば、3つの要素{1,2,3}があったらクラスタリングの仕方は{{1,2,3}},{{1},{2,3}},{{2},{1,3}},{{3},{1,2}},{{1},{2},{3}}の4パターンあって、ハイパーパラメータαを決めると、それぞれの分割の出現確率が決まる。αは集団遺伝学の文脈では、突然変異"率"を表し、大きいほど、沢山のクラスターに分割されやすくなる(α→∞の極限で{{1},{2},{3}}への分割確率は1となる)。


正規分布のパラメータの事前分布としては、正規逆ガンマ分布というものを使う(これが、Dirichlet processの基底分布)。多変数(でfull共分散)の場合は、正規逆Wishart分布。この分布にも、パラメータ(ハイパーパラメータ)が入っており、これも決める必要がある。このへんの事前分布の決め方は、"In general the form of the priors are chosen to have (hopefully) reasonable modelling properties, with an eye to mathematical convenience (through the use of conjugate priors)"とか書いてあって、計算に都合のよさそうな形を選んだということ(?)



実装には、ガンマ分布に従う乱数の生成が必要になるけど、例えば、numpyを使う場合、

numpy.random.gamma

https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.gamma.html

の確率密度関数は、Rasmussenの論文の定義とは異なるものになっているに注意が必要。

p(x|k,¥theta) = x^{k-1} ¥frac{e^{-x/¥theta}}{¥theta^k ¥Gamma(k)}

がnumpyの確率密度関数で、論文のは

G(x|k,¥theta) = ¥frac{x^{k/2-1} e^{-k x/2¥theta}}{¥Gamma(k/2) (2¥theta/k)^{k/2}}

となっている。両者は

G(x|k,¥theta) = p(x|k/2 , 2 ¥theta/k)

という関係で結びついている。

スピン3以上のゲージ場(1)

物理に於ける代数的なテンソル計算の例

http://d.hatena.ne.jp/m-a-o/20170131#p1

の"おまけ"の続き。


狭い意味では、ゲージ場は定義によってスピン1であるけども、higher spin gauge fieldという言葉は(最近では?)一般的に使われているようなので、それを踏襲する。これは1996年のVasilievの論文(arXiv:hep-th/9611024)以降市民権を得たのだと思う(少なくとも、1980年頃の論文で用語自体は登場している)(私はVasilievの理論については、殆ど何も知らない)。この名称の背景には、質量0の(スピン0でない)ボソンは、"ポテンシャル"を持つことができ、スピン1の場合は、ベクトルポテンシャルと一致するという事実があり、これを以って、質量0の(スカラー場でない)ボソン場をhigher spin gauge fieldと呼ぶ。Vasilievの論文で引用されているいくつかの論文では、表現論的な議論は特にされていないが、

[1] Massless Fields as Unitary Representations of the Poincare Group

http://onlinelibrary.wiley.com/doi/10.1002/prop.19790270403/abstract

は表現論的視点から類似の議論を行っているように見える。これは1979年の論文で、大体この頃、higher spin gauge fieldに関するコンセンサスができたっぽい


通常の電磁場と一光子状態は、形式的には同じ形のMaxwell方程式を満たす(従って、共に"ベクトルポテンシャル"を持つ)ものの、前者は実数値の場で特に規格化されていない一方、後者は複素数値の場で自然なユニタリ内積が定まっていて(多分)規格化されているなど、大きな違いがある(規格化されていると積極的に考えるべき理由はないけど、定数倍しても観測量に影響がない)。Vasilievの理論は古典論であるので、前者のタイプの場を扱っている一方、論文[1]で扱っているのは、後者のタイプの場である。この2つの場の関係は、物理の教科書の説明では、以下のような感じだと思う。

『通常の電磁場を第二量子化することによって、光子の状態空間を得ることができ、この中には、光子数が正確にN個の状態が含まれていて、特に光子数が1個の状態からなる部分空間を取ってくることができる。これが、Poincare代数のスピン1のmasslessユニタリ表現(ヘリシティ±1のmassless既約ユニタリ表現の直和)と同値になる(はず)。そして、光子の状態空間は、この一光子状態空間の対称テンソル積として得られる』

この関係は、数学的に証明されている事柄ではないと思うので、どっちのタイプの場を扱ってるかは意識しておいたほうがいい気がする。以下では「相対論的場の理論の一粒子状態はPoincare代数の既約ユニタリ表現で分類される」という、教科書レベルで広く受け入れられてる(と思う)認識を元に、(0でない)整数ヘリシティを持つmassless既約ユニタリ表現に対応する一粒子状態が、"higher spin gauge fieldのsingle particle state"であるということにする。これは、Vasiliev理論で扱われている場とは正確には異なるものだ(と思う)けど、多くの数学的性質を共有しているっぽい



整数スピンsを持つmassless粒子は、Lorentz代数の表現D(s,0)+D(0,s)に値を取る場となっている。これはtwistor理論で完全に確立している事実で、ferimionでも同様に成立し、wave equationも分かっている。複素表現としては、D(s,0)とD(0,s)の2つの既約表現があり、helicityの正負と対応している。s=1の時は、D(1,0)+D(0,1)は2階の反対称テンソル場なので、これらの場がスピンsを持つゲージ粒子のfield strengthの場であると推測するのは自然。Lorentz代数の表現を忘れて単にベクトル空間として見れば

D(s,0) ¥simeq S^{2s}(¥mathbf{C}^2) ¥simeq ¥mathbf{C}^{2s+1}

D(0,s) ¥simeq S^{2s}(¥mathbf{C}^2) ¥simeq ¥mathbf{C}^{2s+1}

という同型があり、twistor理論で"既約なspinor場は対称spinor場となる"という事実(単にsl(2,C)の有限次元既約表現は自然表現の対称テンソル積表現で尽きるという事実の帰結)に関連している。論文[1]では、D(s,0)+D(0,s)が

¥bigotimes^{2s} D(1/2,1/2)

の部分表現として実現できる(D(1/2,1/2)はベクトル表現)という事実を使って、ランク$2s$のテンソル場として、扱っている。普通の(スピン1の)ゲージ場の曲率は2階の反対称テンソルで、Weylテンソルは4階のテンソルであることを思い出せば、自然な見方であるのかもしれない。


疑問)$V=D(1/2,1/2)$として、2階の反対称テンソルの空間は、分割(1,1)に対応するYoung symmetrizerの像として記述される。この空間は、gl(V)の既約表現であると同時に、so(4,C)の既約表現でもある。前回書いたように、分割(2,2)に対応するYoung symmetrizerの像は、代数的曲率の空間と呼ばれ、gl(V)の既約表現であるが、so(4,C)の表現としては可約である。この分解から、WeylテンソルとRicciテンソルが生じる(WeylテンソルはD(2,0)+D(0,2)で、RicciテンソルはD(1,1)+D(0,0)になる)。多分、一般相対論の教科書を見ると、Weylテンソルの名前が出ないことはあってもRicciテンソルの名前が出ないことはないだろうと思われるけど、gl(V)による既約分解を見ることは、意味のあることなのかもしれない。そうすると、スピンsの場合でも、代数的曲率空間に相当する何かがあるだろうか(例えば、分割(s,s)に対応するYoung symmetrizerの像だとか)。論文[1]の式(3.6)(I1)(I2)は、この高階テンソルがS^s(¥Lambda^2(V))に含まれるという条件に見える。



ゲージポテンシャル。スピン1の場合は、ベクトル場であり、D(1/2,1/2)に値を取る場となっていて、スピン2の場合は、D(1,1)場は対称トレースレステンソル場となるので、一般にはLorentz代数の表現D(s/2,s/2)に値を取る場をポテンシャルとして持つのだろうと予想される。s/2は半整数でないといけないので、これは、bosonの場合しか意味がない(massless fermionがポテンシャル場を持たないとは言えないが)。一般に

S^{k+2}(D(1/2,1/2)) ¥simeq D(k/2+1,k/2+1) ¥oplus S^k(D(1/2,1/2)

S^k(D(1/2,1/2)) ¥simeq ¥bigoplus_{j=0}^{¥[k/2¥]} D(k/2-j,k/2-j)

S^2(D(1/2,1/2)) ¥simeq D(1,1) ¥oplus D(0,0)

S^3(D(1/2,1/2)) ¥simeq D(3/2,3/2) ¥oplus D(1/2,1/2)

S^4(D(1/2,1/2)) ¥simeq D(2,2) ¥oplus D(1,1) ¥oplus D(0,0)

etc.

という既約分解が成立するので、物理の文献では、単にhigher spin gauge potential field=高階対称テンソル場として扱われているようである。


(※)D(1/2,1/2)の基底がe_0,e_1,e_2,e_3であるとして、Lorentz代数の標準的な表現を考えるとe_0 ¥otimes e_0 - e_1 ¥otimes e_1 - e_2 ¥otimes e_2 - e_3 ¥otimes e_3は明らかに2階の対称テンソルで、かつLorentz代数の作用は全て0となることは容易に分かる。これが2階の対称テンソルに於けるD(0,0)成分となっている。一般に2つの添字でトレースを取ることで

S^{k+2}(D(1/2),D(1/2)) ¥to S^k(D(1/2,1/2)

A_{¥mu_1 ¥mu_2 ¥cdots ¥mu_{k+2}} ¥mapsto A^{¥lambda}_{¥lambda ¥mu_3¥cdots ¥mu_{k+2}}

のように(k+2)階の対称テンソルから、k階の対称テンソルへの射影が得られる(対称テンソルなので、どの2つの添字を選んでも同じ)


#最近、「古典群の表現論と組み合わせ論(下)」という本を見ていたら、シンプレクティック群と直交群では、(symmetricとは限らない)traceless tensorの空間を、Schur-Weylの時と同様に分解することができて、更に、一般の直交群の対称テンソル積表現の既約分解が、上の分解と同じように与えられていた(108ページ、命題8.26)。最高ウェイトベクトルもYoung symmetrizerを使って計算することができる。


#局所的には、一般相対論の計量は、$GL(4,R)/O(3,1)$に住んでいると言える。この計量の空間には線形性がない一方、上のスピン2のゲージポテンシャルは線形性がある。Lie環のレベルでは

gl(4,¥mathbf{C}) ¥simeq o(4,¥mathbf{C}) ¥oplus S^2(¥mathbf{C}^4)

という分解が成立し、右辺の第二成分は、2階の対称テンソルになっている(複素化しているのは、量子論では場が複素数値になるので、計量なども複素化されるべきだろうという理由から)。というわけで、少なくとも、線形近似の範囲で、計量はスピン2のゲージポテンシャルと見なすことが正当化される(と思われている)



次に、ポテンシャルとfield strengthの関係について。数学的に納得できる記述は見たことがなかったけど、論文[1]の8節に書いてある(完)。これを理解するには、Lorentz代数の表現だけ見ていても不十分で、「場」が必要になる。Lorentz代数の(既約とは限らない)有限次元表現空間$V$に対して、「$V$に値を取る場」というのは、通常

V ¥otimes C^{¥infty}(¥mathbf{R}^4)

の"ようなもの"をイメージしているのだと思う。4次元Euclid空間は、実Minkowski時空の場合もあるし、物理の人は、運動量空間をよく使う。これは論文[1]でもそんな感じ。こういうのは、代数的に扱いづらいので、もっと限定した形の空間として

V ¥otimes H_0

を考える。H_0はPoincare代数のhelicity 0のmassless既約表現。とりあえず、重要そうな話として

(1)この空間には、Poincare代数の表現が定義できる(一般にユニタリーではないが)

(2)この空間では、massless Klein-Gordon方程式が、いつでも自動的に成立している。

の2点が確認できる


H_0

Maxwell方程式の表現論(2)Gupta-Bleuler量子化の表現論的側面

http://d.hatena.ne.jp/m-a-o/20151106#p2

ではV_0と書いているのと同じもので

H_0 = ¥{f ¥in ¥mathbf{C}¥[z_1,z_2,z_3,z_4¥] | (z_1 ¥partial_1 + z_2 ¥partial_2-z_3 ¥partial_3 - z_4 ¥partial_4)f = 0¥}

という実現を持つ。Poincare代数の作用の仕方は、前にも書いたので省略。例えば、2階の反対称テンソル場は、以下のデータ

F_{¥mu ¥nu} ¥in H_0 (¥mu,¥nu=0,1,2,3)

F_{¥mu ¥nu} = -F_{¥nu ¥mu}

と同一視できる。というわけで、場とは、ある条件を満たす4変数多項式の多次元配列である(テストで書いたら0点になりそうな定義)。一般に、helicity hのmassless既約ユニタリ表現は

H_{h} = ¥{f ¥in ¥mathbf{C}¥[z_1,z_2,z_3,z_4¥] | (z_1 ¥partial_1 + z_2 ¥partial_2-z_3 ¥partial_3 - z_4 ¥partial_4)f = 2 h f¥}

で実現できる(Poincare代数の作用の仕方は全て同一)



以上の準備のもとで、論文[1]8節に書いてあることを、そのまま読み替えると、例えば、スピン2のゲージポテンシャルからfield strengthへの写像S^2(D(1/2,1/2)) ¥otimes H_0 ¥to D(1/2,1/2)^{¥otimes 4} ¥otimes H_0

C_{¥mu ¥nu ¥rho ¥sigma} = P_{¥mu} P_{¥rho} h_{¥nu ¥sigma} - P_{¥nu} P_{¥rho} h_{¥mu ¥rho} - P_{¥mu} P_{¥sigma} h_{¥nu ¥rho} + P_{¥nu} P_{¥sigma} h_{¥mu ¥rho}

となり、これは線形化したWeylテンソルとなっている(この写像がPoincare代数の作用と可換であることは別途確かめる必要がある)。これは、S^2(V) ¥otimes H_0全体で定義されるが、h_{ab}はtracelessと仮定することで、D(1,1) ¥otimes H_0上に制限できる。


また、スピン1の時と同様、以下のゲージ変換が存在する

h_{ab} ¥mapsto h_{ab} + P_{a} ¥eta_{b} + P_{b} ¥eta_{a}

ここでh_{ab}に、traceless条件とHilbertゲージ条件を課すことにすると、任意のゲージ変換が許されず¥eta_{a}にも条件が必要となる。traceless条件は単に

P^a ¥eta_{a} = 0

という条件になる。traceless条件がある場合、Hilbertゲージは

P^a h_{ab} = 0

となり、ゲージ変換部分が0となるには

P^a (P_a ¥eta_{b} + P_b ¥eta_{a}) = (P_0^2-P_1^2-P_2^2-P_3^2)¥eta_{b} + P_b(P^a ¥eta_{a}) = 0

で、第一項は、massless Klein-Gordon方程式より0になり、第二項は、さっきと同じ条件

P^a ¥eta_{a} = 0

から0になる(本来のHilbertゲージ条件は、任意のゲージ変換で保たれるので当然。このへんの計算は、重力波の計算で出てくるものと形式的には同一)。ところで、これはLorenzゲージ条件となっている。つまり、

V_{fix}^{(0)} = H_0

V_{fix}^{(1)} = ¥{ A_{a} ¥in V ¥otimes H_0 | P^a A_{a} = 0 ¥}

V_{g}^{(1)} = ¥{ A_{a} ¥in V ¥otimes H_0 | A_a = P_a f , f ¥in V_{fix}^{(0)} ¥}

V_{fix}^{(2)} = ¥{ h_{ab} ¥in S^2(V) ¥otimes H_0 | P^a h_{ab}=0,h^a_{a} = 0 ¥}

V_{g}^{(2)} = ¥{ h_{ab} ¥in S^2(V) ¥otimes H_0 | h_{ab} = P_a v_b + P_b v_a , v_a ¥in V_{fix}^{(1)} ¥}

などと置くと、これらは全部Poincare代数の表現空間となっていて

V_{phys}^{(1)} = V_{fix}^{(1)}/V_{g}^{(1)} ¥simeq H_{+1} ¥oplus H_{-1}

がスピン1のmassless粒子の空間(helicity +1/-1の既約ユニタリ表現の直和)で

V_{phys}^{(2)} = V_{fix}^{(2)}/V_{g}^{(2)} ¥simeq H_{+2} ¥oplus H_{-2}

がスピン2のmassless粒子の空間(helicity +2/-2の既約ユニタリ表現の直和)となる(と思う)。沢山の表現空間と写像があって(これらの写像は、大体、Poincare代数の作用と可換となっているはず)、ややこしいけど、曲率写像

 C: V_{fix}^{(2)} ¥to V^{¥otimes 4} ¥otimes H_0

については

 Ker(C) = V_{g}^{(2)}

が成立するはず。スピン1の時の計算は以前やったけど、Im(C)の具体的な記述も同様にできると思う(が、面倒なのでやらない)。論文[1]にはスピン3以上の場合の条件も書いており、同じような形となっている。


論文[1]にも書いてある通り、V_{fix}^{(s)}は、traceless条件から、D(s/2,s/2) ¥otimes H_0に含まれることが言え、これは[dim D(s/2,s/2)=(s+1)^2]個の独立な成分を持つ場で、ゲージ固定でdim D((s-1)/2,(s-1)/2)=s^2の自由度が消え、2s+1の自由度が残る。ところで、V_{g}^{(s)}V_{fix}^{(s-1)}から来るので、2(s-1)+1=2s-1個分の自由度はゲージ変換の自由度であり、従って(2s+1)-(2s-1)=2成分が残る。これは、電磁波、重力波に共通する自由度で、ヘリシティ正負の2成分があることに対応する。ということから、論文には書いてないけど

V_{fix}^{(s)} ¥simeq V_{g}^{(s+1)}

が成立すると思われる(多分)



それなりに色々な計算をした結果得られるものがH_{+s} ¥oplus H_{-s}なので、苦労した意味はあるのか、という疑問はある(古典的には、ポテンシャルを考えることは、物質場との相互作用を書くのに必須だったけど)


※)QED(Quantum Electrodynamics)の状態空間というのは、(電子と光子のみを考える場合)、Poincare代数のスピンs、質量mのユニタリ表現をV(s,m)とする(m>0の時は、通常の既約ユニタリ表現。m=0の時は、massless既約ユニタリ表現2つの直和)時、

H_{QED} = S^{*}(V(1,0)) ¥otimes ¥Lambda(V(1/2,m))

のはず。V(1,0)の対称テンソル積空間には、V(1,0)の基底ごとにボソンの生成・消滅演算子が作用し、V(1/2,m)の外積空間には、V(1/2,m)の基底ごとにフェルミオンの生成・消滅演算子が作用する。H_{QED}は可算個の基底しか持たないので、電子と光子の相互作用がなければ、H_{QED}にPoincare代数の表現が普通に定義できる。これは自由場の理論なので、特に凄くも何ともないけど、現代の場の量子論の標準的な考えでは、相互作用を入れても、H_{QED}そのものは変わらない(多分)。変わるのは、Poincare代数の作用の仕方で、例えば、Hamiltonian(=Poincare代数の無限小時間並進)には相互作用項が入り、形式的には、その形はよく知られている。Poincare代数では

¥[K_i , P_i¥] = i P_0 (i=1,2,3)

なので(ここで、K_i,P_iはLorentz boostと運動量演算子、P_0はHamiltonian)、Hamiltonian以外にも相互作用項が入らないと、Poincare代数の作用を整合的に変形できない(はず。全く計算してないけど、運動量と角運動量は、相互作用項とかでないので、boost operatorに付加的な項が出て、適切な交換関係を満たすのだと思う。教科書とかで、そんな記述見たことないけど)。そんで、Hamiltonianの相互作用項は、形式的には

H_{int} ¥propto  ¥displaystyle ¥int ¥frac{d ¥mathbf{p} d ¥mathbf{k}}{¥sqrt{|¥mathbf{k}|}} (c_{¥mathbf{k}}^{+} a_{¥mathbf{p}}^{+}a_{¥mathbf{p}+¥mathbf{k}} + c_{¥mathbf{k}} a_{¥mathbf{p}}^{+}a_{¥mathbf{p}-¥mathbf{k}})

という形で、これを

H_{int} = ¥sum_{¥alpha,i,j} p_{¥alpha,i,j} c_{¥alpha}^{+} a_{i}^{+} a_j + ¥sum_{¥alpha,i,j} q_{¥alpha,i,j} c_{¥alpha} a_{i}^{+} a_{j}

というふうに書けるだろうと思う(p,qは係数で、具体的な値を計算して決定できるはず)。但し、c_{¥alpha}^{+}は光子の生成演算子、a_{i}^{+}は電子の生成演算子、a_{j}は電子の消滅演算子で、¥alphaは一光子の状態空間V(1,0)の基底を走り、i,jは一電子の状態空間V(1/2,m)の基底を走る



結局、(少なくともQEDくらいの)場の量子論が数学的に怪しいことになるのは平面波基底を使うからで、平面波基底は、上のQED状態空間の元でないし(勿論、並進演算子の連続スペクトル状態ではある)(そもそも、数学的な意味では厳密には基底ですらない)、上のQED状態空間は、可算個の基底を持つので、それでHamiltonianとか他の演算子を、ちゃんと書けば何も問題は起きないのでないかというようなことを、1.5年前くらいから考えていて、頑張れば、数学的に厳密なQED Hamiltonianを書けるのじゃないかと思うけど、一向に計算を進める気が起きないので、いくつかやるべき計算について書いておく(各問題につき、できたら10万円支払いますので、できた方は教えてください)


問題1:de Sitter代数so(4,1)の適当な既約ユニタリ表現のcontractionから、Poincare代数iso(3,1)のmassive既約ユニタリ表現を具体的に基底を取って書く(スピン0,1/2の場合くらいあれば十分)


問題2:問題1で計算したPoincare代数のmassive既約ユニタリ表現の基底に対応するMinkowski時空上の関数を具体的に決定すること(これも、スピン0と1/2の場合だけでOK)



問題1について。so(4,1)の既約ユニタリ表現は一杯あるけど、

GENERALIZED EIGENVECTORS AND GROUP REPRESENTATIONS - THE CONNECTION BETWEEN REPRESENTATIONS OF SO(4, 1) AND THE POINCARE GROUP

http://link.springer.com/chapter/10.1007%2F978-94-010-2669-7_8

を読むと、できそうな感じ。練習として、so(2,1)=su(1,1)の既約ユニタリ表現から、contractionによって、iso(1,1)の既約ユニタリ表現を作るというようなことをやるといいと思う。iso(1,1)の(contractionによらない)既約ユニタリ表現については

(1+1)-次元Poincare代数の表現論

http://formalgroup.tumblr.com/post/124635771365/1-1-%E6%AC%A1%E5%85%83poincare%E4%BB%A3%E6%95%B0%E3%81%AE%E8%A1%A8%E7%8F%BE%E8%AB%96

などを参考。



問題2は必須というわけでもないけど、masslessの場合の計算は、

Maxwell方程式の表現論

http://d.hatena.ne.jp/m-a-o/20150621#p1

でやっているので、massive版も一応やろうという話。

若干気になるのは、masslessの場合、光円錐上で特異性を持っていて、以前計算した時は、ただの有理関数として代数的に扱ったけれども、massiveの時は、多分Bessel関数とか出てきて代数的に扱えないかもしれない。それでも、やっぱりある超曲面上で特異になるんじゃないかと思うので、そのへんの扱いをどうするか。予想では、基底は全て同じ超曲面上に特異性を持つので、この超曲面を除いた領域で考えればいいと思う。masslessの場合は、どの基底もt^2-x^2-y^2-z^2=0上でのみ値が定まらず発散する。勿論、Poincare群の作用には並進があり、並進作用によって、光円錐は動くので、一見奇妙な気がするけど、並進は、Poincare代数の演算子のexponentialで、基底の(可算)無限個の和となるので、光円錐が動くことに問題はない。有限個の基底の和を考える限りは、同じ光円錐t^2-x^2-y^2-z^2=0上で特異となる


#似たような状況として、1/x,1/x^2,1/x^3,...を考えると、x>aにおいて

1/(x-a) = 1/x ( 1 / (1 - a/x) ) = 1/x(1 + a/x + (a/x)^2 + .... )

という等式が成立する。1/(x-a)は1/xを並進で動かしたものでx=aで特異であるけど、一方、右辺は、1/x,1/x^2,1/x^3,...の無限和であり、個々の関数はx=0で特異となっている

あまり綺麗な解決というわけでもないので、もっといい理解の仕方があるかもしれない。



問題をやるのに、物理の知識は微塵も必要ないけど、数学を多少出来る必要はある。これができれば、QED状態空間の具体的な記述がわかり、相互作用のない場合は、Hamiltonian、運動量、角運動量、boost演算子の全てが具体的に書ける。相互作用がある場合は、上に書いたように、Hamiltonianとboost operatorに相互作用項が入ると思うので、それを決めればいい。まぁダメだったら私の頭が悪かったというだけだし、これで、"QED Hamiltonian"が書けても、それが本当に通常のQEDの結果を再現するのかということは別途示す必要があるので、先は長い。

水素原子の表現論・零(二次元の場合)

水素原子の表現論

http://d.hatena.ne.jp/m-a-o/20140130#p1

で、4次元だと、どうなるか考えたりしたけど、一番簡単な2次元をやってないことに、気付いたので、計算してみた

2次元量子kepler系の束縛状態

http://formalgroup.tumblr.com/post/153902406125/2%E6%AC%A1%E5%85%83%E9%87%8F%E5%AD%90kepler%E7%B3%BB%E3%81%AE%E6%9D%9F%E7%B8%9B%E7%8A%B6%E6%85%8B


普通にやっても面白くないけど、3次元のKustaanheimo-Stiefel変換(KS変換)と類似するLevi-Civita変換があり、調和振動子に帰着させると計算が楽(個人的には、いくつもある計算方法の中で一番簡単な気がする)。Levi-Civita変換では、調和振動子の空間で$O(1)$の作用があって($(u_1,u_2) \to (-u_1 , -u_2)$)、これによって、Levi-Civita変換の行き先は不変に保たれる。そして、O(1)の作用は固有値±1を持ち、シンプレクティック代数の表現空間(振動子表現)も、それに従って既約分解し、その一つ(不変部分を取るので、O(1)作用の固有値+1部分)が水素原子の束縛状態の空間となる。



本編(二次元の場合)終わり。で、3次元に戻ると、KS変換がある。それは例えば

x = 2(u_1 u_3 + u_2 u_4)

y = 2(-u_1 u_4 + u_2 u_3)

z = u_1^2+u_2^2-u_3^2-u_4^2

などと書ける。今度は、4次元空間から3次元空間への写像となっている。KS変換は、3次元球面に制限すると、Hopf写像と同じものになる(S^3のS^2上のHopf fibrationを与える)。

Hopf Map

http://mathworld.wolfram.com/HopfMap.html

また、quaternionで回転行列を書こうとする(つまり、SU(2)をquaternionのノルム1の元と同一視して、SU(2)からSO(3)への準同型を作る)と、行列の成分に、やはり同じ式が出てくるのは、CG関係の仕事とかしてる人は見たことあると思う。他にも

Cayley-Klein parameterの量子化と量子化に関わる諸々

http://d.hatena.ne.jp/m-a-o/20100322#p2

でも、同じ式が出ている。KS変換はどこにでもあるような、ありふれた式と言える



Levi-Civita変換の時と同様、KS変換でも$(u_1,u_2,u_3,u_4)$の取り方には、任意性があり、$u$空間上のある群作用によって、$(u_1,u_2,u_3,u_4)$を動かしても、$(x,y,z)$の値は変わらないということが起きる。二次元Kepler系の場合と違って、今度は、連続群で、$U(1)$と同型になる。連続群なので、無限小変換を考えればいい。生成子は

 L = u_1 ¥frac{¥partial}{¥partial u_2} - u_2 ¥frac{¥partial}{¥partial u_1} - u_3 ¥frac{¥partial}{¥partial u_4} - u_4 ¥frac{¥partial}{¥partial u_3}

となる。Lをx,y,zに作用させて消えることを確認すればいい(これで尽きるということを言うには、多少議論が必要だけども、省略)


Lは、4次元の調和振動子の波動関数に作用する。調和振動子の波動関数の空間というのは、振動子表現のことといっていい。シンプレクティック代数sp(8,R)の中で、$L$が生成するLie環$u(1)$と交換する代数は$u(2,2)$であり、これは水素原子のspectrum generating algebra$so(4,2)$と$L$自身からなる。これは、Howe dualityの簡単な例となっていて、$(u(1),u(2,2))$はdual pairであるという。特に、Lの作用で消える、調和振動子の波動関数全体が、水素原子の束縛状態の空間であり、$so(4,2)$の極小表現である。$L$の固有値で振動子表現を固有空間分解すると、0以外の固有値に属する固有空間にも$so(4,2)$の作用があり、既約表現を得る。これが、ladder representationというものになる(4次元共形代数の整数helicityを持つmassless既約ユニタリ表現でもある。この場合、Lの固有値はヘリシティの二倍と等しい)。


参考)調和振動子の対称性とspectrum generating algebra

http://d.hatena.ne.jp/m-a-o/20161117#p1


古典論に於いて、Lに対応するのは、調和振動子の相空間への$U(1)$作用に関する運動量写像

¥Phi = q_1 p_2 - q_2 p_1 - q_3 p_4 + q_4 p_3 ¥in Map(T^{*}(¥mathbf{R}^4 ¥backslash ¥{0¥}) ¥to ¥mathbf{R})

なので、$L$は量子化運動量写像であるという見方もある。そうすると、水素原子のspectrum generating algebraは、調和振動子の対称性のうち、$U(1)$作用に関する量子化運動量写像と交換するもののみが、(quantum) Hamiltonian reductionの後も残る、という理解になる($L$自身は、水素原子の束縛状態に作用すると0になるので実質的に見えなくなってしまう)

#普通、quantum Hamiltonian reductionという時は、代数の簡約だけを考えて、表現については、何も言わないので、ここの言葉の使い方はあまり正しくはない

#KS変換と、水素原子のspectrum generating algebraは1960年代に見つかっているけど、上のような運動量写像を考えて、Kepler系の相空間を、調和振動子系の相空間のHamiltonian reductionとして理解する、というのは、もうちょっと後の1980年代くらいの論文あたりで普通になっているようである。真面目にsurveyしたわけではないので、もう少し古い論文もあるかもしれないけど、MarsdenとWeinsteinがsymplectic reductionの論文を出したのが1974年らしいから、何となく気付いていた人はいるのかもしれないけど、明確化したのは、それより後じゃないかと思う。振動子表現/Weil表現も、数学の方では、1960年代前半あたりのよう。Howe dualityは1970年代にpreprintを出したけど、完全版が出版されたのは1989年だったらしい(と、Wikipediaに書いてあった)


水素原子のspectrum generating algebraは、原論文を読むと天下り的に出てきた感が強い(とはいえ、本当に何もないところからは出てこない。角運動量とRLベクトルは知られており、一方、動径方向にはsu(1,1)作用があることが知られていたので、それを統一するという視点で得たのでないかと思う)けど、KS変換は、それに比べれば、随分自然にも見える(それで調和振動子に帰着できるというのは、やっぱり不思議なことではあるけど)ので、初学者向けの説明としては、例えば

(1)調和振動子の説明。特に、Heisenberg代数の既約表現とシンプレクティック代数の振動子表現の説明

(2)KS変換による水素原子の束縛状態のエネルギー準位計算

(3)KS変換から量子化運動量写像とspectrum generating algebraを導出

みたいな順序で説明すると、比較的抵抗が少なくて済みそう(表現論前提ではあるので、初めて量子力学を勉強する人向きかどうかは知らないけど)



更に高次元。2次元Kepler系ではLevi-Civita変換ときて、3次元Kepler系ではKS変換と来たけど、4次元Kepler系には、直接の類似物はないっぽい(束縛状態の空間のso(5,2)の既約ユニタリー表現としての構造は別の方法で調べることはできる。任意次元に於けるspectrum generating algebraとKepler系の表現論的記述は、既に知られている)。少し視点を変えると、次に"Hopf fibration"が現れるのはS^7->S^4であり、

参考)Quaternionic Hopf fibrations

https://en.wikipedia.org/wiki/Hopf_fibration#Quaternionic_Hopf_fibrations

5次元Kepler系の相空間が8次元調和振動子のHamiltonian reductionで出ると予想される。ファイバーはS^3であり、これはSp(1)=SU(2)=Spin(3)という群構造を持つ。これは、既に結果が知られている

The geometry of the SU(2) Kepler problem

http://www.sciencedirect.com/science/article/pii/039304409090004M

The Sp(1)-Kepler Problems

https://arxiv.org/abs/0805.0840


試しに少し計算してみると、quaternionic Hopf mapも二次の多項式写像で書けて

x_1 = u_1^2 + u_2^2 + u_3^2 + u_4^2 - u_5^2 - u_6^2 - u_7^2 - u_8^2

x_2 = 2(u_1 u_5 - u_2 u_6 - u_3 u_7 - u_4 u_8)

x_3 = 2(u_1 u_6 + u_2 u_5 + u_3 u_8 - u_4 u_7)

x_4 = 2(u_1 u_7 - u_2 u_8 + u_3 u_5 + u_4 u_6)

x_5 = 2(u_1 u_8 + u_2 u_7 - u_3 u_6 + u_4 u_5)

r=|x| = ¥sum_{i=1}^8 u_i^2

となる。

次の3つの無限小変換で、$x_1,x_2,x_3,x_4,x_5$は消え、これらはLie環$su(2)$と同型な閉じた実Lie環をなす。

D_1 = (u_1 ¥partial_2 - u_2 ¥partial_1) - (u_3 ¥partial_4 - u_4 ¥partial_3) - (u_5 ¥partial_6 - u_6 ¥partial_5) - (u_7 ¥partial_8 - u_8 ¥partial_7)

D_2 = (u_1 ¥partial_3 - u_3 ¥partial_1) + (u_2 ¥partial_4 - u_4 ¥partial_2) - (u_5 ¥partial_7 - u_7 ¥partial_5) - (u_6 ¥partial_8 - u_8 ¥partial_6)

D_3 = (u_1 ¥partial_4 - u_4 ¥partial_1) - (u_2 ¥partial_3 - u_3 ¥partial_2) - (u_5 ¥partial_8 - u_8 ¥partial_5) - (u_6 ¥partial_7 - u_7 ¥partial_6)

確かめていないけど、シンプレクティック代数sp(16,R)の中で、この$D_1,D_2,D_3$と交換するもの全体が、Lie環$so(6,2)$と同型になるのだと思う(上の論文によれば$(su(2),so(6,2))$もdual pair)。一応、Risa/Asirで以下の計算は確認した。


def assert(V){
   if(!V){error(V);}
   return 1;
}


X1 = u1^2 + u2^2 + u3^2 + u4^2 - u5^2 - u6^2 - u7^2 - u8^2;
X2 = 2*(u1*u5 - u2*u6 - u3*u7 - u4*u8);
X3 = 2*(u1*u6 + u2*u5 + u3*u8 - u4*u7);
X4 = 2*(u1*u7 - u2*u8 + u3*u5 + u4*u6);
X5 = 2*(u1*u8 + u2*u7 - u3*u6 + u4*u5);

/*以下は、全部0になる。*/

assert(diff(X1,u1)*u2-diff(X1,u2)*u1-diff(X1,u5)*u6+diff(X1,u6)*u5+( -diff(X1,u3)*u4+diff(X1,u4)*u3-diff(X1,u7)*u8+diff(X1,u8)*u7)==0);
assert(diff(X2,u1)*u2-diff(X2,u2)*u1-diff(X2,u5)*u6+diff(X2,u6)*u5+( -diff(X2,u3)*u4+diff(X2,u4)*u3-diff(X2,u7)*u8+diff(X2,u8)*u7)==0);
assert(diff(X3,u1)*u2-diff(X3,u2)*u1-diff(X3,u5)*u6+diff(X3,u6)*u5+( -diff(X3,u3)*u4+diff(X3,u4)*u3-diff(X3,u7)*u8+diff(X3,u8)*u7)==0);
assert(diff(X4,u1)*u2-diff(X4,u2)*u1-diff(X4,u5)*u6+diff(X4,u6)*u5+( -diff(X4,u3)*u4+diff(X4,u4)*u3-diff(X4,u7)*u8+diff(X4,u8)*u7)==0);
assert(diff(X5,u1)*u2-diff(X5,u2)*u1-diff(X5,u5)*u6+diff(X5,u6)*u5+( -diff(X5,u3)*u4+diff(X5,u4)*u3-diff(X5,u7)*u8+diff(X5,u8)*u7)==0);

assert(diff(X1,u1)*u3-diff(X1,u3)*u1-diff(X1,u5)*u7+diff(X1,u7)*u5-( -diff(X1,u2)*u4+diff(X1,u4)*u2-diff(X1,u6)*u8+diff(X1,u8)*u6)==0);
assert(diff(X2,u1)*u3-diff(X2,u3)*u1-diff(X2,u5)*u7+diff(X2,u7)*u5-( -diff(X2,u2)*u4+diff(X2,u4)*u2-diff(X2,u6)*u8+diff(X2,u8)*u6)==0);
assert(diff(X3,u1)*u3-diff(X3,u3)*u1-diff(X3,u5)*u7+diff(X3,u7)*u5-( -diff(X3,u2)*u4+diff(X3,u4)*u2-diff(X3,u6)*u8+diff(X3,u8)*u6)==0);
assert(diff(X4,u1)*u3-diff(X4,u3)*u1-diff(X4,u5)*u7+diff(X4,u7)*u5-( -diff(X4,u2)*u4+diff(X4,u4)*u2-diff(X4,u6)*u8+diff(X4,u8)*u6)==0);
assert(diff(X5,u1)*u3-diff(X5,u3)*u1-diff(X5,u5)*u7+diff(X5,u7)*u5-( -diff(X5,u2)*u4+diff(X5,u4)*u2-diff(X5,u6)*u8+diff(X5,u8)*u6)==0);

assert(diff(X1,u1)*u4-diff(X1,u4)*u1-diff(X1,u5)*u8+diff(X1,u8)*u5+( -diff(X1,u2)*u3+diff(X1,u3)*u2-diff(X1,u6)*u7+diff(X1,u7)*u6)==0);
assert(diff(X2,u1)*u4-diff(X2,u4)*u1-diff(X2,u5)*u8+diff(X2,u8)*u5+( -diff(X2,u2)*u3+diff(X2,u3)*u2-diff(X2,u6)*u7+diff(X2,u7)*u6)==0);
assert(diff(X3,u1)*u4-diff(X3,u4)*u1-diff(X3,u5)*u8+diff(X3,u8)*u5+( -diff(X3,u2)*u3+diff(X3,u3)*u2-diff(X3,u6)*u7+diff(X3,u7)*u6)==0);
assert(diff(X4,u1)*u4-diff(X4,u4)*u1-diff(X4,u5)*u8+diff(X4,u8)*u5+( -diff(X4,u2)*u3+diff(X4,u3)*u2-diff(X4,u6)*u7+diff(X4,u7)*u6)==0);
assert(diff(X5,u1)*u4-diff(X5,u4)*u1-diff(X5,u5)*u8+diff(X5,u8)*u5+( -diff(X5,u2)*u3+diff(X5,u3)*u2-diff(X5,u6)*u7+diff(X5,u7)*u6)==0);


で、Wikipediaには、もうひとつOctonionic Hopf fibrations S^15->S^8があると書いてある。そうすると、16次元調和振動子から9次元Kepler系への簡約が存在するのか考えたくなる。ファイバーであるところのS^7は、残念ながらコンパクトLie群の構造を持たないので、今までとは違う事情がありそうな感じである。good newsとしては、(9次元Kerpler系のspectrum generating algebraである)so(10,2)がsp(32,R)の部分代数として実現できる(らしい)という事実がある。それで、こういうものも調べた論文があるんじゃないかと思って、ググると、色々見つかったりするもので。

On the SO(10, 2) dynamical symmetry group of the MICZ-Kepler problem in a nine-dimensional space

http://aip.scitation.org/doi/10.1063/1.3606515

Exact analytical solutions of the Schrödinger equation for the nine-dimensional MICZ-Kepler problem

http://aip.scitation.org/doi/abs/10.1063/1.4921171?journalCode=jmp

という論文が出ており、9次元Kepler系を16次元調和振動子に帰着して解けるらしい。ただ、内容自体は簡単な計算をしただけ(端的に言って、殆ど何も新しいことはない)のもので、特に"低次元"のKS変換にいたO(1)=Spin(1),U(1)=Spin(2),Sp(1)=SU(2)=Spin(3)に相当するものが、どうなるのかについて情報がなく、相空間の簡約があるのかどうかは不明


#以下の論文によれば、更に高い次元でも類似の構成が存在するらしい(私は、きちんと確認していないので本当かどうかは知らない)。"一般化KS変換"の作り方は、Clifford代数$Cl_{0,n}$の$2^n$次元表現を利用しているように見える(論文の式(3)には、謎の条件Sp(¥Gamma_{¥lambda})=0というのが付いてるけど、何のことか不明)。論文が正しければ、例えば17次元Kepler系を32次元調和振動子に帰着して解くことができる(so(18,2)がsp(64,R)の部分Lie代数になる?)。しかし、相空間の簡約が存在するのかどうかは分からない。

Theory of the generalized Kustaanheimo-Stiefel transformation

http://www.sciencedirect.com/science/article/pii/037596019390520A


#Hopf fibrationの存在はdivision algebraの存在と関連しているのだと思う(S^15->S^31->S^16というfibrationは存在しないらしい)。一方、一般化KS変換の存在は、もうちょっと広いCayley-Dickson代数というものと関係しているらしい。Cayley-Dickson代数は、通常の実数・複素数・四元数・八元数を含んで、更に16元数(sedenion),etc.と無限に続く。これらの文脈では、Hopf写像と呼ばずに、Hurwitz変換と呼ばれていることが多いっぽい。同じものに、違う名前が一杯ついていて辛い。。とりあえず"Kustaanheimo"のスペルが何回書いても覚えられそうにないので、違う名前を採用して欲しい(まぁ、CFL条件とか、Courant以降誰だっけって感じだけど、CFL条件という名前で覚えているので、KS変換という名前で普及していけば関係ないかもしれない)


#d次元Kepler系のspectrum generating algebraは、D=d+1次元の共形代数で、d=2,3,5,9に対応するのは、D=3,4,6,10である。これらの次元は、SUSY関連ではよく見る。so(10,2)がsp(32,R)の部分Lie環になるというのも、その手の文献に書いてあった(ので、自分では何も確認してない)。これらの次元が重要な理由は、division algebraの存在と関係しているらしい。

super Poincare Li algebra#Lie algebra cohomology

https://ncatlab.org/nlab/show/super+Poincare+Lie+algebra#LieAlgebraCohomology


#KS変換は、(古典力学特に天文周りの文献では)KS正則化という名前で登場することも多い。重力は原点に特異点があり衝突時の計算に問題を引き起こすが、この特異点を何らかの方法で除去することを正則化と呼ぶらしい(量子力学側では、正則化しても表現としては同値になってるので、その観点を強調する意味は薄いかもしれない)。別の方法として、(負エネルギーの)Kepler系を球面上の測地流と"同一視する"というMoserの方法がある(これはKS変換と違って、任意次元で使える)。測地流は量子化すると、ただのラプラシアンである。S^nの共形変換はSO(n+1,1)なので、水素原子(3次元Kepler系)のspectrum generating algebraのうち、so(4,1)までは自然に見える。この方法は量子力学側ではFockによる解法として知られている。Moserの方法で正エネルギー(散乱状態)を扱う場合、hyperboloid上の測地流を扱うことになるようである(量子力学側では、spherical harmonicsの代わりに"hyperbolic harmonics"を考えることになる)

Regularization of kepler's problem and the averaging method on a manifold

http://onlinelibrary.wiley.com/doi/10.1002/cpa.3160230406/abstract


##Moserの正則化/Fockの解法の表現論的理解については、以下の論文などを参照

On the dynamical symmetries of the Kepler problem

http://aip.scitation.org/doi/abs/10.1063/1.524511

Quantization of the Kepler manifold

https://projecteuclid.org/euclid.cmp/1104160353


##so(4,1)の既約表現に実はso(4,2)の作用があるというのは不思議なことで、何でこういうことが起きるのか、よく分からない(似た事例として、Poincare代数のmassless規約ユニタリ表現に共形代数の作用があるというのがある)。Liouville可積分性の立場から見ると、so(4,1)からso(4,2)へ拡大する時、ラプラシアンと交換する作用素が一個増える



最後に、今計算しようという気は全然起きないけど、octonionic Hopf mapの式だけ書いておく。間違えそうなので、Risa/Asirで。

X1 = (u1^2+u2^2+u3^2+u4^2+u5^2+u6^2+u7^2+u8^2)-(u9^2+u10^2+u11^2+u12^2+u13^2+u14^2+u15^2+u16^2);
X2 = 2*(u1*u9-u2*u10-u3*u11-u4*u12-u5*u13-u6*u14-u7*u15-u8*u16);
X3 = 2*(u1*u10+u2*u9+u3*u12-u4*u11+u5*u14-u6*u13-u7*u16+u8*u15);
X4 = 2*(u1*u11-u2*u12+u3*u9+u4*u10+u5*u15+u6*u16-u7*u13-u8*u14);
X5 = 2*(u1*u12+u2*u11-u3*u10+u4*u9+u5*u16-u6*u15+u7*u14-u8*u13);
X6 = 2*(u1*u13-u2*u14-u3*u15-u4*u16+u5*u9+u6*u10+u7*u11+u8*u12);
X7 = 2*(u1*u14+u2*u13-u3*u16+u4*u15-u5*u10+u6*u9-u7*u12+u8*u11);
X8 = 2*(u1*u15+u2*u16+u3*u13-u4*u14-u5*u11+u6*u12+u7*u9-u8*u10);
X9 = 2*(u1*u16-u2*u15+u3*u14+u4*u13-u5*u12-u6*u11+u7*u10+u8*u9);


/*
fctr(X1^2+X2^2+X3^2+X4^2+X5^2+X6^2+X7^2+X8^2+X9^2);
が
[[1,1],[u1^2+u2^2+u3^2+u4^2+u5^2+u6^2+u7^2+u8^2+u9^2+u10^2+u11^2+u12^2+u13^2+u14^2+u15^2+u16^2,2]]
になってれば多分OK
*/