Past exam of the numerical linear algebra and continuous optimisation course of the University of Oxford / 2007 / A

The outer product has entries $\boxed{(x{y}^T{)}_{ij}=x_i{y}_j}$. Every column is a multiple of $x$, and at least one column is nonzero because $y\ne 0$. Its column space is therefore exactly $\mathrm{span}\{x\}$, giving $\boxed{\mathrm{rank}(x{y}^T)=1}$.

A Householder matrix is $H=I-2v{v}^T/(v^{T}v)$, $v\ne 0$. With $P=v{v}^T/(v^{T}v)$ one has $P^T=P$ and $P^2=P$, so $H^{T}H=(I-2P{)}^2=I$. Thus it is an orthogonal matrix and is its own inverse.

To express an arbitrary orthogonal $Q$ as a product of these reflections, argue by induction on its dimension. Its first column $q$ is a unit vector. If $q\ne e_1$, take $v=q-e_1$. Then $v^{T}q=1-q_1$ and $v^{T}v=2(1-q_1)$, so the corresponding reflection satisfies $Hq=q-v=e_1$. If $q=e_1$, no reflection is required. The resulting orthogonal matrix has first column $e_1$; orthogonality of its columns then makes its first row $e_1^T$ as well. It is consequently $\mathrm{diag}(1,Q^{\prime })$ with $Q^{\prime }$ orthogonal. Apply the induction hypothesis to $Q^{\prime }$ and embed its reflections in the full space. In dimension one the only possibilities are $1$, the empty product, and $-1$, itself a reflection. Reversing the performed transformations proves the Householder factorization of an orthogonal matrix, with at most the dimension many reflections.

A full real singular value decomposition iswhere $\Sigma \in {\mathbb{R}}^{m\times n}$ has only diagonal entries, with ${\sigma }_1\ge \cdots \ge {\sigma }_{\min (m,n)}\ge 0$. These are the singular values. The number of positive singular values is the matrix rank. The columns of $V$ are right singular vectors and those of $U$ are left singular vectors, giving orthonormal input and output coordinate systems.

For a rectangular diagonal matrix $D$, define $D^{+}\in {\mathbb{R}}^{n\times m}$ by transposing its shape and takingEntry by entry, $D{D}^{+}D=D$ and $D^{+}D{D}^{+}=D^{+}$. For $A=U\Sigma V^T$, orthogonality therefore givessince $A{A}^{+}A=U\Sigma {\Sigma }^{+}\Sigma V^T=A$, and similarly $A^{+}A{A}^{+}=A^{+}$. This candidate is the Moore-Penrose inverse. It also makes $A{A}^{+}$ and $A^{+}A$ symmetric, although those two extra conditions are not required by the question's stated definition. The two displayed conditions alone need not specify a unique generalized inverse.

Direct solution usually uses a pivoted LU decomposition $PA=LU$: solve $Lz=Pb$ by forward substitution, then $Ux=z$ by back substitution. For a symmetric positive-definite matrix, the Cholesky decomposition iswith $C$ lower triangular and positive diagonal. It exploits symmetry, requiring about $n^3/3$ arithmetic operations rather than the leading $2n^3/3$ for general elimination.

If $A$ has full column rank and its reduced QR decomposition is $A=QR$, $Q^{T}Q=I$, thenChoose the signs of $Q,R$ so that $R$ has positive diagonal. Then $R$ is precisely the upper Cholesky factor of the normal-equation matrix $A^{T}A$, or $R^T$ is its lower factor. This identifies the factors; in computation, obtaining $R$ directly by QR decomposition avoids explicitly forming the normal equations, which square the condition number.

The matrix splitting produces $M{x}^{(k+1)}=N{x}^{(k)}+b$. If the iterates converge to $z$, passing to the limit yields $Mz=Nz+b$, hence $Az=b$. Thus any finite limit solves the system; when $A$ is nonsingular it is the unique solution.

The iteration matrix is $G=M^{-1}N=I-M^{-1}A$. For $M$ symmetric positive definite, its positive-definite square root giveswhich is symmetric when $A$ is symmetric. Let ${\theta }_i$ be the eigenvalues of $\tilde{A}=M^{-1/2}A{M}^{-1/2}$. In its orthonormal eigenbasis the transformed error is multiplied by $1-{\theta }_i$ at each step. Therefore convergence to the unique solution for every starting vector is equivalent toSufficiency follows because each scalar mode tends to zero. Necessity follows by starting with a nonzero error in any mode whose multiplier has modulus at least one. Equivalently the condition says $\boxed{A\text{ and }2M-A\text{ are positive definite}}$, by congruence of their quadratic forms.

Under the usual nonsingularity assumption, a symmetric $A$ that is not positive definite has a negative quadratic direction. Its congruent $\tilde{A}$ then has a negative eigenvalue, making an error multiplier $1-{\theta }_i>1$. A starting error in that mode diverges, proving the final claim in that setting.

If singular matrices are permitted, the last claim needs qualification. A zero eigenvalue gives a persistent error mode, but the iterates may still converge to another solution of a compatible singular system. For example $A=0$, $M=I$, $b=0$ makes every iterate constant and already a solution, despite $A$ not being positive definite. Thus literal nonconvergence cannot be deduced from a zero eigenvalue alone.

