Почему выходные данные inv() и pinv() не равны в Matlab и Octave?

Я заметил, что если A является матрицей NxN и имеет обратную матрицу. Но то, что вывод функции inv() и pinv() отличается.  - Моя среда - Win7x64 SP1, Matlab R2012a, Cygwin Octave 3.6.4, FreeMat 4.2

Взгляните на примеры из Octave:

A = rand(3,3)
A =
0.185987   0.192125   0.046346
0.140710   0.351007   0.236889
0.155899   0.107302   0.300623

pinv(A) == inv(A)
ans =
0 0 0
0 0 0
0 0 0
  • Это все равно ans результат, выполнив ту же самую команду в Matlab.
  • И я вычисляю inv(A)*A или A*inv(A), результатом является матрица тождества 3x3 как в Octave, так и в Matlab.
  • Результатом A*pinv(A) и pinv(A)*A являются матрица тождества 3x3 в Matlab и FreeMat.
  • Результатом A*pinv(A) является тождественная матрица 3x3 в Octave.
  • Результатом pinv(A)*A является не тождественная матрица 3x3 в Octave.

Я не знаю причины, почему inv(A) != pinv(A), я рассмотрел детали элемента в матрице. Это проблема с плавающей точностью, которая вызывает эту проблему.

10-значные цифры после точки точки могут быть разными:

  • 6.65858991579923298331777914427220821380615200000000 элемент inv(A)(1,1) против

  • 6.65858991579923209513935944414697587490081800000000 элемент в pinv(A)(1,1)

Ответы

Ответ 1

Этот вопрос довольно старый, но я все равно отвечу, потому что он появляется почти в верхней части некоторых поисковых запросов Google.

Я использую для своего примера функцию magic (N), которая возвращает магический квадрат N-by-N.

Я создам магический квадрат M3 3x3, возьму псевдоинверсивный PI_M3 и размножу их:

   prompt_$ M3 = magic(3) , PI_M3 = pinv(M3) , M3 * PI_M3
  M3 =

    8   1   6
    3   5   7
    4   9   2

  PI_M3 =

     0.147222  -0.144444   0.063889
    -0.061111   0.022222   0.105556
    -0.019444   0.188889  -0.102778

  ans =

     1.0000e+00  -1.2212e-14   6.3283e-15
     5.5511e-17   1.0000e+00  -2.2204e-16
    -5.9952e-15   1.2268e-14   1.0000e+00

Как вы можете видеть, ответ - это матрица идентичности, сохраняющая некоторые ошибки округления. Я повторю операцию с магическим квадратом 4x4:

   prompt_$ M4 = magic(4) , PI_M4 = pinv(M4) , M4 * PI_M4
 
  M4 =

     16    2    3   13
      5   11   10    8
      9    7    6   12
      4   14   15    1

  PI_M4 =

     0.1011029  -0.0738971  -0.0613971   0.0636029
    -0.0363971   0.0386029   0.0261029   0.0011029
     0.0136029  -0.0113971  -0.0238971   0.0511029
    -0.0488971   0.0761029   0.0886029  -0.0863971

  ans =

     0.950000  -0.150000   0.150000   0.050000
    -0.150000   0.550000   0.450000   0.150000
     0.150000   0.450000   0.550000  -0.150000
     0.050000   0.150000  -0.150000   0.950000

Результат не является тождественной матрицей, это означает, что магический квадрат 4x4 не имеет обратного. Я могу проверить это, попробовав одно из правил псевдорежима Мура-Пенроуза:

   prompt_$ M4 * PI_M4 * M4
  
ans =

   16.00000    2.00000    3.00000   13.00000
    5.00000   11.00000   10.00000    8.00000
    9.00000    7.00000    6.00000   12.00000
    4.00000   14.00000   15.00000    1.00000

Правило A * B * A = A выполняется. Это показывает, что pinv возвращает обратную матрицу, когда она доступна, и псевдообратная, когда инверсная недоступна. Вот почему в некоторых ситуациях вы получаете небольшую разницу, только некоторые ошибки округления, а в других ситуациях вы получаете большую разницу. Чтобы показать это, я получу инверсию обоих магических квадрантов и вычитаю их из псевдоинверса:

   prompt_$ I_M3 = inv(M3) , I_M4 = inv(M4) , DIFF_M3 = PI_M3 - I_M3, DIFF_M4 = PI_M4 - I_M4
  I_M3 =

     0.147222  -0.144444   0.063889
    -0.061111   0.022222   0.105556
    -0.019444   0.188889  -0.102778

  warning: inverse: matrix singular to machine precision, rcond = 1.30614e-17
  I_M4 =

     9.3825e+13   2.8147e+14  -2.8147e+14  -9.3825e+13
     2.8147e+14   8.4442e+14  -8.4442e+14  -2.8147e+14
    -2.8147e+14  -8.4442e+14   8.4442e+14   2.8147e+14
    -9.3825e+13  -2.8147e+14   2.8147e+14   9.3825e+13

  DIFF_M3 =

     4.7184e-16  -1.0270e-15   5.5511e-16
    -9.9226e-16   2.0470e-15  -1.0825e-15
     5.2042e-16  -1.0270e-15   4.9960e-16

  DIFF_M4 =

    -9.3825e+13  -2.8147e+14   2.8147e+14   9.3825e+13
    -2.8147e+14  -8.4442e+14   8.4442e+14   2.8147e+14
     2.8147e+14   8.4442e+14  -8.4442e+14  -2.8147e+14
     9.3825e+13   2.8147e+14  -2.8147e+14  -9.3825e+13

Ответ 2

Мне кажется, что вы ответили на свой вопрос внизу. Причина - арифметика с плавающей запятой. Алгоритмы для inv() и pinv() не совсем совпадают, так как pinv() должен иметь возможность обрабатывать неквадратные матрицы. Поэтому ответы не будут одинаковыми.

Если вы посмотрите на значение pinv(A)*A, вы увидите, что оно очень близко к единице.

Я получаю:

ans =

   1.0000e+00   6.1062e-16  -3.0809e-15
  -5.8877e-15   1.0000e+00   6.3942e-15
   2.4425e-15  -3.0184e-16   1.0000e+00

Вместо сравнения матриц с == используйте < tolerance_limit

c = A*pinv(A);
d = pinv(A)*A;

(c-d) < 1e-10

Sidenote:

x = A^-1*b не следует решать x = inv(A)*b;, а скорее x = A \ b; см. ссылку Shai для объяснений.

Ответ 3

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

Очень простая строка кода в октаве, чтобы продемонстрировать проблемы:

>>> (1/48)*48==(1/49)*49
ans = 0
>>> (1/48)*48-(1/49)*49
ans =  1.1102e-16
>>>