Numerical Methods for Eigenvalues

Section 2.9 Numerical Methods for Eigenvalues

As mentioned above, the eigenvalues and eigenvectors of an \(n \times n\) matrix where \(n \ge 4\) must be found numerically instead of by hand. The numerical methods that are used in practice depend on the geometric meaning of eigenvalues and eigenvectors which is equation (2.7.1). The essence of all these methods is captured in the Power method, which we now introduce.

🔗

Subsection 2.9.1 The Power Method

In the command window of MATLAB enter

>> A = hilb(5)
>> x = ones(5,1)
>> x = A*x
>> el = max(x)
>> x = x/el

Compare the new value of x with the original. Repeat the last three lines (you can use the scroll up button). Compare the newest value of x with the previous one and the original. Notice that there is less change between the second two. Repeat the last three commands over and over until the values stop changing. You have completed what is known as the Power Method. Now try the command

>> [v e] = eig(A)

The last entry in e should be the final el we computed. The last column in v is x/norm(x). Below we explain why our commands gave this eigenvalue and eigenvector.

🔗

For illustration consider a \(2 \times 2\) matrix whose eigenvalues are \(1/3\) and \(-2\) and whose corresponding eigenvectors are \(\vb_{1}\) and \(\vb_{2}\text{.}\) Let \(\xb_{0}\) be any vector which is a combination of \(\vb_{1}\) and \(\vb_{2}\text{,}\) e.g.,

\begin{equation*} \xb_{0} = \vb_{1} + \vb_{2}\text{.} \end{equation*}

Now let \(\xb_{1}\) be \(A\) times \(\xb_{0}\text{.}\) It follows from (2.7.1) that

\begin{equation*} \begin{split}\xb_{1}&= A \vb_{1} + A \vb_{2} \\&= \frac{1}{3}\vb_{1} + -2 \vb_{2}\end{split}\text{.} \end{equation*}

Thus the \(\vb_{1}\) part is shrunk while the \(\vb_{2}\) is stretched. If we repeat this process \(k\) times then

\begin{equation*} \begin{split}\xb_{k}&= A \xb_{k-1}\\&= A^{k} \xb_{0} \\&= \left(\frac{1}{3}\right)^{k} \vb_{1} + (-2)^{k} \vb_{2}\end{split}\text{.} \end{equation*}

Clearly, \(\xb_{k}\) grows in the direction of \(\vb_{2}\) and shrinks in the direction of \(\vb_{1}\text{.}\) This is the principle of the Power Method, vectors multiplied by \(A\) are stretched most in the direction of the eigenvector whose eigenvalue has the largest absolute value.

🔗

The eigenvalue with the largest absolute value is called the dominant eigenvalue. In many applications this quantity will necessarily be positive for physical reasons. When this is the case, the MATLAB code above will work since max(v) will be the element with the largest absolute value. In applications where the dominant eigenvalue may be negative, the program must use flow control to determine the correct number.

🔗

Summarizing the Power Method:

Repeatedly multiply \(\xb\) by \(A\) and divide by the element with the largest absolute value.
🔗

🔗
The element of largest absolute value converges to largest absolute eigenvalue.
🔗

🔗
The vector converges to the corresponding eigenvector.
🔗

🔗

Note that this logic only works when the eigenvalue largest in magnitude is real. If the matrix and starting vector are real then the power method can never give a result with an imaginary part. Eigenvalues with imaginary part mean the matrix has a rotational component, so the eigenvector would not settle down either.

🔗

Try

>> A = randn(15,15);
>> e = eig(A)

You can see that for a random square matrix, many of the eigenvalues are complex. However, matrices in applications are not just random. They have structure, and this can lead to real eigenvalues as seen in the next section.

🔗

Subsection 2.9.2 Symmetric, Positive-Definite Matrices

As noted in the previous paragraph, the power method can fail if \(A\) has complex eigenvalues. One class of matrices that appear often in applications and for which the eigenvalues are always real are called the symmetric matrices. A matrix is symmetric if

\begin{equation*} A' = A\text{,} \end{equation*}

i.e. \(A\) is symmetric with respect to reflections about its diagonal.

🔗

Try

>> A = rand(5,5)
>> C = A'*A
>> e = eig(C)

You can see that the eigenvalues of these symmetric matrices are real.

🔗

Next we consider an even more specialized class for which the eigenvalues are not only real, but positive. A symmetric matrix is called positive definite if for all vectors \(\vb \neq \mathbf{0}\) the following holds:

\begin{equation*} A \vb \cdot \vb > 0\text{.} \end{equation*}

