2016-11-17 5 views
1

Примечание: этот вопрос связан с my previous question here on Matrix right divisionматрица между MATLAB и .NET

Когда я сравниваю конечный результат в следующем примере в обоих MATLAB и C#, я заметил, что есть заметная разница. Почему это так?

Результат поиска матричных инверсий, похоже, подсчитан, но A * A^-1, похоже, ушел.

Пример в MATLAB:

>> a = [1 2 3; 4 5 6; 7 8 10] 

a = 

    1  2  3 
    4  5  6 
    7  8 10 
>> inv(a) 

ans = 

    -0.6667 -1.3333 1.0000 
    -0.6667 3.6667 -2.0000 
    1.0000 -2.0000 1.0000 
>> z = mtimes(a, inv(a)) 
>> z 

z = 

    1.0000e+00 -4.4409e-16 -1.1102e-16 
    1.3323e-15 1.0000e+00 -2.2204e-16 
    2.2204e-15 -2.6645e-15 1.0000e+00 

те же данные в C#: using CSML Matrix Library

//using CSML Matrix Library 
public static Matrix operator *(Matrix A, Matrix B) 
{ 
    if (A.ColumnCount != B.RowCount) 
    throw new ArgumentException("Inner matrix dimensions must agree."); 
    Matrix C = new Matrix(A.RowCount, B.ColumnCount); 

     for (int i = 1; i <= A.RowCount; i++) 
     { 
      for (int j = 1; j <= B.ColumnCount; j++) 
      { 
       C[i, j] = Dot(A.Row(i), B.Column(j)); 
      } 
     } 

    return C; 
} 
Console.WriteLine(e1); 
1;  2;  3;  \ 
4;  5;  6;  \ 
7;  8;  10;  \ 

Console.WriteLine(e1.Inverse()); 
-0.666666666666667;  -1.33333333333333;  1;  \ 
-0.666666666666669;  3.66666666666667;  -2;  \ 
1;  -2;  1;  \ 

Console.WriteLine(e1 * e1.Inverse()); 
0.999999999999999;  1.77635683940025E-15; -8.88178419700125E-16; \ 
-5.32907051820075E-15; 1.00000000000001;  -3.5527136788005E-15; \ 
-1.06581410364015E-14; 3.5527136788005E-15; 0.999999999999996;  \ 
+2

Различия кажутся очень маленькими - по порядку точности машины, 'eps (1)'. Даже если обратный алгоритм был реализован идентично, вы все еще используете два совершенно разных языка/компиляторы и должны ожидать таких различий при работе в плавающей точке. – horchler

ответ

2

Оба результата, кажется разумным. MATLAB вычисляет inv(A) с использованием сокращения строк. Каждый шаг сокращения строк приводит к числовой ошибке (MATLAB интерпретирует 1/3 как десятичную с конечным числом десятичных знаков). Поэтому из-за числовой ошибки я ожидал бы, что элементы inv(A) будут отключены на 10^{-16}. Результат, который

а * INV (а) = 1 +/- 10^{- 16} вдоль диагоналей и

а * INV (а) = +/- 10^{- 16} вдоль off-diagonals

соответствует a*inv(a), равным единичной матрице плюс или минус некоторая численная ошибка.

Смежные вопросы