QR分解

实数域

Hermite矩阵

$满足{\mathrm{P}}^{\mathrm{H}}=\mathrm{P}$
$P=(I-2v{v}^T)为初等正交Hermite（对称）矩阵$
$即：P^T=P,P^{T}P=I$

QR分解

Q为正交矩阵 R为上三角阵

一步

$先考虑Px=ke,P=(I-2v{v}^T),x\in {\mathrm{R}}^d的情况下，求解P。$
$Px=ke$
$k^2=\Vert Px{\Vert }_2^2=x^2=x_1^2+x_2^2+\dots +x_d^2$
$S:=\sqrt{x_1^2+x_2^2+\dots +x_d^2}$
$then,k=\pm S$
$K:=v^{T}x$
$ x-2Kv=ke$
$then, x_1-2Kv_1=k $
$x_i-2K{v}_i=0,i̸=1$
$\therefore v_1=(x_1-k)/2K$
$v_i=x_i/2K$
$\because Px=x-2Kv=ke,2Kv=x-ke$
$\therefore 2K^2=S^2-x_{1}k=S^2\pm x_{1}S$
$ifu:=(x_1\pm S,x_2,x_3,\dots ,x_d)$
$then,P=(I-u{u}^T/2K^2)$

正负号的选择

$k=\pm S$
$2K^2=S^2\pm x_{1}S$
$如何选择正负号呢，为了避免K=0,我们选择符号为x_1的符号即可。$

多步

$假设矩阵A的前r列已经是上三角形，则A_r如下$

$为了把对角化剩下的，我们只需取：$

二者相乘，我们容易发现，右下角变成了第一步的情况，依此知道对角化完全。
可得:
$A=QR$
$Q=P_1{P}_2\dots$

Tips 子空间

$V=QR$
$Q=(q_1,q_2,\dots ,q_r,\dots ,q_d)$
$V=(v_1,v_2,\dots ,v_r)$
$如果v_1,v_2,\dots ,v_r的极大线性无关组为本身，那么$
$q_1,q_2,\dots ,q_r与v_1,v_2,\dots ,v_r张成同一个子空间（互相线性表出）。$

代码

import numpy as np

def step_one_QR(x):
    """
    一步
    x 为ndarray
    :param x:向量
    :return: P
    """
    u = np.array(x, dtype=float)
    sign = lambda x: 1 if x >= 0 else -1
    S_2 = u @ u
    S = np.sqrt(S_2) * sign(u[0])
    K_2 = S_2 + S * u[0] #K_2 = 2K^2 here
    vector_ones = np.ones(len(u), dtype=float)
    u[0] += S
    return np.diag(vector_ones) - np.outer(u, u) / K_2


def step_all_QR(A):

    """
    d:A的行数
    r:A的列数
    P:每步得到的Hermite矩阵
    :param A:d x r的矩阵
    :return:Q, R(A)
    """
    A = np.array(A, dtype=float) #注意经此操作后，后续操作不会改变原来的A
    d, r = A.shape
    Q = np.diag(np.ones(d, dtype=float))
    for i in range(r):
        P = np.diag(np.ones(d, dtype=float))
        P[i:, i:] = step_one_QR(A[i:,i])
        Q = Q @ P
        A = P @ A

    return Q, A

"""
对上面的一个改进，速度明显提高，不过
俩个方法对列非满秩的抗性都不高。
"""
def step_one_QR_2(x):
    """
    一步 第二种
    x 为ndarray
    :param x:向量
    :return:u, K_2
    """
    u = np.array(x, dtype=float)
    sign = lambda x: 1 if x >= 0 else -1
    S_2 = u @ u
    S = np.sqrt(S_2) * sign(u[0])
    K_2 = S_2 + S * u[0] #K_2 = 2K^2 here
    vector_ones = np.ones(len(u), dtype=float)
    u[0] += S
    return u, K_2

def step_all_QR_2(A):
    """
     d:A的行数
     r:A的列数
     :param A:
     :return:
     """
    A = np.array(A, dtype=float)  # 注意经此操作后，后续操作不会改变原来的A
    d, r = A.shape
    Q = np.diag(np.ones(d, dtype=float))
    for i in range(r):
        u, k = step_one_QR_2(A[i:, i])
        if k <= 0.0000001:
            print("Matrix may not full rank...")
            return Q, A  #A非列满秩
        Q[:, i:] -= np.outer(Q[:, i:] @ u, u / k)
        A[i:, :] -= np.outer(u / k, u @ A[i:, :])

    return Q, A