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

An orthogonal projector sends $x$ to the foot $Px$ of its perpendicular onto the image of a linear map $\mathrm{range}P$. The residual $r=x-Px$ is orthogonal to the entire range, not just to the single vector $Px$. Indeed, for any range vector $Pv$,In particular $\boxed{(Px{)}^T(x-Px)=0}$. The orthogonal projection diagram illustrates the resulting right triangle.

The vector Euclidean norm is $\boxed{\Vert x{\Vert }_2=(x^{T}x{)}^{1/2}=(\sum _i{x}_i^2{)}^{1/2}}$, the ordinary length from the origin. The induced matrix norm isGeometrically this is the largest stretch of a unit vector. A linear map sends the unit sphere to an ellipsoid, possibly degenerate; its largest semiaxis has length $\Vert A{\Vert }_2$. Equivalently the matrix 2-norm is the largest singular value.

The perpendicular decomposition in part (a) gives, by the Pythagorean theorem,Thus an orthogonal projection cannot lengthen a vector and $\boxed{\Vert P{\Vert }_2\le 1}$. Algebraically the cross term vanishes because $P^T(I-P)=P-P^2=0$, giving the same identity directly from the hypotheses. If $P\ne 0$, take a nonzero $x$ in its range; then $Px=x$ and $\Vert P{\Vert }_2=1$. For $P=0$ the norm is zero.

The nullspace is the vector subspace $\boxed{\ker A=\{x\in {\mathbb{R}}^n:Ax=0\}}$. A full real singular value decomposition iswhere $U,V$ are orthogonal matrices and ${\Sigma }_{ij}=0$ for $i\ne j$, with ${\Sigma }_{jj}={\sigma }_j\ge 0$ for $1\le j\le n$. We may order ${\sigma }_1\ge \cdots \ge {\sigma }_n\ge 0$. The columns of $V$ form an orthonormal basis of the input space ${\mathbb{R}}^n$, and the columns of $U$ form one of the output space ${\mathbb{R}}^m$. The positive singular values number exactly the matrix rank.

Let $r$ be the matrix rank, and partition the right singular vectors as $V=(V_r{V}_0)$, with $V_0\in {\mathbb{R}}^{n\times (n-r)}$ consisting of the columns associated with zero singular values. From $Ax=0⟺\Sigma V^{T}x=0$, these columns form an orthonormal basis of the nullspace. Therefore its orthogonal projector isIt acts on the input space ${\mathbb{R}}^n$. Symmetry and idempotence follow from $V_0^T{V}_0=I_{n-r}$. At full column rank it is zero, and at rank zero it is the identity.

At elimination column $k$, let $r=m-k$ be the size of the remaining trailing block. Computing its multipliers costs $r$ divisions; updating all $r^2$ entries costs $2r^2$ floating-point operations, one multiplication and one subtraction per entry. Thus the dominant count isUpdating the right side and back substitution add only $O(m^2)$ floating-point operations. Pivot searches involve comparisons and row exchanges, without changing the leading arithmetic count. Hence Gaussian elimination with pivoting has the requested asymptotic cost.

An orthogonal similarity transformation is $H=Q^{T}AQ$ with $Q^{T}Q=I$. An upper Hessenberg matrix has $H_{ij}=0$ whenever $i>j+1$; its first subdiagonal may be nonzero.

At step $k=1,\dots ,m-2$, choose a Householder reflector $Q_k=\mathrm{diag}(I_k,{\hat{Q}}_k)$ whose trailing block sends the current column segment in rows $k+1,\dots ,m$ to a multiple of its first coordinate. The left multiplication zeros positions $(k+2,k),\dots ,(m,k)$. Apply the matching right multiplication to retain similarity: $A\leftarrow Q_k^{T}A{Q}_k$. Since $Q_k$ acts only on rows and columns $k+1,\dots ,m$, earlier Hessenberg zeros are preserved. Finally $Q=Q_1\cdots Q_{m-2}$ gives $H=Q^{T}AQ$.

