Day 22:

More on diagonalizable matrices

The Power Method

Diagonalizable matrices

Theorem. If \(A\) is an \(n\times n\) matrix and \(\{v_{1},\ldots,v_{n}\}\) is a basis of eigenvectors of \(A\) with corresponding eigenvalues \(\{\lambda_{1},\ldots,\lambda_{n}\}\), then

\[A = X\Lambda X^{-1}\]

where

\[X = \begin{bmatrix} | & | & | & & |\\ v_{1} & v_{2} & v_{3} & \cdots & v_{n}\\  | & | & | & & |\end{bmatrix}\]

and

\[\Lambda = \text{diag}(\lambda_{1},\lambda_{2},\ldots,\lambda_{n}) = \begin{bmatrix} \lambda_{1} & 0 & \cdots & 0\\ 0 & \lambda_{2} & & \vdots\\ \vdots & & \ddots & 0\\ 0 & \cdots & 0 & \lambda_{n}\end{bmatrix}\]

Diagonalizable matrices

Proof.  We will show that \(Ax=X\Lambda X^{-1} x\) for an arbitrary \(x\in\R^{n}\). Since \(\{v_{1},\ldots, v_{n}\}\) is a basis, there are scalars \(\alpha_{1},\ldots,\alpha_{n}\) such that

\[x = \alpha_{1}v_{1} + \alpha_{2}v_{2} + \cdots + \alpha_{n}v_{n}\]

Multiplying by \(A\) on the left we have

\[Ax = \alpha_{1}Av_{1} + \alpha_{2}Av_{2} + \cdots + \alpha_{n}Av_{n}\]

\[ = \alpha_{1}\lambda_{1}v_{1} + \alpha_{2}\lambda_{2}v_{2} + \cdots + \alpha_{n}\lambda_{n}v_{n}\]

By the definition of \(X\) we see that \(Xe_{i} = v_{i}\), and hence \(X^{-1}v_{i} = e_{i}\) for \(i=1,\ldots,n\). Hence,

\[X\Lambda X^{-1}x = X\Lambda X^{-1}\big(\alpha_{1}v_{1} + \alpha_{2}v_{2} + \cdots + \alpha_{n}v_{n}\big)\]

\[ = X\Lambda \big(\alpha_{1}X^{-1}v_{1} + \alpha_{2}X^{-1}v_{2} + \cdots + \alpha_{n}X^{-1}v_{n}\big)\]

\[ = X\Lambda \big(\alpha_{1}e_{1} + \alpha_{2}e_{2} + \cdots + \alpha_{n}e_{n}\big)\]

\[ = X\big(\alpha_{1}\lambda_{1} e_{1} + \alpha_{2}\lambda_{2} e_{2} + \cdots + \alpha_{n}\lambda_{n}e_{n}\big)\]

\[ = \alpha_{1}\lambda_{1} v_{1} + \alpha_{2}\lambda_{2} v_{2} + \cdots + \alpha_{n}\lambda_{n}v_{n}.\]

\(\Box\)

Diagonalizable matrices

Example. Consider the matrix

\[A=\left[\begin{array}{rrr} 2 & -1 & 2\\ 0 & 3 & -2\\ 0 & 0 & -1\end{array}\right]\]

We want to find \(\lambda\in\R\) such that

\[A-\lambda I=\left[\begin{array}{rrr} 2-\lambda & -1 & 2\\ 0 & 3-\lambda & -2\\ 0 & 0 & -1-\lambda\end{array}\right]\]

is not invertible. We know that an uppertriangular matrix is invertible if and only if it has no zeros on the diagonal. Hence, the eigenvalues of \(A\) are \(2\), \(3\), and \(-1\).

Hence, we need to find each of the subspace \(N(A-\lambda I)\) for \(\lambda=2,3,\) and \(-1\).

Diagonalizable matrices

\[A-(-1) I=\left[\begin{array}{rrr} 3 & -1 & 2\\ 0 & 4 & -2\\ 0 & 0 & 0\end{array}\right]\]

Row reducing we get

\[\text{rref}(A-(-1)I) = \left[\begin{array}{rrr} 1 &0 & \frac{1}{2}\\ 0 & 1 & -\frac{1}{2}\\ 0 & 0 & 0\end{array}\right]\]

Hence, we see that

\[N(A-(-1)I) = \text{span}\left\{\begin{bmatrix} -\frac{1}{2}\\ \frac{1}{2}\\ 1\end{bmatrix}\right\}\]

