ガウスの消去法とは、行列を対角行列に変換することで、逆行列や連立一次方程式の解を求めるアルゴリズムです。
別名掃き出し法とも言います。
この記事では C++ での実装例を掲載します。
あわせて読みたい
目次
ガウス消去法の実装サンプル
ガウス消去法の実装サンプルを紹介します。
なお、ここでのサンプルは、なるべくコンパクトにまとめることを目的として、単一のソースコードファイルに記述しています。
Matrix クラスの定義
まず行列を表現するための Matrix クラスを定義します。
#include <iostream>
#include <iomanip>
#include <vector>
class Matrix
{
public:
Matrix(size_t n_rows, size_t n_cols, const std::vector<double> &v)
: n_rows_(n_rows),
n_cols_(n_cols)
{
v_ = v;
v_.resize(n_rows * n_cols);
}
size_t n_rows() const { return n_rows_; }
size_t n_cols() const { return n_cols_; }
const double &elem(size_t i, size_t j) const
{
return v_[i * n_cols_ + j];
}
double &elem(size_t i, size_t j)
{
return v_[i * n_cols_ + j];
}
// ２つの要素を交換する
void swap_elems(size_t i1, size_t j1, size_t i2, size_t j2)
{
double tmp = elem(i1, j1);
elem(i1, j1) = elem(i2, j2);
elem(i2, j2) = tmp;
}
// ２つの行を交換する
void swap_rows(size_t i, size_t j)
{
for (size_t k = 0; k < n_cols_; ++k)
{
swap_elems(i, k, j, k);
}
}
// 行列を標準出力する
void print() const
{
std::cout << std::fixed;
std::cout << std::setprecision(2);
for (size_t i = 0; i < n_rows_; ++i)
{
if (i == 0)
{
std::cout << "[ ";
}
else
{
std::cout << " ";
}
for (size_t j = 0; j < n_cols_; ++j)
{
std::cout << std::setw(5) << elem(i, j) << " ";
}
if (i == n_rows_ - 1)
{
std::cout << "]\n";
}
else
{
std::cout << "\n";
}
}
}
private:
std::vector<double> v_;
size_t n_rows_;
size_t n_cols_;
};ガウス消去法に必要な機能のみ実装しています。
追加で、行列の乗算演算子も実装しておきます。
アルゴリズムの結果を検証するために用います。
Matrix operator*(const Matrix &lhs, const Matrix &rhs)
{
size_t n_rows = lhs.n_rows();
size_t n_cols = rhs.n_cols();
size_t n_inner = lhs.n_cols();
Matrix out(n_rows, n_cols, {});
for (size_t i = 0; i < n_rows; ++i)
{
for (size_t j = 0; j < n_cols; ++j)
{
for (size_t k = 0; k < n_inner; ++k)
{
out.elem(i, j) += lhs.elem(i, k) * rhs.elem(k, j);
}
}
}
return out;
}ガウスの消去法
アルゴリズムの本体の実装例です。
ガウスの消去法は、前進消去と後退代入の二つの処理で構成されます。
途中経過を把握するため、要所要所に標準出力を設けています。
void gaussian_elimination(Matrix &a, Matrix &b)
{
std::cout << "starting Gaussian Elimination" << std::endl;
std::cout << "input:" << std::endl;
a.print();
b.print();
size_t n = a.n_rows();
size_t m = b.n_cols();
// 前進消去
for (size_t i = 0; i < n; ++i)
{
double max_abs = 0.0;
size_t max_row = i;
for (size_t k = i; k < n; ++k)
{
double value = std::abs(a.elem(k, i));
if (value > max_abs)
{
max_abs = value;
max_row = k;
}
}
if (i != max_row)
{
a.swap_rows(i, max_row);
b.swap_rows(i, max_row);
}
double denomi = a.elem(i, i);
for (size_t j = i; j < n; ++j)
{
a.elem(i, j) /= denomi;
}
for (size_t j = 0; j < m; ++j)
{
b.elem(i, j) /= denomi;
}
for (size_t k = i + 1; k < n; ++k)
{
double scale = a.elem(k, i);
for (size_t j = i; j < n; ++j)
{
a.elem(k, j) -= scale * a.elem(i, j);
}
for (size_t j = 0; j < m; ++j)
{
b.elem(k, j) -= scale * b.elem(i, j);
}
}
std::cout << "前進消去：ステップ " << i << std::endl;
a.print();
b.print();
}
// 後退代入
for (size_t i = n - 2;; --i)
{
for (size_t j = i + 1; j < n; ++j)
{
double scale = a.elem(i, j);
a.elem(i, j) = 0.0;
for (size_t k = 0; k < m; ++k)
{
b.elem(i, k) -= scale * b.elem(j, k);
}
}
if (i == 0)
{
break;
}
}
std::cout << "output:" << std::endl;
a.print();
b.print();
std::cout << "finished Gaussian Elimination" << std::endl;
}サンプルの実行例
連立一次方程式の求解
サンプルとして、次の一次方程式を解きます。
$$
\begin{pmatrix}
1 & 1 & 1 \\
1 & 2 & 2 \\
2 & 3 & 2
\end{pmatrix}
\begin{pmatrix}
x_0 \\ x_1 \\ x_2
\end{pmatrix}
=
\begin{pmatrix}
2 \\ 3 \\ 2
\end{pmatrix}
$$
gaussian_elimination 関数にインプットする b の内容が、解である x の値に更新されます。
結果として、解である \( x = \begin{pmatrix} 1 \\ -2 \\ 3 \end{pmatrix} \) が求められます。
int main()
{
Matrix a(3, 3, {1, 1, 1, 1, 2, 2, 2, 3, 2});
Matrix b(3, 1, {2, 3, 2});
Matrix a_copied = a;
gaussian_elimination(a, b);
std::cout << std::endl;
std::cout << "結果の検証： A * x = " << std::endl;
Matrix y = a_copied * b;
y.print();
return 0;
}starting Gaussian Elimination
input:
[ 1.00 1.00 1.00
1.00 2.00 2.00
2.00 3.00 2.00 ]
[ 2.00
3.00
2.00 ]
前進消去：ステップ 0
[ 1.00 1.50 1.00
0.00 0.50 1.00
0.00 -0.50 0.00 ]
[ 1.00
2.00
1.00 ]
前進消去：ステップ 1
[ 1.00 1.50 1.00
0.00 1.00 2.00
0.00 0.00 1.00 ]
[ 1.00
4.00
3.00 ]
前進消去：ステップ 2
[ 1.00 1.50 1.00
0.00 1.00 2.00
0.00 0.00 1.00 ]
[ 1.00
4.00
3.00 ]
output:
[ 1.00 0.00 0.00
0.00 1.00 0.00
0.00 0.00 1.00 ]
[ 1.00
-2.00
3.00 ]
finished Gaussian Elimination
結果の検証： A * x =
[ 2.00
3.00
2.00 ]逆行列の計算
ひとつ前の例と同様に、次の行列をインプットとして逆行列を計算します。
$$
A =
\begin{pmatrix}
1 & 1 & 1 \\
1 & 2 & 2 \\
2 & 3 & 2
\end{pmatrix}
$$
この場合は、gaussian_elimination 関数の引数 b に、単位行列 \( \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \) を入力します。
int main()
{
Matrix a(3, 3, {1, 1, 1, 1, 2, 2, 2, 3, 2});
Matrix b(3, 3, {1, 0, 0, 0, 1, 0, 0, 0, 1});
Matrix a_copied = a;
gaussian_elimination(a, b);
std::cout << std::endl;
std::cout << "結果の検証： A * invA = " << std::endl;
Matrix y = a_copied * b;
y.print();
return 0;
}starting Gaussian Elimination
input:
[ 1.00 1.00 1.00
1.00 2.00 2.00
2.00 3.00 2.00 ]
[ 1.00 0.00 0.00
0.00 1.00 0.00
0.00 0.00 1.00 ]
前進消去：ステップ 0
[ 1.00 1.50 1.00
0.00 0.50 1.00
0.00 -0.50 0.00 ]
[ 0.00 0.00 0.50
0.00 1.00 -0.50
1.00 0.00 -0.50 ]
前進消去：ステップ 1
[ 1.00 1.50 1.00
0.00 1.00 2.00
0.00 0.00 1.00 ]
[ 0.00 0.00 0.50
0.00 2.00 -1.00
1.00 1.00 -1.00 ]
前進消去：ステップ 2
[ 1.00 1.50 1.00
0.00 1.00 2.00
0.00 0.00 1.00 ]
[ 0.00 0.00 0.50
0.00 2.00 -1.00
1.00 1.00 -1.00 ]
output:
[ 1.00 0.00 0.00
0.00 1.00 0.00
0.00 0.00 1.00 ]
[ 2.00 -1.00 0.00
-2.00 0.00 1.00
1.00 1.00 -1.00 ]
finished Gaussian Elimination
結果の検証： A * invA =
[ 1.00 0.00 0.00
0.00 1.00 0.00
0.00 0.00 1.00 ]
コメント