量子力学で登場する波動関数の記述に現れる演算子はエルミート行列なのですが、この挙動を調べるためにランダムエルミート行列というものが物理の方では研究されています。行列の性質を知るためには特性方程式やその根である固有値が重要になります。
実はゼータ関数とこのランダムエルミート行列には深い関係があります。具体的に言えば、ゼータ関数の非自明零点の分布とランダムエルミート行列の固有値の分布は一致するのです！！
とか言われると、実際に確認したくなりますよね。
久しぶりに眺めていたら「素数に憑かれた人たち」（ダービーシャー）にはこの辺のことが丁寧に描かれており、Maximaを使ってこの確認計算を再現することができたので、記事にまとめました。この本をお持ちの方は第18章を合わせて読んでいただけると良いと思います。
Maximaの主に数値計算や統計計算に関わる機能を色々と使っていますので、そちらも参考になるかもしれません。
今回は、ランダムエルミート行列の固有値の分布を確認します。適当に（適当とは言っても色々な条件はあるのですが）エルミート行列を作ってその固有値を求めてみます。「素数に憑かれた人たち」に合わせて大きさは4x4にします。これでランダムエルミート行列を実感してみてください。さらに固有値のリストを正規化してから固有値の間隔を求めてみます。
最後に全く同じ手続きを1000x1000の行列で実行してみます。999個の間隔が得られるのでそれをヒストグラムにして固有値の間隔の分布をグラフ化して確認します。
早速始める前に、ランダムエルミート行列の定義を述べておきます。
まずエルミート行列とは、対角要素は実数、それ以外の要素は複素数で、かつij要素とji要素が複素共役になっているような行列のことです。
ランダムエルミート行列とは、エルミート行列で各要素の実部、虚部が正規分布(中心0, 分散1)に従った乱数になっているような行列のことです。
では実際にMaximaで１つ作ってみます。まずは準備。実数の表示を小数点以下４桁に指定し、必要なパッケージを読み込みます。
(%i1) (fpprintprec:4,load(draw),load(descriptive),load(distrib),load(lapack))$
行列のサイズをDという変数で定義します。
(%i2) D:4;
$$ \tag{%o2} 4 $$
後で使う4x4のサイズの配列を用意します。
(%i3) array(RHM,D,D)$
次に4x4の複素行列の要素に必要な実部の乱数を4*(4-1)/2+4=10個生成します。この10個の乱数は中心0, 分散1.0の正規分布に従う必要があります。Distribパッケージのrandom_normal(m,s,N)という関数を使うと平均m, 分散sの正規分布に従う乱数をN個生成してリストにして返します。使う前にload(distrib)が必要です。
(%i4) srel:random_normal(0.0, 1.0, D*(D-1)/2+D);
$$ \tag{%o4} \left[ 0.5488 , 2.298 , 1.205 , 0.3392 , 1.411 , 1.181 , \\ -1.763 , 0.1562 , -1.015 , -0.2761 \right] $$
同じように虚数部の乱数を4*(4-1)/2=6個生成します（対角要素には虚数部はないのでその分少ないのです）。
(%i5) sim:random_normal(0.0, 1.0, D*(D-1)/2);
$$ \tag{%o5} \left[ -0.2082 , 1.063 , -0.199 , -0.159 , -0.6508 , -1.075 \right] $$
ランダムエルミート行列では対角要素に\( \sqrt{2} \)を掛ける必要があるので、そのための定数sqrt2を定義しておきます。
(%i6) sqrt2:ev(sqrt(2),numer);
$$ \tag{%o6} 1.414 $$
では先ほど定義した配列RHMに、要素を詰めていきます。i=jの時は対角要素なので虚数部はありません。
(%i7) for i:0 step 1 while i<D do
for j:i step 1 while j<D do
if i=j then RHM[i,j]:pop(srel)*sqrt2
else (relp:pop(srel),imp:pop(sim),RHM[i,j]:relp+%i*imp, RHM[j,i]:relp-%i*imp);
$$ \tag{%o7} \mathbf{done} $$
この配列RHMから行列を生成し、変数RHMmatに代入します。
(%i8) RHMmat:genmatrix(RHM,D-1,D-1,0,0);
$$ \tag{%o8} $$
$$ \notag \begin{pmatrix}0.7761&2.298-0.2082\,i&1.063\,i+1.205&0.3392-0.199\,i\\ 0.2082\,i+2.298&1.995&1.181-0.159\,i&-0.6508\,i-1.763\\ 1.205-1.063\,i&0.159\,i+1.181&0.2209&-1.075\,i-1.015\\ 0.199\,i+0.3392&0.6508\,i-1.763&1.075\,i-1.015&-0.3905 \end{pmatrix} $$
確かに、対角要素は全て実数ですし、それ以外の要素は全て複素数です。また例えば1,2要素 \( 2.298-0.2082\,i \) と2,1要素\(0.2082\,i+2.298 \)は複素共役になっていますし、他の要素も同様です。これが4x4のランダムエルミート行列です。じっくりと見てみてください。紹介した本の中でも4x4のランダムエルミート行列の例が挙げられていますが、当然、数値は全く異なります。ランダムなので、、、。
次に行列RHMmatの固有値を求めてみます。Maximaにはeigenvalues()という関数があり、結果を数式として求める場合にはそちらを使うのがオススメです。しかし今回は数値計算ですからlapackパッケージの中の関数を使います。この中にも固有値を求める関数は複数あるのですが、今回の目的にぴったりのものは、zheev()です。lapackのマニュアルページの最後の方に乗っています。
(%i9) eiglist:first(zheev(RHMmat));
$$ \tag{%o9} \left[ -2.53 , -1.182 , 1.226 , 5.088 \right] $$
4x4の行列には4つの固有値があります。そしてエルミート行列の固有値は全て実数になります。
ここからこの固有値のリストを正規化します。以下の過程は紹介した本のp343に書かれているのと全く同じです。
(%i10) eig2list:map(lambda([x],x-eiglist[1]),eiglist);
$$ \tag{%o10} \left[ 0.0 , 1.348 , 3.757 , 7.618 \right] $$
(%i11) eig3list:map(lambda([x],x/(eig2list[D]/D)),eig2list);
$$ \tag{%o11} \left[ 0.0 , 0.7077 , 1.972 , 4.0 \right] $$
これでうまく正規化ができたので、隣り合う固有値間の差を計算してみます。
(%i12) eig3difflist:[ ]$
(%i13) for i:1 step 1 while i<D do
eig3difflist:endcons(eig3list[i+1]-eig3list[i],eig3difflist);
$$ \tag{%o13} \mathbf{done} $$
(%i14) eig3difflist;
$$ \tag{%o14} \left[ 0.7077 , 1.265 , 2.028 \right] $$
この数字がどう分布しているのかを見るわけですが、3個ではもちろん何も言えません。そこで上記の過程をD:1000 として実行してみます。ただし途中経過は表示しないようにします。
(%i15) D:1000;
$$ \tag{%o15} 1000 $$
(%i16) array(RHM,D,D)$
(%i17) srel:random_normal(0.0, 1.0, D*(D-1)/2+D)$
(%i18) sim:random_normal(0.0, 1.0, D*(D-1)/2)$
(%i19) sqrt2:ev(sqrt(2),numer)$
(%i20) for i:0 step 1 while i<D do
for j:i step 1 while j<D do
if i=j then RHM[i,j]:pop(srel)*sqrt2
else (relp:pop(srel),imp:pop(sim),RHM[i,j]:relp+%i*imp, RHM[j,i]:relp-%i*imp)$
(%i21) RHMmat:genmatrix(RHM,D-1,D-1,0,0)$
(%i22) eiglist:first(zheev(RHMmat))$
(%i23) eig2list:map(lambda([x],x-eiglist[1]),eiglist)$
(%i24) eig3list:map(lambda([x],x/(eig2list[D]/D)),eig2list)$
(%i25) eig3difflist:[ ]$
(%i26) for i:1 step 1 while i<D do
eig3difflist:endcons(eig3list[i+1]-eig3list[i],eig3difflist)$
さあ、準備完了です。ちなみに固有値の間隔の分布は理論的に求められています。これも合わせてグラフに描いてみましょう。
(%i27) p2(s):=32/%pi^2*s^2*exp(-4/%pi*s^2)$
(%i28) p2(s);
$$ \tag{%o28} \frac{32\,s^2\,e^ {- \frac{4\,s^2}{\pi} }}{\pi^2} $$
(%i29) draw2d(grid=true, title="Spacing distribution of eigenvalues of a Hermite random matrix of size 1000x1000", histogram_description(eig3difflist, nclasses=50, fill_color=blue,frequency=density), explicit(p2(s),s,0,3.0), xrange=[0.0,3.0])$
かなり雰囲気はよく出ています。皆さんがお手元の環境で実行するとヒストグラムの形は多少変わると思います。元がランダムなので、、、。
このグラフから明らかに隣り合う固有値間の差が0に近い、ということは滅多にありません。本の中ではこのことを指して、「固有値の間には斥力が働いている」という表現が使われていました。
