實數域
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 2K^2=S^2-x_1k=S^2\pm 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^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