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

조회 수: 9(최근 30일)
S 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표시숨기기 이전 댓글 수: 2
Walter Roberson 2022년 7월 18일

댓글을 달려면 로그인하십시오.

### 답변(2개)

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표시숨기기 이전 댓글 수: 1
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 2022년 8월 5일
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.

댓글을 달려면 로그인하십시오.

### 범주

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

R2018a

### Community Treasure Hunt

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

Start Hunting!

Translated by