ACM11 Homework 2 Solutions

Problem 1: Taylor series
A "magic" fix to method (b) in Problem 1
Problem 2: Fair dice
Problem 3: RANDU
Problem 4: Vectorization and Anonymous Functions
Problem 5(a): Fourier Interpolation
Problem 5(b): Polynomial Interpolation

Problem 1: Taylor series

This problem shows an example of overflow when we use method (c). Method (b) isn't quite as prone to numerical errors because the numerator (10^x) and the denominator (x!) cancel partially, and then we convert to a single precision. However, method (b) still returns a useless answer that doesn't have any correct digits.

x = 0:100;
v1=(-10).^x./factorial(x); s1=sum(v1);
v2=single( (-10).^x./factorial(x) ); s2=sum(v2);
v3=single( (-10).^x )./factorial(x); s3=sum(v3);
fprintf('Exact value: %.8e\n',exp(-10))
fprintf('Method (a):  %.8e\n',s1)
fprintf('Method (b):  %.8e\n',s2)
fprintf('Method (c):  %.8e\n',s3)

Exact value: 4.53999298e-05
Method (a):  4.53999294e-05
Method (b):  3.40425468e-05
Method (c):  NaN

To explain this, we can look at the limits of single precision floating point numbers.

rMax = realmax('single')
fprintf('The largest number a "single" can hold is %e\n', rMax )
fprintf('log10 of this maximum number is: %f\n', log10(rMax) );

rMax =

  3.4028e+38

The largest number a "single" can hold is 3.402823e+38
log10 of this maximum number is: 38.531841

Hence, on iteration 39, we compute 10^39 and this value exceeds the maximum number, so we get overflow error (and the value is labelled "Inf"). Method (b) doesn't have overflow because, before we convert 10^39 to a single, we first divide by 39!, and this results in a small number.

