C言語で逆行列を求める

C言語で基本変形を利用したn次元の逆行列の求め方を紹介します。

この方法は掃き出し法といい、基本原理はC言語で連立一次方程式を解くと同じになります。

基本変形を用いて逆行列を求めるには、行列Aと単位行列Eを用いて、

この性質を利用します。下にコードを示します。

int i, j, k, n;

/*
  例として
  A = 
  [1  1  1  1]
  [2  1  2  1]
  [1  2  3 -4]
  [1 -1 -1  1]
  の逆行列を求めます。
  [A E] → [E X]
  このときX =A^-1であることを利用します。
*/

// N行 2N列の行列 M = [A E]
double M[N][N * 2] =
{
	{1,	 1,  1,  1,    1, 0, 0, 0},
	{2,	 1,  2,  1,    0, 1, 0, 0},
	{1,	 2,  3, -4,    0, 0, 1, 0},
	{1,	-1, -1,  1,    0, 0, 0, 1}
};

double pivot, mul;

// 対角成分が1で正規化された階段行列を作る
for (i = 0; i < N; ++i)
{
	// 対角成分の選択、この値で行成分を正規化
	pivot = M[i][i];
	for (j = 0; j < N * 2; ++j)
	{
		M[i][j] = (1 / pivot) * M[i][j];
	}

	// 階段行列を作る為に、現在の行より下の行につ
	// i列目の成分が0になるような基本変形をする
	for (k = i + 1; k < N; ++k)
	{
		mul = M[k][i];
		for (n = i; n < N * 2; ++n)
		{
			M[k][n] = M[k][n] -  mul * M[i][n];
		}
	}
}

// 下から上に向かって変数に代入して、左側について
for (i = N - 1; i > 0; --i)
{
	for (k = i - 1; k >= 0; --k)
	{
		mul = M[k][i];
		for (n = i; n < N * 2; ++n)
		{
			M[k][n] = M[k][n] - mul * M[i][n];
		}
	}
}

// ここで、M = [E X] となっている

/* 実行結果:
   M =
   [1 0 0 0    0.5  -0.25   0.25   0.75]
   [0 1 0 0     2   -1.25   0.25   0.25]
   [0 0 1 0   -1.5   1.25  -0.25  -0.75]
   [0 0 0 1     0    0.25  -0.25  -0.25]
   よって
   A^-1 =
   [  0.5  -0.25   0.25   0.75]
   [   2   -1.25   0.25   0.25]
   [ -1.5   1.25  -0.25  -0.75]
   [   0    0.25  -0.25  -0.25]
*/

このように結果は変形後の行列Mの右側に現れ、これによってn行n列の逆行列を得ることができます。基本変形は行列計算をする上で非常に便利で重要な性質であるということがわかりますね。

C言語で連立一次方程式を解く

C言語で掃き出し法を用いて連立一次方程式を解きます。

掃き出し法はガウスの消去法とも呼ばれるもので、線形代数学を学んだことがある人はご存知かと思いますが、簡単に説明すると、連立一次方程式について係数拡大行列を作り、対角成分のみを1とし、他の成分を0とすると対応する変数について解が求まるというものですね。

要するに、

これを行列で表して、

左辺の行列と右辺のベクトルを合わせて拡大係数行列をつくると、

これを変形していきます。今回は行列のrankは3である(x, y, zはそれぞれ独立した解をもつ)場合にのみ限った連立方程式を解いていくと前提にします。すると、基本変形を駆使して最終的には

の形になり、x = l, y = m, z = nが求まるという仕組みです。これをC言語で表現してみましょう。

#define N 4

int i, j, k, n;

/*
  例として、
   x +  y +  z +  w = 10
  2x +  y + 2z +  w = 14
   x + 2y + 3z - 4w = -2
   x -  y -  z +  w =  0
  を解く。
*/

// 拡大係数行列 M
double M[N][N + 1] =
{
	{1,	 1,  1,  1, 10},
	{2,	 1,  2,  1, 14},
	{1,	 2,  3, -4, -2},
	{1,	-1, -1,  1,  0}
};

double pivot, mul;

// 対角成分が1で正規化された階段行列を作る(前進消去)
for (i = 0; i < N; ++i)
{
	// 対角成分の選択、この値で行成分を正規化
	pivot = M[i][i];
	for (j = 0; j < N + 1; ++j)
	{
		M[i][j] = (1 / pivot) * M[i][j];
	}

	// 階段行列を作る為に、現在の行より下の行について
	// i列目の成分が0になるような基本変形をする
	for (k = i + 1; k < N; ++k)
	{
		mul = M[k][i];
		for (n = i; n < N + 1; ++n)
		{
			M[k][n] = M[k][n] -  mul * M[i][n];
		}
	}
}

// 下から上に向かって変数に代入して、独立した解の形にする(後進代入)
// このとき一番下の行はすでに独立した解を得ている
for (i = N - 1; i > 0; --i)
{
	for (k = i - 1; k >= 0; --k)
	{
		mul = M[k][i];
		for (n = i; n < N + 1; ++n)
		{
			M[k][n] = M[k][n] - mul * M[i][n];
		}
	}
}

/* 実行結果:
   M =
   [1 0 0 0 1]
   [0 1 0 0 2]
   [0 0 1 0 3]
   [0 0 0 1 4]
   これより、x = 1, y = 2, z = 3, w = 4
*/