2016-11-17

最近見た面白い論文のコーナー

Special functions and the Mellin transforms of Laguerre and Hermite functions

https://arxiv.org/abs/math-ph/0612085

Hermite多項式とLaguerre多項式のMellin変換は、ゼロ点がRe(s)=1/2上にある(そして、全て単純零点)という事実が、1980年代から知られているらしい。

このMellin変換が超幾何関数の特殊値で表せて(パラメータ(a,b,c)がsの関数として表せる)、関数等式が超幾何関数のKummerの関数等式の帰結として出てくるというのが新しい発見のよう(そして、Kummerの関数等式は、GL(4)のWeyl群作用からの解釈がある)

そしてまた、任意次元の水素原子の固有エネルギー状態のMellin変換(この場合のMellin変換の定義は論文参照)も、Re(s)=1/2上にのみゼロ点を持ち、いずれも単純零点であるとのこと。これの背景には、調和振動子の適当な波動関数($O(n)$対称性のある既約表現に属するという条件がついている)のMellin変換を考えると、やはりRe(s)=1/2上にゼロ点があるという話があるらしい。

A Local Riemann Hypothesis, I

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


水素原子の束縛状態の空間(3次元の場合はso(4,2)の極小表現)は、振動子表現(Weil表現/Segal-Shale-Weil表現)の部分空間としての実現を持つけど、表現の異なる実現に於いて、このような性質を共有するのは意外。KS変換する前と後でMellin変換の結果が同じものになるのかどうかは確認していないけど、もし一致するのであれば、表現の実現の仕方に依存しない形で、Mellin変換が理解できるのか興味がある。わたしは数論には興味がないし、Riemann予想は解けたら何が嬉しいのか、さっぱり分からないので、今まで興味がなかったけど、表現論的に何か面白いことがあるのかもなぁという感じ

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

ただのSegal-Shale-Weil表現(振動子表現という名前もある)の話。無限次元既約ユニタリー表現の研究は、Lorentz群やPoincare群で始まったけど、ある意味では、それより以前から知られていた、最初の無限次元既約ユニタリー表現の例でもある。


水素原子の表現論

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

水素原子の表現論(1.5-2)spectrum generating algebra

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

などで、

・3次元Kepler系の相空間には、群$SO(4,2)$が推移的にsymplectic変換として作用し、symplectic等質多様体であること

・水素原子の束縛状態の空間は、Lie環$so(4,2)$のある既約ユニタリー表現(極小表現)であり、$so(4,2)$はspectrum generating algebraという名前で呼ばれること

を書いた。また、Hamiltonianと交換する保存量として、角運動量とRunge-Lenzベクトルがあり、対応して、古典的にはエネルギー等値面上に、(エネルギーの大きさに応じて)$SO(4)$や$SO(3,1)$などの群が作用している。そして、水素原子に於ける主量子数はLie環$so(4,2)$の既約ユニタリー表現を$so(4)$に制限した時の分岐則として理解できる。というわけで、Kepler系では、相空間の幾何学的対称性とHamiltonianの対称性が、古典力学と量子力学を共に統制している様子を見ることができる。



調和振動子については、Heisenberg群/Heisenberg代数の表現論による記述がよく知られている。調和振動子の相空間$\mathbf{R}^{2n}$はHeisenberg群の余随伴軌道であって、幾何学的量子化が完璧に機能する。一方、(振動数が全て同一の$D$次元)調和振動子が、$su(D)$の(無限小)対称性を持つことも知られている。$su(D)$対称性は、(知らなくてもエネルギー固有値を計算できるし、そもそも一次元しか扱わない場合も多いので)量子力学の教科書には必ずしも書いていないことが多い気がするけど、調和振動子のエネルギーの縮退は、$su(D)$対称性の存在によって説明できる。調和振動子は明らかに回転対称性を持つが、これは、この調和振動子の対称性の一部で、このあたりの事情は、水素原子のエネルギーの縮退が、回転対称性より大きな対称性代数$so(4)$の表現によって説明されるのと似ている。


ところで、調和振動子の対称性$su(D)$はHeisenberg代数の部分Lie環になっていない。これ自体は、珍しいことでもないので、別に構わないという立場もありうるけど、「水素原子の表現論(1.5-2)」で触れたJosephの論文

Minimal realizations and spectrum generating algebras

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

には、symplectic代数$sp(2n,\mathbf{R})$が調和振動子のspectrum generating algebraであると書いてある。シンプレクティック群$Sp(2D,\mathbf{R})$は$\mathbf{R}^{2D}$に推移的には作用してないけど、稠密な開軌道は持つ。


Heisenberg代数とsymplectic群/代数との関係は、symplectic群がWeyl代数の環同型として作用することから来る(Weyl代数が、Euclid空間の余接空間の量子化であることを思い出せば自然)(Weyl代数の自己同型群自体は、もっとずっと大きい。 arXiv:math/0512169などを参照)。調和振動子の状態空間は、Heisenberg代数の既約ユニタリー表現空間であるが、この既約ユニタリー表現は、symplectic代数とHeisenberg代数の半直積の既約ユニタリー表現に一意に拡張でき、Segal-Shale-Weil表現と呼ばれる。この既約ユニタリー表現をsymplectic代数に制限すると、可約になって2つの既約表現に直和分解する。この2つの既約表現は極小表現になっていることが知られている。


注)Segal-Shale-Weil表現は、シンプレクティック群の射影表現にしかならない。量子力学で、幾何学的対称性の群の射影表現を考える必要が最初に生じたのは、スピン表現だと思うけど、事情は似ていて、スピン表現が$SO(3)$の二重被覆$SU(2)$の表現にはなっているように、シンプレクティック群の二重被覆であるメタプレクティック群の表現にはなる。まぁ、群で考えると、計算が面倒になるだけで特にメリットもないので、以下計算する時はLie環の方で考える。



というわけで、調和振動子の状態空間への$sp(2D,\mathbf{R})$の作用や、Hamiltonianの対称性$su(D)$の$sp(2D,\mathbf{R})$の部分Lie代数としての実現を見て、分岐則を調べれば、調和振動子のエネルギー固有値ごとに縮退した状態が、対称性$su(D)$の既約表現として得られるっぽいことを見る。この分岐則は、Segal-Shale-Weil表現の説明に必ず書いてあるようなことではあるし、特に複雑な計算もない(が、よい教科書は知らない)。以下、D>1を仮定する(D=1の時は、$su(1)$は自明になり、エネルギーに縮退もないので)


bosonの生成・消滅演算子と調和振動子のHamiltonianを、通常のように

 ¥[a_i^{-} , a_j^{+}¥] = ¥delta_{ij}

 ¥[a_i^{¥pm} , a_j^{¥pm} ¥] = 0

 H = ¥omega ¥sum_{k=1}^{D} (a_{k}^{+}a_k^{-} + ¥frac{1}{2})

で定める。容易な計算で

 X_{ij} = a_i^{+} a_j^{-}

は全てHamiltonianと可換になることが分かる。また簡単な計算によって

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

なので、これがLie環$gl(D,\mathbf{C})$と同型なLie代数の基底をなすことが分かる。


Hamiltonianの対称性だけ見ていると、これで十分なのだけど、(後でspectrum generating algebraの基底も見ると)これでは(複素化した)spectrum generating algebraの部分Lie環にならないことが分かるので、

 X_{ij} = a_i^{+} a_j^{-} + ¥frac{1}{2} ¥delta_{ij}

とする必要がある。これも、依然としてLie環$gl(D,\mathbf{C})$の基底となるけど、おまけで付いてくる1/2のせいで、この表現は、一般線形群$GL(D,\mathbf{C})$の表現には持ち上がらない(群に拘りたい場合は、一般線形群の二重被覆metalinear群を考えることもできる)。しかし、$sl(D,\mathbf{C})$のCartan部分環は、標準的な基底として

 h_{i} = X_{ii} - X_{i+1,i+1}

を取ることができ、うまいこと1/2が相殺してくれるので、特殊線形群の表現には持ち上がる。$sl(D,\mathbf{C})$に入らない残りの一次元分の対称性は

 ¥sum_{k=1}^D X_{ii} = ¥sum_{k=1}^D (a_k^{+} a_k^{-} + ¥frac{1}{2})

で、定数倍を除けばHamiltonianそのものである。


#群に表現が持ち上がるかどうかという観点では、$SL(D,\mathbf{C})$や$SU(D)$は単連結だが、$GL(1)$や$U(1)$はそうでないというあたりから、面倒さがやってくる。$D$変数多項式環への$gl(D,\mathbf{C})$の作用として、行列単位に対して

 E_{ij} = x_i ¥partial_j + s ¥delta_{ij}

を取れる。対応する一般線形群$GL(D,\mathbf{C})$の表現は形式的には

 ¥pi_{s}(A)f(¥mathbf{x}) = ¥det(A)^{s} f(¥mathbf{x} A)

であるが、$s$が非整数の時、群の表現の方は$(det A)^{s}$が一価でない。しかし、$det A=1$の時に限定すれば問題はなくなり、特殊線形群$SL(D,\mathbf{C})$の表現は定まる。


spectrum generating algebraの実現。まず、複素化した$sp(2D,\mathbf{C})$を考えたほうがわかりやすい。次のように基底を取る

A_{ij} = ¥begin{pmatrix} E_{ij} & O ¥¥ O & -E_{ji} ¥end{pmatrix}

B_{ij} = ¥begin{pmatrix} O & E_{ij}+E_{ji} ¥¥ O & O ¥end{pmatrix}

C_{ij} = ¥begin{pmatrix} O & O ¥¥ E_{ij}+E_{ji} & O ¥end{pmatrix}

対応する演算子をbosonの生成・消滅演算子を使って

X_{ij} = a_i^{+} a_j^{-} + ¥frac{1}{2} ¥delta_{ij}

Y_{ij} = a_i^{+} a_j^{+}

Z_{ij} = a_i^{-} a_j^{-}

とする。これが、symplectic代数と同型なLie環をなすことは計算するのみ。Heisenberg代数にderivationとして作用していることも簡単な計算で分かる。Cartan対合$s$を

s(A_{ij}) = -A_{ji}

s(B_{ij}) = C_{ij}

で定める。このCartan対合に関する実型は$sp(2D,\mathbf{R})$と同型(証明略)。特に、上で作ったLie環$gl(D,\mathbf{C})$は部分Lie環になってることが分かり、このLie環の元のうち、Cartan対合で不変なもの全体は実Lie環$u(D)$をなす。


以上の定義によって、調和振動子の状態空間への$sp(2D,\mathbf{C})$と$sl(D,\mathbf{C})$(あるいは$sp(2D,\mathbf{R})$と$su(D)$)の作用が定まった(ユニタリー表現になることなどは示されてないけど)。symplectic代数の作用を見ると、振動子の個数を変化させないか、2個増減させることが分かるので、少なくとも振動子の個数の偶奇によって、2つの表現に直和分解する。


調和振動子の状態空間を、$su(D)$の表現として既約分解する。試みとして、|0>を基底状態として

 v_i = a_{i}^{+} |0 ¥rangle

を基底とする、ベクトル空間(振動子が一個の状態全体)を考える。そうすると

 a_{i}^{+} a_{j}^{-} (v_{k}) = ¥delta_{jk} v_{i}

なので、$su(D)$が自然表現として作用していることが分かる。一般に振動子$N$個の状態空間は、$su(D)$の自然表現の$N$次対称テンソル積表現となっていて、これらは全て既約表現である($N$次斉次多項式の空間と同一視できるので、適当に基底を取って、最高ウェイトetc.を直接計算すれば示せる。また、この事実から、上で書いた$sp(2D,\mathbf{R})$の2つの表現が既約であることが従う)。$su(D)$を複素化すると$sl(D,\mathbf{C})$なので、振動子の状態空間は、$sl(D,\mathbf{C})$の自然表現の対称テンソル積代数(必要なら、適当に完備化したもの)と同定できる。


Liouville可積分性について。調和振動子の場合、$su(D)$のCartan部分環は、$D-1$個の独立でinvolutiveな保存量を与える。Hamiltonian自身を合わせて、相空間の次元の半分の個数のinvolutieな保存量が得られる。各振動子の振動数ωが異なっていても、Liouville可積分ではあるが、対称性は、$su(D)$の部分代数に落ちていく。最も退化した場合は、$u(1)^{D-1}$となる


Heisenberg vs symplectic。量子Kepler系/水素原子との関連で考えると、Heisenberg代数でなく、symplectic代数をspectrum generating algebraとした方が自然ではある。水素原子については、調和振動子の$U(1)$作用によるHamiltonian reductionとして得られる(散乱状態の場合は、"虚の振動数"を持つHamiltonianを考える)けど、

M = ¥sum_{i=1}^2 X_{ii} - ¥sum_{i=3}^{4} X_{ii}

と可換な$sp(8,\mathbf{R})$の部分Lie代数が、$u(2,2)$と同型であるということが、Kepler系のspectrum generating algebraが$su(2,2)$であることの、1つの説明を与える。この$M$は$U(1)$作用の量子化運動量写像を与えると解釈でき、$u(1)$部分は、Hamiltonian reductionを取る際に0に取るので、結局$su(2,2)$が残る


#古典的な場合、Kustaanheimo-Stiefel変換と、このHamiltonian reductionの関係は

On the regularization of the Kepler problem

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

の3章に簡潔にまとまっている。


#また、Kepler系の相空間は、$su(2,2)$の極小軌道であった。一般に(symplectic群の)dual pairがあった場合、pairの片割れによる簡約相空間は、残りのpairの余随伴軌道(の閉包)と同型になることは、漠然と予想されているっぽい。が、この対応が、いつ正しいのか、正確なところは、よく分からない。

Examples of singular reduction

http://www.math.cornell.edu/~sjamaar/papers/lms.pdf

の4節(特に定理4.3)には、かなり限定的な条件ではあるが、dual pairの簡約相空間と余随伴軌道の対応について、証明がある(1990年頃のものなので、今ではもっと一般的な事実が知られていてもよさそうではある)。Kepler系の相空間が、極小軌道となることも、こうした事実の一類型となっている。



余談。Segal-Shale-Weil表現は、symplectic代数の表現としては既約表現になってない点に不満がある。Heisenberg代数とsymplectic代数の半直積を考えるという道もありそうではある。群としては、シンプレクティック群とHeisenberg群の半直積を考えることになり、この群はJacobi群という名前が付いている。Jacobi群という名前は

Jacobi group

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

によると、EichlerとZagierによる命名らしい。彼らはJacobi形式の研究のために、使ったとか。物理の方では、squeezed stateをJacobi群のgeneralized coherent statesとして理解できる、と書いている人がいる。私は、量子光学を知らないので、squeezed状態が何か知らないし、これが妥当なのかも確認していない(Segal-Shale-Weil表現は、Jacobi群の表現としては当然射影表現にしかならない)

Generalized squeezed states for the Jacobi group

https://arxiv.org/abs/0812.0717

The Jacobi group and the squeezed states - some comments

https://arxiv.org/abs/0910.5563



余談その2。私は全然興味がないけど、テータ関数の変換公式をWeil表現から導出するという話があったりするらしい

リー群の表現から見たテータ関数

http://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/48222

で、主に、$SL(2,R)$の場合に絞って考え方を理解しやすく解説してくれている


調和振動子は量子力学の模型としては簡単で面白くもない部類のものだけど、振動子表現は、表現論的になかなか面白いものである。


その他の初等的な可積分系について。

まず、(massiveな)自由粒子。自由粒子は教科書の最初に出てくる最も簡単な系にも関わらず、表現論的に理解するのは、それほど自明ではない。明らかに重要な代数として、(非相対論的な場合には)ガリレイ代数、(相対論的な場合には)ポアンカレ代数がいる。こいつらは、Hamiltonianと可換でない演算子(Galilean boost,Lorentz boost)を含み、対称性代数をこれらの部分代数と考えることができる。一方、Wignerの仕事によって、粒子状態の空間は、ガリレイ代数/ポアンカレ代数の既約ユニタリー表現と理解できるはずで、これを対称性代数に制限した時の既約分解を考えると、エネルギー固有状態に対応する表現が得られると予想される。エネルギースペクトルが離散的でないことから、分岐則は、直和の形ではなく直積分の形になると思われる。


#ガリレイ代数は、粒子の質量と関連する中心拡大を持つ。相空間に作用する幾何学的変換としては、この項は見えない。波動関数の変換としては、位相の回転と関連し、ガリレイ群の射影表現を考える必要があることの表れでもある。まぁ、群で考えると、計算が面倒になるだけで特にメリットもないのでry


#masslessな自由粒子としては、光子が代表的であるけども、これについては

Maxwell方程式の表現論(1)(2)

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

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

で書いたものに相当する。これらは量子力学系に相当するけど、対応する古典論は"幾何光学"と呼ばれるべきものである。通常の意味の幾何光学と同一なのか知らないけど、正準理論の起源が幾何光学にあった(らしい)ことを思うと、面白い。またmasslessなスカラー粒子は、素粒子として存在するのか知らないけど、対応する既約ユニタリー表現は、水素原子の束縛状態全体の空間と同じである。masslessスカラー粒子の場合と、水素原子の場合では、Hamiltonianは異なるので、対称性と分岐則も異なる


個人的に、よく分からない例として、自由剛体がある。可積分なコマとしては、Euler top,Lagrange top,Kovalevski topがあるけど、特にEuler topを理解したい。この系はHamiltonian以外に1つの保存量を持つLiouville可積分系であり、古典的には、楕円関数を使って解を書ける。量子論の場合、明らかに、$so(3)$の既約ユニタリー表現空間に作用し、それは全て有限次元空間なので、Hamiltonianの行列要素を具体的に書くことができ、原理的には解ける。とはいえ、有限次元空間に作用しているから解けますというのでは、Heisenberg模型(XXX模型)なんかも自明ということになってしまう。

Lamé equation, quantum top and elliptic Bernoulli polynomials

https://arxiv.org/abs/math-ph/0508068

なんかは、quantum Euler topの面白さと非自明さが現れているのでないかと思う。Euler topは、調和振動子や量子Kepler系のように、Lie環の表現の分岐則を見れば殆どの情報が得られるというような、単純な話にはなっていない。

また、以下のようなRisa/Asirコード書いて計算すると、半整数スピン(Nが偶数)の時は、ハミルトニアンの固有多項式がF(x)^2の形に因数分解できるらしい。つまり、理由はわからないけど、半整数スピンの時は、全ての固有状態が二重に縮退している。不思議なことである

N=5;  /*  2s+1  */

Emat = newmat(N,N);
Fmat = newmat(N,N);
Hmat = newmat(N,N);

Umat = newmat(N,N);

