## 1.Gram-Schmidt正交化

$$B=b-Pb=b-\frac{A^Tb}{A^TA}A$$

## 3.Householder矩阵与Householder变换

 1 #coding:utf8
2 import numpy as np
3
4 def gram_schmidt(A):
5     """Gram-schmidt正交化"""
6     Q=np.zeros_like(A)
7     cnt = 0
8     for a in A.T:
9         u = np.copy(a)
10         for i in range(0, cnt):
11             u -= np.dot(np.dot(Q[:, i].T, a), Q[:, i]) # 减去待求向量在以求向量上的投影
12         e = u / np.linalg.norm(u)  # 归一化
13         Q[:, cnt] = e
14         cnt += 1
15     R = np.dot(Q.T, A)
16     return (Q, R)
17
18 def givens_rotation(A):
19     """Givens变换"""
20     (r, c) = np.shape(A)
21     Q = np.identity(r)
22     R = np.copy(A)
23     (rows, cols) = np.tril_indices(r, -1, c)
24     for (row, col) in zip(rows, cols):
25         if R[row, col] != 0:  # R[row, col]=0则c=1,s=0,R、Q不变
26             r_ = np.hypot(R[col, col], R[row, col])  # d
27             c = R[col, col]/r_
28             s = -R[row, col]/r_
29             G = np.identity(r)
30             G[[col, row], [col, row]] = c
31             G[row, col] = s
32             G[col, row] = -s
33             R = np.dot(G, R)  # R=G(n-1,n)*...*G(2n)*...*G(23,1n)*...*G(12)*A
34             Q = np.dot(Q, G.T)  # Q=G(n-1,n).T*...*G(2n).T*...*G(23,1n).T*...*G(12).T
35     return (Q, R)
36
37 def householder_reflection(A):
38     """Householder变换"""
39     (r, c) = np.shape(A)
40     Q = np.identity(r)
41     R = np.copy(A)
42     for cnt in range(r - 1):
43         x = R[cnt:, cnt]
44         e = np.zeros_like(x)
45         e[0] = np.linalg.norm(x)
46         u = x - e
47         v = u / np.linalg.norm(u)
48         Q_cnt = np.identity(r)
49         Q_cnt[cnt:, cnt:] -= 2.0 * np.outer(v, v)
50         R = np.dot(Q_cnt, R)  # R=H(n-1)*...*H(2)*H(1)*A
51         Q = np.dot(Q, Q_cnt)  # Q=H(n-1)*...*H(2)*H(1)  H为自逆矩阵
52     return (Q, R)
53
54 np.set_printoptions(precision=4, suppress=True)
55 A = np.array([[6, 5, 0],[5, -1, 4],[5, 1, -14],[0, 4, 3]],dtype=float)
56
57 (Q, R) = gram_schmidt(A)
58 print(Q)
59 print(R)
60 print np.dot(Q,R)
61
62 (Q, R) = givens_rotation(A)
63 print(Q)
64 print(R)
65 print np.dot(Q,R)
66
67 (Q, R) = householder_reflection(A)
68 print(Q)
69 print(R)
70 print np.dot(Q,R)

posted on 2016-11-18 21:53  1357  阅读(18049)  评论(0编辑  收藏  举报