Past exam of the numerical linear algebra and continuous optimisation course of the University of Oxford / 2007 (C12.1, DMAC 2488)

Original solutions to both sections, all eight questions and all thirty-two terminal parts of the Numerical Linear Algebra and Continuous Optimization paper sat on 5 June 2007. Source: released Mathematics papers, Part C 2007 11.pdf. The relaxed-Jacobi printing error and omitted zero-point qualifications are addressed where they affect the arguments.

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.

At the current point $x$, a trust-region method approximately minimizes a local quadratic modelwhere $B$ is a symmetric Hessian approximation and $\Delta >0$ is the trust-region radius. For a step with positive predicted decrease, compare the true and predicted reductions throughChoose fixed thresholds $0<{\eta }_1<{\eta }_2<1$. If $\rho <{\eta }_1$, the step is unsuccessful: reject it and reduce the radius. If ${\eta }_1\le \rho <{\eta }_2$, accept the successful step and normally retain the radius. If $\rho \ge {\eta }_2$, the step is very successful: accept it and permit an increase in radius, especially when the step reaches its boundary. The ratio tests whether the local model is accurate enough to justify moving. An approximate subproblem solver should deliver at least the model decrease achieved at the Cauchy point.

The Cauchy point minimizes the model along the negative-gradient ray within the ball. For $g\ne 0$, write $s=-tg$, $0\le t\le \Delta /\Vert g{\Vert }_2$. Thenwhose scalar minimum occurs at $t=1$ unless the ball cuts the ray short. ThereforeIf $\Vert g{\Vert }_2\le \Delta$, its predicted decrease is $\frac{1}{2}\Vert g{\Vert }_2^2$. Otherwise it is $\Delta \Vert g{\Vert }_2-\frac{1}{2}{\Delta }^2\ge \frac{1}{2}\Delta \Vert g{\Vert }_2$. Thus every feasible step with $m(s)\le m(s_c)$ satisfies the Cauchy decrease boundFor $g=0$, take $s_c=0$; the bound reduces to zero.

The global optimality conditions for a quadratic trust-region subproblem with Hessian $B$ are: for some $\lambda \ge 0$,Here $B=I$, so the shifted matrix is positive definite. If $\Vert g{\Vert }_2\le \Delta$, choose $\lambda =0$ and $s=-g$. If $\Vert g{\Vert }_2>\Delta$, the active boundary gives $s=-g/(1+\lambda )$ and $1+\lambda =\Vert g{\Vert }_2/\Delta$. HenceThis is exactly the Cauchy point. Because the model is strictly convex, it is the unique global minimizer, including $s_{\ast }=0$ when $g=0$.

The true quadratic has $g(x)=c+Hx$ and the exact expansionSubtracting the model with Hessian $I$ yields $f(x+s)-m(s)=\frac{1}{2}{s}^T(H-I)s$. Since $H$ is symmetric positive definite and $\Vert H{\Vert }_2=1$, its eigenvalues lie in $(0,1]$, so $\Vert H-I{\Vert }_2\le 1$. ThereforeIn fact the stated assumptions give a factor one-half better than the requested bound.

Write $D=f(x)-m(s)$ for the predicted decrease, assuming a nonstationary point and a step at least as good as the Cauchy point. When $\Delta \le \Vert g{\Vert }_2$, part (ii) gives $D\ge \frac{1}{2}\Vert g{\Vert }_2\Delta$. Part (iv), using its weaker requested bound, then givesConsequently $\Delta \le \frac{1}{2}(1-\eta )\Vert g{\Vert }_2$ implies $\boxed{\rho \ge \eta }$.

For the radius conclusion, use the standard update rules from part (i). Suppose no subsequence of the generated gradient norms tends to zero. Then some $\epsilon >0$ bounds those norms below eventually. If ${\Delta }_k\to 0$, eventually ${\Delta }_k\le \frac{1}{2}(1-{\eta }_2)\epsilon$, so every subsequent step is very successful by the preceding estimate. The radius is then maintained or increased, contradicting convergence to zero. Hence a radius tending to zero forces a subsequence of gradient norms to tend to zero. The conclusion uses the trust-region update rules as well as the model estimates.

