C++ NxM 矩阵的高斯消元法

Gauss Elimination for NxM matrix

提问by user645332

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

The above gaussian eleimination algorithm works fine on NxN matrices. But I need it to work on NxM matrix. Can anyone help me to do it? I am not very good at maths. I got this code on some website and i am stuck at it.

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

回答by Alexandre C.

1. (optional) Understand this. Do some examples on paper.
2. Don't write code for Gaussian elimination yourself. Without some care, the naive gauss pivoting is unstable. You have to scale the lines and take care of pivoting with the greatest element, a starting point is there. Note that this advice holds for most linear algebra algorithms.
3. If you want to solve systems of equations, LU decomposition, QR decomposition(stabler than LU, but slower), Cholesky decomposition(in the case the system is symmetric) or SVD(in the case the system is not square) are almost always better choices. Gaussian elimination is best for computing determinants however.
4. Use the algorithms from LAPACKfor the problems which need Gaussian elimination (eg. solving systems, or computing determinants). Really. Don't roll your own. Since you are doing C++, you may be interested in Armadillowhich takes care of a lot of things for you.
5. If you must roll your own for pedagogical reasons, have a look first at Numerical Recipes, version 3. Version 2 can be found online for free if you're low on budget / have no access to a library.
6. As a general advice, don't code algorithms you don't understand.

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

回答by Jean-Marc Valin

You just cannot apply Gaussian elimination directly to an NxM problem. If you have more equations than unknowns, the your problem is over-determined and you have no solution, which means you need to use something like the least squares method. Say that you have A*x = b, then instead of having x = inv(A)*b (when N=M), then you have to do x = inv(A^T*A)*A^T*b.

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

In the case where you have less equations then unknowns, then your problem is underdetermined and you have an infinity of solutions. In that case, you either pick one at random (e.g. setting some of the unknowns to an arbitrary value), or you need to use regularization, which means trying adding some extra constraints.

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

回答by Gregor Bereznicki

You can apply echelon reduction, like in this snippet

您可以应用梯形减少，就像在这个片段中一样

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