机器学习-统计与数学

【分解】一篇入门之-矩阵的QR分解是什么

作者 : 老饼 发表日期 : 2023-02-05 04:26:27 更新日期 : 2025-06-25 00:23:11
本站原创文章,转载请说明来自《老饼讲解-机器学习》www.bbbdata.com



QR分解是一种矩阵分解方法,它把m×n的矩阵分解为列正交矩阵Q与上三角矩阵R的积

本节讲解QR分解的定义,并进一步讲述它的原理与分解流程,最后给出python实现QR分解的代码

通过本文,可以了解什么是QR分解,如何进行QR分解,以及它的分解原理和代码实现方法





     01. 什么是QR分解      




本节讲解QR分解是什么和如何进行QR分解




    QR分解的定义    


对于m×nm\times n列满秩矩阵A
可将其分解为列正交矩阵与上三角矩阵的积,
 
 Amn=QmnRnnA_{m*n} = Q_{m*n}*R_{n*n}
 其中
 
👉 Q列与列之间两两正交  
👉 R为非奇异上三角矩阵





    QR分解的方法    


符号说明
AA的第i列为αi\alpha_i  ,Q的第i列为βi\beta_i

A=[α1,α2,...,αn]A=[\alpha_1,\alpha_2,...,\alpha _n]
 Q=[β1,β2,...,βn]Q=[\beta_1,\beta_2,...,\beta _n] 

 Q的求法
 
Q中的β1β2β3...βn\beta_1、\beta_2、\beta_3、...、\beta_n求解如下
 

 β1=α1β2=α2α2β1β1β1β1β3=α3α3β1β1β1β1α3β2β2β2β2...βn=αnαnβ1β1β1β1αnβ2β2β2β2...αnβn1βn1βn1βn1\begin{aligned} \beta_1 =& \alpha_1 \\\beta_2 =& \alpha_2 - \dfrac{\alpha_2\cdot \beta_1}{\beta_1\cdot \beta_1} \beta_1 \\\beta_3 =& \alpha_3 - \dfrac{\alpha_3 \cdot \beta_1 }{\beta_1\cdot\beta_1}\beta_1 -\dfrac{\alpha_3 \cdot \beta_2 }{\beta_2\cdot\beta_2}\beta_2 \\... \\\beta_n=& \alpha_n - \dfrac{\alpha_n \cdot \beta_1 }{\beta_1\cdot\beta_1}\beta_1 -\dfrac{\alpha_n \cdot \beta_2 }{\beta_2\cdot\beta_2}\beta_2 -...-\dfrac{\alpha_n \cdot \beta_{n-1} }{\beta_{n-1}\cdot\beta_{n-1}}\beta_{n-1} \end{aligned}
 
 PASS
实际就是新的
α\alpha减去在各个之前已经正交化的β\beta上的投影分量,
从而得到由新 
α\alpha生成的新正交量β\beta

R的求法
在求得Q后,R可如下算得
R=[1α2β1β1β1αnβ1β1β101αnβ2β2β2001]R= \begin{bmatrix} 1 & \dfrac{\alpha_2 \cdot \beta_1}{\beta_1 \cdot \beta_1} & \dots & \dfrac{\alpha_n \cdot \beta_1}{\beta_1 \cdot \beta_1} \\ 0& 1 & \cdots& \dfrac{\alpha_n \cdot \beta_2}{\beta_2 \cdot \beta_2}\\ \vdots & \vdots& \ddots & \vdots\\ 0 & 0 & \cdots &1 \end{bmatrix}
✍️补充:关于Q的单位化
如果需要进一步单位化Q,
使Q的每一条列向量βi\beta_i 为单位向量,
则可以在以上结果进一步如下单位化Q
 
Q的单位化公式
 Q=[β1,β2,...,βn][1β10001β20000001βn]Q=[\beta_1,\beta_2,...,\beta_n]*\begin{bmatrix} \dfrac{1}{|\beta_1|} & 0 &\dots & 0\\ 0 & \dfrac{1}{| \beta_2|} &\dots & 0\\ 0 & 0 & \ddots & \vdots \\ 0 & 0 & 0&\dfrac{1}{|\beta_n|} \end{bmatrix}
 
 对应的R为
 

 R=[β1000β2000000βn]RR = \begin{bmatrix} |\beta_1| & 0 &\dots & 0\\ 0 &| \beta_2| &\dots & 0\\ 0 & 0 & \ddots & \vdots \\ 0 & 0 & 0&| \beta_n| \end{bmatrix} * R              








    02. QR分解的原理与公式推导    




本节讲解上节所提供的QR分解公式的原理推导





    相关引理    


 引理1:两条向量的正交
