Gauss Elimination for NxM matrix

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

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.

Answer 1

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

Answer 2

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.

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.

Answer 3

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 LAPACK for 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 Armadillo which 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.