实数域
Hermite矩阵
满
足
P
H
=
P
满足\mathrm{P^{H}=P}
满足PH=P
P
=
(
I
−
2
v
v
T
)
为
初
等
正
交
H
e
r
m
i
t
e
(
对
称
)
矩
阵
P = (I-2vv^{T})为初等正交Hermite(对称)矩阵
P=(I−2vvT)为初等正交Hermite(对称)矩阵
即
:
P
T
=
P
,
P
T
P
=
I
即:P^{T}=P,P^TP=I
即:PT=P,PTP=I
QR分解
Q为正交矩阵 R为上三角阵
一步
先
考
虑
P
x
=
k
e
,
P
=
(
I
−
2
v
v
T
)
,
x
∈
R
d
的
情
况
下
,
求
解
P
。
先考虑Px=ke,P=(I-2vv^{T}),x\in\mathrm{R}^{d}的情况下,求解P。
先考虑Px=ke,P=(I−2vvT),x∈Rd的情况下,求解P。
P
x
=
k
e
Px=ke
Px=ke
k
2
=
∥
P
x
∥
2
2
=
x
2
=
x
1
2
+
x
2
2
+
…
+
x
d
2
k^2=\|Px\|_2^2=x^2=x_1^2+x_2^2+\ldots+x_d^2
k2=∥Px∥22=x2=x12+x22+…+xd2
S
:
=
x
1
2
+
x
2
2
+
…
+
x
d
2
S:=\sqrt{x_1^2+x_2^2+\ldots+x_d^2}
S:=x12+x22+…+xd2
t
h
e
n
,
k
=
±
S
then, k=\pm S
then,k=±S
K
:
=
v
T
x
K:=v^Tx
K:=vTx
$ x-2Kv=ke$
$then, x_1-2Kv_1=k $
x
i
−
2
K
v
i
=
0
,
i
≠
1
x_i-2Kv_i=0, \quad i \ne1
xi−2Kvi=0,i̸=1
∴
v
1
=
(
x
1
−
k
)
/
2
K
\therefore v_1=(x_1-k)/2K
∴v1=(x1−k)/2K
v
i
=
x
i
/
2
K
v_i = x_i/2K
vi=xi/2K
∵
P
x
=
x
−
2
K
v
=
k
e
,
2
K
v
=
x
−
k
e
\because Px=x-2Kv=ke, 2Kv=x-ke
∵Px=x−2Kv=ke,2Kv=x−ke
∴
2
K
2
=
S
2
−
x
1
k
=
S
2
±
x
1
S
\therefore 2K^2=S^2-x_1k=S^2\pm x_1S
∴2K2=S2−x1k=S2±x1S
i
f
u
:
=
(
x
1
±
S
,
x
2
,
x
3
,
…
,
x
d
)
if \quad u:=(x_1 \pm S, x_2, x_3,\ldots, x_d)
ifu:=(x1±S,x2,x3,…,xd)
t
h
e
n
,
P
=
(
I
−
u
u
T
/
2
K
2
)
then, P = (I-uu^T/2K^2)
then,P=(I−uuT/2K2)
正负号的选择
k
=
±
S
k=\pm S
k=±S
2
K
2
=
S
2
±
x
1
S
2K^2=S^2\pm x_1S
2K2=S2±x1S
如
何
选
择
正
负
号
呢
,
为
了
避
免
K
=
0
,
我
们
选
择
符
号
为
x
1
的
符
号
即
可
。
如何选择正负号呢,为了避免K=0,我们选择符号为x_1的符号即可。
如何选择正负号呢,为了避免K=0,我们选择符号为x1的符号即可。
多步
假
设
矩
阵
A
的
前
r
列
已
经
是
上
三
角
形
,
则
A
r
如
下
假设矩阵A的前r列已经是上三角形,则A_r如下
假设矩阵A的前r列已经是上三角形,则Ar如下
为
了
把
对
角
化
剩
下
的
,
我
们
只
需
取
:
为了把对角化剩下的,我们只需取:
为了把对角化剩下的,我们只需取:
二者相乘,我们容易发现,右下角变成了第一步的情况,依此知道对角化完全。
可得:
A
=
Q
R
A=QR
A=QR
Q
=
P
1
P
2
…
Q=P_1P_2\ldots
Q=P1P2…
Tips 子空间
V
=
Q
R
V=QR
V=QR
Q
=
(
q
1
,
q
2
,
…
,
q
r
,
…
,
q
d
)
Q=(q_1, q_2, \ldots, q_r, \ldots, q_d)
Q=(q1,q2,…,qr,…,qd)
V
=
(
v
1
,
v
2
,
…
,
v
r
)
V=(v_1, v_2, \ldots, v_r)
V=(v1,v2,…,vr)
如
果
v
1
,
v
2
,
…
,
v
r
的
极
大
线
性
无
关
组
为
本
身
,
那
么
如果v_1, v_2, \ldots,v_r的极大线性无关组为本身,那么
如果v1,v2,…,vr的极大线性无关组为本身,那么
q
1
,
q
2
,
…
,
q
r
与
v
1
,
v
2
,
…
,
v
r
张
成
同
一
个
子
空
间
(
互
相
线
性
表
出
)
。
q_1,q_2,\ldots,q_r与v_1, v_2, \ldots, v_r张成同一个子空间(互相线性表出)。
q1,q2,…,qr与v1,v2,…,vr张成同一个子空间(互相线性表出)。
代码
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