前回の記事ではガウス消去法のアルゴリズムを紹介しました。
当アルゴリズムの派生で、ガウス・ジョルダン消去法というのがあります。
ガウス消去法では、前進消去と後退代入と２つの処理ステップで構成されます。
ガウス・ジョルダン消去法では、前進消去と後退代入を同時に進めるようなイメージのアルゴリズムです。
ガウス消去法と比較すると、計算量が大きくなってしまうデメリットがあります。
一方で、入力行列を直接逆行列に変換できることから、メモリ節約のメリットがあります。
この記事では、ガウス・ジョルダン消去法の実装サンプルを掲載します。
がりテック
C++ | ガウス消去法で逆行列や連立一次方程式の解を求める | がりテック
ガウスの消去法とは、行列を対角行列に変換することで、逆行列や連立一次方程式の解を求めるアルゴリズムです。別名掃き出し法とも言います。 この記事では C++ での実装例…
目次
ガウス・ジョルダン消去法
Matrix クラス
Matrix クラスについては、前回定義したものを踏襲します。
今回のアルゴリズムでは、行列の列交換も必要であるため、swap_columns メソッドを追加しています。
#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 swap_columns(size_t i, size_t j)
{
for (size_t k = 0; k < n_rows_; ++k)
{
swap_elems(k, i, k, j);
}
}
// 行列を標準出力する
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;
}ガウス・ジョルダン消去法
ガウス・ジョルダン消去法の実装例です。
前回記事のガウス消去法では、入力行列の A を単位行列となるよう変換しました。
以下のアルゴリズムでは、A が直接逆行列となるよう変換します。
void gauss_jordan_elimination(Matrix &a)
{
std::cout << "Starging ガウス・ジョルダン消去法" << std::endl;
std::cout << "入力行列：" << std::endl;
a.print();
size_t n = a.n_rows();
std::vector<size_t> col_indices;
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;
}
}
col_indices.push_back(max_row);
if (i != max_row)
{
a.swap_rows(i, max_row);
}
double denomi = a.elem(i, i);
a.elem(i, i) = 1.0;
for (size_t j = 0; j < n; ++j)
{
a.elem(i, j) /= denomi;
}
for (size_t k = 0; k < n; ++k)
{
if (k == i)
{
continue;
}
double scale = a.elem(k, i);
a.elem(k, i) = 0.0;
for (size_t j = 0; j < n; ++j)
{
a.elem(k, j) -= scale * a.elem(i, j);
}
}
std::cout << "ガウス・ジョルダン消去：ステップ " << i << std::endl;
a.print();
}
for (size_t i = n - 1;; --i)
{
if (col_indices[i] != i)
{
a.swap_columns(col_indices[i], i);
}
if (i == 0)
{
break;
}
}
std::cout << "出力行列：" << std::endl;
a.print();
std::cout << "Finished ガウス・ジョルダン消去法" << std::endl;
}実行サンプル
前回記事同様に、次の行列をインプットとして逆行列を計算します。
$$
A =
\begin{pmatrix}
1 & 1 & 1 \\
1 & 2 & 2 \\
2 & 3 & 2
\end{pmatrix}
$$
int main()
{
Matrix a(3, 3, {1, 1, 1, 1, 2, 2, 2, 3, 2});
Matrix a_copied = a;
gauss_jordan_elimination(a);
std::cout << std::endl;
std::cout << "結果の検証： A * invA = " << std::endl;
Matrix y = a_copied * a;
y.print();
return 0;
}Starging ガウス・ジョルダン消去法
入力行列：
[ 1.00 1.00 1.00
1.00 2.00 2.00
2.00 3.00 2.00 ]
ガウス・ジョルダン消去：ステップ 0
[ 0.50 1.50 1.00
-0.50 0.50 1.00
-0.50 -0.50 0.00 ]
ガウス・ジョルダン消去：ステップ 1
[ 2.00 -3.00 -2.00
-1.00 2.00 2.00
-1.00 1.00 1.00 ]
ガウス・ジョルダン消去：ステップ 2
[ 0.00 -1.00 2.00
1.00 0.00 -2.00
-1.00 1.00 1.00 ]
出力行列：
[ 2.00 -1.00 0.00
-2.00 0.00 1.00
1.00 1.00 -1.00 ]
Finished ガウス・ジョルダン消去法
結果の検証： A * invA =
[ 1.00 0.00 0.00
0.00 1.00 0.00
0.00 0.00 1.00 ]
コメント