1. 程式人生 > >冪迭代和逆冪迭代求最大最小特徵值

冪迭代和逆冪迭代求最大最小特徵值

參考連結 https://wenku.baidu.com/view/ee7ecbeca98271fe910ef9fc.html?from=search

冪迭代演算法:

逆冪迭代演算法:

實際在使用中A可以先進行LU分解

 

無論是冪迭代,還是逆冪迭代,只能求出最大和最小特徵值與對應的特徵向量。

具體問題

 

python實現

 1 import numpy as np
 2 import matplotlib.pyplot as plt
 3 
 4 eigenval1_ = []
 5 eigenval2_ = []
 6
eigenvec1_ = [] 7 eigenvec2_ = [] 8 cnt_ = [] 9 x0 = np.array([1 / np.sqrt(2), 1 / np.sqrt(2)]) 10 x0.resize((2, 1)) 11 dis = 10 ** -12 12 13 14 def normalized_power_iteration(A): 15 v0 = x0 16 distance = 1 17 times = 0.0 18 19 while distance > dis: 20 v1 = np.dot(A, v0)
21 v1 = v1 / np.linalg.norm(v1) 22 distance = np.linalg.norm(v1 - v0) 23 v0 = v1 24 times += 1.0 25 eig = np.linalg.norm(np.dot(A, v1)) / np.linalg.norm(v1) 26 return (v1, eig, times) 27 28 29 def normalized_inverse_power_iteration(A): 30 # 2x2矩陣 進行LU分解,在這裡進行手動的LU分解求逆,因為題目找那個說了是一個2x2的矩陣,實際上當矩陣的維數比較高的時候,應該用LU分解然後解線性方程組來求解。
31 L = np.array([1, 0, A[1, 0] / A[0, 0], 1]) 32 L.resize((2, 2)) 33 LI = np.array([1, 0, -L[1, 0], 1]) 34 LI.resize((2, 2)) 35 36 U = np.array([A[0, 0], A[0, 1], 0, A[1, 1] - A[0, 1] * A[1, 0] / A[0, 0]]) 37 U.resize((2, 2)) 38 UI = np.array([1 / U[0, 0], -U[0, 1] / (U[0, 0] * U[1, 1]), 0, 1 / U[1, 1]]) 39 UI.resize((2, 2)) 40 lu = np.dot(UI, LI) 41 42 v0 = x0 43 distance = 1 44 times = 0 45 46 while distance > dis: 47 u1 = np.dot(lu, v0) 48 m1 = np.linalg.norm(u1) 49 v1 = u1 / m1 50 distance = np.linalg.norm(v1 - v0) 51 v0 = v1 52 times += 1 53 eig = np.linalg.norm(np.dot(lu, v1)) / np.linalg.norm(v1) 54 return (v1, 1 / eig, times) 55 56 57 for i in np.arange(0, n): 58 (vec, val, round1) = normalized_power_iteration(As[i, :, :]) 59 eigenval1_.append(val) 60 eigenvec1_.append(vec[0, 0]) 61 eigenvec1_.append(vec[1, 0]) 62 cnt_.append(round1) 63 64 (vec, val, round2) = normalized_inverse_power_iteration(As[i, :, :]) 65 eigenvec2_.append(vec[0, 0]) 66 eigenvec2_.append(vec[1, 0]) 67 eigenval2_.append(val) 68 69 eigenval1 = np.array(eigenval1_) 70 eigenvec1 = np.array(eigenvec1_) 71 eigenvec1.resize((n, 2)) 72 73 cnt = np.array(cnt_) 74 75 eigenval2 = np.array(eigenval2_) 76 eigenvec2 = np.array(eigenvec2_) 77 eigenvec2.resize((n, 2)) 78 79 plt.xlabel("bilv") 80 plt.ylabel("round times") 81 plt.title("result") 82 plt.show()