題目回顧:
一般解線性方程Ax=b有哪幾種做法?你能在Eigen中實現嗎?
解: 線性方程組Ax = b的解法 :
1、直接法：(1,2,3,4,5）
2、迭代法:如Jacobi迭代法(6)
其中只有2 3方法不要求方程組個數與變量個數相等
下面簡略說明下Jacobi迭代算法:
由迭代法求解線性方程組的基本思想是將聯立方程組的求解歸結為重復計算一組彼此獨立的線性表達式，這就使問題得到了簡化，類似簡單迭代法轉換方程組中每個方程式可得到雅可比迭代式
迭代法求解方程組有一定的局限性，比如下面Jacobi函數程序實現的過程中，會判斷迭代矩陣的譜半徑是不是小於1，如果小於1表示Jacobi迭代法收斂，如果求不出來迭代矩陣即D矩陣不可逆的話，無法通過收斂的充要條件來判斷是不是收斂，此時可以試着迭代多次，看看輸出結果是否收斂，此時Jacobi迭代算法並不一定收斂，只能試着做下，下面的程序實現過程僅僅處理了充要條件，迭代法同時有十分明顯的優點——算法簡單，因而編制程序比較容易，所以在實際求解問題中仍有非常大利用價值。
具體代碼實現 如下:
1 #include<iostream> 2 #include<ctime> 3 #include <cmath> 4 #include <complex> 5 /* 線性方程組Ax = b的解法 ( 直接法（1,2,3,4,5）+迭代法(6) ) 其中只有2 3方法不要求方程組個數與變量個數相等 */ 6 7 //包含Eigen頭文件 8 //#include <Eigen/Dense> 9 #include<Eigen/Core> 10 #include<Eigen/Geometry> 11 #include <Eigen/Eigenvalues> 12 13 //下面這兩個宏的數值一樣的時候 方法1 4 5 6才能正常工作 14 #define MATRIX_SIZE 3 //方程組個數 15 #define MATRIX_SIZE_ 3 //變量個數 16 //using namespace std; 17 typedef Eigen::Matrix<double,MATRIX_SIZE, MATRIX_SIZE_> Mat_A; 18 typedef Eigen::Matrix<double ,MATRIX_SIZE,1 > Mat_B; 19 20 //Jacobi迭代法的一步求和計算 21 double Jacobi_sum(Mat_A &A,Mat_B &x_k,int i); 22 23 //迭代不收斂的話 解向量是0 24 Mat_B Jacobi(Mat_A &A,Mat_B &b, int &iteration_num, double &accuracy ); 25 26 int main(int argc,char **argv) 27 { 28 //設置輸出小數點后3位 29 std::cout.precision(3); 30 //設置變量 31 Eigen::Matrix<double,MATRIX_SIZE, MATRIX_SIZE_> matrix_NN = Eigen::MatrixXd::Random(MATRIX_SIZE,MATRIX_SIZE_); 32 Eigen::Matrix<double ,MATRIX_SIZE,1 > v_Nd = Eigen::MatrixXd::Random(MATRIX_SIZE,1); 33 34 //測試用例 35 matrix_NN << 10,3,1,2,-10,3,1,3,10; 36 v_Nd <<14,-5,14; 37 38 //設置解變量 39 Eigen::Matrix<double,MATRIX_SIZE_,1>x; 40 41 //時間變量 42 clock_t tim_stt = clock(); 43 44 /*1、求逆法 很可能沒有解 僅僅針對方陣才能計算*/ 45 #if (MATRIX_SIZE == MATRIX_SIZE_) 46 x = matrix_NN.inverse() * v_Nd; 47 std::cout<<"直接法所用時間和解為："<< 1000*(clock() - tim_stt)/(double)CLOCKS_PER_SEC 48 <<"MS"<< std::endl << x.transpose() << std::endl; 49 #else 50 std::cout<<"直接法不能解!(提示:直接法中方程組的個數必須與變量個數相同，需要設置MATRIX_SIZE == MATRIX_SIZE_)"<<std::endl; 51 #endif 52 53 /*2、QR分解解方程組 適合非方陣和方陣 當方程組有解時的出的是真解，若方程組無解得出的是近似解*/ 54 tim_stt = clock(); 55 x = matrix_NN.colPivHouseholderQr().solve(v_Nd); 56 std::cout<<"QR分解所用時間和解為："<<1000*(clock() - tim_stt)/(double)CLOCKS_PER_SEC 57 << "MS" << std::endl << x.transpose() << std::endl; 58 59 /*3、最小二乘法 適合非方陣和方陣，方程組有解時得出真解，否則是最小二乘解(在求解過程中可以用QR分解 分解最小二成的系數矩陣) */ 60 tim_stt = clock(); 61 x = (matrix_NN.transpose() * matrix_NN ).inverse() * (matrix_NN.transpose() * v_Nd); 62 std::cout<<"最小二乘法所用時間和解為:"<< 1000*(clock() - tim_stt)/(double)CLOCKS_PER_SEC 63 << "MS" << std::endl << x.transpose() << std::endl; 64 65 /*4、LU分解方法 只能為方陣（滿足分解的條件才行） */ 66 #if (MATRIX_SIZE == MATRIX_SIZE_) 67 tim_stt = clock(); 68 x = matrix_NN.lu().solve(v_Nd); 69 std::cout<< "LU分解方法所用時間和解為:" << 1000*(clock() - tim_stt)/(double)CLOCKS_PER_SEC 70 << "MS" << std::endl << x.transpose() << std::endl; 71 #else 72 std::cout<<"LU分解法不能解!(提示:直接法中方程組的個數必須與變量個數相同，需要設置MATRIX_SIZE == MATRIX_SIZE_)"<<std::endl; 73 #endif 74 75 /*5、Cholesky 分解方法 只能為方陣 (結果與其他的方法差好多)*/ 76 #if (MATRIX_SIZE == MATRIX_SIZE_) 77 tim_stt = clock(); 78 x = matrix_NN.llt().solve(v_Nd); 79 std::cout<< "Cholesky 分解方法所用時間和解為:" << 1000*(clock() - tim_stt)/(double)CLOCKS_PER_SEC 80 << "MS"<< std::endl<< x.transpose()<<std::endl; 81 #else 82 std::cout<< "Cholesky法不能解!(提示:直接法中方程組的個數必須與變量個數相同，需要設置MATRIX_SIZE == MATRIX_SIZE_)"<<std::endl; 83 #endif 84 85 /*6、Jacobi迭代法 */ 86 #if (MATRIX_SIZE == MATRIX_SIZE_) 87 int Iteration_num = 10 ; 88 double Accuracy =0.01; 89 tim_stt = clock(); 90 x= Jacobi(matrix_NN,v_Nd,Iteration_num,Accuracy); 91 std::cout<< "Jacobi 迭代法所用時間和解為:" << 1000*(clock() - tim_stt)/(double)CLOCKS_PER_SEC 92 << "MS"<< std::endl<< x.transpose()<<std::endl; 93 #else 94 std::cout<<"LU分解法不能解!(提示:直接法中方程組的個數必須與變量個數相同，需要設置MATRIX_SIZE == MATRIX_SIZE_)"<<std::endl; 95 #endif 96 97 return 0; 98 } 99 100 //迭代不收斂的話 解向量是0 101 Mat_B Jacobi(Mat_A &A,Mat_B &b, int &iteration_num, double &accuracy ) { 102 Mat_B x_k = Eigen::MatrixXd::Zero(MATRIX_SIZE_,1);//迭代的初始值 103 Mat_B x_k1; //迭代一次的解向量 104 int k,i; //i,k是迭代算法的循環次數的臨時變量 105 double temp; //每迭代一次解向量的每一維變化的模值 106 double R=0; //迭代一次后，解向量每一維變化的模的最大值 107 int isFlag = 0; //迭代要求的次數后，是否滿足精度要求 108 109 //判斷Jacobi是否收斂 110 Mat_A D; //D矩陣 111 Mat_A L_U; //L+U 112 Mat_A temp2 = A; //臨時矩陣獲得A矩陣除去對角線后的矩陣 113 Mat_A B; //Jacobi算法的迭代矩陣 114 Eigen::MatrixXcd EV;//獲取矩陣特征值 115 double maxev=0.0; //最大模的特征值 116 int flag = 0; //判斷迭代算法是否收斂的標志 1表示Jacobi算法不一定能收斂到真值 117 118 std::cout<<std::endl<<"歡迎進入Jacobi迭代算法！"<<std::endl; 119 //對A矩陣進行分解 求取迭代矩陣 再次求取譜半徑 判斷Jacobi迭代算法是否收斂 120 for(int l=0 ;l < MATRIX_SIZE;l++){ 121 D(l,l) = A(l,l); 122 temp2(l,l) = 0; 123 if(D(l,l) == 0){ 124 std::cout<<"迭代矩陣不可求"<<std::endl; 125 flag =1; 126 break; 127 } 128 } 129 L_U = -temp2; 130 B = D.inverse()*L_U; 131 132 //求取特征值 133 Eigen::EigenSolver<Mat_A>es(B); 134 EV = es.eigenvalues(); 135 // cout<<"迭代矩陣特征值為:"<<EV << endl; 136 137 //求取矩陣的特征值 然後獲取模最大的特徵值 即爲譜半徑 138 for(int index = 0;index< MATRIX_SIZE;index++){ 139 maxev = ( maxev > __complex_abs(EV(index)) )?maxev:(__complex_abs(EV(index))); 140 } 141 std::cout<< "Jacobi迭代矩陣的譜半徑為："<< maxev<<std::endl; 142 143 //譜半徑大於1 迭代法則發散 144 if(maxev >= 1){ 145 std::cout<<"Jacobi迭代算法不收斂！"<<std::endl; 146 flag =1; 147 } 148 149 //迭代法收斂則進行迭代的計算 150 if (flag == 0 ){ 151 std::cout<<"Jacobi迭代算法譜半徑小於1,該算法收斂"<<std::endl; 152 std::cout<<"Jacobi迭代法迭代次數和精度： "<< std::endl << iteration_num<<" "<<accuracy<<std::endl; 153 154 //迭代計算 155 for( k = 0 ;k < iteration_num ; k++ ){ 156 for(i = 0;i< MATRIX_SIZE_ ; i++){ 157 x_k1(i) = x_k(i) + ( b(i) - Jacobi_sum(A,x_k,i) )/A(i,i); 158 temp = fabs( x_k1(i) - x_k(i) ); 159 if( fabs( x_k1(i) - x_k(i) ) > R ) 160 R = temp; 161 } 162 163 //判斷進度是否達到精度要求 達到進度要求后 自動退出 164 if( R < accuracy ){ 165 std::cout <<"Jacobi迭代算法迭代"<< k << "次達到精度要求."<< std::endl; 166 isFlag = 1; 167 break; 168 } 169 170 //清零R，交換迭代解 171 R = 0; 172 x_k = x_k1; 173 } 174 if( !isFlag ) 175 std::cout << std::endl <<"迭代"<<iteration_num<<"次后仍然未達到精度要求，若不滿意該解，請再次運行加大循環次數！"<< std::endl; 176 return x_k1; 177 } 178 //否則返回0 179 return x_k; 180 } 181 182 //Jacobi迭代法的一步求和計算 183 double Jacobi_sum(Mat_A &A,Mat_B &x_k,int i) { 184 double sum; 185 for(int j = 0; j< MATRIX_SIZE_;j++){ 186 sum += A(i,j)*x_k(j); 187 } 188 return sum; 189 }
歡迎大家關注我的微信公眾號「佛系師兄」，里面有關於 Ceres 以及 OpenCV 庫等一些技術文章。
比如
「反復研究好幾遍，我才發現關於 CMake 變量還可以這樣理解！」
更多好的文章會優先在里面不定期分享！打開微信客戶端，掃描下方二維碼即可關注！
