c++ - 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;
}

上面的高斯消除算法在 NxN 矩阵上运行良好。但我需要它来处理 NxM 矩阵。任何人都可以帮我做吗？我数学不是很好。我在某个网站上得到了这个代码，我被困在它上面。

Answer 1

1. （可选）了解这一点。在纸上做一些例子。
2. 不要自己编写高斯消元代码。如果不小心，天真的高斯旋转是不稳定的。您必须缩放线条并注意使用最大元素进行旋转，起点就在那里。请注意，此建议适用于大多数线性代数算法。
3. 如果要求解方程组，LU 分解、QR 分解（比 LU 更稳定，但更慢）、Cholesky 分解（在系统对称的情况下）或SVD（在系统不是正方形的情况下）几乎总是更好选择。然而，高斯消元法最适合计算行列式。
4. 使用LAPACK的算法来解决需要高斯消除的问题（例如求解系统或计算行列式）。真的。不要自己滚。由于您正在使用 C++，您可能会对Armadillo感兴趣，它可以为您处理很多事情。
5. 如果您出于教学原因必须自己动手，请先查看Numerical Recipes, version 3。如果您预算不足/无法访问图书馆，则可以在网上免费找到第 2 版。
6. 作为一般建议，不要编写您不了解的算法。

Answer 2

您不能将高斯消除直接应用于 NxM 问题。如果您的方程多于未知数，则您的问题是超定的并且您没有解决方案，这意味着您需要使用最小二乘法之类的方法。假设你有 A*x = b，那么你必须做 x = inv(A^T*A)*A^T*b 而不是 x = inv(A)*b（当 N=M 时） .

如果你的方程比未知数少，那么你的问题是不确定的，你有无限的解决方案。在这种情况下，您要么随机选择一个（例如将一些未知数设置为任意值），要么您需要使用正则化，这意味着尝试添加一些额外的约束。

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