HW FOUR

1. 4.2

   (a) Eigenvalue $3$ with eigenvector $(1,1)$ and eigenvalue $-2$ with eigevector $(-1,4)$

   (b) Eigenvalue $5$ with eigenvector $(1,-1)$ and eigenvalue $-1$ with eigenvector $(7,-1)$.

   (c) The matrix is deffective with eigenvalue $1$ and eigenvector $(1,1)$.

   (d) The eigenvalues are $(7+ sqrt{13})/2$, $(7- sqrt{13})/2$ and $1$ with corresponding eigenvectors $(0,(-3+ sqrt{13})/2)$, $(0,(-3- sqrt{13})/2)$ and $(1,0,0)$.

   (e) Eigenvalues are $(29+\sqrt{65})/2$, $16$ and $(29-\sqrt{65})/2$ with corresponding eigenvectors $((-7+\sqrt{65})/4,1,0)$, $(0,0,1)$ and $((-7-\sqrt{65})/4,1,0)$.

   (f)Eigenvalues are $(5+\sqrt{5})/2$, $(5-\sqrt{5})/2$ and 0 with eigenvectors $(0, (-1+\sqrt{5})/2,1)$, $(0, (-1-\sqrt{5})/2,1)$ and $(1,0,0)$.

2. 4.4

   It is easier to work with $4\times 4$ or $6\times 6$ matrix first to gain intuition. For example, type in matlab

   In what follow, we assume that $n$ is even.

   To find eigenvalues, write $det(A - \lambda I) =0.$. To calculate this determinant, expand it in the first raw. The second $(n-1)\times(n-1)$ determinant should be expanded in the last raw. The result is

   $$(\lambda-a)^n - (\lambda-a)^{n-2}=0.$$
   The roots of this are $\lambda=a$, and $\lambda=a\pm 1$.

   To calculate eigenvectors, write $A x = \lambda x$ in components.

   The eigenvectors corresponding to $\lambda=a$ are

   $$\begin{bmatrix}0 \\ 1 \\ \dots \\ 0 \\ 0\\ 0 \end{bmatrix},$$
   $$\begin{bmatrix}0 \\ \\ 0 \\ 1 \\ \dots \\ 0 \\ 0\\ 0 \end{bmatrix},$$
   $$\dots$$
   $$\begin{bmatrix}0 \\ \\ 0 \\ \dots \\ 0 \\ 1\\ 0 \end{bmatrix}.$$

   The eigenvector corresponding to $\lambda = a+1$ is

   $$\begin{bmatrix}1 \\ 0 \\ 0 \\ \dots 1 \end{bmatrix}.$$

   Similarly, the eigenvector corresponding to $\lambda = a-1$ is

   $$\begin{bmatrix}1 \\ 0 \\ 0 \\ \dots \\ 0\\ -1 \end{bmatrix}.$$

