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

Original solutions to both sections, all six questions and all twenty-nine terminal parts of the Numerical Linear Algebra and Continuous Optimisation paper sat on 8 June 2010. Source: released Mathematics papers, Part C 2007 11.pdf. Statements are not reproduced.

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.

At an unconstrained local minimizer of a twice continuously differentiable objective function, the first-order optimality condition is $\boxed{\nabla f(x_{\ast })=0}$. The second-order necessary condition is $d^T{\nabla }^{2}f(x_{\ast })d\ge 0$ for every direction $d$, so the Hessian matrix is positive semidefinite.

For the Euclidean norm, the unit steepest descent direction is $-g/\Vert g{\Vert }_2$ when $g=\nabla f(x)\ne 0$; the algorithm usually uses the unnormalized direction $p=-g$. It is a descent direction because $\boxed{g^{T}p=-\Vert g{\Vert }_2^2<0}$.

Two limitations of unpreconditioned steepest descent method in high-accuracy smooth optimization are its slow zigzagging on ill-conditioned valleys and its sensitivity to scaling of the coordinates. It uses no curvature information to correct those directions, so it can require many function and gradient evaluations compared with a suitable Newton or quasi-Newton method. These are reasons to prefer curvature-aware methods in that setting, rather than a claim that gradient-based methods are never useful.

The descent lemma follows by integrating the gradient along the line and using its Lipschitz continuity:This upper bound is at most the right side of the Armijo condition when $\alpha [(1-\beta )g^{T}p+\gamma \alpha \Vert p{\Vert }_2^2/2]\le 0$. Since $g^{T}p<0$, it holds forThe zero step gives equality; every positive step in this interval satisfies the prescribed sufficient-decrease test.

The quadratic model has gradient $g(x_k)+B_{k}p$ and positive-definite Hessian matrix $B_k$. Its unique minimizer is thereforeIf $g(x_k)\ne 0$, then $p_k\ne 0$ and $g(x_k{)}^T{p}_k=-p_k^T{B}_k{p}_k<0$, proving that it is a descent direction.

Backtracking line search starts with a trial step and repeatedly reduces it until the Armijo condition holds. Part (ii) guarantees finite termination because sufficiently small steps are acceptable. The test guarantees a specified positive fraction of the decrease predicted by the linear model, namely $f(x_k+\alpha p_k)\le f(x_k)+\alpha \beta g(x_k{)}^T{p}_k<f(x_k)$ for an accepted positive step. It does not require solving the line-minimization problem exactly.

If the first trial is accepted, ${\alpha }_k={\alpha }_{\mathrm{init}}$. Otherwise the preceding trial was rejected and so exceeded ${\alpha }_{\max }$, since part (ii) guarantees acceptance below that threshold. The accepted step, reduced by the factor $\tau$, consequently exceeds $\tau {\alpha }_{\max }$. In either case,The angle condition $-g_k^T{p}_k\ge \epsilon \Vert g_k{\Vert }_2\Vert p_k{\Vert }_2$ gives ${\alpha }_{\max }\ge 2\epsilon (1-\beta )\Vert g_k{\Vert }_2/(\gamma \Vert p_k{\Vert }_2)$. Hence the backtracking line search has the lower bound

From $B_k{p}_k=-g_k$ and the matrix 2-norm bound, $\Vert g_k{\Vert }_2\le {\kappa }_b\Vert p_k{\Vert }_2$, so $\Vert p_k{\Vert }_2\ge \Vert g_k{\Vert }_2/{\kappa }_b$. Let $C=2\tau \epsilon (1-\beta )/\gamma >0$. The summand in the supplied decrease bound satisfiesThus, with the fixed positive constant $c=\beta \epsilon \min ({\alpha }_{\mathrm{init}}/{\kappa }_b,C)$,The nonnegative series therefore converges, forcing its terms to zero:The argument uses the given angle condition and upper bound on $\Vert B_k{\Vert }_2$; it does not need an additional uniform lower eigenvalue bound on $B_k$.

Use the optimization Lagrangian $\mathcal{L}(x,y)=f(x)-y^{T}c(x)$, with $A(x)=Dc(x)$ and $H(x,y)={\nabla }_{xx}^2\mathcal{L}(x,y)$. Under the usual linear independence constraint qualification, a local minimizer has multipliers $y_{\ast }$ satisfying the equality-constrained Karush-Kuhn-Tucker conditionsThe second-order necessary condition isOnly feasible tangent directions enter this last test; the Hessian matrix need not be positive semidefinite on the entire ambient space.

With $y=-c/\mu$, differentiation of the quadratic penalty method objective givesThe primal Newton system is thereforeFor bounded multipliers and regular derivatives, $H$ remains of ordinary size while ${\mu }^{-1}{A}^{T}A$ becomes large in directions normal to the constraints. If the feasible tangent space is nontrivial and retains ordinary-size curvature, this produces a large condition number as $\mu \downarrow 0$.

