提綱
- 背景介紹
- 三角方程組
- Gauss消去法
- 附錄
一、背景介紹
1.1 線性方程組的相關概念
線性方程組在解決現實實際問題中直接產生，最小二乘資料擬合、微分方程邊值問題和初邊值問題的數值解產生了大量的線性方程組。
線性方程組係數矩陣的型別分別有
- 稠密型(dense)：幾乎所有元素都是非零的
- 稀疏型(sparse)：有大量零元素
- 帶狀的(banded)
- 三角狀(triangular)
- 塊狀的(block structure)
解線性方程組的方法可以分為兩類
- 直接法(direct method)
經過有限步四則運算可球的方程組準確解的方法 - 迭代法(iterative method)
從一個近似值出發，構造某種演算法，使其逐步接近準確解
大多科學計算應用經過建模和數值離散後，都可歸結為如下兩種形式方程組的求解：
方程組形式
矩陣形式
\(Ax=b\)有唯一解\(\iff A\)非奇異
C++中的線性方程組
線上性代數中，一矩陣的尺寸通常稱為階數(order)或維度(dimension)。以下示例程式碼在主函式中定義了稀疏矩陣\(A\)，常向量\(b\)和解向量\(x\)。
在Eigen庫中，可以採用Eigen::MatrixXd表示矩陣型別，採用Eigen::VectorXd表示向量型別。矩陣和向量的尺寸可以在建立時進行設定。
需要注意的是，
Eigen庫中Eigen::VectorXd預設為列向量，如果需要將其作為行向量進行運算，需要在使用時進行轉置，例如：X.transpose()
即使沒有硬性的要求，但還是建議讀者使用
const size_t型別的變數單獨儲存矩陣的尺寸，這將使得程式碼維護變得更容易。
#include <iostream>
#include <Eigen/Dense>
int main() {
// 矩陣的階數
const size_t order = 6;
// 定義係數矩陣 A
Eigen::MatrixXd A(order, order); // 指定尺寸為 order * order
// 定義常向量 b
Eigen::VectorXd b(order); // 指定尺寸為 order * 1
// 定義解向量 x
Eigen::VectorXd x(order); // 指定尺寸為 order * 1
}
採用直接法求解線性方程組的求解器通常包含三個輸入，即：係數矩陣\(A\)、常向量\(b\)和解向量\(x\)。
在進行求解前，首先應當檢查輸入是否符合求解器要求，例如針對上三角矩陣的求解器需要檢查係數矩陣是否為上三角矩陣；一般的，輸入應滿足三個要求：
- 係數矩陣\(A\)為方陣
- 係數矩陣\(A\)的行數等於常向量\(b\)的行數
- 係數矩陣\(A\)的列數等於解向量\(x\)的行數
矩陣的行數可以透過
.rows()方法得到
矩陣的列數可以透過.cols()方法得到
該方法對於向量同樣適用，特別的，
Eigen庫中向量的列數總是1
以下給出參考的實現：
void size_check(const Eigen::MatrixXd& A,
const Eigen::VectorXd& b,
const Eigen::VectorXd& x)
{
// 檢查A是否為方陣
if (A.rows() != A.cols()) {
throw std::invalid_argument("Error: The coefficient matrix of the system of equations is not a square matrix.");
}
// 檢查係數矩陣A的尺寸是否與常向量b的尺寸匹配
if (A.rows() != b.rows()) {
throw std::invalid_argument("Error: The order of the coefficient matrix A does not match the order of the constant vector b.");
}
// 檢查係數矩陣A的尺寸是否與解向量x的尺寸匹配
if (A.cols() != x.rows()) {
throw std::invalid_argument("Error: The order of the coefficient matrix A does not match the order of the solution vector x.");
}
}
void solve(const Eigen::MatrixXd& A,
const Eigen::VectorXd& b,
Eigen::VectorXd& x)
{
// 檢查尺寸是否合適
size_check(A, b, x);
// 求解
// ...
}
在實際實現時有幾個應注意的細節
為什麼不將解向量\(x\)作為輸出？
將解向量\(x\)作為輸出的函式的使用方式為：ans=solve(A,b)，如果返回值的尺寸與變數ans的尺寸不一致則會導致程式錯誤。為了避免該問題，必須在建立變數ans時設定尺寸，並在求解前檢查尺寸，虛擬碼如下
Eigen::VetorXd x(order);
if (x.rows() == A.cols()) { // 尺寸檢查
x = solve(A,b);
}
顯然，形如ans=solve(A,b,x)的求解器更為易用，其型別檢查可以在函式內部完成，這帶來了更好的封裝性、可維護性。
在必要的時候新增&和const關鍵字
在傳遞函式引數時，&關鍵字表明瞭該傳參方式為引用傳參，區別於普通傳參，引用傳參方式使得函式無需在其內部複製一個副本，而是可以直接在原變數上進行操作。無需複製副本顯著降低了程式的效能開銷。
對於普通傳參，const關鍵字表明內部複製的副本為常變數。對於引用傳參，const關鍵字表明該函式不具有修改該變數的許可權，只具備讀取（訪問）的許可權。
三角方程組
下三角方程組
解法：前代法(Forward substitution)
下三角矩陣判斷
Eigen庫並沒有提供直接的判斷矩陣是否為下三角矩陣的方法，因此採用瞭如下的判斷方法：
- 首先提取矩陣的嚴格上三角部分（不包含對角線）
- 判斷其是否全部為零，如果嚴格上三角部分全部為零，那麼其為下三角矩陣
前代法求解
- 檢查輸入尺寸是否匹配
- 判斷係數矩陣是否為下三角矩陣
- 採用前代法求解。
外層迴圈用於遍歷解向量\(x\)的每個元素，從下標0開始，遍歷至下標n-1結束。迴圈內部分佈實現式\((2.1)\)的計算，對於求和部分，巢狀內層迴圈實現。
矩陣/向量元素訪問
在訪問矩陣/向量的元素時元素，採用括號運算子進行訪問。
#include "check.h"
bool isLowerTriangular(const Eigen::MatrixXd& A) {
// 獲取矩陣的嚴格上三角部分（不包括對角線）
Eigen::MatrixXd upperTriangularPart = A.triangularView<Eigen::StrictlyUpper>();
// 檢查嚴格上三角部分是否全為零
return upperTriangularPart.isZero();
}
void forward_substitution(const Eigen::MatrixXd& A,
const Eigen::VectorXd& b,
Eigen::VectorXd& x)
{
// 檢查尺寸是否匹配
size_check(A, b, x);
// 判斷係數矩陣是否為下三角矩陣
if (!isLowerTriangular(A)) {
throw std::invalid_argument("Error: The matrix is not lower triangular.");
}
for (size_t i = 0; i < A.rows(); ++i) {
x(i) = b(i);
for (size_t j = 0; j + 1 <= i; ++j) { // j < i - 1
x(i) -= A(i, j) * x(j);
}
x(i) /= A(i, i);
}
}
注意事項
應當注意C++中的陣列索引是從
0開始的，Eigen庫也沿用了這一習慣。
在求和\(\sum_{j=1}^{i-1}a_{ij}x_j\)的實現中，很容易錯誤的使用
j<=i-1作為迴圈的終止條件，這實際上有一個風險，當i=0的時候，i-1並不是-1，而是最大的size_t型別的數，這將導致終止條件錯誤，因此，應當用j+1<=i
上三角方程組
解法：回代法(Back substitution)
上三角矩陣判斷
Eigen庫並沒有提供直接的判斷矩陣是否為上三角矩陣的方法，因此採用瞭如下的判斷方法：
- 首先提取矩陣的嚴格下三角部分（不包含對角線）
- 判斷其是否全部為零，如果嚴格下三角部分全部為零，那麼其為上三角矩陣
回代法求解
- 檢查輸入尺寸是否匹配
- 判斷係數矩陣是否為上三角矩陣
- 採用回代法求解。
外層迴圈用於遍歷解向量\(x\)的每個元素，從下標n-1開始，遍歷至下標0結束。迴圈內部分佈實現式\((2.2)\)的計算，對於求和部分，巢狀內層迴圈實現。
bool isUpperTriangular(const Eigen::MatrixXd& A) {
// 獲取矩陣的嚴格下三角部分（不包括對角線）
Eigen::MatrixXd lowerTriangularPart = A.triangularView<Eigen::StrictlyLower>();
// 檢查嚴格下三角部分是否全為零
return lowerTriangularPart.isZero();
}
void back_substitution(const Eigen::MatrixXd& A,
const Eigen::VectorXd& b,
Eigen::VectorXd& x)
{
// 檢查尺寸是否匹配
size_check(A, b, x);
// 判斷係數矩陣是否為上三角矩陣
if (!isUpperTriangular(A)) {
throw std::invalid_argument("Error: The matrix is not upper triangular.");
}
size_t n = A.rows();
for (size_t i = n - 1; i != size_t(-1); --i) { // i != -1
x(i) = b(i);
for (size_t j = i + 1; j <= n - 1; ++j) {
x(i) -= A(i, j) * x(j);
}
x(i) /= A(i, i);
}
}
注意事項
外層迴圈的遍歷是從下標
n-1開始，遍歷至下標0結束；一般習慣性的寫法是，以i>=0作為截止條件，但應當注意，size_t型別是非負的，事實上，對於size_t型別的變數，當其值為0時再做-1，其值為size_t(-1)，因此，可以採用i!=size_t(-1)作為截止條件
高斯消元法
一般高斯消元法
高斯消元法（Gaussian Elimination）是一種用於求解線性方程組的經典方法。它透過逐步消去未知數，將方程組化為上三角形式，然後透過回代法求解未知數。高斯消元法主要分為兩個步驟：前向消元和後向回代，本文中將以前向消元為例展開討論。
前向消元（Forward Elimination）
前向消元法是從第一列開始，透過一些列的行變換，逐漸將原矩陣變換為一個上三角矩陣。假定矩陣的尺寸為\(N*N\)，那麼高斯消元法需要進行\(N-1\)次，在第\(i\)時執行如下操作：
- 選擇主元：選擇第\(i\)列的元素\(A_{i,i}\)作為主元
- 消去操作：透過將第\(i\)行的適當倍數加到其他行，使得當前列的其它元素變為零。
消去操作的公式如下：
矩陣的第一步消元過程可以參考以下公式：
在下述程式中，取樣行向量相減的方式實現高斯消元法，相較於逐個元素相減，程式碼更簡潔易懂，易維護。
void simple_gauss_elimination(Eigen::MatrixXd& A, Eigen::VectorXd& b) {
// 檢查尺寸是否匹配
size_check(A, b);
size_t n = A.rows();
// 逐步消元為上三角矩陣
for (size_t k = 0; k < n - 1; ++k) {
// 提取矩陣的第k行
Eigen::VectorXd temp = A.row(k);
// 將第i列索引大於i的元素消為0
for (size_t i = k + 1; i < n; ++i) {
// 計算比值
double m = A(i, k) / A(k, k);
// 消元
A.row(i) -= m * temp;
b(i) -= m * b(k);
}
}
}
改進的高斯消元法
若\(a^{(k)}_{kk}\to 0\)，則\(m=a_{ik}^{(k)}/a_{kk}^{(0)}\to\infty\)，此時直接用高斯消元法求解線性方程組是會由於舍入誤差的擴大，而導致解失真。
因此在原高斯消元法的基礎上，可以做改進，新增主元的選擇過程，該方法稱為列主元法，具體流程如下：
- 尋找第\(k\)列中第\(k\)行到第\(n\)行最大的元素，記為\(a_{jk}\)
- 將第\(j\)行與第\(k\)行交換
- 進行高斯消元法
void gauss_elimination(Eigen::MatrixXd& A, Eigen::VectorXd& b) {
// 檢查尺寸是否匹配
size_check(A, b);
size_t n = A.rows();
// 逐步消元為上三角矩陣
for (size_t k = 0; k < n; ++k) {
// 選擇主元
size_t j = k;
double max = abs(A(j, k));
for (size_t i = k + 1; i < n; ++i) {
double d = abs(A(i, k));
if (d > max) { // 選擇絕對值最大的元素
j = i; max = d;
}
}
// 交換主元
if (j != k) {
Eigen::VectorXd temp = A.row(j);
A.row(j) = A.row(k);
A.row(k) = temp;
double temp_b = b(j);
b(j) = b(k);
b(k) = temp_b;
}
// 將第i列索引大於i的元素消為0
for (size_t i = k + 1; i < n; ++i) {
// 計算比值
double m = A(i, k) / A(k, k);
// 消元
A.row(i) -= m * A.row(k);
b(i) -= m * b(k);
}
}
}
注意事項
對方程\(Ax=b\)的係數矩陣\(A\)和常向量\(b\)同時做行變換時，方程的解\(x\)不變。
基於高斯消元法的一般線性方程求解
對於一般的線性方程組，可以先用高斯消元法將係數矩陣轉化為上三角矩陣，再透過回代法求解。
void gauss_solve(Eigen::MatrixXd A,
Eigen::VectorXd b,
Eigen::VectorXd& x)
{
// 檢查尺寸是否匹配
size_check(A, b, x);
// 高斯消元法轉為上三角矩陣
gauss_elimination(A, b);
// 透過回代法求解
back_substitution(A, b, x);
}
注意事項
切忌捨本逐末，雖然新增引用修飾符可以一定程度上提升效能，但是這會導致稀疏矩陣\(A\)和常向量\(b\)被修改，而使用者往往容易忽略這一點，因此為了保證安全性，此處不使用引用傳參。
截止到目前，對係數矩陣\(A\)為下三角形矩陣的線性方程組有兩種求解方法，一種是採用前代法，一種是採用高斯消元結合回代法，在附錄中我們對同一組資料採用兩種方法分別計算結果，進行交叉驗證。
附錄
功能測試方法
構建函式（方法）的測試程式流程如下：
- 從函式（方法）的名稱中提取縮寫，作為名聲空間的字首
- 定義測試函式，命名為
test()，如果需要可以設計多個，例如：test1(),test2() - 實現測試函式，一般來說，有以下步驟：①生成資料，②呼叫方法，③列印資料以及結果
- 在主函式中，呼叫該名聲空間下的測試函式
test()，一般需要使用try-catch結構
示例程式碼如下：
namespace SMP{
void test() {
std::cout << "Hello World!";
}
}
int main() {
try{
SMP::test();
}
catch (const std::exception& e) {
std::cerr << "Error: " << e.what() << std::endl;
}
}
在後續的附錄內容中，將省略main函式的設計，讀者只需按照上述方法呼叫即可。
前代法測試
namespace FWD{
// test for forward_substitution()
void test() { // 矩陣的階數
const size_t order = 5;
// 定義係數矩陣 A
Eigen::MatrixXd A(order, order);
// 定義常向量 b
Eigen::VectorXd b(order);
// 定義解向量 x
Eigen::VectorXd x(order);
// 設定矩陣為隨機數
A.setRandom();
b.setRandom();
// 處理為方便手算的數字
A = (1.5 + A.array()) * 2;;
b *= 10;
A = A.array().round().matrix();
b = b.array().round().matrix();
// 將嚴格上三角部分設定為零，使其成為下三角矩陣
A.triangularView<Eigen::StrictlyUpper>().setZero();
// 前代法
forward_substitution(A, b, x);
// 輸出結果
std::cout << "A=\n" << A << "\n";
std::cout << "b=\n" << b << "\n";
std::cout << "x=\n" << x << "\n";
}
}
效果展示
程式的輸出如下圖所示（經過拼接），經檢驗，該計算結果正確（讀者感興趣的可以手算一下試試）。
回代法測試
namespace BCK{
// test for back_substitution()
void test() { // 矩陣的階數
const size_t order = 5;
// 定義係數矩陣 A
Eigen::MatrixXd A(order, order);
// 定義常向量 b
Eigen::VectorXd b(order);
// 定義解向量 x
Eigen::VectorXd x(order);
// 設定矩陣為隨機數
A.setRandom();
b.setRandom();
// 處理為方便手算的數字
A = (1.5 + A.array()) * 2;;
b *= 10;
A = A.array().round().matrix();
b = b.array().round().matrix();
// 將嚴格下三角部分設定為零，使其成為上三角矩陣
A.triangularView<Eigen::StrictlyLower>().setZero();
// 回代法
back_substitution(A, b, x);
// 輸出結果
std::cout << "A=\n" << A << "\n";
std::cout << "b=\n" << b << "\n";
std::cout << "x=\n" << x << "\n";
}
}
效果展示
程式的輸出如下圖所示（經過拼接），經檢驗，該計算結果正確.
一般高斯消元法測試
namespace S_GSE {
// test for simple_gauss_elimination
void test() { // 矩陣的階數
const size_t order = 5;
// 定義係數矩陣 A
Eigen::MatrixXd A(order, order);
// 定義常向量 b
Eigen::VectorXd b(order);
// 設定矩陣為隨機數
A.setRandom();
b.setRandom();
// 調整顯示精度為小數點後兩位
std::cout << std::fixed << std::setprecision(2);
// 輸出消元前矩陣
std::cout << "A=\n" << A << "\n";
std::cout << "b=\n" << b << "\n";
// 前代法
simple_gauss_elimination(A, b);
// 輸出消元后矩陣
std::cout << "A=\n" << A << "\n";
std::cout << "b=\n" << b << "\n";
}
}
效果展示
程式的輸出如下圖所示（經過拼接），顯示精度為小數點後兩位；經檢驗，該計算結果正確.
列主元法改進的高斯消元法測試
namespace GSE {
// test for simple_gauss_elimination
void test() { // 矩陣的階數
const size_t order = 5;
// 定義係數矩陣 A
Eigen::MatrixXd A(order, order);
// 定義常向量 b
Eigen::VectorXd b(order);
// 設定矩陣為隨機數
A.setRandom();
b.setRandom();
// 調整顯示精度為小數點後兩位
std::cout << std::fixed << std::setprecision(2);
// 輸出消元前矩陣
std::cout << "A=\n" << A << "\n";
std::cout << "b=\n" << b << "\n";
// 前代法
gauss_elimination(A, b);
// 輸出消元后矩陣
std::cout << "A=\n" << A << "\n";
std::cout << "b=\n" << b << "\n";
}
}
程式的輸出如下圖所示（經過拼接），顯示精度為小數點後兩位；經檢驗，該計算結果正確.
高斯+回代法求解
namespace GS_SOLVE{
void test1() {
const size_t order = 5;
Eigen::MatrixXd A(order, order);
Eigen::VectorXd b(order);
Eigen::VectorXd x(order);
// 設定矩陣為隨機數
A.setRandom();
b.setRandom(); b = (1.0 + b.array()) * 5;
// 前代法
gauss_solve(A, b, x);
// 輸出結果
std::cout << std::fixed << std::setprecision(2);
std::cout << "A=\n" << A << "\n";
std::cout << "b=\n" << b << "\n";
std::cout << "x=\n" << x << "\n";
}
void test2() {
const size_t order = 5;
Eigen::MatrixXd A(order, order);
Eigen::VectorXd b(order);
Eigen::VectorXd x1(order);
Eigen::VectorXd x2(order);
// 設定矩陣為隨機數
A.setRandom();
b.setRandom(); b = (1.0 + b.array()) * 5;
// 將上三角部分設定為零，使其成為下三角矩陣
A.triangularView<Eigen::StrictlyUpper>().setZero();
// 高斯
gauss_solve(A, b, x1);
// 前代法
forward_substitution(A, b, x2);
// 輸出結果
std::cout << std::fixed << std::setprecision(2);
std::cout << "GS_solve:\n" << "x1=\n" << x1 << "\n";
std::cout << "back_stt:\n" << "x2=\n" << x2 << "\n";
}
}
測試1
函式GS_SOLVE::test1()用於測試高斯求解是否能夠正常工作，該程式的輸出如下圖所示（經過拼接），顯示精度為小數點後兩位；經檢驗，該計算結果正確.
測試2
函式GS_SOLVE::test2()採用交叉驗證法，分別採用前代法，一種是採用高斯消元結合回代法求解係數矩陣\(A\)為下三角矩陣的線性方程組，並對比計算結果；經檢驗，結果各方面功能正常。
