Singular value decomposition and inverse of square matrix

I've previously touched the subject in this question. There (and subsequently on other places), I've learned that if a SVD is applied to a square matrix $M$, $M=USV^T$, then the inverse of $M$ is relatively easy to calculate as $M^{-1}=V S^{-1}U^T$. I've implemented the SVD algorithm and began to receive wrong results, so I fed my test examples to Matlab and was surprised to find that $M^{-1}=V S^{-1}U^T$ apparently doesn't hold.

So, my question is am I calculating the inverse of a matrix based on its SVD correctly? Am I missing something there?

Please note that I'm not asking for debugging help, seeking bugs in Matlab, etc. I'm just perplexed that the equation $M^{-1}=V S^{-1}U^T$ doesn't hold for some reason, as, when calculated separately, $M^{-1}$ and $V S^{-1}U^T$ give different results.

Here is an example (the results come from Matlab and have nothing to do with any implementation of mine):

 M = 32.7276   -5.0470   -5.3461   -1.7619
-5.0470   10.1665   -5.1195   -2.0058
-5.3461   -5.1195   38.7891   10.4173
1.7619    2.0058  -10.6087   38.5192


$M=USV^T$ (obtained with the command [U S V] = svd(M))

U = -0.4313   -0.0317    0.8703    0.2355
-0.0785   -0.0293   -0.2974    0.9511
0.8860    0.1502    0.3905    0.1999
-0.1509    0.9877   -0.0402    0.0054

S = 43.3263         0         0         0
0   39.9753         0         0
0         0   31.9654         0
0         0         0    7.8516

V = -0.4321    0.0012    0.8705    0.2356
-0.0799    0.0269   -0.2971    0.9511
0.8927   -0.1084    0.3893    0.1996
0.1001    0.9937    0.0495   -0.0042


Now, $V S^{-1}U^T$ (obtained with the command V.' * inv(S) * U.', and I am aware that this is a non-conjugate transpose, however the case is a real matrix) yields:

0.0317    0.0047    0.0043   -0.0015
0.0268    0.1214    0.0241    0.0015
0.0037    0.0010    0.0227   -0.0108
0.0022   -0.0035    0.0107    0.0224


While a direct inverse of $M$ (command inv(M)) yields:

 0.0351    0.0212    0.0078    0.0006
0.0212    0.1181    0.0191    0.0020
0.0078    0.0190    0.0277   -0.0061
-0.0006   -0.0019    0.0063    0.0241


The two should be the same, but clearly are not.

• Check the line "Now, $VS^{-1}U^T$ (...)", either you wrote it wrong here or there is (at least of) your mistake(s). – Nigel Overmars Sep 24 '16 at 20:33
• I haven't checked the arithmetic but I wonder if the discrepancy might be explained by rounding too early. $\qquad$ – Michael Hardy Sep 24 '16 at 20:42
• @NigelOvermars Yes, there was a typo in the command, I've edited it to correct – user3209815 Sep 24 '16 at 22:22
• @MichaelHardy That comment is actually very relevant to me, as the precision of the results can have impact, especially in large matrices. Would you care to elaborate or to point me in some reading direction? – user3209815 Sep 24 '16 at 22:24

Now, $VS^{−1}U^T$ (obtained with the command U.' * inv(S) * V.'

This is wrong, the command should be V*inv(S)*U', which yields the answer you are looking for.

>> [U,S,V] = svd(M)
inv(M) - V*inv(S)*U'

ans =

1.0e-15 *

0.0625    0.0243   -0.0954    0.0181
0.0243    0.1388    0.0451   -0.0017
0.0408    0.0486   -0.0278    0.0052
0.0035   -0.0004   -0.0026         0

• I don't understand why isn't V transposed again? SVD gives the already transposed $V^T$ into the variable V, so to invert it you have to transpose the variable V (technically $(V^T)^T$). – user3209815 Sep 24 '16 at 22:20
• You are just making a mistake in your code, please check the highlighted area in my answer again. @user3209815 – Nigel Overmars Sep 24 '16 at 22:24
• I see the error, I just don't understand why it is V*inv(S)*U' instead of V'*inv(S)*U'. You have to transpose the V part just as the U part, right? – user3209815 Sep 24 '16 at 22:26
• Matlab returns $V$, not $V^T$ if [U,S,V] = svd(M). See also here:nl.mathworks.com/help/matlab/ref/svd.html – Nigel Overmars Sep 24 '16 at 22:31