for(I=0;I<N-1;I++){
   Emat[I][I+1] = I+1;
   Fmat[I+1][I] = N-I-1;
   Hmat[I][I] = (N-1-2*I);
   Umat[I][I] = 1;
}
Hmat[N-1][N-1] = -(N-1);
Umat[N-1][N-1] = 1;

S1 = Emat+Fmat;
S2 = @i*(-Emat+Fmat);
S3 = Hmat;

H = a*S1*S1+b*S2*S2+c*S3*S3;

fctr(det(H-4*x*Umat));

(量子力学の入門的な教科書には載っていないという意味で)初等的でない例として、rational Calogero-Moser模型に興味がある(trigonometricなSutherland模型も殆ど同じことになると思うけど)。この模型は、一次元空間上の相互作用する多体粒子系と説明される。表現論的には、粒子が何次元空間に住んでるとか、多体系か否かは副次的な解釈の問題に過ぎないので、どうでもいいが。1970年代はじめに、量子可積分系として発見され、1978年には、Kazhdan-Kostant-Sternbergによって、Calogero-Moser spaceという"真の相空間"が構成されている。量子系の方は、Dunkl作用素という、調和振動子に於ける生成・消滅演算子のようなものを導入して解かれる。Dunkl operatorは1980年代終わり頃発見されたものだけど、90年代以降も、非常に多くの研究があり、多くの一般化や様々な代数との関連が議論されている。

The Quantum Calogero-Moser Model: Algebraic Structures

https://www.jstage.jst.go.jp/article/jpsj1946/62/9/62_9_3035/_article

Integrability and algebraic structure of the quantum Calogero-Moser model

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

Exact Solvability of the Calogero and Sutherland Models

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

Lectures on Calogero-Moser systems

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

Dunkl operators: Theory and applications

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

SBMLによる生化学ネットワークの記述と計算

SBMLは、Systems Biology Markup Languageの略。シグナル伝達や代謝系などの生化学ネットワークのモデルを形式的に記述するXML仕様で、2000年頃から継続的に策定と改訂が行われているらしい。生物分野で古典的な数理モデルといえば、酵素反応速度論(Michaelis-Menten式etc.)、Hodgkin-Huxley方程式、Turingの反応拡散系などがあるけど、前の2つは、SBMLで記述されるモデルのプロトタイプと言っていいと思う。これらのモデルは、数個の変数しかない単純なモデルだけど、最近では細胞内の生化学ネットワークを調整する遺伝子たちが沢山同定されて、複雑な生化学ネットワークを、そのままモデル化して調べるという道ができた。


その極致として、2012年に、500個の遺伝子を持つ細菌マイコプラズマの全細胞シミュレーションが行われている

A whole-cell computational model predicts phenotype from genotype.

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

いずれ、線虫(細胞数1000、遺伝子数20000弱)くらいの生物で全細胞シミュレーションができれば面白いと思うけど、全ての要素をモデル化するだけでも大変そう(マイコプラズマの場合でも、シミュレーションの結果から新しい酵素が見つかったとかあったようだし)


#関連する話題として、最近、RIKENで同じくマイコプラズマの全原子分子動力学シミュレーションを行ったという報告があった

バクテリア細胞質の全原子分子動力学計算(プレスリリース)

http://www.riken.jp/pr/press/2016/20161101_3/

こっちは、20nsを数カ月かけて計算するというものらしいので、計算量的に何かへ応用するというには、難しいけれど。


More Detailed Summary of SBML

http://sbml.org/More_Detailed_Summary_of_SBML

に、もう少し詳しいSBMLの説明と、ミカエリス・メンテン型酵素反応の例が載っている。


既に知られているモデルの場合は、SBMLを

BioModels Database

http://www.ebi.ac.uk/biomodels-main/

とかから見つけてこれるかもしれない。Curatedでないモデルは色々ダメなことが多いっぽいので、使えるモデルは600ちょっとと思っておいたほうがよさげではある。まぁ、この手のシミュレータを実装すること自体は、難しいことでもないと思うけど、何気にパラメータが多くて、論文などから適切な値を探すのが面倒だったりするので、モデルをまとめて公開してくれているとありがたい。


SBMLのシミュレータはいくつも公開されている。

SBML Software Matrix

http://sbml.org/SBML_Software_Guide/SBML_Software_Matrix

にソフトウェアのリストがある。libSBMLとかあるけど、こいつにはシミューレション機能とかはなく、SBMLファイルの読み書きをするだけで、割とどうでもいい。対応しているSBMLのバージョン、実装されている機能や計算速度、継続的にメンテナンスがされているのかなどを考慮して決めたいのだけど、いまいちよく分からないので、自分で実装した方が早そうであった



最も単純な利用法は、常微分方程式系に変換してソルバーに投げるというのだと思う。細胞は非常に小さく、分子種によっては分子数が非常に少ない場合があり、ゆらぎが無視できない場合もあると考えられており、その場合は確率モデルを用いる必要があると考えられている。その他にも、目的に応じて、適用したい解析手法は多数あると思う。

自分でシミュレータを実装する場合は

SBML Test Suite

http://sbml.org/Software/SBML_Test_Suite

から、テストデータやらを引っ張ってこれるっぽいので、確認用に使える。



例として、カルシウム振動のモデル

http://www.ebi.ac.uk/compneur-srv/biomodels-main/BIOMD0000000184

を使う。アストロサイトにおけるカルシウム振動の(多分最も単純な)モデル。変数としては細胞質のカルシウム濃度とIP3(イノシトール3リン酸)濃度と、ERのカルシウムの濃度。

(1)細胞外から細胞内へのCa2+の流入及び流出

(2)ERから細胞質へのカルシウムの放出(IP3受容体を介するもの)

(3)細胞質からERへのカルシウムの浸透

(4)ERから細胞質へのカルシウムの漏れ

(5)IP3の合成(カルシウム濃度依存。ホスホリパーゼCの活性がカルシウムによる調節を受けるらしい)

(6)IP3の分解

が考慮される。3つの変数が振動を起こすことが確認できる。計算結果では、振動周期は3分程度である。これが、カルシウム振動のモデルとして適切なのかどうかは不明


とりあえず、これくらいのモデルがPython(+numpy/scipy/matplotlib)で動くようにした

・多分、2.7と3.5の両方で動く

・Test Suiteは通してないし、仕様の確認もしてないので、仕様を網羅もしてないと思う

・凄く遅い


いうて、常微分方程式やし、簡単にできると思ったが、色々面倒くさいことがわかった。とりあえずミカエリス・メンテン型反応とカルシウム振動の計算結果を貼っておく

f:id:m-a-o:20161117232853p:image:w360

f:id:m-a-o:20161117232851p:image:w360

from scipy.integrate import ode
import xml.etree.ElementTree as et
import math
from numpy import prod


def sbml_root_fn(xs):
    if len(xs)==1:
        return math.sqrt(xs[0])
    elif len(xs)==2:
        return math.pow(xs[1] , 1.0/xs[0])
    else:
        assert(False),"invalid argument in sbml_root_fn:{}".format(str(xs))

class Closure:
    def __init__(self , _vars , fnbody):
        self.vars = _vars
        self.fndef = fnbody
        self.env = {}
    def __call__(self , args):
        assert(len(args)==len(self.vars)),"invalid closure call"
        saved_env = {}
        for (key,val) in zip(self.vars , args):
            if key in self.env:
                saved_env[key] = self.env[key]
            self.env[key] = val
        ret = eval_math(self.fndef , self.env)
        for key in self.vars:
            del self.env[key]
        for (key,val) in saved_env.items():
            self.env[key] = val
        return ret


def parse_math(e):
    def parse_cn(e):
       assert(e.attrib.get('type',None)!="rational"),"rational cn@parse_cn"
       return eval(e.text)
    def parse_apply(es):
       fn = None
       args = []
       for e0 in es:
          _ , tagname  = e0.tag.split('}')
          if tagname=="plus":
              fn = (lambda x:x[0]+x[1])
          elif tagname=="minus":
              fn = (lambda x:x[0]-x[1])
          elif tagname=="divide":
              fn = (lambda x:x[0]/x[1])
          elif tagname=="times":
              fn = (lambda x:prod(x))
          elif tagname=="power":
              fn = (lambda x:math.pow(x[0],x[1]))
	  elif tagname=="exp":
	      fn = (lambda x:math.exp(x[0]))
          elif tagname=="root":
              fn = sbml_root_fn
          elif tagname=="apply":
              args.append( parse_math(e0) )
          elif tagname=="degree":
              assert(fn==sbml_root_fn)
              assert(len(e0.getchildren())==1),"error in parse_apply@degree"
              e1 = e0.getchildren()[0]
              _ , tagname1  = e1.tag.split('}')
              if tagname1=="ci":
                 args.append( e1.text.strip() )
              elif tagname1=="cn":
                 assert(e0.attrib.get('type',None)!="rational"),"rational cn@parse_math"
                 args.append( eval(e1.text) )
              else:
                 assert(False),"failed to parse degree"
          elif tagname=="ci":
              if fn!=None:
                  args.append( e0.text.strip() )
              else:
                  fn = e0.text.strip()
          elif tagname=="cn":
              assert(e0.attrib.get('type',None)!="rational"),"rational cn@parse_math"
              args.append( eval(e0.text) )
          else:
              assert(False),"@parse_math:unknown tag {}".format(tagname)
       assert(fn!=None),et.tostring(es)
       return (fn,args)
    def parse_lambda(es):
       varnames = []
       for e0 in es:
          _ , tagname  = e0.tag.split('}')
          if tagname=="bvar":
              for it in e0.getchildren():
                  varnames.append( it.text.strip() )
          elif tagname=="apply":
              fnbody = parse_apply(e0.getchildren())
              return Closure(varnames , fnbody)
          else:
              assert(False),"unknown tag name@parse_lambda:{}".format(tagname)
    _ , tagname  = e.tag.split('}')
    if tagname=="math":
         assert(len(e.getchildren())==1),"error in parse_math"
         return parse_math(e.getchildren()[0])
    elif tagname=="lambda":
         return parse_lambda(e.getchildren())
    elif tagname=="apply":
         return parse_apply(e.getchildren())
    elif tagname=="cn":
         return parse_cn(e)
    else:
         assert(False),"unknown tagname@parse_math:{}".format(tagname)


def eval_math(e , env):
   if type(e)==tuple:
      fn,_args = e
      args = [eval_math(x,env) for x in _args]
      if type(fn)!=str:
          return fn(args)
      else:
          closure = env[fn]
          closure.env = env
          return closure(args)
   elif type(e)==str:
      if e in env:
          return env[e]
      else:
          assert(False),"Unknown variable in env:{}".format(e)
   else:
      return e



class Model:
   def __init__(self):
      self.args = []
      self.species = []
      self.namemap = {}
      self.env = {}
      self.reactions = []
      self.compartment = {}
      self.rules = {}
   def values(self):
      ret = []
      for key in self.args:
          v = self.env[key]
	  if type(v)==float or type(v)==int:
             ret.append( self.env[key] )
          else:
	     ret.append( eval_math(v ,self.env) )
      return ret
   def __call__(self , t , y):
      assert(len(y)==len(self.args)),"invalid input"
      ret = [0.0 for _ in y]
      for (name,val) in zip(self.args , y):
          self.env[name] = val
      rest = self.rules.keys()
      Nmax = len(rest)
      #-- process assignment rules
      for _ in range(Nmax):
          if len(rest)==0:break
          for (key,sm) in self.rules.items():
	      if key in rest:
	         try:
	            self.env[key] = eval_math(sm , self.env)
		    rest.remove(key)
		 except:
		    pass
      #-- process reactions and rate rules
      for (reactants , products , sm , local_env) in self.reactions:
          saved_env = {}
          for key,val in local_env.items():
              if key in self.env:
                 saved_env[key] = self.env[key]
              if type(val)==str:
                 self.env[key] = self.env[val]
              else:
                 self.env[key] = val
          try:
              v = eval_math(sm , self.env)
          except:
              assert(False),("failed to evaluate math expresion",sm)
          for r in reactants:
	      try:
                 idx = self.args.index(r)
	         if r in self.compartment:
                     ret[idx] -= v/self.env[self.compartment[r]]
	         else:
	             ret[idx] -= v
	      except:
	         pass
          for r in products:
	      try:
                 idx = self.args.index(r)
	         if r in self.compartment:
                     ret[idx] += v/self.env[self.compartment[r]]
	         else:
	             ret[idx] += v
	      except:
	         pass
          for key in local_env:
              del self.env[key]
          for key,val in saved_env.items():
              self.env[key] = val
      return ret


def readsbml(root):
   ns = {'sbml' : root.tag[1:-5] , 'mml':"http://www.w3.org/1998/Math/MathML"}
   ret = Model()
   for it in root.find('.//sbml:listOfCompartments' , namespaces=ns).getchildren():
       key,value = it.attrib['id'] ,float( it.attrib['size'] )
       ret.env[key] = value
   for it in root.find('.//sbml:listOfSpecies' , namespaces=ns).getchildren():
       if it.attrib.get('boundaryCondition' , 'false')!='true':
           ret.args.append( it.attrib['id'] )
       if 'name' in it.attrib:
           ret.species.append( it.attrib['name'] )
	   ret.namemap[ it.attrib['id'] ] = it.attrib['name']
       else:
           ret.species.append( it.attrib['id'] )
	   ret.namemap[ it.attrib['id'] ] = it.attrib['id']
       ret.compartment[ it.attrib['id'] ] = it.attrib['compartment']
       if "initialAmount" in it.attrib:
            vol = ret.env[it.attrib["compartment"]]
            key,val = it.attrib['id'] , float(it.attrib['initialAmount'])/vol
            ret.env[key] = val
       elif "initialConcentration" in it.attrib:
            vol = ret.env[it.attrib["compartment"]]
            key,val = it.attrib['id'] ,  float(it.attrib['initialConcentration'])
            ret.env[key] = val
       else:
            pass
   for it in root.findall('.//sbml:listOfParameters' , namespaces=ns):
       for r in it.getchildren():
           if 'value' in r.attrib:
               ret.env[ r.attrib['id'] ] = eval(r.attrib['value'])
   for it in root.findall('.//sbml:functionDefinition' , namespaces=ns):
       fndef = it.find('.//mml:math' , namespaces=ns)
       key = it.attrib['id']
       if fndef!=None:
           fn = parse_math(fndef)
           ret.env[key] = fn
   listOfInitialAssignments = root.find('.//sbml:listOfInitialAssignments' ,  namespaces=ns)
   if listOfInitialAssignments!=None:
       for it in listOfInitialAssignments.getchildren():
           varname = it.attrib['symbol']
	   ret.env[varname] = eval_math(parse_math(it.getchildren()[0]) , ret.env)
   listOfRules = root.find('.//sbml:listOfRules' , namespaces=ns)
   if listOfRules!=None:
       for it in listOfRules.getchildren():
           _ , tagname  = it.tag.split('}')
           key = it.attrib['variable']
	   if tagname=="assignmentRule":
	       rule = parse_math(it.getchildren()[0])
	       ret.rules[key] = rule
	   else:
	       ret.args.append( key )
	       ret.reactions.append( ([] , [key] , parse_math(it.getchildren()[0]) ,{}) )
   for r in root.find('.//sbml:listOfReactions' , namespaces=ns).getchildren():
       reactants = []
       products = []
       local_env = {}
       listOfReactants = r.find('.//sbml:listOfReactants' , namespaces=ns)
       if listOfReactants!=None:
           for it in listOfReactants.getchildren():
               reactants.append( it.attrib['species'] )
               if 'id' in it.attrib:
                    local_env[ it.attrib['id'] ] = it.attrib['species']
       listOfProducts = r.find('.//sbml:listOfProducts' , namespaces=ns)
       if listOfProducts!=None:
           for it in listOfProducts.getchildren():
               products.append( it.attrib['species'] )
               if 'id' in it.attrib:
                    local_env[ it.attrib['id'] ] = it.attrib['species']
       listOfModifiers = r.find('.//sbml:listOfModifiers' , namespaces=ns)
       if listOfModifiers!=None:
           for it in listOfModifiers.getchildren():
               if 'id' in it.attrib:
                    local_env[ it.attrib['id'] ] = it.attrib['species']
       listOfLocalParameters = r.find('.//sbml:listOfLocalParameters' , namespaces=ns)
       if listOfLocalParameters!=None:
           for it in listOfLocalParameters.getchildren():
               local_env[ it.attrib['id'] ] = eval(it.attrib['value'])
       sm = r.find('.//mml:math' , namespaces=ns)
       ret.reactions.append( (reactants , products , parse_math(sm) , local_env) )
   return ret


try:  #-- python2
   from urllib2 import urlopen
except: #-- python3
   from urllib.request import urlopen


import matplotlib.pyplot as plt


def runsbml(path , t0 , t1, dt):
   try:
      root = et.parse(path).getroot()
   except:
      req = urlopen(path)
      s = req.read()
      root = et.fromstring(s)
   m = readsbml(root)
   r = ode(m, None).set_integrator('vode', method='bdf', nsteps=1500 , with_jacobian=False)
   r.set_initial_value(m.values() , t0)
   res = []
   times = []
   while r.successful() and r.t < t1:
      _ = r.integrate(r.t+dt)
      times.append( r.t )
      res.append( r.y )
   for n,name in enumerate(m.args):
      dat = [x[n] for x in res]
      plt.plot(times , dat , label=m.namemap[name])
   plt.legend()
   plt.show()


if __name__=="__main__":
 #-- calcium oscillation
  runsbml("http://www.ebi.ac.uk/compneur-srv/biomodels-main/download?mid=BIOMD0000000184" , 0.0 , 2000.0 , 0.5)

2016-10-10

水素原子の表現論(1.5-2)spectrum generating algebra

前回の補足など。


前回、極小表現上で消えるU(so(4,2))の元全体をJosephイデアルと呼んでしまったけど、Josephイデアルの(数学者が一般的に採用する)定義は、複素単純Lie代数の普遍展開環のcompletely prime primitive idealであって、so(4,2)は複素単純ではないし、Josephイデアルは、A型"以外の"複素単純Lie環では一意に決まるが、A型の複素単純Lie環ではwell-definedではない。so(4,2)は複素化すると、so(6,C)であるが、これはsl(4,C)と同型なので(実形で見れば、so(4,2)とsu(2,2)の同型がある)、U(so(6,C))のJosephイデアルは、上の意味では定義されないことになる。

Josephイデアルのこの定義は、Josephのオリジナルの論文で採用されている定義でもあり、

The minimal orbit in a simple Lie algebra and its associated maximal ideal

https://eudml.org/doc/81975

前回の最後に書いたGarfinkleの論文の結果

A new construction of the Joseph ideal

https://dspace.mit.edu/handle/1721.1/15620#files-area

も、今回の問題では使えない。


