实数域
Hermite矩阵
(满足mathrm{P^{H}=P})
(P = (I-2vv^{T})为初等正交Hermite（对称）矩阵)
(即：P^{T}=P,P^TP=I)
QR分解
Q为正交矩阵 R为上三角阵
一步
(先考虑Px=ke,P=(I-2vv^{T}),xinmathrm{R}^{d}的情况下，求解P。)
(Px=ke)
(k^2=|Px|_2^2=x^2=x_1^2+x_2^2+ldots+x_d^2)
(S:=sqrt{x_1^2+x_2^2+ldots+x_d^2})
(then, k=pm S)
(K:=v^Tx)
$ x-2Kv=ke(
)then, x_1-2Kv_1=k (
)x_i-2Kv_i=0, quad i
e1(
) herefore v_1=(x_1-k)/2K(
)v_i = x_i/2K(
)ecause Px=x-2Kv=ke, 2Kv=x-ke(
) herefore 2K^2=S^2-x_1k=S^2pm x_1S(
)if quad u:=(x_1 pm S, x_2, x_3,ldots, x_d)(
)then, P = (I-uu^T/2K^2)$
正负号的选择
(k=pm S)
(2K^2=S^2pm x_1S)
(如何选择正负号呢，为了避免K=0,我们选择符号为x_1的符号即可。)
多步
(假设矩阵A的前r列已经是上三角形，则A_r如下)
(为了把对角化剩下的，我们只需取：)
二者相乘，我们容易发现，右下角变成了第一步的情况，依此知道对角化完全。
可得:
(A=QR)
(Q=P_1P_2ldots)
Tips 子空间
(V=QR)
(Q=(q_1, q_2, ldots, q_r, ldots, q_d))
(V=(v_1, v_2, ldots, v_r))
(如果v_1, v_2, ldots,v_r的极大线性无关组为本身，那么)
(q_1,q_2,ldots,q_r与v_1, v_2, ldots, v_r张成同一个子空间（互相线性表出）。)
代码
import numpy as np
def step_one_QR(x):
"""
一步
x 为ndarray
:param x:向量
:return: P
"""
u = np.array(x, dtype=float)
sign = lambda x: 1 if x >= 0 else -1
S_2 = u @ u
S = np.sqrt(S_2) * sign(u[0])
K_2 = S_2 + S * u[0] #K_2 = 2K^2 here
vector_ones = np.ones(len(u), dtype=float)
u[0] += S
return np.diag(vector_ones) - np.outer(u, u) / K_2
def step_all_QR(A):
"""
d:A的行数
r:A的列数
P:每步得到的Hermite矩阵
:param A:d x r的矩阵
:return:Q, R(A)
"""
A = np.array(A, dtype=float) #注意经此操作后，后续操作不会改变原来的A
d, r = A.shape
Q = np.diag(np.ones(d, dtype=float))
for i in range(r):
P = np.diag(np.ones(d, dtype=float))
P[i:, i:] = step_one_QR(A[i:,i])
Q = Q @ P
A = P @ A
return Q, A
"""
对上面的一个改进，速度明显提高，不过
俩个方法对列非满秩的抗性都不高。
"""
def step_one_QR_2(x):
"""
一步 第二种
x 为ndarray
:param x:向量
:return:u, K_2
"""
u = np.array(x, dtype=float)
sign = lambda x: 1 if x >= 0 else -1
S_2 = u @ u
S = np.sqrt(S_2) * sign(u[0])
K_2 = S_2 + S * u[0] #K_2 = 2K^2 here
vector_ones = np.ones(len(u), dtype=float)
u[0] += S
return u, K_2
def step_all_QR_2(A):
"""
d:A的行数
r:A的列数
:param A:
:return:
"""
A = np.array(A, dtype=float) # 注意经此操作后，后续操作不会改变原来的A
d, r = A.shape
Q = np.diag(np.ones(d, dtype=float))
for i in range(r):
u, k = step_one_QR_2(A[i:, i])
if k <= 0.0000001:
print("Matrix may not full rank...")
return Q, A #A非列满秩
Q[:, i:] -= np.outer(Q[:, i:] @ u, u / k)
A[i:, :] -= np.outer(u / k, u @ A[i:, :])
return Q, A