Similarly

\[N(A-3I) = \text{span}\left\{\begin{bmatrix} -1\\ 1\\ 0\end{bmatrix}\right\}\quad\text{and  }\quad N(A-2I) = \text{span}\left\{\begin{bmatrix} 1\\ 0\\ 0\end{bmatrix}\right\}\]

Diagonalizable matrices

Selecting a nonzero vector from each eigenspace and forming the matrix with those columns:

\[X=\begin{bmatrix} 1 & -1 & -\frac{1}{2}\\ 0 & 1 & \frac{1}{2}\\ 0 & 0 & 1\end{bmatrix}\]

Since the columns of \(X\) form a basis for \(\R^{3}\), we see that \(X\) is invertible. Indeed,

\[X^{-1} = \begin{bmatrix} 1 & 1 & 0\\ 0 & 1 & -\frac{1}{2}\\ 0 & 0 & 1\end{bmatrix}.\]

Finally,

\[\begin{bmatrix} 1 & -1 & -\frac{1}{2}\\ 0 & 1 & \frac{1}{2}\\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix} 2 & 0 & 0\\ 0 & 3 & 0\\ 0 & 0 & -1\end{bmatrix}\begin{bmatrix} 1 & 1 & 0\\ 0 & 1 & -\frac{1}{2}\\ 0 & 0 & 1\end{bmatrix} = \left[\begin{array}{rrr} 2 & -1 & 2\\ 0 & 3 & -2\\ 0 & 0 & -1\end{array}\right]\]

Hence, \(A\) is diagonalizable!

\(X\)

\(\Lambda\)

\(X^{-1}\)

\(A\)

Diagonalizable matrices

Another Example. Consider the matrix

\[A = \left[\begin{array}{rrr} -1 &\ \ 0 & 0\\ -8 & 3 & 8\\ 0 & 0 & -1\end{array}\right].\]

The eigenvalues of \(A\) are \(3\) and \(-1\). Diagonalize \(A\), that is, find an invertible matrix \(X\) and a diagonal matrix \(\Lambda\) such that \(A=XAX^{-1}\). 

\[\text{rref}(A-3I) = \begin{bmatrix} 1 & 0 & 0\\ 0 & 0 & 1\\ 0 & 0 & 0\end{bmatrix}\]

Therefore

\[\left\{\begin{bmatrix} 0\\ 1\\ 0\end{bmatrix}\right\}\]

is a basis for the eigenspace of \(A\) associated with the eigevalue \(3\).

Diagonalizable matrices

Another Example. Next we note that 

\[\text{rref}(A-(-1)I) = \begin{bmatrix} 1 & -\frac{1}{2} & -1\\ 0 & 0 & 0\\ 0 & 0 & 0\end{bmatrix}.\]

From this we see that

\[\left\{\begin{bmatrix} 1/2\\ 1\\ 0\end{bmatrix},\begin{bmatrix} 1\\ 0\\ 1\end{bmatrix}\right\}\]

is a basis for the eigenspace of \(A\) associated with the eigenvalue \(-1\).

Diagonalizable matrices

Another Example. If we set

\[X = \left[\begin{array}{rrr} 0 & \frac{1}{2} & 1\\ 1 & 1 & 0\\ 0 & 0 & 1\end{array}\right]\]

\[\Lambda = \left[\begin{array}{rrr} 3 & 0 & 0\\ 0 & -1 & 0\\ 0 & 0 & -1\end{array}\right]\]

and

then

\[X^{-1} = \left[\begin{array}{rrr} -2 &\ \ 1 & 2\\ 2 & 0 & -2\\ 0 & 0 & 1\end{array}\right]\]

and

\[X\Lambda X^{-1} = \left[\begin{array}{rrr} 0 & \frac{1}{2} & 1\\ 1 & 1 & 0\\ 0 & 0 & 1\end{array}\right]\left[\begin{array}{rrr} 3 & 0 & 0\\ 0 & -1 & 0\\ 0 & 0 & -1\end{array}\right]\left[\begin{array}{rrr} -2 &\ \ 1 & 2\\ 2 & 0 & -2\\ 0 & 0 & 1\end{array}\right]\]

\[= \left[\begin{array}{rrr} -1 &\ \ 0 & 0\\ -8 & 3 & 8\\ 0 & 0 & -1\end{array}\right] = A\]

Finding Eigenvalues

