基于Householder变换实现超定矩阵的QR分解的算法

实现矩阵的QR分解有多种算法，除了利用Gram-Schmidt正交化之外，另一种经典的方法是基于Householder变换的。Householder 变换在算法优化中扮演着重要的角色，其定义如下：

Householder 变换有诸多性质，如对称性、正交性等。事实上，存在 Householder 变换使得任一非零向量的第一个分量变为原向量的二范数或二范数的相反数，其余分量变为 0，这条性质在实现 QR 分解的算法中发挥着至关重要的作用。

基于这条性质，可以借助Householder变换获取上三角矩阵，在第i次迭代时，选取第i个对角元以及其正下方的元素作为第i个待处理的向量，基于此向量，获取相应的Householder向量并对原始矩阵的子矩阵做相应的Householder变换，可以使得原始矩阵的第i个对角元正下方的元素全部变为0，经过反复迭代，就可以得到上三角矩阵。

正交矩阵的获取则是基于Householder变换的正交性，设置初始矩阵为单位矩阵(阶数为m,即原始矩阵的行数)，每次迭代的时候都右乘以新构造的Householder矩阵，经反复迭代，最终可得正交矩阵。

对于一个超定矩阵(m>n，即矩阵行数大于列数)，基于Householder变换实现QR分解的思想如下：

基于上述过程实现QR分解的python代码参考如下：

import numpy as np
#获取未归一化的Householder向量，相应的系数和调节符号的变量sign
def house(x):
    N=len(x)#向量x的维数
    u=np.array(x)#没有归一化的Householder向量，后仅需修改第一个分量
    v0=0
    #sign可以将对角元进行调整，若x[0]>0，初始对角元算出为负值，可乘以sign调整为正值
    if x[0]<0:
        sign=1
    elif x[0]>=0:
        sign=-1
    for i in range(1,N):
        v0=v0+x[i]**2
    if v0==0:#原始向量后面的分量全部为0的特殊情况
        corr=0#系数归0
    else:
        alpha=np.sqrt(v0+x[0]**2)
        if x[0]<0:
            u[0]=u[0]-alpha #更换第一个分量
        elif x[0]>=0:
            u[0]=u[0]+alpha
        corr=2/(u[0]**2+v0)#确定系数
    return u,corr,sign
#Householder变换作用于给定的矩阵A，得出最后的上三角矩阵
def House_tran(A,vec,beta,t,sign):
    m=A.shape[0]#矩阵A的行数，也是vec的维数
    n=A.shape[1]
    Q=A[t:m,t:n]#第(t+1)次迭代，仅仅取第(t+1)个对角元到最后一行和最后一列构造的子矩阵
    Q_row=Q.shape[0]
    Q_col=Q.shape[1]
    sub_new_A=np.zeros((Q_row,Q_col))
    for i in range(0,Q_row):
        for j in range(0,Q_col):
            sum=0
            for k in range(0,Q_row):
                sum+=vec[i]*vec[k]*Q[k][j]
            sub_new_A[i][j]=Q[i][j]-beta*sum
    #直接在上一次循环的第二层赋值会改变Q矩阵(Q来自于A,A的变动会更新Q)
    for i in range(0,Q_row):
        for j in range(0,Q_col):
            A[t+i][t+j]=sign*sub_new_A[i][j]
    return A
#计算正交矩阵的程序
def ortho_comput(A,vec,beta,t,sign):
    m=A.shape[0]
    n=A.shape[1]
    Q=A[:,t:n]
    Q_row=Q.shape[0]
    Q_col=Q.shape[1]
    sub_new_B=np.zeros((Q_row,Q_col))
    for i in range(Q_row):
        for j in range(Q_col):
            temp=0
            for k in range(Q_col):
                temp+=Q[i][k]*vec[k]*vec[j]
            sub_new_B[i][j]=Q[i][j]-beta*temp
    for i in range(Q_row):
        for j in range(Q_col):
            A[i][t+j]=sign*sub_new_B[i][j]
    return A
#基于Householder变换实现QR分解的算法
def House_QR(A):
    n=A.shape[1]#矩阵A的列数
    m=A.shape[0]#矩阵A的行数
    Q_0=np.eye(m)#构造m阶的单位矩阵，参与后续正交矩阵的迭代
    R_0=A#利用Householder变换反复迭代，最终可得上三角矩阵
    if m<n:
        print("行数小于列数，为欠定矩阵，无法进行QR分解！")
    elif m>n:
        for i in range(0,n):
            #若行数等于列数，迭代到第n-1列，或行数大于列数，迭代到最后一列
            A_vec=R_0[i:,i]#第i+1个对角元正下方的元素组成的向量
            v_i,b,sign_1=house(A_vec)#获取相应的Householder向量和系数
            R_0=House_tran(R_0,v_i,b,i,sign_1)
            Q_0=ortho_comput(Q_0,v_i,b,i,sign_1)
            Q = Q_0[0:m, 0:n]
            R = R_0[0:n, 0:n]
    return Q,R
M=np.random.rand(7,5)
M_copy = np.copy(M)#保持原始矩阵的初值
Q,R=House_QR(M)
print("原始矩阵：\n",M_copy)
print("列正交矩阵：\n",Q)
print("上三角矩阵：\n",R)
#程序计算得到的结果
M_val=np.dot(Q,R)
U=np.transpose(Q)
I_val=np.dot(U,Q)
#验证结果
def valid(A,B):
    m=A.shape[0]
    n=A.shape[1]
    k=0
    for i in range(m):
        for j in range(n):
            if abs(A[i][j]-B[i][j])>1e-15:
                print("两矩阵不相等")
                return
    print("两矩阵相等")
    return
print("验证QR分解计算的结果:")
valid(M_copy,M_val)
print("验证Q是否为列正交矩阵：")
valid(I_val,np.eye(5))

经过验证，基于此算法算出的正交矩阵和上三角矩阵的通常积，其每个元素和原始矩阵的相应元素之间的差值不超过1e-15,这个精度是比较高的。