针对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; 
}