最終更新: 2017-06-11 00:20
連立一次方程式の精度保証付き解法を実装しました。
行列と次元ベクトルがあるときという連立一次方程式を考えることができます。が正則だとすると、を計算すれば解を求めることができるのですが、コンピュータで計算を行う場合丸め誤差が存在するので得られる解は真の解とは異なります。
計算機で求めた解を、真の解をとすると、相対誤差は
最後の不等式はという評価を使いました。
したがって、相対誤差は残差が倍拡大されて反映されることになります。を条件数と言います。
倍精度浮動小数点数の精度は程度なので、残差がどんなに小さくても条件数がより大きい場合は解の精度は1桁も存在しないことになります。したがって解の精度がどれくらいあるかを評価することが非常に重要だということが分かります。
ここではLU分解を用いた連立一次方程式の精度保証をLU分解1回程度の手間で行う方法を実験します。詳細は精度保証付き数値計算 (大石 進一)の2章を読んでください。
おおよそ同じ内容が線形方程式はもう一度数値解を計算する手間で精度保証できるに書いてあります。
要点は連立一次方程式のLU分解による近似解と真の解の間に
という評価が成り立つということです。
ここで
double
なら)、はの次元この評価で一番時間がかかるのは, の計算で、それぞれ回の計算が必要です。LU分解が程度なのでLU分解を使って方程式を解く場合の倍程度の計算量で精度保証を行えることになります。
精度保証付き数値計算 (大石 進一)の2.3.2にあるMATLABコードを参考にEigenを用いたC++のコードを書きました。
#include <iostream>
#include <limits>
#include <cfenv>
#include <cmath>
#include <algorithm>
#include <Eigen/Core>
#include <Eigen/LU>
using namespace Eigen;
#pragma STDC FENV_ACCESS ON
int main() {
const int N = 1000;
const MatrixXd A = MatrixXd::Random(N, N);
const VectorXd b = VectorXd::Random(N);
fesetround(FE_UPWARD);
const double eps = std::numeric_limits<double>::epsilon()/2.0;
double nu = N * eps;
fesetround(FE_DOWNWARD);
double den = 1 - nu;
fesetround(FE_UPWARD);
double gamma = nu / den;
fesetround(FE_TONEAREST);
PartialPivLU<MatrixXd> LU(A);
// 対角成分が0になったLが得られる
MatrixXd L = LU.matrixLU().triangularView<StrictlyLower>();
// 対角成分を1にするために単位行列を足す
L = L + MatrixXd::Identity(N, N);
MatrixXd U = LU.matrixLU().triangularView<Upper>();
MatrixXd P = LU.permutationP();
MatrixXd XL = L.inverse();
MatrixXd XU = U.inverse();
VectorXd x = XU * (XL * (P * b));
fesetround(FE_UPWARD);
VectorXd rd = (A*x - b).cwiseAbs();
fesetround(FE_UPWARD);
VectorXd ru = (A*x - b).cwiseAbs();
ru = rd.cwiseMax(ru);
double r = ru.lpNorm<Infinity>();
VectorXd e = VectorXd::Ones(N);
VectorXd d1 = XU.cwiseAbs() * (XL.cwiseAbs() * (L.cwiseAbs() * (U.cwiseAbs() * e)));
VectorXd d2 = XU.cwiseAbs() * (U.cwiseAbs() * e);
double d = 2 * gamma * d1.lpNorm<Infinity>() + gamma * d2.lpNorm<Infinity>();
double R_norm = (XU.cwiseAbs() * (XL.cwiseAbs() * e)).lpNorm<Infinity>();
fesetround(FE_DOWNWARD);
double D = 1 - d;
if(D > 0) {
fesetround(FE_UPWARD);
double error = R_norm * r / D;
std::cout << "||Ax - b||:" << (A * x - b).lpNorm<Infinity>() << std::endl;
std::cout << "||x - x*||:" << error << std::endl;
} else {
std::cerr << "false" << std::endl;
}
}
実行結果は
||Ax - b||:3.26639e-12
||x - x*||:1.44821e-08
のようになり残差の他に解の誤差も評価できました。