nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

MKL11.3：逆行列計算ルーチンの罠

プログラミング

自分の今の環境では、MKL11.3を使ってfortran90でプログラムを弄っている。違う環境では、ここで言うことは当てはまらないかも知れない。

逆行列を計算したい。私の願いはそれだけであった。

一般の行列に対して対応可能なsubroutineは、GETRF (LU分解)とGETRI(それを踏まえて逆行列計算)である。頭文字のGはGeneralのGであろう。両方セットで使う。引数は、変換したい行列a(*)、変換するときの情報が入っているipiv (integer)、変換が上手く行ったか教えてくれるinfo (integer)、の三つのみである。aは実数でも倍精度実数でも複素数でも倍精度複素数でも、GETRFやGETRIは内部でそれぞれの型に対応した関数に渡して処理してくれる。
実際、GETRFやGETRIはLAPACK95のモジュールを使用しているのだが、これはInterfaceの集まりであり、計算そのものを実行しているのはF95_PRECISIONの中のCGETRF（単精度複素数用）、ZGETRF（倍精度複素数用）といった、型がキッチリ指定されたやつである。使うときには、LAPACK95のみをuseすると宣言してGETRF, GETRIをただ使えば十分である。

罠からは外れるが、使うときにはlapack95.modがあれば良く、f95_precision.modは無くて大丈夫。
ややこしい設定ファイルは、Windowsだと/Qmklとか/QMkl:sequentialとか、Linuxだと-mklで通ると書いてあるが、僕の環境では通らなかったので、自力でファイルを持って来たりした。使ったファイルだけ列挙すると、
libiomp5md.lib
mkl_core.lib
mkl_intel_ilp64.lib
mkl_intel_thread.lib
mkl_lapack95_ilp64.lib
である。libiomp5md.libは mkl/lib/intel64_win/ の中ではなく、compiler/lib/intel64_win/ に入っているため、探すのに苦労した気がする。
後、ilp64のファイルを使う場合には、コンパイルオプションに /4I8 を追加する必要がある。
本来は、Makefileをキチっと作って使い回すものなんでしょうけれども、Makefileはまだ勉強中であります。

んで、何が罠だったかというと、自分がやりたいのはエルミート行列の逆行列計算なので、エルミート行列用のHETRFとHETRIが使いたかった。頭文字のHはもちろんHermitianのHであろう。これまでは適当にGETRF,GETRIを使っていて、せっかくaがエルミートだから速く計算できる方が良かろうということであるのだが、引数の並びに問題があった。

call GETRF( a, [ipiv], [info] )
call GETRI( a, [ipiv], [info] )

に対して、

call HETRF( a, [uplo], [ipiv], [info] )
call HETRI( a, ipiv, [uplo], [info] )

である。uploは、Bunch-Kaufman diagonal pivoting methodにおいて上三角行列と下三角行列どちらを使うかというもので、'U'か'L'が入る。[...]はOPTIONAL引数で、省略可能である。
この引数の順番の違いで、コンパイルが全然通らなくてだいぶ困った。多分エラーメッセージが良くない。
「この汎用サブルーチン呼び出しと一致する特殊サブルーチンがありません。」
って言われたら、真っ先に設定ファイルが足りてないことを疑うのではないのだろうか？設定ファイルの多さと、GETRF,GETRIの流れから引数を疑うまでになかなか至らなかった。。。
もしMKLで詰まったら、大元のlapack.f90を見てみて下さい（というか、本来はそうするのが当たり前なんだろうけれども。。。）。

逐次代入法を用いた超シンプルな例と収束因子

数学

x = 1-x

を考える。もちろん答えは x = 0.5である。
これを敢えて逐次代入法で解いてみよう。この方法は固体物理や第一原理計算で用いるSCF法と同じノリである。

最初に初期値を設定する。答えに近ければ近い程良いが、とりあえず x_0 = 0.7とおく。これを右辺に代入したものが左辺を与えるとして、

x_1 = 1 - x_0 = 0.3
x_2 = 1 - x_1 = 0.7

と次々に項が求まる。が、よく見ると全然0.5に近づいておらず、このままだとずっと振動して xが求まらない。

そこで収束因子 \alphaを導入する。 \alphaは新しい xをどれくらい古い xに混ぜるかというもので、例えば x_1は次のようになる。

x'_0 = 1 - x_0
x_1 = \alpha x'_0 + ( 1 - \alpha) x_0

\alphaが小さい程、古い xの情報が残り、 xの変化が緩やかになる。
具体的に値を入れてその振舞を見てみると、

\alpha = 0.2
x_0 = 0.7
x'_0 = 0.3
x_1 = 0.2 \cdot 0.3 + 0.8 \cdot 0.7 = 0.06 + 0.56 = 0.62
x_2 = 0.2 \cdot 0.38 + 0.8 \cdot 0.62 = 0.072 + 0.496 = 0.568

