Numerical Methods Qualification Exam Problems and Solutions (University of Maryland)/January 2008

Problem 1

Consider the system $Ax=b$. The GMRES method starts with a point $x_0$ and normalizes the residual $r_0=b-A{x}_0$ so that $v_1=\frac{r_0}{\nu }$ has 2-norm one. It then constructs orthonormal Krylov bases $V_k=(v_1{v}_2\cdots v_m)$ satisfying $A{V}_k=V_{k+1}{H}_k$ where $H_k$ is a $(k+1)\times k$ upper Hessenberg matrix. One then looks for an approximation to $x$ of the form $x(c)=x_0+V_{k}c$ choosing $c_k$ so that $\Vert r(c)\Vert =\Vert b-Ax(c)\Vert$ is minimized, where $\Vert \cdot \Vert$ is the usual Euclidean norm.

Problem 1a

Show that $c_k$ minimizes $\Vert \nu e_1-H_{k}c\Vert$.

Solution 1a

We wish to show that $\Vert b-Ax(c)\Vert =\Vert \nu e_1-H_{k}c\Vert$

$\begin{array}{cc}\Vert b-Ax(c)\Vert & =\Vert b-A(x_0+V_{k}c)\Vert \\ & =\Vert b-A{x}_0-A{V}_{k}c\Vert \\ & =\Vert r_0-A{V}_{k}c\Vert \\ & =\Vert r_0-V_{k+1}{H}_{k}c\Vert \\ & =\Vert \nu v_1-V_{k+1}{H}_{k}c\Vert \\ & =\Vert V_{k+1}\underset{h_c}{\underbrace{(\nu e_1-H_{k}c)}}\Vert \\ & =(V_{k+1}{h}_c,V_{k+1}{h}_c{)}^{\frac{1}{2}}\\ & =((V_{k+1}{h}_c{)}^T{V}_{k+1}{h}_c{)}^{\frac{1}{2}}\\ & =(h_c^T{V}_{k+1}^T{V}_{k+1}{h}_c{)}^{\frac{1}{2}}\\ & =(h_c^T{h}_c{)}^{\frac{1}{2}}\\ & =\Vert h_c\Vert \\ & =\Vert \nu e_1-H_{k}c\Vert \end{array}$

Problem 1b

Suppose we choose to solve the least squares problem in part a for $c_k$ by the method of orthogonal triangularization (QR). What is the order of the floating point operation for this method. Give reasons.

Solution 1b

We would like to transform $H_k$, a $(k+1)\times k$ upper Hessenberg matrix, into QR form.

The cost is on the order of $4k^2+\frac{1}{2}{k}^2=4.5k^2$ from the cost of Given's Rotations and backsolving.

Given's Rotations Cost

We need $k$ Given's Rotation multiplies to zero out each of the $k$ subdiagonal entries and hence transform $H_k$ into upper triangular form $R$,

Each successive Given's Rotations multiply on an upper Hessenberg matrix requires four fewer multiplies because each previous subdiagonal entry has been zeroed out by a Given's Rotation multiply.

Hence the cost of $k$ Given's Rotations multiplies is

$4k+4(k-1)+\dots +4\cdot 1<4k\cdot k=4k^2$

Back Solving Cost

$R$ is a $(k+1)\times k$ upper triangular matrix with last row zero. Hence, we need to backsolve $k$ upper triangular rows.

$1+2+\dots +k=\frac{k(k+1)}{2}=\frac{k^2}{2}+\frac{k}{2}$

Problem 2

We wish to approximate the integral $I(f)\equiv \underset{a}{\overset{b}{\int }}f(x)dx$

Problem 2a

State the composite trapezoidal rule $Q_{T,n}$ for approximating $I(f)$ on a uniform partitioning of width $h=\frac{b-a}{n}$, and give a formula for the error $I(f)-Q_{T,n}$ that is in a form suitable for extrapolation.

Solution 2a

The composite trapezoidal rule is

$Q_{T,n}(f)=\frac{h}{2}(f(a)+f(b))+h\sum _{k=1}^{n-1}f(x_k)$

The error is

$I(f)-Q_{T,n}(f)=-\frac{(b-a)f^{\prime \prime }(\xi )}{12}{h}^2$

where $\xi \in [a,b]$.

Derivation of Composite Trapezoidal Error

