C++でガウス・ジョルダン法を使って行列の逆行列を求める方法

この問題では、2次元行列 mat[][] が与えられ、ガウス・ジョルダン法(Gauss-Jordan Method)を用いて行列の逆行列を求めることが課題となります。

まず、問題を理解するための基礎知識から確認していきましょう。

行列(MATRIX)とは、数値を2次元配列の形に並べたものです。

例

$\begin{bmatrix}2&5&4 \\1&6&7 \\9&3&8\end{bmatrix}$

逆行列 [A^-1]とは、正方行列に対して行われる演算のひとつです。行列が逆行列を持つためには、次の条件を満たす必要があります。

• 元の行列が正方行列であること。

• 正則(非特異)行列であること。

• 行列Aに対する単位行列Iが存在し、次の式が成り立つこと。

$$AA^{-1} = A^{-1}.A = I$$

任意の行列の逆行列を求めるために使える公式として、次のものが知られています。

$A^{-1}\:=\:\left(\frac{adj(A)}{\det(A)}\right)$

adj(A) は行列Aの余因子行列(随伴行列)、det(A) は行列Aの行列式を表します。

逆行列を求める方法は複数ありますが、本記事では基本変形(Elementary Row Operation)とも呼ばれるガウス・ジョルダン法について解説します。

これは、行列の逆行列を段階的に求めていく手法です。具体的な手順は以下の通りです。

• 単位行列を用いて拡大行列を作成する。

• 手順1で作成した拡大行列に行基本変形(掃き出し操作)を施し、階段(エシェロン)形式へと変形していく。

• 変形の過程で拡大行列に施すことができる操作には、次のようなものがあります。

  • 行の入れ替え(任意の2つの行を交換できます)

  • スカラー倍(行の各要素を0以外の定数倍できます)

  • 行の加算(ある行を「その行」と「別の行の定数倍」の和で置き換えます)

実装例

この解法の動作を示すサンプルプログラムです。

#include <iostream>
#include <vector>
using namespace std;
void printMatrixValues(float** arr, int n, int m){
   for (int i = 0; i < n; i++) {
      for (int j = 0; j < m; j++) {
         cout<<arr[i][j]<<"\t";
      }
      cout<<endl;
   }
   return;
}
void printInverseMatrix(float** arr, int n, int m){
   for (int i = 0; i < n; i++) {
      for (int j = n; j < m; j++) {
         printf("%.3f\t", arr[i][j]);
      }
      cout<<endl;
   }
   return;
}
void findInvMatGaussJordan(float** mat, int order){
   float temp;
   printf("The inverse of matrix : A = \n");
   printMatrixValues(mat, order, order);
   for (int i = 0; i < order; i++) {
      for (int j = 0; j < 2 * order; j++) {
         if (j == (i + order))
            mat[i][j] = 1;
      }
   }
   for (int i = order - 1; i > 0; i--) {
      if (mat[i - 1][0] < mat[i][0]) {
         float* temp = mat[i];
         mat[i] = mat[i - 1];
         mat[i - 1] = temp;
      }
   }
   for (int i = 0; i < order; i++) {
      for (int j = 0; j < order; j++) {
         if (j != i) {
            temp = mat[j][i] / mat[i][i];
            for (int k = 0; k < 2 * order; k++) {
               mat[j][k] -= mat[i][k] * temp;
            }
         }
      }
   }
   for (int i = 0; i < order; i++) {
      temp = mat[i][i];
      for (int j = 0; j < 2 * order; j++) {
         mat[i][j] = mat[i][j] / temp;
      }
   }
   cout<<"A' =\n";
   printInverseMatrix(mat, order, 2 * order);
   return;
}
int main(){
   int order = 3;
   float** mat = new float*[20];
   for (int i = 0; i < 20; i++)
   mat[i] = new float[20];
   mat[0][0] = 6; mat[0][1] = 9; mat[0][2] = 5;
   mat[1][0] = 8; mat[1][1] = 3; mat[1][2] = 2;
   mat[2][0] = 1; mat[2][1] = 4; mat[2][2] = 7;
   findInvMatGaussJordan(mat, order);
   return 0;
}

アルゴリズムのポイント

プログラムでは、入力行列の右側に単位行列を結合した n×2n の拡大行列を作成し、左半分が単位行列になるまで掃き出し処理を繰り返します。すべての処理が完了した時点での右半分が、求める逆行列 A⁻¹ となります。計算量は O(n³) であり、行列のサイズがある程度大きくなっても効率的に動作します。なお、対角成分(ピボット)が0に近い場合は行の入れ替えが必要になるなど、数値的な安定性への配慮も重要です。

出力

The inverse of matrix : A =
6 9 5
8 3 2
1 4 7
A' =
-0.049  0.163  -0.011
0.205  -0.141  -0.106
-0.110  0.057  0.205