Матрицы

Для нахождения обратной матрицы нужно найти `n^2` чисел `x_{ij}`, таких что
`[(a_{11},...,a_{1n}),(vdots,ddots,vdots),(a_{n1},...,a_{n n})] [(x_{11},...,x_{1n}),(vdots,ddots,vdots),(x_{n1},...,x_{n n})]=[(1,0,...,0,0),(vdots,,ddots,,vdots),(0,0,...,0,1)]`

Мы можем найти неизвестные числа, решая `n` систем линейных уравнений с одной и той же матрицей коэффициентов `A`, у которых вектор неизвестных `x^j` представляет собой `j`-ый столбец обратной матрицы, а вектор свободных членов `e^j` — `j`-ый столбец единичной матрицы:
`A*x^j=e^j`

Для вычислений можно использовать преобразование матрицы `A` в форму, именуемую `LU`-разложением, или `LU`-декомпозицией матрицы коэффициентов.

Данная форма может быть получена с помощью алгоритма исключения Гаусса за `Theta(n^3)`.

Верхнетреугольная матрица `U` представляет собой результат исключения, а нижнетреугольную матрицу `L` образована единицами на главной диагонали и множителями, вычисляемыми в процессе исключения Гаусса для обнуления соответствующих коэффициентов в строках матрицы.

Произведение `LU` этих матриц равно исходной матрице `А`.

Следовательно, решение системы линейных уравнений `Ax=b` эквивалентно решению системы `LUx = b`. Решить ее можно следующим образом. Обозначим `y = Ux`, тогда `Ly = b`.

Сначала решим систему `Ly =b`, что очень просто сделать, поскольку `L` — нижнетреугольная матрица. Затем решим систему `Ux = y`, что опять же несложно, поскольку `U` — верхнетреугольная матрица.

Таким образом, с помощью `LU`-разложения матрицы `A`, мы можем решать системы линейных уравнений `Ax = b` для разных векторов свободных членов `b` за `Theta(n^2)`.

Так как каждое решение можно найти за `Theta(n^2)`, то общее время для нахождения обратной матрицы равно `Theta(n^3)`.

bool LUdecomp(vector<vector<double>>&A, vector<int> &P /*, bool& inv */) {
  int n=A.size();
//bool inv=0; для определителя
  P.resize(n);
  for(int i=0;i<n;++i) P[i]=i;
  for(int i=0; i<n-1; ++i)
  { // находим маскимальное значение в столбце для уменьшения ошибки вычислений
    int imax=i;
    for(int j=i+1; j<n; ++j)
      if(fabs(A[j][i])>fabs(A[imax][i])) imax=j;
    if(fabs(A[imax][i])<1e-12) return 0;
    if(imax!=i)
    { swap(A[imax],A[i]);
      swap(P[imax],P[i]);
//    inv=!inv;
    }
    for(int j=i+1; j<n; ++j)
    {
      double t=A[j][i]/A[i][i];
      for(int k=i; k<n; ++k)
        A[j][k]-=A[i][k]*t;
      A[j][i]=t;
    }
  }
  return 1;
}
// решение системы LUx=b
vector<double> LUsolve(const vector<vector<double>>&LU, 
    const vector<int> &P, const vector<double>& b) {
  int n=LU.size();
  vector<double> x(n),y(n);
// Ly=b
  for (int i = 0; i < n; i++) {
        y[i] = b[P[i]];
        for (int j = 0; j < i; j++)
            y[i] -= LU[i][j] * y[j];
    }
// Ux=y
    for (int i = n - 1; i >= 0; i--) {
        for (int j = i + 1; j < n; j++)
            y[i] -= LU[i][j] * x[j];
        x[i] = y[i]/LU[i][i];
    }
  return x;
}
bool LUinverse(vector<vector<double>> A, vector<vector<double>>&IA)
{ int n=A.size();
  vector<int> P;
  if(! LUdecomp(A,P)) return 0;
  IA.assign(n,vector<double>(n));
  for(int i=0;i<n;++i)
  { vector<double> b(n,0.0);
    b[i]=1;
    auto x=LUsolve(A,P,b);
    for(int j=0;j<n;++j)
      IA[j][i]=x[j];
  }
  return 1;
}

Метод исключения Гаусса помогает и при вычислении определителя.

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

// определитель
double LUdeterm(const vector<vector<double>>&LU, bool inv) {
    double det =inv?-1:1;
    for (int i = 0; i < LU.size(); i++)
        det *= LU[i][i];
    return det;
}