pinv failing on single precision matrices

I just noticed that pinv is not giving the expected result on single precision matrices.

A

A =

  4×4 single matrix

   1.0e+03 *

    0.0010    0.0000    0.0000    0.0537
   -0.0000    0.0010   -0.0000   -0.0714
   -0.0000    0.0000    0.0010   -2.4149
         0         0         0    0.0010        

Its inverse and pseudo inverses are computed as

inv(A)

ans =

  4×4 single matrix

   1.0e+03 *

    0.0010   -0.0000   -0.0000   -0.0563
    0.0000    0.0010    0.0000    0.0724
    0.0000   -0.0000    0.0010    2.4148
         0         0         0    0.0010
         
pinv(A)

ans =

  4×4 single matrix

    0.9995    0.0003    0.0222   -0.0000
    0.0010    0.9991   -0.0295   -0.0000
    0.0233   -0.0300    0.0014   -0.0000
    0.0000   -0.0000   -0.0004    0.0000

The pseudo inverse is not matching the inverse and hence incorrect - as A is an invertible matrix.

Making A a double precision matrix the pseudo inverse gives the expected result

Adouble = double(A)

Adouble =

   1.0e+03 *

    0.0010    0.0000    0.0000    0.0537
   -0.0000    0.0010   -0.0000   -0.0714
   -0.0000    0.0000    0.0010   -2.4149
         0         0         0    0.0010
         
pinv(Adouble)

ans =

   1.0e+03 *

    0.0010   -0.0000   -0.0000   -0.0563
    0.0000    0.0010    0.0000    0.0724
    0.0000   -0.0000    0.0010    2.4148
   -0.0000    0.0000   -0.0000    0.0010

This is not expected behaviour of pinv I suppose?

[John D'Errico]

The thing is, we can't really know what the expected result should have been, since we don't know the true matrix. You want to make it possible to get help, by posting the matrix in a .mat file, not as only what was displayed by default. pinv is a pretty simple code, so the difference between what you got and what should have happened will be easy to track down.

I would guess the matrix has a singular value very near the cutoff point (but just below it), and depending on the tolerance for pinv relative to that singular value, then pinv effectively decides the matrix is singular in single precision. In that case, pinv zeros out that singular value when computing the pseudo-inverse.

The above scenario does not happen when you convert to double though, so everything works. But this is not really a bug in pinv at all, merely a reflection of the tolerance chosen.

Answer 1

[Matt J]

It doesn't look like a bug. Contrary to what you claim, the matrix is not very invertible, showing a condition number of ~1e7,

format long

A =[ 0.001000000000000                   0                   0  -0.056300000000000
                   0   0.001000000000000                   0   0.072400000000000
                   0                   0   0.001000000000000   2.414800000000000
                   0                   0                   0   0.001000000000000];

cond(A)

ans = 
     5.839672489999828e+06

svd(A)

ans = 4×1

2.416541431467673 0.001000000000000 0.001000000000000 0.000000413814548

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Bob]

Right! I did not check the condition number assuming that is would not be that large. Then the behaviour of pinv seems consistent to expected numerical errors using single precision values.

Thank you for the clarification!

Answer 2

[Harald]

Hi,

I find the display to be confusing here. Try

format shortG

and you may find the display of the correct inverse to be easier on the eyes than the exponential notation.

Still, there is a big deviation in the last column and the bottom right. This might be due to the somewhat large norm of the matrix and the resulting default tolerance. Try lowering the tolerance, and results will still not be the same but much more comparable:

pinv(A, 1e-7)

Best wishes,

Harald

[Bob]

Thank you for your reply!

But actually I already found as solution to use double precision matrices. My concern was to report this unexpected behaviour of pinv to create awareness for other users and allowing the Matlab developers to fix it if they agree to my judgement.

[Harald]

Thanks Bob, I misunderstood and apologize for that.

If using doubles is acceptable to you, then this is certainly a good idea.

If you suspect a bug, my advice would generally be to contact the Technical Support team.

https://www.mathworks.com/support/contact_us.html

The colleagues will then assess and contact the proper development team if needed.

Best wishes,

Harald