Исключение Гаусса для матрицы NxM - PullRequest
Новая
       4

Исключение Гаусса для матрицы NxM

1 голос
/ 05 марта 2011 
/* 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. Может ли кто-нибудь помочь мне сделать это? Я не очень хорош в математике. Я получил этот код на каком-то сайте, и я застрял на нем.

Ответы [ 3 ]

17 голосов
/ 10 марта 2011
  1. (необязательно). Поймите это .Сделайте несколько примеров на бумаге.
  2. Не пишите код для исключения Гаусса самостоятельно .Без некоторой заботы, наивный поворот Гаусса является неустойчивым.Вы должны масштабировать линии и позаботиться о повороте с наибольшим элементом, отправная точка которого там .Обратите внимание, что этот совет применим для большинства алгоритмов линейной алгебры.
  3. Если вы хотите решить системы уравнений, Разложение LU , Разложение QR (стабильнее, чем LU, но медленнее), Разложение Холецкого (в случае, когда система симметрична) или SVD (в случае, когда система не квадратная) почти всегда являются лучшим выбором.Однако гауссовское исключение лучше всего подходит для вычисления определителей.
  4. Используйте алгоритмы из LAPACK для задач, требующих устранения Гаусса (например, решение систем или вычисление определителей).В самом деле.Не катай свои собственные.Так как вы занимаетесь C ++, вас может заинтересовать Armadillo , который позаботится о многих вещах для вас.
  5. Если вам приходится кататься самостоятельно по педагогическим причинам, сначала посмотрите на Числовые рецепты, версия 3 .Версию 2 можно найти в Интернете бесплатно, если у вас мало средств или у вас нет доступа к библиотеке.
  6. В качестве общего совета не пишите алгоритмы, которые вы не понимаете.
3 голосов
/ 10 марта 2011

Вы просто не можете применить исключение Гаусса непосредственно к проблеме NxM. Если у вас больше уравнений, чем неизвестных, ваша проблема переопределена, и у вас нет решения, что означает, что вам нужно использовать что-то вроде метода наименьших квадратов. Скажем, что у вас есть A * x = b, тогда вместо x = inv (A) * b (когда N = M) вы должны выполнить x = inv (A ^ T * A) * A ^ T * b .

В случае, когда у вас меньше уравнений, чем неизвестных, тогда ваша проблема недоопределена, и у вас есть бесконечное количество решений. В этом случае вы либо выбираете одно случайным образом (например, устанавливаете для некоторых из неизвестных произвольное значение), либо вам необходимо использовать регуляризацию, что означает попытку добавления некоторых дополнительных ограничений.

0 голосов
/ 14 июня 2017

Вы можете применить уменьшение эшелона, как в этом фрагменте 

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