However, method (b) has its own problems. The factorial grows faster than the exponential (roughly, x! grows like x^x -- think of Stirling's approximation). So, the quantity 10^x/x! gets very small, and in fact we have underflow! The quantity gets so small that we can't represent it as a single precision number.

rMin = realmin('single')
fprintf('The smallest number a "single" can hold is %e\n',rMin)
fprintf('When x = 87, the value of 10^x/x! is %e\n', 10^87/factorial(87) );
fprintf('   but if we store it as a single, Matlab treats it as: %e\n', single(10^87/factorial(87)) );
% It's unlikely that this underflow contributed to the loss of precision in
% (b), because errors from loss of relative precision crept in earlier.
% These errors are from the summing!

rMin =

  1.1755e-38

The smallest number a "single" can hold is 1.175494e-38
When x = 87, the value of 10^x/x! is 4.744379e-46
   but if we store it as a single, Matlab treats it as: 0.000000e+00

A "magic" fix to method (b) in Problem 1

Earlier, we summed from beginning to end, i.e. each new term that we added to the sum was smaller than the previous. This is the worst way to sum! If we instead sum from smaller to larger numbers, we can keep much more precision, at no extra cost. To flip the order of a vector in Matlab, use the command "fliplr" (the "lr" is for "left-right"), or "flipud"( the "ud" is for "up-down"), depending on if it is a row or column vector.

fprintf('Exact value: %.8e\n',exp(-10))
fprintf('Method (b), summing from large to small: %.3e\n', sum(v2) );
fprintf('Method (b), summing from small to large: %.3e\n', sum(fliplr(v2)) );
% We still don't have much accuracy, but we at least have 2 correct digits,
% instead of being off by a factor of 10.

Exact value: 4.53999298e-05
Method (b), summing from large to small: 3.404e-05
Method (b), summing from small to large: 3.404e-05

Problem 2: Fair dice

In general, we expect a fair die to have a mean value of 3.5. Of course, if we roll a die several thousand times, we won't find an empirical mean of exactly 3.5 -- we expect this to fluctuate a bit about the true mean (think of the Central Limit Theorem). If the two die have similar means that are barely different, this isn't evidence that one is fair. If one die had a mean of 5, then yes, this is evidence, but that's not the case. And of course, to be tricky, we ensured that the unfair die actually had a mean closer to 3.5 than the fair die.

load diceData
fprintf('The mean of die 1 is %.4f\n', mean(die1) )
fprintf('The mean of die 2 is %.4f\n', mean(die2) )
% Is a fluctuation of .02 from the true mean, after 5000 rolls, telling?
% We could answer precisely using a statistical test, but that's not
% necesssary.  Intuitively, even though it's a bit suspicious, it doesn't
% really raise the alarm.  We should find a better test.  A good test is to
% look at the frequency for which the numbers 1:6 show up:
figure(1);
subplot(1,2,1); hist(die1,6); title('Die 1'); xlim([1,6])
subplot(1,2,2); hist(die2,6); title('Die 2'); xlim([1,6])

The mean of die 1 is 3.5218
The mean of die 2 is 3.5048

Judging by the histogram, there is very clear evidence that die1 is relatively fair, and die2 is certainly unfair. Below is how we generated the two die:

rounding = @(x) floor(x)+1;
roundingFunction = @(x) min( max(rounding(x),1),6);
N = 5000;
die1 = roundingFunction( 6*rand(N,1)+.03  );     % a fair die
die2 = roundingFunction( 1.5*randn(N,1)+3 );     % loaded die
% The fair die was generated from a uniform distribution (and the mean was
% shifted a bit to account for the systematic errors from the rounding
% operation), while the unfair die was generated from a normal
% distribution.

Problem 3: RANDU

Taking a look at Moler's book or the wikipedia page in RANDU gives us a 3-term recurrence relation for a RANDU sequence.

Below, we plot many triplets x_k, x_{k+1} and x_{k+2}, and we can see graphically that these are not isotropically distributed. Basically, any kind of structure in psuedo-random numbers is unwanted (unless you're doing something like Quasi-Monte Carlo).

Cleve Moler's chapter in his online textbook talks about the old generators that Matlab used, and about the modern generators (which themselves are not perfect).

load randomNumbers
figure(1); clf;
subplot(1,2,1);
plot3( rand1(1:3:end), rand1(2:3:end), rand1(3:3:end ), '.',...
    'markersize',2)
title('Sequence 1: generated by RANDU'); axis square
subplot(1,2,2);
plot3( rand2(1:3:end), rand2(2:3:end), rand2(3:3:end ), '.',...
    'markersize',2)
title('Sequence 2: generated by Matlab''s "rand"'); axis square

Problem 4: Vectorization and Anonymous Functions

There are several ways you can write this. Don't forget to put periods before operations like multiplication and exponentiation. Here are two possible ways:

f_a = @(x) (x>=0).*(x<2).*(  sin(2*x).*(x<1) + exp(-x.^2).*(x>=1) );
f_b = @(x) (x>=0 & x<2).*(  sin(2*x).*(x<1) + exp(-x.^2).*(x>=1) );
x = [-.3, .3, 1.5, 2, 200];
disp( f_a(x) )
disp( f_b(x) )

         0    0.5646    0.1054         0         0

         0    0.5646    0.1054         0         0

Problem 5(a): Fourier Interpolation

Part (i). We find below that it takes only 10 points to acheive near perfect reconstruction. This means our 512-point sample can be compressed to only 10 points. This is not surprising, because our function f(x) is a trigonometric polynomial itself, and it only has 2 frequencies (or 4, if you count negatives). For more information, read about the Nyquist–Shannon sampling theorem. Note: not all of the plots shown below were needed for the homework; they are included just to explain what's going on.

a = 0; b = 2;       % the domain
f = @(x) sin( 4*pi*x ) .* sin(x*pi) .*(x>=a).*(x<b);
NN = 512;
dX = (b-a)/NN;
grid = a:dX:(b-dX);
signal = f(grid);

figure(1); clf
plot(grid,signal,'y-','linewidth',6); hold on;

Colors = get(gca,'colororder'); ci = 1; Labels = {'True signal'};
err = Inf; N = 0;
while err > 1e-14
    N = N + 1;
    dX = (b-a)/N;
    sampleGrid = a:dX:(b-dX);
    samples = f( sampleGrid );
    estimate = interpft(samples,512);
    err = norm(estimate-signal,'inf');
    fprintf('For N = %d, l_inf error is %.2e\n',N,err)

    % do some fancy stuff so that we get a new color every time we plot
    colr  = Colors(ci,:); ci = ci+1; if ci > size(Colors,1), ci = 1; end
    if N < 6
        plot(grid,estimate,':','color',colr);
    else
        plot(grid,estimate,'color',colr);
    end
    Labels{N+1} = sprintf('with N = %d',N');
end
legend(Labels,'location','bestoutside');
title('Reconstructions via Fourier Interpolation for various grid sizes');
fprintf('Problem 5(a) part i: it take %d points for near perfect reconstruction\n',N);

For N = 1, l_inf error is 9.29e-01
For N = 2, l_inf error is 9.29e-01
For N = 3, l_inf error is 1.38e+00
For N = 4, l_inf error is 9.29e-01
For N = 5, l_inf error is 1.81e+00
For N = 6, l_inf error is 9.08e-01
For N = 7, l_inf error is 1.00e+00
For N = 8, l_inf error is 9.29e-01
For N = 9, l_inf error is 1.00e+00
For N = 10, l_inf error is 2.23e-15
Problem 5(a) part i: it take 10 points for near perfect reconstruction

We see that the spectrum is very sparse. There's not much information in this signal!

figure(2);
F = fft(signal);
plot( abs(F(1:256+1)),'o:','markersize',4 );
title('Spectrum of the signal f(x)');
xlabel('frequency'); xlim([0,257])

Part (ii). With b = 2.2, we don't get such good results. Even though the signal still contains the exact same frequencies, the spectrum, via a discrete Fourier transform, will not look as sparse. This is because we have boundary effects. The Fourier transform works best on periodic functions -- in essence, it makes the function periodic -- so the fact that our signal f(x) is not zero at the boundary makes for discontinuities. When we look at the error in interpolation, we see that almost all the error is made at the boundary. The functions we use in interpolation -- Fourier polynomials on [0,2.2] -- all vanish at the endpoint, so we have a hard time reconstructing a function that doesn't vanish at the endpoint.

a = 0; b = 2.2;       % the domain
f = @(x) sin( 4*pi*x ) .* sin(x*pi);
NN = 512;
dX = (b-a)/NN;
grid = a:dX:(b-dX);
signal = f(grid);

figure(1); clf; subplot(1,2,1);
plot(grid,signal,'b-','linewidth',3); hold on;

N=50;
dX = (b-a)/N;
sampleGrid = a:dX:(b-dX);
samples = f( sampleGrid );
estimate = interpft(samples,512);
err = norm(estimate-signal,'inf');
fprintf('For N = %d, l_inf error is %.2e\n',N,err)
plot(grid,estimate,'r');
title('Actual signal and interpolation, b = 2.2');
xlim([0,2.2]);
legend('Actual signal','Interpolation with 50 points');

subplot(1,2,2);
plot(grid, abs(estimate-signal),'k*' )
xlim([0,2.2]); title('Error in interpolation');

For N = 50, l_inf error is 3.32e-01

Part (iii). We don't expect such dramatic results using normal polynomial interpolation because sines and cosines are not particularly well approximated by polynomials. Of course, as we increase the number of sample points, there is improvement. We show only N = 20 below.

a = 0; b = 2;       % the domain
f = @(x) sin( 4*pi*x ) .* sin(x*pi);
NN = 512;
dX = (b-a)/NN;
grid = a:dX:(b-dX);
signal = f(grid);
figure(1); clf;
plot(grid,signal,'b','linewidth',3); hold on

N = 20;
dX = (b-a)/N;
sampleGrid = a:dX:(b-dX);
samples = f(sampleGrid);

% Use spline interpolation
est_spline = spline(sampleGrid,samples,grid);
plot(grid,est_spline,'r')
err_spline = est_spline-signal;
fprintf('Error via spline is %.2e\n',norm(err_spline,'inf'))

% Use pchip interpolation
est_pchip = pchip(sampleGrid,samples,grid);
plot(grid,est_pchip,'g')
err_pchip = est_pchip-signal;
fprintf('Error via pchip is  %.2e\n',norm(err_pchip,'inf'))

legend('Actual signal','spline interpolation','pchip interpolation',...
    'location','southwest');
title('Polynomial Interpolation with N = 20');

Error via spline is 1.69e+00
Error via pchip is  1.02e-01

Both pchip and spline use cubic polynomials. They don't use a single cubic, but rather have many different cubic polynomials and combine them (e.g. they use a certain cubic for the [a,b] region, and another for the [b,c] region, etc.). We can verify this by looking at the third derivative of the interpolations. For a cubic, we expect this to be a constant. In Matlab, we can use the differences, using the "diff" command, to estimate the derivative.

figure(1); clf;
% the command diff(x,n) finds the nth-order discrete differences
plot( diff(est_pchip,3),'r'); hold on
plot( diff(est_spline,3),'b');
ylim(2e-3*[-1,1]);
legend('3rd derivative of pchip','3rd derivative of spline','location',...
    'southeast');
title('Derivatives of interpolations via various cubic methods');

Part (iv): Bonus -- padding with zeros. We keep the same setup for the signal. We can redefine f(x) such that f(x) is zero for x outside of [a,b]. This trick will make the zero-padding easier. As we see from the error results below, the zero-padding helps enormously, especially at the end-points.

f = @(x) sin( 4*pi*x ) .* sin(x*pi) .*(x>=a).*(x<b);

d = (a+b)/2;
extendedSampleGrid = (a-d):1/N:(b-1/N+d);
extendedSamples = f(extendedSampleGrid);

figure(1); clf;
plot(extendedSampleGrid, extendedSamples,'o:','markersize',4 );
hold on; plot(grid,signal,'r'); title('Zero-padded samples');
legend('Samples (with zeros at front and back)','Actual signal');
figure(2); clf; plot(grid,signal,'b','linewidth',3); hold on

% Use spline interpolation
est_spline = spline(extendedSampleGrid,extendedSamples,grid);
plot(grid,est_spline,'r')
err_spline = est_spline-signal;
fprintf('Error via zero-padded spline is %.2e\n',norm(err_spline,'inf'))

% Use pchip interpolation
est_pchip = pchip(extendedSampleGrid,extendedSamples,grid);
plot(grid,est_pchip,'g')
err_pchip = est_pchip-signal;
fprintf('Error via zero-padded pchip is  %.2e\n',norm(err_pchip,'inf'))

legend('Actual signal','spline interpolation','pchip interpolation',...
    'location','southwest');
title('Zero-padded Polynomial Interpolation with N = 20');

Error via zero-padded spline is 5.60e-03
Error via zero-padded pchip is  3.55e-02

Problem 5(b): Polynomial Interpolation

Part (i). The setup is very similar to problem 5(a); we just need to change our definition of f(x). We see that, for even just N = 4, the spline interpolation has virtually no error. This isn't too surprising, since our function f(x) is a polynomial with degree less than 3, so we only need 4 or fewer points to specify it.

a = 0; b = 1;
p = poly ( [.4, .7] );
fprintf('The polynomial looks like: %s\n',poly2str(p,'x') );
f = @(x) polyval(p,x);
NN = 512;
dX = (b-a)/NN;
grid = a:dX:(b-dX);
signal = f(grid);
figure(1); clf;
plot(grid,signal,'b','linewidth',3); hold on

N = 4;
dX = (b-a)/N;
sampleGrid = a:dX:(b-dX);
samples = f(sampleGrid);

% Use spline interpolation
est_spline = spline(sampleGrid,samples,grid);
plot(grid,est_spline,'r')
err_spline = est_spline-signal;
fprintf('Error via spline is %.2e\n',norm(err_spline,'inf'))

% Use pchip interpolation
est_pchip = pchip(sampleGrid,samples,grid);
plot(grid,est_pchip,'g')
err_pchip = est_pchip-signal;
fprintf('Error via pchip is  %.2e\n',norm(err_pchip,'inf'))
legend('Actual signal','spline interpolation','pchip interpolation',...
    'location','southwest');
title('Polynomial Interpolation of a quadratic polynomial with N = 4');

The polynomial looks like:    x^2 - 1.1 x + 0.28
Error via spline is 2.50e-16
Error via pchip is  4.90e-02

Part (ii). The setup is just like part (i). We don't expect to do as well with just 4 points, since now our actual signal is a quartic.

a = 0; b = 1;
p = poly( [.2, .3, .8, .88] );
fprintf('The polynomial looks like: %s\n',poly2str(p,'x') );
f = @(x) polyval(p,x);
NN = 512;
dX = (b-a)/NN;
grid = a:dX:(b-dX);
signal = f(grid);
figure(1); clf;
plot(grid,signal,'b','linewidth',3); hold on

N = 10;
dX = (b-a)/N;
sampleGrid = a:dX:(b-dX);
samples = f(sampleGrid);

% Use spline interpolation
est_spline = spline(sampleGrid,samples,grid);
plot(grid,est_spline,'r')
err_spline = est_spline-signal;
fprintf('Error via spline is %.2e\n',norm(err_spline,'inf'))

% Use pchip interpolation
est_pchip = pchip(sampleGrid,samples,grid);
plot(grid,est_pchip,'g')
err_pchip = est_pchip-signal;
fprintf('Error via pchip is  %.2e\n',norm(err_pchip,'inf'))
legend('Actual signal','spline interpolation','pchip interpolation',...
    'location','southwest');
title('Polynomial Interpolation of a quartic polynomial with N = 10');

The polynomial looks like:    x^4 - 2.18 x^3 + 1.604 x^2 - 0.4528 x + 0.04224
Error via spline is 1.81e-03
Error via pchip is  6.50e-03