下面的代碼提供了兩個函數
import numpy as np
# coefficients of givens matrix
def Givens(alpha: np.float64, beta: np.float64):
eta = np.sqrt(alpha * alpha + beta * beta)
if eta != np.float64(0):
c = alpha / eta
s = beta / eta
else:
c = 1
s = 0
return c, s, eta
# modify hess_mat in place, return Q
# Q = G1^T * G2^T * ... * Gn^T
def upperHessenbergQR(hess_mat: np.array):
size = np.size(hess_mat, 0)
Q = np.identity(size)
for i in range(size - 1):
c, s, eta = Givens(hess_mat[i][i], hess_mat[i + 1][i])
if eta != 0:
hess_mat[i][i] = eta
hess_mat[i + 1][i] = 0
tmp_1 = c * hess_mat[i, i + 1:] + s * hess_mat[i + 1, i + 1:]
tmp_2 = -s * hess_mat[i, i + 1:] + c * hess_mat[i + 1, i + 1:]
hess_mat[i, i + 1:] = tmp_1
hess_mat[i + 1, i + 1:] = tmp_2
tmp_3 = np.matmul(Q[:, i:i + 2], np.array([[c, -s], [s, c]]))
Q[:, i:i + 2] = tmp_3
return Q
def HouseHolder(u: np.array):
gamma = np.sqrt(np.sum(u * u))
w = np.zeros(np.size(u), dtype=np.float64)
if gamma != 0 and gamma != u[0]:
detla = np.sqrt(gamma * (gamma - u[0]))
w[0] = (u[0] - gamma) / detla
w[1:] = u[1:] / detla
else:
gamma = 0
return gamma, w
def QRHouseHolder(mat: np.array):
size = np.size(mat, 0)
Q = np.identity(size, dtype=np.float64)
for i in range(size - 1):
gamma, u = HouseHolder(mat[i:, i])
temp = np.identity(len(u), dtype=np.float64) - np.matmul(u.reshape(-1, 1), u.reshape(1, -1))
H = np.identity(size, dtype=np.float64)
H[i:, i:] = temp
if gamma != 0:
Q = np.matmul(H, Q)
v = np.matmul(u, mat[i:, i + 1:])
mat[i][i] = gamma
mat[i + 1:, i] = 0
for j in range(i, size):
mat[j, i + 1:] = mat[j, i + 1:] - u[j - i] * v
else:
mat[i + 1:, i] = 0
return Q.T
# solve linear equation by QR decomposition
def solve_linear_equ(mat_copy: np.array, constants: np.array):
R = mat_copy.copy()
Q = QRHouseHolder(R)
right = np.matmul(Q.T, constants.reshape(-1, 1)).reshape(-1)
size = np.size(constants)
roots = np.zeros(size)
roots[size - 1] = 1 if R[size - 1][size - 1] == 0 else right[size - 1] / R[size - 1][size - 1]
for i in range(size - 1, -1, -1):
tmp = right[i] - np.sum(R[i, i + 1:] * roots[i + 1:])
roots[i] = 1 if R[i][i] == 0 else tmp / R[i][i]
return roots
# transform the matrix into hessenberg matrix
def HouseholderReduction(mat: np.array):
size = np.size(mat, 0)
for i in range(size - 2):
gamma, u = HouseHolder(mat[i + 1:, i])
if gamma != 0:
v = np.matmul(u, mat[i + 1:, i + 1:])
mat[i + 1][i] = gamma
mat[i + 2:, i] = 0
for j in range(i + 1, size):
mat[j, i + 1:] = mat[j, i + 1:] - u[j - i - 1] * v
v = np.matmul(mat[:, i + 1:], u)
for j in range(i + 1, size):
mat[:, j] = mat[:, j] - u[j - i - 1] * v
else:
mat[i + 2:, i] = 0
""" def HouseholderReduction2(mat: np.array): size = np.size(mat, 0) First = 1 for i in range(size - 2): gamma, u = HouseHolder(mat[i + 1:, i]) if gamma != 0: temp = np.identity(len(u)) - u.reshape(-1, 1) * u.reshape(1, -1) H = np.identity(size) H[i + 1:, i + 1:] = temp if First == 0: print(np.matmul(H, mat)) mat = np.matmul(np.matmul(H, mat), H) mat[i + 2:, i] = 0 if First == 0: print(mat) First = 1 return mat """
def check_converge(mat: np.array, ord: int):
return np.abs(mat[ord][ord - 1]) < 1.0e-10
def solve_eigen_vec(mat: np.array, size: int, ord: int):
roots = np.zeros(size, dtype=np.float64)
roots[ord] = 1
for inx in range(ord - 1, -1, -1):
v = np.sum(mat[inx, inx + 1:] * roots[inx + 1:])
coeff = mat[ord][ord] - mat[inx][inx]
if coeff == 0:
roots[inx] = 1
else:
roots[inx] = v / coeff
return roots / np.sqrt(np.sum(roots * roots))
def eigen(m: np.array):
mat = m.copy()
size = np.size(mat, 0)
ord = size - 1
HouseholderReduction(mat)
mat_copy = mat.copy()
cnt = 0
while ord:
u = mat[ord][ord] if cnt > 100 else 0
cnt += 1
mat -= u * np.identity(size, dtype=np.float64)
q = upperHessenbergQR(mat)
mat = np.matmul(mat, q) + u * np.identity(size, dtype=np.float64)
if check_converge(mat, ord):
mat[ord][ord - 1] = 0
ord -= 1
eigen_value = np.array([mat[i][i] for i in range(size)])
# print(eigen_value)
ord = size - 1
ret_q = np.identity(size)
while ord:
mat_copy -= eigen_value[ord] * np.identity(size)
q = upperHessenbergQR(mat_copy)
ret_q = np.matmul(ret_q, q)
mat_copy = np.matmul(mat_copy, q) + eigen_value[ord] * np.identity(size)
if check_converge(mat_copy, ord):
mat_copy[ord][ord - 1] = 0
ord -= 1
eigen_value = np.array([mat_copy[i][i] for i in range(size)])
eigen_vec = np.array(
[np.matmul(ret_q, solve_eigen_vec(mat_copy, size, i).reshape(-1, 1)).reshape(-1) for i in
range(size)])
"""e, v = np.linalg.eig(mat_copy) #print(e) #print(v) eigen_vec = np.array( [solve_linear_equ(mat_copy, size, i) for i in range(size)]) #print(eigen_value) #print(eigen_vec)"""
return eigen_value, eigen_vec
""" 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 check_converge(mat: np.array, size: int): for i in range(1, size): for j in range(i - 1): if mat[i][j] > 1.0e-14: return False return True def eigen(mat: np.array): size = np.size(mat, 0) while 1: q, r = HouseHolder(mat) new_mat = r * q if check_converge(new_mat, size): print(new_mat) return np.array([new_mat[i][i] for i in range(size)]), np.array([]) mat = new_mat y = np.array([[2, -1, 0, -1], [-1, 2, -1, 0], [0, -1, 2, -1], [-1, 0, -1, 2]]) h = y[1:, 1:] print(h) h[0][0] = 100 """
""" c, _ = eigen(y) print(c) t, _ = np.linalg.eig(y) print(t) """
""" y = np.array([[1, 2, 3] for i in range(3)]) q, r = HouseHolder(y) print(q) print(r) print(np.matmul(q, r)) """
# y = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]], dtype=np.float64)
# H, Q = hessenberg(y, calc_q=True)
# H_copy = H.copy()
# print(H)
# q, r = qr(H)
# mat = HouseholderReduction2(y)
# print(mat)
# print(q)
# print(r)
# print(np.matmul(q, r))
# my_q = upperHessenbergQR(H)
# print(my_q)
# print(q)
# print(my_q)
# print(np.matmul(my_q, H))
if __name__ == '__main__':
y = np.array([[2, -1, 0, -1],
[-1, 2, -1, 0],
[0, -1, 2, -1],
[-1, 0, -1, 2]], dtype=np.float64)
y = np.array([[1, 1, 1, -2], [1, 1, 0, 3], [1, 2, 1, 2.7], [-1, 3, 2, 0.9]], dtype=np.float64)
ret = np.array([1.3, 2.7, 1, -1])
print(y)
print(solve_linear_equ(y, ret))
print(np.matmul(np.linalg.inv(y), ret.reshape(-1, 1)).reshape(-1))
注：