#結局今回は不要になったので、ちゃんと読んでないのだけど、Josephイデアルの一意性証明にはgapがあると、以下の論文に書いてある(新しい証明も)

UNIQUENESS OF JOSEPH IDEAL

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.153.6344&rank=1


いつからminimal representationという名前が一般的になったのか分からないけど、極小表現の"極小"という名前の由来は、極小冪零軌道の量子化と見なせる/見なしたいという気持ちから来ているのだと思う(Gelfand-Kirillov次元が最小という意味もあるかもしれない)。このへんの問題意識と事情は、Josephイデアルの論文のIntroductionにも書かれている。タイトルの"minimal orbit"=極小軌道というのは、(自明な軌道以外で)次元が最小という意味で、A型以外の場合は、半単純軌道も含めて一つしかなく、極小冪零軌道と同じもので、A型の場合は、冪零軌道の中で次元が最小のもの=極小冪零軌道が一意に存在するが、極小冪零軌道と次元が等しい半単純軌道が無数に存在する。


現在、A型以外の複素単純Lie環では、annihilatorがJosephイデアルであることを以て、極小表現を定義することも、しばしばある。この定義だと、A型の場合は、極小表現が定義されないことになってしまうが、A型の場合もこれが極小表現だというのはあって、ただ、統一的かつ自然に見える定義というのがない(極小冪零軌道の量子化とは何かというのが、きちんと定義されて、それで定義するというのが王道だろうとは思う)

The Uniqueness of the Joseph Ideal for the Classical Groups

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

という論文では、(A_1以外の場合に)Josephイデアルを定義している(Josephイデアルが定義できれば、そのannihilatorが極小表現であるという定義を、拡張できる)。生成元を明示的に与えるということをやってのけていて、スマートな定義とはいいがたいけど、計算には使いやすい。極小冪零軌道の変形量子化を考えている人もいて、その場合も、$A_1$型を除く$sl(n,C)$で変形量子化は一意であるらしい(例えば、arXiv:math/0010257)


#極小表現が極小冪零軌道の量子化という時は、"幾何学的量子化"と同じ意味で、極小表現は量子力学の波動関数の空間に相当する。極小冪零軌道の変形量子化を考えるときは、与えられるのは、物理的な言い方では、波動関数ではなく、物理量の代数の方になる


水素原子の束縛状態を考えている今は、具体的な表現は分かっていて、そのannihilatorの生成元を知りたいだけなので、Josephイデアルとか極小表現の定義は重要ではないけれど、用語としては、そんな感じで、ちょっと揺れがある。



本題。色々論文を眺めていたら、Josephの論文

Minimal realizations and spectrum generating algebras

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

は、Josephイデアルや極小表現と深く関連する、"量子化運動量写像"(今の場合、局所化したWeyl代数への環準同型写像)を構成していて面白かった。abstractには、"与えられた半単純Lie環が入るような自由度最小の量子力学系の決定"云々と書かれていて、タイトルの"minimal realizations"とは、このことを指す。


この論文の5章に、「運動量空間に於ける水素原子のspectrum generating algebraの実現」というものが与えられている。論文通り計算すると、Chevalley生成元は、以下の式で与えられる

$e_1 = q_1$

$e_2 = q_2 p_1$

$e_3 = q_3 p_2$

$h_1 = -2 q_1 p_1 - q_2 p_2 - q_3 p_3$

$h_2 = q_1 p_1 - q_2 p_2$

$h_3 = q_2 p_2 - q_3 p_3$

$f_1 = -(q_1 p_1 + q_2 p_2 + q_3 p_3)p_1$

$f_2 = q_1 p_2$

$f_3 = q_2 p_3$

論文に書いてある通り、これによって、sl(4,C)の(3変数)多項式環への作用が与えられるけども、極小表現は、この多項式環の中で実現されるわけではない。


上記の量子化運動量写像(Lie環の普遍展開環からWeyl代数への環準同型。普通の運動量写像がLie環の対称代数からPoisson代数へのPoisson準同型を与えることの類似)には、当然、対応して、古典的な運動量写像も存在し、その像は極小冪零軌道の閉包になってるのじゃないかと期待したい。次元を数えると、sl(3,C)の極小冪零軌道は4次元であり、相空間の方の次元(=Weyl代数の生成元の数)も4だから一致している。同様に、$sl(4,C)$の極小冪零軌道は6次元であり、相空間の方の次元も6で一致する。これらは、次元最小の軌道ではあるが、sl(n,C)の場合は、半単純な次元最小の軌道(極小軌道)が存在するので、次元の勘定だけで、像が極小冪零軌道であるとは言えない。


しかし、変数$q_i,p_j$を消去してやれば、像(の閉包)の定義方程式は得ることができるし、極小冪零軌道の定義方程式も群作用をexplicitに書いて、変数消去すれば直接計算で得られる。闇雲に変数を消去しても、同じ定義式が出るとは限らないので、被約グレブナー基底を計算して比較するのがいい。グレブナー基底による変数消去計算については、グレブナー基底の演習問題


#冪零軌道の次元については、例えば、$sl(4)$の場合、Collingwood,McGovern "Nilpotent Orbits In Semisimple Lie Algebra"の4.3節などに、そのまま書いてある


まず、Chevalley生成元以外の基底、他のCartan-Weyl basisを、以下のように

$e_4 = [e_1 , e_2]$

$e_5 = [e_2 , e_3]$

$e_6 = [e_1 , e_5]$

$f_4 = [f_1 , f_2]$

$f_5 = [f_2 , f_3]$

$f_6 = [f_1 , f_5]$

で定義する。それぞれ対応するCartan部分環の元は

$h_4 = [e_4 , f_4] = -(h_1+h_2)$

$h_5 = [e_5 , f_5] = -(h_2+h_3)$

$h_6 = [e_6 , f_6] = h_1+h_2+h_3$

のようになる。


これらの準備の下で、以下のRisa/Asirコードによって、$sl(3)$と$sl(4)$のそれぞれの場合で、運動量写像の像(の閉包)の定義方程式を具体的に得ることができる(この計算は、手元の貧弱なノートPCでも一瞬で終わる)

/* Josephの運動量写像のkernelの計算 */

load("noro_pd.rr");
/* sl(3)-case */
B = [e1-q1 , e2-q2*p1 , e3-q2 , h1+2*q1*p1+q2*p2 , h2-q1*p1+q2*p2 , 
     f1+q1*p1*p1+q2*p1*p2 , f2-q1*p2 , f3+(q1*p1+q2*p2)*p2];
V = [q1,q2,p1,p2,e1,e2,e3,h1,h2,f1,f2,f3];
G = nd_gr_trace(B,V,1,1,[[0,4],[0,8]]);
D = noro_pd.elimination(G,[e1,e2,e3,h1,h2,f1,f2,f3]);


/* sl(4)-case */
B = [e1-q1 , e2-q2*p1 , e3-q3*p2 , e4-q2 , e5-q3*p1 , e6-q3 ,
     h1+2*q1*p1+q2*p2+q3*p3 , h2-q1*p1+q2*p2 , h3-q2*p2+q3*p3,
     f1+(q1*p1+q2*p2+q3*p3)*p1,f2-q1*p2,f3-q2*p3 , f4-(q1*p1+q2*p2+q3*p3)*p2 , f5+p3*q1 , 
     f6+(q1*p1+q2*p2+q3*p3)*p3];
V = [q1,q2,q3,p1,p2,p3,e1,e2,e3,e4,e5,e6,h1,h2,h3,f1,f2,f3,f4,f5,f6];
G = nd_gr_trace(B,V,1,1,[[0,6],[0,15]]);
D = noro_pd.elimination(G,[e1,e2,e3,e4,e5,e6,h1,h2,h3,f1,f2,f3,f4,f5,f6]);

$sl(3)$の場合は11本の斉次式が、そして$sl(4)$の場合は44本の斉次式が、被約グレブナー基底として出てくる。極小冪零軌道の定義方程式を計算しなくても、斉次式になってる時点で、冪零だって分かるし、後は上に書いたように次元勘定から、極小冪零軌道だと分かる


一応、念のため、$sl(3)$の場合、極小冪零軌道(の閉包)の定義方程式(の被約グレブナー基底)は、以下のRisa/Asirコードで得られる。これは、上の計算結果と一致する

/* defining equations for a minimal nilpotent orbit of sl(3) */
load("noro_pd.rr");

A = newmat(3,3,[[a11,a12,a13],[a21,a22,a23],[a31,a32,a33]]);   /* SL(3) */

M = (invmat(A)[0])*newmat(3,3,[[0,1,0],[0,0,0],[0,0,0]])*A;   /* orbit */

B = [e1 - 6*M[1][0] , e2 - 6*M[2][1] , e3-6*M[2][0] ,
     h1-12*M[0][0]-6*M[2][2] , h2+12*M[2][2]+6*M[0][0] , 
     f1-6*M[0][1] , f2-6*M[1][2] , f3-6*M[0][2] , det(A)-1];
V = [a11,a12,a13,a21,a22,a23,a31,a32,a33,e1,e2,e3,h1,h2,f1,f2,f3];
G = nd_gr_trace(B,V,1,1,[[0,9],[0,8]]);
D2 = noro_pd.elimination(G,[e1,e2,e3,h1,h2,f1,f2,f3]);

つまり、Josephの論文のminimal realizationの自由度というのは、Lie環の極小軌道の次元の半分となる。半分なのは、自由度=配位空間の次元で、極小軌道の次元=相空間の次元だから、ということで、極小軌道を量子化すれば、求めるminimalな"量子力学系"が得られる。よかった。


Hamiltonianが得られるわけではないので、"量子力学系"というのは、やや違和感があるけど、古典的には、相空間の構造が明らかになるという状況に対応する。冪零軌道はsymplectic等質多様体となっているという意味で、かなり特殊な相空間であるといえる(一般のsymplectic等質多様体と余随伴軌道の関係は、 Kirillov-Kostant-Souriau Theorem)


実際のKepler問題でも、群$SO(4,2)$はKepler多様体に推移的に作用して、相空間はsymplectic等質多様体となっている。spectrum generating algebraは相空間のsymplectic構造を保つ幾何学的対称性だというのが、一つの解釈になる。実際には、symplectic等質多様体でない相空間の方が普通なので、この解釈は狭すぎる。spectrum generating algebraは相空間上の関数環を量子化した代数(短くて的確な呼称が欲しい)で、相空間がsymplectic等質多様体の場合は、その代数は、Lie環の普遍展開環の商として得られると理解するほうが一般性があるのだと思う。Hamiltonianの与え方は色々な可能性があり可積分になるとは限らない


# 元々、素朴にはKepler問題は、3次元Euclid空間を配位空間としている(原点では、ポテンシャルが発散するので、原点を除きたければ除く)。相空間は、素朴には、これの余接空間。$SO(4,2)$をこの相空間に作用させた場合、行き先がwell-definedでないことは起こりうる(実際、どうなってるのか知らないけど)。似たような状況として、共形変換群のEuclid空間やMinkowski空間への作用で、特殊共形変換の行き先は"無限遠点"の場合がある。共形変換群は、数学的に厳密にはEuclid空間やMinkowski空間に推移的に作用しない。が、原点を動かさない部分群全体は分かるので、共形変換群を、そのような部分群で割ってやると、等質空間が得られる。こうした操作は、共形コンパクト化と呼ばれ、4次元Euclid空間の共形コンパクト化はS^4となる。Kepler問題の相空間の場合も、事情は同じようなことと思う。可積分系界隈では、物理とかで扱う素朴な相空間と、数学的に自然な相空間が、若干異なるということは、時々ある(けど、明示的に事情を説明してくれていることは少ない)


# ところで、水素原子の場合、包合的な第一積分は、丁度3つある($so(4)$のCartan部分環から2つとHamiltonian自身)。これらは独立なので、Kepler系はLiouvilleの意味で可積分。ところが、水素原子が3自由度という結論は古典的な場合からの推測であって、$so(4,2)$の極小表現とか一般の既約ユニタリー表現を、ぽっと与えても、その自由度は定義から明らかなものではない。妄想の一つとして、Gelfand-Kirillov次元で自由度を定義するということが考えられる。定義するのは自由で、Liouville-Arnoldの定理の量子力学版を示すというようなことをやらないと意味がないけど


などと考えると、以前書いた

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

『主量子数n+1の電子殻のヒルベルト空間は、so(4)の最高ウェイト(n,0)を持ついくつかの既約表現の直和で書けることが分かったけど、各既約表現が何個出てくるかということは、表現論的に知る方法はなくて、Schrodinger方程式を解く(あるいは、Pauliの時は、まだSchrodinger方程式がなかったのでHeisenberg方程式でやっていたようだけど)しかないように思う(状態のHilbert空間を、どう取るかは物理的な理由によって決まるべきもので、表現論の立場から決める原理はないから)』

については、水素原子の状態空間は、かなり表現論的な理由のみによって決められるのかもしれない。とはいえ、一方で、状態が束縛状態か散乱状態かは、少なくとも一見ではHamiltonianに依存する概念なので、そのへん、どう理解すればいいのか分からない



おまけ。Josephイデアルは、大雑把に言って、極小冪零軌道の定義方程式の非可換版。極小冪零軌道の定義方程式の計算は、グレブナー基底でできたので、Josephイデアルの生成元を決めるのは、非可換グレブナー基底を使えばできるだろうと考える。残念ながら、この計算は、$sl(3)$の場合ですら、3日くらい待っても終わらなかったので、断念した。

非可換グレブナー基底の実装は

[1] Bergman

http://servus.math.su.se/bergman/

[2] GAPのGBNP

http://www.gap-system.org/Packages/gbnp.html

[3] Singular

A.6.1 Left and two-sided Groebner bases

https://www.singular.uni-kl.de/Manual/4-0-3/sing_883.htm

などがある。Singularは定番。GAPは有限群メインと思ってたら、最近は半単純Lie代数計算用のパッケージとかあって、冪零軌道を扱うための諸々が入ってたりする。Bergmanは、未だに使い方が謎でよく分からないのだけど、以前どこかで見たベンチマークでは結構早かった気がする(他と違うアルゴリズムが使われていたはず)。Bergmanで計算する場合、以下のような感じ(終わらないけど)

[1]> (purelexify)
NIL
[2]> (noncommify)
NIL
[3]> (algforminput)
algebraic form input> vars q1,q2,p1,p2,h1,h2,e1,e2,e3,f1,f2,f3;
e1-q1,
e2-q2*p1,
e3-q2,
h1+2*q1*p1+q2*p2,
h2-q1*p1+q2*p2,
f1+q1*p1*p1+q2*p1*p2,
f2-q1*p2,
f3+q1*p1*p2+q2*p2*p2,
p1*q1-q1*p1-1,
p2*q2-q2*p2-1;
T
[4]> (homogeniseinput)
NIL
[5]> (groebnerinit)
NIL
[6]> (groebnerkernel)

余談と妄想。

$SU(2,2)$に付随する力学系

http://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/83876

では、4自由度の力学系の$U(1)$作用による、Hamiltonian reductionとして、Kepler系を得ている。簡約前の相空間は、"4次元Euclid空間から原点を除いた空間の余接空間で、従って、symplectic群Sp(8,R)が自然に作用している。一方、$su(2,2)$の極小表現は、Metaplectic群$Mp(8,R)$のSegal-Shale-Weil表現の部分空間上に実現されるのだった。この部分空間を定める拘束条件が、上記$U(1)$作用から定まるmoment mapが0となる条件の量子化と思える形になっている(Sp(8,R)の自然な作用に対応する運動量写像を正準量子化することによって、Segal-Shale-Weil表現は得られる)。きちんと確認したわけではないけれど、古典的には、極小冪零軌道と$U(1)$-Hamiltonian reductionの2通りの方法で、Kepler問題の相空間が得られ、表現論では、後者の構成が極小表現のladder representationによる実現と対応している形になると推測するのは自然


このHamiltonian reductionによる構成は、もうちょっと一般化できそうな気がする。多分、きちんと読んでないけど

The U(1)-Kepler Problems

https://arxiv.org/abs/0805.0833

The Sp(1)-Kepler Problems

https://arxiv.org/abs/0805.0840

あたりが、そういうことを考えているらしい。可積分系としては(初等的な表現論のみで解けてしまうので)それほど面白いわけでもないが。論文では相空間は明示的に出てこないが、(群$G$が空間$X$に自由に作用する場合)簡約相空間$T^{*}X // G$は$T^{*}(X/G)$に同型なので同じことをやっていると言える。もっと一般のdual pairに、このような構成がないかは分からない

HOL Lightの実装を理解したい(0)

理解したい。問題は、OCaml読めないことと、HOL Lightが、どういうものか微塵も分からんので、とりあえず、HOL Lightを動かしてみたメモ。ソースコード落とすと、中にTutorialがあるので適当に見れば良い気がする。


Coqと双璧をなすかどうかは知らないけど、HOL系のproof assistantは、よく見かける。Isabelle/HOLというのが一番(?)有名っぽい。こいつは

seL4という、マイクロカーネルタイプOSの形式検証

https://github.com/seL4/l4v

とかいう実績があるらしい。


HOL Lightは、

Kepler予想の証明を形式化したflyspeckプロジェクト

https://github.com/flyspeck/flyspeck

で使われたらしい。


HOL Light開発の経緯として

https://en.wikipedia.org/wiki/HOL_(proof_assistant)

には、'The second current implementation is HOL Light. This started as an experimental "minimalist" version of HOL. Although it has subsequently grown into another mainstream HOL variant, its logical foundations remain unusually simple.'と書いてある。Isabelleで書かれた証明とHOL Lightの証明は、全く似てなかったりするけど、できることは、基本的に同じと考えて、差し支えなさそう?


Towards self-verification of HOL Light

https://www.cl.cam.ac.uk/~jrh13/papers/holhol.html

の冒頭には"The HOL Light prover is based on a logical kernel consisting of about 400 lines of mostly functional OCaml, whose complete formal verification seems to be quite feasible."という一文がある。別にformal verificationしたいとは思わないのだけど、400行なら、なんか理解する気になると思った(あと、apt-get install hol-lightってしたらインストールされたので)


# HOL Lightのソースは

The HOL Light theorem prover (moved from Google code)

https://github.com/jrh13/hol-light

で公開されている。この中のfusion.mlが、『400行程度のlogical kernel』に相当する部分っぽい?コメント入れて670行くらい。冒頭で、lib.mlというのが、参照されているっぽいけど、curryだのuncurryだのの関数とか、リスト操作関連関数が定義されているだけ。Wikipediaの

HOL Light#Logical Foundations

https://en.wikipedia.org/wiki/HOL_Light#Logical_foundations

