私は（学習演習として）行列クラスを作成していて、逆関数をテストしているときに出くわして問題が発生しました。
私はそのように任意の行列を入力します：
2 1 1
1 2 1
1 1 2
そして逆数を計算するためにそれを手に入れました、そして私は正しい結果を得ました：
0.75 -0.25 -0.25
-0.25 0.75 -0.25
-0.25 -0.25 0.75
しかし、2つを掛け合わせて、次の単位行列を確実に取得しようとすると、次のようになります。
1 5.5111512e-017 0
0 1 0
-1.11022302e-0.16 0 1
なぜこれらの結果が得られるのですか？いくつかの丸め誤差を理解できる奇妙な数値を乗算している場合は理解できますが、その合計は次のようになります。
2 * -0.25 + 1 * 0.75 + 1 * -0.25
これは明らかに0であり、5.111512e-017ではありません
手動で計算を行うようにした場合。例えば：
std::cout << (2 * -0.25 + 1 * 0.75 + 1 * -0.25) << "\n";
期待通り0になりますか？
すべての数値はdoubleとして表されます。これが私の乗算の過負荷です：
Matrix operator*(const Matrix& A, const Matrix& B)
{
if(A.get_cols() == B.get_rows())
{
Matrix temp(A.get_rows(), B.get_cols());
for(unsigned m = 0; m < temp.get_rows(); ++m)
{
for(unsigned n = 0; n < temp.get_cols(); ++n)
{
for(unsigned i = 0; i < temp.get_cols(); ++i)
{
temp(m, n) += A(m, i) * B(i, n);
}
}
}
return temp;
}
throw std::runtime_error("Bad Matrix Multiplication");
}
およびアクセス機能：
double& Matrix::operator()(unsigned r, unsigned c)
{
return data[cols * r + c];
}
double Matrix::operator()(unsigned r, unsigned c) const
{
return data[cols * r + c];
}
逆関数を見つける関数は次のとおりです。
Matrix Inverse(Matrix& M)
{
if(M.rows != M.cols)
{
throw std::runtime_error("Matrix is not square");
}
int r = 0;
int c = 0;
Matrix augment(M.rows, M.cols*2);
augment.copy(M);
for(r = 0; r < M.rows; ++r)
{
for(c = M.cols; c < M.cols * 2; ++c)
{
augment(r, c) = (r == (c - M.cols) ? 1.0 : 0.0);
}
}
for(int R = 0; R < augment.rows; ++R)
{
double n = augment(R, R);
for(c = 0; c < augment.cols; ++c)
{
augment(R, c) /= n;
}
for(r = 0; r < augment.rows; ++r)
{
if(r == R) { continue; }
double a = augment(r, R);
for(c = 0; c < augment.cols; ++c)
{
augment(r, c) -= a * augment(R, c);
}
}
}
Matrix inverse(M.rows, M.cols);
for(r = 0; r < M.rows; ++r)
{
for(c = M.cols; c < M.cols * 2; ++c)
{
inverse(r, c - M.cols) = augment(r, c);
}
}
return inverse;
}