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

8 次查看（过去 30 天）
N/A
N/A 2022-7-18
回答： Christine Tobler 2022-8-5
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 个评论
显示 1更早的评论
N/A
N/A 2022-7-18
编辑：N/A 2022-7-18
@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).

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Bruno Luong
Bruno Luong 2022-7-18
编辑：Bruno Luong 2022-7-18
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 个评论
N/A
N/A 2022-7-18
编辑：N/A 2022-7-18
@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
Bruno Luong 2022-7-18
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.

请先登录，再进行评论。

Christine Tobler
Christine Tobler 2022-8-5
Ran in:
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.

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by