にあるのと、同じ名前の関数が定義されている(REFLとかMK_COMBとかetc.)。大体わかった(適当


HOL LightとCoqで、どういう差があるのか、よく分かってないけど、HOL Lightでは

・命題はただのブール値であり排中律が成立する

・functional exntensionalityが成立する

・選択公理も成立する(?)

・依存型はない

という感じであるらしい。



HOL Lightの例1。

# let thp= ASSUME `p:bool`;;
val thp : thm = p |- p
# let thq = ASSUME `q:bool`;;
val thq : thm = q |- q
# let thpq = CONJ thp thq;;     
val thpq : thm = p, q |- p /\ q

何をしてるかは見ればわかる。`...`は、何かtermというものになるっぽい

# let t0 = `x+1`;;
val t0 : term = `x + 1`
# `1+2`;;
val it : term = `1 + 2`

実数(というか浮動小数点数?)はダメらしい

# `3.4`;;
Exception: Failure "Unparsed input following term".

何もつけないと

# 3;;
val it : int = 3

となる。簡単な計算もできる

# (NUM_REDUCE_CONV `2 + 2`);;
val it : thm = |- 2 + 2 = 4

fusion.ml見ても、ASSUMEはあるけど、CONJはないし、NUM_REDUCE_CONVもない。CONJはbool.mlに定義されている

let CONJ =
  let f = `f:bool->bool->bool`
  and p = `p:bool`
  and q = `q:bool` in
  let pth =
    let pth = ASSUME p
    and qth = ASSUME q in
    let th1 = MK_COMB(AP_TERM f (EQT_INTRO pth),EQT_INTRO qth) in
    let th2 = ABS f th1 in
    let th3 = BETA_RULE (AP_THM (AP_THM AND_DEF p) q) in
    EQ_MP (SYM th3) th2 in
  fun th1 th2 -> substitute_proof (
    let th = INST [concl th1,p; concl th2,q] pth in
    PROVE_HYP th2 (PROVE_HYP th1 th))
  (proof_CONJ (proof_of th1) (proof_of th2));;

なるほど、よく分からん。これの少し上の方に

let AND_DEF = new_basic_definition
 `(/\) = \p q. (\f:bool->bool->bool. f p q) = (\f. f T T)`;;

というのがある。やりたいことはわかる


termの型はtype_ofで取れる

# type_of `1`;;
val it : hol_type = `:num`
# type_of `1=2`;;
val it : hol_type = `:bool`
# type_of `[1,2,3]`;;
val it : hol_type = `:(num#num#num)list`

hol_typeはfusion.mlで定義されている

  type hol_type = Tyvar of string
                | Tyapp of string * hol_type list

だと思う


自然数は最初から組み込まれているが、自分で定義する場合は、以下のようにするとできる。

# let nat_ind , nat_rec = define_type "Nat = Zero | Suc Nat";;
val nat_ind : thm = |- !P. P Zero /\ (!a. P a ==> P (Suc a)) ==> (!x. P x)
val nat_rec : thm =
  |- !f0 f1. ?fn. fn Zero = f0 /\ (!a. fn (Suc a) = f1 a (fn a))

これで、

# type_of `Zero`;;
val it : hol_type = `:Nat`
# type_of `Suc Zero`;;
val it : hol_type = `:Nat`

とか、使えるようになる。


一方、組みこまれてる奴は

# type_of `SUC`;;
val it : hol_type = `:num->num`
# mk_comb(`SUC`,`0`);;
val it : term = `SUC 0`

という感じ。


組み込みのnumの方は

# NUM_REDUCE_CONV `(SUC 0)`;;  
val it : thm = |- SUC 0 = 1

勝手に、human-friendlyな方に正規化してくれる。


帰納法。HOL LightにもCoq同様、tacticがある。forall(x:num),x+0=xは

# let thm0 = prove(`!x.x+0=x`,INDUCT_TAC THEN ASM_REWRITE_TAC[ADD_CLAUSES]);;
val thm0 : thm = |- !x. x + 0 = x

のように書ける(!はforallを表す)。proveは以下の型を持つ関数

# prove;;
val it : term * tactic -> thm = <fun>

ほとんどの(多分全ての)人類は極めてIQが低く愚かなので、複雑なケースではtacticを一発で書くのは難しい。Coqのように、subgoalを確認しながら進める場合

# g `!x.x+0=x`;;
val it : goalstack = 1 subgoal (1 total)

`!x. x + 0 = x`

# e INDUCT_TAC;;
val it : goalstack = 2 subgoals (2 total)

  0 [`x + 0 = x`]

`SUC x + 0 = SUC x`

`0 + 0 = 0`

# e (ASM_REWRITE_TAC[ADD_CLAUSES]);;
val it : goalstack = 1 subgoal (1 total)

  0 [`x + 0 = x`]

`SUC x + 0 = SUC x`

# e (ASM_REWRITE_TAC[ADD_CLAUSES]);;
val it : goalstack = No subgoals
# let thm1 = top_thm();;
val thm1 : thm = |- !x. x + 0 = x



numはnums.mlで定義されている

let NUM_REP_RULES,NUM_REP_INDUCT,NUM_REP_CASES =
  new_inductive_definition
   `NUM_REP IND_0 /\
    (!i. NUM_REP i ==> NUM_REP (IND_SUC i))`;;

let num_tydef = new_basic_type_definition
  "num" ("mk_num","dest_num")
    (CONJUNCT1 NUM_REP_RULES);;

let ZERO_DEF = new_definition
 `_0 = mk_num IND_0`;;

let SUC_DEF = new_definition
`SUC n = mk_num(IND_SUC(dest_num n))`;;

だと思う。

# num_tydef;;
val it : thm * thm =
  (|- mk_num (dest_num a) = a, |- NUM_REP r <=> dest_num (mk_num r) = r)

nums.mlのもうちょっと下の方を見ると

# prove(`_0 = 0`,REWRITE_TAC[NUMERAL]);;
val it : thm = |- _0 = 0

みたくできるのがわかる。NUMERALの定義は

# NUMERAL;;
val it : thm = |- !n. NUMERAL n = n
# type_of `NUMERAL`;;
val it : hol_type = `:num->num`

となっている。紛らわしいけど、項としてのNUMERALは、numからnumへの関数で、定理NUMERALは、単にそれが恒等関数であることを主張している。

let NUMERAL =
  let num_ty = type_of(lhand(concl ZERO_DEF)) in
  let funn_ty = mk_fun_ty num_ty num_ty in
  let numeral_tm = mk_var("NUMERAL",funn_ty) in
  let n_tm = mk_var("n",num_ty) in
new_definition(mk_eq(mk_comb(numeral_tm,n_tm),n_tm));;

なので、 `!n. NUMERAL n = n`は実際にはNUMERALの定義なのだと思う。で、`0`は内部的には`NUMERAL _0`として扱われているっぽい


tactic嫌いな人間は

# INST [`_0` ,`n:num`] (List.hd(mk_rewrites false NUMERAL []));;
val it : thm = |- 0 = _0

のようにやる。



線形代数をやることを考える。Multivariate/vector.mlに色々定義がある。realという型は、起動時点で作られている。

# needs "Multivariate/vectors.ml";;

とか打つと、何か一杯出てきて使えるようになる。で

# `!x:real^N. x dot x >= real(0)`;;
val it : term = `!x. x dot x >= real 0`

と書けるようになる。しかしNってなんだ。どっから来たんだって感じである

# types();;
val it : (string * int) list =
  [("4", 0); ("3", 0); ("2", 0); ("finite_sum", 2); ("cart", 2);
   ("finite_image", 1); ("int", 0); ("real", 0); ("hreal", 0); ("nadd", 0);
   ("char", 0); ("list", 1); ("option", 1); ("sum", 2); ("recspace", 1);
   ("num", 0); ("ind", 0); ("prod", 2); ("1", 0); ("bool", 0); ("fun", 2)]

を見てもいない。


別に、Nじゃなくてもいいらしい

# `!x:real^M. x dot x >= real(0)`;; 
val it : term = `!x. x dot x >= real 0`
# `!x:real^hoge. x dot x >= real(0)`;;
val it : term = `!x. x dot x >= real 0`

以下は、正しくない命題だが

# `!x:real^N.x dot x < real (dimindex(:N)) `;;
val it : term = `!x. x dot x < real (dimindex (:N))`

とか書ける。HOL Lightは依存型とかがないので、"forall (n:num) (u v:real^n),...."みたいなことが書けない。代わりに、(実)n次元ベクトル空間を{1,2,...,n} -> realみたいな関数と同一視する。"forall (X:Type) (u v:X->real),..."みたいには書けないけど、HOL Lightでは、forall (u v :X->real),...みたいに書いてあると、Xは型なんだなと推論してくれる仕組みがあるっぽい?(嘘かも)。考えてみると、リストを操作する関数とかで、こういう多相関数を許す仕組みがないと、辛いことになる


fusion.mlにあるmk_vartypeというのを使うと

# mk_vartype "X";;
val it : hol_type = `:X`

というのが作れる。内部的には、(Tyvar "X")という風になっているはず。"HOL Light: an overview"という論文に"HOL Light's logic is simple type theory with polymorphic type variables"と書いてあって、polymorphic type variablesというのは、このTyvarを指しているっぽい。


dimindexはcart.mlで定義されているけど

let dimindex = new_definition
`dimindex(s:A->bool) = if FINITE(:A) then CARD(:A) else 1`;;

みたいに書いてある。dimindexとかFINITEとか、直接、型を引数に取ることはできないのだけど、sets.mlで

# UNIV;;
val it : thm = |- (:A) = (\x. T)

みたいのが定義されている。で、printer.ml見ると

(* Print the universal set UNIV:A->bool as "(:A)".                           *)

とか書いてあったり。なんか、そういうtrickが使われている


そんなわけで

# `!f:A->B x:A.f(x)=f(x)`;;
val it : term = `!f x. f x = f x`

もOK.

自明だけど証明してみる

# g`!f:A->B x:A.f(x)=f(x)`;;
val it : goalstack = 1 subgoal (1 total)

`!f x. f x = f x`

# e STRIP_TAC;;
val it : goalstack = 1 subgoal (1 total)

`!x. f x = f x`

# e STRIP_TAC;;
val it : goalstack = 1 subgoal (1 total)

`f x = f x`

# e REFL_TAC;;
val it : goalstack = No subgoals

まとめ。Coq使ったことある人は、あまり戸惑うこともなさそうな感じ。実装も、かなり単純で、読めば大体分かる気がする(読んだことあるproof assistantの実装が、HOL Lightとcubical interpreterしかないけど、cubical interpreterは、それほど自明じゃない感じがある)。proof assistant使ったことない人には、あまり親切とは言えないと思う(そもそも、情報が少ない)。flyspeckの人は、なぜHOL Lightを選んだのか謎である。

2016-06-14

水素原子の表現論(1.5)spectrum generating algebra

水素原子の束縛状態全体は、Lie代数so(4,2)の極小表現になっていることが知られている。このso(4,2)は、dynamical symmetryと呼ばれたり、spectrum generating algebraと呼ばれていたりする。個人的には、ハミルトニアンと可換でないものを"symmetry"と呼ぶのは紛らわしいと思うので、spectrum generating algebraという名前を使用する。spectrum generatingというのも、いまいちな名前とは思う。そして、伝説によれば、散乱状態の全体も、同様にso(4,2)の既約ユニタリー表現になっているらしいのだけど、散乱状態の方は、出典不明で、本当かどうかも、よく分からない。


#表現論的に理解したいこととして、一つには、散乱状態の全体を、何かよい基底を見つけて、so(4,2)の作用を書きたいというのがある。もう一つは、極小表現がso(4)の既約ユニタリ表現の直和に分解できるのと同様、散乱状態の方も、エネルギー一定の状態はso(3,1)の既約ユニタリ表現をなすので、散乱状態全体は、これらの表現の直積分で書けるというようなことが期待される。知る限り、そのような説明を見たことがない


束縛状態に関しては、spectrum generating algebraの具体的な実現は分かっているし、一方で、全ての束縛状態も分かっていて、極小表現は最低ウェイト表現であるから、基底状態(1s軌道解)に対して、最低ウェイト条件と特異ベクトルが消えるという条件の成立を確認すれば、水素原子の束縛状態がso(4,2)の極小表現であるということは証明できる。so(4,2)の極小表現を、so(4)に制限した時の既約分解の仕方などは、表現論の知識のみで分かる


ところで以前、

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

で『主量子数n+1の電子殻のヒルベルト空間は、so(4)の最高ウェイト(n,0)を持ついくつかの既約表現の直和で書けることが分かったけど、各既約表現が何個出てくるかということは、表現論的に知る方法はなくて、Schrodinger方程式を解く(あるいは、Pauliの時は、まだSchrodinger方程式がなかったのでHeisenberg方程式でやっていたようだけど)しかないように思う(状態のHilbert空間を、どう取るかは物理的な理由によって決まるべきもので、表現論の立場から決める原理はないから)』と書いた


けれどもし、

(1)spectrum generating algebraの生成元が分かって

(2)束縛状態はspectrum generating algebraの何らかの既約ユニタリー表現であるということは分かっているが

(3)にも関わらず、シュレディンガー方程式の解き方は分からない

というような人がいたとしたら、束縛状態の全体について、どのくらいのことが分かるか、考えたい。現在のところ、条件(2)とか、束縛状態の全体が何故so(4,2)の極小表現となる理由を、シュレディンガー方程式を解かずに理解する方法はないと思うけど


とりあえず考えることとして、spectrum generating algebraの生成元の間には、いくつかの非自明な代数的関係式(角運動量とRunge-Lenzベクトルが直交するというような)が成立し、それらの関係式によって、U(so(4,2))の両側イデアル$J$が決定される。この両側イデアルは、極小表現上では、当然0になるはずである。$J$が消えるような、既約ユニタリー表現は、比較的限定されるだろうと期待したい。


一方、極小表現上で消えるU(so(4,2))の元全体を考えることができ、やはり両側イデアルをなす(これにはJosephイデアルという名前が付いている)。まずは、この2つのイデアルに、どれくらいの差があるか、ということを調べたい



Kepler問題に於けるspectrum generating algebraの生成元の具体的な式は、例えば、以下の論文に書いてある

SO(4,2)‐Formulation of the Symmetry Breaking in Relativistic Kepler Problems with or without Magnetic Charges

http://scitation.aip.org/content/aip/journal/jmp/12/5/10.1063/1.1665653


簡単のため、長さをBohr半径で規格化して、プランク定数も1とすると、生成元は、以下の15個

角運動量

L_1 = -y ¥frac{¥partial}{¥partial z} + z ¥frac{¥partial}{¥partial y}

L_2 = -z ¥frac{¥partial}{¥partial x} + x ¥frac{¥partial}{¥partial z}

L_3 = -x ¥frac{¥partial}{¥partial y} + y ¥frac{¥partial}{¥partial x}

scaled Runge-Lenz vector

A_1 = -¥frac{1}{2} x (¥frac{¥partial^2}{¥partial x^2} + ¥frac{¥partial^2}{¥partial y^2} + ¥frac{¥partial^2}{¥partial z^2}) + ¥frac{¥partial}{¥partial x} ¥cdot ( x¥frac{¥partial}{¥partial x} + y ¥frac{¥partial}{¥partial y} + ¥frac{¥partial}{¥partial z}) - ¥frac{1}{2} x

A_2 = -¥frac{1}{2} y (¥frac{¥partial^2}{¥partial x^2} + ¥frac{¥partial^2}{¥partial y^2} + ¥frac{¥partial^2}{¥partial z^2}) + ¥frac{¥partial}{¥partial y} ¥cdot ( x¥frac{¥partial}{¥partial x} + y ¥frac{¥partial}{¥partial y} + ¥frac{¥partial}{¥partial z}) - ¥frac{1}{2} y

A_3 = -¥frac{1}{2} z (¥frac{¥partial^2}{¥partial x^2} + ¥frac{¥partial^2}{¥partial y^2} + ¥frac{¥partial^2}{¥partial z^2}) + ¥frac{¥partial}{¥partial z} ¥cdot ( x¥frac{¥partial}{¥partial x} + y ¥frac{¥partial}{¥partial y} + ¥frac{¥partial}{¥partial z}) - ¥frac{1}{2} z


dual Runge-Lenz vector

B_1 = A_1 + x

B_2 = A_2 + y

B_3 = A_3 + z


other generators

¥Gamma_1 = r ¥frac{¥partial}{¥partial x}

¥Gamma_2 = r ¥frac{¥partial}{¥partial y}

¥Gamma_3 = r ¥frac{¥partial}{¥partial z}

T_1 = -¥frac{r}{2}(¥frac{¥partial^2}{¥partial x^2} + ¥frac{¥partial^2}{¥partial y^2} + ¥frac{¥partial^2}{¥partial z^2}) - ¥frac{r}{2}

T_2 = x ¥frac{¥partial}{¥partial x} + y ¥frac{¥partial}{¥partial y} + z ¥frac{¥partial}{¥partial z} + 1

T_3 = -¥frac{r}{2}(¥frac{¥partial^2}{¥partial x^2} + ¥frac{¥partial^2}{¥partial y^2} + ¥frac{¥partial^2}{¥partial z^2}) + ¥frac{r}{2}


そして、Chevalley生成元を以下のように決める。

e_1 = (A_1 + L_2 + i(L_1 - A_2))/2

f_1 = (A_1 - L_2 + i(L_1 + A_2))/2

h_1 = -(A_3 + i L_3)

e_2 = i(T_1+T_2+B_3-¥Gamma_3)/2

f_2 = i(T_1-T_2+B_3+¥Gamma_3)/2

h_2 = A_3+T_3

e_3 = (A_1 + L_2 -i(L_1 -A_2))/2

f_3 = (A_1 - L_2 -i(L_1+A_2))/2

h_3 = -A_3 + i L_3

末尾に、Risa/Asirでの確認用コードを載せてある。ただし、多項式係数でない微分作用素は扱いづらいので、付加的な変数$r,r'$を使って

p_x = ¥frac{¥partial}{¥partial x} + r’ x ¥frac{¥partial}{¥partial r} - r’^3 x ¥frac{¥partial}{¥partial r’}

etc.

などとして、関係式r’ = r^{-1} , r^2 = x^2+y^2+z^2によるmoduloを取ることで、計算している


長さをBohr半径で規格化しているので、1s軌道解は、規格化定数を無視すれば、単にexp(-r)である。最低ウェイト条件と特異ベクトルが消えることを確認するのは容易い。確認すべき条件は

Maxwell方程式の表現論

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

で書いた(極小表現は、helicity0のmassless既約表現と同じ)。

e_1 ¥phi = e_3 ¥phi = 0

f_1 ¥phi = f_2 ¥phi = f_3 ¥phi = 0

h_1 ¥phi = (h_2-1)¥phi = h_3 ¥phi = 0


ここまでが準備。それで両側イデアル$J$を計算したいし、多分、機械的に計算していくことは可能だけど、直接計算は面倒そうなので、Josephイデアルの生成元を調べて、それらは勿論極小表現の上では消える(はず)だけど、微分作用素として恒等的に消えるかどうかチェックするという手順を取ろうと思う(続く)


Josephイデアルの構造については、Garfinkleによる結果が知られている

A new construction of the Joseph ideal

https://dspace.mit.edu/handle/1721.1/15620#files-area


その他、本題とは関係ない補足

(1)Josephイデアルのassociated varietyは、極小冪零軌道の閉包になっていて、極小冪零軌道は、Kepler問題の古典的相空間とみなすことができる(散乱状態は無視している)。これには、Kepler多様体という名前がついているらしい。また、表現空間の中で、generalized coherent states全体を考えることができ、Kepler多様体と同型になると思う(要確認)

(2)Coulomb problemは、(一種の拘束条件付きの)4次元調和振動子へ変換できることが知られている(この座標変換にはKustaanheimo-Stiefel変換という名前がついているらしい。この変換は水素原子の経路積分でも使われたらしい(要確認))。数学的にみると、極小表現の4変数多項式環での実現(Segal-Shale-Weilによる?)と同じもの。この実現は

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

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

にも書いてある。次数$n$の多項式が(n+1)次元で、量子数n+1の電子軌道と対応する。発見の順番とか、歴史的経緯とかは知らない(この実現自体は、sp(4,R)のmetaplectic表現を、su(2,2)に制限して得られ、metaplectic表現は、物理にinspireされて生まれたと見たことがある気がする)。Kustaanheimo-Stiefel変換はある種のmagicのように私には見えるけど、表現論的には、同じ表現を二通りの方法で実現して、その間のintertwinnerを作っているのだ、と理解できる。そのような理解を得るためには、so(4)対称性だけ見ていたのでは不十分である


(3)普遍展開環$U(so(4,2))$をJosephイデアルで割った代数$hs(4,2)$は、higher spin algebraという名前で、最近(というのは主に21世紀になってから)、物理の人が調べたりしているらしい


/**** common stuffs ****/
load("noro_pd.rr");
load("noro_matrix.rr");

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


def appdiff_rat(P , F , V){
   N = length(V)/2;
   if(type(P)==9){  /* 分散多項式表現の時 */
      DP = dp_dtop(P, V);
   }else if(type(P)<=2){
      DP = P;
   }else{
      return;
   }
   Terms = p_terms(DP , V , 2);
   for(Ret = 0 , Terms = Terms; Terms!=[] ; Terms=cdr(Terms)){
      T = car(Terms);
      if(type(F)==9){
           Tmp = dp_dtop(F , V);
      }else{
           Tmp = F;
      }
      if(T==1){
         Ret += p_nf(DP , V , V ,0)*Tmp;
      }else{
         C = p_nf(sdiv(DP-p_nf(DP,[T],V,0) , T) , V , V , 0);  /* 単項式Tの係数 */
         for(I = 0 ; I < N ; I++){
            Deg = deg(T , V[N+I]);
            for(J = 0 ; J < Deg ; J++){
               Tmp = red(diff(Tmp , V[I]));
            }
         }
         for(I = 0 ; I < N ; I++){
            Deg = deg(T , V[I]);
            Tmp *= V[I]^Deg;
         }
         Ret += C*Tmp;
      }
   }
   return red(Ret);
}


def appdiff_complex(P , F , V){
    if(type(P)==9){
        DP = dp_dtop(P,V);
    }else{
        DP = P;
    }
    RP = real(DP);
    IP = imag(DP);
    if(type(F)<=2){
        RF = real(F);
        IF = imag(F);
    }else if(type(F)==3){ /*有理式*/
        DN_conj = conj(dn(F));
        if(DN_conj != dn(F)){
           TF = red(DN_conj*nm(F)/(DN_conj*dn(F)));
        }else{
           TF = F;
        }
        RF = red(real(nm(TF))/dn(TF));
        IF = red(imag(nm(TF))/dn(TF));
    }
    Ret = red( appdiff_rat(RP , RF , V) - appdiff_rat(IP,IF,V) );
    Ret += @i*red( appdiff_rat(IP,RF,V) + appdiff_rat(RP,IF,V) );
    return Ret;
}

/* 
微分作用Pを式Fに作用させる 
Vは変数のリスト(dp_dtopなどで使われるのと同様の形式)
*/
def appdiff(P,F,V){
    if(type(P)!=6){
        if(type(F)<=3){
           return appdiff_complex(P,F,V);
        }else if(type(F)==5){   /* vector */
           Ret = newvect(length(F));
           for(I=0;I<length(F);I++){
              Ret[I] = appdiff_complex(P, F[I] , V);
           }
           return Ret;
        }
    }else{ /*行列*/
        Sz = size(P);
        N = Sz[0];
        M = Sz[1];
        if(type(F)==5){
           Ret = newvect(N);
           for (I = 0 ; I < N ; I++){
              for(J = 0 ; J < M ; J++){
                Ret[I] += appdiff_complex(P[I][J] , F[J] , V);
              }
              Ret[I] = red(real(nm(Ret[I]))/dn(Ret[I])) + @i*red(imag(nm(Ret[I]))/dn(Ret[I]));
           }
           return Ret;
        }else if(type(F)==6){
            Sz2 = size(F);
            N2 = Sz2[0];
            M2 = Sz2[1];
            Ret = newmat(N,M2);
            for (I = 0 ; I < N ; I++){
                for(J = 0 ; J < M2 ; J++){
                    for(K = 0 ; K < M ; K++){
                       Ret[I][J] += appdiff_complex(P[I][K] , F[K][J] , V);
                    }
                }
            }
            return Ret;
        }
    }
}


def dp_weyl_matmul(A,B){
    Sz1 = size(A);
    Sz2 = size(B);
    L = Sz1[0];
    M = Sz1[1];
    N = Sz2[0];
    Ret = newmat(L,N);
    for(I = 0; I < L ; I++){
       for(J = 0 ; J < N ; J++){
          for(K = 0 ; K < M ; K++){
               Ret[I][J] += dp_weyl_mul(A[I][K] , B[K][J]);
          }
       }
    }
    return Ret;
}

/*===== end of common stuffs =====*/


/* stuffs for 3D hydrogen atom */
V=[x1,x2,x3,r,ri,d1,d2,d3,dr,dri];
X1=dp_ptod(x1,V);
X2=dp_ptod(x2,V);
X3=dp_ptod(x3,V);
R=dp_ptod(r,V);
RI=dp_ptod(ri,V);
D1=dp_ptod(d1,V);
D2=dp_ptod(d2,V);
D3=dp_ptod(d3,V);
DR = dp_ptod(dr,V);
DRI = dp_ptod(dri,V);
P1 = D1 + dp_weyl_mul(RI*X1 , DR) - dp_weyl_mul(RI*RI*RI*X1 , DRI);
P2 = D2 + dp_weyl_mul(RI*X2 , DR) - dp_weyl_mul(RI*RI*RI*X2 , DRI);
P3 = D3 + dp_weyl_mul(RI*X3 , DR) - dp_weyl_mul(RI*RI*RI*X3 , DRI);
I=gr([r*ri-1,r*r-x1*x1-x2*x2-x3*x3],V,0);
Delta = dp_weyl_mul(P1,P1)+dp_weyl_mul(P2,P2)+dp_weyl_mul(P3,P3); /* Laplacian */


/* angular momentum */
L1 = -dp_weyl_mul(X2,P3) + dp_weyl_mul(X3,P2);
L2 = -dp_weyl_mul(X3,P1) + dp_weyl_mul(X1,P3);
L3 = -dp_weyl_mul(X1,P2) + dp_weyl_mul(X2,P1);

/* check [L1,L2]=L3 , [L2,L3] = L1 , [L3,L1] = L2 */
assert(dp_weyl_mul(L1,L2)-dp_weyl_mul(L2,L1)==L3);
assert(dp_weyl_mul(L2,L3)-dp_weyl_mul(L3,L2)==L1);
assert(dp_weyl_mul(L3,L1)-dp_weyl_mul(L1,L3)==L2);

T1 = (1/2)*(-R*Delta-R);
T2 = X1*P1 + X2*P2 + X3*P3 + 1;
T3 = (1/2)*(-R*Delta+R);
/* check [T1,T2]=T3 , [T2,T3]=-T1 , [T3,T1]=T2 */
/* T1,T2,T3 generates so(2,1) */
assert(p_nf(dp_dtop(dp_weyl_mul(T1,T2)-dp_weyl_mul(T2,T1)-T3,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T2,T3)-dp_weyl_mul(T3,T2)+T1,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,T1)-dp_weyl_mul(T1,T3)-T2,V),I,V,0)==0);



