pinv failing on single precision matrices

1 次查看(过去 30 天)
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?
  1 个评论
John D'Errico
John D'Errico 2025-3-12
编辑:John D'Errico 2025-3-12
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.

请先登录,再进行评论。

采纳的回答

Matt J
Matt J 2025-3-12
编辑:Matt J 2025-3-13
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>
  1 个评论
Bob
Bob 2025-3-13
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!

请先登录,再进行评论。

更多回答(1 个)

Harald
Harald 2025-3-12
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
  2 个评论
Bob
Bob 2025-3-12
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
Harald 2025-3-12
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.
The colleagues will then assess and contact the proper development team if needed.
Best wishes,
Harald

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Operating on Diagonal Matrices 的更多信息

产品


版本

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by