求b对a的正交向量d:
先求出b在a的投影向量c
c=(bcos<a,b>)aa=(babab)aa=(ab)a2a=(ab)aaac =(|b|*cos<a,b> )* \dfrac{a}{|a|}=\left (|b|*\dfrac{a\cdot b}{|a||b|} \right )*\dfrac{a}{|a|}=\dfrac{(a\cdot b)}{|a|^2}*a =\dfrac{(a\cdot b)}{ a\cdot a}*a

则b对a的正交向量 为
d=bc=b(ab)aaad=b-c=b-\dfrac{(a\cdot b)}{ a\cdot a}*a
引理2:与正交向量集的正交量
求 b对a=[a1,a2,...,an]a=[a_1,a_2,...,a_n]的正交量 d,使 d与a1,a2...ana_1,a_2...a_n都正交
其中, aa 是正交向量集,它的任意两条向量ai,aja_i,a_j 正交

将 b可以减去在a1,a2...ana_1,a_2...a_n上的投影分量分量即可
d=bk1a1k2a2....knand = b - k_1a_1-k_2a_2-....k_na_n
其中
kiaik_ia_i为b在aia_i上的投影分量

对于任意一个aia_i, d 都与aia_i正交,因为dd 中除了 bb和aia_i项,其余项都与aia_i正交,而b-k_ia_i又与a_i正交,所以d与a_i正交 
即:
aid=ai[(bkiai)(k1a1+k2a2+...ki1ai1+ki+1ai+1+....knan)]=ai[(bkiai)]ai[(k1a1+k2a2+...ki1ai1+ki+1ai+1+....knan)]=00=0\begin{aligned} a_i*d &= a_i* [(b-k_ia_i)-(k_1a_1+k_2a_2+...k_{i-1}a_{i-1}+k_{i+1}a_{i+1}+....k_na_n)] \\&= a_i* [(b-k_ia_i)]-a_i* [(k_1a_1+k_2a_2+...k_{i-1}a_{i-1}+k_{i+1}a_{i+1}+....k_na_n)] \\&= 0-0 = 0 \end{aligned}





     问题结论推导    


Q的求法
因此,要用A=[α1,α2,...,αn]A=[\alpha_1,\alpha_2,...,\alpha _n]生成正交向量 Q=[β1,β2,...,βn]Q=[\beta_1,\beta_2,...,\beta _n] , 
只需令β1=α1\beta_1=\alpha_1,然后不断的用新的α\alpha来生成与已有β\beta正交向量集的正交量即可
即可得到
β1=α1\beta_1 =\alpha_1                                                                                   
 β2=α2α2β1β1β1β1\beta_2 = \alpha_2 - \dfrac{\alpha_2\cdot \beta_1}{\beta_1\cdot \beta_1} \beta_1                                                                
 β3=α3α3β1β1β1β1α3β2β2β2β2\beta_3= \alpha_3 - \dfrac{\alpha_3 \cdot \beta_1 }{\beta_1\cdot\beta_1}\beta_1 -\dfrac{\alpha_3 \cdot \beta_2 }{\beta_2\cdot\beta_2}\beta_2                                           
....                                                                               
 βn=αnαnβ1β1β1β1αnβ2β2β2β2...αnβn1βn1βn1βn1\beta_n= \alpha_n - \dfrac{\alpha_n \cdot \beta_1 }{\beta_1\cdot\beta_1}\beta_1 -\dfrac{\alpha_n \cdot \beta_2 }{\beta_2\cdot\beta_2}\beta_2 -...-\dfrac{\alpha_n \cdot \beta_{n-1} }{\beta_{n-1}\cdot\beta_{n-1}}\beta_{n-1}
R的求法
又由上,可得
β1=α1\beta_1 =\alpha_1                                                                                     
 α2β1β1β1β1+β2=α2\dfrac{\alpha_2\cdot \beta_1}{\beta_1\cdot \beta_1} \beta_1+\beta_2 = \alpha_2                                                                  
 α3β1β1β1β1+α3β2β2β2β2+β3=α3\dfrac{\alpha_3 \cdot \beta_1 }{\beta_1\cdot\beta_1}\beta_1 +\dfrac{\alpha_3 \cdot \beta_2 }{\beta_2\cdot\beta_2}\beta_2 +\beta_3= \alpha_3                                              
...
 αnβ1β1β1β1+αnβ2β2β2β2+...+αnβn1βn1βn1βn1+βn=αn\dfrac{\alpha_n \cdot \beta_1 }{\beta_1\cdot\beta_1}\beta_1 +\dfrac{\alpha_n \cdot \beta_2 }{\beta_2\cdot\beta_2}\beta_2 +...+\dfrac{\alpha_n \cdot \beta_{n-1} }{\beta_{n-1}\cdot\beta_{n-1}}\beta_{n-1} +\beta_n= \alpha_n    