/* scaled Runge-Lenz vector */
A1 = -(1/2)*X1*Delta + dp_weyl_mul(P1,X1*P1+X2*P2+X3*P3) - X1/2;
A2 = -(1/2)*X2*Delta + dp_weyl_mul(P2,X1*P1+X2*P2+X3*P3) - X2/2;
A3 = -(1/2)*X3*Delta + dp_weyl_mul(P3,X1*P1+X2*P2+X3*P3) - X3/2;

/* dual Runge-Lenz vector */
B1 = -(1/2)*X1*Delta + dp_weyl_mul(P1,X1*P1+X2*P2+X3*P3) + X1/2;
B2 = -(1/2)*X2*Delta + dp_weyl_mul(P2,X1*P1+X2*P2+X3*P3) + X2/2;
B3 = -(1/2)*X3*Delta + dp_weyl_mul(P3,X1*P1+X2*P2+X3*P3) + X3/2;
   

assert(dp_weyl_mul(A1,A2) - dp_weyl_mul(A2,A1) == -L3);
assert(dp_weyl_mul(A2,A3) - dp_weyl_mul(A3,A2) == -L1);
assert(dp_weyl_mul(A3,A1) - dp_weyl_mul(A1,A3) == -L2);
assert(dp_weyl_mul(A1,L2) - dp_weyl_mul(L2,A1)==A3);
assert(dp_weyl_mul(@i*A1+L1,@i*A2+L2)-dp_weyl_mul(@i*A2+L2,@i*A1+L1)==2*(@i*A3+L3));


assert(dp_weyl_mul(B1,B2)-dp_weyl_mul(B2,B1)==L3);
assert(dp_weyl_mul(B2,B3)-dp_weyl_mul(B3,B2)==L1);
assert(dp_weyl_mul(B3,B1)-dp_weyl_mul(B1,B3)==L2);
assert(dp_weyl_mul(B1,L2) - dp_weyl_mul(L2,B1)==B3);
assert(dp_weyl_mul(B2,L3) - dp_weyl_mul(L3,B2)==B1);
assert(dp_weyl_mul(B3,L1) - dp_weyl_mul(L1,B3)==B2);
assert(dp_weyl_mul(B1+L1,B2+L2)-dp_weyl_mul(B2+L2,B1+L1)==2*(B3+L3));


assert(dp_weyl_mul(A1,B1)-dp_weyl_mul(B1,A1)==T2);
assert(dp_weyl_mul(A2,B2)-dp_weyl_mul(B2,A2)==T2);
assert(dp_weyl_mul(A3,B3)-dp_weyl_mul(B3,A3)==T2);
assert(dp_weyl_mul(T2,A1)-dp_weyl_mul(A1,T2)+B1==0);
assert(dp_weyl_mul(T2,A2)-dp_weyl_mul(A2,T2)+B2==0);
assert(dp_weyl_mul(T2,A3)-dp_weyl_mul(A3,T2)+B3==0);
assert(dp_weyl_mul(T2,B3)-dp_weyl_mul(B3,T2)+A3==0);


G1 = R*P1;
G2 = R*P2;
G3 = R*P3;
assert(p_nf(dp_dtop(dp_weyl_mul(G1,T1)-dp_weyl_mul(T1,G1)-A1,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(G2,T1)-dp_weyl_mul(T1,G2)-A2,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(G3,T1)-dp_weyl_mul(T1,G3)-A3,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(G1,T3)-dp_weyl_mul(T3,G1)-B1,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T1,A2)-dp_weyl_mul(A2,T1)-G2,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T1,A3)-dp_weyl_mul(A3,T1)-G3,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T1,B1)-dp_weyl_mul(B1,T1),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T1,B2)-dp_weyl_mul(B2,T1),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T1,B3)-dp_weyl_mul(B3,T1),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,A1)-dp_weyl_mul(A1,T3),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,A2)-dp_weyl_mul(A2,T3),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,B1)-dp_weyl_mul(B1,T3)+G1,V),I,V,0)==0)
assert(p_nf(dp_dtop(dp_weyl_mul(T3,B2)-dp_weyl_mul(B2,T3)+G2,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,B3)-dp_weyl_mul(B3,T3)+G3,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,G3)-dp_weyl_mul(G3,T3)+B3,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(A3,G3)-dp_weyl_mul(G3,A3)+T1,V),I,V,0)==0);


assert(p_nf(dp_dtop(dp_weyl_mul(A3,T3)-dp_weyl_mul(T3,A3),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(L3,T3)-dp_weyl_mul(T3,L3),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(L3,A3)-dp_weyl_mul(A3,L3),V),I,V,0)==0);



/*
$T_1 (e^{-r}) = (1-r)e^{-r}$
$T_2 (e^{-r}) = (1-r)e^{-r}$
$T_3 (e^{-r}) = e^{-r}$

$A_1 (e^{-r}) = 0$
$A_2 (e^{-r}) = 0$
$A_3 (e^{-r}) = 0$

$L_1 (e^{-r}) = 0$
$L_2 (e^{-r}) = 0$
$L_3 (e^{-r}) = 0$

$\phi$: lowest weight vector for minrep of so(4,2)

$f_1(\phi) = f_2(\phi) = f_3(\phi) = 0$
$h_1 \phi = (h_2-1)\phi = h_3 \phi = 0$
$e_1 \phi = e_3 \phi = 0$

$h_4 = h_1+h_2 , h_5 = h_2+h_3 , h_6 = h_1+h_2+h_3$
*/

/* Chevalley generators */
E1 = (A1+L2+@i*L1-@i*A2)/2;
F1 = (A1-L2 +@i*L1+@i*A2)/2;
H1 = -(A3+@i*L3);

E2 = @i*(T1+T2+B3-G3)/2;
F2 = @i*(T1-T2+B3+G3)/2;
H2 = A3+T3;

E3 = -(@i*(@i*A1+L1)-(@i*A2+L2))/2;
F3 = (A1-L2-@i*L1-@i*A2)/2;
H3 = -A3+@i*L3;

assert(dp_weyl_mul(E1,F1)-dp_weyl_mul(F1,E1)==H1);
assert(dp_weyl_mul(H1,E1)-dp_weyl_mul(E1,H1)==2*E1);
assert(dp_weyl_mul(H1,F1)-dp_weyl_mul(F1,H1)==-2*F1);

assert(p_nf(@i*dp_dtop(dp_weyl_mul(H2,E2)-dp_weyl_mul(E2,H2)-2*E2,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(E2,F2)-dp_weyl_mul(F2,E2)-H2,V),I,V,0)==0);
assert(p_nf(@i*dp_dtop(dp_weyl_mul(H2,F2)-dp_weyl_mul(F2,H2)+2*F2,V),I,V,0)==0);

assert(dp_weyl_mul(E3,F3)-dp_weyl_mul(F3,E3)==H3);
assert(dp_weyl_mul(H3,E3)-dp_weyl_mul(E3,H3)==2*E3);
assert(dp_weyl_mul(H3,F3)-dp_weyl_mul(F3,H3)==-2*F3);

assert(dp_weyl_mul(H1,E3)-dp_weyl_mul(E3,H1)==0);
assert(dp_weyl_mul(H3,E1)-dp_weyl_mul(E1,H3)==0);
assert(dp_weyl_mul(E1,F3)-dp_weyl_mul(F3,E1)==0);
assert(dp_weyl_mul(E3,F1)-dp_weyl_mul(F1,E3)==0);

/* [H2 , E1] = -E1 */
assert(p_nf(dp_dtop(dp_weyl_mul(A3+T3,A1+L2)-dp_weyl_mul(A1+L2,A3+T3)+(A1+L2),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(A3+T3,L1-A2)-dp_weyl_mul(L1-A2,A3+T3)+(L1-A2),V),I,V,0)==0);

/* [H2 , F1] = F1 */
assert(p_nf(dp_dtop(dp_weyl_mul(A3+T3,L1+A2)-dp_weyl_mul(L1+A2,A3+T3)-(L1+A2),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(A3+T3,A1-L2)-dp_weyl_mul(A1-L2,A3+T3)-(A1-L2),V),I,V,0)==0);

/* [H1,E2] = -E2 */
assert(p_nf(dp_dtop(@i*(dp_weyl_mul(H1,E2)-dp_weyl_mul(E2,H1)+E2),V),I,V,0)==0);

/* [H1,F2] = F2 */
assert(p_nf(dp_dtop(@i*(dp_weyl_mul(H1,F2)-dp_weyl_mul(F2,H1)-F2),V),I,V,0)==0);



/* some vanishing conditions */ 
assert(dp_weyl_mul(L1,A1)+dp_weyl_mul(L2,A2)+dp_weyl_mul(L3,A3)==0);

/* A1*A1+A2*A2+A3*A3-L1*L1-L2*L2-L3*L3==T3*T3-1 */
assert(p_nf(dp_dtop(dp_weyl_mul(A1,A1)+dp_weyl_mul(A2,A2)+dp_weyl_mul(A3,A3)-dp_weyl_mul(L1,L1)-dp_weyl_mul(L2,L2)-dp_weyl_mul(L3,L3)-dp_weyl_mul(T3,T3)+1,V),I,V,0)==0);

/* B1*B1+B2*B2+B3*B3+L1*L1+L2*L2+L3*L3==T1*T1+1 */
assert(p_nf(dp_dtop(dp_weyl_mul(B1,B1)+dp_weyl_mul(B2,B2)+dp_weyl_mul(B3,B3)+dp_weyl_mul(L1,L1)+dp_weyl_mul(L2,L2)+dp_weyl_mul(L3,L3)-dp_weyl_mul(T1,T1)-1,V),I,V,0)==0);


三角形の合同類のなす空間と有限変形理論

三角形の合同条件と不変式論

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

みたいなことを昔書いた。ここで書いた考え方だと、三角形が退化した場合(三点が一直線上にある場合とか、一点に集中している場合)の扱いが面倒で、そのへんのことは、詳しく書いてない。が、最近、三角形の合同類のなす空間は、もっと単純な記述ができることに気付いた。


