Numerical linear algebra course of the University of Oxford (C6.1)

Solutions to all three questions of the University of Oxford numerical linear algebra paper, sat on 2 June 2025. Source: released examination PDF, parent PDF pages 1–3. The question statements are omitted. The Gaussian-matrix estimates in parts 2(b)(v) and 3(b)(iv) are used as the paper explicitly permits; their limitations are explained in those solutions.

The orthonormal columns give $A^{T}A=I_n$. Thus $P=A{A}^T$ is the orthogonal projection onto the column space of $A$, andThe two terms are orthogonal, so the Pythagorean identity givesOnly the first term depends on $x$. Consequently the linear least-squares problem has the unique answer

For every vector $z\in {\mathbb{R}}^n$, splitting the rows givesBoth summands are nonnegative. Taking the supremum over unit vectors in the definition of the spectral norm proves

Put $B=A_1^T{A}_1$. It is a symmetric matrix and a positive semidefinite matrix, while $A_2^T{A}_2=I_n-B$ is also a positive semidefinite matrix. By the spectral theorem, choose an orthogonal matrix $V=[v_1,\dots ,v_n]$ such thatThe bounds follow directly from $z^{T}Bz=\Vert A_{1}z{\Vert }_2^2\in [0,\Vert z{\Vert }_2^2]$. Define the diagonal matricesThe columns of $A_{1}V$ are mutually orthogonal, since $(A_1{v}_i{)}^T(A_1{v}_j)={\lambda }_j{\delta }_{ij}$. For each ${\lambda }_j>0$, set $u_{1j}=A_1{v}_j/\sqrt{{\lambda }_j}$. These columns are orthonormal. At a zero ${\lambda }_j$, one has $A_1{v}_j=0$; choose $u_{1j}$ by completing the already prescribed columns to an orthonormal set of $n$ columns in ${\mathbb{R}}^k$. Such a completion exists because $k\ge n$. This gives $A_{1}V=U_1{\Sigma }_1$, including the zero columns.

Likewise, the Gram matrix of $A_{2}V$ is $\mathrm{diag}(1-{\lambda }_j)$. Normalize its nonzero columns and complete at its zero columns in ${\mathbb{R}}^{m-k}$, using $m-k\ge n$. The resulting $U_2$ has orthonormal columns and $A_{2}V=U_2{\Sigma }_2$. Combining the identities givesThis is a form of the cosine-sine decomposition. Completing the zero columns is essential: neither block was assumed to have full column rank.

Using the orthonormal columns of $U_1,U_2$ in the factorization givesMultiplying by the orthogonal matrices $V^T$ and $V$ gives

Choose the nonnegative square roots of the eigenvalues of $A^{T}A$ for the diagonal entries of $\Sigma$, in nonincreasing order. Its assumed matrix inverse exists, so these entries are positive. The eigendecomposition givesHence $U$ has orthonormal columns. AlsoTogether with the orthogonal matrix $V$ and the nonnegative diagonal matrix $\Sigma$, these are precisely the properties of a thin singular value decomposition.

Let $\tilde{A}=\hat{U}\hat{\Sigma }{\hat{V}}^T$. Interpret the diagonal of $\hat{\Sigma }$ as nonnegative computed singular values and sort them consistently. If signed diagonal entries are used instead, take their absolute values and absorb signs into the orthonormal columns. Because the two outer factors have the stated orthonormality, these sorted numbers are exactly ${\sigma }_j(\tilde{A})$.

For completeness, the singular value perturbation bound follows from the Eckart-Young-Mirsky theorem. If $B_{j-1}$ is the truncated singular value decomposition of a matrix $B$, the triangle inequality givesExchange $A$ and $B$ to get $|{\sigma }_j(A)-{\sigma }_j(B)|\le \Vert A-B{\Vert }_2$. Applying this to $B=\tilde{A}$ and the given backward error provesThus backward stability of the singular value decomposition gives an absolute forward error bound for every singular value, even a very small one.

