martes, 16 de mayo de 2023

The Power Method eigenvalyues Python

 

https://pythonnumericalmethods.berkeley.edu/notebooks/chapter15.02-The-Power-Method.html

The Power Method

Find the largest eigenvalue

In some problems, we only need to find the largest dominant eigenvalue and its corresponding eigenvector. In this case, we can use the power method - a iterative method that will converge to the largest eigenvalue. Let’s see the following how the power method works.

Consider an × matrix  that has  linearly independent real eigenvalues 1,2,, and the corresponding eigenvectors 1,2,,. Since the eigenvalues are scalars, we can rank them so that |1|>|2|>>|| (actually, we only require |1|>|2|, other eigenvalues may be equal to each other).

Because the eigenvectors are independent, they are a set of basis vectors, which means that any vector that is in the same space can be written as a linear combination of the basis vectors. That is, for any vector 0, it can be written as:

0=11+22++

where 10 is the constraint. If it is zero, then we need to choose another initial vector so that 10.

Now let’s multiply both sides by :

0=11+22++

Since =, we will have:

0=111+222++

We can change the above equation to:

0=11[1+21212++11]=111

where 1 is a new vector and 1=1+21212++11.

This finishes the first iteration. And we can multiply  to 1 to start the 2nd iteration:

1=11+212212++121

Similarly, we can rearrange the above equation to:

1=1[1+2122122++1212]=12

where 2 is another new vector and 2=1+2122122++1212.

We can continue multiply  with the new vector we get from each iteration  times:

1=1[1+21212++11]=1

Because 1 is the largest eigenvalue, therefore, the ratio 1<1 for all >1. Thus when we increase  to sufficient large, the ratio of (1) will be close to 0. So that all the terms that contain this ratio can be neglected as  grows:

1=11

Essentially, as  is large enough, we will get the largest eigenvalue and its corresponding eigenvector. When implementing this power method, we usually normalize the resulting vector in each iteration. This can be done by factoring out the largest element in the vector, which will make the largest element in the vector equal to 1. This normalization will get us the largest eigenvalue and its corresponding eigenvector at the same time. Let’s take a look of the following example.

You may ask when should we stop the iteration? The basic stopping criteria should be one of the three: in the consecutive iterations, (1) the difference between eigenvalues is less than some specified tolerance; (2) the angle between eigenvectors is smaller than a threshold ; or the norm of the residual vector is small enough.

TRY IT! We know from last section that the largest eigenvalue is 4 for matrix =[0223], now use the power method to find the largest eigenvalue and the associated eigenvector. You can use the initial vector [1, 1] to start the iteration.

1st iteration: $$

(7)
[0223]
(8)
[11]

=\begin{bmatrix} 2\5\ \end{bmatrix} =5\begin{bmatrix} 0.4\1\ \end{bmatrix} $$

2nd iteration: $$

(9)
[0223]
(10)
[0.41]

=\begin{bmatrix} 2\3.8\ \end{bmatrix} =3.8\begin{bmatrix} 0.5263\1\ \end{bmatrix} $$

3rd iteration: $$

(11)
[0223]
(12)
[0.52631]

=\begin{bmatrix} 2\ 4.0526\ \end{bmatrix} = 4.0526\begin{bmatrix} 0.4935\1\ \end{bmatrix} $$

4th iteration: $$

(13)
[0223]
(14)
[0.49351]

=\begin{bmatrix} 2\ 3.987\ \end{bmatrix} = 3.987\begin{bmatrix} 0.5016\1\ \end{bmatrix} $$

5th iteration: $$

(15)
[0223]
(16)
[0.50161]

=\begin{bmatrix} 2\ 4.0032\ \end{bmatrix} = 4.0032\begin{bmatrix} 0.4996\1\ \end{bmatrix} $$

6th iteration: $$

(17)
[0223]
(18)
[0.49961]

=\begin{bmatrix} 2\ 3.9992\ \end{bmatrix} = 3.9992\begin{bmatrix} 0.5001\1\ \end{bmatrix} $$

7th iteration: $$

(19)
[0223]
(20)
[0.50011]

=\begin{bmatrix} 2\ 4.0002\ \end{bmatrix} = 4.0002\begin{bmatrix} 0.5000\1\ \end{bmatrix} $$

We can see after 7 iterations, the eigenvalue converged to 4 with [0.5, 1] as the corresponding eigenvector.

TRY IT! Implement the power method in Python

import numpy as np
Copy to clipboard
def normalize(x):
    fac = abs(x).max()
    x_n = x / x.max()
    return fac, x_n
Copy to clipboard
x = np.array([1, 1])
a = np.array([[0, 2], 
              [2, 3]])

for i in range(8):
    x = np.dot(a, x)
    lambda_1, x = normalize(x)
    
print('Eigenvalue:', lambda_1)
print('Eigenvector:', x)
Copy to clipboard
Eigenvalue: 3.999949137887188
Eigenvector: [0.50000636 1.        ]
Copy to clipboard

The inverse power method

The eigenvalues of the inverse matrix 1 are the reciprocals of the eigenvalues of . We can take advantage of this feature as well as the power method to get the smallest eigenvalue of , this will be basis of the inverse power method. The steps are very simple, instead of multiplying  as described above, we just multiply 1 for our iteration to find the largest value of 11, which will be the smallest value of the eigenvalues for . As for the inverse of the matrix, in practice, we can use the methods we covered in the previous chapter to calculate it. We won’t got to the details here, but let’s see an example.

TRY IT! Find the smallest eigenvalue and eigenvector for =[0223].

from numpy.linalg import inv
Copy to clipboard
a_inv = inv(a)

for i in range(8):
    x = np.dot(a_inv, x)
    lambda_1, x = normalize(x)
    
print('Eigenvalue:', lambda_1)
print('Eigenvector:', x)
Copy to clipboard
Eigenvalue: 0.20000000000003912
Eigenvector: [1. 1.]
Copy to clipboard

The shifted power method

In some cases, we need to find all the eigenvalues and eigenvectors instead of the largest and smallest. One simple but inefficient way is to use the shifted power method (we will introduce you an efficient way in next section). Given =, and 1 is the largest eigenvalue obtained by the power method, then we can have:

[1]=

where ’s are the eigenvalues of the shifted matrix 1, which will be 0,21,31,,1.

Now if we apply the power method to the shifted matrix, then we can determine the largest eigenvalue of the shifted matrix, i.e. . Since =1, we can get the eigenvalue  easily. We can repeat this process many times to find the all the other eigenvalues. But you can see that, it involves a lot of work! A better method for finding all the eigenvalues is to use the QR method, let’s see the next section how it works!

15.1 Mathematical Characteristics of Eigen-problems | Contents | 15.3 The QR Method >

No hay comentarios:

Publicar un comentario