Past exam of the numerical linear algebra course of the University of Oxford / 2025 / 3

Write $x=Q_{k}y$ and $\beta =\Vert b{\Vert }_2$. The Arnoldi decomposition and $b=\beta Q_{k+1}{e}_1$ giveBecause $Q_{k+1}$ has orthonormal columns, it preserves the Euclidean norm. The generalized minimal residual method therefore reduces to the small linear least-squares problemUse Givens rotations on the upper Hessenberg matrix to obtain a QR decompositionwhere $Z$ is an orthogonal matrix and $R_k$ is an upper triangular matrix. Then solve $R_k{y}_k=g_{1:k}$ by backward substitution in a triangular system and form $Q_k{y}_k$. The resulting minimum linear-system residual has Euclidean norm $|g_{k+1}|$.

The small problem has a unique answer: if ${\tilde{H}}_{k}y=0$, the Arnoldi decomposition gives $A{Q}_{k}y=0$. The invertible matrix $A$ implies $Q_{k}y=0$, and the orthonormality of $Q_k$ then gives $y=0$. Thus ${\tilde{H}}_k$ has full column rank and $R_k$ is nonsingular.

For the small QR decomposition, each of the $k$ Givens rotations eliminates one subdiagonal entry of the upper Hessenberg matrix. At column $j$, apply the rotation to two rows in columns $j,\dots ,k$, costing $O(k-j+1)$ arithmetic operations. Summing gives $O(k^2)$, and applying the rotations to the right-hand side costs $O(k)$. Backward substitution in a triangular system costs $O(k^2)$.

Computing $\Vert b{\Vert }_2$ costs $O(n)$ and the final dense matrix multiplication $Q_k{y}_k$ costs $O(nk)$. Thus the total, with the Arnoldi decomposition already supplied, isHere $1\le k<n$, since $Q_{k+1}$ has orthonormal columns in ${\mathbb{R}}^n$. In the standard dense arithmetic model, forming the general matrix–vector product gives the sharp $\Theta (nk)$ cost. The pair above records that two-parameter bound; looser bounds such as $O(n^2)$ are also possible because $k<n$. Treating ${\tilde{H}}_k$ as an unstructured dense matrix would unnecessarily replace its $O(k^2)$ factorization by an $O(k^3)$ computation.

The claim is true, withIt is a symmetric positive-definite matrix: for every nonzero vector $z$, the fact that $A$ is an invertible matrix gives $z^{T}Mz=\Vert Az{\Vert }_2^2>0$. Since $A{x}_{\ast }=b$, every candidate $x$ satisfiesThe generalized minimal residual method therefore minimizes precisely this weighted norm of the error over the linear span of $Q_k$. The same matrix $M$ works for every $k$ and does not depend on $b$.

Multiplying the Arnoldi decomposition by $Q_k^T$ shows that the leading square block isIt is an upper Hessenberg matrix by the Arnoldi iteration. Because $A=A^T$, it is also a symmetric matrix. Entries below the first subdiagonal vanish by the Hessenberg structure, and their reflected entries above the first superdiagonal vanish by symmetry. ConsequentlyThe extra last row has only its final entry possibly nonzero. With the usual choice of positive normalization factors, the individual column relations are the three-term Lanczos algorithm recurrencewhere the first term is omitted for $j=1$. This reduction is an exact-arithmetic statement; round-off error can affect the orthogonality in an implementation.

For any nonzero vector $y\in {\mathbb{R}}^k$, orthonormality gives $\Vert Q_{k}y{\Vert }_2=\Vert y{\Vert }_2>0$. Since $A$ is a symmetric positive-definite matrix,The preceding part already gives symmetry, so $T_k$ is a symmetric positive-definite matrix. In fact the spectral theorem gives the useful boundswhich place its eigenvalues between the smallest and largest eigenvalues of $A$.

The least-squares normal equations areThe full column rank of $G$ gives $z^{T}Hz=\Vert Gz{\Vert }_2^2>0$ for nonzero $z$, so $H$ is a symmetric positive-definite matrix. Apply the Conjugate gradient method to this linear system. To verify that its solution is the desired one, let $x_{\ast }$ solve the displayed equations. Then $G^T(G{x}_{\ast }-b)=0$, and for any $x$ the Pythagorean identity givesThus $x_{\ast }$ is the unique minimizer of the linear least-squares problem.

One can apply $H$ without forming the Gram matrix: compute $Hp=G^T(Gp)$. More explicitly, initialize $x_0$, $r_0=c-H{x}_0$, $p_0=r_0$; until the linear-system residual vanishes, use the Conjugate gradient method recurrenceEach dense application of $H$ this way costs $O(mn)$ operations. The condition number satisfies ${\kappa }_2(H)={\kappa }_2(G{)}^2$, where ${\kappa }_2(G)={\sigma }_{\max }(G)/{\sigma }_{\min }(G)$; this is the conditioning cost of the least-squares normal equations.

Under the permitted singular value bounds, $m=4n$ givesFor the symmetric positive-definite matrix $H=G^{T}G$ in part (iii), its eigenvalues are the squared singular values of $G$, soApply the supplied conjugate gradient convergence bound to $H$. For $e_k=x_k-x_{\ast }$ and nonzero $e_0$, the function $(\sqrt{\kappa }-1)/(\sqrt{\kappa }+1)$ increases with $\kappa \ge 1$, givingThus the error in this weighted norm contracts at a guaranteed exponential factor of at most $1/2$ per iteration, up to the prefactor $2$. Here $\Vert e{\Vert }_H=\Vert Ge{\Vert }_2$. The Pythagorean identity in part (iii) also identifies $\Vert e_k{\Vert }_H^2$ with the excess squared linear least-squares problem objective, which consequently decreases by a bound proportional to $4^{-k}$. For the Euclidean norm, the same eigenvalue bounds give $\Vert e_k{\Vert }_2/\Vert e_0{\Vert }_2\le 62^{-k}$.

For relative accuracy $0<ϵ<1$ in the $H$-weighted norm, it is sufficient to takeCompute $c=G^{T}b$ once in $O(mn)$ operations. Each iteration applies $p\mapsto G^T(Gp)$ with two dense matrix multiplications and uses $O(n)$ further arithmetic operations; when $m=4n$ this costs $O(n^2)$. Hence the total cost isUsing finite termination of the conjugate gradient method in exact arithmetic refines this to $O(n^2\min \{n,\lceil {\log }_2(2/ϵ)\rceil \})$. Explicitly constructing $G^{T}G$ first would add $O(m{n}^2)=O(n^3)$ setup work and is unnecessary. These accuracy estimates concern relative error; an absolute tolerance replaces $2/ϵ$ by $2\Vert e_0{\Vert }_H/ϵ$ when $\Vert e_0{\Vert }_H>ϵ$. If $e_0=0$, the initial guess already solves the problem.

As in part 2(b)(v), the Gaussian assumption alone does not make the printed singular value interval a certainty. For example, the event that every entry of $G$ has absolute value below a sufficiently small positive number has positive probability and forces $\Vert G{\Vert }_F<\sqrt{n}$, contradicting the asserted lower endpoint. The displayed rate is therefore conditional on the singular-value bounds authorized by the paper. The cost and convergence derivation are exact under that condition; an unconditional probabilistic statement would need an additional failure-probability estimate.