Preliminary Examination, Numerical Analysis, Fall 1998

Linear Least Squares The Linear Least Squares problem for an \(m\times n\) real matrix \(A\) and \(b\in \mathbb {R}^m\) is the problem: \[ \text {Find }x\in \mathbb {R}^n\text { such that }\lVert Ax-b\rVert _2\text { is minimized.} \]

a) Suppose that you have data \(\{(t_j,y_j)\}\), \(j=1,2,\ldots ,m\) that you wish to approximate by an expansion \[ p(t)=\sum _{k=1}^n x_k\phi _k(t). \] Here, the functions \(\phi _k(t)\) are given functions. Which norm on the difference between the approximation function \(p\) and the data gives rise to a linear least squares problem for the unknown expansion coefficients \(x_k\)? What is the matrix \(A\) in this case, and what is the vector \(b\)?

b) Suppose that \(A\) is a real \(m\times n\) matrix of full rank and let \(b\in \mathbb {R}^m\). Show that the Least Squares problem has a unique solution and describe the geometry of the Least Squares problem (in terms of vectors and/or subspaces of \(\mathbb {R}^m\) and/or \(\mathbb {R}^n\)).

c) For the matrix \(A\) and vector \(b\) of part b), answer the following questions: What are the ‘normal equations’ for the Least Squares problem? How can they be used to solve the Least Squares problem? What is the QR factorization of \(A\) and how can it be used to solve the Least Squares problem? How can the SVD of \(A\) be used to solve the Least Squares problem? Compare and contrast these approaches for numerically solving the Least Square problem.

Singular Value Decomposition Suppose that \(A\) is a real \(m\times n\) matrix, \(m\ge n\).

a) What is the singular value decomposition (SVD) of \(A\)? Does every matrix \(A\) of this type have an SVD? If your answer is yes, prove the assertion. If your answer is no, give an example of a matrix \(A\) which has no SVD.

b) For a matrix \(A\) which does have an SVD, how would having the SVD give easy access to information about the rank of \(A\)? If \(A\) were \(m\times m\) and nonsingular, how would knowing \(A\)’s SVD help you determine the sensitivity of the solution of a linear system of the following form \[ Ax=b? \]

c) For a nonsingular matrix \(A\) which does have an SVD, express the relative distance from \(A\) to the nearest rank deficient matrix in terms of \(A\)’s singular values.

Iterative Methods for Linear Systems a) Consider the fixed-point iteration \[ u^{(k+1)}=Tu^{(k)}+c \] for finding a solution of the problem \[ u=Tu+c \] where \(T\) is an \(m\times m\) real matrix and \(c\) is a real \(m\)-vector. Show that the fixed point iteration will converge for an arbitrary initial guess \(u^{(0)}\) if and only if the spectral radius of \(T\), \(\rho (T)\), is less than 1.

b) Consider the boundary value problem \[ -u''(x)=f(x),\qquad \text {for }0\le x\le 1 \] with \(u(0)=u(1)=0\), and the following discretization of it: \[ -U_{j-1}+2U_j-U_{j+1}=F_j \] for \(j=1,2,\ldots ,N-1\) where \(Nh=1\) and \(F_j=h^2f(jh)\).

Analyze the convergence properties of the Jacobi iterative method for this problem. In particular, express the speed of convergence as a function of the discretization stepsize \(h\). How does the number of iterations required to reduce the initial error by a factor of \(\delta \) depend on \(h\)? In practice, would you use this method to solve the given problem? If so, explain why this is a good idea? If not, how would you solve it in practice?

Eigenvalue Problems

a) Consider the problem of computing an eigenvalue/eigenvector pair of a real symmetric \(m\times m\) matrix \(A\) with eigenvalues of distinct magnitudes. Two algorithms for trying to do this are Power Iteration and Rayleigh Quotient Iteration:

i) Analyze the convergence properties of the Power Iteration algorithm for a generic guess \(v^0\).

ii) Suppose that the Rayleigh Quotient Iteration iterates \(\lambda ^{(k)}\) converge to a simple eigenvalue \(\lambda \) of \(A\). Analyze the speed of convergence of \(\lambda ^{(k)}\) to \(\lambda \).

b) Many algorithms for finding eigenvalues of a complex \(m\times m\) matrix \(A\) are based on the existence of a Schur decomposition of \(A\), namely, \(A=QTQ^*\) where \(Q\) is unitary and \(T\) is upper triangular. Prove that every matrix \(A\) has a Schur decomposition and explain how knowing the Schur decomposition of \(A\) would help one determine the eigenvalues of \(A\).

