87 Eigen應用：解線性方程組的古典迭代法

0 引言

線性方程組的迭代法就是用某種極限過程逐步逼近線性方程組精確解的方法。迭代法具有需要的存儲空間少、程序設計簡單、原始系數矩陣在計算過程中始終不變等優點，但有收斂性或收斂速度的問題。迭代法是解大型稀疏矩陣方程組的重要方法。迭代法的基本思想是構造一串收斂到解的序列，即建立一種從已有近似解計算新的近似解的規則。由不同的計算規則得到不同的迭代法。

本文的主要思想來自於中南大學 鄭洲順 教授在中國大學MOOC上的 科學計算與數學建模 課程， 第六章 回歸問題-線性方程組的迭代解法，鏈接如下。

https://www.icourse163.org/learn/CSU-1001985002?tid=1206784225#/learn/content?type=detail&id=1211843352&sm=1

鄭洲順 高等工程數學：
鏈接：https://pan.baidu.com/s/1F1iQ5r0z8RUUOGkKV8sUiw 
提取碼：9d9p
https://www.bilibili.com/video/BV1dE411j7G6?from=search&seid=1204602480490752881

1 問題的引入與定義

        

  求回歸系數的本質是解線性方程組。迭代法解線性方程組的基本思路如下。

      

為了應用迭代法解線性方程組，首先應該找到一個迭代矩陣B，當k->無窮大時，誤差->0，此時稱該線性方程組迭代法收斂。值得一提的是，線性方程組迭代法的斂散性與具體的迭代矩陣有關。迭代矩陣不同，線性方程組迭代法的斂散性也有可能不同。 

2 線性方程組迭代法的收斂性

 

 直接給出結論如上。即當迭代矩陣的譜半徑 < 1時，線性方程組時收斂的。由於譜半徑時矩陣范數的下限，故存在如下迭代法收斂的充分條件。   

3 評價迭代法的優良

迭代矩陣譜半徑可以衡量迭代矩陣的收斂速度，譜半徑越小，迭代速度越快。

另外一個維度是精度。在保證精度的條件下，迭代的速度越快，認為算法的效果越好。

4 解線性方程組的古典迭代法-eigen實現

#include <iostream>
#include <Eigen/Dense>  // eigen library
#include <ctime>        // clock, etc
#include <iomanip>      // std::setprecision, etc
#include <cmath>        // use fabs, pow, etc
#include <climits>


/* 迭代法解線性方程組
* Ax = b
* 該迭代法收斂的必要條件有兩條
*** 條件1：對角線元素不等於零
*** 條件2：0 < w < 2
* 迭代法解線性方程組的充要條件有一個
*** 譜半徑 < 1 
* 
*/
enum IterativeSolverMethod{
    JACOBI,               // 雅可比迭代法
    GAUSS_SEIDEL,         // 高斯迭代法
    SOR                   // 超松弛迭代法
};
/* JacobiIterative: 雅可比迭代法解線性方程組
*/
void JacobiIterative(const Eigen::MatrixXd& A,     // 系數矩陣
                     const Eigen::VectorXd& b,      // 常數項 
                     Eigen::VectorXd& x){            // 解
    // 行遍歷，每次迭代計算一個解
    Eigen::VectorXd x_temp = Eigen::VectorXd(x.size());
    for(int i = 0; i < A.rows(); ++ i){
        x_temp[i] = (b[i] - (A.row(i).dot(x) - x[i]*A(i, i)))/A(i, i);
    }
    x = x_temp;   // 全部算完再更新
}

/* GaussSeidelIterative: 高斯-賽達爾迭代法，雅克比迭代法的改進版本
* 區別在於，每次迭代計算一個新的解，均采用當前最新的結果
*/
void GaussSeidelIterative(const Eigen::MatrixXd& A,   // 系數矩陣
                          const Eigen::VectorXd& b,        // 常數項 
                          Eigen::VectorXd& x){             // 解
    // 行遍歷，每次迭代計算一個解
    for(int i = 0; i < A.rows(); ++ i){
        x[i] = (b[i] - (A.row(i).dot(x) - x[i]*A(i, i)))/A(i, i);   // 實時更新
    }
}

/* SorIterative: 給定超松弛系數的超松弛迭代法
* 當w == 1時，就是Gauss-Seidel迭代法
*/
void SorIterative(const Eigen::MatrixXd& A,      // 系數矩陣
                  const Eigen::VectorXd& b,      // 常數項 
                  const double& w,               // 松弛因子w, 通常滿足  0 < w < 2
                  Eigen::VectorXd& x){           // 解
    // 行遍歷，每次迭代計算一個解
    for(int i = 0; i < A.rows(); ++ i){
        double x_gauss_seidel = (b[i] - (A.row(i).dot(x) - x[i]*A(i, i)))/A(i, i);  // Gauss-Seidel 迭代的結果
        x[i] = x[i] + w * (x_gauss_seidel - x[i]);
    }    
}
 
