大学の演習で逆行列を求めるプログラムを作る課題があったんですが、「2次」と決められていたのでめんどくさくてたすき掛けで作って提出してしまいました。

それだけではつまらないので、3次以上でも適用できるメソッドを実装してみることにしました。

候補は次の2つです。

• LU分解
• 掃き出し法

今回はLU分解を選択し、これを実装することにします。掃き出し法はまた今度。

[12/29追記] 掃き出し法についても記事を書きました。
逆行列を求めるプログラムの考察 (2) — 掃き出し法

ソースを示す前にLU分解について軽く説明します。

原理

LU分解とは、\(n\) 次正方行列 \(A\) があったときに、適当な上三角行列 \(U\) と下三角行列 \(L\) があって \(A=LU\) と分解できるという方法です。

\( \left( \begin{array}{ccc} a_{11} & \ldots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{n1} & \ldots & a_{nn} \end{array} \right) = \left( \begin{array}{ccc} l_{11} & \ldots & 0 \\ \vdots & \ddots & \vdots \\ l_{n1} & \ldots & l_{nn} \end{array} \right) \left( \begin{array}{ccc} u_{11} & \ldots & u_{1n} \\ \vdots & \ddots & \vdots \\ 0 & \ldots & u_{nn} \end{array} \right) \)

もちろん \(n\) 次の行列なので、式は \(n^2\) 個、未知係数は \(n(n+1)\) 個あるのでこのままでは解けません。そこで、次の2つの方法を利用します。

Doolittle 法

\(L\) の対角成分すべてを 1 とおく。

Crout 法

\(U\) の対角成分すべてを 1 とおく。

Doolittle 法を例にとって具体的な操作手順を紹介します。

1. \(L\) の対角成分すべてを 1 とおく。 \( L = \left( \begin{array}{ccc} l_{11} & \ldots & 0 \\ \vdots & \ddots & \vdots \\ l_{n1} & \ldots & l_{nn} \end{array} \right) \)
2. \(A=LU\) の両辺の係数を比較する。 \( \begin{cases} a_{1k} = u_{1k} \\ \vdots \end{cases} \) めんどくさいから省略しました。
3. 上から順番に解いていくと、\(a_{ij}\) が決定される。

しかし、逆行列の計算をする上で、\(LU\) の対角成分に 0 があるととても困ります。そこで、対角成分に 0 があるときは 0 のない行と入れ替え、置換行列 \(P\) にその情報を覚えさせておきます。そうすればもとの行列はあとから復元可能です（注：今回は逆行列の計算しか行わないので不要ですが、連立1次方程式の解を求める場合は \(P\) を持っておくことで解が入れ替わってしまう可能性を排除できます）。

なお、以下では特に注記がなければ行列の添字（1, 2, …, \(n\)）をもとに話を進めます。配列の添字については i というように等幅フォントで明示します。

ソース

C言語で実装した例を次に挙げます。なお、main 関数以外はほとんど C# – 行列分解 をC#から移植しただけです。

C#と違い、C言語は多次元配列をスムーズに扱えるように設計されていません。そのため、\(n\) 次正方行列を Matrix、\(n\) 次列ベクトルを bvector という構造体を定義して扱っています。

それでは各部分の動作を見てみましょう。

MatrixCreate で行列を初期化します。A は入力行列、I は逆行列です。 MatrixInput でキーボードから行列 A を入力させます。

MatrixInverse で逆行列を求め、MatrixPrint で表示しています。なお、MatrixInverse に渡した変数 err が 1 になっている場合は \(\det A = 0\) ですから、逆行列が存在しない旨を表示して終了します。

LUP分解と逆行列の計算

MatrixInverse とそれに関係するコードを深く掘り下げてみましょう。

MatrixDecompose でLUP分解を行います。