The local error, the error on one interval, is

$-\frac{h^3}{12}{f}^{\prime \prime }(\eta )$.

Observe that

$n\underset{\eta \in [a,b]}{\min }{f}^{\prime \prime }({\eta }_i)\le \sum _{i=1}^n{f}^{\prime \prime }({\eta }_i)\le n\underset{\eta \in [a,b]}{\max }{f}^{\prime \prime }({\eta }_i)$

which implies

$\underset{\eta \in [a,b]}{\min }{f}^{\prime \prime }({\eta }_i)\le \sum _{i=1}^n\frac{f^{\prime \prime }({\eta }_i)}{n}\le \underset{\eta \in [a,b]}{\max }{f}^{\prime \prime }({\eta }_i)$

Hence, the Intermediate Value Theorem implies there exists a $\xi \in [a,b]$ such that

$f^{\prime \prime }(\xi )=\sum _{i=1}^n\frac{f^{\prime \prime }({\eta }_i)}{n}$.

Multiplying both sides of the equation by $n$,

$n{f}^{\prime \prime }(\xi )=\sum _{i=1}^n{f}^{\prime \prime }({\eta }_i)$

Using this relationship, we have

$\begin{array}{cc}I(f)-Q_{T,n}(f)& =\sum _{i=1}^n{I}_i(f)-Q_i(f)\\ & =\sum _{i=1}^n[-\frac{h^3}{12}{f}^{\prime \prime }({\eta }_i)]\\ & =-\frac{h^3}{12}{f}^{\prime \prime }(\xi )n\\ & =-\frac{h^3}{12}{f}^{\prime \prime }(\xi )(\frac{b-a}{h})\\ & =-\frac{h^2}{12}{f}^{\prime \prime }(\xi )(b-a)\end{array}$

Derivation of Local Error

The error in polynomial interpolation can be found by using the following theorem:

 Assume $f^{n+1}$ exists on $[a,b]$ and $\{x_0,x_1,\dots ,x_n|x\in [a,b]\}.$ $P_n$ interpolates $f$ at $\{x_j{\}}_{j=0}^n$.  
 Then there is a ${\xi }_x$ ($\xi$ is dependent on $x$) such that 
 
 $f(x)-p_n(x)=\frac{(x-x_0)(x-x_1)...(x-x_n)}{(n+1)!}{f}^{(n+1)}({\xi }_x)$
 
 where ${\xi }_x$ lies in $(m,M)$,  $m=\min \{x_0,x_1,\dots ,x_n,X\}$,  $M=\max \{x_0,x_1,\dots ,x_n,X\}$

Applying the theorem yields,

$f(x)-p_1(x)=(x-a)(x-b)\frac{f^{\prime \prime }({\xi }_x)}{2}$

Hence,

$\begin{array}{cc}E(f)& =\underset{a}{\overset{b}{\int }}f(x)-p_1(x)dx\\ & =\underset{a}{\overset{b}{\int }}\underset{w(x)}{\underbrace{\underset{>0}{\underbrace{(x-a)}}\underset{<0}{\underbrace{(x-b)}}}}\frac{f^{\prime \prime }({\xi }_x)}{2}dx\end{array}$

Since $a$ is the start of the interval, $x-a$ is always positive. Conversely, since $b$ is the end of the interval, $x-b$ is always negative. Hence, $w(x)<0$ is always of one sign. Hence from the mean value theorem of integration, there exists a $\zeta \in [a,b]$ such that

$E(f)=\frac{f^{\prime \prime }(\zeta )}{2}\underset{a}{\overset{b}{\int }}(x-a)(x-b)dx$

Note that $f^{\prime \prime }(\zeta )$ is a constant and does not depend on $x$.

Integrating $\underset{a}{\overset{b}{\int }}(x-a)(x-b)dx$, yields,

$\begin{array}{cc}E(f)& =\frac{f^{\prime \prime }(\zeta )}{2}(-\frac{(b-a)}{6})\\ & =-\frac{f^{\prime \prime }(\zeta )(b-a)}{12}\end{array}$

Problem 2b

Use the error formula to derive a new quadrature rule obtained by performing one step of extrapolation on the composite trapezoidal rule. What is this rule, and how does its error depend on $h$? You may assume here that $f$ is as smooth as you need it to be.

Solution 2b

STOER AND BUESCH pg 162