A multigrid method alternates local error damping with correction on coarser grids. For clarity write $A_h=h^{-2}$ times the five-point stencil and $A_h{u}_h=f_h$, incorporating Dirichlet boundary data into the right side. One two-grid cycle proceeds as follows.

The coarse solve is performed recursively, with a direct solve on the smallest grid, giving a V-cycle. Smoothing treats oscillations that are cheap for local relaxation to remove; coarse-grid correction treats errors that appear smooth on the fine grid. The scaling of the residual and coarse operator must agree if one instead uses the paper's unscaled stencil.

For homogeneous error boundary values, the tensor-product sine vectors $U_{j,k}=\sin (r\pi jh)\sin (s\pi kh)$ form a basis. The sine addition formula gives $U_{j+1,k}+U_{j-1,k}=2\cos (r\pi h)U_{j,k}$, and similarly in the other coordinate. Hence the five-point matrix hasThere is a printing error in the displayed iteration. Taken literally, its error multiplier is $1/20-(2/5)(\cos (r\pi h)+\cos (s\pi h))$, and its fixed-point equation is $(19/4)U_{j,k}+{\sum }_{\mathrm{neighbors}}U=h^2{f}_{j,k}$, rather than the stated Poisson equation. Thus the printed affine update cannot be used as a consistent relaxation for this system.

The intended weighted Jacobi iteration with $\omega =4/5$ isIts iteration matrix is $I-(\omega /4)A$, whose mode multiplier isIf at least one of the two cosines is negative, their sum lies strictly between $-2$ and $1$, so $-3/5<{\mu }_{r,s}<3/5$. Thus every such high-frequency mode contracts by at most $\boxed{3/5}$ per sweep, uniformly in grid spacing. This weighted Jacobi smoothing bound for the square Poisson problem justifies the corrected method as a smoother. The literal printed update also damps modes, but toward the wrong fixed-point equation.

The Krylov subspace is $K_k(A,r)=\mathrm{span}\{r,Ar,\dots ,A^{k-1}r\}$. Write $x_k-x_0=q(A)r_0$ with $\mathrm{deg}q\le k-1$. For the exact solution, $r_0=A(x-x_0)$, soThis polynomial has degree at most $k$ and $p(0)=1$. It is the error-polynomial representation for a Krylov subspace method.

Every point in the affine search space is uniquely $x_k=x_0+V_{k}y$, since the columns of $V_k$ form a basis. Its linear-system residual is $r_k=r_0-A{V}_{k}y$, so minimizing the residual is exactly the linear least-squares problem ${\min }_y\Vert r_0-A{V}_{k}y{\Vert }_2$.

Use the standard starting choice $v_1=r_0/\Vert r_0{\Vert }_2$ and an orthonormal extended basis $V_{k+1}$. If $A{V}_k=V_{k+1}{L}_k$, then $r_0=V_{k+1}(\Vert r_0{\Vert }_2{e}_1)$ and orthonormality givesThe starting-vector choice is needed for this particular $e_1$ representation; an arbitrary orthonormal basis need not have those initial coordinates.

The Arnoldi iteration constructs the basis. At column $j$, orthogonalize $A{v}_j$ against $v_1,\dots ,v_j$, recording coefficients $h_{ij}=v_i^{T}A{v}_j$; normalize the remaining vector with $h_{j+1,j}$ to get $v_{j+1}$. This makes $L_k$ an upper Hessenberg matrix of size $(k+1)\times k$. The resulting GMRES algorithm starts from $x_0,r_0$, grows this Arnoldi basis, solves the small least-squares problem, and sets $x_k=x_0+V_{k}y$. Its small residual norm provides the stopping test. Incremental Givens rotations efficiently update the small QR decomposition instead of solving it anew. Zero initial residual means the system is already solved; Arnoldi breakdown gives an invariant subspace and, for nonsingular $A$, permits exact solution within it.

One restarted step chooses $x_j=x_{j-1}+{\alpha }_j{r}_{j-1}$ to minimize the new residual norm. Since $r_{j-1}=A(x-x_{j-1})$,and induction givesThese polynomials have degree at most one; a zero step gives the constant polynomial one.

In the specified example the first search space is spanned by $b=(1,1{)}^T$. Writing $x_1=\alpha b$, its squared residual norm is $(1-\alpha {)}^2+(1+\alpha {)}^2=2+2{\alpha }^2$, uniquely minimized at $\alpha =0$. Thus restarted GMRES gives $\boxed{x_1=x_2=(0,0{)}^T}$ and stagnates forever, repeatedly facing the same minimization problem.

Without restarting, the first iterate is still zero, but $K_2(A,b)=\mathrm{span}\{(1,1{)}^T,(1,-1{)}^T\}={\mathbb{R}}^2$. The second minimization finds the exact solution $\boxed{x_2=(1,-1{)}^T}$. This one-step restarted GMRES stagnation example shows how discarding the accumulated search space can prevent convergence even for a nonsingular two-dimensional system.