Use the optimization Lagrangian $L(x,y)=f(x)-y^{T}c(x)$ and let $A=Dc$ be the $m\times n$ constraint Jacobian matrix. Under the standard regularity hypothesis $\mathrm{rank}A(x_{\ast })=m$, a local constrained minimizer satisfies, for some multiplier $y_{\ast }$,The second-order necessary condition isRegularity makes each such tangent vector the velocity of a feasible smooth curve. Twice differentiating the objective along that curve and using stationarity and the differentiated constraints gives the displayed nonnegative second derivative.

The regularity hypothesis is needed for these multiplier conditions to be necessary in this form. For example minimizing $f(x)=x$ subject to $c(x)=x^2=0$ has the unique feasible minimizer zero, but $f^{\prime }(0)-y{c}^{\prime }(0)=1$ for every $y$. The question does not explicitly state the qualification, so it must be supplied when using the usual KKT conditions.

The sequential quadratic programming step minimizes a quadratic local objective model while enforcing the first-order constraint expansion $c(x_k+s)\approx c_k+A_{k}s=0$. When the iterate and multipliers are close to a regular solution and $B_k$ approximates the Lagrangian Hessian, the resulting correction approximates a Newton step for the nonlinear KKT conditions.

Full row rank of $A_k$ ensures that $A_{k}s=-c_k$ is solvable for every $c_k$ and that the multiplier of the quadratic subproblem is unique. A positive-definite $B_k$ makes the objective strictly convex and coercive on the affine feasible set, giving a unique global subproblem minimizer. In fact positivity on $\ker A_k$, rather than on every vector, is enough for uniqueness and boundedness on that set; full positive definiteness is a convenient stronger assumption.

Stationarity of the equality-constrained quadratic subproblem gives $B_k{s}_k+g_k-A_k^T{y}_{k+1}=0$. Together with feasibility, this isMultiplying the bottom block row by $-1$ gives a symmetric indefinite saddle-point matrix, with right side $(-g_k,c_k{)}^T$.

One solution strategy uses the Schur complement. From $s_k=B_k^{-1}(A_k^T{y}_{k+1}-g_k)$ obtainThis is an $m\times m$ positive-definite system when $A_k$ has full row rank; factor $B_k$ and solve with it rather than forming its inverse.

A second strategy uses a nullspace basis $Z$ for $\ker A_k$. Choose one particular $s_p$ with $A_k{s}_p=-c_k$, write $s_k=s_p+Zz$, and solve the reduced systemThe reduced Hessian is positive definite. These range-space and nullspace reductions exploit different relative sizes of $m$ and $n-m$.

Since $p$ is even, differentiating $c_i^p$ givesAt a stationary point define $y_i=-c_i(x{)}^{p-1}/\mu$. Then $\nabla f(x)-Dc(x{)}^{T}y=0$, precisely the stationarity equation for the optimization Lagrangian. Bounded multipliers giveThus along any subsequence of stationary points with $x_{\mu }\to x_{\ast }$, a further subsequence of bounded multipliers converges to $y_{\ast }$, and continuity yields $c(x_{\ast })=0$ and $\nabla f(x_{\ast })-Dc(x_{\ast }{)}^T{y}_{\ast }=0$. These are the equality-constrained first-order conditions. Bounded multipliers alone do not establish convergence of the primal points; the precise conclusion concerns their finite cluster points. This is the power penalty multiplier limit.