というように振動が抑えられ、0.5に近づいていくのが分かる。
では本当にこれを続けて0.5になるかどうかを解析的に解いてみよう。

x_1 = \alpha ( 1 - x_0) + ( 1 - \alpha) x_0 = x_0 (1 - 2\alpha) + \alpha
x_2 = x_1 (1 - 2\alpha) + \alpha = x_0 ( 1 - 2\alpha)^2 + \alpha (1 - 2\alpha) + \alpha
\vdots
x_n = x_0 ( 1 - 2\alpha)^n + \alpha ( 1 + (1 - 2\alpha) + (1 - 2\alpha)^2 + \cdots + (1 - 2\alpha)^{n-1})

よって、これを無限に繰り返すと、 |1 - 2\alpha| < 1 であれば、

x_{\infty} = \alpha \frac{1}{1-(1-2\alpha)} = \frac{1}{2}

となり、初期値にも収束因子にも依存せず答えがちゃんと求まることが示せた。
|1 - 2\alpha| < 1 、つまり 0 < \alpha < 1 という条件は、元々が新しい xと古い xの割合ということを思い出せば全く自然な範囲である。
ちなみに、最初に x_0=0.5を代入すると一発で求まるし、 \alpha=0.5でも一発で求まるので、問題によって適切な初期値と収束因子を選ぶことが重要であることがわかる。しかし、一般に収束因子に関しては物理的なものではないので類推は困難であり、経験的に適切だと思われる値を入れることが多いと個人的には思う。
最適な収束因子を求める問題とか数値計算の方でありそうな気がするが、不用意に飛び込むと全ての数値計算法に満足出来なくなりそうで怖いのでやめておこう。。。

初等的な縮退の話

量子力学

サクライの章末問題に載っててフムフムとなった縮退の話。

演算子 A_1, A_2が互いに交換せず( [ A_1, A_2 ] \neq 0)、かつそれぞれがハミルトニアンと同時固有状態を作るとする [ A_1, H ] = 0, [ A_2, H ] = 0 。
このときに、一般に縮退が存在することが証明出来る。

これを真っ向から挑む場合に、必要十分であることを言うのはなかなかしんどい気がする。ホワンとした証明しか思いつかん。
これは背理法を使うと結構スッキリ話を持って行けて気持ちいい。

背理法として、縮退が無いとまず仮定する。あるエネルギー固有状態 |E_i>を A_1の固有値 a^{(1)}_iでラベル付けるか、 A_2の固有値 a^{(2)}_iでラベル付けるか好きな方を選択出来るが、 [ A_1, A_2 ] \neq 0の前提条件により a^{(1)}_iと a^{(2)}_iを一つの状態に同時に指定することは出来ない。しかし、今の仮定では縮退が無いため、

|a^{(1)}_i, E_i > = |a^{(2)}_i, E_i >

である。これは次の式で見るように、 [ A_1, A_2 ] \neq 0に抵触している。

( AB - BA ) |a^{(1)}_i, E_i > = (AB-Ba^{(1)}_i) |a^{(1)}_i, E_i > \\ = (AB-Ba^{(1)}_i)|a^{(2)}_i, E_i > = (A a^{(2)}_i - a^{(2)}_i a^{(1)}_i) |a^{(2)}_i, E_i > \\ = (A a^{(2)}_i - a^{(2)}_i a^{(1)}_i) |a^{(1)}_i, E_i > = (a^{(1)}_i a^{(2)}_i - a^{(2)}_i a^{(1)}_i) |a^{(1)}_i, E_i > \\ = 0

よって、縮退が無いという仮定が誤りであり、実際には縮退していることが言える。

これの具体例で良さそうだなぁと思ったのが点群で、例えば C_{3v}を考えて、 A_1 = C_3、 A_2 = \sigma_vと対応付けられて、図を描くとわかるが、 C_3 \sigma_v = \sigma''_v、 \sigma_v C_3 = \sigma'_v（鏡映面の異なる鏡映操作）となり、交換しない。加えて、系が C_{3v}に属していればハミルトニアンは C_{3v}の対称操作に対して変わらない。よって、エネルギー固有状態は縮退している。
実際、 C_{3v}は二次元の既約表現 Eを持ち、二重に縮退する。もちろん、一次元の全対称既約表現 A_1（どの対称操作に対しても不変）しかないような部分空間を用意するならば、話は別で縮退は無い（この部分空間においては [ C_3, \sigma_v ] = 0 が成立）。
でもやっぱ点群とかは図が無いと説明が絶望的に伝わらんなぁ。。。でもサクライだと角運動量 L_zと L_xで説明してて、交換しないのがただ数式的な感じがして個人的には縮退のイメージがピンとこなかったんだよなぁ。。。悩ましい。。。