Each reflector is applied through a rank-one update, costing $O(m^2)$ rather than a full matrix multiplication. Across $m-2$ steps, the Hessenberg reduction costs $O(m^3)$.

Solving the $m$ unrelated shifted systems by separate Gaussian elimination costs $\boxed{\frac{2}{3}{m}^4+O(m^3)}$ floating-point operations.

Instead compute $H=Q^{T}AQ$ once in $O(m^3)$ work. For each shift, solveThe shift changes only the diagonal, so $H-{\lambda }_{j}I$ remains an upper Hessenberg matrix. In each elimination column only the entry immediately below the pivot needs removal. With partial pivoting the only possible nonzero candidate rows are the current row and the next, and each update touches $O(m)$ entries. Factorization therefore costs $O(m^2)$ per shift; the triangular solve and both orthogonal-vector transformations also cost $O(m^2)$ each.

Thus $m$ shifted solves plus the initial Hessenberg reduction cost $\boxed{O(m^3)}$ in total. Pivoting is still needed in general: nonsingularity excludes a singular system but does not by itself guarantee that every unpivoted elimination pivot is nonzero.

The off-diagonal Toeplitz matrix $T=A-I$ is symmetric and has row and column absolute sums bounded byFor an elementary Euclidean bound, apply the weighted Cauchy-Schwarz inequality in each row:Sum over $i$ and use both absolute-sum bounds to get $\Vert Tx{\Vert }_2^2\le r^2\Vert x{\Vert }_2^2$. ThereforeThe entry $t_1$ is unused because the diagonal of $A$ is fixed to one; including it only weakens the bound. This is the row-column bound for the matrix 2-norm.

A real symmetric matrix is positive-definite when $x^{T}Ax>0$ for every nonzero real vector $x$. The preceding matrix 2-norm bound givesThus $A$ is a symmetric positive-definite matrix, independently of the detailed positions of the nonzero entries of $t$.

An eigenvalue of $A$ is a scalar $\lambda$ for which $Az=\lambda z$ has a nonzero eigenvector. Allowing initially $z\in {\mathbb{C}}^m$, symmetry of the real matrix gives $A^{\ast }=A$, so $z^{\ast }Az$ is real. Since $z^{\ast }z>0$,This proves reality without presuming it when choosing the eigenvector. Once $\lambda$ is real, the real or imaginary part of $z$ supplies a nonzero real eigenvector.

For a unit eigenvector of $A$, $(A-I)x=(\lambda -1)x$. The matrix 2-norm bound from part (a) gives $|\lambda -1|\le 2/3$. Hence every eigenvalue lies inIn particular the condition number in the Euclidean operator norm satisfies ${\kappa }_2(A)={\lambda }_{\max }/{\lambda }_{\min }\le 5$.

The energy norm associated with this symmetric positive-definite matrix isIt is the norm induced by the inner product $⟨u,v{⟩}_A=u^{T}Av$. The conjugate gradient method minimizes this error norm over its successive affine Krylov subspaces.

The standard conjugate gradient convergence bound isHere $\kappa \le 5$, so the factor in parentheses is at most $q=(3-\sqrt{5})/2\simeq 0.381966$. Since $2q^{25}<10^{-10}$, 25 iterations suffice to reduce the energy-norm error by a factor exceeding $10^{10}$.

There are at most $m/20$ nonzero entries of $t$, and each contributes to at most two off-diagonals. Thus $A-I$ has at most $m^2/10$ nonzero entries. Initialize a matrix-vector product with the diagonal contribution $p$, which costs no arithmetic because the diagonal is one. Each off-diagonal term needs one multiplication and one addition, giving at most $m^2/5$ floating-point operations per product.

A standard conjugate gradient step adds two inner products, three vector updates and two scalar divisions, at most $10m$ further operations. Starting from $x_0=0$ makes $r_0=b$ and avoids an initial matrix-vector product; its squared norm needs $2m-1$ operations. The total for 25 steps is at mostThis supplies an explicit meaning of the paper's large-$m$ qualification. The quoted conjugate gradient convergence bound is the usual exact-arithmetic result.