вычислить обратную матрицу

[d.dragon.n76]

Возникла необходимость вычислить обратную матрицу. Вспомнил метод Гаусса и набросал код. Вроде работает без ошибок

Вопросы такие:
1. когда применим этот метод (код) /желательно указать точный критерий типа норма, собственные числа, число обусловленности и отослать к конкретной книге/
2. можно ли быстрее или рациональнее /со ссылкой на книжку/

Цитата:

void inverse_matrix(std::vector<std::vector<double> > & A)
{
int k,i,j;
double koef;
int n = A.size();

std::vector<std::vector<double> > a(n,std::vector<double>(n+n));
for( i=0; i<n; ++i){
for( j=0; j<n; ++j){
a[i][j] = A[i][j];
}// for_j
for( j=n; j<n+n; ++j){
if(i+n==j) a[i][j] = 1.0;
else a[i][j] = 0.0;
}
} // for_i

for( k=0; k<n; ++k ) {
for( i=k+1; i<n; ++i){
koef = -a[i][k]/a[k][k];
for( j=k; j<n+n; ++j ){
a[i][j] += koef*a[k][j];
}// for_j
}// for_i
} // for_k



for(int i=0; i<n; ++i){
koef = a[i][i];
for(int j=0; j<n+n; ++j){
a[i][j] /= koef;
}
}

for(int k=n-1; k>0; --k) {
for(int i=k-1; i>=0; --i) {
koef = a[i][k];
for(int j=0; j<n+n; ++j) {
a[i][j] -= koef*a[k][j];
} // for_j
} // for_i
} // for_k


for( i=0; i<n; ++i){
for( j=0; j<n; ++j){
A[i][j] = a[i][j+n];
} // for_i
} // for_j

}

[Brain13]

Метод Гаусса(как и метод Гаусса-Жордана) применим практически везде, например, его применяют тогда, когда не знают особенностей матрицы(симметричность, положительность, разряженность,...), потому что в общем случае это самый эффективный метод. И в общем случае лучше него нет.
Есть еще такие методы как LU-факторизации(очень эффективен, если у Вас матрица постоянна, а меняется только правая часть), QR методы, метод Холесского(для сверхбольших, разряженных, симметричных матриц), но все они являются модификациями метода Гаусса и дают выигрыш к производительности в определенных случаях.
Есть книжка на англ. с исходниками в нормальном pdf:
NUMERICAL
RECIPES
The Art of Scientific Computing
Third Edition
Поищите, или могу скинуть.

[Cave]

Ваша реализация также не сработает, если в матрице на диагонали есть нули.
Offtop: вижу, что вы относительно заботитесь о понятности кода (комментарии "// for i" - это хорошо). Так вот, называть две различные переменные с матрицами 'a' и 'A' - это не очень хорошо.

[d.dragon.n76]

Brain13 в сообщении #289038 писал(а):

Метод Гаусса(как и метод Гаусса-Жордана) применим практически везде, например, его применяют тогда, когда не знают особенностей матрицы(симметричность, положительность, разряженность,...), потому что в общем случае это самый эффективный метод. И в общем случае лучше него нет.
Есть еще такие методы как LU-факторизации(очень эффективен, если у Вас матрица постоянна, а меняется только правая часть), QR методы, метод Холесского(для сверхбольших, разряженных, симметричных матриц), но все они являются модификациями метода Гаусса и дают выигрыш к производительности в определенных случаях.
Есть книжка на англ. с исходниками в нормальном pdf:
NUMERICAL
RECIPES
The Art of Scientific Computing
Third Edition
Поищите, или могу скинуть.

подскажите страничку с критериями применимости метода Гаусса или Ж.Г. (спрашиваю из-за очень плохого знания английского)

[d.dragon.n76]

Cave в сообщении #289112 писал(а):

Ваша реализация также не сработает, если в матрице на диагонали есть нули.

Проблема возникает и в случае если на диагонали не нули, например

A =

  -0.2993330   0.1142370   0.0015756  -0.0516769   0.2351970
   0.0530580  -0.1895260   0.0115484   0.1561540  -0.0312342
   0.0065948   0.5349510  -1.0830900   0.5349510   0.0065948
  -0.0312342   0.1561540   0.0115484  -0.1895260   0.0530580
   0.2351970  -0.0516769   0.0015756   0.1142370  -0.2993330

octave:4> inv(A)
ans =

   1.3254e+06   3.8689e+06   8.6360e+04   3.8689e+06   1.3254e+06
   1.3254e+06   3.8689e+06   8.6360e+04   3.8689e+06   1.3254e+06
   1.3254e+06   3.8689e+06   8.6359e+04   3.8689e+06   1.3254e+06
   1.3254e+06   3.8689e+06   8.6360e+04   3.8689e+06   1.3254e+06
   1.3254e+06   3.8689e+06   8.6360e+04   3.8689e+06   1.3254e+06

в моей процедуре результат совсем другой

4.35869e+15 1.27235e+16 2.84009e+14 1.27235e+16 4.35869e+15 
5.11576e+15 1.49335e+16 3.3334e+14 1.49335e+16 5.11576e+15 
4.55463e+15 1.32955e+16 2.96777e+14 1.32955e+16 4.55463e+15 
4.76076e+15 1.38972e+16 3.10208e+14 1.38972e+16 4.76076e+15 
4.66593e+15 1.36204e+16 3.04029e+14 1.36204e+16 4.66593e+15 

Что в этом конкретном случае можно сделать, если

eig(A)
ans =
  -5.9144e-01
  -2.8877e-01
   9.5467e-08
  -8.5874e-02
  -1.0947e+00

cond(A)
ans =  1.7099e+07

[ewert]

Ваша матрица очень плохо обусловлена:

eig(A*A')

ans =

    4.710648587799392e-015
    7.709622222254331e-003
    8.137039870563619e-002
    3.584791817884138e-001
    1.747695310808111e+000

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

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

[Zuzia]

Не хочу создавать отдельную тему из-за пустякового вопроса. Вопрос такой, я написала программу на VBA, но для того что бы она заработала в эксель надо ввести матрицу. Допустим я пишу так Call InputMatrix("лист1", 5, 1, 2, 3, A) и матрицу надо ввести в определенные ячейки листа, вот только в какие? Подскажите пожалуста

[Melevir]

Обратная матрица существует, если определитель не равен нулю.
Чтобы решить проблему с нулями на диагонали обратитесь к модификации метода Гаусса - с выбором ведущего элемента по столбцу или по всей матрице.
Для решения вопроса с обусловленностью матрицы также есть готовые и несложные в реализации методы.