`svd` sometimes blows up - how to fix it?

I am having trouble with svd in pinv.m on r2018a, the version of Matlab that I have access to. I was trying to use this workaround. Also, see attached. Usually, it works, but if you repeatedly apply it, e.g.,

trace(pinv_modified(test_matrix))

then sometimes it blows up. How can I get by this? The modified line in pinv is replacing

[U,s,V] = svd(A,'econ');

by

[U,s,V] = svd(A + norm(A)*1e-14*randn(size(A)));

Attached is test_matrix as well.

[KSSV]

I am getting both the values same:

[trace(pinv_modified(test_matrix)) trace(pinv(A))]

ans =

  123.8666  123.8666

[N/A]

@KSSV If you do trace(pinv_modified(test_matrix)) repeatedly, sometimes it blows up and the values are different by 12 orders of magnitude and not even the right sign. This is an issue because I'm doing this calculation many times (on many different matrices).

[Walter Roberson]

original discussion https://www.mathworks.com/matlabcentral/answers/1761920-error-in-svd-in-r2018a?s_tid=srchtitle

Answer 1

[Bruno Luong]

The thread subject is missleading, SVD does NOT blow up, what blowed up is PINV and it's normal given :

>> cond(test_matrix)

ans =

2.7819e+16

Well if you expect that you can get any stable solution without any careful regularized strategy, you must read lterature about "ill-posd problem".

It might be also that your matrix is wrongly built thus you get a singular system.

[N/A]

@Bruno Luong So what specifically would you recommend here? The matrix is calculated from a physical system, so I'm not sure what can be done to modify that.

[Bruno Luong]

I recommend you read about ill-posed system and check your matrix. "The matrix is calculated from a physical system" is not a proof that your matrix is correct : the physics is wrong, the units used is not normalized, one forget to takeinto account the boundary condition or extra equaltion that makes the system wll-posed, a bug in the code, ... or simply what you are trying to solve is not physical possible. The math doesn't lie.

Answer 2

[Christine Tobler]

In short, the problem is that pinv_modified is based on a misunderstanding of the workaround here. The idea is to check if SVD fails and, in that case, compute the svd of a modified matrix (since it's unlikely to fail again there). Always modifying the matrix doesn't prevent SVD from failing on its own, and it introduces some randomness which is causing the problem here.

Here's what's going on: The matrix test_matrix has one singular value which is nearly 0. The pinv function will see this as below its tolerance for a non-zero singular value, and will not invert that last singular value when computing the inverse.

load test.mat
semilogy(svd(test_matrix), 'x')

rank(test_matrix)

ans = 136

trace(pinv(test_matrix))

ans = 123.8666

However, if we now modify this matrix by a random perturbation as done in pinv_modified, sometimes that perturbation will cause the smallest singular value to be too large for pinv to ignore. Basically, pinv will then detect that matrix as being full rank, and invert the smallest singular value. Since this singular value is very small, its inverse is very large and will dominate the output of pinv.

rng default;
for ii=1:10
    trace(pinv_modified(test_matrix))
end

ans = 123.8666

ans = 123.8666

ans = 123.8666

ans = 123.8666

ans = 123.8666

ans = 3.8543e+12

ans = 123.8666

ans = 123.8666

ans = 123.8666

ans = 123.8666

Now how to prevent this? The pinv and rank functions have a tolerance input, telling the functions singular values beneath which value should be treated as zero. To avoid the random perturbation introduced in pinv_modified from being perceived as an important part of your original matrix, set this tolerance to be larger than that perturbation:

rng default;
for ii=1:10
    trace(pinv_modified(test_matrix, 1e-12))
end

ans = 123.8666

ans = 123.8666

ans = 123.8666

ans = 123.8666

ans = 123.8666

ans = 123.8666

ans = 123.8666

ans = 123.8666

ans = 123.8666

ans = 123.8666

Take-away:

1. Copy over the new proposed workaround into pinv_modified to avoid perturbing every matrix instead of only the very rare matrices where SVD otherwise fails in non-updated R2018a.
2. Make sure to adjust the tolerance used in pinv_modified in case the perturbation was necessary, so that the random perturbation doesn't end up dominating the values of the result in the end.