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

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

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

LU分解
掃き出し法

があり、このうち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式の－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式の－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式の \(\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} \]
これで前進消去は終わりましたが、第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} \]

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

第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} \]
第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} \]
第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} \]
解が求まりました。
\[ \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\) 列を \(x\)、\(E\) の第 \(i\) 列を \(b_i\) とすると、\(Ax=b_i\) を解けはいいことになります。

コード

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

#include <stdio.h>
#include <math.h>

#define SIZE 4

// 掃き出し法によって逆行列を求めてみる

typedef struct {
  double data[SIZE][SIZE];
} Matrix;

typedef struct {
  double data[SIZE];
} Vector;

typedef struct {
  double data[SIZE][2*SIZE];
} Gauss;

Matrix MatrixInput();
Matrix MatrixInverse(Matrix X, int *error);
Gauss GaussCreate(Matrix X);
Gauss GaussPivot(Gauss X, int *error);
Gauss GaussElimination(Gauss X);
Gauss GaussSubstitution(Gauss X);
int MatrixPrint(Matrix X);

int main() {
  Matrix A, I;
  int err = 0;

  A = MatrixInput();

  I = MatrixInverse(A, &err);

  if (err) {
    printf("\n逆行列は存在しません\n");
  } else {
    printf("\n逆行列は\n");
    MatrixPrint(I);
  }
}

Matrix MatrixInput() {
  int i, j;
  Matrix X;

  for (i=0;i<SIZE;i++) {
    for (j=0;j<SIZE;j++) {
      printf("%d 次正方行列の (%d, %d) 成分を入力してください\n", SIZE, i+1, j+1);
      scanf("%lf", &X.data[i][j]);
    }
  }

  return X;
}

Matrix GaussExtract(Gauss X) {
  Matrix R;
  int i, j;

  for (i=0;i<SIZE;i++) {
    for (j=0;j<SIZE;j++) {
      R.data[i][j] = X.data[i][j+SIZE];
    }
  }

  return R;
}

Gauss GaussCreate(Matrix X) {
  Gauss R;
  int i, j;

  for (i=0;i<SIZE;i++) {
    for (j=0;j<SIZE;j++) {
      R.data[i][j] = X.data[i][j];
      R.data[i][j+SIZE] = (i==j ? 1 : 0);
    }
  }

  return R;
}

Gauss GaussPivot(Gauss X, int *error) {
  int i, j, k, row;
  double max, ptr;

  for (j=0;j<SIZE-1;j++) {
    max = fabs(X.data[j][j]);
    row = j;
    for (i=j+1;i<SIZE;i++) {
      if (X.data[i][j] > max) {
	max = X.data[i][j];
	row = i;
      }
    }

    if (row != j) {
      for (k=0;k<2*SIZE;k++) {
	ptr = X.data[row][k];
	X.data[row][k] = X.data[j][k];
	X.data[j][k] = ptr;
      }
    }

    if (fabs(X.data[j][j]) < 1e-10) {
      *error = 1;
    }
  }

  return X;
}

Gauss GaussElimination(Gauss X) {
  int i, j, k;
  double tmp;

  for (i=0;i<SIZE;i++) {
    tmp = X.data[i][i];
    for (j=0;j<2*SIZE;j++) {
      X.data[i][j] /= tmp;
    }
    for (j=i+1;j<SIZE;j++) {
      tmp = X.data[j][i];
      for (k=0;k<2*SIZE;k++) {
	X.data[j][k] -= X.data[i][k] * tmp;
      }
    }
  }

  return X;
}

Gauss GaussSubstitution(Gauss X) {
  int i, j, k;
  double tmp;

  for (i=SIZE-1;i>=0;i--) {
    for (j=i-1;j>=0;j--) {
      tmp = X.data[j][i];
      for (k=0;k<2*SIZE;k++) {
	X.data[j][k] -= X.data[i][k] * tmp;
      }
    }
  }

  return X;
}

Matrix MatrixInverse(Matrix X, int *error) {
  Gauss G;
  Matrix R;
  int err;

  G = GaussCreate(X);

  G = GaussPivot(G, &err);
  if (!err) {
    G = GaussElimination(G);
    G = GaussSubstitution(G);
  } else {
    *error = 1;
  }

  return GaussExtract(G);
}

