針對NxM矩陣的高斯消除

Question

/* Program to demonstrate gaussian <strong class="highlight">elimination</strong> 
    on a set of linear simultaneous equations 
*/ 

#include <iostream> 
#include <cmath> 
#include <vector> 

using namespace std; 

const double eps = 1.e-15; 

/*Preliminary pivoting strategy 
    Pivoting function 
*/ 
double pivot(vector<vector<double> > &a, vector<double> &b, int i) 
{ 
    int n = a.size(); 
    int j=i; 
    double t=0; 

    for(int k=i; k<n; k+=1) 
    { 
     double aki = fabs(a[k][i]); 
     if(aki>t) 
     { 
      t=aki; 
      j=k; 
     } 
    } 
    if(j>i) 
    { 
     double dummy; 
     for(int L=0; L<n; L+=1) 
     { 
      dummy = a[i][L]; 
      a[i][L]= a[j][L]; 
      a[j][L]= dummy; 
     } 
     double temp = b[j]; 
     b[i]=b[j]; 
     b[j]=temp; 
    } 
    return a[i][i]; 
}   



/* Forward <strong class="highlight">elimination</strong> */ 
void triang(vector<vector<double> > &a, vector<double> &b) 
{ 
    int n = a.size(); 
    for(int i=0; i<n-1; i+=1) 
    { 
     double diag = pivot(a,b,i); 
     if(fabs(diag)<eps) 
     { 
      cout<<"zero det"<<endl; 
      return; 
     } 
     for(int j=i+1; j<n; j+=1) 
     { 
      double mult = a[j][i]/diag; 
      for(int k = i+1; k<n; k+=1) 
      { 
       a[j][k]-=mult*a[i][k]; 
      } 
      b[j]-=mult*b[i]; 
     } 
    } 
} 

/* 
DOT PRODUCT OF TWO VECTORS 
*/ 
double dotProd(vector<double> &u, vector<double> &v, int k1,int k2) 
{ 
    double sum = 0; 
    for(int i = k1; i <= k2; i += 1) 
    { 
    sum += u[i] * v[i]; 
    } 
    return sum; 
} 

/* 
BACK SUBSTITUTION STEP 
*/ 
void backSubst(vector<vector<double> > &a, vector<double> &b, vector<double> &x) 
{ 
    int n = a.size(); 
    for(int i = n-1; i >= 0; i -= 1) 
    { 
    x[i] = (b[i] - dotProd(a[i], x, i + 1, n-1))/ a[i][i]; 
    } 
} 
/* 
REFINED GAUSSIAN <strong class="highlight">ELIMINATION</strong> PROCEDURE 
*/ 
void gauss(vector<vector<double> > &a, vector<double> &b, vector<double> &x) 
{ 
    triang(a, b); 
    backSubst(a, b, x); 
}        

// EXAMPLE MAIN PROGRAM 
int main() 
{ 
    int n; 
    cin >> n; 
    vector<vector<double> > a; 
    vector<double> x; 
    vector<double> b; 
    for (int i = 0; i < n; i++) { 
     vector<double> temp; 
     for (int j = 0; j < n; j++) { 
      int no; 
      cin >> no; 
      temp.push_back(no); 
     } 
     a.push_back(temp); 
     b.push_back(0); 
     x.push_back(0); 
    } 
    /* 
    for (int i = 0; i < n; i++) { 
     int no; 
     cin >> no; 
     b.push_back(no); 
     x.push_back(0); 
    } 
    */ 

    gauss(a, b, x); 
    for (size_t i = 0; i < x.size(); i++) { 
     cout << x[i] << endl; 
    } 
    return 0; 
} 

上述高斯消除算法在N×N矩陣上正常工作。但我需要它在NxM矩陣上工作。任何人都可以幫我做到嗎？我不擅長數學。我在一些網站上得到了這個代碼，我堅持了。

Answer 1

1. （可選）瞭解this。在紙上做一些例子。
2. 不要編寫自己的高斯消除代碼。沒有多少照顧，天真的高斯擺動是不穩定的。您必須縮放線條並採用最大的元素進行旋轉，起點爲​​。請注意，這個建議適用於大多數線性代數算法。
3. 如果你想解決方程組，LU decomposition，QR decomposition（穩定性優於LU，但速度較慢），Cholesky decomposition（在情況下，系統是對稱的）或SVD（在系統中是不是正方形的情況下）幾乎總是更好選擇。然而，高斯消除對於計算行列式來說是最好的。
4. 將LAPACK中的算法用於需要高斯消除（例如求解系統或計算決定因素）的問題。真。不要推出自己的。既然你在做C++，你可能會對Armadillo感興趣，它會爲你處理很多事情。
5. 如果您出於教學原因必須推出自己的產品，請首先看Numerical Recipes，版本3。如果您的預算不足/無法使用圖書館，則可免費在線查找版本2。
6. 作爲一般建議，不要編碼算法，你不明白。

Answer 2

你不能直接應用高斯消除NxM問題。如果你有更多的方程而不是未知數，那麼你的問題就超出了決定性，而你沒有解決方案，這意味着你需要使用像最小二乘法這樣的方法。假設你有A * x = b，那麼你不得不做x = inv（A）* b（當N = M時），那麼你必須做x = inv（A^T * A）* A^T * b 。