3. 4.5

   (a) Let ${\bf A} = (a_{ij}), i = 1\dots n, j=1\dots m$. Then ${\bf A}^T = (a_{ji})$.

   The product of matrices ${\bf C}= (C_{ij})$ and $D=(d_{ij})$ is given by

   $${\bf C}\cdot {\bf D} = ( \sum\limits_{k=1}^n a_{i k} b_{k j}).$$
   Then
   $${\bf A} \cdot {\bf A}^T =(\sum\limits_{k=1}^n a_{i k} a_{j k}),$$
   which is $i\to j$ symmetric.

   (b) Denote the eigenvector $x$ to be $x=(x_i)$ and using Rayleigh quontient, we write

   $$\lambda=\frac{ x*T {\bf {A}} x}{x^T x} = \frac{\sum\limits{ i j k... ...} = \frac{\sum\limits_{k=1}^n (\sum_i( a_{i k} x_i)^2)} {\sum\limits_{i}x_i^2}.$$ (44)

   Since the numerator and denominator are the sum of squares, they both are positive, and thus $\lambda$ is also positive.

   (c) The eigenvalues of ${\bf A}+\mu{\bf {I}}$ are $\lambda_i + \mu$. We showed in (b) that $\lambda_i$ are positive, and since $\mu$ is assumed to be postive, the sum is positive. Therefore the eigenvalues of ${\bf A}+\mu{\bf {I}}$ are positive.

   (d) We now write the Rayleigh quontient in vector form, substituting ${\bf A} = {\bf B} {\bf {B}}^T$ so that

   $$0=\lambda=\frac{ x*T {\bf {A}} x}{x^T x} = \frac{ x*T {\bf {B}^T} {\bf {B}} x}{x^T x} .$$ (45)

   Therefore, since ${x^T x}>0$ for nonzero $x$, we obtain that either ${\bf {B}} x=0$ or $x^T {\bf {B}}^T =0$.

   (e) If columns of ${\bf B}$ are linearly independent, then there are no ${\bf x}$ such that ${bf B}{\bf x}=0$. Then due to (b) above all eigenvalues are positive, and thus the ${\bf A}$ are positive definite.

4. 4.19

   (a) The power will fail since it will not be able to “decide” whether to converge to eigenvector corresponding to eigenvalue $2$ or $-2$. The resulting vector will therefore converge to the linear combination of the two eigenvectors corresponding to eigenvalue $2$ and $-2$.

   (b) Consider the matrix $B = A+10 I$, and use power method and inverse power method to compute eigenvalues $2$ and $-2$.

   (c) Consider the eigenvalues $1$ and $2$. Their ratio is $1/2$ and adding any positive $s$ would decrease the ratio to $(1+s)/(2+s)$. Now consider the eigenvalues $-2$ and $2$. Their absolute ratio is $1$. Adding any positive $s$ will make the ratio to be $\vert s-2\vert/(s+2)$. Any $s$ smaller than $2$ would decrease the ratio from one to zero. So adding $s$ makes one ratio worse, and another better. The optimum is where both ratios are as small as possible, i.e. they are equal, that is

   $$(1+s)/(2+s) = (2-s)/(s+2),\ \ 0<s<2.$$
   The solution to this is $s=1/2$.

5. 4.27

   This is a cool problem. The code that implements power iterations, inverse power iterations and power iterations with a shift is given below.

   Note that for the million times million matrix one does not need to form the matrix $A$. Instead, it is enough to store the vector $X$ and perform operations on the $k$ and $k\pm 1$ elements.

   %First Question of this problem is to write power method for this matrix
   
   n=100;format long;
   MainDiagonal=linspace(1,n,n);
   MainDiagonal(n)=pi*n;
   SuperDiagonal=ones(n-1,1);
   A = diag(MainDiagonal,0)+diag(SuperDiagonal,1)+diag(SuperDiagonal,-1);
   
   x =rand(n,1);
   
   yold=ones(n,1);
   ynew=10*yold;
   Counter=1;
   
   while norm(yold-ynew)>1e-100
   yold=ynew;Counter=Counter+1;
   ynew=A*yold;
   lambda=ynew'*A*ynew/(ynew'*ynew);
   %fprintf("Counter=%d Lambda= %f old = %f new=%f \n",Counter,lambda,norm(yold),norm(ynew))
    
   ynew=ynew/norm(ynew,2);
   end
   
   fprintf("Largest Eigenvalue is given by = %f\n",lambda);
   
   %Second Question of this problem is to write INVERSE power method for this matrix
   
   
   ynew=10*yold;
   
   while norm(yold-ynew)>1e-100
   yold=ynew;Counter=Counter+1;
   ynew=A\yold;
   lambda=ynew'*A*ynew/(ynew'*ynew);
   %fprintf("Counter=%d Lambda= %f old = %f new=%f \n",Counter,lambda,norm(yold),norm(ynew))
    
   ynew=ynew/norm(ynew,2);
   end
   
   fprintf("Smallest Eigenvalue is given by = %f\n",lambda);
   
   
   yold=ones(n,1);
   ynew=10*yold;
   B = A -100*eye(n,n);
   Counter=1;
   while Counter<10
   yold=ynew;Counter=Counter+1;
   ynew=B\yold;
   lambda=ynew'*B*ynew/(ynew'*ynew);
   %fprintf("Counter=%d Lambda= %f old = %f new=%f \n",Counter,lambda,norm(yold),norm(ynew))
   ynew=ynew/norm(ynew,2);
   end
   
   fprintf("Eigenvalue closest to 100 is given by = %f\n",lambda+100);
   
   
   %Fourth question is to consider 10^6 times 10^6 matrix
   
   n = 10^6;
   y0 = rand(10^6,1);
   
   yold=ones(n,1);
   ynew=10*yold;
   Counter=1;
   
   while norm(yold-ynew)>1e-100
   yold=ynew;Counter=Counter+1;
   
   ynew(1) = yold(1)+yold(2);
   ynew(n) = yold(n-1)+pi*n*yold(n);
   for k = 2:n-1
   	  ynew(k) = k*yold(k) + yold(k-1)+yold(k+1);
   end
   
   lambda=yold'*ynew/(yold'*yold);
   %fprintf("Counter=%d Lambda= %f old = %f new=%f \n",Counter,lambda,norm(yold),norm(ynew))
   ynew=ynew/norm(ynew,2);
   end
   
   fprintf("Million Over Million Largest Eigenvalue is given by = %f\n",lambda);
   

6. 4.23

   (a)

   $$\pm 10 \sqrt{10405.0} \simeq \pm 1020.049018429996823846313$$
   $$510+100\sqrt{26.0}\simeq 1019.9019513592784830028224,$$
   $$1020=1020.$$
   So there are four eigenvalues close by absolute value to $1020$.

   (b) Power method will fail since four eigenvalues are very close to each other, and $(1020/ 1020.05)^2 \simeq 0.999902$, which is very close to 1

   (c) If you do a shift, and $\omega=-1$, then the dominant eigenvalue is

   $$10 \sqrt{10405.0} +1 \simeq \pm 1021.05.$$

   (d) The next biggest eigenvalue in the shifted matrix will be $1021$, so the error will be reduced each time by

   $$(1021/1021.05)^2 \simeq 0.999902.$$

   Now solve

   $$0.999902^n = \frac{1}{10},$$
   which gives $n=23494$.

   (e) To converge to $10 \sqrt{10405.0}$ one can choose $\omega=-1$. Then the eigenvalue $10 \sqrt{10405.0}+1$ will be the dominant one, so power method will converge to it.