1. 置換行列 \(P\) は理論上は \(n\) 次正方行列ですが、ここでは入れ替えた行が追跡できれば目的を満たすので、第 \(i\) 成分が \(i\) である \(n\) 次ベクトルになっています。
2. toggle は線形代数学でいうところの「置換の符号」です。1 ならば偶置換（すなわち、行入れ替えが偶数回行われた）、-1 ならば奇置換です。これは行列式を求める場合に使います。
3. 各列について、対角成分の絶対値が最大となるように入れ替えを行います。colMax に初めの対角成分の絶対値を格納し、以降の行でそれより大きい成分があればその行を row に記録しておいて、探索が終わったら row 行と入れ替えます。
4. 行を入れ替えても対角成分が 0 ならば逆行列は存在しないので、error フラグを立てて終了します。単に行列分解を行う（逆行列を求めることを目的としない）ならばこの部分は必要ありません。
5. 対角成分が 0 でないならば、分解を行います。ここから先は少々複雑です。
   1. \(A=LU\) の係数を比較することで、次の方程式を得ます。\( a_{ij} = \sum_{k=1}^n l_{ik} u_{kj} \)
   2. いま、\(j=1\) とします。\(a_{11} = l_{11} u_{11} = u_{11}\) で、\(i \geq 2\) に対し \(a_{i1} = l_{i1} u_{11} = l_{i1} a_{11}\) ですから、\(a_{i1} / a_{11}\) を実行すれば A の第1列は L の第1列と（(1,1) 成分以外は）等しくなります。
   3. さらに、\(a_{1k} = l_{11} u_{1k} = u_{1k}\) を用いると、\(a_{2k} = l_{21} u_{1k} + u_{2k}\) ですから、\(a_{2k} – (a_{i1} / a_{11}) a_{1k}\) を実行すれば、Aの第2列は U の第2列と（(2,1) 成分以外は）等しくなります。
   4. これらを順次実行していけば、\( \left( \begin{array}{cccc} u_{11} & u_{12} & \cdots & u_{1n} \\ l_{21} & u_{22} & \cdots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & u_{nn} \end{array} \right) \) のような行列が完成します。

逆行列の計算――連立1次方程式の解

つづいて、HelperSolve で連立1次方程式の解を求めます。逆行列を求めるのに連立1次方程式？と思うかもしれませんが、置換行列の第 \(i\) 列を \(\vec{b}\) として、\(A\vec{x}=\vec{b}\) なる \(\vec{b}\) を求め、これを逆行列 \(I\) の第 \(i\) 列とするということです。

\(\vec{b}\) を作るには、置換行列 \(P\) から第 \(i\) 列を取り出します。つまり、perm[j] が i ならば、\(\vec{b}\) の第 \(j\) 成分は 1、そうでなければ 0 とすればよいのです。たとえば i=0 で、 perm[j] = {2, 1, 0, 3} なら、\( \vec{b} = \left( \begin{array}{c} 0 \\ 0 \\ 1 \\ 0 \end{array} \right) \) です。

\(A=LU\) ですから、\(A\vec{x}=\vec{b}\) はすなわち \(LU\vec{x}=\vec{b}\) です。ここで、\(U\vec{x}=\vec{c}\) とすると、\(L\vec{c}=\vec{b}\) ですから、まず \(\vec{c}\) を求め、次に \(\vec{c}\) を用いて \(\vec{x}\) を求めます。これが HelperSolve の主な処理です。

詳しく見てみましょう。

1. \(\vec{x}\) に \(\vec{b}\) を代入します。
2. まず \(L\vec{c}=\vec{b}\) を解きます。\(l_{11} c_1 = c_1 = b_1\) です。次に、\(l_{21} c_1 + l_{22} c_2 = l_{21} c_1 + c_2 = b_2\) ですから、\(c_2 = b_2 – A_{21} b_1\) として計算できます。以下、前の計算で求まった \(c_i\) を次の方程式に代入していきます（これを前進代入といいます）。
3. 前進代入が終わったら、\(U\vec{x} = \vec{c}\) を解きます。まず、\(u_{nn} x_n = c_n\) ですから、\(x_n = c_n / A_{nn}\) です。次に、\(u_{n-1,n-1} x_{n-1} + u_{n-1,n} x_n = c_{n-1}\) ですから、\(x_{n-1} = (c_{n-1} – A_{n-1,n} x_n) / A_{n-1,n-1}\) として計算できます。以下、前の計算で求まった \(x_i\) を次の方程式に代入していきます（これを後退代入といいます）。

以上の処理を終えると、MathInverse 内で得られた R が求める逆行列ということになります。

感想

私がなぜ演習の課題としてこれを提出するのを諦めたかというと、これだけ洗練されたプログラムは自分一人では到底書けなかったからです。実際、このプログラムを移植したはいいものの、仕組みを理解するのにかなり時間がかかりました。

本記事は線形代数学（特に行列の演算）についてある程度の知識を必要とします。詳しく学びたい方は 線型代数のビブリオグラフィー [数学についてのwebノート] を参考に、大きな書店や図書館で参考書を探してみてください。ちなみに私は『理系のための線型代数の基礎』で学びました。

掃き出し法の方も出来上がったら記事を書きます。