二次元ハミルトニアンとパウリ行列

量子力学

H = a ( |1><1| - |2><2| + |1><2| + |2><1| )の固有値を求める問題をパウリ行列使って解く方法が割と楽しい（サクライでは章末問題になっている）。
これを行列で書くと、

\displaystyle H = a \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} = a \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} + a \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \\ = a ( \sigma_z + \sigma_x ) = a \sqrt{2}( \sigma \cdot \hat{n} )

ただし、 \hat{n} = \frac{1}{\sqrt{2}}(1,0,1)の単位ベクトルである。
せっかくなので、パウリ行列をスピン演算子に直してあげれば、 S_i = \frac{\hbar}{2}\sigma_iより

\displaystyle H = \sqrt{2} a \frac{2}{\hbar} ( {\bf S} \cdot \hat{n} )

で、変形して何なんだという話であるが、 {\bf S} \cdot \hat{n}の固有値は \pm \frac{\hbar}{2}で物理的に明らか（ \hat{n}方向を向いたスピンの固有値を表している。なぜそうなるのは自明というよりかは自然（実験事実）からの要請に近い。）であるため、このハミルトニアンの固有値は E = \pm \sqrt{2}aであることが行列式を使わずにわかってしまう。

ちなみに固有ベクトルだが、違う章末問題で {\bf S} \cdot \hat{n}の固有値を回転行列を使わないで真面目に固有値方程式を解く趣旨のものがあり、そっちも教育的で楽しい。固有値はもうわかっているから展開係数を直接求めるのだが、そういうプロセスはあんまりないなぁとボンヤリ思ったりもした。いずれそれもまとめたい。

ハミルトニアンの行列表示

量子力学

J.J. Sakuraiのゼミをしていたとき、ハミルトニアンの行列表示について議論になった。

\displaystyle H \doteq \begin{pmatrix} H_{11} & H_{12} & \cdots & \\ H_{21} & H_{22} & & \\ \vdots & & \ddots & \end{pmatrix}

Diracのbraketを用いた表示で書けば、

\displaystyle H \doteq \begin{pmatrix} <1|H|1> & <1|H|2> & \cdots & \\ <2|H|1> & <2|H|2> & & \\ \vdots & & \ddots & \end{pmatrix}

となる。
Sakuraiだと、イコールではないことを強調して一応 \doteqを用いてはいるものの、結構ポッと出な感もあるし、何よりも <1|H|1>の中の Hは外の Hとどう違うの？とか、もう一回中に代入するってこと？みたいな誤解が起きることが分かった。確かに H \psi(x) = E \psi(x) に直接行列表現ぶち込もうとすると訳分からんことになる。

ポイントは一本の二階微分方程式を解くモードか、微分を含まない連立方程式で解くモードかということだと思っている。

シュレーディンガー方程式から出発して、全然分からん固有状態 \phiを素性が良くわかっている固有状態 \phi_iで展開する（例えば何にも相互作用の無いときの固有状態である平面波とか）。

\displaystyle H |\psi> = E |\psi> \\ H ( C_1 |\phi_1> + C_2 |\phi_2> + \cdots ) = E ( C_1 |\phi_1> + C_2 |\phi_2> + \cdots )

ここで、例えば <\psi_1|を左から掛けることで、 H|\psi>にどれだけ \phi_1の情報が含まれているかを見る様なことをすると、

\displaystyle <\phi_1| H |\psi> = <\phi_1| E |\psi> \\ C_1 <\phi_1|H|\phi_1> + C_2 <\phi_1|H|\phi_2> + \cdots = E C_1

これを他の状態 <\phi_i|でもやってあげれば、展開係数 C_iに関する連立方程式が出来る。

\displaystyle C_1 <\phi_1|H|\phi_1> + C_2 <\phi_1|H|\phi_2> + \cdots = E C_1 \\ C_1 <\phi_2|H|\phi_1> + C_2 <\phi_2|H|\phi_2> + \cdots = E C_2 \\ \vdots

ここまで、 Hは微分方程式（シュレーディンガー方程式）のハミルトニアンのままであり、別に変なことはしていない。この連立方程式を行列でまとめて書くと、

\displaystyle \begin{pmatrix} <\phi_1|H|\phi_1> & <\phi_1|H|\phi_2> & \cdots & \\ <\phi_2|H|\phi_1> & <\phi_2|H|\phi_2> & & \\ \vdots & & \ddots & \end{pmatrix} \begin{pmatrix} C_1 \\ C_2 \\ \vdots \end{pmatrix} = \begin{pmatrix} E & & & O \\ & E & & \\ O & & \ddots & \end{pmatrix} \begin{pmatrix} C_1 \\ C_2 \\ \vdots \end{pmatrix}

