Let x = {x1,…,xn} and y = {y1,…,xm} be n and m equispaced gridpoints on [-1, 1] respectively.

(a) Let p(x) be a n-1-th order polynomial. Derive a formula for the m x n matrix A that maps the n-vector px = (p(x0), … p(xn)) to the m-vector py = (p(y0), … p(yn)). (Read py as a vector marked with a subscript y). We can interpret this as forming a polynomial interpolant p(x) to the “data” px and evaluating the interpolating polynomial at the points {y1,…,xm} to get py = (p(y0), … p(yn)).

An n-1-th order polynomial p(x) = c0 + c1 x + … c{n-1} x^{n-1} is defined by n coefficients c = (c0, c1, …, c{n-1}). Evaluating a polynomial at a given set of gridpoints is equivalent to multiplying c by the Vandermonde matrix Vx for the gridpoints x, as defined in example 1.1 (and where Vx is a matrix with a x subscript. Computing the coefficients from the data is equivalent to solving c = Vx^{-1} px.

Now to evaluate the interpolating polynomial at the gridpoints {y1,…,xm}, we would like to multiply c by the Vandermonde matrix for the points y, i.e. py = Vy c. However, note that the Vy will be an m x m matrix, whereas c is n-dimensional. Let's assume m > n. We can either define c_j = 0 for n-1<j<m-1 and compute py = Vy c, or, equivalently, remove the last m-n columns of Vy so that it is m x n and then compute py = Vy c.

Assuming that Vy is so truncated, we have py = Vy c = Vy Vx^{-1} px, so the linear map from px to py is Vy Vx^{-1}.

(b) Matlab code to compute this matrix is

function A = ex12_2a(n)
m = 2*n-1;

X = vander(linspace(-1, 1, n));
Y = vander(linspace(-1, 1, m));
Y = Y(:,(m-n+1):m);  % argh, matlab defines Vandermonde with reversed column order

A = Y*inv(X);

(b) The inf-norm of A is plotted as a function of n on the left. (c) The inf-norm condition number of interpolating the constant function 1 from n to m equispaced points is the condition number of the matrix multiplication A px

The simplification occurs because both px = (1 1 1 1 …) and A px have unit inf-norms.

(d) On the left is a reproduction of fig 11.1 using A for m x n = 21 x 11 and px = [0 0 0 1 1 1 0 0 0 0 0]. The interpolant looks the same as in the book (good!) and the inf-norm (max) amplification for this value of px appears to be about 4. This is somewhat smaller than 24.7, the inf-norm condition number of interpolating the constant function 1, but that's ok. The condition number of interpolation is the sup over the space of data px; the particular px above does not achieve the sup. On the right is a similar plot for px = [1 1 1 1 1 1 1 1 1 1 1] , though on the vertical axis we plot p(x)-1 so that small deviations from unity are visible. I would expect the errors in the interpolating polynomial to be O(kappa(A) epsilon_machine) = 24 * 2^-52 = 5e-15 –that is, we think of the constant function 1 and the matrix multiplication as being perturbed by machine-epsilon rounding errors, which are then amplified by a factor kappa(A). But the errors are two orders of magnitude larger! Why is this?

Error analysis shows that the numerical estimate of e computed from

is a stable whether summed rightwards (largest terms first) or leftwards (smallest terms first). The above plots conforms this, by showing the absolute value of the error as a function of simulated machine precision. Both rightwards and leftwards sums have error that scale as machine precision.

Calculations at a given machine precision eps were simulated by multiplying each floating point operation in the formula by (1 + rand*eps), where rand is a random number uniformly distributed on [-1, 1]. The sums were truncated at the smallest K for which 1/k! < eps.

Matlab code for leftwards summation

function x = eleftwards(eps);

x = 0;

n = 0;
for k=0:100;
  if 1/factorial(k) > eps
    n = k;
  else
    break;
  end
end 

for k=n:-1:0
  factk = 1;
  for i=1:k
    factk = factk*i*(1 + randu*eps);
  end
  invfactk = 1/factk*(1 + randu*eps);
  x = (x + invfactk)*(1 + randu*eps);
end

the function randu produces a random uniform on [-1,1] from matlab's rand on [0,1]

function r = randu 
r = 2*rand - 1;

the driver program for plotting results

% plot error of sum 1/k! computation of e as a function of 
% simulated machine precision eps, using both leftards and rightwards 
% summation

hold off
e = exp(1);
mf = 'markerface';

for n=0:19;
    eps = 10^-n;
    loglog(eps, abs(e-eleftwards(eps)), 'ro', mf, 'r');
    hold on
    loglog(eps, abs(e-erightwards(eps)), 'bo', mf, 'b');
end
plot([10^-20 1], [10^-20 1], 'k') 

xlabel('machine epsilon')
ylabel('| fl(e) - sum 1/k! |');
legend('leftwards sum', 'rightwards sum', 'machine epsilon', 'Location', 'northwest')
legend boxoff

Note that the structure of the driver program and eleftwards/rightwards functions ends up calculating the same partial sums many times over, but I'm optimizing for my time, not CPU time.

(b) These plots show (left) the normwise errors of a numerically computed SVD of randomly constructed matrices A with known SVD factors (with the signs of the columns of U and V fixed), 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.

% 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')

To modify this program to do similar tests on QR decomposition, substitute in this code segment

  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

and appropriate changes in the plotting and labeling.