You are about to erase your work on this activity. Are you sure you want to do this?
Updated Version Available
There is an updated version of this activity. If you update to the most recent version of this activity, then your current progress on this activity will be erased. Regardless, your record of completion will remain. How would you like to proceed?
Mathematical Expression Editor
The Power Method and the Dominant Eigenvalue
We know eigenvalues as the roots of a characteristic polynomial of a matrix. In
practice, the problem of finding eigenvalues is virtually never solved by finding roots
of the characteristic polynomial, as this task is computationally prohibitive for large
matrices. Iterative methods are used instead. In this section we describe the
Power Method and some variations. In QR Factorization we will touch upon
a better iterative method for approximating eigenvalues of a non-singular
matrix called the QR-algorithm. The methods we present in this text are a
mere introduction to the iterative methods used to compute eigenvalues and
eigenvectors.
Finding the Dominant Eigenvalue using the Power Method
In Exploration ?? of Diagonalizable Matrices and Multiplicity, our initial rationale
for diagonalizing matrices was to be able to compute the powers of a square matrix,
and the eigenvalues were needed to do this. In this section, we are interested in
efficiently computing eigenvalues, and it may come as no surprise that the first
method we discuss uses the powers of a matrix.
An eigenvalue \(\lambda \) of an \(n \times n\) matrix \(A\) is called a dominant eigenvalue if \(\lambda \) has multiplicity \(1\),
and
\begin{equation*} |\lambda | > |\mu | \quad \mbox { for all eigenvalues } \mu \neq \lambda \end{equation*}
Any corresponding eigenvector is called a dominant eigenvector of \(A\).
When such an eigenvalue exists, one technique for finding it is as follows: Let \(\vec {x}_{0}\) in \(\RR ^n\) be
a first approximation to a dominant eigenvector \(\lambda \), and compute successive
approximations \(\vec {x}_{1}, \vec {x}_{2}, \dots \) as follows:
\begin{equation*} \vec {x}_{k+1} = A\vec {x}_{k} \quad \mbox { for each } k \geq 0 \end{equation*}
We will see in this section that if
the first estimate \(\vec {x}_{0}\) is good enough, these vectors \(\vec {x}_{n}\) will approximate the dominant
eigenvector \(\lambda \). This technique is called the Power Method (because \(\vec {x}_{k} = A^{k}\vec {x}_{0}\) for each \(k \geq 1\)). Observe
that if \(\vec {z}\) is any eigenvector corresponding to \(\lambda \), then
Because the vectors \(\vec {x}_{1}, \vec {x}_{2}, \dots , \vec {x}_{n}, \dots \)
approximate dominant eigenvectors, this suggests that we define the Rayleigh
quotients as follows:
Then the numbers \(r_{k}\) approximate the dominant eigenvalue
\(\lambda \).
Use the power method to approximate a dominant eigenvector and eigenvalue of
\(A = \left [ \begin{array}{rr} 1 & 1 \\ 2 & 0 \end{array}\right ]\).
The eigenvalues of \(A\) are \(2\) and \(-1\), with eigenvectors \(\left [ \begin{array}{rr} 1 \\ 1 \end{array}\right ]\) and \(\left [ \begin{array}{rr} 1 \\ -2 \end{array}\right ]\), as we can verify using
Octave.
For more information on using Octave to find eigenvalues, please refer to the Octave
Tutorial section of this text.
If we take \(\vec {x}_{0} = \left [ \begin{array}{rr} 1 \\ 0 \end{array}\right ]\) as the first approximation and compute \(\vec {x}_{1}, \vec {x}_{2}, \dots ,\) successively, from \(\vec {x}_{1} = A\vec {x}_{0}, \vec {x}_{2} = A\vec {x}_{1}, \dots \), the result is
These vectors are approaching scalar multiples of the dominant eigenvector \(\left [ \begin{array}{rr} 1 \\ 1 \end{array}\right ]\).
Moreover, we will find that the Rayleigh quotients are
and these are approaching the
dominant eigenvalue \(2\).
(a)
Verify the above calculations using the following Octave code.
(b)
Try a different starting vector. Do we still get convergence?
(c)
What happens if the starting vector is changed to \(\vec {x}_{0} = \left [ \begin{array}{rr} 1 \\ -2 \end{array}\right ]\)? Why doesn’t this
method work in this case?
(d)
Try using \(\vec {x}_{0} = \left [ \begin{array}{rr} 1 \\ -2.01 \end{array}\right ]\). You may need to increase the number of iterations to see
convergence.
% The Power Method
A=[1 1; 2 0]; %Input a square matrix
x_0=[1;0]; %Input an initial vector
maxit =5; %Choose number iterations
x = x_0;
count = 0;
r=0;
while count < maxit
y=x;
x = A*y;
if count > 0;
r=(transpose(y)*x)/(transpose(y)*y);
RQUOTIENTS(count)=r;
end
count = count + 1;
VECTORS(:,count)=x;
end
VECTORS
RQUOTIENTS
To see why the power method works, let \(\lambda _{1}, \lambda _{2}, \dots , \lambda _{m}\) be eigenvalues of \(A\) with \(\lambda _{1}\) dominant and let \(\vec {y}_{1}, \vec {y}_{2}, \dots , \vec {y}_{m}\)
be corresponding eigenvectors. What is required is that the first approximation \(\vec {x}_{0}\)
be a linear combination of these eigenvectors:
The right side approaches \(a_{1}\vec {y}_{1}\) as \(k\) increases because \(\lambda _{1}\) is
dominant
\(\left ( \left |\frac {\lambda _{i}}{\lambda _{1}} \right | < 1 \mbox { for each } i > 1 \right )\). Because \(a_{1} \neq 0\), this means that \(\vec {x}_{k}\) approximates the dominant eigenvector \(a_{1}\lambda _{1}^k\vec {y}_{1}\).
The power method requires that the first approximation \(\vec {x}_{0}\) be a linear combination of
eigenvectors. In Example , the eigenvectors form a basis of \(\RR ^2\). But even in this case,
the method fails if \(a_{1} = 0\), where \(a_{1}\) is the coefficient of the dominant eigenvector (as we saw
in part (c)). In general, the rate of convergence is quite slow if any of the ratios \(\left | \frac {\lambda _{i}}{\lambda _{1}} \right |\) is
near \(1\) (as we saw in part (d)).
Use the Power Method to approximate a dominant eigenvector and eigenvalue of
\(A = \left [ \begin{array}{rrr} 9 & 2 & 8 \\ 2 & -6 & -2 \\ -8 & 2 & -5 \end{array}\right ]\)
You can use Octave to employ the power method.
% The Power Method
A=[9 2 8; 2 -6 -2; -8 2 -5]; %Input a square matrix
xinit=[1;1;1];
%Input an initial vector
maxit =15;
%Choose number iterations
x = xinit;
count = 0;
r=0;
while count < maxit
y=x;
x = A*y;
if count > 0;
r=(transpose(y)*x)/(transpose(y)*y);
RQUOTIENTS(count)=r;
end
count = count + 1;
VECTORS(:,count)=x;
end
VECTORS
RQUOTIENTS
[EIGENVECTORS, EIGENVALUES]=eig(A)
VECTORS(:,maxit)/norm(VECTORS(:,maxit))
We observe that the Rayleigh quotients approach \(-3\), which is the dominant eigenvalue.
We can confirm this using the command ‘eig’. To see that our vectors are approaching
the dominant eigenvector, in the last step we divide our final approximate vector
by its norm to obtain a unit vector. Note that we get the same result as
the ‘eig’ command. In practice, we achieve greater accuracy if we normalize
the approximate eigenvector at each iteration rather than waiting until the
end.
Finding Other Eigenvalues
If the Power Method converges to the dominant eigenvalue, how can we determine
other eigenvalues of a matrix? This section is devoted to some variations of the Power
Method that address this problem. Many of the ideas in this section can be employed
in other settings.
The Inverse Power Method
The Inverse Power Method is based upon the following theorem.
Let \(M\) be an invertible \(n \times n\) matrix. Then the eigenvalues of \(M^{-1}\) are reciprocals of the
eigenvalues of \(M\).
The proof is left as an exercise. See Practice Problem .
We return to the matrix in Exploration . Let \(M=\left [ \begin{array}{rrr} 9 & 2 & 8 \\ 2 & -6 & -2 \\ -8 & 2 & -5 \end{array}\right ]\) and let \(A=M^{-1}\).
Now use Octave to verify that \(A=\left [ \begin{array}{rrr} 17/3 & 13/3 & 22/3 \\ 13/3 & 19/6 & 17/3 \\ -22/3 & -17/3 & -29/3 \end{array}\right ]\) and that the eigenvalues of \(A\) are indeed the reciprocals
of the eigenvalues of \(M\).
% Eigenvalues of the inverse are reciprocals of eigenvalues of the matrix
M=[9 2 8; 2 -6 -2; -8 2 -5];
eig(M)
A=inv(M)
eig(A)
If we now apply the Power Method to \(A\), the Rayleigh quotients converge toward the
dominant eigenvalue of \(A\), which is \(-1\). If we take the reciprocal of the dominant
eigenvalue of \(A\), we will have found the “least dominant” eigenvalue of \(M\), i.e., the
eigenvalue closest to zero. Sure enough, the reciprocal of \(-1\) is \(-1\), which is the eigenvalue
of \(M\) closest to zero.
The Shifted Power Method
Again, we begin with a theorem.
Let \(M\) be an invertible \(n \times n\) matrix and let \(z\) be any number (real or complex). If \(\lambda \) is an
eigenvalue of \(M\), then \(\lambda - z\) is an eigenvalue of \(M - cI\).
The proof is left as an exercise. See Practice Problem .
Once again, let \(M=\left [ \begin{array}{rrr} 9 & 2 & 8 \\ 2 & -6 & -2 \\ -8 & 2 & -5 \end{array}\right ]\) and let \(A=M - (-3)I\).
Using Octave, you can verify that \(A=\left [ \begin{array}{rrr} 12 & 2 & 8 \\ 2 & -3 & -2 \\ -8 & 2 & -2 \end{array}\right ]\) and that the eigenvalues of \(A\) are obtained by
subtracting \(c=-3\) from (i.e. by adding 3 to) each eigenvalue of \(M\). The dominant eigenvalue
of \(M\) gives us an eigenvalue of zero for \(A\) (you may find it gives a number very close
to zero, due to rounding error). So if we apply the Power Method to \(A\) to
find its dominant eigenvalue (which is 5), we can add \(-3\) to get the “second
largest” eigenvalue of \(M\) (i.e., the eigenvalue second-farthest from zero), which is
2.
To see this, we will use Octave.
% The shifted power method
M=[9 2 8;
2 -6 -2;
-8 2 -5];
eig(M)
c=-3;
A=M-c*eye(3)
eig(A)
xinit=[1;1;1]; %Input an initial vector
maxit =15; %Choose number iterations
x = xinit;
count = 0;
r=0;
while count < maxit
y=x;
x = A*y;
if count > 0;
r=(transpose(y)*x)/(transpose(y)*y);
RQUOTIENTS(count)=r;
end
count = count + 1;
VECTORS(:,count)=x;
end
VECTORS
RQUOTIENTS
RQUOTIENTS(maxit-1)+c
The Shifted Inverse Power method combines the two ideas we just studied. Using the
Inverse Power method, we know how to find the eigenvalue closest to zero. We also
know that we can change which eigenvalue is the dominant eigenvalue using the
Shifted Power method. With the Shifted Inverse Power method, we are able to find
the eigenvalue closest to a target \(\tau \), where \(\tau \) can be any number, real or complex. We
begin with the following theorem.
Let \(M\) be an \(n \times n\) matrix and let \(\lambda \) be an eigenvalue of \(M\). If \(\tau \) is any number (real or complex)
other than an eigenvalue of \(M\), then \(\frac {1}{\lambda - \tau }\) is an eigenvalue of \((M - \tau I)^{-1}\).
Suppose \(M\vec {x}=\lambda \vec {x}\). By Theorem 4 we have
Since \(\tau \) is not an eigenvalue, \(\det (M - \tau I) \ne 0\), so \(M - \tau I\) is invertible. Also, since \(\tau \) is not an eigenvalue, \(\lambda \ne \tau \). We
multiply both sides on the left by \((M - \tau I)^{-1}\) and divide both sides by the nonzero scalar \(\lambda - \tau \). We
arrive at
Again we let \(M=\left [ \begin{array}{rrr} 9 & 2 & 8 \\ 2 & -6 & -2 \\ -8 & 2 & -5 \end{array}\right ]\). We will use the Shifted Inverse Power Method to compute the
eigenvalue of \(M\) closest to the target \(\tau = 1\).
Using Octave, you can verify that \(A=\left [ \begin{array}{rrr} 5.75 & 3.5 & 6.5 \\ 3.5 & 2 & 4 \\ -6.5 & -4 & -7.5 \end{array}\right ]\). When we apply the Power Method to \(A\) to find its
dominant eigenvalue, our Rayleigh quotients approach 1. This means that the
eigenvalue of \(M\) closest to our target \(\tau = 1\) is \(1+\tau \), or 2.
We will use Octave.
% The shifted inverse power method
M=[9 2 8; 2 -6 -2; -8 2 -5];
eig(M)
tau=1;
A=inv(M-tau*eye(3)) %I love the fact that eye(n) is the command for the identity matrix
eig(A)
xinit=[1;1;1]; %Input an initial vector
maxit =15; %Choose number iterations
x = xinit;
count = 0;
r=0;
while count < maxit
y=x;
x = A*y;
if count > 0;
r=(transpose(y)*x)/(transpose(y)*y);
RQUOTIENTS(count)=r;
end
count = count + 1;
VECTORS(:,count)=x;
end
VECTORS
RQUOTIENTS
RQUOTIENTS(maxit-1)+tau
Let \(\lambda \) be an eigenvalue of \(M\). Write down what this means.
Then try to show \(M^{-1}\vec {x} = \frac {1}{\lambda }\vec {x}\) for some vector \(\vec {x}\).
If \(\vec {x}\) is the eigenvector of \(M\) corresponding to \(\lambda \), it will also be
the eigenvalue of \(M-cI\) corresponding to \(\lambda - c\).
In each case, find the exact eigenvalues and determine corresponding eigenvectors.
Then start with \(\vec {x}_{0} = \left [ \begin{array}{rr} 1 \\ 1 \end{array}\right ]\) and compute \(\vec {x}_{4}\) and \(r_{3}\) using the power method.