Introduce $w=Ap/\mu$. The same step is obtained from the symmetric augmented systemEliminating $w$ recovers the primal system exactly. This matrix contains no explicit ${\mu }^{-1}$ entries and tends to a saddle-point matrix. Under full row rank of the limiting constraint Jacobian matrix and nonsingularity of the limiting reduced Hessian, that limiting matrix is nonsingular, so this formulation avoids the artificial penalty-induced separation of scales.

These conditioning statements require regularity beyond bounded $y$ alone. A rank-deficient constraint Jacobian matrix or singular reduced Hessian can still make the augmented system ill-conditioned. If $m=n$ and $A$ is nonsingular, all directions may scale together in the primal system, so a large condition number is not automatic there either.

Complete the square in the augmented Lagrangian method objective:For fixed $u,\mu$, the last term is constant in $x$. Thus its minimization is the quadratic penalty method applied to the shifted constraint $c(x)=\mu u$.

Its practical advantage is that multiplier updates, typically $u^{+}=u-c(x)/\mu$, can drive the original constraint residual to zero while $\mu$ stays moderately sized. The plain quadratic penalty method generally needs $\mu$ very small to produce the same feasibility accuracy, exposing the conditioning problem in part (ii).

The linearized constraints may be inconsistent, for example when $-c(x_k)$ is not in the image of a linear map defined by $A(x_k)$. Even if feasible, the quadratic objective may have negative curvature along its feasible tangent space, making it unbounded below and leaving no local minimizer.

Fletcher's Sℓ1QP replaces the hard linearized constraints by an L1 exact penalty and bounds the step inside a trust region. One usual form isThe zero step is always admissible. The trust region is compact, so a minimum exists even when the model Hessian matrix is indefinite; the constraint penalty measures infeasibility rather than making the subproblem impossible to solve. Adjusting its weight promotes feasibility, while the radius keeps the local model within a controlled neighbourhood. Introducing nonnegative positive- and negative-part slack variables converts the absolute values to an equivalent smooth quadratic program with box constraints.^[ref]

Stationarity of the equality-constrained quadratic model gives $g_k+H_k{s}_k-A_k^T{y}_{k+1}=0$, while feasibility gives $A_k{s}_k=-c_k$. A symmetric system linking the step and new multiplier isLet $Z$ have columns forming a basis of $\ker A_k$. The saddle-point matrix is nonsingular exactly when $A_k$ has full row rank and $Z^T{H}_{k}Z$ is nonsingular. Full row rank is necessary because a nonzero vector in $\ker A_k^T$ would give a null vector of the block matrix. Under full row rank, solving a homogeneous system forces $s=Zq$ and then $Z^T{H}_{k}Zq=0$; this proves sufficiency and the converse. Positive definiteness of the reduced Hessian is sufficient for a constrained minimum but is stronger than the nonsingularity condition asked here.

A null-space method first finds $s_p$ with $A_k{s}_p=-c_k$, then solvesand recovers $y_{k+1}$ from $A_k^T{y}_{k+1}=g_k+H_k{s}_k$. It needs the stated full row rank and nonsingular reduced matrix; Cholesky decomposition is available when the latter is positive definite.

A range-space method additionally assumes $H_k$ is invertible. It solves the Schur complement equationsIf $H_k$ is positive definite and $A_k$ has full row rank, the Schur matrix is positive definite and both solves can use Cholesky decomposition. An indefinite but invertible $H_k$ instead requires a suitable indefinite factorization; the Schur matrix must remain nonsingular.

For $\mathcal{L}=f-y^{T}c$, the inequality-constrained Karush-Kuhn-Tucker conditions areUnder an appropriate constraint qualification, the second-order necessary condition is $d^{T}H(x_{\ast },y_{\ast })d\ge 0$ on the critical coneEquivalently active constraints with positive multiplier require $A_{i}d=0$, while active zero-multiplier constraints retain $A_{i}d\ge 0$. Inactive constraints impose no first-order tangent restriction. This distinction is why the relevant cone need not be the nullspace of every active row.

On the strictly feasible region $c_i(x)>0$, the logarithmic barrier objective isWith $y_i(x)=\mu /c_i(x)$, its gradient isThe positive coefficients $y_i$ are the barrier's multiplier estimates.

The equation $F_{\mu }(x,y)=0$, together with $c(x)>0$ and $y>0$, is the perturbed Karush-Kuhn-Tucker conditions: stationarity is unchanged, while exact complementary slackness becomes $c_i(x)y_i=\mu$ for each $i$. A finite limit of such points as $\mu \downarrow 0$ satisfies the original conditions by continuity.

