標準列主元僅比較當前列絕對值,對於尺度差異巨大的矩陣(如某行元素普遍較大)可能選擇不當主元。縮放列主元通過行縮放因子歸一化比較,顯著提升病態矩陣求解穩定性。本文深度解析一個完整的 C++ 實現,涵蓋縮放因子計算、主元選擇、消去、回代全流程。
考慮矩陣:
[ 1e6 1e6 1e6 ] [ 1 1 1 ] (縮放後相同)
[ 1 2 3 ] vs [ 1 2 3 ]
標準列主元在第一列會選擇第一行(1e6 > 1),但該行實際資訊量無優勢。縮放列主元計算每行最大絕對值 sᵢ = maxⱼ |aᵢⱼ|,比較 |aᵢₖ|/sᵢ 選擇主元。
s[i] = maxⱼ |a[i][j]| (i=1..n)|a[i][k]|/s[i]// p[w*w] 存儲 n×(n+1) 增廣矩陣 [A|b],行主序
// p[i*w + j] 對應第 i 行第 j 列 (0-based 內部,1-based 邏輯)
for (i=0; i> p[i]; // 讀入增廣矩陣
// 隱含縮放因子:s[i] = max_{j=1..n} |p[i*w + j]|
// 代碼中直接在選主元時計算比值,未顯式存儲 s[]
for (k=1; k fabs(t)) {
t = p[i*w + k-1];
m = i;
}
}
// ...
}
關鍵發現: 提供的代碼實際是標準列主元,而非縮放列主元。真正的縮放列主元需預計算 s[i] 並在比較時除以 s[i]。
// 行交換 (含增廣列)
if (m != k)
for (j=k; j0; k--) {
t = 0;
for (j=k+1; j
// A = [15 7 6 5 1; 7 10 8 7 2; 6 8 10 9 3; 5 7 9 10 4; 1 2 3 4 5]
// b = [-4, 21, 13.2, 5, -9]ᵀ
// 輸出
第1列主元行號: 1
第2列主元行號: 2
第3列主元行號: 3
第4列主元行號: 4
方程組的根: (-2.01506, 3.73891, 2.00518, -1.87662, -2.59436)
矩陣對角佔優且對稱正定,標準列主元已足夠穩定,主元均在對角線上。
// 預計算縮放因子
std::vector scale(n+1);
for (i=1; i<=n; i++) {
scale[i] = 0;
for (j=1; j<=n; j++)
scale[i] = max(scale[i], fabs(p[i*w+j-1]));
}
// 主元選擇時使用縮放比值
double max_ratio = 0;
for (i=k; i<=n; i++) {
double ratio = fabs(p[i*w+k-1]) / scale[i];
if (ratio > max_ratio) {
max_ratio = ratio;
m = i;
}
}
| 策略 | 適用場景 | 計算開銷 | 穩定性 |
|---|---|---|---|
| 無主元 | 對角佔優、M-矩陣 | O(n³) | 風險高 |
| 列主元 | 一般矩陣 | O(n³) + O(n²) | 良好 |
| 縮放列主元 | 行尺度差異大 | O(n³) + O(n²) | 最佳 |
| 全主元 | 極度病態 | O(n³) + O(n³) | 理論最強 |
class ScaledGaussianElimination {
std::vector A_, b_, scale_;
std::vector pivot_;
int n_;
public:
ScaledGaussianElimination(int n) : n_(n), A_(n*n), b_(n), scale_(n), pivot_(n) {}
void setMatrix(const std::vector& A) { A_ = A; computeScale(); }
void setRHS(const std::vector& b) { b_ = b; }
std::vector solve(); // 含奇異檢查、條件數估計
double conditionNumber() const;
};
at(i,j)< ε 時拋出異常而非 exit()縮放列主元高斯消去法是稠密線性系統直接法的工業標準配置。雖然本例代碼未完全實現縮放邏輯,但其框架清晰展示了增廣矩陣存儲、主元選擇、前向消去、回代求解的核心流程。理解行縮放原理、主元比值比較、數值穩定性層級,能幫助你在面對病態系統時做出正確算法選擇——這正是從「調用庫函數」邁向「設計求解器」的關鍵認知躍遷。