scilabによるSU(9)の実装

昨日のブログに書いた、9次の特殊ユニタリSU(9)の実装を、scilabで行ってみます。
scilabは、信号処理に特化したオープンソースのプログラミング言語なので
複素行列の扱いが得意です。


まず80個の自由度を持つので、0から1までの間の一様（？）乱数を80個生成しましょう。
1行80列のベクトルsとして生成します。

s=rand(1,80)

次に、このベクトルのノルム(長さ)を1にします。
s=s/norm(s)

できあがったら一応、長さが1になっていることを確認しましょう。

norm(s)


9つの対角成分の入力をします。
x1～x9の9つです。

x1=s(3)+s(8)/sqrt(3)+s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x2=-s(3)+s(8)/sqrt(3)+s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x3=-2*s(8)/sqrt(3)+s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x4=-3*s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x5=-4*s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x6=-5*s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x7=-6*s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x8=-7*s(63)/sqrt(28)+s(80)/sqrt(36)

x9=-8*s(80)/sqrt(36)

ここで、計算ミスがないように、検算をしてみましょう。
x1からs9までの和はゼロになるはずです。



次に、非対角成分のうちの上三角部分を定義しましょう。

a12=s(1)+%i*s(2)

a13=s(4)+%i*s(5)
a23=s(6)+%i*s(7)

a14=s(9)+%i*s(10)
a24=s(11)+%i*s(12)
a34=s(13)+%i*s(14)

a15=s(16)+%i*s(17)
a25=s(18)+%i*s(19)
a35=s(20)+%i*s(21)
a45=s(22)+%i*s(23)

a16=s(25)+%i*s(26)
a26=s(27)+%i*s(28)
a36=s(29)+%i*s(30)
a46=s(31)+%i*s(32)
a56=s(33)+%i*s(34)

a17=s(36)+%i*s(37)
a27=s(38)+%i*s(39)
a37=s(40)+%i*s(41)
a47=s(42)+%i*s(43)
a57=s(44)+%i*s(45)
a67=s(46)+%i*s(47)

a18=s(49)+%i*s(50)
a28=s(51)+%i*s(52)
a38=s(53)+%i*s(54)
a48=s(55)+%i*s(56)
a58=s(57)+%i*s(58)
a68=s(59)+%i*s(60)
a78=s(61)+%i*s(62)

a19=s(64)+%i*s(65)
a29=s(66)+%i*s(67)
a39=s(68)+%i*s(69)
a49=s(70)+%i*s(71)
a59=s(72)+%i*s(73)
a69=s(74)+%i*s(75)
a79=s(76)+%i*s(77)
a89=s(78)+%i*s(79)

要素を書きだしたら、これを行列に入れます。

いっぺんにやると凡ミスをしやすいので、1行ずつ入れることにします。
B1=[0,a12,a13,a14,a15,a16,a17,a18,a19]
B2=[0,0,a23,a24,a25,a26,a27,a28,a29]
B3=[0,0,0,a34,a35,a36,a37,a38,a39]
B4=[0,0,0,0,a45,a46,a47,a48,a49]
B5=[0,0,0,0,0,a56,a57,a58,a59]
B6=[0,0,0,0,0,0,a67,a68,a69]
B7=[0,0,0,0,0,0,0,a78,a79]
B8=[0,0,0,0,0,0,0,0,a89]
B9=[0,0,0,0,0,0,0,0,0]

A=[B1;B2;B3;B4;B5;B6;B7;B8;B9]

こうして、上三角部分は埋まりました。

エルミート行列なので、下三角は上三角の複素共役です。

A=A+A'

ここまでで、非対角成分すべてが埋まりました。

念のため、ここでエルミートであることを確認してみましょう。
A-A'=0

対角要素を入れます。
A(1,1)=x1
A(2,2)=x2
A(3,3)=x3

A(4,4)=x4
A(5,5)=x5
A(6,6)=x6
A(7,7)=x7
A(8,8)=x8
A(9,9)=x9


これでエルミート行列Aが完成したので
今一度A-A'=0が成立するのを確かめておきましょう。


ここで、いきなりexpm(%i*A)とやってもいいのですが、せっかくなので寄り道して

Aの固有値方程式と固有値を見てみましょう。

変数xを定義します。
x=poly(0,"x")

次に、固有値方程式Bを見てみましょう。

B=det(A-x*eye(9,9))
detは行列式、eye(n,m)はn行m列の単位行列を意味します。

 0.0000013 + 0.0000969x + 0.0011569x^2 - 0.0025147x^3 - 0.0637862x^4 - 0.1498150x^5 + 0.2886920x^6 + x^7 - x^9

しっかりと、7次の符号が9次の符号の逆相になっていて、8次の係数がゼロになっていますね。
係数がちゃんと実数だけでできていることも確認しておきましょう。

この9次方程式の解はすべて実数です。エルミート行列であることからそういえます。

Aの固有値がこの9次方程式の解なのですが

spec(A)を計算してみるとわかります。

  - 0.6196769  

  - 0.4362829  

  - 0.2728734  

  - 0.2080910  

  - 0.0761649  

  - 0.0165342  

    0.1421505  

    0.4476562  

    1.0398167 

エルミート行列をただの数に例えると実数
歪エルミート行列は純虚数に対応するため
expm(%i*A)は、複素平面上の単位円に相当します。

これがユニタリ＝ユニットな行列と呼ばれる所以です。

C=expm(%i*A)
を計算してみて(行列指数関数などの場合は、関数名の後ろに行列(matrix)を意味するmを付ける必要があります)

spec(C)を計算してみますと

    0.5063784 + 0.8623114i  

    0.9014641 + 0.4328538i  

    0.8140662 - 0.5807721i  

    0.9899136 + 0.1416723i  

    0.9063287 - 0.4225735i  

    0.9630005 - 0.2694997i  

    0.9784271 - 0.2065925i  

    0.9998633 - 0.0165335i  

    0.9971009 - 0.0760913i  
このようになり、固有値がすべて複素平面における単位円周上にあることが確認できます。


また、このユニタリ行列は「特殊」ユニタリ行列といって、
abs(det(C))とやるまでもなく
det(C)そのものが1になるという性質があります。

つまり個の場合は9つすべての固有値の積が、ちょうど実正数の1になるということです。



なお、今回乱数にすべて正の数を用いたのは、
複素共役を見やすくするためでした。
Aの上三角の虚部がすべて正の数で、下三角の虚部がすべて負の数になります。