QR分解
实数域
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}),x\in\mathrm{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 \ne1\(
\)\therefore v_1=(x_1-k)/2K\(
\)v_i = x_i/2K\(
\)\because Px=x-2Kv=ke, 2Kv=x-ke\(
\)\therefore 2K2=S2-x_1k=S^2\pm x_1S\(
\)if \quad u:=(x_1 \pm S, x_2, x_3,\ldots, x_d)\(
\)then, P = (I-uuT/2K2)$
正负号的选择
\(k=\pm S\)
\(2K^2=S^2\pm x_1S\)
\(如何选择正负号呢，为了避免K=0,我们选择符号为x_1的符号即可。\)
多步
\(假设矩阵A的前r列已经是上三角形，则A_r如下\)
\(为了把对角化剩下的，我们只需取：\)
二者相乘，我们容易发现，右下角变成了第一步的情况，依此知道对角化完全。
可得:
\(A=QR\)
\(Q=P_1P_2\ldots\)
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
