Computational Algorithms for Fitting a Regression

Write a program to implement the Gram-Schmidt algorithm (Algorithm 11.3 in Section 11.3). Try out your program on the Longley data and the polynomial approximations described in Example 11.3 in Section 11.8.4. Compare your results with those in Examples 11.2 and 11.3.

Verify the formulas for the updating algorithm in Section 11.7.

Try some experiments to compare the execution time of a regression performed using the Cholesky method with one using the SVD. Do the flop counts given in the text give a reasonable indication of the relative speed of the two methods?

Calculate algebraically the solution to the system of linear equations corresponding to the matrix

using Gaussian elimination with and without partial pivoting. Show that if partial pivoting is not used, the solution $x_{3}$ is computed as

and as $1-2 \epsilon$ if partial pivoting is used. Comment on the accuracy of these two methods if $\epsilon$ is very small.

Show that applying one Gaussian elimination step to the matrix $\mathbf{X}^{T} \mathbf{X}$ yields the matrix

where $\widetilde{\mathbf{X}}$ is the centered version of $\mathbf{X}$.

Verify (11.15) and hence show that $\sqrt{\mathrm{RSS}}$ is given by the $(p+1, p+1)$ th element of $\mathbf{R}_{A}$.

If $\mathbf{R}^{T} \mathbf{R}$ is the Cholesky decomposition of $\mathbf{X}^{T} \mathbf{X}$, show that

Let $\mathbf{q}_{1}, \ldots , \mathbf{q}_{p}$ be the vectors constructed in Algorithm 11.3. Show that $\mathbf{q}_{1}, \ldots , \mathbf{q}_{j}$ is an orthonormal basis for $\mathcal{C}\left(\mathbf{a}^{(1)}, \ldots , \mathbf{a}^{(j)}\right)$ for $j = 1,2, \ldots , p$.

Show that an upper triangular matrix with positive diagonal elements is nonsingular. Also, prove that the inverse of an upper triangular matrix with unit diagonal elements is also upper triangular with unit diagonal elements.

Derive formulas (11.23), (11.24), and (11.25).

Prove that at the beginning of stage $j$ of the MGSA, the columns $\mathbf{a}^{(1)}, \ldots , \mathbf{a}^{(j)}$ are orthogonal, are orthogonal to $\mathbf{a}^{(j+1)}, \ldots , \mathbf{a}^{(p)}$, and span the first $j$ columns of the original matrix.

Show that a product of orthogonal matrices is orthogonal.

Show that if $\mathbf{Q}_{p}, \mathbf{Q}_{n-p}$ and $\mathbf{r}_{1}$ are as defined in (11.36), the fitted values $\widehat{\mathbf{Y}}$ are given by

Let $\mathbf{G}_{i h}$ be given by $(11.40)$ and $\widetilde{\mathbf{D}} = \operatorname {diag}\left(\tilde{d}_{1}, \ldots , \tilde{d}_{n}\right)$. Show that there is a diagonal matrix $\mathbf{D} = \operatorname {diag}\left(d_{1}, \ldots , d_{n}\right)$ and a fast rotator of the form (11.45) such that $\mathbf{G}_{i h} = \mathbf{D M} \tilde{\mathbf{D}}^{-1}$. Show that

where $t = s / c$ and $r = \tilde{d}_{h} / \tilde{d}_{i}$.

Prove that the condition number of $\mathbf{X}_{w}$ is given by (11.61).

Repeat the calculations in Example 11.3 using the function $f(t) = \cos 4 t$ in place of $\exp (\sin 4 t)$. Does the ranking of the methods still hold good? Sample R-code is given after Exercise 3 below.

Modify the code given below to count the number of flops that are required to execute each algorithm. How accurate are the formulas given in Table 11.2? ####################################################### # MGSA function: calculates Q and R^{-1} of QR decomp mgsa<-function(A){ n<-dim(A)[1] p<-dim(A)[2] GG<-diag(p) for(i in 1:(p-1)){ a<-A[,i] denom<-sum(aa) G<-numeric(p) for(j in (i+1):p){ num<-sum(A[,j]a) G[j]<- -num/denom A[,j]<-A[,j]+aG[j] } GG<-GG + outer(GG[,i],G) } list(W=A,G=GG) } ####################################################### #Householder function: calculates R of QR decomp of A house<-function(A){ n<-dim(A)[1] p<-dim(A)[2] for(j in 1:p){ indices<-j:n # construct Householder vector u<-numeric(n) u[indices]<-A[indices,j] unorm1<-sqrt(sum(u[indices]^2)) unorm2<-sum(u[indices[-1]]^2) uu<-u[j] u[j]<-if(uu<0)uu-unorm1 else -unorm2/(uu+unorm1) gamma<-0.5(unorm2 + u[j]^2) # multiply by householder matrix A[indices,j]<-0 A[j,j]<-unorm1 if(j==p)return(A) for(l in (j+1):p){ k<-sum(A[indices,l]u[indices])/gamma A[indices,l]<-A[indices,l] - ku[indices] } } A } ####################################################### # Givens function: calculates R of QR decomp of A givens<-function(A){ n<-dim(A)[1] p<-dim(A)[2] for(j in 1:(p-1)){ indices<-j:p for(i in (j+1):n){ sqr<-sqrt(A[j,j]^2+A[i,j]^2) cc<-A[j,j]/sqr ss<-A[i,j]/sqr temp<-ccA[j,indices] + ssA[i,indices] A[i,indices]<- -ssA[j,indices] + ccA[i,indices] A[j,indices]<-temp } } A }

Suppose that we fit a least squares regression to a subset $J$ of the $n$ available cases. Let $\mathbf{X}_{J}$ be the submatrix of $\mathbf{X}$ corresponding to the cases in $J$, let $e_{i}$ be the $i$ th residual from this fit, $i = 1, \ldots , n$, and let $h_{i j} = \mathbf{x}_{i}^{T}\left(\mathbf{X}_{J}^{T} \mathbf{X}_{J}\right)^{-1} \mathbf{x}_{j}$. Show that if we drop case $i \in J$ and add case $j \notin J$, the change in the residual sum of squares is

Using the result in Exercise 1 above, write an R function implementing Hawkins's method of calculating an approximation to the LTS estimate.

Write an R function to implement Ruppert's algorithm to calculate an approximation to the LTS estimate. You will need to refer to Ruppert [1992] for the details of the algorithm. Devise a small simulation to compare Ruppert's algorithm to Hawkins's.