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

Problem 1

Let $x_1,x_2,\dots ,x_n$ be distinct real numbers, and consider the matrix $V=\left(\begin{array}{ccccc}1& x_1& x_1^2& \cdots & x_1^{n-1}\\ 1& x_2& x_2^2& \cdots & x_2^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_n& x_n^2& \cdots & x_n^{n-1}\end{array}\right)$ Prove that $V^{-1}$ can be expressed as $V^{-1}=\left(\begin{array}{c}\\ \\ {\alpha }^{(1)},{\alpha }^{(2)},{\alpha }^{(3)},\cdots ,{\alpha }^{(n)}\\ \\ \end{array}\right)$ where the column vector ${\alpha }^{(j)}=({\alpha }_1^{(j)},{\alpha }_2^{(j)},\cdots ,{\alpha }_n^{(j)}{)}^T$ generates a polynomial $p_j(x)=\sum _{i=1}^n{\alpha }_i^{(j)}{x}^{i-1}$ satisfying $p_j(x)=\frac{\prod _{i=1,i\ne j}^n(x-x_i)}{\prod _{i=1,i\ne j}^n(x_j-x_i)}$ You may cite anything you know about polynomial interpolation to justify that $V$ is nonsingular.

Solution 1

We want to show that

$V{V}^{-1}=I$

or equivalently the $i$th, $j$th entry of $V{V}^{-1}$ is $1$ for $i=j$ and $0$ for $i\ne j$ i.e.

