LU分解を用いた連立一次方程式の精度保証付き解法

  1. Home
  2. 数値計算
  3. LU分解を用いた連立一次方程式の精度保証付き解法

最終更新: 2017-06-11 00:20

連立一次方程式の精度保証付き解法を実装しました。

動機

n×nn\times n行列AAとnn次元ベクトルxxがあるときAx=bAx = bという連立一次方程式を考えることができます。AAが正則だとすると、x=A−1bx = A^{-1}bを計算すれば解を求めることができるのですが、コンピュータで計算を行う場合丸め誤差が存在するので得られる解は真の解とは異なります。

計算機で求めた解をxx、真の解をx∗x^*とすると、相対誤差∣∣x−x∗∣∣/∣∣x∗∣∣||x-x^*||/||x^*||は

∣∣x−x∗∣∣∣∣x∗∣∣=∣∣A−1(Ax−b)∣∣∣∣x∗∣∣≤∣∣A∣∣ ∣∣A−1∣∣∣∣Ax−b∣∣∣∣b∣∣\frac{||x-x^*||}{||x^*||} = \frac{||A^{-1}(Ax - b)||}{||x^*||}\leq ||A||\ ||A^{-1}||\frac{||Ax-b||}{||b||}

最後の不等式は∣∣b∣∣=∣∣Ax∗∣∣≤∣∣A∣∣∣∣x∗∣∣||b||=||Ax^*||\leq||A||||x^*||という評価を使いました。

したがって、相対誤差は残差∣∣Ax−b∣∣||Ax-b||が∣∣A∣∣ ∣∣A−1∣∣/∣∣b∣∣||A||\ ||A^{-1}|| / ||b||倍拡大されて反映されることになります。∣∣A∣∣ ∣∣A−1∣∣||A||\ ||A^{-1}||を条件数と言います。

倍精度浮動小数点数の精度は10−1610^{-16}程度なので、残差がどんなに小さくても条件数が101610^{16}より大きい場合は解の精度は1桁も存在しないことになります。したがって解の精度がどれくらいあるかを評価することが非常に重要だということが分かります。

理論

ここではLU分解を用いた連立一次方程式の精度保証をLU分解1回程度の手間で行う方法を実験します。詳細は精度保証付き数値計算 (大石 進一）の2章を読んでください。

おおよそ同じ内容が線形方程式はもう一度数値解を計算する手間で精度保証できるに書いてあります。

要点は連立一次方程式Ax=bAx = bのLU分解による近似解x~\tilde{x}と真の解x∗x^*の間に

∣∣x∗−x~∣∣∞≤(∣XU∣⋅∣XL∣)∣∣∞⋅∣∣Ax~−b∣∣∞1−(2γn∣ ∣∣XU∣ ∣XL∣ ∣L~∣ ∣U~∣ ∣∣∞+γn∣∣ ∣XU∣ ∣U~∣ ∣∣∞)||x^*- \tilde{x}||_{\infty} \leq \frac{(|X_U| \cdot |X_L|)||_{\infty} \cdot ||A\tilde{x} - b||_{\infty}}{1-(2 \gamma_n |\ ||X_U|\ |X_L|\ |\tilde{L}|\ |\tilde{U}|\ ||_{\infty} + \gamma_n ||\ |X_U|\ |\tilde{U}|\ ||_{\infty})}

という評価が成り立つということです。
ここで

この評価で一番時間がかかるのはXLX_L, XUX_Uの計算で、それぞれn3/6n^3/6回の計算が必要です。LU分解がn3/3n^3/3程度なので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

のようになり残差の他に解の誤差も評価できました。

参考ページ