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

1 行列の計算

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

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

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

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

1.1 逆行列の求め方

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

• 余因子行列を使う方法

• Gauss-Jordan の掃き出し法を使う方法 (線形代数の本 [1] の p47)

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

などがある。実際のプログラムやアルゴリズムが簡単であるので、ここでは 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 行列の基本変形

行列の (左) 基本変形は次のようなものである。

1. 2 つの行の入れ換え。

2. ある行に 0 でない数を掛ける。

3. ある行に別の行の定数倍を加える。

このうち、 2、 3 を使うと、単位行列への変換 (+逆行列の導出) が一応出来る。基本変形 1. は後出の pivoting で使う。

1.2.2 プログラムの指針 (A を E にする)

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

• \underline{ステップ 1.}

以下のステップ 2,3 を k = 1 から k = n まで順次実行する。

• \underline{ステップ 2.}

中間変数 wa := A_{k,k} を導入し、第 k 行を wa で割ったものに置き換える。

A_{k,j} := A_{k,j}/wa (j = 1, ⋯, n)  (1)
• \underline{ステップ 3.}

第 k 行以外の i = 1 ∼ n 行(i = 1, ⋯, k − 1, k + 1, ⋯, n) について以下の操作を行なう。中間変数 wb := A_{i,k} を導入し、第 i 行から第 k 行を wb 倍したものを引く。

A_{i,j} := A_{i,j} − A_{k,j} ∗ wb (j = 1, ⋯, n)  (2)

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

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

∘ 確認

この操作で、行列 A を単位行列に変形できることを確かめよ。 1.2.3 行列出力の方法 (do 形反復)

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

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

とした方がスマート。行列の入力にも同じ方法が使える。

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

A から E へ変形した操作を単位行列にもほどこしてやれば、 E を B = A^{−1} にすることができる。操作としては、始めに (A : E) からなる n 行 2n 列の行列を作り、これに対し掃き出し法の操作 (但し、列の操作を j = 1, ⋯, n から j = 1, ⋯, 2n へ拡張) を行なう。

実際に

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

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

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

これで A の逆行列が求まったが、正しい答えになってるかを確認する。実際に A と B の行列の積を計算して単位行列になるかを確かめればいい。

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) について計算すると行列の積が求まる。課題先ほど求めた B と 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 列目を処理する時に、

A_{k,k}, A_{k+1,k}, ⋯, A_{n−1,k}, A_{n,k}

の中から絶対値の 1 番大きな成分を探す。それを A_{l,k} とした時に、第 l 行と第 k 行とを入れ替えて計算する。第 k 行より下の要素のみ扱ったのは、既に半ば出来上がっている単位行列を壊さないようにするためである。

1.5 その他

1.5.1 計算量について

• 余因子行列法と比較 (O(N^{3})vs.O(N!))

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

1.5.2 メモリ節約型

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

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

参考文献

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

[2] 「数値計算の常識」伊理正夫・藤野和建著 (共立出版)

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