Error in svd in r2018a

4 views (last 30 days)
S
S on 17 Jul 2022
I only have access to r2018a, where there is an error in svd. I tried the recommended fix
[U,s,V] = svd(A + 1e-14*randn(size(A)));
but I still getting the same error every now and then. How do you recommend getting around this error without increasing the size of the perturbation? Could I increase the number of iterations to test convergence, if so how? Or something else?
  1 Comment
Walter Roberson
Walter Roberson on 17 Jul 2022
The problem would continue if eps(A) > 1e-14 which would be the case if abs(A) > 100 or so

Sign in to comment.

Accepted Answer

Bruno Luong
Bruno Luong on 17 Jul 2022
Edited: Bruno Luong on 17 Jul 2022
As Walter has pointed out, the workwaround given by TMW might be still flawed because it ignores the scale of A.
A better solution would be
A = randi(100,4,6)
A = 4×6
82 8 13 18 98 82 77 8 10 61 29 85 28 48 21 23 86 96 70 94 98 63 60 33
[U,s,V] = svd(A + norm(A)*1e-14*randn(size(A)))
U = 4×4
-0.5147 0.4562 -0.1289 -0.7144 -0.4295 0.2837 0.7830 0.3493 -0.4856 0.2094 -0.6069 0.5932 -0.5610 -0.8170 0.0443 -0.1255
s = 4×6
275.2436 0 0 0 0 0 0 110.2492 0 0 0 0 0 0 66.4140 0 0 0 0 0 0 42.8050 0 0
V = 6×6
-0.4656 0.0719 0.5394 -0.5573 -0.0481 -0.4174 -0.3037 -0.5517 -0.2972 0.3213 0.2364 -0.5966 -0.2767 -0.6068 -0.0339 -0.1317 -0.6095 0.4065 -0.2978 -0.1917 0.5160 0.3314 0.5339 0.4619 -0.5026 0.1989 -0.5943 -0.3829 0.3450 0.2984 -0.5226 0.4959 -0.0125 0.5589 -0.4078 -0.0469
Or might be this cumbersome code that uses EIG (hopefully not buggy in your version) instead of SVD.
[m,n] = size(A);
AAc = A*A'; [U,s2u] = eig(1/2*(AAc+AAc'));
AcA = A'*A; [V,s2v] = eig(1/2*(AcA+AcA'));
[~, isu] = sort(sqrt(max(diag(s2u),0)),'descend');
[s, isv] = sort(sqrt(max(diag(s2v),0)),'descend');
U = U(:,isu);
S = diag(s);
if n < m
S(m,1) = 0;
isv = isv(1:n);
else
S = S(1:m,:);
end
V = V(:,isv);
U
U = 4×4
0.5147 0.4562 0.1289 -0.7144 0.4295 0.2837 -0.7830 0.3493 0.4856 0.2094 0.6069 0.5932 0.5610 -0.8170 -0.0443 -0.1255
S
S = 4×6
275.2436 0 0 0 0 0 0 110.2492 0 0 0 0 0 0 66.4140 0 0 0 0 0 0 42.8050 0 0
V % Note the two last columns are arbitrary provided they span the same subspace and are orthogonal
V = 6×6
0.4656 0.0719 -0.5394 -0.5573 -0.4065 -0.1064 0.3037 -0.5517 0.2972 0.3213 -0.6239 0.1500 0.2767 -0.6068 0.0339 -0.1317 0.4883 -0.5462 0.2978 -0.1917 -0.5160 0.3314 0.3822 0.5936 0.5026 0.1989 0.5943 -0.3829 0.2469 0.3836 0.5226 0.4959 0.0125 0.5589 0.0110 -0.4103
  1 Comment
Christine Tobler
Christine Tobler on 5 Aug 2022
Thank you for pointing this out, I have included scaling by norm(A) into the proposed workaround.

Sign in to comment.

More Answers (1)

Steven Lord
Steven Lord on 17 Jul 2022
Are you using Update 3 of release R2018a or a later Update, or are you using release R2018a (no Update) or Updates 1 or 2?
If you have not installed Update 3 or a later Update, please do.
  1 Comment
S
S on 17 Jul 2022
Edited: S on 17 Jul 2022
Thanks @Steven Lord. It's the version R2018a (9.4.0.813654) 64-bit (glnxa64) from Feb. 23, 2018. It's on a cluster where I don't have admin access (I'm not able to install a new version).

Sign in to comment.

Tags

Products


Release

R2018a

Community Treasure Hunt

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

Start Hunting!