Python解线性方程组的直接法（基于吉文斯变换的QR分解）

程序员文章站 2022-04-18 09:26:10

基于吉文斯变换的QR分解# 基于吉文斯变换的QR分解import numpy as npdef givens(A): m = A.shape[0] n = A.shape[1] P = np.mat(np.eye(n)) p1 = P for i in range(0, n-1): for j in range(i+1, n): if A[i, j] == 0: continue...

基于吉文斯变换的QR分解

#  基于吉文斯变换的QR分解 import numpy as np def givens(A): m = A.shape[0] n = A.shape[1] P = np.mat(np.eye(n)) p1 = P for i in range(0, n-1): for j in range(i+1, n): if A[i, j] == 0: continue else: c1 = np.sqrt(np.sum(np.multiply(A[i:j, i], A[i:j, i]))) d = np.sqrt(np.sum(np.multiply(A[i:j+1, i], A[i:j+1, i]))) s1 = A[j, i] c = c1/d
                s = s1/d
                P[i, i] = c
                P[i, j] = s
                P[j, i] = -s
                P[j, j] = c
                A = P*A
                p1 = P * p1
                P = np.mat(np.eye(n)) for i in range(m): for j in range(n): if abs(A[i, j]) < 10**-10: A[i, j] = 0 return A, p1.T if __name__ == "__main__": ai = np.mat([[3, 5, 5], [0, 3, 4], [4, 0, 5]], dtype=float) aa, po = givens(ai) print('基于吉文斯变换的QR分解:') print('R:') print(aa) print('Q:') print(po)