INTRODUCTION TO APPLIED NUMERICAL ANALYSIS by RICHARD HAMMING pg 178

The error for the composite trapezoidal rule at $n$ points in $[a,b]$ is

$E_n=I(f)-Q_{T,n}=\frac{-(b-a)h^2}{12}+𝒪(h^3)$

With $2n$ points, there are twice as many intervals, but the intervals are half as wide. Hence, the error for the composite trapezoidal rule at $2n$ points in $[a,b]$ is

$E_{2n}=I(f)-Q_{T,2n}=2\cdot \frac{-(b-a)(\frac{h}{2}{)}^2}{12}+𝒪(h^3)=\frac{-(b-a)h^2}{24}+𝒪(h^3)$

We can eliminate the $h^2$ term by choosing using an appropriate linear combination of $E_n$ and $E_{2n}$. This gives a new error rule with $h^3$error.

$\begin{array}{cc}{E}_n-2E_{2n}& =\frac{-(b-a)h^2}{12}-2\cdot \frac{-(b-a)h^2}{24}+𝒪(h^3)\\ & =𝒪(h^3)\end{array}$

Substituting our equations for $E_n$ and $E_{2n}$ on the left had side gives

$I(f)-Q_{T,n}-2I(f)-2Q_{T,2n}=𝒪(h^3)$

$I(f)=2Q_{T,2n}-Q_{T,n}+𝒪(h^3)$

If we call our new rule $\hat{Q}$ we have

$\hat{Q}=2Q_{T,2n}-Q_{T,n}$ whose error is on the order of $𝒪(h^3)$

Problem 3

Consider the shifted QR iteration for computing a 2 x 2 matrix $A$: starting with $A_0=A$, compute $A_k-{\sigma }_{k}I=Q_k{R}_k{A}_{k+1}=R_k{Q}_k+{\sigma }_{k}I$ where ${\sigma }_k$ is a scalar shift.

Problem 3a

If $A_k=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right],$ specify the orthogonal matrix $Q_k$ used to perform this step when Givens rotations are used. The matrix should be described in terms of the entries of the shifted matrix.

Solution 3a

Let ${\tilde{A}}_k$ be the ${\sigma }_k$-shifted matrix of $A_k$ i.e.

${\tilde{A}}_k=A_k-{\sigma }_{k}I=\left[\begin{array}{cc}{a}_{11}-{\sigma }_k& a_{12}\\ a_{21}& a_{22}-{\sigma }_k\end{array}\right],$

$Q_k^T$ is an orthogonal 2x2 Given's rotations matrix. $Q_k^T$'s entries are chosen such that when $Q_k^T$ is pre-multiplied against the 2x2 matrix ${\tilde{A}}_k$, $Q_k^T$ will zero out the (2,1) entry of ${\tilde{A}}_k$ and scale ${\tilde{A}}_k$'s remaining three entries i.e.

$Q_k^T{\tilde{A}}_k=\left[\begin{array}{cc}*& *\\ 0& *\end{array}\right]=R_k$

where $*$ denotes calculable scalar values we are not interested in and $R_k$ is our desired upper triangular matrix sought by the QR algorithm.

Since $Q_k^T$ is orthogonal, the above equation implies

$\begin{array}{cc}{Q}_k^T{\tilde{A}}_k& =R_k\\ Q_k{Q}_k^T{\tilde{A}}_k& =Q_k{R}_k\\ {\tilde{A}}_k& =Q_k{R}_k\end{array}$

The Given's rotation $Q_k^T$ is given by

$Q_k^T=\left[\begin{array}{cc}c& s\\ -s& c\end{array}\right]$

where

$\begin{array}{cc}c& =\frac{a_{11}-{\sigma }_k}{r}\\ s& =\frac{a_{21}}{r}\\ r& =\sqrt{(a_{11}-{\sigma }_k{)}^2+a_{21}^2}\end{array}$

Taking the transpose of $Q_k^T$ yields

$Q_k=\left[\begin{array}{cc}c& -s\\ s& c\end{array}\right]$

Problem 3b

Suppose $a_{21}=\delta$, a small number, and $|a_{12}|\le (a_{11}-a_{22}{)}^2$. Demonstrate that with an appropriate shift ${\sigma }_k$, the (2,1)-entry of $A_{k+1}$ is of magnitude $O({\delta }^2)$. What does this suggest about the convergence rate of the QR iteration?

