Preliminary Examination, Numerical Analysis, January 2005

Sensitivity Analysis In this problem, we consider a real, non-singular \(n\times n\) matrix \(A\) and vectors \(b,x\in \mathbb {R}^n\). We are concerned with the error in solving the linear system \(Ax=b\). It turns out that the numerical solution \(\hat {x}\) is only an approximate solution to this equation. However, it satisfies a perturbed equation \(\hat {A}\hat {x}=b\) exactly, where \(\hat {A}=A+\delta A\). Let the error be defined as \(\delta x=\hat {x}-x\); we are interested in bounding \(\lVert \delta x\rVert \) in terms of \(\lVert \delta A\rVert \lVert A^{-1}\rVert \), or \(\lVert \delta A\rVert \lVert \hat {A}^{-1}\rVert \), and \(\lVert x\rVert \), or \(\lVert \hat {x}\rVert \), where \(\lVert \cdot \rVert \) is some suitable norm, and we proceed as follows:

a) First prove the identity \[ A^{-1}-\hat {A}^{-1}=A^{-1}\cdot \delta A\cdot \hat {A}^{-1}, \] and hence deduce that \[ \lVert A^{-1}-\hat {A}^{-1}\rVert \le \lVert A^{-1}\rVert \cdot \lVert \delta A\rVert \cdot \lVert \hat {A}^{-1}\rVert . \]

b) Assuming that \(\hat {\delta }=\lVert \delta A\rVert \lVert \hat {A}^{-1}\rVert <1\), show that \[ \lVert A^{-1}\rVert \le \frac {1}{1-\hat {\delta }}\lVert \hat {A}^{-1}\rVert ,\qquad \lVert A^{-1}-\hat {A}^{-1}\rVert \le \frac {\hat {\delta }}{1-\hat {\delta }}\lVert \hat {A}^{-1}\rVert . \]

c) By comparing the equations \(Ax=b\) and \((A+\delta A)(x+\delta x)=b\), show that \[ \lVert \delta x\rVert \le \frac {\hat {\delta }}{1-\hat {\delta }}\lVert \hat {x}\rVert . \] Similarly, show that if \(\delta =\lVert \delta A\rVert \lVert A^{-1}\rVert <1\), \[ \lVert \delta x\rVert \le \frac {\delta }{1-\delta }\lVert x\rVert . \]

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:

Find \(x\in \mathbb {R}^n\) such that \(\lVert Ax-b\rVert _2\) 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\)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\)th column of \(V\) and \(\sigma _j\) is the \(j\)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\)?

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 \): \[ \text {Solve }(A-\lambda ^{(k-1)}I)w=v^{(k-1)}\text { for }w, \] \[ v^{(k)}=\frac {w}{\lVert w\rVert _2}, \] \[ \lambda ^{(k)}=(v^{(k)})^TAv^{(k)}. \] 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 \).

Iterative Methods for Linear Systems 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 \), \[ \alpha _n=\frac {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=\frac {r_n^T r_n}{r_{n-1}^T r_{n-1}}, \] \[ p_n=r_n+\beta _n p_{n-1}. \] End.

The vectors generated by this algorithm have the property that, provided \(r_{n-1}\ne 0\), \[ \begin {aligned} \mathcal {K}_n &= \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^T Ax)^{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\to 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 ODEs

a) The difference equation \[ y_{n+k}+a_1y_{n+k-1}+\cdots +a_ky_n=0 \] can, after the substitution \(\mathbf {y}_j=(y_j,y_{j+1},\ldots ,y_{j+k-1})^T\), be written in the form \(\mathbf {y}_{n+1}=A\mathbf {y}_n\). What does the matrix \(A\) look like? Describe the condition for \(y\) to remain bounded in terms of \(A\).

b) The differential equation \[ y^{(k)}+a_1y^{(k-1)}+\cdots +a_ky=0 \] can, after the substitution \(\mathbf {y}=(y,y',\ldots ,y^{(k-1)})^T\), be written in the form \(\mathbf {y}'=A\mathbf {y}\). What does the matrix \(A\) look like? Describe the condition for \(y\) to remain bounded in terms of \(A\).

c) Show that all solutions of the difference equation \[ y_{n+1}-2\lambda y_n+y_{n-1}=0 \] are bounded, as \(n\to \infty \), if \(\lambda \) is real and \(\lvert \lambda \rvert <1\).

Finite-Difference Methods Consider the initial boundary value problem \[ v_t=\beta v_{xx},\qquad 0<x<1,\quad t>0,\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 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 \(\ell _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\). (Note that von Neumann analysis is not enough to prove stability for the initial boundary value problem.) What other issues are important to consider in deciding whether this scheme is a good choice for the problem?

Fact 1: A real symmetric \(n\times n\) matrix \(A\) can be diagonalized by an orthogonal similarity transformation, and \(A\)’s eigenvalues are real.

Fact 2: The \((N-1)\times (N-1)\) matrix \(M\) defined by \[ M_{ij}= \begin {cases} -2, & i=j,\\ 1, & \lvert i-j\rvert =1,\\ 0, & \text {otherwise} \end {cases} \qquad (1\le i,j\le N-1) \] has eigenvalues \[ \mu _l=-4\sin ^2\left (\frac {\pi l}{2N}\right ),\qquad l=1,2,\ldots ,N-1. \]

Fact 3: The \((N+1)\times (N+1)\) matrix \(M\) defined by \[ M_{ij}= \begin {cases} -1, & i=j\in \{0,N\},\\ -2, & i=j,\ 1\le i\le N-1,\\ 1, & \lvert i-j\rvert =1,\\ 0, & \text {otherwise} \end {cases} \qquad (0\le i,j\le N) \] has eigenvalues \[ \mu _l=-4\sin ^2\left (\frac {\pi l}{2(N+1)}\right ),\qquad l=0,1,\ldots ,N. \]

Fact 4: For a real \(n\times n\) matrix \(A\), the Rayleigh quotient of a vector \(x\in \mathbb {R}^n\) is the scalar \[ r(x)=\frac {x^T Ax}{x^T x}. \] The gradient of \(r(x)\) is \[ \nabla r(x)=\frac {2}{x^T x}(Ax-r(x)x). \] If \(x\) is an eigenvector of \(A\) then \(r(x)\) is the corresponding eigenvalue and \(\nabla r(x)=0\).