$\{V{V}^{-1}{\}}_{ij}={\delta }_{ij}=\{\begin{array}{cc}1& \text{ if }i=j\\ 0& \text{ if }i\ne j\end{array}$

First notice that

$\{V{V}^{-1}{\}}_{ij}=\left[\begin{array}{c}1+x_i^1+x_i^2+\cdots +x_i^{n-1}\end{array}\right]\left[\begin{array}{c}{\alpha }_1^{(j)}\\ {\alpha }_2^{(j)}\\ ⋮\\ {\alpha }_n^{(j)}\end{array}\right]=\sum _{k=1}^n{\alpha }_k^{(j)}{x}_i^{k-1}$

Also notice that

$p_j(x_i)=\sum _{k=1}^n{\alpha }_k^{(j)}{x}_i^{k-1}={\delta }_{ij}$

Hence,

$\{V{V}^{-1}{\}}_{ij}=p_j(x_i)={\delta }_{ij}$

Problem 2

Let $A$ be a real matrix of order $n$ with eigenvalues ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$ and $n$ linearly independent eigenvectors $v_1,v_2,\dots ,v_n$. Suppose you want to find an eigenvector $v_j$ corresponding to the eigenvalue ${\lambda }_j$ and you are given $\sigma$ such that $|{\lambda }_j-\sigma |<|{\lambda }_i-\sigma |\text{ for all }i\ne j$ Specify a numerical algorithm for finding $v_j$ and give a convergence proof for this algorithm demonstrating that it is convergent under appropriate circumstances

Solution 2

Shifted Inverse Power Method

Let

$\begin{array}{cc}\tilde{A}& =(A-\sigma I{)}^{-1}\\ \end{array}$

Then,

$\widetilde{{\lambda }_i}=\frac{1}{|{\lambda }_i-\sigma |}$

which implies

$\widetilde{{\lambda }_j}>{\tilde{\lambda }}_i\text{ for all }i\ne j$

Since shifting the eigenvalues and inverting the matrix does not affect the eigenvectors,

$\tilde{A}{v}_i=\widetilde{{\lambda }_i}{v}_{i}i=1,2,\dots ,n$

Assume $\Vert v_i\Vert =1$ for all $i$. Generate $w_0,w_1,w_2,\dots$ to find $v_j$. Start with arbitrary $w_0$ such that $\Vert w_0\Vert =1$.

For $k=0,1,2,\dots$
  ${\hat{w}}_{k+1}=\tilde{A}{w}_k$
  $w_{k+1}=\frac{{\hat{w}}_{k+1}}{\Vert {\hat{w}}_{k+1}\Vert }$
  ${\lambda }_j^{(k+1)}=\frac{({\hat{w}}_{k+1},w_k)}{(w_k,w_k)}=\frac{(\tilde{A}{w}_k,w_k)}{(w_k,w_k)}$ (Rayleigh quotient)
End

Convergence of Power Method

If $|{\lambda }_1|>|{\lambda }_i|$, for all $i\ne 1$, then $A^k{w}_0$ will be dominated by $v_1$.

Since $v_1,v_2,\dots ,v_n$ are linearly independent, they form a basis of $R^n$. Hence,

$w_0={\alpha }_1{v}_1+{\alpha }_2{v}_2+\dots +{\alpha }_n{v}_n$

From the definition of eigenvectors,

$\begin{array}{cc}A{w}_0& ={\alpha }_1{\lambda }_1{v}_1+{\alpha }_2{\lambda }_2{v}_2+\dots +{\alpha }_n{\lambda }_n{v}_n\\ A^2{w}_0& ={\alpha }_1{\lambda }_1^2{v}_1+{\alpha }_2{\lambda }_2^2{v}_2+\dots +{\alpha }_n{\lambda }_n^2{v}_n\\ & ⋮\\ A^k{w}_0& ={\alpha }_1{\lambda }_1^k{v}_1+{\alpha }_2{\lambda }_k^2{v}_2+\dots +{\alpha }_n{\lambda }_n^k{v}_n\\ \end{array}$

To find a general form of $w_k$, the approximate eigenvector at the kth step, examine a few steps of the algorithm:

$\begin{array}{cc}{\hat{w}}_1& =A{w}_0\\ w_1& =\frac{{\hat{w}}_1}{\Vert {\hat{w}}_1\Vert }=\frac{A{w}_0}{\Vert A{w}_0\Vert }\\ {\hat{w}}_2& =A{w}_1=\frac{A^2{w}_0}{\Vert A{w}_0\Vert }\\ w_2& =\frac{{\hat{w}}_2}{\Vert {\hat{w}}_2\Vert }=\frac{\frac{A^2{w}_0}{\Vert A{w}_0\Vert }}{\frac{\Vert A^2{w}_0\Vert }{\Vert A{w}_0\Vert }}=\frac{A^2{w}_0}{\Vert A^2{w}_0\Vert }\\ {\hat{w}}_3& =\frac{A^3{w}_0}{\Vert A^2{w}_0\Vert }\\ w_3& =\frac{\frac{A^3{w}_0}{\Vert A^2{w}_0\Vert }}{\frac{\Vert A^3{w}_0\Vert }{\Vert A^2{w}_0\Vert }}=\frac{A^3{w}_0}{\Vert A^3{w}_0\Vert }\\ \end{array}$

From induction,

$w_k=\frac{A^k{w}_0}{\Vert A^k{w}_0\Vert }$

Hence,

$\begin{array}{cc}{w}_k& =\frac{A^k{w}_0}{\Vert A^k{w}_0\Vert }\\ & =\frac{{\alpha }_1{\lambda }_1^k{v}_1+{\alpha }_2{\lambda }_2^k{v}_2+\dots +{\alpha }_n{\lambda }_n^k{v}_n}{\Vert A^k{w}_0\Vert }\\ & =\frac{{\lambda }_1^k({\alpha }_1{v}_1+{\alpha }_2(\frac{{\lambda }_2}{{\lambda }_1}{)}^k{v}_2+\dots +{\alpha }_n(\frac{{\lambda }_n}{{\lambda }_1}{)}^k{v}_n)}{\Vert A^k{w}_0\Vert }\end{array}$

Comparing a weighted $w_k$ and $v_1$,

$\begin{array}{cc}\Vert \frac{\Vert A^k{w}_0\Vert }{{\lambda }_1^k}{w}_k-{\alpha }_1{v}_1\Vert & =\Vert {\alpha }_2{\left(\frac{{\lambda }_2}{{\lambda }_1}\right)}^k{v}_2+\dots +{\alpha }_n{\left(\frac{{\lambda }_n}{{\lambda }_1}\right)}^k{v}_n\Vert \\ & \le {\alpha }_2{\left|\frac{{\lambda }_2}{{\lambda }_1}\right|}^k+\dots +{\alpha }_n{\left|\frac{{\lambda }_n}{{\lambda }_1}\right|}^k\end{array}$

since $\Vert v_i\Vert =1$ by assumption.

The above expression goes to $0$ as $k\to \mathrm{\infty }$ since $|{\lambda }_1|>|{\lambda }_i|$, for all $i\ne 1$. Hence as $k$ grows, $w_k$ is parallel to $v_1$. Because $\Vert w_k\Vert =1$, it must be that $w_k\to \pm v_1$.

Problem 3

Suppose A is an upper triangular, nonsingular matrix. Show that both Jacobi iterations and Gauss-Seidel iterations converge in finitely many steps when used to solve $A𝐱=𝐛$.

Solution 3

Derivation of Iterations

Let $A=D+L+U$ where $D$ is a diagonal matrix, $L,$ is a lower triangular matrix with a zero diagonal and $U$ is an upper triangular matrix also with a zero diagonal.

The Jacobi iteration can be found by substituting into $Ax=b$, grouping $L+U$, and solving for $x$ i.e.

$\begin{array}{cc}Ax& =b\\ (D+L+U)x& =b\\ Dx+(L+U)x& =b\\ Dx& =b-(L+U)x\\ x& =D^{-1}(b-(L+U)x)\end{array}$

Since $L=0$ by hypothesis, the iteration is

$x^{(i+1)}=D^{-1}(b-U{x}^{(i)})$

Similarly, the Gauss-Seidel iteration can be found by substituting into $Ax=b$, grouping $D+L$, and solving for $x$ i.e.

$\begin{array}{cc}Ax& =b\\ (D+L+U)x& =b\\ (D+L)x+Ux& =b\\ (D+L)x& =b-Ux\\ x& =(D+L{)}^{-1}(b-Ux)\end{array}$

Since $L=0$ by hypothesis, the iteration has identical from as Jacopi:

$x^{(i+1)}=D^{-1}(b-U{x}^{(i)})$

Convergence in Finite Number of Steps

Jacobi and Gauss-Seidel are iterative methods that split the matrix $A\in {ℝ}^{n\times n}$ into $D$, $U$, and $L$: Diagonal, Upper (everything above the diagonal), and Lower (everything below the diagonal) triangular matrices, respectively. Their iterations go as such

$x^{(i+1)}=(D+L{)}^{-1}(b-U{x}^{(i)})$ (Gauss-Seidel)

$x^{(i+1)}=D^{-1}(b-(U+L)x^{(i)})$ (Jacobi)

In our case $A$ is upper triangular, so $L$ is the zero matrix. As a result, the Gauss-Seidel and Jacobi methods take on the following identical form

$x^{(i+1)}=D^{-1}(b-U{x}^{(i)})$

Additionally, $x$ can be written

$x=D^{-1}(b-Ux)$

Subtracting $x$ from $x^{(i+1)}$ we get

$\begin{array}{cc}{e}^{(i+1)}& =x^{(i+1)}-x\\ & =D^{-1}(b-U{x}^{(i)})-D^{-1}(b-Ux)\\ & =D^{-1}U(-x^{(i)}+x)\\ & =D^{-1}U{e}^{(i)}\\ & =D^{-1}U(D^{-1}U{e}^{(i-1)})=(D^{-1}U{)}^2{e}^{(i-1)}\\ & ⋮\\ & =(D^{-1}U{)}^{i+1}{e}^{(0)}\end{array}$

In our problem, $D^{-1}$ is diagonal and $U$ is upper triangular with zeros along the diagonal. Notice that the product $D^{-1}U$ will also be upper triangular with zeros along the diagonal.

Let $R=D^{-1}U$, $R\in {ℝ}^{n\times n}$

$R=\left(\begin{array}{ccccc}0& *& *& \cdots & *\\ 0& 0& *& \cdots & *\\ ⋮& ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & \cdots & 0& *\\ 0& \cdots & \cdots & 0& 0\end{array}\right)$

Also, let $\tilde{R}\in {ℝ}^{n\times n}$ be the related matrix

$\begin{array}{cc}\tilde{R}& =\left(\begin{array}{ccccccc}0& \cdots & 0& *& *& \cdots & *\\ 0& \cdots & 0& 0& *& \cdots & *\\ ⋮& ⋮& ⋮& ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & 0& 0& \cdots & 0& *\\ 0& \cdots & 0& 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \cdots & ⋮& ⋮\\ 0& \cdots & 0& 0& \cdots & 0& 0\end{array}\right)\\ & =\left(\begin{array}{cc}𝐚& T\\ \mathrm{𝟎}& {𝐚}^T\end{array}\right)\end{array}$

Where $𝐚$ is $k\times (n-k)$, $T$ is $k\times k$, and $\mathrm{𝟎}$ is $(n-k)\times (n-k)$.

Finally, the product $R\tilde{R}$ (call it (1)) is

$\begin{array}{cc}R\tilde{R}& =\left(\begin{array}{ccccc}0& *& *& \cdots & *\\ 0& 0& *& \cdots & *\\ ⋮& ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & \cdots & 0& *\\ 0& \cdots & \cdots & 0& 0\end{array}\right)\left(\begin{array}{ccccccc}0& \cdots & 0& *& *& \cdots & *\\ 0& \cdots & 0& 0& *& \cdots & *\\ ⋮& ⋮& ⋮& ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & 0& 0& \cdots & 0& *\\ 0& \cdots & 0& 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \cdots & ⋮& ⋮\\ 0& \cdots & 0& 0& \cdots & 0& 0\end{array}\right)\\ & =\left(\begin{array}{cccccccc}0& \cdots & 0& 0& *& \cdots & \cdots & *\\ 0& \cdots & 0& 0& 0& *& \cdots & *\\ ⋮& ⋮& ⋮& ⋮& \cdots & \cdots & \ddots & ⋮\\ 0& \cdots & 0& 0& \cdots & \cdots & 0& *\\ 0& \cdots & 0& 0& \cdots & \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \cdots & \cdots & ⋮& ⋮\\ 0& \cdots & 0& 0& \cdots & \cdots & 0& 0\end{array}\right)\\ & =\left(\begin{array}{cc}𝐚& \tilde{T}\\ \mathrm{𝟎}& {𝐚}^T\end{array}\right)\end{array}$

Here $\tilde{T}$ is almost identical in structure to $T$, except that its diagonal elements are zeros.

At this point the convergence in $n$ steps (the size of the starting matrix) should be apparent since $R$ is just $\tilde{R}$ where $k=0$ and each time $R$ is multiplied by $\tilde{R}$, the k-th super-diagonal is zeroed out (where k=0 is the diagonal itself). After $n-1$ applications of $R$, the result will be the zero matrix of size $n$.

In brief, $R^n$ is the zero matrix of size $n$.

Therefore $e^{(n)}=(D^{-1}U{)}^n{e}^{(0)}\equiv 0$, i.e. the Jacobi and Gauss-Seidel Methods used to solve $A𝐱=b$ converge in $n$ steps when $A\in {ℝ}^{n\times n}$ is upper triangular.

Examples of (1)

$n=3$

$R=\left(\begin{array}{ccc}0& *& *\\ 0& 0& *\\ 0& 0& 0\end{array}\right)$

$R^2=\left(\begin{array}{ccc}0& *& *\\ 0& 0& *\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccc}0& *& *\\ 0& 0& *\\ 0& 0& 0\end{array}\right)=\left(\begin{array}{ccc}0& 0& *\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)$