答えとしては、退化してない三角形の合同類の空間は$GL_{+}(2,\mathbf{R})/SO(2)$になる。実際、三角形の一点を原点に持ってきて、残りの二点の座標が、(a,b),(c,d)とすると、非退化な三角形であれば、ad-bc≠0が成立するから、そのような2点の配置の空間は$GL(2,\mathbf{R})$そのものである。一点を原点に置いてるので、これで並進の自由度は消えていて、残りは回転(と鏡映)の自由度で、商空間を考えれば、退化してない三角形の合同類の空間を与える。一般の次元でも、$GL_{+}(n,\mathbf{R})$(行列式が正のn行n列の行列全体)が$n$-単体の合同類の空間に推移的に作用して、固定部分群が$SO(n)$になるので、非退化な$n$-単体の合同類の空間は、$GL_{+}(n,\mathbf{R})/SO(n)$になる。この空間は、正定値対称行列全体と同一視できる($GL_{+}(n,\mathbf{R})$の元が、正定値対称行列と直交行列の積に一意に分解できるというのが、極分解定理)


#$SL(2,\mathbf{R})/SO(2)$は有名なポアンカレの上半平面になる。$GL_{+}(1,\mathbf{R})$は定数倍であるから、単なる拡大・縮小であり、$SO(n)$と合わせて相似な三角形全体を通る。$GL_{+}(2,\mathbf{R})$は$GL_{+}(1,\mathbf{R})$と$SL(2,\mathbf{R})/SO(2)$の直積であるから、$SL(2,\mathbf{R})/SO(2)$は三角形の相似類の空間である。$SL(n,\mathbf{R})/SO(n)$も対称空間の代表的な例である。数学的には、合同類より相似類の方が自然なのかもしれない


並進の自由度をad hocに消しているけど、$n$-単体の一様並進も入れると、$n$-単体の合同類の空間に、アフィン変換群が推移的に作用し、その場合、固定部分群は、合同変換群$Iso(n,\mathbf{R})$となる。回転と鏡映を考えると$GL(n,\mathbf{R})$が作用し、固定部分群は$O(n)$になる。


これを考えたのは、有限要素法の大変形解析のことを考えていた時。有限要素法では物体の変形がアフィン変換(並進の自由度を消せば一般線形変換)で線形補間すると十分よい近似になる微小領域を使う。この微小領域の初期状態(平衡状態)として、三角形や四面体を取ると、線形変換やアフィン変換によって、任意の三角形や四面体が得られる。そして、微小領域の形状のみを考えるなら、合同変換の自由度(剛体並進・剛体回転)は邪魔である。これは、$n$-単体の合同類の空間を考えた時と同じ考え方になっている。結局、(局所的な)有限歪みは、$GL(n,\mathbf{R})/O(n)$で分類され、超弾性体などの場合、弾性エネルギー/ひずみエネルギーは、この空間上の関数とみなすことができる。次元は、n^2-n(n-1)/2=n(n+1)/2で、n=2なら3、n=3なら6で、これは三角形や四面体の辺の数と等しい


※)三角形や四面体の定ひずみ要素を考える場合、平衡状態(要素の歪みエネルギーが最小となる配置)に於ける節点の座標と、変形後の節点の座標から、前者を後者へ変換するアフィン変換が一意に定まり、これの線形変換部分が変形勾配テンソル(テンソルと名前が付いてるけど、単なる行列)となる。変形勾配テンソルは直接、極分解するか、あるいは、転置行列との積を取ることにより、剛体回転成分が除かれ、有限歪みを得ることができる。応力と歪みの構成則があれば、応力成分と歪み成分の積を足し合わせて、弾性エネルギー密度を得られる。こうして得られた弾性エネルギーは、変形後の節点座標の関数と見なすことができるので、変形後の各節点の座標成分で微分すると、節点力が定まる


連続体が非圧縮性条件を満たす場合、局所的には体積が保存するので、$GL_{+}(n,\mathbf{R})$の代わりに$SL(n,\mathbf{R})$を考えることになる。従って、非圧縮性有限歪みの分類空間は、対称空間$SL(n,\mathbf{R})/SO(n)$で記述される(という事実が、何かの役に立つわけではない)

日本語/英語の組み合わせ範疇文法パーサーを作った

toyccg

https://github.com/vertexoperator/toyccg


英語については、CCG(Combinatory Categorial Grammar)のparserは公開されているものが、いくつか存在する(疑問文とか命令文は、あんまり対応してなかったりする。toyccgも疑問文や命令文は、あんまり対応できてない。やればできるとは思う)けど、日本語については、知る限りではない(公開せずに実装している人はいると思うけど)ので。原理的に新しいことは何もないけど、実用的に使うのに、必要そうな機能は、一通り実装してあると思う


自然言語用のパーサーとしては

・pythonで1500〜2000行程度((自然)言語非依存のコア部分は、1000行くらい。外部ライブラリは使ってないし他(プログラミング)言語への移植は楽と思う)で比較的小さい

・日本語と英語同時対応

・機械学習フリー

などの特徴がある。機械学習を使用してないのは、そもそも、日本語の学習用データが存在しない(多分)(英語用では、CCGBankがあって、普通これを使うと思うけど、どうやって入手するのか謎なので、やっぱり使えない)のもあるし、機械学習頼りだと、手で修正するということがやりづらくなるという面もある。手でルール作るとか、時代の流れに逆行している気もするけど、実用上は悪いことばかりでもない(CCGBankはTreeBankからの自動生成っぽいけど、学習用データを自分で作る場合、やっぱり手間がかかるわけだし)

#CFGを確率化した確率文脈自由文法(PCFG)と同様、CCGを確率化した、確率的組み合わせ範疇文法(PCCG)のような方法も存在する


nltkのCCGパーサーとの違いは

・デフォルトで辞書を用意している

・normal form parsingをサポートしている

・nltkのパーサーで日本語を解析する場合、単語単位への分かち書きは別に行う必要があるけど、toyccgは、分かち書きと構文解析を同時にやる

というあたり。分かち書き対応したおかげで、英語でも、New YorkとかWhite Houseみたいなのを一単語のようにして扱うことができるようになった(やろうと思えば、"This is a pen"のような文単位で統語範疇を辞書に登録することもできる)。あとまぁ、長単位・短単位みたいなのを気にする必要も、あまりない


#"数学的構造"みたいなのを辞書登録すると、数学+的+構造で分解するのと、"数学的構造"で単語扱いする結果が、ダブって出現することになる。こういう単語が多い文では、計算時間の増大にも繋がるし、本質的に同じ解析結果が沢山出現することにもなる。toyccgで、このパターンで計算量が増えるケースは結構ある


#CCGでは、長単位・短単位以外に、品詞についても、考える必要がない(これをメリットと捉えるかは人によるかもしれないけど)。名詞・動詞は分かりやすいけど、実際に作業していて、日本語文法の知識が殆どない私には、助詞・助動詞の区別は迷うことがあるし、例えば、もう死語のような気がするけど

「なう」の品詞ってなんですか?

http://togetter.com/li/258854

という話でも、CCGだと、"S\N"みたいな統語範疇を設定してやれば"渋谷なう"みたいな文章を解析できる。これに類似した統語範疇(つまり、"S[*]\N[*]")を持つ単語は、今のところ他にはない。動詞や形容詞は統語範疇として"S\NP"を持つので、近いと言えるかもしれない。"ご飯食べてるなう"だと、"S[now]\S"とするとか。統語素性[now]を付けているのは、"ご飯食べてるなうなう"みたいな文章を是としないため。



統語素性の設計について。日本語の統語素性は、ベースになるようなものがないので、勝手に作ったものが沢山ある(一応、「日本語文法の形式理論」という本が一冊出ていて、それは参考にした)。英語の場合は、CCGBankがデファクトスタンダードになっているので、なるべくそれに従っているけれど、改変してある部分もある(例えば、toyccgは前置詞句PPを使わない)。これらは、辞書の問題なので、辞書を差しかえれば、実装の方を変えることなく、他の統語範疇・統語素性設計に対応できる(はず)

#CCGBankは、結構謎のunary rulesがあって、本来のCCGからは、はみ出しているので、これが気に入らない。toyccgも、このへんの問題はクリアできてなくて、適切に統語範疇・統語素性を設計すれば解決するのか、それともどうしようもないのかは、不明。英語ではCCGBankが標準のようになっているけど、本当に、これでいいのかはよく検討するべきことと個人的には思う


unary rulesとして、どういうものがあるかは、

https://github.com/mikelewis0/easyccg/blob/e42d58e08eb2a86593d52f730c5afe222e939781/training/unaryRules

などを見ると分かる



機械学習ベースじゃないことのデメリットとして、構文解析できないケースが、増えてしまう(実用上は、何でもいいから、とりあえず結果を出してくれた方が嬉しい時もあれば、ダメな時には、あからさまに失敗してくれた方がいい時もあると思うので、一概にデメリットとも言い切れない。toyccgの場合、日本語として、文法が崩壊している文章は、解析されない可能性が高い)。一応、それなりには頑張っていて

>>> import toyccg.japanese as jpn
>>> jpn.run(u"すももももももももの内。")
test run : sentence=すももももももももの内。
すもも	N
も	(NP[sbj]/NP[nom-enum])\N
もも	N
も	NP[nom-enum]\N
もも	N
の	((S[nom]\NP[sbj])/(S[nom]\NP[sbj]))\N
内	S[nom]\NP[sbj]
。	ROOT\S[nom]

>>> jpn.run(u"象は鼻が長い。")
test run : sentence=象は鼻が長い。
象	N
は	NP[sbj]\N
鼻	N
が	NP[ga-acc]\N
長い	(S\NP[sbj])\NP[ga-acc]
。	ROOT\S

くらいの解析はできる。sentences.ja.txtには、解析できる文章として、どれくらいのものがあるか、100以上の例がある。一応、ニュースや論文に出てくるような、比較的固めの文章なら、そこそこ解析できるように頑張っている(解析できるとは言っていない)。twitterや話し言葉で出てくるような文章となると、なかなか難しい。これらの文章では、格助詞の省略がよく起きて、難しい原因の一つとなる。その他、未知語が含まれている場合(これは、ニュースや論文でも、よくある。一応、文字種に基づく未知語推定を行えるしているが)や、倒置法とかは全く対応されてない


助詞が省略されていても、

>>> jpn.run(u"お腹すいた。")
test run : sentence=お腹すいた。

>>> jpn.parser.lexicon.setdefault(u"お腹",[]).append( "NP[sbj]" )
>>> jpn.run(u"お腹すいた。")
test run : sentence=お腹すいた。
お腹	NP[sbj]
すい	IV[euph]
た	(S\NP[sbj])\IV[euph]
。	ROOT\S

みたいに、適切な統語範疇を辞書登録してやれば、対応できないわけではないけど。安易にこういうことをして、間違った構文解析結果が増えてしまうと、それはそれで困る。


#既存部分への影響を、頭を使わずに、確実に抑えたいなら、NP[x]とか適当な新しい素性を導入して、動詞や形容詞のような単語にも、対応する統語範疇("S\NP[x]"など)を追加していくなどの手もある。スマートな解決策とは言えないけど、実用的な解析を行うには、そういうことも必要かもしれない。日本語は、話し言葉レベルでは、「私、昨日、学校行った」みたいな格助詞が全部落ちたような文章も可能であるから、そこまで行くと、現時点では、どれが主語か正しく判定するのは難しい



jpn.runは、形態素解析的な結果しか出さないけど、内部的には、構文木が作られている

>>> import toyccg.japanese as jpn
>>> r = jpn.parser.parse(u"あらゆる現実を、全て自分の方へ捻じ曲げたのだ。")
>>> t = r.next()
>>> print(t.show())
(LApp (LApp (SkipCommaJP (LApp (RApp [あらゆる:N/N] [現実:N]) [を:NP[obj]\N]) [、:COMMA]) (RBx [全て:S[null]/S[null]] (RBx (LApp (RApp (LApp [自分:N] [の:(N/N[base])\N]) [方:N[base]]) [へ:(S[null]/S[null])\N]) (LB (LApp [捻じ曲げ:TV[euph]] [た:(S[null]\NP[obj])\TV[euph]]) [のだ:S[null]\S[null]])))) [。:ROOT\S[null]])

あんまり見やすい出力ではないけど。多くの場合、構文解析結果は一意でないので、結果はジェネレータで、一つずつ返ってくる。

#デフォルトでは、ソースコード中で、shortcutと書いてある"近似"が有効になっていて、例えば、"かわいい私の妹"は、二種類の係り受け構造("かわいい"が"私"に係る場合と、"妹"に係る場合)が存在するが、どちらの場合も統語範疇は"N"となる(これは、本質的に構造が違うので、normal form parsingでは対処できない)ので、それ以降の解析には影響しない。こういうのが沢山あると、長文では、解析結果が沢山出てきて、計算時間が膨大になる。こういう状況を回避するため、先に出てきた結果だけを取って、それ以降のは使わないという処理をしている。最終的な構文木が一つしか必要ない場合は、shortcutをTrueにしておく方が、計算が速く済む(大体、ニュースとか論文とかの文章だと、これがないと一文の計算が全然終わらないことが多い)けど、本当に、全ての結果が必要な場合は、shortcutをFalseにする必要がある


分かち書きしたいだけなら、

>>> print [c.token for c in t.leaves()]
[u'\u3042\u3089\u3086\u308b', u'\u73fe\u5b9f', u'\u3092', u'\u3001', u'\u5168\u3066', u'\u81ea\u5206', u'\u306e', u'\u65b9', u'\u3078', u'\u637b\u3058\u66f2\u3052', u'\u305f', u'\u306e', u'\u3060', u'\u3002']

とかできる



感想。色々と実験するプロトタイプ的なつもりだったのだけど、思ったよりも、よく解析できるな、という印象(自画自賛)。


課題

・複合名詞対応

・格助詞落ち対応


今後。もう少し改善しようかと思ってたけど、飽きてしまったので、当分いじることはなさげ(完)

2016-01-03

非可換Groebner基底と表現論(1)概要メモ

Groebner基底は、多項式環のイデアルの"よい"生成元のことで、Groebner基底を計算するBuchbergerアルゴリズムなどの存在によって、計算機代数で広く使われている。Groebner基底を非可換多項式環に拡張しようというのは、誰でも考えるようなことではあると思うけど、1986年に、Moraという人が、非可換多項式環の両側イデアルに対する非可換Groebner基底を計算するように、Buchbergerアルゴリズムを拡張したらしい(これは、停止するとは限らない半アルゴリズム)。


ところで、非可換代数では、両側イデアルと左イデアル/右イデアルの区別が生じる。Lie環の表現論の文脈では、普遍展開環の定義関係式は、非可換多項式環の両側イデアルを与え、普遍展開環の左イデアルによる商加群を考えると、Lie環の表現が得られる(抽象論としては、unitalな非可換結合代数$R$の単純加群は、ある極大左イデアル$I$によって、$R/I$の形で書ける)。そういうわけで、表現論の主役は、むしろ左イデアルであるから、両側イデアルの非可換Groebner基底が分かっても、特別嬉しいことはない。


などと思っていたところ

GROEBNER-SHIRSHOV BASES FOR REPRESENTATION THEORY

http://www.math.uconn.edu/~khlee/Papers/GS-pair.pdf

という論文を読んだ。もしかしたら、これ以前に同じことを考えた人はいるかもしれないけど、確認はしてない(求情報)。タイトルにあるGroebner-Shirshov basisというのは、非可換Groebner基底と同じものだ(と思う)けど、上記論文では、Groebner-Shirshov pairという両側イデアルの"よい"生成元と左イデアルの"よい"生成元の組を定義している。論文では、Lie環sl(3)の有限次元既約表現について、Groebner-Shirshov pairを与えている(全ての有限次元既約表現に対する一般形を知るため、証明を行っているけど、Buchbergerのアルゴリズムの拡張として、Groebner-Shirshov pairを計算する半アルゴリズムも与えている)


#その後の論文で、Lie環sl(n)の有限次元既約表現についても具体形を与えている

Grobner-Shirshov Bases for Irreducible sl_{n + 1}-Modules

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

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



Groebner-Shirshov pairが分かると、何ができるのかというと、イデアルメンバーシップ問題を解いたり、イデアルが等しいかどうかの判定を行ったりできるはずだけど、それは表現論的には嬉しい場面がなさそうな気がする。応用の一つとして、標準単項式の計算がある。これは、Groebner基底の時にもあったものであるけど、標準単項式は表現の基底を与える(Groebner"基底"とかGroebner-Shirshov"基底"は、生成元であるけど、こっちは表現空間の基底を指す)。ある代数の表現論で、具体的な計算を実行しようと思うと、表現の"よい"基底を取るということは、重要な問題となる。この種の問題に対する答えとして、Gelfand-Tsetlin基底や結晶基底があるけども、その存在や構成は、代数や表現の性質に強く依存している。標準単項式は、より一般的な代数に対して定義可能なものとなる


#Groebner-Shirshov pairの前に、左イデアルの生成元が分かっている必要はあるわけだけど、例えば、複素半単純Lie環の有限次元既約表現の場合、これは既約最高ウェイト表現であり、特異ベクトルがどこに出るかも全部分かるので、左イデアルの生成元は分かる。既約表現の構成の仕方によっては、左イデアルの生成元を調べること自体が難しい場合もあるとは思うし、最高ウェイト表現でない場合、一般にどれくらいのことが知られているのかよく分からない(Casimir作用素=普遍展開環の中心の元は、既約表現に定数倍で作用するので、それによって、左イデアルの生成元がいくつか得られる場合もある。それ以外に、どれくらいの付加的な条件があるのかは、分からない)



とりあえず、証明とかは抜きで、ざっと流れを把握する。可換Groebner基底の時と殆ど同じではあるけど、detailは、非可換性のために多少複雑になっている


まず単項式順序。非可換の時は、length-lexicographic order(長いのでdeglex順序と呼ぶことにする)というもの以外使われてるのを見ない(これしか適切な単項式順序がないのかどうかは知らない。要検討)。ちゃんとした定義は上論文参照。2変数x,yでx<yの場合、1 < x < y < x^2 < xy < yx < y^2 < x^3 < x^2y < xyx < xy^2 < yx^2 < ....などとなる。それで、可換の時と同様、leading monomialとかleading coefficientが定まる。leading coefficientが1の時、monicであるという


normal formもしくは剰余。可換の時と同様、剰余によって、normal formを計算/定義する。これが一意に定まるのが、非可換Groebner基底やGroebner-Shirshov pairで重要な性質(定義としてもいい)。$S$と$T$をそれぞれ、monicな非可換多項式の有限集合として($S$を両側イデアルの生成元、$T$を左イデアルの生成元と思う)、非可換多項式$p$が

p = ¥sum ¥alpha_i (a_i s_i b_i) + r

と書けるとする。但し

(1)¥alpha_iは定数(複素数とか有理数)