Theorem. Let \(A\) be an \(n\times n\) matrix, and let \(v\) be a nonzero vector. Define the sequence \((v_{n})\) by \(v_{1} = \frac{1}{\|v\|}v\), and

\[v_{k+1} = \frac{Av_{k}}{\|Av_{k}\|}\quad\text{for all }k\in\N.\]

If \(A\) is diagonalizable with eigenvalues \(\lambda_{1},\ldots,\lambda_{n}\), and \(|\lambda_{1}|>|\lambda_{2}|\geq|\lambda_{3}|\geq\cdots\geq |\lambda_{n}|\), then the sequence

\[v_{1}^{\top}Av_{1}, v_{2}^{\top}Av_{2}, v_{3}^{\top}Av_{3},\ldots,v_{k}^{\top}Av_{k},\ldots\]

converges to an eigenvalue of \(A\) for almost all initial vectors \(v\).

Example. Let

\[A = \left[\begin{array}{rrrr} -1 & -3 & 7 & 5\\ -3 & -1 & 5 & 7\\ 7 & 5 & -1 & -3\\ 5 & 7 & -3 & -1\end{array}\right].\]

Set \(v_{1} = \begin{bmatrix} 1 & 0 & 0 & 0\end{bmatrix}^{\top}\). Then,

\[v_{1} = \begin{bmatrix}1\\ 0\\ 0\\ 0\end{bmatrix},\ v_{2} = \begin{bmatrix}-0.1091\\ -0.3273\\ 0.7638\\ 0.5455\end{bmatrix}, v_{3} = \begin{bmatrix} 0.6355\\ 0.5750\\ -0.3329\\ -0.3934 \end{bmatrix}, v_{4} = \begin{bmatrix} -0.4263 \\ -0.4418\\ 0.5658\\ 0.5503 \end{bmatrix}, \]

\[v_{5} = \begin{bmatrix} 0.5322 \\ 0.5283 \\ -0.4659\\ -0.4698 \end{bmatrix},\ldots,v_{12} = \begin{bmatrix} -0.4998\\ -0.4998 \\ 0.5002\\ 0.5002 \end{bmatrix}, v_{13} = \begin{bmatrix} 0.5001\\ 0.5001\\ -0.4999\\ -0.4999 \end{bmatrix},\ldots\]

Example Continued. However,

\[v_{1}^{\top}Av_{1} = -1,\quad v_{2}^{\top}Av_{2} = -10.4762,\quad v_{3}^{\top}Av_{3} = -14.5201,\]

\[v_{4}^{\top}Av_{4}= -15.6261,\quad  v_{5}^{\top}Av_{5} = -15.9063, \quad \ldots\]

\[ v_{12}^{\top}Av_{12} = -15.999994277954329,\quad\ldots\quad,v_{40}^{\top}Av_{40} = -16\]

Why Does this Power Method Work?

If \(A\) is diagonaliable, then there is an invertible matrix \(X\) and a diagonal matrix \(\Lambda\) such that \(A = X\Lambda X^{-1}.\)

\[A^2 = (X\Lambda X^{-1})(X\Lambda X^{-1}) = X\Lambda\Lambda X^{-1} = X\Lambda^{2} X^{-1}\]

\[A^3 = AA^2 = (X\Lambda X^{-1})(X\Lambda^2 X^{-1}) = X\Lambda\Lambda^2 X^{-1} = X\Lambda^{3} X^{-1}\]

\[A^4 = AA^3 = (X\Lambda X^{-1})(X\Lambda^3 X^{-1}) = X\Lambda\Lambda^3 X^{-1} = X\Lambda^{4} X^{-1}\]

\[\vdots\]

In general,

\[A^{k} = X\Lambda^{k} X^{-1}.\]

Why Does this Power Method Work?

Since \(\Lambda\) is a diagonal matrix, we have

\[\Lambda = \begin{bmatrix} \lambda_{1} & 0 & \cdots & 0\\ 0 & \lambda_{2} & & \vdots\\ \vdots & & \ddots & \vdots\\ 0 & \cdots & \cdots & \lambda_{n}\end{bmatrix}\quad\text{and}\quad \Lambda^{k} = \begin{bmatrix} \lambda_{1}^{k} & 0 & \cdots & 0\\ 0 & \lambda_{2}^{k} & & \vdots\\ \vdots & & \ddots & \vdots\\ 0 & \cdots & \cdots & \lambda_{n}^{k}\end{bmatrix}\]

If we assume \(|\lambda_{1}|> |\lambda_{2}|\geq\cdots\geq |\lambda_{n}|\), then we have

