Past exam of the numerical linear algebra and approximation course of the University of Oxford

Original solutions to all four questions and all thirteen terminal parts, including the two inline definitions in 1(a), of the Numerical Linear Algebra and Approximation paper sat on 3 June 2008. Source: released Mathematics papers, Part C 2007 11.pdf.

A real square orthogonal matrix $Q$ satisfies $\boxed{Q^{T}Q=Q{Q}^T=I}$. Equivalently its columns form an orthonormal basis. In particular it preserves the Euclidean norm, since $\Vert Qx{\Vert }_2^2=x^T{Q}^{T}Qx=\Vert x{\Vert }_2^2$.

A Householder matrix has the formIt reflects in the hyperplane perpendicular to $v$. Put $P=v{v}^T/(v^{T}v)$, an orthogonal projector satisfying $P^T=P$ and $P^2=P$. Then $H^{T}H=(I-2P{)}^2=I-4P+4P^2=I$, proving that every Householder matrix is orthogonal.

If $x\ne 0$, choose $\alpha =-\mathrm{sgn}(x_1)\Vert x{\Vert }_2$, with $\mathrm{sgn}(0)=1$ for this choice, and put $v=x-\alpha e_1$. This nonzero vector satisfiesThus the Householder matrix $H=I-2v{v}^T/(v^{T}v)$ gives $Hx=x-v=\alpha e_1$. Symmetry gives $\boxed{x^{T}H=\alpha e_1^T}$. If $x=0$, any Householder matrix works with $\alpha =0$.

For an LQ decomposition, process rows $j=1,\dots ,\min (m,n)$. At step $j$, postmultiply by an orthogonal block matrix $H_j=\mathrm{diag}(I_{j-1},{\hat{H}}_j)$ that zeros entries to the right of the diagonal in the current row. If its row tail is already suitable, use the identity. Earlier processed rows have zero entries throughout this trailing block, so their zeros are preserved. Consequently $L=A{H}_1\cdots H_r$ has $l_{ij}=0$ for $i<j$. Setting $\boxed{Q=H_r^T\cdots H_1^T}$ gives $\boxed{A=LQ}$. The argument handles both rectangular shapes; $L$ is lower trapezoidal when it is not square.

The singular value decomposition $A=U\Sigma V^T$ givesTherefore $\boxed{{\sigma }_i(A)=\sqrt{{\lambda }_i(A^{T}A)}}$, counting multiplicities. These are the nonnegative square roots: zero singular values must be included if $A$ is singular. Strict positivity requires full rank, which the question does not assume.

Write the lower bidiagonal matrix entries as $B_{ii}=d_i$ and $B_{i+1,i}=e_i$. Its columns $i$ and $j$ have disjoint supports if $|i-j|>1$, so $T=B^{T}B$ has no entries outside the three central diagonals. More explicitly,This proves both symmetry and the tridiagonal matrix structure.

A generic lower triangular matrix has a dense first column. The first left Householder transformation in its QR decomposition mixes those rows, usually creating entries throughout the rows to the right. Subsequent eliminations impose upper triangularity but no further systematic zeros; thus $R$ is generically a full upper triangular matrix.

For Householder bidiagonalization use two-sided operations instead. At stage $j$, a right Householder transformation acting on columns $j,\dots ,n$ clears the current row $j$ beyond its diagonal. Then a left Householder transformation acting on rows $j+1,\dots ,n$ clears column $j$ below its first subdiagonal. The left operation leaves row $j$ unchanged; both operations preserve the already completed lower-bidiagonal rows and columns. Continue for $j=1,\dots ,n-1$, omitting a final unnecessary left operation. The result isIt follows that $B^{T}B=Q^T{A}^{T}AQ$. Hence the two symmetric matrices are orthogonally similar and have the same eigenvalues. By part (c), compute the eigenvalues of the symmetric tridiagonal matrix $B^{T}B$ and take their nonnegative square roots to obtain all singular values of $A$.

For a tridiagonal matrix with subdiagonal $a_i$, diagonal $b_i$ and superdiagonal $c_i$, banded Gaussian elimination introduces no fill beyond the original band. The Thomas algorithm usesUpdate the right side by $d_1=f_1$, $d_i=f_i-{\ell }_i{d}_{i-1}$. Then back substitution givesEach row requires only a fixed number of arithmetic operations, so both work and storage are $\boxed{O(n)}$. The assumed absence of a need for pivoting includes nonzero elimination pivots $p_i$.

Write $G=M^{-1}N$. Subtracting the exact equation $Mx=Nx+b$ from the iteration gives $e_{\ell +1}=G{e}_{\ell }$, where $e_{\ell }=x-x^{(\ell )}$. Hence $e_{\ell }=G^{\ell }{e}_0$. Alsoso the stated scalars are precisely the eigenvalues of the iteration matrix $G$.

The spectral radius criterion for matrix powers states that $G^{\ell }\to 0$ if and only if all its eigenvalues have modulus strictly less than one. Applying that theorem gives $e_{\ell }\to 0$ for every initial vector. Therefore $\boxed{x^{(\ell )}\to A^{-1}b}$.

Order points along each $j$-line before moving to the next $k$-line. The splitting hasand $N$ couples the two neighboring $k$-lines. Factor $T_n$ once with the Thomas algorithm; each iteration then needs $n$ tridiagonal solves and the neighboring-line right sides, costing $\boxed{O(n^2)}$ for $n^2$ unknowns. Known boundary values are moved to the right side; the error eigenmodes have homogeneous boundary values.