(2)a_i,b_iは任意の単項式、s_i ¥in Sで、a_i LT(s_i) b_iLT(p)を超えない($LT(f)$は$f$のleading monomial)

(3)$r$の任意の単項式は、任意の$S$の元$s$について$a LT(s) b$と書けないとする

このような表示は、剰余の手続きによって得ることができる($a_i$,$b_i$は単項式なので、和は有限和であれば、異なる$i,j$について$s_i=s_j$となってもよい)。この時$r$を$p$の$S$に関するnormal formとする。$S$に関するnormal formが任意の非可換多項式$p$について一意に定まる時、$S$の元は、非可換Groebner基底をなすという。


同様に$(S,T)$に関するnormal formを定める。Groebner-Shirshov pairの場合、両側剰余と左剰余が混ざることに注意する必要がある。非可換多項式$p$が

p = ¥sum ¥alpha_i (a_i s_i b_i) + ¥sum ¥beta_j c_j t_j + r

と書けるとする。但し

(1)¥alpha_i,¥beta_jは定数

(2)a_i,b_i,c_jは任意の単項式、s_i ¥in S , t_j ¥in T

(3)$r$の任意の単項式は、任意の$S$の元$s$についてa LT(s) bと書けず、任意の$T$の元$t$についてc LT(t)とも書けないとする

この時$r$を$p$の$(S,T)$に関するnormal formとする。$(S,T)$に関するnormal formが任意の非可換多項式$p$について一意に定まる時、$(S,T)$はGroebner-Shirshov pairであるという。


ちゃんと書くと式は何か複雑であるけど、$S$は両側イデアルの生成元であるので、両側剰余を取っている。$T$は($S$が生成する両側イデアルによる剰余代数の)左イデアルの生成元となるものなので、左側からの掛け算しか許していない


非可換Groebner基底もGroebner-Shirshov pairも、Buchbergerアルゴリズムを拡張した完備化手続きによって計算できる。注意することとして、非可換Groebner基底の場合、可換の場合のS-polynomialに相当するものが一つとは限らないという事情がある。例えば、f = xy - x , g = yx - yでx<yとするdeglex順序が入っているとすると、fとgのleading monominalはそれぞれxyとyxになる。可換の時、これらが消えるようにLCMを選んでいたけども、fx - xg = -x^2 + xyで(-x^2+xy=-x^2+(f+x)なので-x^2+xにreductionされることに注意)、一方、yf - gy = -yx + y^2でもある(-yx+y^2=-(g+y)+y^2なのでy^2-yにreductionされる)。こうして、次は4つの多項式{f,g,x^2-x,y^2-y}の各ペアに対して、同様の計算を繰り返す(この例では、これで終わりとなる)。ついでに停止性も保証されない


$(S,T)$のGroebner-Shirshov pairを計算する場合は、まず$S$の非可換Groebner基底を計算して、その後、$T$の方を完備化していくことになる。$T$の方の完備化は、少し込み入っていて、2ステップの繰り返しとなる。まず、$T$同士の元で、right-justified composition(これも、S-polynomialの計算の一般化)というものを計算して簡約、新しい元の追加を行う。次に、完備化した$S$の元とupdateした$T$の元の対で、$T$に元を追加していく。これを繰り返す(詳細は論文参照。言葉だけで説明するのが辛くなってきた)。


アルゴリズムの正しさは、非可換Groebner基底/Groebner-Shirshv pairでも、剰余の一意性と"closed under composition"という性質(可換の場合の、全てのS-polynomialを簡約すると0になるという性質に対応)が同値であることを示せば自明となる


例。Lie環$sl(3)$の最高ウェイト¥Lambda_1+¥Lambda_2の既約表現を考える。$sl(3)$の下三角部分しか重要でないので、Chevalley生成元$f_1,f_2$のみを考える。単項式順序は$f_1 < f_2$で定まるdeglex順序とする。これらには2つのSerre関係式がある

p = f_2^2f_1 - 2 f_2f_1f_2 + f_1f_2^2

q = f_2f_1^2 - 2 f_1f_2f_1 + f_1^2f_2

が両側イデアルの方の生成元。唯一の"S-polynomial"p f_1 - f_2 qが存在するけども、$\{p,q\}$による剰余を計算すると0になるので、$S=\{p,q\}$は非可換Groebner基底をなしている。


一方、片側イデアルの方は、2つの非可換多項式

T=¥{t_1 = f_1^2 , t_2 = f_2^2¥}

で生成される。$T$の完備化の第一段階であるright-justified compositionは、今回することがない。次に、$S$と$T$の間の完備化計算を行うと、次の3つのS-polynomialが出る($S$と$T$による剰余も計算している)

 t_3 = q - f_2 t_1 = -2 f_1 f_2 f_1 + f_1^2 f_2

 t_4 = p f_1 - f_2^2 t_1 = -2 f_2 f_1 f_2 f_1 + 2 f_1 f_2 f _1 f_2 - f_1^2 t_2 ¥equiv -2 f_2 f_1 f_2 f_1 + 2 f_1 f_2 f _1 f_2

 t_5 = q f_1 - f_2 f_1 t_1 = -3 f_1 f_2 t_1/2 - f_1 q/2 + f_1^3 f_2/2 ¥equiv f_1^3 f_2/2

これら3つの多項式+元の$T$の2つの多項式について、"right-justified completion"を行う。のだけど、その前に、2つ目の多項式は、

 f_2 t_3 - q f_2 + f_1^2 t_2

なので、除いてしまっても問題ない。ということを繰り返すのだけど、手計算では、予想より大分辛かったので以下略


というわけで、何はともあれ、多くのことが可換Groebner基底の時とパラレルに進む


TODO:あとで、Groebner-Shirshov pairを計算するコードを追加する

2015-11-06

Combinatory Categorical Grammar

というものを最近調べた。前身として、Categorical Grammar(範疇文法)というものがあったらしい。


歴史的には、範疇文法は文脈自由文法(CFG)より古いっぽいが、範疇文法は、文脈自由文法と同等の計算能力しかないことが知られているらしい。基本的には、

S -> NP VP

のような生成規則があったら、VP=S\NPないしNP=S/VPとして、前者であれば、VPは、"左に型NPの引数を受け取って、型Sの値を返す関数"と見なし、後者なら、NPは、"右に型VPの引数を受け取って型Sの値を返す関数"と見なす。一般的なプログラム言語で、普通の中置きされた(int型同士の)"+"演算子であれば、(int\int)/intのような型になる感じだと理解すると、分かりやすい気がする。


CFGと比較した場合

・CFGと比べるとラムダ計算で意味付けする方法が明確なので、意味論に繋げやすい

・各単語ごとに、"型"が付いて、統語規則の方は変更しない語彙化文法というものになっている(実用上のメリットとして、CFGで文法ルールを作って、生成規則を増やす場合、影響範囲を予測するのが難しいけど、語彙化文法なら、追加した語彙が出てこない文には影響しない)

というメリットがある。CFGから範疇文法への変換は、一回標準形にする必要があり、例えば、二重他動詞DTVは、CFGでは生成規則VP->DTV NP NPを持つと思われるけど、直接対応する範疇文法は作れない(範疇文法なら、DTV=(VP/NP)/NPという形にするのが普通)。


NPやVPやSは、型とは呼ばずに、(統語)範疇と呼ぶのが普通ぽい。時制や人称のような情報は、統語素性として扱って、単数形名詞をN[sg]、複数形名詞をN[pl]などと書いたりする(場合もある)。この場合、実用上は、N[sg]とN[pl]は、特に何の関係もない、別の素性という扱いになる。現実的に、文法的に許される文章を過不足なく記述しようと思うと、かなり多くの統語素性は必要になると思う(時制や人称がないと、"He meet ..."とか"him looks ..."みたいな文章も合法となったりする)。



組み合わせ範疇文法(CCG)。範疇文法には関数適用しかないけど、ラムダ計算で意味付けされることを考えると、関数合成を許すのは自然に思える。関数適用が左と右の2つあるので、合成のバリエーションも多い。関数合成があると、

He conjectured and might prove completeness.

のような文章で、

conjectured |- (S\NP)

might |- (S\NP)/(S\NP)

prove |- (S\NP)

から、

might prove |- S\NP

とできて、conjecturedとmight provedが同じ型になって対をなすようにできる。andは、forall x,(x\x)/xという多相関数みたいな型として実現されるとする。他に、

(X/Y)/Z Y/Z => X/Z

(X/Y)\Z (Y\Z) => X\Z

Y\Z (X\Y)\Z => X\Z

Y/Z (X\Y)/Z => X/Z

のような代入規則も、CCGでは使われる


ここまでは、いわば"二項演算"であるけども、CCGには、もう一つ、type raisingというのもあって、これは、一項演算みたいなもの。インフォーマルには、S->NP VPをCG/CCGに変換するとき、VP=S\NPとしても、NP=S/VPとしてもよいという任意性があり(意味論の自然さから、普通は、NPを基本範疇にして、VP=S\NPとする)、VP=S\NPと、NP=S/VPを組み合わせると、NP=S/(S\NP)と"解ける"。一般に、XとY=T\Xがあれば、X=T/Yでもあり、X=T/(T\X)と解ける。これがtype raisingの由来と思う

#左右の区別を忘れると、Type raisingは、イメージ的には、型Aから、多相型forall X,(A->X)->Xを作り出すような操作と見なせる。一方、多相型ラムダ計算では、最小不動点のChurch Encodingとして、forall X,(F(X) -> X) -> Xというのがあって、その特殊な場合と思える。

#日本語(トルコ語とかもそうらしい)などは、語順の自由度が英語より高い(と一般的に言われる)が、それらを系統的に扱うために、scramblingというものも追加する場合があるらしい


type raisingの恩恵として、関係代名詞の扱いが自然になるという説明が、よく書いてある。例えば、the company which Microsoft bought

だと、(boughtは他動詞なので)

Microsoft |- NP

bought |- (S\NP)/NP

だけど、type raisingによって

Microsoft |- S/(S\NP)

となって、関数合成によって

Microsoft bought |- S/NP

となって、

which |- (N\N)/(S\NP)

が適用できる形になる



意味が扱えるという話。例として、"0"と"1"と"+"だけで構成される算術式を、範疇文法(<=>文脈自由文法)で考える(個人的な感覚としては、CCGは、数式やプログラミング言語のような形式言語を扱うには、自由度が高すぎて、制御しづらいという気がする)。この場合、"0","1","+"の範疇として

"0" |- S : 0

"0" |- S\S : (\x -> 10*x)

"1" |- S : 1

"1" |- S\S : (\x -> 10*x + 1)

"+" |- (S\S)/S : (\x y -> x+y)

があるとする("001"のような表記も許容する)。':'の後ろにあるのが、意味。範疇文法なので、右適用と左適用だけ考えればよく、右適用を>で、左適用を<で表すと、"110+1"の導出木は、例えば

(< (< (< "1"|S "1"|(S\S)) "0"|(S\S)) (> "+"|( (S\S)/S ) "1"|S))

みたいになる。対応する意味を

(< (< (< 1 (\x -> 10*x + 1)) (\x -> 10*x)) (> (\x y -> x+y) 1))

のように置き換えて、簡約化していくと

(< (< (< 1 (\x -> 10*x + 1)) (\x -> 10*x)) (> (\x y -> x+y) 1))

=> (< (< 11 (\x -> 10*x)) (> (\x y -> x+y) 1))

=> (< 110 (> (\x y -> x+y) 1))

=> (< 110 (\y -> 1+y))

=> 111

という風に計算できる。というわけで、"110+1"の意味は、今の場合、111となる。


どういう意味を付けるかは、CGやCCGとは独立の話で、JSONパーサーなんかだと、対応するJSONオブジェクトが意味ということになるし、英語の意味を表すのに、日本語を使ってもいい。自然言語形式意味論としては、discourse representation theory(DRT)というものが比較的よく使われているっぽい(あまり調べていない)。一階述語論理との違いがよく分からないけど、東ロボくんの数学部分は、問題を計算機代数システムで処理できる形に変換するのに、CCG+DRTを使っているらしい(と書いてある)

cf.深い言語理解と数式処理の接合による入試数学問題解答システム

https://kaigi.org/jsai/webprogram/2013/pdf/622.pdf



CCGはCFGより計算能力が高いことが知られている(文脈依存文法よりは弱いらしい)けど、CCGまでいくと、主に関数合成を許しているせいで、導出木は、一般に非常に沢山ある。が、意味の方で見ると同じになる擬似的曖昧性(spurious ambiguity。導出木は違うけど、意味は同じと見なせる)と、そもそも、異なる意味に対応する、本質的な曖昧性が存在する

(本質的曖昧性の例1)有名な例

I saw a girl with a telescope

で、withが動詞句"saw a girl"にかかるか、名詞句"a girl"にかかるか、で、withの範疇は、( (S\NP)\(S\NP) )/NPと(NP\NP)/NPになって、意味的に異なる

(本質的曖昧性の例2)

The fish twisted and turned on the bent hook.

で、"on"の範疇が、( (S\NP)\(S\NP) )/NPであるとすると、"on the bent hook"は、turnedのみにかかると見ることもできるし、"twisted and turned"にかかると見ることもできる(あくまで、文法的には)。というわけで、全ての範疇が決まっても、、意味の異なる導出木が複数存在する場合がある



慣用表現の扱い。特に体系的な扱いは存在しない。英語には、色々な慣用表現がある。特に、句動詞(phrasal verb)というものが多い。look at ...みたいなの。こういう時には、at ...がlookにかかるというより、look atが他動詞であるかのように見たい。こういう時は、"at ..."の範疇をPP(前置詞句)とする。atの範疇はPP/NPとなる。lookの範疇は(S\NP)/PPとなる。look atの範疇は合成によって(S\NP)/NPとなるので、目的は果たせられる。ただ、どこからが句動詞とするのか、明確な線引きはない気がする。take ... apartみたいなパターンを慣用表現として捉えたければ、apartの範疇を何か新しいXPとして、takeの範疇を(VP/XP)/NPとすれば、一応、そういう扱いもできる


別の慣用表現の例。"It is easy to show that."とか"There is ...."みたいな文で、ItやThereに、もはや意味はない。これらは、CCGBankなどでは、それぞれNP[expl]とか、NP[thr]という範疇を持ち、isは((S\NP[expl])/VP[to])/VP[adj]とか(S\NP[thr])/NPみたいな特殊な範疇となる。


実際の構文解析は、例えばCYK法で行うことができる。英語のCCG parserを手っ取り早く試すには、

EasyCCG

http://homepages.inf.ed.ac.uk/s1049478/easyccg.html

とかがいい気がする(多分、これはProbabilistic CCG parserで、CCGBankで学習してるのだと思う。命令文は、実質対応できてないっぽい?)(CCGBankはPenn Treebankから機械的に生成されたものだと思う)(何故か、自然言語parser業界では、Javaが使われることが多く、これもJavaで書かれている)。python nltkにもCCG parserはあるけど、lexiconを自分で定義する必要がある。


日本語のCCG parserは、自分で作るしかない気がする

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

Photons and gravitons in conformal field theory

http://link.springer.com/chapter/10.1007%2F3540171630_73

という論文の2ページ目の冒頭を見たら、次のようなことが書いてあった。



$V$を、Poincare代数のhelicity 0のmassless既約ユニタリ表現$V_0$と、Poincare代数の自然な4次元表現(並進成分は自明に働くので、実質的に$so(3,1)$のベクトル表現)のテンソル積表現空間とする。$V$は、ユニタリー表現ではないし、既約でもないが、直既約表現になっているらしい。$V$の重要な部分表現として、

$V_L$:Lorenz gauge条件を満たす$V$の部分空間

$V_g$:A_{¥mu}=¥partial_{¥mu} ¥Lambdaと書けるような$V$の元全体

があって、$V_g$は$V_L$の部分空間でもある。そして、以下の3つの命題が成立する

(1)商表現$V/V_L$は、Poincare代数のhelicity 0のmassless既約ユニタリ表現と同値

(2)商表現$V_L/V_g$は、Poincare代数のhelicity +1/-1の2つのmassless既約ユニタリ表現の直和表現と同値

(3)$V_g$は、Poincare代数のhelicity 0のmassless既約ユニタリ表現と同値


Gupta-Bleuler量子化の時は、非負ノルムを持つ状態の部分空間を0ノルム空間で割って、正ノルムを持つ状態を取り出していた。同じことを、1粒子状態に対してやっているのが、上の操作と考えられ、$(V,V_L,V_g)$を指して、Gupta-Bleuler tripletと呼んでいる文献もある。$V$には、平面波のような、ノルムが定まらない状態は含まれておらず、高々可算個の基底を持つ、"小さな"空間になっている


#重力波は、Klein-Gordon方程式を満たす、10成分の場である。Poincare代数の自然表現の2次の対称テンソル積表現は10次元表現で、これとhelicity 0のmassless既約ユニタリ表現のテンソル積によって、新しい表現を作ることができる。この表現から出発して、上記と類似の方法で、helicity +2/-2のmassless既約ユニタリ表現を得ることができるっぽい(要確認)


Poincare代数のhelicity +1/-1の2つのmassless既約ユニタリ表現の直和表現は、Maxwell方程式の解空間と同型になっているというのが、twistor理論の教えだった(多分、Maxwell方程式のノルム有限な解は、これ(の完備化したもの)で尽きると思うのだけど、証明とかは知らない)。この解空間は、$V_L/V_g$と同型であるが、一方、$V_L$の元は、電磁ポテンシャルなので、標準的な電磁気学の知識から、対応する電磁場を直接作ることができ、これはMaxwell方程式を満たす。真空中の電磁場は、波動方程式を満たすので、helicity 0のmassless既約ユニタリ表現と、Poincare代数の6次元表現(自然表現の2次の反対称テンソル積表現)のテンソル積空間$V_2$に含まれる解を考えるのは自然。そうすると、Poincare代数の作用と可換な写像

d_1 : V_L ¥to V_2

があって、その像$Im(d_1)$は、$V_2$の中で、Maxwell方程式を満たす元全体に等しいと期待される。そして

Ker(d_1) ¥simeq V_g

であれば、$Im(d_1)$は、(2)の表現空間$V_L/V_g$と同型になるので、辻褄が合う。


$V_0$は、共形代数$su(2,2)$の最低ウェイト表現に一意に拡張できるので、具体的な基底を取って、Poincare代数の作用を明示的に書くことができる。同様に、$V_L$や$V_g$の基底も具体的に記述することができるはず。ここでは、$V_0$を、以下のように、4変数多項式環の部分空間として実現する

V_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¥}

$su(2,2)$の複素化sl(4,¥mathbf{C})のChevalley生成元の作用は、

e_1 = z_1 ¥partial_2

f_1 = z_2 ¥partial_1

h_1 = z_1 ¥partial_1 - z_2 ¥partial_2

