Why is the output of inv() and pinv() not equal in Matlab and Octave?

Question

I have noticed that if A is a NxN matrix and it has the inverse matrix. But what the inv() and pinv() function output is different. - My environment is Win7x64 SP1, Matlab R2012a, Cygwin Octave 3.6.4, FreeMat 4.2

Have a look at the examples from 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

It's all the same ans result by running the same command above in Matlab.

And I calculate inv(A)*A or A*inv(A), the result is identity 3x3 matrix in both Octave and Matlab.
The result of A*pinv(A) and pinv(A)*A are identity 3x3 matrix in Matlab and FreeMat.
The result of A*pinv(A) is identity 3x3 matrix in Octave.
The result of pinv(A)*A is not identity 3x3 matrix in Octave.

I don't know the reason why inv(A) != pinv(A), I have considered the details of the element in the matrix. It seems to be the floating accuracy problem which causes this problem.

The 10+ digits after the dot point may be different like this:

6.65858991579923298331777914427220821380615200000000 element in inv(A)(1,1) against
6.65858991579923209513935944414697587490081800000000 element in pinv(A)(1,1)

possible duplicate of [Why is Matlab's inv slow and inaccurate?](http://stackoverflow.com/questions/1419580/why-is-matlabs-inv-slow-and-inaccurate) — Shai, Oct 17 '13 at 10:26
@Shai, I believe OP might benefit from reading the answers to the question you linked to (at least if OP is using `inv` for solving `x = A^-1*b`), but IMO this is not a duplicate. — Stewie Griffin, Oct 17 '13 at 10:42

goth · Answer 1 · 2015-08-10T19:05:24.707

This question is quite old, but I'll answer it anyway because it appears almost on top in some google searches.

I'll use for my example the magic(N) function which returns an N-by-N magic square.

I'll create a 3x3 magic square M3, take the pseudoinverse PI_M3 and multiply them:

   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

As you can see the answer is the identity matrix save for some rounding errors. I'll repeat the operation with a 4x4 magic square:

   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

The result is not the identity matrix, it means that the 4x4 magic square does not have an inverse. I can verify this by trying one of the rules of the Moore-Penrose pseudoinverse:

   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

The rule A*B*A = A is satisfied. This shows that pinv returns the inverse matrix when it is available and the pseudoinverse when the inverse is not available. This is the reason why in some situations you get a small difference, just some rounding errors, and in other situations you get a bigger difference. To show it I'll get the inverse of both magic quadrants and subtract them from the pseudoinverse:

   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

score 9 · Answer 2 · edited May 23 '17 at 12:34

9

It seems to me like you answered your own question in the bottom here. The reason is floating point arithmetic. The algortihms for inv() and pinv() are not exactly the same, as pinv() must be able to handle non-square matrices. Therefore the answers will not be exactly the same.

If you look at the value of pinv(A)*A, you will see that it's very very close to the identity matrix.

I get:

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

Instead of comparing the matrices with ==, use < tolerance_limit

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

(c-d) < 1e-10

Sidenote:

x = A^-1*b should not be solved x = inv(A)*b;, but rather x = A \ b; See the link Shai posted for explanations.

edited May 23 '17 at 12:34

Community

1
1

answered Oct 17 '13 at 10:34

Stewie Griffin

14,889
11
39
70

1

What I meant to be is that pinv(A)*A should be identity matrix while now the Octave don't give me this answer, I know why this tiny result comes out, but shouldn't we always get matrix multipy it's inverse matrix equal to identity matrix just like that the Matlab and Freemat now is doing? – myme5261314 Oct 17 '13 at 11:44
What does Octave give you? I guess a result _very close_ to the identity matrix (as my example). To understand why this is, Daniel's answer is quite explanatory. Octave first calculates the inverse than multiplies with the original matrix. If you get something completely different, please give an example. – Stewie Griffin Oct 17 '13 at 12:16
I get the ans matrix just similar to your ans in your answer. What I mean is that as a product or say as a library, the Octave should let pinv(A)*A=identity matrix which make sense, but now the ans matrix it returns back doesn't make sense, and pinv(A)*A equals to identity matrix in Matlab and FreeMat. – myme5261314 Oct 17 '13 at 13:05
1

@myme5261314 The inverse is probably calculated differently in Octave, Matlab and FreeMat, therefore the answers differ. It makes sense if you consider it as two separate calculations, first calculate the inverse of `A`, then multiply it by `A`. As the first calculation is not accurate, the answer won't either. Octave doesn't recognize `pinv(A)*A` as the identity matrix (as we do mathematically), it actually does the calculations. – Stewie Griffin Oct 17 '13 at 13:19

score 2 · Answer 3 · answered Oct 17 '13 at 10:06

Floating point arithmetic has a certain precision, you can not rely on equality. To avoid such errors, you could try to work with the symbolic toolbox of matlab.

Very simple line of code in octave to demonstrate the problems:

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

user2874626 · Answer 4 · 2021-11-04T14:45:06.993

0

I computed the pseudo-inverse matrix of A=[1,1,0;1,0,1;2,1,1] (which is of rank 2) by using pinv(A). A*pinv(A) gives a non-identity matrix, which is A*pinv(A)=[0.667, -0.333, 0.333; -0.333, 0.667,0.333; 0.333, 0.333, 0.667]. I think for singular matrices, it is better to compute the pseudo-inverse matrix using svd() manually.

Here is some update: A*pinv(A) itself may be different to identity matrix because of it is non-full rank. It may not depends on the Matlab algorithm of pinv(A).

edited Nov 04 '21 at 14:45

answered Nov 03 '21 at 15:46

user2874626

13
4

`A` has no inverse, which means you cannot find a matrix `B` such that `A*B` is the identity matrix. So it is no wonder than `A*pinv(A)` is not the identity matrix. – Cris Luengo Nov 04 '21 at 15:06

Why is the output of inv() and pinv() not equal in Matlab and Octave?

4 Answers4

Linked

Related