/* IterativeSolver: 選擇迭代方法，根據閾值迭代求解線性方程組的解
*/
void IterativeSolver(const Eigen::MatrixXd& A,   // 系數矩陣
                     const Eigen::VectorXd& b,      // 常數項 
                     const int& threshold_iteration_times,            // 迭代次數
                     const double& threshold__iteration_error,        // 迭代誤差
                     const IterativeSolverMethod& solver,                   // 迭代方法
                     Eigen::VectorXd& x){            // 解
    // x = Eigen::VectorXd::Zero(A.rows(), 1);  // 初始解
    x = Eigen::VectorXd::Zero(A.rows());  // 初始解
    Eigen::VectorXd last_x = x;   // 保存上一次的解用於計算兩次迭代之間的誤差

    // 判斷是否滿足迭代條件
    int current_iteration_times = 1;
    double current_iteration_error = static_cast<double>(INT_MAX);      // 借用INT_MAX對初始誤差進行初始化

    const double w = 1.3;   // 在switch case外定義
    while(current_iteration_times < threshold_iteration_times &&  current_iteration_error > threshold__iteration_error){
        switch(solver){
            case JACOBI:
                JacobiIterative(A, b, x);
                break;
            case GAUSS_SEIDEL:
                GaussSeidelIterative(A, b, x);
                break;
            case SOR:
                SorIterative(A, b, w, x);
                break;
            default:
                break;
        }
        
        current_iteration_error = (x-last_x).norm();
        last_x = x;
        current_iteration_times ++;
    }
    std::cout << "current_iteration_error = (solution-last_solution).norm() = " << current_iteration_error << std::endl;
    std::cout << "current_iteration_times = " << current_iteration_times << std::endl;
}

void TestIterativeSolver(){
    const int dim = 4;
    Eigen::MatrixXd A(dim, dim);
    Eigen::VectorXd b(dim, 1);
    A << -4, 1, 1, 1,
        1, -4, 1, 1,
        1, 1, -4, 1,
        1, 1, 1, -4;
    b << 1, 1, 1, 1;

    // ldl分解求解精確值
    Eigen::VectorXd x_qr(dim, 1);
    x_qr = A.colPivHouseholderQr().solve(b);

    // 迭代法求解
    const int threshold_iteration_times = 100;        // 迭代次數
    const double threshold__iteration_error = std::pow(0.1, 4);    // 迭代誤差

    Eigen::VectorXd x_jacobi_iteration;
    IterativeSolverMethod solver = JACOBI;
    IterativeSolver(A, b, threshold_iteration_times, threshold__iteration_error, solver, x_jacobi_iteration);

    Eigen::VectorXd x_gauss_seidel_iteration;
    solver = GAUSS_SEIDEL;
    IterativeSolver(A, b, threshold_iteration_times, threshold__iteration_error, solver, x_gauss_seidel_iteration);

    Eigen::VectorXd x_sor_iteration;
    solver = SOR;
    IterativeSolver(A, b, threshold_iteration_times, threshold__iteration_error, solver, x_sor_iteration); 

    std::cout << "x_qr = " << x_qr << std::endl;
    std::cout << "x_jacobi_iteration = " << x_jacobi_iteration << std::endl;
    std::cout << "absolute error of jacobi iteration is: " << (x_jacobi_iteration - x_qr).norm() << std::endl; 

    std::cout << "x_gauss_seidel_iteration = " << x_gauss_seidel_iteration << std::endl;
    std::cout << "absolute error of gauss seidel iteration is: " << (x_gauss_seidel_iteration - x_qr).norm() << std::endl; 

    std::cout << "x_sor_iteration = " << x_sor_iteration << std::endl;
    std::cout << "absolute error of SOR iteration is: " << (x_sor_iteration - x_qr).norm() << std::endl;    
}

int main()
{
    TestIterativeSolver();
    system("pause");
    return 0;
}

5 運行結果及分析

 如上圖所示，采用Jacobi迭代法，Gauss-Seidel迭代法和超松弛迭代法解線性方程組的結果評價如下。

（1）從速度上來說，Jacobi迭代法 < Gauss-Seidel迭代法 < 超松弛迭代法

（2）從精度上來說，Jacobi迭代法 < Gauss-Seidel迭代法 < 超松弛迭代法

 因此，可以認為當w（松弛因子）選擇合適的情況下，超松弛迭代法是最好的古典迭代法。