Начало
Учебные материалы
Матрицы
Представление матриц Операции над матрицами LU-разложение и вычисление обратной матрицы Умножение матриц Алгоритм умножения матриц Штрассена
LU-разложение и вычисление обратной матрицы
`[(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`-декомпозицией матрицы коэффициентов.
Данная форма может быть получена с помощью [алгоритма исключения Гаусса](51057-3.html) за `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)`.
--
```c++
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; ++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;
}
```
--
Метод исключения Гаусса помогает и при вычислении определителя.
Определитель верхнетреугольной матрицы равен произведению ее элементов на главной диагонали,
а элементарные операции, выполняемые алгоритмом исключения Гаусса либо оставляют его значение неизменным,
либо меняют знак, либо приводят к умножению его на константу. В результате мы можем вычислить
определитель матрицы за кубическое время.
```c++
// определитель
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;
}
```