$R^3=\left(\begin{array}{ccc}0& *& *\\ 0& 0& *\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccc}0& 0& *\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)$

$n=4$

$R=\left(\begin{array}{cccc}0& *& *& *\\ 0& 0& *& *\\ 0& 0& 0& *\\ 0& 0& 0& 0\end{array}\right)$

$R^2=\left(\begin{array}{cccc}0& *& *& *\\ 0& 0& *& *\\ 0& 0& 0& *\\ 0& 0& 0& 0\end{array}\right)\left(\begin{array}{cccc}0& *& *& *\\ 0& 0& *& *\\ 0& 0& 0& *\\ 0& 0& 0& 0\end{array}\right)=\left(\begin{array}{cccc}0& 0& *& *\\ 0& 0& 0& *\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)$

$R^3=\left(\begin{array}{cccc}0& *& *& *\\ 0& 0& *& *\\ 0& 0& 0& *\\ 0& 0& 0& 0\end{array}\right)\left(\begin{array}{cccc}0& 0& *& *\\ 0& 0& 0& *\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)=\left(\begin{array}{cccc}0& 0& 0& *\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)$

$R^4=\left(\begin{array}{cccc}0& *& *& *\\ 0& 0& *& *\\ 0& 0& 0& *\\ 0& 0& 0& 0\end{array}\right)\left(\begin{array}{cccc}0& 0& 0& *\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)$