An error of order $ϵ$ in a squared singular value can become an error of order $\sqrt{ϵ}$ after taking its square root. To keep $A$ nonsingular as required, take $0<ϵ<1/2$ andThen $\Vert A{\Vert }_2=1$, and the computed eigendecomposition has exactly the permitted backward error:Nevertheless the second computed singular value has errorBy the singular value perturbation bound, every product with these computed singular values and exactly orthonormal outer factors has spectral norm error at least this large. In particular, even the favorable choice $\hat{U}=I_2$ giveswhich has no uniform $O(ϵ)$ bound in big O notation. Therefore backward stability for the eigendecomposition of the Gram matrix does not imply backward stability for the singular value decomposition.

There is a related loss of orthonormality in the raw formula: evaluating $A\hat{V}{\hat{\Sigma }}^{-1}$ exactly for these perturbed factors gives $\mathrm{diag}(1,ϵ/\sqrt{{ϵ}^2+ϵ})$, whose second column has length tending to zero. Its product with $\hat{\Sigma }{\hat{V}}^T$ equals $A$, but its diagonal factors are consequently not a singular value decomposition. The paper's idealized orthonormality assumption lets us isolate the unavoidable singular value error even if this additional problem is repaired.

No uniform absolute error bound of order $ϵ$ follows. Take the Jordan blockHere $\Vert B{\Vert }_2=1$ and $\Vert E{\Vert }_2=ϵ$. The exact eigenvalues of $B$ are both zero, whereasThese are admissible outputs under a backward error guarantee of size $ϵ$, yet their forward errors are $\sqrt{ϵ}$. This is the square-root splitting of a defective double eigenvalue. More generally, eigenvalue sensitivity depends on the matrix; backward stability alone does not control it. The singular value perturbation bound used in part (ii) has no equally uniform counterpart for the eigenvalues of arbitrary nonsymmetric matrices.

By the Eckart-Young-Mirsky theorem, truncate the singular value decomposition:The same $A_r$ minimizes the spectral norm and Frobenius norm errors over matrices of matrix rank at most $r$. This conventional meaning of rank-$r$ low-rank approximation includes the case $\mathrm{rank}A<r$, when $A_r=A$.

Choose the first $r$ left singular vectors:Their orthonormality gives $Q_{\ast }^T{Q}_{\ast }=I_r$. The orthogonal projection $Q_{\ast }{Q}_{\ast }^T$ fixes $u_j$ for $j\le r$ and annihilates $u_j$ for $j>r$. Applying it to the singular value decomposition therefore gives

Put $P=Q{Q}^T$. Since $Q$ has $r$ columns, the rank bound for a matrix product gives $\mathrm{rank}(PA)\le r$. The Eckart-Young-Mirsky theorem and part (a) then giveOne can also see the lower bound directly. The restriction of $PA$ to the $(r+1)$-dimensional linear span of $v_1,\dots ,v_{r+1}$ has a nonzero null space, because its matrix rank is at most $r$. For a unit vector $z$ in that null space, $\Vert (A-PA)z{\Vert }_2=\Vert Az{\Vert }_2\ge {\sigma }_{r+1}$ by the singular value decomposition.

The orthogonal projection $P=Q{Q}^T$ givesThe first summand has every column in the column space of $Q$, and the second has every column in its orthogonal complement. They therefore have zero Frobenius inner product. The Pythagorean identity, together with $Q^{T}Q=I_r$, gives the stronger exact formulaConsequentlywith equality exactly when $M=Q^{T}A$. Thus solving the linear least-squares problem for each column chooses the best coefficients in the prescribed column space.

The thin QR decomposition gives $AX=QR$. Its columns already belong to the column space of $Q$, so the orthogonal projection $Q{Q}^T$ leaves them fixed:This identity uses only $AX=QR$ and $Q^{T}Q=I_r$; it remains valid if $R$ is singular, provided such a factorization with orthonormal $Q$ is used.