即:  

 [β1,β2,...,βn][1α2β1β1β1αnβ1β1β101αnβ2β2β2001]=[α1,α2,...,αn][\beta_1,\beta_2,...,\beta_n] \begin{bmatrix} 1 & \dfrac{\alpha_2 \cdot \beta_1}{\beta_1 \cdot \beta_1} & \dots & \dfrac{\alpha_n \cdot \beta_1}{\beta_1 \cdot \beta_1} \\ 0& 1 & \cdots& \dfrac{\alpha_n \cdot \beta_2}{\beta_2 \cdot \beta_2}\\ \vdots & \vdots& \ddots & \vdots\\ 0 & 0 & \cdots &1 \end{bmatrix} = [\alpha_1,\alpha_2,...,\alpha_n]

可知

 R=[1α2β1β1β1αnβ1β1β101αnβ2β2β2001]R= \begin{bmatrix} 1 & \dfrac{\alpha_2 \cdot \beta_1}{\beta_1 \cdot \beta_1} & \dots & \dfrac{\alpha_n \cdot \beta_1}{\beta_1 \cdot \beta_1} \\ 0& 1 & \cdots& \dfrac{\alpha_n \cdot \beta_2}{\beta_2 \cdot \beta_2}\\ \vdots & \vdots& \ddots & \vdots\\ 0 & 0 & \cdots &1 \end{bmatrix}
3、单位化QR
QR=[β1,β2,...,βn][1β10001β20000001βn][β1000β2000000βn]RQR=[\beta_1,\beta_2,...,\beta_n]*\begin{bmatrix} \dfrac{1}{|\beta_1|} & 0 &\dots & 0\\ 0 & \dfrac{1}{| \beta_2|} &\dots & 0\\ 0 & 0 & \ddots & \vdots \\ 0 & 0 & 0&\dfrac{1}{|\beta_n|} \end{bmatrix}*\begin{bmatrix} |\beta_1| & 0 &\dots & 0\\ 0 &| \beta_2| &\dots & 0\\ 0 & 0 & \ddots & \vdots \\ 0 & 0 & 0&| \beta_n| \end{bmatrix} * R
可得,单位化的Q为:
Q=[β1,β2,...,βn][1β10001β20000001βn]Q=[\beta_1,\beta_2,...,\beta_n]*\begin{bmatrix} \dfrac{1}{|\beta_1|} & 0 &\dots & 0\\ 0 & \dfrac{1}{| \beta_2|} &\dots & 0\\ 0 & 0 & \ddots & \vdots \\ 0 & 0 & 0&\dfrac{1}{|\beta_n|} \end{bmatrix}
对应的R为:
R=[β1000β2000000βn]RR = \begin{bmatrix} |\beta_1| & 0 &\dots & 0\\ 0 &| \beta_2| &\dots & 0\\ 0 & 0 & \ddots & \vdots \\ 0 & 0 & 0&| \beta_n| \end{bmatrix} * R







     03. QR分解的代码实现    





本节展示如何用python代码实现QR分解,并展示运行结果





     python实现矩阵的QR分解     


使用python实现QR分解只需按上述理论步骤进行即可
具体代码实现如下:
"""
代码说明:本代码用于实现-矩阵的QR分解
本代码来自《老饼讲解-机器学习》www.bbbdata.com
"""
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---")                       # 打印矩阵A
    print('A=\n',A)                                             # 打印矩阵Q
    print('Q=\n',Q)                                             # 打印矩阵
    print('R=\n',R)                                            
    print("\n-----用QTQ验证Q各个向量正交--")                    # 验证Q各个向量正交 
    QTQ= Q.T@Q                                                  # 计算QTQ
    QTQ[np.abs(QTQ)<0.00001]=0                                  # 将过小的元素置0
    print('Q^T*Q=\n',QTQ)                                       # 打印QTQ
    print("\n-----用A-QR验证A=QR--")                            # 验证A=QR
    err = A-Q@R                                                 # 计算A与QR的误差
    err[np.abs(err)<0.00001]=0                                  # 将过小的元素置0
    print('A-QR=\n',err)                                        # 打印误差

    #----------单位化------------------------------------
    print("\n----单位化后的QR:nQ,nR---")                        # 打印单位化后的QR
    print('nQ=\n',nQ)                                           # 打印单位化的Q
    print('nR=\n',nR)                                           # 打印单位化的R
    print("\n-----用nQ^T*nQ验证nQ各个向量正交--")               # 验证单位化的结果
    nQTnQ= nQ.T@nQ                                              # 计算nQTnQ
    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)                                    # 
代码运行结果如下:
 
代码运行结果-分解后的QR 
 代码运行结果-单位化后的QR







好了,以上就是矩阵的QR分解以及代码实现了~









 End 






图标 评论
添加评论