Suppress the iteration subscript. From the quadratic subproblem, $As=-c$ and $g=-Bs+A^{T}y$. ThereforeWhen $|c^{T}y|>0$, the stated bound on $\mu$ makes ${\mu }^{-1}\Vert c{\Vert }_p^p\ge |c^{T}y|$. Consequentlyusing positive definiteness of $B$. Thus every nonzero such SQP step is a descent direction for the penalty function. If $c^{T}y=0$, the quotient bound is unnecessary and the same conclusion follows for every $\mu >0$. A zero step has zero directional derivative and is not a strict descent direction; it corresponds to an already feasible stationary point. This nonzero-step qualification is implicit in the requested conclusion.

For $g=\nabla f(x)\ne 0$, the unit steepest descent direction is $-g/\Vert g{\Vert }_2$ by the Cauchy-Schwarz inequality. The usual unnormalized direction is $\boxed{p=-g}$; choosing its scale separately as a line-search step gives the same ray. It satisfies $\boxed{g^{T}p=-\Vert g{\Vert }_2^2<0}$, so it is a descent direction. The first-order Taylor expansion then guarantees decrease for sufficiently small positive steps. The subsequent bounds use this unnormalized negative-gradient convention.

A Lipschitz gradient with constant $\gamma >0$ gives the descent lemmaFor $\alpha >0$, this upper bound is at most $f(x)+\alpha \beta g^{T}p$ whenever $\frac{\gamma }{2}\alpha \Vert p{\Vert }_2^2\le (\beta -1)g^{T}p$. Since $p$ is a descent direction, the right side is positive. Thus the Armijo condition holds throughoutAt $\alpha =0$ it holds by equality. If the gradient Lipschitz constant can be taken as zero, the objective is affine and every positive step along a descent direction passes this sufficient-decrease test.

The backtracking line search starts with ${\alpha }_{\mathrm{init}}>0$, tests the Armijo condition, and repeatedly replaces a failed step by $\tau \alpha$, with fixed $\tau \in (0,1)$. It accepts the first step that passes. Backtracking makes an overambitious trial step small enough for a local linear model to predict the change reliably; part (ii) ensures finite termination. The Armijo condition guarantees a decrease of at least $-\alpha \beta g^{T}p>0$ along a nonzero descent direction. Merely obtaining some decrease without this quantitative lower bound would not provide the convergence estimate in the following parts.

For $p=-g\ne 0$, part (ii) simplifies to ${\alpha }_{\max }=2(1-\beta )/\gamma$. If the initial step is accepted, ${\alpha }_k={\alpha }_{\mathrm{init}}$. Otherwise the immediately preceding trial ${\alpha }_k/\tau$ failed and hence exceeded ${\alpha }_{\max }$. ThereforeThe fixed global Lipschitz constant and search parameters make this lower bound uniform across iterations. At $g=0$, the zero update passes the first trial and the sequence can be kept stationary.

The Armijo condition with $p_k=-g(x_k)$ gives $f(x_k)-f(x_{k+1})\ge \beta {\alpha }_k\Vert g(x_k){\Vert }_2^2$. Summing these inequalities telescopes toIf $f\ge f_{\inf }$, these partial sums of nonnegative terms are bounded above by $(f(x_0)-f_{\inf })/\beta$. The resulting series converges, so its terms tend to zero: $\boxed{{\alpha }_k\Vert g(x_k){\Vert }_2^2\to 0}$.

Use the uniform lower step bound from part (iv):Thus $\boxed{g(x_k)\to 0}$. This is convergence to stationarity in the usual global convergence theorem for the steepest descent method. It does not itself prove that the iterates converge, or that any stationary limit is a global minimizer of a nonconvex objective.

Two reasons to prefer more advanced methods for smooth, high-accuracy optimization are:

The question's blanket wording is too strong: the steepest descent method and related gradient algorithms are used in practice, especially when problem size, cheap gradients or simplicity make curvature-based alternatives costly. The objections describe limitations, not universal nonuse.