Because $r<n$, choose a unit vector $v$ in the null space of $X^T$, so that $v^{T}X=0$. Let $J\in {\mathbb{R}}^{m\times n}$ have orthonormal columns, for example the first $n$ columns of $I_m$, and choose a real number $\alpha >1$. SetIts singular value along $v$ is $\alpha$, and its other singular values are $1$. Thus ${\sigma }_1(A)=\alpha >1={\sigma }_2(A)$. Meanwhile $AX=JX$. If $X$ has full column rank, the column space of the thin QR decomposition is exactly $\mathrm{range}(JX)$, which is orthogonal to $Jv$. Hence $PJv=0$ for $P=Q{Q}^T$, andConversely, an orthogonal projection is a contraction in the Euclidean norm, so $\Vert (I_m-P)A{\Vert }_2\le \Vert A{\Vert }_2=\alpha$. ThereforeIf $X$ is rank deficient, complete the QR decomposition within $J(v^{\perp })$: this space has dimension $n-1\ge r$ and contains all columns of $JX$, so the required orthonormal completion exists there. The same proof applies. For rank-deficient data the choice of unused columns of $Q$ matters; the construction specifies a valid choice. If $Q$ instead contains only a basis of $\mathrm{range}(AX)$, the conclusion holds without any completion. A fixed sampling matrix can thus miss the unique largest singular value completely.

Write $\ell =r/2$; the requested expression requires the positive even integer $r\ge 2$. Split the thin singular value decomposition after $\ell$ columns:Here $\Vert D_b{\Vert }_2={\sigma }_{\ell +1}=\Vert A-A_{\ell }{\Vert }_2$. Each column of a Gaussian random matrix is an orthogonally invariant random vector: its probability density function is proportional to $\exp (-\Vert z{\Vert }_2^2/2)$, and an orthogonal matrix preserves both $\Vert z{\Vert }_2$ and volume. Thus $V^{T}X$, where $V=[V_a,V_b]$, again has independent and identically distributed standard normal entries. In particular both displayed blocks have the Gaussian distributions needed below.

First prove a deterministic error bound for a randomized range finder. Put $P=Q{Q}^T$ and suppose ${\Omega }_a$ has full row rank. Its Moore-Penrose pseudoinverse is ${\Omega }_a^{†}={\Omega }_a^T({\Omega }_a{\Omega }_a^T{)}^{-1}$, so ${\Omega }_a{\Omega }_a^{†}=I_{\ell }$. The QR decomposition implies $(I_m-P)AX=0$, henceMultiplication by ${\Omega }_a^{†}$ expresses the missed leading columns in terms of the trailing columns. Substitution into $(I_m-P)A$ gives the exact identityFor any matrix $F$, the spectral norm of $[-FI]$ is $\sqrt{1+\Vert F{\Vert }_2^2}$, because its product with its transpose is $F{F}^T+I$. Using submultiplicativity of the operator norm, the orthonormality of $U_b$, and the contraction property of the orthogonal projection, we obtainThis estimate is deterministic and also covers rank-deficient $A$; if ${\sigma }_{\ell +1}=0$, its right-hand side vanishes.

Now use the paper's permitted singular value estimates. Apply the lower estimate to the tall matrix ${\Omega }_a^T$ of size $r\times \ell$, rather than to its wide transpose. It givesThe upper estimate for ${\Omega }_b$, transposing if necessary, givesThus, with $K=2/(1-1/\sqrt{2})$,Since $n/r>1$, the deterministic estimate yieldsThis proves the requested conclusion under the explicitly authorized estimates.

The probabilistic qualification is necessary. A Gaussian random matrix does not always have all its singular values between the two endpoints printed in the hint. Already a $1\times 1$ example would require $|Z|\le 2$ for every standard normal random variable $Z$, which fails with positive probability. The proof above holds on the event where the two block estimates hold. In a probabilistic formulation one supplies estimates for that event, usually allowing a positive margin in both endpoints; no failure probability is specified by the paper.

Nor can the requested inequality with one universal $C$ hold almost surely for every $A$. For example, take $n=m=3$, $r=2$, $A=\mathrm{diag}(L,1,1)$ with $L>1$, and a sampling matrixAt $X_0$, the sampled column space misses the first coordinate, so $\Vert A-\hat{A}{\Vert }_2=L$, whereas $\Vert A-A_1{\Vert }_2=1$. Choose $L>C\sqrt{3/2}$. Because $A{X}_0$ has full column rank, its orthogonal projection depends continuously on $X$ near $X_0$ through $AX(X^T{A}^{T}AX{)}^{-1}{X}^T{A}^T$. The strict violation persists in an open neighborhood of $X_0$. The joint Gaussian probability density function is positive there, so that neighborhood has positive probability. This identifies the missing qualification without weakening the conditional calculation.

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.