正規直交化をした
ちょっと必要だったんで行列の正規直交化しました．
すっごいつかいふるされた話題だと思うんですけど，でも思ったより参考になるサイトとか見つけられなかったので，メモしときます．
ちなみに数学っぽい話のようでプログラムの話です，しかもFortranです．さらに言うとあんまりレベル高い話じゃないです．
Gram-Schmidtの直交化
有名な直交化の方法の一つ．数学的な話はGoogle先生に聞けばわかると思います．google:Gram-Schmidt 直交化
直交化でググったら最初に出てきたのでまずはこれでやってみました．最初に言っておくとうまくいってません．うまくいきました！
プログラム的には
- 行列を格納した配列の列の数(N)だけループ回す(index I ≤ N)
- そのなかで未満の数だけループ回す(index J < I )
- I番目とJ番目の行の内積をとる
- I番目の行から(↑の内積)×(J番目の行)を引く
- 各行を正規化する．つまりノルムで割る．
こんな感じです．わかりづらかったら他のサイトを探してください．
さらに上の繰り返しを数回繰り返すと誤差が小さくなっていって良いってどこかで見ました．
というわけでそれをfortranで書いたコードが下のやつです．
program GramSchmidtOrtho use matrix_io implicit none integer :: Ortho_Loop, I, J, K, N double precision :: Ortho(6,6), Test(6,6), InnerP data Ortho/ 5.d0, 3.d0, 2.d0, 1.d0, 0.d0, 9.d0, & & 3.d0, 1.d0, 1.d0, 2.d0, 5.d0, 10.d0, & & 2.d0, 1.d0, 9.d0, 1.d0, 3.d0, 1.d0, & & 1.d0, 2.d0, 1.d0, 4.d0, 1.d0, 5.d0, & & 0.d0, 5.d0, 3.d0, 1.d0, 0.d0, 2.d0, & & 9.d0, 10.d0, 1.d0, 5.d0, 2.d0, 4.d0/ External dgemm double precision, External :: ddot, dnrm2 N=6 print *, "Gram-Shmidt Orthogonalization." do Ortho_Loop=1, 2 do I=1, N if (I .ne. 1) then do J=1, I-1 InnerP = ddot(N, Ortho(:,I), 1, Ortho(:,J), 1) do K=1, N Ortho(K,I) = Ortho(K,I) - (InnerP * Ortho(K,J)) end do end do end if Ortho(:,I) = Ortho(:,I) / dnrm2(N,Ortho(:,I),1) end do end do call dgemm('T', 'N', N, N, N, 1.d0, Ortho, N, Ortho, N, 0.d0, Test, N) call matrix_print(N, N, Test, "I?") end program GramSchmidtOrtho
行列は適当な実対称行列です．
ddot, dnrm2, dgemmはblasの関数とサブルーチンです．
- ddot
- 内積を求める関数
- dnrm2
- ノルムを求める関数
- dgemm
- 行列の積を求めるサブルーチン．今回は trans(Ortho)*Ortho = I を計算するために使用
あと，matrix_printって言うのは行列を表示させるために自作したサブルーチンなので気にしないでください．ただ行列表示するだけです．matrix_ioっていうのはそれが入ってるモジュールです．
合ってるのか少々不安ですけど，とりあえず実行するとこうなります．
Gram-Shmidt Orthogonalization.
<< I? >>
1 2 3 4 5
1 1.00000 0.00000 0.00000 0.00000 0.00000
2 0.00000 1.00000 0.00000 0.00000 0.00000
3 0.00000 0.00000 1.00000 0.00000 0.00000
4 0.00000 0.00000 0.00000 1.00000 0.00000
5 0.00000 0.00000 0.00000 0.00000 1.00000
6 0.00000 0.00000 0.00000 0.00000 0.00000
6
1 0.00000
2 0.00000
3 0.00000
4 0.00000
5 0.00000
6 1.00000
おおー！すげー！
と思ったんですが，肝心の自分の仕事で使ってみたらなぜか非対角要素にゴミが残ります．
1 2 3 4 5
1 1.00000 0.00000 0.00000 0.00000 0.00000
2 0.00000 1.00000 0.00000 0.00000 0.00000
3 0.00000 0.00000 1.00000 0.00000 0.00000
4 0.00000 0.00000 0.00000 1.00000 0.00000
5 0.00000 0.00000 0.00000 0.00000 1.00000
6 0.00000 0.00000 0.00000 0.00000 0.00000
7 0.00000 0.00000 -0.01951 0.00000 0.00000
8 0.00000 0.00000 0.00000 0.00000 0.00000
9 0.00000 0.00000 0.00000 0.00000 0.02435
10 0.00000 0.00000 0.01045 0.00000 0.00000
11 0.00000 0.00000 0.00000 0.00000 0.00000
12 0.00000 0.00000 0.00000 0.00000 0.04901
13 0.00000 0.00000 0.00000 0.00000 0.00000
14 0.00000 0.00000 0.00000 0.00000 0.00000
15 0.05799 0.00000 0.00000 0.00000 0.00000
16 0.00000 0.00000 0.39492 0.00000 0.00000
17 0.00000 0.00000 0.00000 0.00000 0.02842
18 0.00000 0.00000 0.00000 0.00000 0.00000
19 0.00000 0.00000 0.00000 0.10721 0.00000
20 0.00000 0.00000 0.42938 0.00000 0.00000
21 0.00000 0.00000 0.00000 0.00000 0.00000
22 0.00000 0.00000 -0.17855 0.00000 0.00000
ちなみにこの場合は102*102の行列を使っていてそこからの抜粋です．
さっきの誤差が少なくなるっていう繰り返し回数もアホみたいに増やしたりしてみたんですが，ほとんど効果ありませんでした．
うーん，よくわからん…．何をミスったのかわかりません．とりあえずGram-Schmidtはあきらめました．
原因がわかりました．自分のお仕事のプログラムではいろいろ試行錯誤してる途中に型宣言を直し忘れてdouble precisionをひとつintegerのままにしてました…．凡ミスにもほどがあるぜ！気付けよ！なんてこったい／(^o^)＼
LAPACKガシガシ使うHouseholder分解からの直交化
ここら辺で結構時間が経過して自分でやる気を失い，LAPACKで何とかやる方法はないもんかと探し始めます．なんか負けた気がする…，と思う心を押さえつけて…．
するとdgeqrfというのが見つかりました．どうやらこれでQR分解ができる模様．
そしてさらにdorgqrというのを使うと，そのQ行列から直交行列が求められるみたいです．
というわけでやってみました．内容はさっきと一緒．
program HouseholderOrtho use matrix_io implicit none integer :: Ortho_Loop, I, J, K, N integer :: lwork, info double precision :: Ortho(6,6), Test(6,6), TAU(6), lw double precision, allocatable :: work(:) data Ortho/ 5.d0, 3.d0, 2.d0, 1.d0, 0.d0, 9.d0, & & 3.d0, 1.d0, 1.d0, 2.d0, 5.d0, 10.d0, & & 2.d0, 1.d0, 9.d0, 1.d0, 3.d0, 1.d0, & & 1.d0, 2.d0, 1.d0, 4.d0, 1.d0, 5.d0, & & 0.d0, 5.d0, 3.d0, 1.d0, 0.d0, 2.d0, & & 9.d0, 10.d0, 1.d0, 5.d0, 2.d0, 4.d0/ External dgemm, dgeqrf, dorgqr double precision, External :: ddot, dnrm2 N=6 print *, "Orthogonalization using LAPACK." !!! Householder QR decomposition call dgeqrf(N, N, Ortho, N, TAU, lw, -1, info) lwork = int(lw) allocate(work(lwork)) call dgeqrf(N, N, Ortho, N, TAU, work, lwork, info) deallocate(work) !!! Get orthogonalized matrix from QR matrix call dorgqr(N, N, N, Ortho, N, TAU, lw, -1, info) lwork = int(lw) allocate(work(lwork)) call dorgqr(N, N, N, Ortho, N, TAU, work, lwork, info) call dgemm('T', 'N', N, N, N, 1.d0, Ortho, N, Ortho, N, 0.d0, Test, N) call matrix_print(N, N, Test, "I?") end program HouseholderOrtho
実行結果はもちろん…
Orthogonalization using LAPACK. << I? >> 1 2 3 4 5 1 1.00000 0.00000 0.00000 0.00000 0.00000 2 0.00000 1.00000 0.00000 0.00000 0.00000 3 0.00000 0.00000 1.00000 0.00000 0.00000 4 0.00000 0.00000 0.00000 1.00000 0.00000 5 0.00000 0.00000 0.00000 0.00000 1.00000 6 0.00000 0.00000 0.00000 0.00000 0.00000 6 1 0.00000 2 0.00000 3 0.00000 4 0.00000 5 0.00000 6 1.00000
OK!
肝心の自分のお仕事では…
1 2 3 4 5
1 1.00000 0.00000 0.00000 0.00000 0.00000
2 0.00000 1.00000 0.00000 0.00000 0.00000
3 0.00000 0.00000 1.00000 0.00000 0.00000
4 0.00000 0.00000 0.00000 1.00000 0.00000
5 0.00000 0.00000 0.00000 0.00000 1.00000
6 0.00000 0.00000 0.00000 0.00000 0.00000
7 0.00000 0.00000 0.00000 0.00000 0.00000
8 0.00000 0.00000 0.00000 0.00000 0.00000
9 0.00000 0.00000 0.00000 0.00000 0.00000
10 0.00000 0.00000 0.00000 0.00000 0.00000
11 0.00000 0.00000 0.00000 0.00000 0.00000
12 0.00000 0.00000 0.00000 0.00000 0.00000
13 0.00000 0.00000 0.00000 0.00000 0.00000
14 0.00000 0.00000 0.00000 0.00000 0.00000
15 0.00000 0.00000 0.00000 0.00000 0.00000
16 0.00000 0.00000 0.00000 0.00000 0.00000
17 0.00000 0.00000 0.00000 0.00000 0.00000
18 0.00000 0.00000 0.00000 0.00000 0.00000
19 0.00000 0.00000 0.00000 0.00000 0.00000
20 0.00000 0.00000 0.00000 0.00000 0.00000
21 0.00000 0.00000 0.00000 0.00000 0.00000
22 0.00000 0.00000 0.00000 0.00000 0.00000
わーい(^o^)LAPACKすごいなー！
でもやっぱり何か負けた気がする．
