Let
\(A\) be a symmetric matrix. According to
Real Spectral Theorem, there exists an orthogonal matrix
\(P\) such that
\(P^TAP\) is diagonal, with the entries on the diagonal given by the eigenvalues of
\(A\text{.}\)
In this worksheet, we will work through the process of finding the eigenvalues of a symmetric matrix, along with the orthogonal matrix \(P\text{.}\)
We will also look at the method described in
Section 4.7 for finding dominant eigenvalues.
Subsection 4.3.1 Orthogonal diagonalization
Consider the following \(5\times 5\) matrix \(A\text{:}\)
\begin{equation*}
A = \bbm 97\amp -23 \amp -23\amp -3\amp -46\\
-23\amp 49\amp 1\amp 21\amp 2\\
-23\amp 1\amp 49\amp 21\amp 2\\
-3\amp 21\amp 21\amp 9\amp 42\\
-46\amp 2\amp 2\amp 42\amp 52\ebm\text{.}
\end{equation*}
Exercise 4.3.1.
(a)
First, enter the matrix and confirm that it’s symmetric. You can compute the transpose of a SymPy matrix A using A.T.
(b)
Now, find the eigenvalues of
\(A\text{.}\) Yes, we could jump straight to the answer, but let’s go through the steps anyway. First, compute the characteristic polynomial of
\(A\text{.}\) You may with to refer to
Subsection B.3.2 for the correct syntax to use.
Now, factor the polynomial:
Finally, based on the output, input the eigenvalues into the cell below.
(c)
Recall that if \(\vv\) is an eigenvector of \(A\) corresponding to the eigenvalue \(\lambda\text{,}\) then \(\vv\) belongs to the nullspace of \(\lambda I-A\text{,}\) where \(I\) is the identity matrix. You can create an \(n\times n\) identity matrix in SymPy using the syntax eye(n).
For each eigenvalue \(\lambda\) found in part (b), compute the nullspace of \(\lambda I-A\text{.}\) You will want to refer to these nullspaces later, so give each one a name. For example, if your first eigenvalue was 7, you could enter something like
E1 = (7*sy.eye(5)-A).nullspace()
E1
to get the first eigenspace. Three code cells are provided below. If you are in Jupyter, you can add more cells by clicking the \(+\) button in the toolbar.
Finally, let’s check our work. Use the command A.eigenvects() to compute the eigenvalues and eigenvectors in one step, and confirm that the results match what you found above.
Next, we want to diagonalize \(A\text{.}\) There is a diagonalize command, but it won’t do orthogonal diagonalization, which is what we want.
Recall that the eigenvectors corresponding to distinct eigenvalues are automatically orthogonal. But this is not true of repeated eigenvalues, and there’s a good chance you’ve found that one of your eigenvalues has multiplicity greater than 1.
A matrix \(P\) is orthogonal if \(P^TP=I\text{,}\) in which case te columns of \(P\) form an orthonormal basis. (So we are looking for an orthonormal basis of eigenvectors.)
For any 1-dimensional eigenspaces, you will just need to find a unit eigenvector. For example, if \(\vv = (2,1,0,3,1)\) is an eigenvector for an eigenvalue of multiplicity 1, then you will want to use the eigenvector \(\uu = \frac{1}{\sqrt{15}}(2,1,0,3,1)\text{,}\) since \(\len{\vv} = \sqrt{15}\text{.}\)
For any higher-dimensional eigenspaces, you can use the Gram-Schmidt algorithm. Recall that the optional argument true will produce unit vectors: given a basis B for an eigenspace, the command sy.GramSchmidt(B,true) will produce an orthonormal basis.
Exercise 4.3.2.
(a)
Determine an orthonormal basis of eigenvectors of \(A\text{.}\)
(b)
Once you have found five orthonormal eigenvectors create a matrix \(P\) whose columns are those eigenvectors:
(c)
Confirm that your matrix is orthogonal by computing \(P^TP\text{.}\) You should get the identity matrix at the result.
(d)
Finally, compute the matrix \(P^TAP\text{,}\) and confirm that the result is a diagonal matrix whose entries are the eigenvalues of \(A\text{.}\)
Subsection 4.3.2 The power method for dominant eigenvalues
For the problem you just completed, you may have noticed that there was one eigenvalue that was larger in absolute value than all the other eigenvalues. Such an eigenvalue is called a dominant eigenvalue.
If you complete
Section 4.6, you will see that the state of such systems is ultimately given by a sum of certain geometric series (involving eigenvalues of a matrix), and that the long-term behaviour of the system is approximately geometric, and governed by the dominant eigenvalue.
The power method states that, given an initial state vector \(\xx_0\text{,}\) the sequence of vectors
\begin{equation*}
\xx_0, A\xx_0, A^2\xx_0,\ldots, A^k\xx_0,\ldots
\end{equation*}
will converge to an eigenvector for the dominant eigenvalue.
Exercise 4.3.3.
Let’s see if this is true. Below, you are given an initial guess
x0 and an empty list
L. You want to populate
L with vectors of the form
\(A^k\xx_0\text{,}\) for
\(0\leq k\leq 10\text{.}\) This is most easily done using a
for loop. For details on syntax, see
Subsubsection 4.7.2.1
(a)
Let’s print the last vector in the list, which will probably have some pretty big entries:
(b)
Now let’s check our work. How can we tell if this vector is (approximately) a multiple of the dominant eigenvector? One option is to divide each entry in L[10] by the smallest (in absolute value) non-zero entry in L[10]. How does this compare to the eigenvector you originally found for the dominant eigenvalue?
If you find that some entries look right, but not others, see the last part of the worksheet.
(c)
Finally, suppose we didn’t know what the dominant eigenvalue was, and we wanted to find it. Note that if \(\xx\) is a dominant eigenvector, then \(A\xx = \lambda\xx\text{,}\) where \(\lambda\) is the dominant eigenvalue. Then
\begin{equation*}
\lambda = \frac{\xx\cdot (\lambda\xx)}{\xx\cdot \xx} = \frac{\xx\cdot A\xx}{\xx\cdot \xx}\text{.}
\end{equation*}
It seems reasonable, then, that if \(\xx_k\) is approximately equal to a dominant eigenvector, then the numbers
\begin{equation*}
r_k = \frac{\xx_k\cdot A\xx_k}{\xx_k\cdot \xx_k}=\frac{\xx_k\cdot \xx_{k+1}}{\lVert \xx_k\rVert^2}
\end{equation*}
should be approximately equal to the dominant eigenvalue. (These numbers are the so-called Rayliegh quotients.)
Following the approach in
Subsubsection 4.7.2.1, compute the Rayleigh quotients
\(r_k\) for
\(1\leq k\leq 10\text{,}\) and comment on how well they approximate the dominant eigenvalue of
\(A\text{.}\)
(d)
Does it look like the numbers are approaching the dominant eigenvalue? Or do they seem to be getting near one of the other eigenvalues?
If your numbers seem wrong, it might be for the following reason: when we write our initial guess x0 as a linear combination of the eigenvectors, the coefficient of the dominant eigenvector has to be nonzero. Is that the case here? To check, note that if \(\vec{c}\) is the vector of coefficients, then we must have \(P\vec{c}=\xx_0\text{.}\) Can you solve this equation for \(\vec{c}\text{?}\)
If the first entry in c is zero, enter a new value for \(x0\text{,}\) and try the above steps again in the cell below.
