Error in svd in r2018a

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?

[Walter Roberson]

The problem would continue if eps(A) > 1e-14 which would be the case if abs(A) > 100 or so

Answer 1

[Bruno Luong]

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

[Christine Tobler]

Thank you for pointing this out, I have included scaling by norm(A) into the proposed workaround.

Answer 2

[Steven Lord]

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.

[N/A]

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).