channelflow.org

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== exercise 16.2 ======

{{:unh2010:iam931:svderrs.png?400}} {{:unh2010:iam931:qrerrs.png?400}}

These plots show (left) the normwise errors of a numerically computed SVD of randomly constructed matrices A with known SVD factors, as a function of the condition number of A, (right) ditto for QR decomposition. For the SVD, USV' = A, the errors in the computed unitary factors U and V appear to be bounded above linearly with condition number, suggesting that the accuracy of the SVD computation follows the error scaling law for backward stable algorithms (Trefethen eqn 15.1) --even though on dimensionality grounds the computation of f: A -> U,S,V cannot backward stable! Also, the singular values are computed to nearly machine precision, even when the condition number is O(10^18).

The error scalings of the computed Q,R factors of the QR decomposition are similar to those of the SVD, though they are an order of magnitude or so worse. Also, whereas the errors in the SVD factors are sometimes spread well below the apparent upper bound, the errors in Q and R are almost always within an order of magnitude of linear-in-condition number

For both, the normwise error of the product of the factors (QR or USV') stays right down near machine epsilon, independently of condition number.


<code>
% A matlab program that constructs a matrix via SVD factors and then compares those to the computed SVD

m = 50;  % dimension of matrices
S = 100; % number of tests 

for s = 1:S 

  % construct A from random SVD factors
  [U, tmp] = qr(randn(m,m));
  [V, tmp] = qr(randn(m,m));
  S = diag(sort(rand(m,1).^6, 'descend'));
  A = U*S*V';
  
  % compute SVD of A
  [U2,S2,V2] = svd(A);

  % change signs of cols of U2 if they don't match U, ditto V
  for k=1:m
      uk = U(:,k);
      if uk'*U2(:,k) < 0
         U2(:,k) = -1*U2(:,k);
         V2(:,k) = -1*V2(:,k);
      end
  end
  kappa = cond(A);
  loglog(kappa, norm(U-U2), 'r.')
  hold on
  loglog(kappa, norm(V-V2), 'b.')
  loglog(kappa, norm(S-S2), 'k.')
  loglog(kappa, norm(A-U2*S2*V2'), 'g.');
end

legend('norm(U-U2)','norm(V-V2)','norm(S-S2)', 'norm(A-U2*S2*V2)''', ...
       'Location', 'NorthWest');
legend boxoff
xlabel('cond(A)')
axis([1 10^20 10^-16 1])
set(gcf, 'color', [1 1 1])
set(gcf,'inverthardcopy','off')

title('SVD errors, A has singvals (0<random<1).^6')
</code>

To modify this program to do similar tests on QR decomposition, substitute in this code segment
<code>
  R = triu(randn(m,m));
    [Q, tmp] = qr(randn(m,m));

    for k=1:m
      if R(k,k) < 0
         R(k,k) = -1*R(k,k);
      end
    end
</code>
and appropriate changes in the plotting and labeling.
channelflow.org

User Tools

Site Tools

Page Tools