With $L(x,y)=f(x)-y^{T}c(x)$, the Karush-Kuhn-Tucker conditions areUnder a suitable constraint qualification, such as the linear independence constraint qualification on active constraint gradients, they are necessary for a local minimum. Smoothness alone is insufficient. The sign convention here corresponds to constraints written as $c_i\ge 0$ and a minus sign in the optimization Lagrangian.

Write the constraint as $c(x)=1-x^{T}x\ge 0$ and its multiplier as $y\ge 0$. The KKT conditions becomeIf $A$ is positive definite, $A+yI$ is positive definite for every $y\ge 0$, so $\boxed{x=0,y=0}$ is the only critical point.

For indefinite $A$, the origin remains a critical point, contrary to an unqualified reading of the question. Every nonzero critical point is an eigenvector with $\boxed{Ax=-yx}$, so its eigenvalue is nonpositive. A negative-eigenvalue vector must lie on the unit sphere, with multiplier equal to the negative of that eigenvalue. A nonzero vector in $\ker A$ may have any norm at most one and has multiplier zero. Positive-eigenvalue directions do not satisfy the multiplier sign condition. These first-order critical points include saddles; the global minimum selects a unit vector in the lowest-eigenvalue eigenspace.

The logarithmic barrier isIts gradient is $2Ax+2\mu x/(1-x^{T}x)$. Thus every nonzero stationary point satisfiesA fixed negative-eigenvalue branch exists for $0<\mu <-{\lambda }_1$. The zero point is stationary for every $\mu$ as well, and it is not an eigenvector; at zero the displayed scalar $-\mu$ need not be an eigenvalue of $A$. Thus the claimed eigenvalue conclusion requires the nonzero-point qualification. This distinction is essential for the following Hessian analysis.

Put $d=1-x^{T}x>0$. Differentiating the barrier gradient givesAt a nonzero stationary point on the branch $Ax={\lambda }_{1}x$, $\mu /d=-{\lambda }_1$. Choose an orthonormal eigenbasis of the symmetric matrix $A$ whose first vector is $x/\Vert x{\Vert }_2$. On each of the remaining eigenvectors the rank-one term vanishes, so the corresponding Hessian eigenvalues areAlong $x$, the first two terms cancel and the remaining eigenvalue isIf ${\lambda }_1$ repeats, a transverse eigenvalue can be zero. If a lower eigenvalue than ${\lambda }_1$ exists, a transverse eigenvalue is negative, identifying this stationary branch as a saddle rather than a minimum.

Along a fixed nonzero branch, $\Vert x(\mu ){\Vert }_2^2=1+\mu /{\lambda }_1\to 1$ and the radial Hessian eigenvalue grows like $4{\lambda }_1^2/\mu$. The nonzero transverse eigenvalues $2({\lambda }_i-{\lambda }_1)$ stay fixed. If ${\lambda }_1$ is simple, the ratio of the largest to smallest absolute eigenvalue therefore tends to infinity. If it repeats, the Hessian is already singular because of the zero transverse directions; there is no finite spectral condition number to begin with. These conclusions concern the nonzero branches: at the stationary origin $H=2(A+\mu I)$, and the same blowup is not implied.

Ill-conditioning does not invalidate barrier minimizers as approximations to the constrained solution. Let ${\lambda }_{\min }<0$. At radius $r$, the smallest possible quadratic value is ${\lambda }_{\min }{r}^2$, so global barrier minimization reduces to minimizing ${\lambda }_{\min }t-\mu \log (1-t)$ for $0\le t<1$. For $0<\mu <-{\lambda }_{\min }$ its unique minimizing radius satisfiesThus every limiting unit vector in that eigenspace solves the original problem, with objective ${\lambda }_{\min }$. For $\mu \ge -{\lambda }_{\min }$ the barrier minimizer is the origin instead. The increasing condition number concerns the numerical difficulty of small-parameter subproblems, which can call for scaling or a primal-dual interior-point method; it does not prevent the mathematical convergence just established. This is the quadratic-ball barrier eigenvalue branch.