\[\frac{A^{k}}{\lambda_{1}^{k}} = X\begin{bmatrix} 1 & 0 & \cdots & 0\\ 0 & (\lambda_{2}/\lambda_{1})^{k} & & \vdots\\ \vdots & & \ddots & \vdots\\ 0 & \cdots & \cdots & (\lambda_{n}/\lambda_{1})^{k}\end{bmatrix} X^{-1} \approx X\begin{bmatrix} 1 & 0 & \cdots & 0\\ 0 & 0 & & \vdots\\ \vdots & & \ddots & \vdots\\ 0 & \cdots & \cdots & 0\end{bmatrix} X^{-1}\]

Why Does this Power Method Work?

Note that

\[v_{k+1} = \frac{A^{k}v}{\|A^{k}v\|}\approx \pm \frac{x_{1}}{\|x_{1}\|}\]

Hence,

\[v_{k+1}^{\top}Av_{k+1} \approx \left( \pm \frac{x_{1}}{\|x_{1}\|}\right)^{\top}A\left( \pm \frac{x_{1}}{\|x_{1}\|}\right) = \frac{1}{\|x_{1}\|^{2}} x_{1}^{\top}\lambda_{1}x_{1} = \lambda_{1}\]

What can go wrong?

If \(\lambda\) and \(-\lambda\) are both eigenvalues:

Example. Consider the matrix \(A = \begin{bmatrix} 1 & 0\\ 0 & -1\end{bmatrix}.\)

If we start with a vector \(v_{1}=\begin{bmatrix} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\end{bmatrix}\), then

\[v_{2} = \frac{Av_{1}}{\|Av_{1}\|^{2}} = \begin{bmatrix} \frac{1}{\sqrt{2}}\\ -\frac{1}{\sqrt{2}}\end{bmatrix},\ v_{3} = \begin{bmatrix} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\end{bmatrix},\ v_{4} = \begin{bmatrix} \frac{1}{\sqrt{2}}\\ -\frac{1}{\sqrt{2}}\end{bmatrix},\ v_{5}  = \begin{bmatrix} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\end{bmatrix},\ldots\]

How can we fix it?

Example. Add a bit \(A+I = \begin{bmatrix} 2 & 0\\ 0 & 0\end{bmatrix}.\)

If we start with a vector \(v_{1}=\begin{bmatrix} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\end{bmatrix}\), then

\[v_{2} = \frac{Av_{1}}{\|Av_{1}\|^{2}} = \begin{bmatrix} 1\\ 0\end{bmatrix},\ v_{3} = \begin{bmatrix} 1\\ 0\end{bmatrix},\ v_{4} = \begin{bmatrix} 1\\ 0\end{bmatrix},\ v_{5}  = \begin{bmatrix} 1\\ 0\end{bmatrix},\ldots\]

And

\[v_{2}^{\top}Av_{2} = 2,\ v_{3}^{\top}Av_{3} = 2,\ v_{4}^{\top}Av_{4} = 2,\ldots\]

How can we fix it?

Proposition. If \(A\) is a matrix with eigenvalues \(\lambda_{1},\ldots,\lambda_{n}\) and associated eigenvectors \(v_{1},\ldots,v_{n}\), then \(A+cI\) has eigenvalues \(\lambda_{1}+c,\ldots,\lambda_{n}+c\) with the same associated eigenvectors.

How do we code it?

Input: Matrix \(A\)

1. Let \(v\) be a random vector: \(\texttt{v=randn(n,1);}\)
2. Normalize \(v\): \(\texttt{v=v/norm(v);}\)
3. Compute an approximate eigenvalue: \(\texttt{la = v'*A*v;}\)
4. Check whether \(\texttt{la}\) is an eigenvalue and \(v\) is the eigenvector: \(\texttt{if norm(la*v-A*v)<thresh, done; end;}\)
5. If not, update \(v\): \(\texttt{v = Av/norm(Av);}\) and go back to step \(3\).

Once this algorithm terminates with \(\texttt{done}\), then \(\texttt{la}\) is (approximately) an eigenvalue and \(\texttt{v}\) is a normalized eigenvector.

Notes: If the algorithm takes too long, you should start over with a new random vector.

If \(A\) has two opposite eigenvalues, then the algorithm may not work. So, we start by adding \(\texttt{c*eye(n)}\) to our matrix \(A\) where \(c\) is a random number, then subtract \(c\) from the eigenvalue we get at the end.