The generalized eigenvalue problem iswith $U=0$ at the boundary indices. Put ${\theta }_r=r\pi /(n+1)$, ${\theta }_s=s\pi /(n+1)$. The sine addition formula gives $U_{j+1,k}+U_{j-1,k}=2\cos {\theta }_r{U}_{j,k}$ and the corresponding identity in $k$. ThusThese tensor-product sine vectors form a basis, so they give the whole spectrum. With $c=\cos (\pi /(n+1))\in [0,1)$,The iteration therefore converges.

It is not a uniformly effective multigrid smoother. Take $r=1$ and $s=n$: the error is slowly varying along a line but highly oscillatory between lines, and its multiplier is $-c/(2-c)\to -1$ as the grid is refined. It changes sign at each step while retaining almost its full amplitude. This line-Jacobi Fourier spectrum distinguishes convergence on each fixed grid from the high-frequency damping needed by the multigrid method.

Let $e_j=x-x^{(j)}$. The underlying Richardson iteration gives $e_j=(I-A{)}^j{e}_0$. The polynomials $1,(1-t),\dots ,(1-t{)}^k$ form a basis for polynomials of degree at most $k$. Thus writeDefine $y^{(k)}={\sum }_{j=0}^k{c}_j{x}^{(j)}$. Since the coefficients sum to one,The spectral theorem for real symmetric matrices supplies an orthonormal basis of eigenvectors. In that basis $p_k(A)$ has diagonal entries $p_k({\lambda }_i)$, soThe algebraic recombination is valid in exact arithmetic even if the original unaccelerated iterates themselves fail to converge.

Map $[a,b]$ to $[-1,1]$ by $z(r)=(2r-a-b)/(b-a)$. Since $a>0$, the normalization point $z(0)=-(a+b)/(b-a)$ lies outside that interval. The proposed Chebyshev polynomial ratio has degree $k$ and equals one at zero. On the spectral interval its maximum absolute value isThe Chebyshev minimax polynomial with exterior normalization minimizes this maximum among degree-at-most-$k$ polynomials with $p(0)=1$. One may quote the minimax property. Its reason is the alternating $k+1$ extremal values: a strictly better competitor would force their difference to have at least $k$ roots in $[a,b]$ and the additional root at zero outside it, impossible for degree at most $k$.

The paper uses a monic normalization of $T_k$ for $k\ge 1$, rather than the usual $T_k(1)=1$ convention. The common factor cancels in the ratio. With $\kappa =b/a>1$ and $q=(\sqrt{\kappa }-1)/(\sqrt{\kappa }+1)$, the hyperbolic-cosine formula gives the error factor $\boxed{2q^k/(1+q^{2k})\le 2q^k}$. This tends rapidly to zero when the spectral interval is well-conditioned. If $a=b$, use $p(r)=1-r/a$ for exact one-step termination instead of the displayed interval map.

For $[a,b]=[1,2]$, $z(0)=-3$ and $\cosh \theta =3$ gives $e^{\theta }=3+2\sqrt{2}$. The normalization factors cancel, soIf the eigenvalues are known to be exactly $1,3/2,2$, use the more targeted annihilating polynomialIt equals one at zero and vanishes on the entire spectrum, so the spectral theorem gives $p_3(A)=0$. Thus $y^{(3)}=x$ in exact arithmetic. Degree three is the least degree that can vanish at all three distinct eigenvalues; a particular initial error supported on fewer eigenspaces can terminate earlier.

The Krylov subspace is $\boxed{K_k(A,b)=\mathrm{span}\{b,Ab,\dots ,A^{k-1}b\}}$. If $y\in K_k$, write $y=q(A)b$ with $\mathrm{deg}q\le k-1$. Since $b=Ax$,which has degree at most $k$ and $p(0)=1$. Conversely, if $x-y=p(A)x$ with these properties, then $q(t)=(1-p(t))/t$ is a polynomial of degree at most $k-1$ andThis is the polynomial description of the admissible errors in a Krylov subspace method.

For a symmetric positive-definite matrix, the conjugate gradient method chooses the unique point in $x_0+K_k(A,r_0)$ minimizing the energy norm of the error. Equivalently it minimizes $\frac{1}{2}{y}^{T}Ay-b^{T}y$ on that affine space. Here $x_0=0$ and $r_0=b$.

Choose the eigenbasis orthonormally, as the spectral theorem permits, and write $x={\sum }_i{\alpha }_i{v}_i$. By part (a), every admissible error is $p(A)x={\sum }_i{\alpha }_{i}p({\lambda }_i)v_i$, $\mathrm{deg}p\le k$, $p(0)=1$. Orthogonality then givesThus the implicit CG polynomial minimizes exactly the displayed weighted sum. Strictly, the printed hypothesis that the eigenvectors are individually normalized must be supplemented by mutual orthogonality when eigenvalues repeat. Otherwise cross terms within a repeated-eigenvalue eigenspace need not vanish. An orthonormal eigenbasis always supplies the intended formulation.

If at most $\ell$ eigenvectors contribute to $x$, a polynomial with one factor $1-t/\lambda$ for each distinct contributing eigenvalue is admissible by step $\ell$ and makes that error zero. CG minimizes a nonnegative quantity, so $\boxed{\text{exact termination occurs in at most }\ell \text{ steps}}$.

For the last question, write $x=\alpha v_1+\beta v_2$ with eigenvalues $1,L$, $L=10^6$. The first polynomial has the form $p_1(t)=1-\eta t$. Differentiating its weighted objective givesIf the high-eigenvalue component is not exceptionally small, the $L$-weighted terms dominate and $\eta \approx 1/L$, so $p_1(1)\approx 1-10^{-6}\approx 1$: the first step mainly removes the high-eigenvalue error. This is the spectral weighting of conjugate gradient error. The expectation is conditional on the coefficients, not guaranteed by the two eigenvalues alone; if $\beta =0$, the first step instead solves the low-eigenvalue component exactly.