Preliminary Examination, Numerical Analysis, Fall 2001

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? Compare and contrast these approaches for numerically solving the Least Square problem.

Singular Value Decomposition The singular value decomposition (SVD) of a real \(m\times n\) matrix \(A\) (you may assume \(m\ge n\)) is defined by \[ A=U\Sigma V^T \] where \(U\) is \(m\times m\) and orthogonal, \(\Sigma \) is \(m\times n\) and diagonal, \(V\) is \(n\times n\) and orthogonal.

a) Prove that every matrix \(A\) has an SVD.

b) Use the SVD of \(A\) to prove the famous result from linear algebra that the null space of \(A\) is the orthogonal complement of the range of \(A^T\) (the adjoint of \(A\)).

c) For a nonsingular \(m\times m\) matrix \(A\), the solution to the system \(Ax=b\) can be expressed in terms of \(A\)’s SVD. Write \(b=\sum _{j=1}^m c_j u_j\) where \(u_j\) is the \(j^{\text {th}}\) column of \(U\). Then the solution \(x\) is \[ x=\sum _{j=1}^m \frac {c_j}{\sigma _j}v_j \] where \(v_j\) is the \(j^{\text {th}}\) column of \(V\) and \(\sigma _j\) is the \(j^{\text {th}}\) singular value. Verify this formula. Consider the special case \(b=u_1\). What perturbation \(\delta b\) in \(b\), with \(\lVert \delta b\rVert _2=1\), causes the largest change \(\delta x\) in the solution \(x\)?

Iterative Methods for Linear Systems 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.

a) 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) Let \(A\) be a real \(n\times n\) matrix with simple eigenvalue \(\lambda \), right eigenvector \(x\) (i.e., \(Ax=\lambda x\)), and left eigenvector \(y\) (i.e., \(y^*A=\lambda y^*\)). Assume that \(\lVert x\rVert _2=\lVert y\rVert _2=1\). Show that the condition number for \(\lambda \) to changes in \(A\) is given by \(s(\lambda )=1/\lvert y^*x\rvert \). Give an example of a matrix \(A\) which has a badly conditioned eigenvalue, and another matrix \(A\) which has only well-conditioned eigenvalues.

b) Consider the problem of computing an eigenvalue/eigenvector pair of a real symmetric \(m\times m\) matrix \(A\) with eigenvalues of distinct magnitudes. One algorithm for trying to do this is Rayleigh Quotient Iteration:

Guess \(v^{(0)}\) with \(\lVert v^{(0)}\rVert _2=1\).

Set \(\lambda ^{(0)}=(v^{(0)})^TAv^{(0)}\).

For \(k=1,2,3,\ldots \) \[ \begin {aligned} &\text {Solve }(A-\lambda ^{(k-1)}I)w=v^{(k-1)}\text { for }w,\\ &v^{(k)}=w/\lVert w\rVert _2,\\ &\lambda ^{(k)}=(v^{(k)})^TAv^{(k)}. \end {aligned} \] End

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 \).

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:

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+2}-\frac {4}{3}y_{n+1}+\frac {1}{3}y_n=\frac {2}{3}hf_{n+2} \] 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.