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

7 views (last 30 days)
S
S on 18 Jul 2022
Answered: Christine Tobler on 5 Aug 2022
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.
  3 Comments

Sign in to comment.

Answers (2)

Bruno Luong
Bruno Luong on 18 Jul 2022
Edited: Bruno Luong on 18 Jul 2022
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.
  2 Comments
Bruno Luong
Bruno Luong on 18 Jul 2022
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.

Sign in to comment.


Christine Tobler
Christine Tobler on 5 Aug 2022
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.

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Tags

Products


Release

R2018a

Community Treasure Hunt

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

Start Hunting!