Scientific Computing, Spring 1999 (http://www.math.nyu.edu/faculty/goodman/teaching/sci comp 99/)

Assignment 3.

Given February 9, due February 23.

1. Write a program that computes the Choleski factorization of a symmetric positive definite pentadiagonal matrix. The program should return an error message if it is unable to do this, for example, if the matrix is not positive definite or a parameter is out of range.

2. Write a program that takes a lower triangular matrix, L, and computes LL^{t}. Use this to check the L computed by 1.

3. Write a program to do forward and back substitution to solve LL^{t}u = f. It should use an externally supplied L. Write a tester that multiplies out LL^{t}u, given L and u.

4. Write a procedure that takes as input x, x_{0}, Δx, and computes the value at x of the symmertic cubic b-spline basis function centered at x_{0} with uniform knot spacing Δx. The answer is b(x, x_{0}, Δx). Use finite differences (with a step size much smaller than Δx) to check that b is correct. That is, that b, together with its first and second derivatives is continuous, and that the third derivative is constant between knots.

5. We are trying to reconstruct a function y(x) from noisy data. The sample points are uniformly spaced: x_{k} = k, k = 1, …, n. the approximating function will be a cubic spline with uniformly spaced knots at r_{j} = s · j. Here s is a positive integer, the knot spacing. The approximating curve will have the form

ŷ(x) = ∑_{j=0}^{n/s} u_{j}b(x, r_{j}, s) .

Statisticians use the hat to indicate a hat to indicate a statistical estimator, so ŷ is a statistical estimate of the unknown function y(x). Write a program that takes as input the data^{1} Y = (Y_{0}, …, Y_{n}), and the parameters n and s, and computes the u_{j} to minimize

R^{2} = ∑_{k=0}^{n} (ŷ(x_{k}) − Y_{k})^{2} .

As one check, note that ŷ should interpolate the data if s = 1. Experiment with other powers of 2 for s and with the three data sets (which have n = 1024). Try to get an s that gives an accurate picture of function from the noisy data set. Using least squares splines in this way is one method for smoothing noisy data.

1Statisticians often (but not always) use capital letters to indicate sample data that include random noise.

6. It is hard to know in advance what s to use. Taking s = 1 fits the data best but has little predictive power for noisy data. We will choose s using a simple form of cross validation. The idea of cross validation is to fit the model with part of the data and see how well it predicts the other part. In our situation, we can fit the model by minimizing the “in sample” sum of squares,

RIS^{2} = ∑_{k=0}^{n/2} (ŷ(x_{2k}) − Y_{2k})^{2} ,

and then checking the predictive power by looking the “out of sample” sum of squares,

ROS^{2} = ∑_{k=0}^{n/2−1} (ŷ(x_{2k+1}) − Y_{2k+1})^{2} .

The in sample sum of squares will always be smaller for smaller s, but the out of sample sum of squares behaves differently. Use this to select a good s value for each of the three data sets.