という具合になる。
なので、 E,\psiを求める問題が、 E, C_iを求める問題にすり替わり、特に Eを求めたければ C_iは必要なく、

\displaystyle \begin{vmatrix} <\phi_1|H|\phi_1> - E & <\phi_1|H|\phi_2> & \cdots & \\ <\phi_2|H|\phi_1> & <\phi_2|H|\phi_2> - E & & \\ \vdots & & \ddots & \end{vmatrix} = 0

を解きさえすれば良い（これがゼロでない場合には、上の行列に逆行列が存在して、 C_i = 0（つまり \psi = 0）という全然興味のない答えしか得られないため、ゼロ出ないことを条件付ける訳である）。
そうすると、 Hの微分方程式における具体的な表式は、「行列要素の値が求まっていれば」どうでも良くなり、行列形式で問題を扱って行くときにいちいち違う記号を用いるのは煩雑になる（人によって違う記号を使いだすと多分ヤバい）ので、同じ Hを使っているのだと理解している。
最初の話に戻ると、行列の中の Hは実際のシュレーディンガー方程式の Hで、 \doteqの左側にいる Hはただ Hという記号を使っている行列ということになる。

fortranとpython（読み込み）

プログラミング

fortranユーザーからpythonに移って来ると、ただテキスト編集したいのにも関わらず、fortranの癖に慣れてしまったが故にpythonがなかなか使えるようにならんかったのでちょっと整理。

「ファイルを開いて中身を読んで閉じる」までをfortranとpythonで比較。

implicit none
integer, parameter :: INPUT = 10
real(kind(0d0))    :: InputData( 50, 50 ), ...
!----------
InputData
open( INPUT, 'test.inp', status='old' )
do i = 1, 50
  read(INPUT, *) InputData( i, : )
enddo
close( INPUT )

test.inpが 50 x 50 以上のデータを持っていると仮定して、InputDataの配列サイズを指定。
'test.inp'に対応する番号として10(=INPUT)を割り当ててあげてファイルを開き、i行目の各列のデータをreadでInputDataに流し込む。
それを50行目まで繰り返してファイルを閉じる。

f = open( 'test.inp', 'r' )
lines = f.read().split('\n')
f.close()
InputData = []
for i, line in enumrate( lines ):
    temps = line.split(' ')
    while '' in temps: temps.remove('')
    InputData += temps

読み込み形式('r')で'test.inp'を開き、中身をlinesにいきなりリストとして突っ込みファイルは閉じる。.split('\n')によって改行ごとに（つまり各行が）違う要素としてlinesの中に収納されている。lines[0]は一行目全部、lines[1]は二行目全部、といった具合。
それを値として全部ぶつ切りに区別したデータとするため、linesの中の要素line（つまりある一つの行）を空白で区切って（.split(' ')）、空白が連続している場合には空の要素''が生まれてしまうためwhile文を使って除いている。
綺麗に区切ったデータをInputDataに突っ込む作業を全ての行に対して行う。


おそらくpythonでfortranと同じように読ませようとするとこんな感じになるのではないかと思う。
fortranは「行列」を最初から想定しているような節を感じる。型が決まっているなら、その型にハマる様に最初に設定してしまえばあんまり何も考えなくて済む（その代わり型から外れると酷いが）。
一方でpythonは、一本の長いリストに収納してからグニグニ変形させる。かなり自由が利く分、自分がやりたいことを実践するために出来る方法も多く、結果どうしたら良いかわからなくなることも多かった（ファイルを読み込むだけなのにいくつか方法があるというのも新参者には割とストレスになった。。。）。
とりあえず今はこんな感じで落ち着いているが、もっとスマートな方法が多分あると思う。先は蛇のように長い。。。

奇数次元の虚数について

数学

虚数 iは二次元の正方行列を用いて、

i_2 = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}

と書ける。
では、三次元の正方行列で二乗したら単位行列にマイナスの掛かったものが得られるかをボンヤリ考えたところ、

i_3 = \begin{pmatrix} 0 & 0 & -1 \\ 0 & i & 0 \\ 1 & 0 & 0 \end{pmatrix}

しか思いつかんかった。この iを二次元 i_2と思えば、super matrixと解釈して各行列要素が二次元正方行列に拡張出来るから、

\begin{pmatrix} 0 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \end{pmatrix}

ちゃんと二乗するとマイナス単位行列になる。もちろん、基底を組み替えたものに対応した、

\begin{pmatrix} 0 & 0 & 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}

もちゃんとなるし、どちらかというとこっちの方が帰納的でしっくり来るだろう。

でも、最初の三次元の iを自分自身の i_3で置き換えると無限に行列の次元が増えていく。
世界の自己相似性に奇数次元であることが絡んでいると面白そうだなぁと思った次第である。