For a strictly feasible $x$, the second block uniquely specifies $y_i=\mu /c_i(x)$. Substitution into the first block gives precisely $\nabla B_{\mu }(x)=0$. Thuswithin the stated strictly positive domain. These are stationary barrier points; they need not all be local minimizers in a nonconvex problem.

The Jacobian matrix of $F_{\mu }$ in $(x,y)$ isThe top-left block is the Hessian matrix of $f-y^{T}c$ with $y$ held fixed. In the second block, differentiating $C(x)y$ in $x$ gives $YA$ and in $y$ gives $C$. The explicit derivative with respect to $\mu$ is $(0,-e{)}^T$. Differentiate $F_{\mu }(v(\mu ))=0$ by the chain rule:All blocks are evaluated on the current central path point.

At the current exact path point $v(\mu )$, stationarity is already zero and $Cy=\mu e$. Therefore $F_{\bar{\mu }}(v(\mu ))=(0,(\mu -\bar{\mu })e{)}^T$. Its Newton correction solvesThe matrix has no explicit dependence on the barrier parameter, so comparison with part (iv) giveswhen $J$ is nonsingular. If it is singular, the displayed vector is still a solution of the correction equations, but uniqueness of the correction cannot be inferred from differentiability of the path alone.

This step is exactly the tangent predictor in the stated first-order Taylor expansion. For nearby positive $\bar{\mu }$, it approximates the new path point with error $o(|\bar{\mu }-\mu |)$; a bounded second path derivative improves this to $O((\bar{\mu }-\mu {)}^2)$. Continuity also preserves strict positivity for sufficiently small parameter changes. Larger changes may require damping to keep both the constraint values and multiplier estimates positive.

Original solutions to all six questions and all thirteen printed terminal parts of the Numerical Linear Algebra and Continuous Optimisation paper sat on 9 June 2009. Source: released Mathematics papers, Part C 2007 11.pdf.

A real square orthogonal matrix satisfies $Q^{T}Q=Q{Q}^T=I$, so it preserves the Euclidean norm. The columns of $Z_h$ are orthonormal, hence $\Vert Z_{h}v{\Vert }_2=\Vert v{\Vert }_2$. The induced matrix norm therefore givesFor the left truncation, $\Vert Q_{\ell }^{T}A{\Vert }_2\le \Vert A{\Vert }_2$, and equality is guaranteed for arbitrary $A$ only when $\boxed{\ell =m}$. Indeed, if $\ell <m$, take $A=q_{\ell +1}{e}_1^T$, where $q_{\ell +1}$ is an omitted column of $Q$. Then $\Vert A{\Vert }_2=1$ but $Q_{\ell }^{T}A=0$. Merely taking $\ell \ge n$ does not guarantee preservation for an arbitrary $A$.

A full QR decomposition of a real $m\times n$ matrix is $A=QR$, with $Q$ orthogonal and $R$ upper trapezoidal; for a square matrix, $R$ is upper triangular. Fix the following sign convention for a Givens rotation: it is the identity except for its $(i,k)$ blockThe block transpose times the block is $I_2$ because $c^2+s^2=1$, proving orthogonality. Writing $r=\sqrt{x_i^2+x_k^2}$, the choice $\boxed{c=x_i/r,s=x_k/r}$ gives $(Jx{)}_k=0$ and $(Jx{)}_i=r$. To zero the other coordinate take $\boxed{c=x_k/r,s=-x_i/r}$, giving $(Jx{)}_i=0$. If $r=0$, both entries are already zero and the identity suffices.

For a tridiagonal matrix, successively rotate adjacent row pairs $(1,2),(2,3),\dots ,(n-1,n)$, each time zeroing the current entry $(k+1,k)$. Earlier columns remain triangular because both affected rows are already zero there. The result $R=J_{n-1}\cdots J_{1}A$ is upper triangular, so $Q=J_1^T\cdots J_{n-1}^T$ gives its QR decomposition. Only $n-1$ rotations are required; the resulting $R$ has at most two superdiagonals.

For a full square matrix, process columns $j=1,\dots ,n-1$. Within column $j$, for $i=n,n-1,\dots ,j+1$, rotate rows $(i-1,i)$ to eliminate the current $(i,j)$ entry. Previously cleared columns stay clear. With $J_1,\dots ,J_s$ in application order, $s=n(n-1)/2$, one has $\boxed{R=J_s\cdots J_{1}A,Q=J_1^T\cdots J_s^T}$. This is QR decomposition by Givens rotations.

The Gauss-Seidel method updates coordinates successively using the newest available entries:It is well-defined for every initial vector exactly when every $a_{ii}\ne 0$, equivalently when the triangular matrix $D+L$ is invertible. Nonsingularity of $A$ alone does not ensure this. The matrix form isIf the iterates converge to $z$, passing to the limit gives $(D+L)z=b-Uz$, hence $Az=b$. Nonsingularity identifies $z$ as the unique solution $x$.