Geometrically, \(A\) does not rotate any vector by more than \(\pi/2\text{.}\) In summary:

If \(A\) is symmetric then its eigenvalues are real.
🔗

🔗
If \(A\) is symmetric positive definite, then its eigenvalues are positive numbers.
🔗

🔗

🔗

Notice that the \(B\) matrices in the previous section were symmetric and the eigenvalues were all real. Notice that the Hilbert and Pascal matrices are symmetric.

🔗

Subsection 2.9.3 The residual of an approximate eigenvector-eigenvalue pair

If \(\vb\) and \(\lambda\) are an eigenvector-eigenvalue pair for \(A\text{,}\) then they are supposed to satisfy the equations: \(A \vb = \lambda \vb\text{.}\) Thus a scalar residual for approximate \(\vb\) and \(\lambda\) would be

\begin{equation*} r = \| A \vb - \lambda \vb \|\text{.} \end{equation*}

🔗

Subsection 2.9.4 The Inverse Power Method

In the application of vibration analysis, the mode (eigenvector) with the lowest frequency (eigenvalue) is the most dangerous for the machine or structure. The Power Method gives us instead the largest eigenvalue, which is the least important frequency. In this section we introduce a method, the Inverse Power Method which produces exactly what is needed.

🔗

The following facts are at the heart of the Inverse Power Method:

If \(\lambda\) is an eigenvalue of \(A\) then \(1/\lambda\) is an eigenvalue for \(A^{-1}\text{.}\)
🔗

🔗
The eigenvectors for \(A\) and \(A^{-1}\) are the same.
🔗

🔗

🔗

Thus if we apply the Power Method to \(A^{-1}\) we will obtain the largest absolute eigenvalue of \(A^{-1}\text{,}\) which is exactly the reciprocal of the smallest absolute eigenvalue of \(A\text{.}\) We will also obtain the corresponding eigenvector, which is an eigenvector for both \(A^{-1}\) and \(A\text{.}\) Recall that in the application of vibration mode analysis, the smallest eigenvalue and its eigenvector correspond exactly to the frequency and mode that we are most interested in, i.e. the one that can do the most damage.

🔗

Here as always, we do not really want to calculate the inverse of \(A\) directly if we can help it. Fortunately, multiplying \(\xb_{i}\) by \(A^{-1}\) to get \(\xb_{i+1}\) is equivalent to solving the system \(A \xb_{i+1}= \xb_{i}\,\text{,}\) which can be done efficiently and accurately. Since iterating this process involves solving a linear system with the same \(A\) but many different right hand sides, it is a perfect time to use the LU decomposition to save computations. The following function program does \(n\) steps of the Inverse Power Method.

function [x e] = myipm(A,n)
  % Performs the inverse power method.
  % Inputs: A -- a square matrix   
  %         n --  number of iterations.
  % Outputs: x -- estimated eigenvector
  %          e -- estimated smallest eigenvalue. 
  [L U P] = lu(A);  % LU decomposition of A with pivoting
  m = size(A,1); % determine the size of A
  x = ones(m,1); % make an initial vector with ones
  for i = 1:n
      px = P*x;  % Apply pivot
      y = L\px;  % solve via LU
      x = U\y;
      % find the maximum entry in absolute value, retaining its sign
      M = max(x);  
      m = min(x);
      if abs(M) >= abs(m)
        el = M;
      else
        el = m;
      end
      x = x/el; % divide by the estimated eigenvalue of the inverse of A
  end
  e = 1/el; % reciprocate to get an eigenvalue of A
end

Create a \(5 \times 5\) Hilbert matrix \(H\) and use the program to find its smallest eigenvalue and corresponding eigenvector. Also do [V, E] = eig(H) and compare.

🔗

Exercises 2.9.5 Exercises

1.

For each of the matrices

\begin{equation*} A =\left[\begin{array}{cc}1&2\\ 3&4\end{array}\right] \quad \text{and}\quad B =\left[\begin{array}{ccc}-2&1&0\\ 1&-2&1\\ 0&1&-3\end{array}\right]\text{,} \end{equation*}

perform two iterations of the power method by hand starting with a vector of all ones. State the resulting approximations of the eigenvalue and eigenvector.

🔗

2.

Write a well-commented MATLAB function program mypm.m that inputs a matrix and a tolerance, applies the power method until the scalar residual is less than the tolerance, and outputs the estimated eigenvalue and eigenvector, the number of steps, and the scalar residual.
🔗

🔗
Test your program on the matrices \(A\) and \(B\) in the previous exercise and check that you get the same results as your hand calculations for the first 2 iterations. Turn in your program only.
🔗

🔗

🔗

Prev Top Next