长矩阵（宽矩阵） 基于Householder变换的QR分解 Python源码

---

前言

例题：随机产生一个11 * 8的实质矩阵A,并且分别实现基于Givens变换和Householder变换的QR分解
网上关于QR分解的讲解非常充分，但大多是对方阵进行的。本次例题作用于一个11 * 8的长矩阵的QR分解，查阅网上的代码发现Householder变换的代码在方阵变为长、宽矩阵时无法运行，故按照计算方式写一种计算方法。

一、基于Householder变换的QR分解 原理分析

任何的实方阵都可以做QR分解：其中 Q 是正交阵， R 是上三角矩阵.。

参考图片与原理详见

https://zhuanlan.zhihu.com/p/112327923
https://blog.youkuaiyun.com/m0_37604894/article/details/123360420

根据以上算法，在对实值矩阵A进行Householder变换时，按列依次计算

• 每次需要对矩阵一列进行反射变换，得到上三角矩阵的一列。
• 每次处理完一列，将处理好R矩阵第一行第一列去掉（可以使用矩阵切片操作），其余继续进行QR分解。
• 当R矩阵只剩1行或1列时，中止。
• 将每列的正交变换矩阵合在一起，得到最终的变换矩阵Q，R=Q逆乘A。

注：当实值矩阵不是方阵时（m * n），正交变换矩阵为方阵（m * m），左乘代表行变换。

二、Python源码

代码如下（示例）：

# 2.householder进行QR分解，并输出QR
import numpy as np
def qr_householder(A):
    h, l = A.shape  # h,l分别为输入矩阵的行数和列数
    q = np.identity(h)  # q矩阵h阶单位矩阵
    R = np.copy(A)  # r矩阵为A的副本
    QQQ = np.eye(h)  # 分解后的正交矩阵Q 方阵 大小为原矩阵行数
    for i in range(h-1):
        # 取矩阵的第一列 注：每次处理完一列，将处理好R矩阵第一行第一列去掉
        x = R[:, 0]
        e = np.zeros_like(x)
        e[0] = 1  # 方向矩阵
        x_ = np.linalg.norm(x)  # x模长
        v = (x-x_*e) / np.linalg.norm(x - x_*e)
        H = np.eye(len(x)) - 2 * np.outer(v, v)  # 反射变换矩阵
        R_ = np.dot(H, R)
        # 还原维数 迭代求得分解后的正交矩阵QQQ
        QQQ_new = np.eye(h)  # 单位矩阵
        QQQ_new[i:,i:] = H.T
        QQQ = np.dot(QQQ, QQQ_new)

        # 注：每次处理完一列，将处理好R矩阵第一行第一列去掉
        R = R_[1:, 1:]
        # 当处理到R矩阵只剩1行或1列时中止
        if R_.shape[1] == 1 or R_.shape[1] == 1:
            R = np.dot(QQQ.T, A)
            print('R', R, R.shape)
            break
    return QQQ, R


A = np.random.randint(-10, 10, [11, 8])  # 形成元素取值在-10到10之间的随机矩阵A
np.set_printoptions(suppress=True, precision=2)  # 使用浮点数显示，设置小数位置2位
Q, R = qr_householder(A)  # q和r等于QR分解函数的返回值
print('A', A, A.shape)
print('Q', Q, Q.shape)
print('R', R, R.shape)

三、结果验证

生产单位元素取值在-10到10之间的随机矩阵A（11 * 8），得到分解后的Q矩阵和R矩阵如下图所示：