Solution 3b

From hypothesis and part (a),

$\begin{array}{cc}{A}_{k+1}& =R_k{Q}_k+{\sigma }_{k}I\\ & =\left[\begin{array}{cc}*& *\\ 0& r_{22}\end{array}\right]\left[\begin{array}{cc}c& -s\\ s& c\end{array}\right]+\left[\begin{array}{cc}{\sigma }_k& 0\\ 0& {\sigma }_k\end{array}\right]\end{array}$

Let $A_{k+1}(2,1)$ be the (2,1) entry of $A_k$. Using the above equation, we can find $A_{k+1}(2,1)$ by finding the inner product of the second row of $R_k$ and first column $Q_k$ and adding the (2,1) entry of ${\sigma }_{k}I$ i.e.

$\begin{array}{cc}{A}_{k+1}(2,1)& =(0\cdot c+r_{22}\cdot s)+0\\ & =r_{22}s\end{array}$

We need to find the value of $r_{22}$ so we need to calculate $R_k$.

From hypothesis and the orthogonality of $Q_k$, we have

$\begin{array}{cc}{A}_k-{\sigma }_{k}I& =Q_k{R}_k\\ Q_k^T(A_k-{\sigma }_{k}I)& =Q_k^T{Q}_k{R}_k\\ Q_k^T(A_k-{\sigma }_{k}I)& =R_k\\ R_k& =Q_k^T(A_k-{\sigma }_{k}I)\\ & =\left[\begin{array}{cc}c& s\\ -s& c\end{array}\right]\left[\begin{array}{cc}{a}_{11}-{\sigma }_k& a_{12}\\ \delta & a_{22}-{\sigma }_k\end{array}\right]\end{array}$

From $R_k$, we can find its (2,2) entry $r_{2}2$ by using inner products.

$r_{22}=-s\cdot a_{12}+c\cdot (a_{22}-{\sigma }_k)$

Now that we have $r_{22}$ we can calculate $A_{k+1}(2,1)$ by using appropriate substitutions

$\begin{array}{cc}{A}_{k+1}(2,1)& =r_{22}\cdot s\\ & =(-s{a}_{12}+c(a_{22}-{\sigma }_k))s\\ & =-s^2{a}_{12}+cs(a_{22}-{\sigma }_k)\\ & =\frac{-{\delta }^2{a}_{12}}{(a_{11}-{\sigma }_k{)}^2+{\delta }^2}+\frac{(a_{11}-{\sigma }_k)(a_{22}-{\sigma }_k)\delta }{(a_{11}-{\sigma }_k{)}^2+{\delta }^2}\\ & =\frac{-{\delta }^2{a}_{12}+\delta (a_{11}-{\sigma }_k)(a_{22}-{\sigma }_k)}{(a_{11}-{\sigma }_k{)}^2+{\delta }^2}\\ & \approx \frac{-{\delta }^2{a}_{12}+\delta (a_{11}-{\sigma }_k)(a_{22}-{\sigma }_k)}{(a_{11}-{\sigma }_k{)}^2}\end{array}$

since $\delta$ is a small number.

Let our shift ${\sigma }_k=a_{22}$. Then the above equation yields,

$\begin{array}{cc}{A}_{k+1}(2,1)& \approx \frac{-{\delta }^2{a}_{12}+\delta (a_{11}-{\sigma }_k)(a_{22}-{\sigma }_k)}{(a_{11}-{\sigma }_k{)}^2}\\ & \approx \frac{-{\delta }^2{a}_{12}+\delta (a_{11}-a_{22})(a_{22}-a_{22})}{(a_{11}-a_{22}{)}^2}\\ & =\frac{-{\delta }^2{a}_{12}}{(a_{11}-a_{22}{)}^2}\end{array}$

Hence,

$\begin{array}{cc}|A_{k+1}(2,1)|& \approx \left|\frac{{\delta }^2{a}_{12}}{(a_{11}-a_{22}{)}^2}\right|\\ & \le |{\delta }^2|\end{array}$

since $|a_{12}|\le (a_{11}-a_{22}{)}^2$

Hence we have shown that $A_{k+1}(2,1)$ is of order ${\delta }^2$.

If $\delta$ is small, the QR convergence rate will be quadratic.