使用Gram-Schmidt正交化，Householder变换，Given旋转三种方法实现了QR分解，但是感觉好像并无卵用，貌似实际生产有更好的改进方法？
矩阵分析或者矩阵论的课都有介绍，网上资料也很多，在此不贴了。
import numpy as np
def SchmitOrth(mat:np.array):
cols = mat.shape[1]
Q = np.copy(mat)
R = np.zeros((cols, cols))
for col in range(cols):
for i in range(col):
k = np.sum(mat[:, col] * Q[:, i]) / np.sum( np.square(Q[:, i]) )
Q[:, col] -= k*Q[:, i]
Q[:, col] /= np.linalg.norm(Q[:, col])
for i in range(cols):
R[col, i] = Q[:, col].dot( mat[:, i] )
return Q, R
def HouseHolder(mat:np.array):
cols = mat.shape[1]
Q = np.eye(cols)
R = np.copy(mat)
for col in range(cols-1):
a = np.linalg.norm(R[col:, col])
e = np.zeros((cols- col))
e[0] = 1.0
num = R[col:, col] -a*e
den = np.linalg.norm(num)
u = num / den
H = np.eye((cols))
H[col:, col:] = np.eye((cols- col))- 2*u.reshape(-1, 1).dot(u.reshape(1, -1))
R = H.dot(R)
Q = Q.dot(H)
return Q, R
def GivenRot(mat:np.array):
rows, cols = mat.shape
R = np.copy(mat)
Q = np.eye(cols)
for col in range(cols):
for row in range(col+1, rows):
if abs(R[row, col]) < 1e-6:
continue
f = R[col, col]
s = R[row, col]
den = np.sqrt( f*f+ s*s)
c = f / den
s = s / den
T = np.eye(rows)
T[col, col], T[row, row] = c, c
T[row, col], T[col, row] = -s, s
R = T.dot(R)
Q = T.dot(Q)
return Q.T, R
mat = np.array( [ [1.0, 2.0, 2.0],
[2.0, 1.0, 2.0],
[1.0, 2.0, 1.0], ])
q, r = SchmitOrth(mat)
print("Gram-Schmidt Orthogonal")
print("Q: \n", q)
print("R: \n", r)
print("Q x R: \n", q.dot(r))
#
q, r = HouseHolder(mat)
print("HouseHolder")
print("Q: \n", q)
print("R: \n", r)
print("Q x R: \n", q.dot(r))
#
q, r = GivenRot(mat)
print("Givens")
print("Q: \n", q)
print("R: \n", r)
print("Q x R: \n", q.dot(r))
