% Test accuracy of svd() by a hard problem % % With svd_driver('gesvd') % svd(A) is good upto N==14, with max relative error 1.4e-6. % at N==19, svd(A) max relative error 0.34, start to return wrong answer. % MATLAB behaviours exactly the same as 'gesdd'. %sumsq = @(x) sum(x.*x); % for MATLAB^TM svd_driver('gesvd') %svd_driver('gesdd') N = 26 A = compan (fliplr ([1, 1 ./ cumprod(1:N)])); [u, compan_svd, v] = svd(A); compan_svd = diag(compan_svd); a = A(1,:); % Characteristic polynomial of A*A', note svd(A)=eig(A*A')=roots(c). c = [1, -(1-sumsq(a)), -sumsq(a(1:end-1)), zeros(1, length(a)-2)]; % Numerical stable answer. compan_svd_formula = sort(1-[roots(c(1:3)); zeros(length(a)-2,1)], 'descend') .^0.5; % Solve svd(A) by eig(A*A'). compan_svd_eig = sort(eig(A*A'), 'descend') .^0.5; disp('[compan_svd, compan_svd_eig, compan_svd_formula]'); fprintf('%14.6g, %14.6g, %14.6g\n', ... [compan_svd, compan_svd_eig, compan_svd_formula]'); max_relative_error = max(abs(compan_svd - compan_svd_formula)./compan_svd_formula) % Running result % Octave-4.0.3 (lapack-3.7.0, openblas-0.2.19) %{ N = 26 [compan_svd, compan_svd_eig, compan_svd_formula] 6.08901e+26, 6.08901e+26, 6.08901e+26 3.76175e+09, 1.03585, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 1, 1, 1 0.662326, 0.956064, 1 2.65834e-10, 8.48585e-06, 0.662326 max_relative_error = 3.7617e+09 %}