matlabmatrixsvdinverse

Why pinv answer not equal with svd method answer in Matlab？

a = [1 2 3
     2 4 6
     3 6 9];
b = pinv(a);

[U,S,V] = svd(a); 
T = S;
T(find(S~=0)) = 1./S(find(S~=0));
svda = V * T' * U';

I find that the pinv method in Matlab uses the SVD decomposition to calculate pseudo-inverse, so I tried to solve the matrix a.

And as shown above, theoretically the b should be equal with svda, but the Matlab result said they are totally different. Why?

b is

0.00510204081632653 0.0102040816326531  0.0153061224489796
0.0102040816326531  0.0204081632653061  0.0306122448979592
0.0153061224489796  0.0306122448979592  0.0459183673469388

svda is

-2.25000000000000      -5.69876639328585e+15   3.79917759552390e+15
-2.14021397132170e+15   1.33712246709292e+16  -8.20074512351222e+15
 1.42680931421447e+15  -7.01456098285751e+15   4.20077088383351e+15

How does pinv get to its result?

REASON:

Thanks to Cris, I check my S, and it do have 2 very large number, and that is the source of this strange result.

S:

14.0000000000000    0                       0
0                   1.00758232556386e-15    0
0                   0                       5.23113446604828e-17

By pinv method and Cris method, this 2 latter numbers should set to 0 which I didnt do. So here is the reason。

Solution

  • pinv doesn't just invert all the non-zero values of S, it also removes all the nearly zero values first. From the documentation:

    pinv(A,TOL) treats all singular values of A that are less than TOL as zero. By default, TOL = max(size(A)) * eps(norm(A)).

    pinv more or less does this:

    [U,S,V] = svd(a); 
I = find(abs(S) > max(size(a)) * eps(norm(a)));
T = zeros(size(S));
T(I) = 1 ./ S(I);
svda = V * T.' * U';

    On my machine, isequal(svda,b) is true, which is a bit of a coincidence because the operations we're doing here are not exactly the same as those done by pinv, and you could expect rounding errors to be different. You can see what pinv does exactly by typing edit pinv in MATLAB. It's a pretty short function.

    Note that I used T.', not T'. The former is the transpose, the latter is the complex conjugate transpose (Hermitian transpose). We're dealing with real-valued matrices here, so it doesn't make a difference in this case, but it's important to use the right operator.