int MatrixPrint(Matrix X) {
  int i, j;

  for (i=0;i<SIZE;i++) {
    for (j=0;j<SIZE;j++) {
      printf("% lf\t", X.data[i][j]);
    }
    printf("\n");
  }

  return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

#include <stdio.h>

#include <math.h>

#define SIZE 4

// 掃き出し法によって逆行列を求めてみる

typedef struct {

double data[SIZE][SIZE];

} Matrix;

typedef struct {

double data[SIZE];

} Vector;

typedef struct {

double data[SIZE][2*SIZE];

} Gauss;

Matrix MatrixInput();

Matrix MatrixInverse(Matrix X, int *error);

Gauss GaussCreate(Matrix X);

Gauss GaussPivot(Gauss X, int *error);

Gauss GaussElimination(Gauss X);

Gauss GaussSubstitution(Gauss X);

int MatrixPrint(Matrix X);

int main() {

Matrix A, I;

int err = 0;

A = MatrixInput();

I = MatrixInverse(A, &err);

if (err) {

printf("\n逆行列は存在しません\n");

} else {

printf("\n逆行列は\n");

MatrixPrint(I);

}

Matrix MatrixInput() {

int i, j;

Matrix X;

for (i=0;i<SIZE;i++) {

for (j=0;j<SIZE;j++) {

printf("%d 次正方行列の (%d, %d) 成分を入力してください\n", SIZE, i+1, j+1);

scanf("%lf", &X.data[i][j]);

}

return X;

}

Matrix GaussExtract(Gauss X) {

Matrix R;

int i, j;

for (i=0;i<SIZE;i++) {

for (j=0;j<SIZE;j++) {

R.data[i][j] = X.data[i][j+SIZE];

}

return R;

}

Gauss GaussCreate(Matrix X) {

Gauss R;

int i, j;

for (i=0;i<SIZE;i++) {

for (j=0;j<SIZE;j++) {

R.data[i][j] = X.data[i][j];

R.data[i][j+SIZE] = (i==j ? 1 : 0);

}

return R;

}

Gauss GaussPivot(Gauss X, int *error) {

int i, j, k, row;

double max, ptr;

for (j=0;j<SIZE-1;j++) {

max = fabs(X.data[j][j]);

row = j;

for (i=j+1;i<SIZE;i++) {

if (X.data[i][j] > max) {

max = X.data[i][j];

row = i;

}

if (row != j) {

for (k=0;k<2*SIZE;k++) {

ptr = X.data[row][k];

X.data[row][k] = X.data[j][k];

X.data[j][k] = ptr;

}

if (fabs(X.data[j][j]) < 1e-10) {

*error = 1;

}

return X;

}

Gauss GaussElimination(Gauss X) {

int i, j, k;

double tmp;

for (i=0;i<SIZE;i++) {

tmp = X.data[i][i];

for (j=0;j<2*SIZE;j++) {

X.data[i][j] /= tmp;

}

for (j=i+1;j<SIZE;j++) {

tmp = X.data[j][i];

for (k=0;k<2*SIZE;k++) {

X.data[j][k] -= X.data[i][k] * tmp;

}

return X;

}

Gauss GaussSubstitution(Gauss X) {

int i, j, k;

double tmp;

for (i=SIZE-1;i>=0;i--) {

for (j=i-1;j>=0;j--) {

tmp = X.data[j][i];

for (k=0;k<2*SIZE;k++) {

X.data[j][k] -= X.data[i][k] * tmp;

}

return X;

}

Matrix MatrixInverse(Matrix X, int *error) {

Gauss G;

Matrix R;

int err;

G = GaussCreate(X);

G = GaussPivot(G, &err);

if (!err) {

G = GaussElimination(G);

G = GaussSubstitution(G);

} else {

*error = 1;

}

return GaussExtract(G);

}

int MatrixPrint(Matrix X) {

int i, j;

for (i=0;i<SIZE;i++) {

for (j=0;j<SIZE;j++) {

printf("% lf\t", X.data[i][j]);

}

printf("\n");

}

return 0;

}

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

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

感想

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

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

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

Share on Tumblr

Web being Web

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

原理

コード

感想

関連記事

2024年4月
日	月	火	水	木	金	土
« 1月
	1	2	3	4	5	6
7	8	9	10	11	12	13
14	15	16	17	18	19	20
21	22	23	24	25	26	27
28	29	30