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 exp:motivate_diagonalization 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 of an matrix is called a dominant eigenvalue if has multiplicity ,
and \begin{equation*} |\lambda | > |\mu | \quad \mbox{ for all eigenvalues } \mu \neq \lambda \end{equation*}
Any corresponding eigenvector is called a dominant eigenvector of .
When such an eigenvalue exists, one technique for finding it is as follows: Let in be
a first approximation to a dominant eigenvector , and compute successive
approximations as follows: \begin{equation*} \vec{x}_{1} = A\vec{x}_{0} \quad \vec{x}_{2} = A\vec{x}_{1} \quad \vec{x}_{3} = A\vec{x}_{2} \quad \cdots \end{equation*}
In general, we define \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 is good enough, these vectors
will approximate the dominant eigenvector . This technique is called the Power
Method (because for each ). Observe that if is any eigenvector corresponding to ,
then \begin{equation*} \frac{\vec{z} \dotp (A\vec{z})}{\norm{ \vec{z} }^2} = \frac{\vec{z} \dotp (\lambda \vec{z})}{\norm{ \vec{z} }^2} = \lambda \end{equation*}
Because the vectors approximate dominant eigenvectors, this suggests that we define
the Rayleigh quotients as follows: \begin{equation*} r_{k} = \frac{\vec{x}_{k} \dotp \vec{x}_{k+1}}{\norm{ \vec{x}_{k} }^2} \quad \mbox{ for } k \geq 1 \end{equation*}
Then the numbers approximate the dominant eigenvalue .
Use the power method to approximate a dominant eigenvector and eigenvalue of
.
The eigenvalues of are and , with eigenvectors and , 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 as the first approximation and compute successively, from , the result is
\begin{equation*} \vec{x}_{1} = \left [ \begin{array}{rr} 1 \\ 2 \end{array}\right ], \
\vec{x}_{2} = \left [ \begin{array}{rr} 3 \\ 2 \end{array}\right ], \
\vec{x}_{3} = \left [ \begin{array}{rr} 5 \\ 6 \end{array}\right ], \
\vec{x}_{4} = \left [ \begin{array}{rr} 11 \\ 10 \end{array}\right ], \
\vec{x}_{3} = \left [ \begin{array}{rr} 21 \\ 22 \end{array}\right ], \ \dots \end{equation*}
These vectors are approaching scalar multiples of the dominant eigenvector .
Moreover, we will find that the Rayleigh quotients are \begin{equation*} r_{1} = \frac{7}{5}, r_{2} = \frac{27}{13}, r_{3} = \frac{115}{61}, r_{4} = \frac{451}{221}, \dots \end{equation*}
and these are approaching the dominant eigenvalue .
(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 ? Why doesn’t this
method work in this case?
(d)
Try using . 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 be eigenvalues of with dominant and let
be corresponding eigenvectors. What is required is that the first approximation be a
linear combination of these eigenvectors: \begin{equation*} \vec{x}_{0} = a_{1}\vec{y}_{1} + a_{2}\vec{y}_{2} + \dots + a_{m}\vec{y}_{m} \quad \mbox{ with } a_{1} \neq 0 \end{equation*}
If , the fact that and for each gives \begin{equation*} \vec{x}_{k} = a_{1}\lambda _{1}^k\vec{y}_{1} + a_{2}\lambda _{2}^k\vec{y}_{2} + \dots + a_{m}\lambda _{m}^k\vec{y}_{m} \quad \mbox{ for } k \geq 1 \end{equation*}
Hence \begin{equation*} \frac{1}{\lambda _{1}^k}\vec{x}_{k} = a_{1}\vec{y}_{1} + a_{2}\left (\frac{\lambda _{2}}{\lambda _{1}}\right )^k\vec{y}_{2} + \dots + a_{m}\left (\frac{\lambda _{m}}{\lambda _{1}}\right )^k\vec{y}_{m} \end{equation*}
The right side approaches as increases because is dominant
. Because , this means that approximates the dominant eigenvector .
The power method requires that the first approximation be a linear combination of
eigenvectors. In Example exp:2x2PowerMethod, the eigenvectors form a basis of . But even in this case,
the method fails if , where is the coefficient of the dominant eigenvector (as we saw
in part exp:2x2PowerMethod_c). In general, the rate of convergence is quite slow if any of the ratios is near
(as we saw in part exp:2x2PowerMethod_d).
Use the Power Method to approximate a dominant eigenvector and eigenvalue of
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 , 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 be an invertible matrix. Then the eigenvalues of are reciprocals of the
eigenvalues of .
If we now apply the Power Method to , the Rayleigh quotients converge toward the
dominant eigenvalue of , which is . If we take the reciprocal of the dominant
eigenvalue of , we will have found the “least dominant” eigenvalue of , i.e., the
eigenvalue closest to zero. Sure enough, the reciprocal of is , which is the eigenvalue
of closest to zero.
The Shifted Power Method
Again, we begin with a theorem.
Let be an invertible matrix and let be any number (real or complex). If is an
eigenvalue of , then is an eigenvalue of .
Using Octave, you can verify that and that the eigenvalues of are obtained by
subtracting from (i.e. by adding 3 to) each eigenvalue of . The dominant eigenvalue
of gives us an eigenvalue of zero for (you may find it gives a number very close
to zero, due to rounding error). So if we apply the Power Method to to
find its dominant eigenvalue (which is 5), we can add to get the “second
largest” eigenvalue of (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 , where can be any number, real or complex. We
begin with the following theorem.
Let be an matrix and let be an eigenvalue of . If is any number (real or complex)
other than an eigenvalue of , then is an eigenvalue of .
Proof
Suppose . By Theorem th:eig_shifted we have
Since is not an eigenvalue, , so is invertible. Also, since is not an eigenvalue,
. We multiply both sides on the left by and divide both sides by the nonzero
scalar . We arrive at
as desired.
Again we let . We will use the Shifted Inverse Power Method to compute the
eigenvalue of closest to the target .
Using Octave, you can verify that . When we apply the Power Method to to find its
dominant eigenvalue, our Rayleigh quotients approach 1. This means that the
eigenvalue of closest to our target is , 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