物理学特別講義V(物理情報処理論):No.8

1 行列の計算

様々ある行列の計算の中から、ここでは逆行列を求める事を行う。逆行列を求めると、

• 線形の連立方程式が解ける。

• AA^{−1} = E となってうれしい。

などの利点がある。 (E は単位行列を表す)

1.1 逆行列の求め方

逆行列を求める方法としては、

• 余因子行列を使う方法

• Gauss-Jordan の掃き出し法を使う方法 (数学の本の p47)

• 三角行列 (LU 行列) にわけて計算する方法

などがある。実際のプログラムやアルゴリズムが簡単であるので、ここでは Gauss-Jordan の掃き出し法を使って、逆行列を求める。

1.2 Gauss-Jordan の掃き出し法

A の逆行列を求めたい。まずは、 A と E を並べて (A : E) とします。基本操作 (線形代数の本 [1] p4) を行って、 A が単位行列になるように変形します。その変

形により E が、 B になったとすると、 B = A^{−1} となります。

なぜならば基本操作を繰り返す事は、基本操作の行列 (数学の本 p87) の積を左からかける事である。 A を E にした時に使った基本操作の行列の積を B として、

B(A : E) = (BA : B)
= (E : B)

となるので、 B = A^{−1} が得られる。 1.2.1 プログラムの指針 (A を E にする)

掃き出し法の指針としては、 1 列目から右に向かって順に掃き出していく (1 ∼ n 列まで)] 具体的には

\underline{ステップ 1.} k 行を A_{kk} で割ったものに置き換える。つまり j = 1, ⋯, n まで、行列の成分を次のよう

に書き換える。

A_{kj} := A_{kj}/A_{kk}  (1)

\underline{ステップ 2.} k 行以外の 1 ∼ n 行について以下の操作を行なう。

A_{ij} := A_{ij} − A_{kj}A_{ik}  (i = 1, ⋯, k − 1, k + 1, ⋯, n, (j = 1, ⋯, n))  (2)

ステップ 1、 2 により、 k 列の非対角成分 0、対角成分 1 になる。

( \begin{matrix}? & 0 & ? \\ ? & : & ? \\ ? & 0 & ? \\ ? & 1 & ? \\ ? & 0 & ? \\ ? & : & ? \\ ? & 0 & ? \end{matrix} )

\underline{ステップ 3.} ステップ 1、 2 を k = 1 から k = n まで順次実行する。

∘ 確認この操作で、行列 A を単位行列に変形できることを確かめよ。

1.2.2 行列出力の方法 (do 形反復)

行列を出力する時、 do loop 文で各成分を順に書き出しても良いが、

      do i=1,L
        write (*,*) (A(i,j), j=1,N)
      end do

とした方がスマート。

1.2.3 (E : B) を作るプログラム

A から E を作る時のプログラム上の操作を単位行列にもほどこしてやれば、 E を B にすることができる。操作としては、始めに (A : E) からなる n 行 2n 列の行列を作り、これに対し掃き出し法の操作 (但し、列の操作を j = 1, ⋯, n から j = 1, ⋯, 2n へ拡張) を行なう。実際に

( \begin{matrix}1 & 2 \\ 3 & 4 \end{matrix} )

などの逆行列を求めてみましょう。

1.3 行列の積を計算するプログラム

これで A の逆行列がもとまったのであるが、正しい答えになってるかを確かめておいた方が安心である。逆行列になってるかどうかをチェックするばいいのだから、 A と、 A^{−1} の行列の積を計算して実際に単位行列になるかを確かめればいい。

L ＊ N 行列のものと、 N ＊ M 行列のものの積を計算する。

C(L,M)=A(L,N)*B(N,M)

C の各成分について書いてやると、

C_{ij} = ∑_{k=1}^{N} A_{ik}B_{kj}

となる。各成分を (i = 1 ∼ L) と、 (j = 1 ∼ M) について計算してやれば行列の積が求まる。課題先ほど求めた A^{−1} と A との積を計算して単位行列になっていることを確かめましょう。

1.4 Pivoting

では、次に

( \begin{matrix}3 & 2 & 5 \\ 6 & 4 & 2 \\ 1 & 1 & 1 \end{matrix} )

の逆行列を求めてみましょう。失敗したのではないでしょうか？なぜだ？（答え： A_{kk} = 0 だから）この行列の行列式は 0 ではないので逆行列は存在する。逆行列を求めようと思ったら、 0

で割り算をしないようにすればいい。その方法が pivoting である。 pivoting には、

• full pivoting

• partial pivoting

がある。簡単な partial pivoting を使う。 1.4.1 Partial Pivoting の指針

掃き出し法のプログラムで k 列目の処理する時に、 k 列の中から一番大きい成分を探してその行と k 行とを入れ替えて計算する。

1.5 その他

1.5.1 計算量について • 余因子行列法と比較 (O(N^{3})vs.O(N!))

• LU 行列を求める方法よりは、 1.5 倍多い

1.5.2 メモリ節約型

(Fortran の本 [2] の p73) に Gauss-Jordan の掃き出し法が載っている。この方法だと、行列を 1 つだけ用意すればいいので、メモリの節約になる (但しプログラムはやや複雑になる)。

このプログラムでは、元の行列 A から逆行列を計算する時に、使ってしまったデータのところに計算していった逆行列のデータを入れている。

参考文献

[1] 「基礎線形代数」押川元重・坂口紘治 (培風館)

[2] 「Fortran90 入門」新井親夫 (森北出版株式会社)