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

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.

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。

CodePudding user response：

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.