Krylov Subspace Methods Consider the system \(Ax=b\), in which \(A\) is a symmetric positive definite real \(m\times m\) matrix, and denote by \(x_*\) its true solution. The conjugate gradient (CG) algorithm for finding the solution to the system is:

Set \(x_0=0\), \(r_0=b\), and \(p_0=r_0\).

For \(n=1,2,3,\ldots \) \[ \begin {aligned} \alpha _n &= r_{n-1}^T r_{n-1}/p_{n-1}^T A p_{n-1},\\ x_n &= x_{n-1}+\alpha _n p_{n-1},\\ r_n &= r_{n-1}-\alpha _n A p_{n-1},\\ \beta _n &= r_n^T r_n/r_{n-1}^T r_{n-1},\\ p_n &= r_{n-1}+\beta _n p_{n-1}. \end {aligned} \] End

The vectors generated by this algorithm have the property that, provided \(r_{n-1}\ne 0\), \[ \begin {aligned} \mathcal {K}_n &\equiv \operatorname {span}\{b,Ab,A^2b,\ldots ,A^{n-1}b\} =\operatorname {span}\{x_1,x_2,\ldots ,x_n\}\\ &=\operatorname {span}\{r_0,r_1,\ldots ,r_{n-1}\} =\operatorname {span}\{p_0,p_1,\ldots ,p_{n-1}\} \end {aligned} \] and \(r_n^T r_j=0\) and \(p_n^T A p_j=0\) for all \(j<n\).

Let \(e_n=\lVert x_n-x_*\rVert _A\), where the \(A\)-norm of a vector \(x\) is defined by \(\lVert x\rVert _A=(x^TAx)^{1/2}\).

a) Prove the following:

Theorem: Suppose the CG algorithm is applied to the symmetric positive definite system \(Ax=b\) and that it has not yet converged (\(r_{n-1}\ne 0\)). Then \(e_n<\lVert y-x_*\rVert _A\) for all \(y\in \mathcal {K}_n\), \(y\ne x_n\); \(\lVert e_n\rVert _A\le \lVert e_{n-1}\rVert _A\); and \(e_N=0\) for some \(N\le m\).

b) In practice this result is not very useful. Why is this, and what determines the speed of convergence of \(x_n\) to \(x_*\)? What is ‘preconditioning’ and how is it related to the speed of convergence of the CG algorithm?

Numerical Methods for ODEs Consider the Linear Multistep Method \[ y_{n+1}-y_n=\frac {1}{2}h(f_{n+1}+f_n) \] for solving an initial value problem \(y'=f(y,x)\), \(y(0)=\eta \). You may assume that \(f\) is Lipschitz continuous with respect to \(y\) for all \(x\).

a) Analyze the consistency, stability, accuracy, and convergence properties of this method.

b) Would it be more reasonable to use this method or Euler’s method for the initial value problem: \[ y'=-10^8[y-\cos (x)]-\sin (x),\qquad y(0)=1? \] Justify your answer.

Finite-Difference Methods Consider the initial boundary value problem \[ v_t=\beta v_{xx},\qquad 0<x<1,\qquad 0<t,\qquad v(x,0)=f(x) \] with Dirichlet boundary conditions \[ v(0,t)=0,\qquad v(1,t)=0. \]

Suppose that the interval \([0,1]\) is partitioned by meshpoints \(x_j=jh\) for \(j=0,1,\ldots ,N\) and \(h=1/N\). A possible discretization of the above problem uses the Crank-Nicolson scheme for the differential equation: \[ \left (I-\frac {\beta k}{2}D_+D_-\right )u_j^{n+1} = \left (I+\frac {\beta k}{2}D_+D_-\right )u_j^n \] for \(j=1,\ldots ,N-1\), combined with the boundary conditions: \[ u_0^n=0\qquad \text {and}\qquad u_N^n=0 \] for each time level \(n\). Initial data is obtained from by setting \(u_j^0=f_j\) for \(j=0,1,\ldots ,N\). (Note that \(D_+\) and \(D_-\) are the usual forward and backward difference operators.)

Set up the matrix system that must be solved in order to advance the numerical solution from time level \(n\) to time level \(n+1\). Show that the suggested discretization is stable in the \(l_2\) norm \[ \lVert u\rVert _2= \left (\sum _{j=0}^N h\lvert u_j\rvert ^2\right )^{1/2} \] for all choices of timestep \(k\) and spacestep \(h\). What other issues should be considered in deciding whether the suggested scheme is a good choice for the problem?