本站原创文章，转载请说明来自《老饼讲解-机器学习》www.bbbdata.com
本节讲解QR分解是什么，然后展示如何进行QR分解和分解的原理
最后，文未给出python实现QR分解的代码，进一步展示QR分解的细节
本节讲解QR分解是什么和如何进行QR分解
QR分解的定义
对于m*n列满秩矩阵A
可将其分解为列正交矩阵与上三角矩阵的积,
即
其中
👉 Q列与列之间两两正交
👉 R为非奇异上三角矩阵
QR分解的方法
符号说明
记的第i列为 ,Q的第i列为
即
Q的求法
Q中的求解如下
PASS
实际就是新的减去在各个之前已经正交化的上的投影分量,
从而得到由新 生成的新正交量
R的求法
在求得Q后，R可如下算得
✍️补充：关于Q的单位化
如果需要进一步单位化Q,
使Q的每一条列向量 为单位向量,
则可以在以上结果进一步如下单位化Q
Q的单位化公式
对应的R为
本节讲解上节所提供的QR分解公式的原理推导
相关引理
引理1：两条向量的正交
求b对a的正交向量d:
先求出b在a的投影向量c
则b对a的正交向量 为
引理2：与正交向量集的正交量
求 b对的正交量 d,使 d与都正交
其中, 是正交向量集,它的任意两条向量 正交
将 b可以减去在上的投影分量分量即可
其中
为b在上的投影分量
对于任意一个, d 都与正交,因为 中除了 和项,其余项都与正交,而b-k_ia_i又与a_i正交,所以d与a_i正交
即：
问题结论推导
Q的求法
因此,要用生成正交向量 ,
只需令,然后不断的用新的来生成与已有正交向量集的正交量即可
即可得到
....
R的求法
又由上,可得
...
即：
可知
3、单位化QR
可得,单位化的Q为：
对应的R为：
本节展示如何用python代码实现QR分解，并展示运行结果
python实现矩阵的QR分解
# -*- coding: utf-8 -*-
"""
矩阵的QR分解
"""
import numpy as np
# 将矩阵进行QR分解
def QrDecompose(A):
#------------QR分解----------------------------
Q = np.zeros(A.shape) # 初始化Q
Q[:,0] = A[:,0] # A的第0列作为Q的第0列
R=np.identity(A.shape[1]) # 将R初始化为单位矩阵
for i in range(1,A.shape[1]): # 逐列计算Q、R
cur_beta = A[:,i].copy() # 初始化beta为Ai
for j in range(i): # 将Ai与Q[:,i-1]进行正交
k = A[:,i]@Q[:,j]/(Q[:,j].T@Q[:,j]) # 计算Ai在Qj的投影系数
R[j,i]= k # Rji就是Ai在Qj的投影系数
cur_beta -= k*Q[:,j] # 更新beta,令beta与Qj正交，即减去Ai在Qj的投影
Q[:,i] =cur_beta # 将最终的beta(即与Q[:,i-1]正交后的Ai)作为Qj
return Q,R
# 将QR矩阵进行单位化
def QrNormal(Q,R):
c = Q.shape[1] # Q的列数
D1 = np.zeros((c,c)) # 初始化对角矩阵D1
D2 = np.zeros((c,c)) # 初始化对角矩阵D2
for i in range(c): # 计算D的对角元素
D1[i,i] = np.sqrt(Q[:,i].T@Q[:,i]) # 计算D1的第i个对角元素,即|Q第i列|
D2[i,i] = 1/D1[i,i] # 计算D2的第i个对角元素,即1/|Q第i列|
nQ = Q@D2 # 计算单位化后的Q
nR = D1@R # 计算单位化后的R
return nQ,nR
# 测试样例
if __name__ == "__main__":
A = np.array([[1.,2.,5,8],[3.,5.,4,2],[6.,4,3,1]]).T # 生成需要QR分解的矩阵A
Q,R = QrDecompose(A) # 将A进行QR分解，得到列正交矩阵Q和上三角矩阵R
nQ,nR= QrNormal(Q,R) # 将QR进行标准化
#------------打印单位化前QR的信息-------------------------
print("----原始数据A与分解后的QR---")
print('A=\n',A)
print('Q=\n',Q)
print('R=\n',R)
print("\n-----用QTQ验证Q各个向量正交--")
QTQ= Q.T@Q
QTQ[np.abs(QTQ)<0.00001]=0
print('Q^T*Q=\n',QTQ)
print("\n-----用A-QR验证A=QR--")
err = A-Q@R
err[np.abs(err)<0.00001]=0
print('A-QR=\n',err)
#----------单位化---------------------------------------------
print("\n----单位化后的QR:nQ,nR---")
print('nQ=\n',nQ)
print('nR=\n',nR)
print("\n-----用nQ^T*nQ验证nQ各个向量正交--")
nQTnQ= nQ.T@nQ
nQTnQ[np.abs(nQTnQ)<0.00001]=0
print('nQ^T*nQ=\n',nQTnQ)
print("-----用A-QR验证A=QR--")
nerr = A-nQ@nR
nerr[np.abs(nerr)<0.00001]=0
print('A-nQ*nR=\n',nerr)代码运行结果
----原始数据A与分解后的QR---
A=
[[1. 3. 6.]
[2. 5. 4.]
[5. 4. 3.]
[8. 2. 1.]]
Q=
[[ 1. 2.4787234 2.75700935]
[ 2. 3.95744681 -1.3364486 ]
[ 5. 1.39361702 -0.57009346]
[ 8. -2.17021277 0.34579439]]
R=
[[1. 0.5212766 0.39361702]
[0. 1. 1.14953271]
[0. 0. 1. ]]
-----用QTQ验证Q各个向量正交--
Q^T*Q=
[[94. 0. 0. ]
[ 0. 28.45744681 0. ]
[ 0. 0. 9.8317757 ]]
-----用A-QR验证A=QR--
A-QR=
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
----单位化后的QR:nQ,nR---
nQ=
[[ 0.10314212 0.46465445 0.87927 ]
[ 0.20628425 0.74185174 -0.42622241]
[ 0.51571062 0.26124349 -0.18181515]
[ 0.825137 -0.40682192 0.11028132]]
nR=
[[9.69535971 5.05396411 3.81625861]
[0. 5.33455217 6.13224221]
[0. 0. 3.13556625]]
-----用nQ^T*nQ验证nQ各个向量正交--
nQ^T*nQ=
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
-----用A-QR验证A=QR--
A-nQ*nR=
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]End