e_2 = - z_2 z_3

f_2 = ¥partial_2 ¥partial_3

h_2 = 1 + z_2 ¥partial_2 + z_3 ¥partial_3

e_3 = - z_4 ¥partial_3

f_3 = - z_3 ¥partial_4

h_3 = - z_3 ¥partial_3 + z_4 ¥partial_4

のようになる。1が最低ウェイトベクトルとなり、これに上記の元を作用させていくことで、$V_0$の基底が得られる。別に、この実現を使うことに特に意味はないけど、コンピュータに計算させるときに、ちょっと楽。変数の数はたまたま4つだけど、これ自体は、通常の時空の座標とは関係ない


#$V_0$に含まれる同次多項式の次数は、全て偶数で、次数2n(n=0,1,2,...)の同次多項式が張る部分空間の次元は、$(n+1)^2$になる。これは、$V_0$の定義によって、z1,z2の部分の次数と、z3,z4の部分の次数が等しくないといけないことから、容易に分かる。

e_4 = ¥[e_1,e_2¥] = -z_1 z_3

e_5 = ¥[e_2,e_3¥] = -z_2 z_4

e_6 = ¥[e_1,e_5¥] = -z_1 z_4

P_0 = -(e_2+e_6) = z_2 z_3 + z_1 z_4

P_1 = -¥sqrt{-1}(e_4 - e_5) = ¥sqrt{-1}(z_1z_3 - z_2 z_4)

P_2 = -(e_4+e_5) = z_1 z_3 + z_2 z_4

P_3 = e_6 - e_2 = z_2 z_3 - z_1 z_4

となる。P_0^2 - P_1^2 - P_2^2 - P_3^2は恒等的に0なので、質量0のKlein-Gordon方程式が満たされる。空間回転$M_1,M_2,M_3$及びLorentzブースト$L_1,L_2,L_3$は、

M_1 = -(z_1 ¥partial_2 - z_2 ¥partial_1 - z_3 ¥partial_4 + z_4 ¥partial_3)/2

M_2 = ¥sqrt{-1}(z_1 ¥partial_2 - z_3¥partial_4 - z_4 ¥partial_3 + z_2 ¥partial_1)/2

M_3 = -¥sqrt{-1}(z_1 ¥partial_1 - z_2 ¥partial_2 + z_3 ¥partial_3 - z_4 ¥partial_4)/2

L_1 = -¥sqrt{-1}(z_1 ¥partial_2 + z_3 ¥partial_4 - z_2 ¥partial_1 - z_4 ¥partial_3)/2

L_2 = (-z_4¥partial_3 - z_1 ¥partial_2 - z_2 ¥partial_1 - z_3 ¥partial_4)/2

L_3 = (z_1 ¥partial_1 - z_2 ¥partial_2 - z_3 ¥partial_3 + z_4 ¥partial_4)/2



以降は、普通の電磁気学でやる計算と同じになる。

V = V_0 ¥otimes ¥mathbf{C}^4として、

$w_0 = (1,0,0,0)$

$w_1 = (0,1,0,0)$

$w_2 = (0,0,1,0)$

$w_3 = (0,0,0,1)$

を$\mathbf{C}^4$の基底とする

$d_0 : V_0 \to V$を

d_0(f) = (e_2+e_6)(f) ¥otimes w_0 - ¥sqrt{-1} (e_4-e_5)(f) ¥otimes w_1 - (e_4+e_5)(f) ¥otimes w_2 -(e_2-e_6)(f) ¥otimes w_3 = (-P_0(f) , P_1(f) , P_2(f) , P_3(f))

で定義し、

$V_g = Im(d_0)$とする。$d_0$は、Poincare代数の作用と可換なはずで、従って、$V_g$は$V_0$と同値な表現空間となる。$V_g$の最低ウェイトベクトルは

(z_2 z_3+z_1z_4 , ¥sqrt{-1}(z_1z_3 - z_2 z_4) , z_1 z_3 + z_2 z_4 , z_2 z_3 - z_1 z_4)

になる。第一成分がスカラーポテンシャル、残りがベクトルポテンシャルに相当する


$d_0$とPoincare代数の作用の可換性は、例えば

(1 ¥otimes M_1 + M_1 ¥otimes 1) d_0 = (-M_1 P_0 , M_1 P_1 , M_1 P_2-P_3 , M_1 P_3+P_2) = (-P_0 M_1, P_1 M_1 , P_2 M_1 , P_3 M_1)

などで、この場合、

¥[M_1,P_2¥] = P_3

¥[M_1,P_3¥] = -P_2

より従う。他の生成元についても、同様に計算すればいい


次に、d_1 : V_0 ¥otimes ¥mathbf{C}^4 ¥to V_2は、通常の外微分と同様、A = (A_0,A_1,A_2,A_3) ¥in V_0 ¥otimes ¥mathbf{C}^4に対して

d_1(A) = d_0(A_0) ¥wedge w_0 + d_0(A_1) ¥wedge w_1 + d_0(A_2) ¥wedge w_2 + d_0(A_3) ¥wedge w_3

で定義する。6成分出てくるけど、当然、電場と磁場に対応する


Mikowski計量に対する、Hodgeのスター作用素* : ¥wedge^2 ¥mathbf{C}^4 ¥to ¥wedge^2 ¥mathbf{C}^4

*(w_0 ¥wedge w_1) = -w_2 ¥wedge w_3

*(w_0 ¥wedge w_2) = -w_3 ¥wedge w_1

*(w_0 ¥wedge w_3) = -w_1 ¥wedge w_2

*(w_1 ¥wedge w_2) = w_0 ¥wedge w_3

*(w_3 ¥wedge w_1) = w_0 ¥wedge w_2

*(w_2 ¥wedge w_3) = w_0 ¥wedge w_1

で定義すると、$*^2=-1$なので、¥wedge^2 ¥mathbf{C}^4を$*$の固有値+i/-iの空間に固有分解できる。$*$は、Poincare代数の作用と可換なので、$V_2$も2つの表現空間の直和に分解する。この分解は、$Im(d_1)$に制限すると、Maxwell方程式の解空間が、正helicity/負helicityの解に分解する(それぞれ、自己双対・反自己双対Maxwell方程式の解空間)のと対応している。


自己双対・反自己双対Maxwell方程式の解空間は、共形代数の$su(2,2)$の最低ウェイト表現空間となっているので、最低ウェイトベクトルを計算するのが、構造を知る上で、手っ取り早い。$V_2$には、$su(2,2)$自体は作用しないが、Poincare代数の作用のみでも、最低ウェイトベクトル条件のいくつかは書くことができ、かつ、最低ウェイトベクトルを決定するにはそれで十分である。このへんは、表現論の知識が必要となる。最低ウェイト条件の詳細は

Maxwell方程式の表現論

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

で書いたものと同じ。


具体的に、自己双対Maxwell方程式の"最低ウェイト解"を計算してみる。まず、$*$について、固有値+iの空間の基底として、

w_0 ¥wedge w_1 + ¥sqrt{-1} w_2 ¥wedge w_3

w_0 ¥wedge w_2 + ¥sqrt{-1} w_3 ¥wedge w_1

w_0 ¥wedge w_3 + ¥sqrt{-1} w_1 ¥wedge w_2

を取ることができる。この基底に対する、成分を(¥psi_1,¥psi_2,¥psi_3) ¥in V_0 ¥otimes ¥mathbf{C}^3とする。$V_0$の定義より

(z_1 ¥partial_1 + z_2 ¥partial_2 - z_3 ¥partial_3 - z_4 ¥partial_4)¥psi_i = 0

が分かる。$su(2,2)$作用のChevalley生成元のいくつかは、Poincare代数の元のみで書けて、例えば

h_1 = L_3 - ¥sqrt{-1} M_3 = z_1 ¥partial_1 - z_2 ¥partial_2

の場合、

h_1(¥psi_i) + 2 ¥psi_i = 0

が最低ウェイトベクトルの満たす条件となる。また、

f_1 = (M_1 - ¥sqrt{-1}M_2)/2 - ¥sqrt{-1}(L_1 - ¥sqrt{-1}L_2)/2 = z_2 ¥partial_1

に対して、

f_1(¥psi_i) = 0

なども条件としてある。ここまでの条件で、¥psi_1,¥psi_2,¥psi_3

z_2^2 z_3 z_4 , z_2^2z_3^2, z_2^2z_4^2

の3つの元の一次結合で書けることが分かる。他に、$h_2,h_3,f_3$などに関する条件も使うと、定数倍を除いて

¥psi_1 = -¥frac{¥sqrt{-1}}{2}(z_2^2 z_4^2 + z_2^2 z_3^2)

¥psi_2 = ¥frac{1}{2}(z_2^2 z_4^2 - z_2^2 z_3^2)

¥psi_3 = z_2^2 z_3 z_4

という解を得る。これらの解を与える電磁ポテンシャルとして、例えば

(¥phi,A_x,A_y,A_z) = (-z_2 z_4/2 , ¥sqrt{-1}z_2 z_3/2 , z_2 z_3/2 , -z_2 z_4/2)

などが取れる。反自己双対Maxwell方程式の最低ウェイト解も同様に計算できる。


光子の状態空間は、光子の1粒子状態の対称テンソル積空間として得られる。物理の教科書では、波数ごとに、(非可算無限個の)生成消滅演算子があって、それを真空に作用させて...ということになるが、平面波は、Maxwell方程式のノルム有限な解ではなく、$P_0,P_1,P_2,P_3$の連続スペクトルに対する一般化固有状態として得られるはず。そして、($V_0$の各基底に対応して)高々可算無限個の生成消滅演算子があって、それを真空に作用させていくことで、多光子状態が得られるという形になる。通常、物理の教科書に書いてある、量子電気力学の生成消滅演算子による光子の状態空間と、表現論的("群論的")な光子の一粒子状態空間の説明が、等価だというのは、長い間納得できなかったけど、数学的に正当化できそうではある



電磁気学をやるには、電子の状態空間も知りたい。Poincare代数のmassless既約表現については、共形代数の表現に拡張することで、その構造が数学的によく分かる(通常の群の誘導表現による構成では、詳しい構造を調べるのは結構難しい)のだけど、massive表現については、その手は使えない。が、de Sitter代数$so(4,1)$の主系列表現は、$so(4,1)$からPoincare代数へのcontraction(半径→∞の極限を取る)を行うことによって、Poincare代数のmassive既約ユニタリ表現2つの直和に分解するらしい。

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

http://link.springer.com/chapter/10.1007/978-94-010-2669-7_8

$so(4,1)$の主系列表現は、最高ウェイトも最低ウェイトもないので、Lie環の表現論としてはまだ難しい部類ではある(というのは、現時点では、系統的な手法が使えず、個別に対処する必要があるという意味で)けども、Gelfand-Zetlin basisなど、具体的な基底を構成することができる(はず←やったことない)。これによって、Poincare代数のmassive既約表現の具体的な記述を手にすることができるはず(続く)


#contractionによる構成自体は、Euclid群の既約ユニタリ表現を作る場合にも、同じような方法が取られるので、さほど問題なく出来るはずではある。masslessの方の基底は分かっても、massiveの方がどうにもならなくて、ふと、Eudlic群なら、誰かやってる人いそうだ、と思って調べたら、q-Euclid代数でやってる人がいて、Poincare代数でも、同じことを考えた人が1960年代にいたという...


#http://www.kyoritsu-pub.co.jp/90thmath/math03.html

によると『リー群のユニタリ表現論』という本が刊行予定らしく、"非コンパクトリー群の例として一般化ローレンツ群SO(n,1)と擬ユニタリ群SU(n,1)をとり,それらの(非ユニタリ)主系列表現の構造を,上のGZ基底を拡張した基底を作って詳細に調べることにより"とあり、欲しい情報を書いてくれてそうな感じであるが、いつ出版されるか不明なので、自分で計算する必要がありそう。ぐぬぬ


以下は、計算の一部をRisa/Asirで確かめた時のコード。文章と式が合ってなかったら、文章の方が間違ってると思う

/* Risa/Asir code */

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


V=[z1,z2,z3,z4,d1,d2,d3,d4];
Z1=dp_ptod(z1,V);
Z2=dp_ptod(z2,V);
Z3=dp_ptod(z3,V);
Z4=dp_ptod(z4,V);
D1=dp_ptod(d1,V);
D2=dp_ptod(d2,V);
D3=dp_ptod(d3,V);
D4=dp_ptod(d4,V);

P0 = dp_ptod(z2*z3+z1*z4 , V);
P1 = dp_ptod(@i*(z1*z3-z2*z4) , V);
P2 = dp_ptod(z1*z3 + z2*z4 , V);
P3 = dp_ptod(z2*z3 - z1*z4 , V);
M1 = -(dp_weyl_mul(Z1,D2) - dp_weyl_mul(Z2,D1) - dp_weyl_mul(Z3,D4) + dp_weyl_mul(Z4,D3))/2;
M2 = @i*(dp_weyl_mul(Z1,D2) - dp_weyl_mul(Z3,D4) - dp_weyl_mul(Z4,D3) + dp_weyl_mul(Z2,D1))/2;
M3 = -@i*(dp_weyl_mul(Z1,D1) - dp_weyl_mul(Z2,D2) + dp_weyl_mul(Z3,D3) - dp_weyl_mul(Z4,D4))/2;

L1 = -@i*(dp_weyl_mul(Z1,D2) + dp_weyl_mul(Z3,D4) - dp_weyl_mul(Z2,D1) - dp_weyl_mul(Z4,D3))/2;
L2 = (-dp_weyl_mul(Z4,D3) - dp_weyl_mul(Z1,D2) - dp_weyl_mul(Z2,D1) - dp_weyl_mul(Z3,D4))/2;
L3 = (dp_weyl_mul(Z1,D1) - dp_weyl_mul(Z2,D2) - dp_weyl_mul(Z3,D3) + dp_weyl_mul(Z4,D4))/2;


RM1 = newmat(4,4,[[0,0,0,0],[0,0,0,0],[0,0,0,-1],[0,0,1,0]]);
RM2 = newmat(4,4,[[0,0,0,0],[0,0,0,1],[0,0,0,0],[0,-1,0,0]]);
RM3 = newmat(4,4,[[0,0,0,0],[0,0,-1,0],[0,1,0,0],[0,0,0,0]]);
RL1 = newmat(4,4,[[0,1,0,0],[1,0,0,0],[0,0,0,0],[0,0,0,0]]);
RL2 = newmat(4,4,[[0,0,1,0],[0,0,0,0],[1,0,0,0],[0,0,0,0]]);
RL3 = newmat(4,4,[[0,0,0,1],[0,0,0,0],[0,0,0,0],[1,0,0,0]]);

assert( dp_weyl_mul(P0,M1)==dp_weyl_mul(M1,P0) );
assert( dp_weyl_mul(P1,M1)==dp_weyl_mul(M1,P1) );
assert( dp_weyl_mul(L1,P2)==dp_weyl_mul(P2,L1) );
assert( dp_weyl_mul(L1,P3)==dp_weyl_mul(P3,L1) );
assert( dp_weyl_mul(P0,L1)-dp_weyl_mul(L1,P0)==P1 );
assert( dp_weyl_mul(P1,L1)-dp_weyl_mul(L1,P1)==P0 );
assert( dp_weyl_mul(M1,P2)-dp_weyl_mul(P2,M1)==P3 );
assert( dp_weyl_mul(M1,P3)-dp_weyl_mul(P3,M1)==-P2 );

assert( dp_weyl_mul(M1,M2)-dp_weyl_mul(M2,M1)==M3 );
assert( dp_weyl_mul(M2,M3)-dp_weyl_mul(M3,M2)==M1 );
assert( dp_weyl_mul(M3,M1)-dp_weyl_mul(M1,M3)==M2 );
assert( dp_weyl_mul(L1,L2)-dp_weyl_mul(L2,L1)==-M3 );
assert( dp_weyl_mul(L2,L3)-dp_weyl_mul(L3,L2)==-M1 );
assert( dp_weyl_mul(L3,L1)-dp_weyl_mul(L1,L3)==-M2 );
assert( dp_weyl_mul(M1,L2)-dp_weyl_mul(L2,M1)==L3 );
assert( dp_weyl_mul(L1,M2)-dp_weyl_mul(M2,L1)==L3 );
assert( dp_weyl_mul(M2,L3)-dp_weyl_mul(L3,M2)==L1 );
assert( dp_weyl_mul(L2,M3)-dp_weyl_mul(M3,L2)==L1 );
assert( dp_weyl_mul(L3,M1)-dp_weyl_mul(M1,L3)==L2 );
assert( dp_weyl_mul(M3,L1)-dp_weyl_mul(L1,M3)==L2 );
assert( RM1*RM2-RM2*RM1==RM3 );
assert( RM2*RM3-RM3*RM2==RM1 );
assert( RM3*RM1-RM1*RM3==RM2 );
assert( RL1*RL2-RL2*RL1==-RM3 );
assert( RL2*RL3-RL3*RL2==-RM1 );
assert( RL3*RL1-RL1*RL3==-RM2 );
assert( RM1*RL2-RL2*RM1==RL3 );
assert( RL1*RM2-RM2*RL1==RL3 );
assert( RM2*RL3-RL3*RM2==RL1 );
assert( RL2*RM3-RM3*RL2==RL1 );
assert( RM3*RL1-RL1*RM3==RL2 );
assert( RL3*RM1-RM1*RL3==RL2 );

/* ---------------------------------- */
P0 = z2*z3+z1*z4;
P1 = @i*(z1*z3-z2*z4);
P2 = z1*z3 + z2*z4;
P3 = z2*z3 - z1*z4;

/* lowest weight solution of SD-Maxwell equation */
Ex = -@i*(z2^2*z4^2 + z2^2*z3^2)/2;
Ey = (z2^2*z4^2 - z2^2*z3^2)/2;
Ez = z2^2*z3*z4;
Bx = -@i*Ex;
By = -@i*Ey;
Bz = -@i*Ez;

/* half-part of Maxwell equation */
assert( P1*Ex + P2*Ey + P3*Ez==0 );  /* div E = 0 */
assert( P2*Ez-P3*Ey==-P0*Bx );
assert( P3*Ex-P1*Ez==-P0*By );
assert( P1*Ey-P2*Ex==-P0*Bz );

/* lowest weight potential of SD-Maxwell equation */
A0 = -z2*z4/2;
A1 = @i*z2*z3/2;
A2 = z2*z3/2;
A3 = -z2*z4/2;

assert( -P1*A0 - P0*A1 == Ex );
assert( -P2*A0 - P0*A2 == Ey );
assert( -P3*A0 - P0*A3 == Ez );
assert( P2*A3-P3*A2 == Bx );
assert( P3*A1-P1*A3 == By );
assert( P1*A2-P2*A1 == Bz );