逆行列を求めるプログラムの考察 (2) — 掃き出し法

[latexpage]
聖なる夜は 聖夜に薪を燃やし続けるだけ【薪が燃える27時間生放送】 を見て過ごしました。

さて、本題です。逆行列を求める手法について考察するシリーズです。逆行列を求める手法には

  • LU分解
  • 掃き出し法

があり、このうちLU分解を利用したプログラムは前回紹介しました。まだお読みでない方はぜひどうぞ。

逆行列を求めるプログラムの考察 (1) — LU分解

今回は掃き出し法を利用します。

原理

掃き出し法はガウスの消去法ともいい、多元連立1次方程式を求める方法として有名です。次の連立1次方程式を考えてみましょう。

\[ \begin{cases} \phantom{1}x_1+3x_2-2x_3+2x_4=9 \\ 2x_1-\phantom{1}x_2+5x_3-3x_4=3 \\ 3x_1+2x_2-6x_3+5x_4=9 \\ 5x_1+\phantom{1}x_2-7x_3+5x_4=6 \end{cases} \]

上から順に第1式、第2式、第3式、第4式と呼ぶことにします。

まず、前進消去を行います。

  1. 第1式の-2倍、-3倍、-5倍をそれぞれ第2式、第3式、第4式に代入します。
    \[ \begin{cases} \phantom{1}x_1+\phantom{1}3x_2-2x_3+2x_4=\phantom{-1}9 \\ \phantom{2x_1}-\phantom{1}7x_2+9x_3-3x_4=-15 \\ \phantom{3x_1}-\phantom{1}7x_2\mathbin{\phantom{-}}\phantom{2x_3}-\phantom{1}x_4=-18 \\ \phantom{5x_1}-14x_2+3x_3-5x_4=-39 \end{cases} \]
  2. 第2式の-1倍、-2倍をそれぞれ第3式、第4式に代入します。
    \[ \begin{cases} \phantom{1}x_1+3x_2-\phantom{1}2x_3+2x_4=\phantom{-1}9 \\ \phantom{2x_1}-7x_2+\phantom{1}9x_3-3x_4=-15 \\ \phantom{3x_1-7x_2}-\phantom{1}9x_3+6x_4=-\phantom{1}3 \\ \phantom{5x_1-4x_2}-15x_3+9x_4=-\phantom{1}9 \end{cases} \]
  3. 第3式の $\displaystyle -\frac{5}{3}$ 倍を第4式に代入します。
    \[ \begin{cases} \phantom{1}x_1+3x_2-2x_3+2x_4=\phantom{-1}9 \\ \phantom{2x_1}-7x_2+9x_3-3x_4=-15 \\ \phantom{3x_1-7x_2}-9x_3+6x_4=-\phantom{1}3 \\ \phantom{5x_1-1x_2-1x_3}-\phantom{1}x_4=-\phantom{1}4 \end{cases} \]
  4. これで前進消去は終わりましたが、第3式を $\displaystyle\frac{1}{3}$ 倍、第4式を-1倍しておきます。
    \[ \begin{cases} \phantom{1}x_1+3x_2-2x_3+2x_4=\phantom{-1}9 \\ \phantom{2x_1}-7x_2+9x_3-3x_4=-15 \\ \phantom{3x_1-7x_2}-3x_3+2x_4=-\phantom{1}1 \\ \phantom{5x_1-1x_2-1x_3}\mathbin{\phantom{-}}\phantom{1}x_4=\phantom{-1}4 \end{cases} \]

次に、後退代入を行います。

  1. 第4式を第3式に代入します。
    \[ \begin{cases} \phantom{1}x_1+3x_2-2x_3+2x_4=\phantom{-1}9 \\ \phantom{2x_1}-7x_2+9x_3-3x_4=-15 \\ \phantom{3x_1-7x_2}-3x_3\mathbin{\phantom{-}}\phantom{2x_4}=-\phantom{1}9 \\ \phantom{11x_1-1x_2-1x_3}\mathbin{\phantom{-}}x_4=\phantom{-1}4 \end{cases} \]
  2. 第3式、第4式を第2式に代入します。
    \[ \begin{cases} \phantom{1}x_1+3x_2-2x_3+2x_4=\phantom{-1}9 \\ \phantom{2x_1}-7x_2\mathbin{\phantom{-}}\phantom{9x_3-3x_4}=-14 \\ \phantom{3x_1-7x_2}-3x_3\mathbin{\phantom{-}}\phantom{2x_4}=-\phantom{1}9 \\ \phantom{11x_1-1x_2-1x_3}\mathbin{\phantom{-}}x_4=\phantom{-1}4 \end{cases} \]
  3. 第2式、第3式、第4式を第1式に代入します。
    \[ \begin{cases} \phantom{1}x_1\mathbin{\phantom{-}}\phantom{1x_2-1x_3+1x_4}=\phantom{-1}1 \\ \phantom{1x_1}-7x_2\mathbin{\phantom{-}}\phantom{1x_3-1x_4}=-14 \\ \phantom{1x_1-1x_2}-3x_3\mathbin{\phantom{-}}\phantom{1x_4}=-\phantom{1}9 \\ \phantom{11x_1-1x_2-1x_3}\mathbin{\phantom{-}}x_4=\phantom{-1}4 \end{cases} \]
  4. 解が求まりました。
    \[ \begin{cases} \phantom{1}x_1\mathbin{\phantom{-}}\phantom{1x_2-1x_3+1x_4}=\phantom{-1}1 \\ \phantom{1x_1-1}x_2\mathbin{\phantom{-}}\phantom{1x_3-1x_4}=\phantom{-1}2 \\ \phantom{1x_1-1x_2-1}x_3\mathbin{\phantom{-}}\phantom{1x_4}=\phantom{-1}3 \\ \phantom{1x_1-1x_2-1x_3-1}x_4=\phantom{-1}4 \end{cases} \]

元の連立1次方程式は次のように行列を使って表すことができます。
\[ \left( \begin{array}{rrrr} 1 & 3 & -2 & 2 \\ 2 & -1 & 5 & -3 \\ 3 & 2 & -6 & 5 \\ 5 & 1 & -7 & 5 \end{array} \right) \left( \begin{array}{r} x_1 \\ x_2 \\ x_3 \\ x_4 \end{array} \right) = \left( \begin{array}{r} 9 \\ 3 \\ 9 \\ 6 \end{array} \right) \]
ですから、逆行列を求めるためには、$AI=E$ となる行列 $I$ の第 $i$ 列を $\bvec{x}$、$E$ の第 $i$ 列を $\bvec{b}_i$ とすると、$A\bvec{x}=\bvec{b}_i$ を解けはいいことになります。

コード

やることはお分かりになったと思うので、一気にソースをご覧ください。

1ヶ所だけ上では説明していない操作がありますので、ご紹介します。

GaussPivot では、対角成分に 0 が含まれている場合は、他の行と入れ替えて対角成分から 0 を排除します。入れ替えても 0 が排除できない場合は、error フラグを立て、逆行列が存在しないことを知らせます。前回のLU分解のプログラムで実装した部分をそのまま移植したようなものなので、詳しく知りたい方は前回の記事を参照してください。

感想

わりと簡単なロジックで出来ているので、およそ中級者とはいえない私でもわりとスイスイと書けました。

このコードはご自由に使ってもらって構いません。必要ならば改変して頂いても結構です。コード内にコメントでひとこと「この記事を参考にした」などと書いていただけるとありがたいです。

そのうち行列の階数とかを計算するプログラムも考えるかもしれません。考えないかもしれません。