QR分解迭代求特征值——原生python实现(不使用numpy)

时间:2023-03-09 21:18:44
QR分解迭代求特征值——原生python实现(不使用numpy)

QR分解:

QR分解迭代求特征值——原生python实现(不使用numpy)

有很多方法可以进行QR迭代,本文使用的是Schmidt正交化方法

具体证明请参考链接 https://wenku.baidu.com/view/c2e34678168884868762d6f9.html

迭代格式

QR分解迭代求特征值——原生python实现(不使用numpy)

实际在进行QR分解之前一般将矩阵化为上hessnberg矩阵(奈何这个过程比较难以理解,本人智商不够,就不做这一步了哈哈哈)

迭代终止条件

看了很多文章都是设置一个迭代次数,感觉有些不是很合理,本来想采用A(k+1)-A(k)的对角线元素的二范数来作为误差的,但是我有没有一些严格的证明,所以本文也采用比较大众化的思路,设置迭代次数。

Python实现

 M = [[2, 4, 2], [-1, 0, -4], [2, 2, 1]]

 import copy
import math class QR(object): def __init__(self, data):
self.M = data
self.degree = len(data) def get_row(self, index):
res = []
for i in range(self.degree):
res.append(self.M[i][index])
return res def get_col(self, index):
res = []
for i in range(self.degree):
res.append(self.M[i][index])
return res @staticmethod
def dot(m1, m2):
res = 0
for i in range(len(m1)):
res += m1[i] * m2[i]
return res @staticmethod
def list_multi(k, lt):
res = []
for i in range(len(lt)):
res.append(k * lt[i])
return res @staticmethod
def one_item(x, yArr):
res = [0 for i in range(len(x))]
temp_y_arr = [] n = len(yArr)
if n == 0:
res = x
else:
for item in yArr:
k = QR.dot(x, item) / QR.dot(item, item)
temp_y_arr.append(QR.list_multi(-k, item))
temp_y_arr.append(x) for item in temp_y_arr:
for i in range(len(item)):
res[i] += item[i]
return res @staticmethod
def normal(matrix):
yArr = []
yArr.append(matrix[0]) for i in range(1, len(matrix)):
yArr.append(QR.one_item(matrix[i], yArr))
return yArr @staticmethod
def normalized(lt):
res = []
sm = 0
for item in lt:
sm += math.pow(item, 2)
sm = math.sqrt(sm)
for item in lt:
res.append(item / sm)
return res @staticmethod
def matrix_T(matrix):
mat = copy.deepcopy(matrix)
m = len(mat[0])
n = len(mat)
for i in range(m):
for j in range(n):
if i < j:
temp = mat[i][j]
mat[i][j] = mat[j][i]
mat[j][i] = temp
return mat @staticmethod
def matrix_multi(mat1, mat2):
res = []
rows = len(mat1[0])
cols = len(mat1)
for i in range(rows):
temp = [0 for i in range(cols)]
res.append(temp) for i in range(rows):
for j in range(cols):
sm = 0
for k in range(cols):
sm += (mat1[k][i] * mat2[j][k])
res[j][i] = sm
return res def execute(self):
xArr = []
for i in range(self.degree):
xArr.append(self.get_col(i))
yArr = QR.normal(xArr)
self.Q = []
for item in yArr:
self.Q.append(QR.normalized(item)) self.R = QR.matrix_multi(QR.matrix_T(self.Q), xArr)
return (self.Q, self.R) # A = [
# [1, 0, -1, 2, 1],
# [3, 2, -3, 5, -3],
# [2, 2, 1, 4, -2],
# [0, 4, 3, 3, 1],
# [1, 0, 8, -11, 4]
# ]
# A = [
# [1, 2, 2],
# [2, 1, 2],
# [2, 2, 1]
# ]
A = [
[3, 2, 4],
[2, 0, 2],
[4, 2, 3]
] temp = copy.deepcopy(A)
val = [] # 特征值
times = 20 # 迭代次数
for i in range(times):
qr = QR(temp)
(q, r) = qr.execute()
temp = QR.matrix_multi(r, q)
temp = QR.matrix_T(temp) for i in range(len(temp)):
for j in range(len(temp[0])):
if i == j:
val.append(temp[i][j])
# 特征值
print(val)

结果展示

QR分解迭代求特征值——原生python实现(不使用numpy)

总结

使用QR分解迭代求特征值,收敛的比较快,也可以求出所有的特征值,但是如果要求特征向量的话,还是需要求解线性方程组(感觉很麻烦)