如果您有較少的方程，那麼未知，那麼你的問題是欠定的，你有一個無限的解決方案。在這種情況下，您可以隨機選擇一個（例如，將某些未知數設置爲任意值），或者您需要使用正則化，這意味着嘗試添加一些額外的約束。

Answer 3

您可以在這個片段中

#include <iostream> 
#include <algorithm> 
#include <vector> 
#include <iomanip> 

using namespace std; 
/* 
A rectangular matrix is in echelon form(or row echelon form) if it has the following 
    three properties : 
1. All nonzero rows are above any rows of all zeros. 
2. Each leading entry of a row is in a column to the right of the leading entry of 
    the row above it. 
3. All entries in a column below a leading entry are zeros. 

If a matrix in echelon form satisfies the following additional conditions, 
    then it is in reduced echelon form(or reduced row echelon form) : 
4. The leading entry in each nonzero row is 1. 
5. Each leading 1 is the only nonzero entry in its column.        
*/ 

template <typename C> void print(const C& c) { 
    for (const auto& e : c) { 
     cout << setw(10) << right << e; 
    } 

    cout << endl; 
} 
template <typename C> void print2(const C& c) { 
    for (const auto& e : c) { 
     print(e); 
    } 

    cout << endl; 
} 

// input matrix consists of rows, which are vectors of double 
vector<vector<double>> Gauss::Reduce(const vector<vector<double>>& matrix) 
{ 
    if (matrix.size() == 0) 
     throw string("Empty matrix"); 
    auto A{ matrix }; 
    auto mima = minmax_element(A.begin(), A.end(), [](const vector<double>& a, const vector<double>& b) {return a.size() < b.size(); }); 
    auto mi = mima.first - A.begin(), ma = mima.second - A.begin(); 
    if (A[mi].size() != A[ma].size()) 
     throw string("All rows shall have equal length"); 
    size_t height = A.size(); 
    size_t width = A[0].size(); 
    if (width == 0) 
     throw string("Only empty rows"); 

    for (size_t row = 0; row != height; row++) { 
     cout << "processing row " << row << endl; 

     // Search for maximum below current row in column row and move it to current row; skip this step on the last one 
     size_t col{ row }, maxRow{ 0 }; 
     // find pivot for current row (partial pivoting) 
     while (col < width) 
     { 
      maxRow = distance(A.begin(), max_element(A.begin() + row, A.end(), [col](const vector<double>& rowVectorA, const vector<double>& rowVectorB) {return abs(rowVectorA[col]) < abs(rowVectorB[col]); })); 
      if (A[maxRow][col] != 0) // nonzero in this row and column or below found 
       break; 
      ++col; 
     } 
     if (col == width) // e.g. in current row and below all entries are zero 
      break; 
     if (row != maxRow) 
     { 
      swap(A[row], A[maxRow]); 
      cout << "swapped " << row << " and " << maxRow; 
     } 
     cout << " => leading entry in column " << col << endl; 

     print2(A); 
     // here col >= row holds; col is the column of the leading entry e.g. first nonzero column in current row 
     // moreover, all entries to the left and below are zeroed 
     if (row+1 < height) 
      cout << "processing column " << col << endl; 

     // Make in all rows below this one 0 in current column 
     for (size_t rowBelow = row + 1; rowBelow < height; rowBelow++) { 
      // subtract product of current row by factor 
      double factor = A[rowBelow][col]/A[row][col]; 
      cout << "processing row " << rowBelow << " below the current; factor is " << factor << endl; 
      if (factor == 0) 
       continue; 
      for (size_t colRight{ col }; colRight < width; colRight++) 
      { 
       auto d = A[rowBelow][colRight] - factor * A[row][colRight]; 
       A[rowBelow][colRight] = abs(d) < DBL_EPSILON ? 0 : d; 
      } 
      print(A[rowBelow]); 
     } 
    } 
    // the matrix A is in echelon form now 
    cout << "matrix in echelon form" << endl; 
    print2(A); 
    // reduced echelon form follows (backward phase) 
    size_t row(height-1); 
    auto findPivot = [&row, A]() -> size_t { 
     do 
     { 
      auto pos = find_if(A[row].begin(), A[row].end(), [](double d) {return d != 0; }); 
      if (pos != A[row].end()) 
       return pos - A[row].begin(); 
     } while (row-- > 0); 
     return A[0].size(); 
    }; 
    do 
    { 
     auto col = findPivot(); 
     if (col == width) 
      break; 
     cout << "processing row " << row << endl; 
     if (A[row][col] != 1) 
     { 
      //scale row row to make element at [row][col] equal one 
      auto f = 1/A[row][col]; 
      transform(A[row].begin()+col, A[row].end(), A[row].begin()+col, [f](double d) {return d * f; });    
     } 
     auto rowAbove{ row}; 
     while (rowAbove > 0) 
     { 
      rowAbove--; 
      double factor = A[rowAbove][col]; 
      if (abs(factor) > 0) 
      { 
       for (auto colAbove{ 0 }; colAbove < width; colAbove++) 
       { 
        auto d = A[rowAbove][colAbove] - factor * A[row][colAbove]; 
        A[rowAbove][colAbove] = abs(d) < DBL_EPSILON ? 0 : d;     
       } 
       cout << "transformed row " << rowAbove << endl; 
       print(A[rowAbove]); 
      } 
     } 
    } while (row-- > 0); 

    return A; 
}