【介绍】矩阵分解之- QR分解

QR分解也是矩阵中的一种常用分解，它可能将矩阵A分解为一个正交矩阵与上三角矩阵的积，即A=Q*R，这不仅有利于矩阵的计算，在一些理论推导中也常扮演着重要的角色，下面我们就一起来看看QR分解是什么，以及QR分解的方法、具体代码实现例子。

一、QR分解

QR分解是指：对于列满秩矩阵A，可以将其分解为列正交矩阵与上三角矩阵的积，即：

其中

1. Q列与列之间两两正交

2. R为非奇异上三角矩阵

更形象地，QR分解示图如下：

如图所示，QR分解将A分解为Q和R，Q列列正交，而R则是上三角矩阵。

二、QR分解的方法

好了，我们都知道QR分解是怎么一回事了，那怎么才能将矩阵QR分解呢，这里我们先讲QR分解的公式，下面再慢慢来说公式是怎么来的。

2.1. 符号说明

记的第列为  ,Q的第列为，那么，A、Q分别可以记为：



 

2.2. Q的求法

Q中的求解如下

初看公式可能有点摸不着头脑，其实仔细观察一下就会发现，它只不过是每次用新的减去在各个之前已经正交化的上的投影分量，从而得到由新 生成的新正交量。

2.3. R的求法

在求得Q后，R可如下算得

2.4. 关于Q的单位化

上面求解得到的Q只是列列正交，虽然也满足QR分解的要求，但有时我们希望Q是单位正交矩阵，那么可以进一步将Q单位化，使Q的每一条列向量 为单位向量。

好了，下面我们来看看怎么单位化Q，以及Q正交化后，R对应所做的修改，闲话少说，直接上公式。

Q的单位化公式如下：

对应的R为

就是这么简单，对着公式来单位化就可以了。

三、QR分解的原理与公式推导

好了，下面我们再来说说QR分解所用到的公式是怎么推导出来的。

3.1. QR分解-公式推导-相关引理

在说推导之前，我们先来说一下相关引理。

引理1：两条向量的正交向量

假设现在有两条向量，现在要求对的正交向量d：

先求出b在a的投影向量c

则对的正交向量为

引理2：与正交向量集的正交量

求 b对的正交量 d，使 d与都正交，其中,  是正交向量集，它的任意两条向量 正交。

只需要将 b减去在上的投影分量即可

其中，为b在上的投影分量

对于任意一个,  都与正交，因为 中除了 和项，其余项都与正交，而又与正交，所以d与正交，即：

3.2. QR分解-Q的公式推导

我们先来推导QR分解中Q的公式。

要用生成正交向量  ，

只需令，然后不断的用新的来生成与已有正交向量集的正交量即可。

即可得到

3.3. QR分解-R的公式推导

接下来推导QR分解中R的公式。

由上，可得

即：  

可知

3.3. QR分解-QR单位化-公式推导

由于：

可得，单位化的Q为：

对应的R为：

四、QR分解-代码实现

好了，下面我们根据之前所说的QR分解公式，用代码来实现QR分解。

具体python代码实现如下：

"""
代码说明：本代码用于实现-矩阵的QR分解
本代码来自《老饼讲解-机器学习》www.bbblearn.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)                                   # 打印单位化后的QTQ
    print("-----用A-QR验证A=QR--")                              # 打印
    nerr = A-nQ@nR                                              # 计算A与还原后的误差
    nerr[np.abs(nerr)<0.00001]=0                                # 屏蔽掉一些太小的值
    print('A-nQ*nR=\n',nerr)                                    # 打印A与还原后的误差

代码运行结果如下：

可以看到，按照上述公式，就能对A进行QR分解，并验证Q*R能还原出A。

总结

总的来说，QR分解就是将矩阵A分解成一个正交矩阵和上三角矩阵的积，而分解的过程也是比较简单的，直接使用公式来进行分解就可以了，QR分解可以简化计算，同时也可以帮助一些矩阵问题进行理论推导，熟悉它有时可以化腐朽为神奇。