householder进行矩阵QR分解

1. HouseHolder 理论分析

HouseHolder变换的目的是为了将向量$x_1$通过镜面反射的方式，在保持大小不变的情况下，变成在x轴上的向量$x_2$具体如下图所示：

HouseHolder变换和Givens变换的目的是一致的，都是通过反射，旋转的方式，将任意向量$x_1$变到x轴上，形成向量$x_2$

• 由图可得向量之间关系如下
  $$x_1=x_2+2w\text{\hspace{1em}(1)}$$
• 做一个辅助单位向量$w_0$，使得$w_0$平行于w,$|w_0|=1$
  
• $w$的长度可以计算如下：
  $$w_0^T{x}_1=|w_0||x_1|\cos {\theta }_2\text{\hspace{1em}(2)}$$
  $$|w|=|x_1|\cos {\theta }_2\text{\hspace{1em}(3)}$$
• 方程1,3结合可得如下：
  $$x_1=x_2+2\frac{w_0^T{x}_1}{|w_0|}{w}_0\text{\hspace{1em}(4)}$$
• 因为我们定义$|w_0|=1$,并且$w_0^T{x}_1$为一个数，所以在矩阵中可以任意调位置，即
  $$x_1=x_2+2w_0{w}_0^T{x}_1\text{\hspace{1em}(5)}$$
• 整理可得如下
  $$x_2=x_1-2w_0{w}_0^T{x}_1=(I-2w_0{w}_0^T)x_1\text{\hspace{1em}(6)}$$
• 变成常规大家熟悉的公式如下,u为单位向量
  $$u=\left[\begin{array}{c}1\\ \\ 0\\ \\ 0\end{array}\right]\text{\hspace{1em}(7)}$$
  $$y=(I-2u{u}^T)x\text{\hspace{1em}(8)}$$

2. Householder 案例

2.1 Householder进行QR分解

$$A=QR\text{\hspace{1em}(9)}$$
$$A=\left[\begin{array}{cccc}1& 2& 0& 1\\ & & & \\ 1& 0& 3& 1\\ & & & \\ 1& 0& 3& 2\\ & & & \\ 1& 2& 0& 2\end{array}\right]\text{\hspace{1em}(10)}$$
从矩阵A中取出第一列 $\overrightarrow{a_1}=[1,1,1,1{]}^T$ ，$\Vert {\alpha }_1{\Vert }_2=\sqrt{1^2+1^2+1^2+1^2}=2$
$$\overrightarrow{u_1}=\frac{\overrightarrow{a_1}-{\alpha }_1\overrightarrow{e_1}}{\Vert \overrightarrow{a_1}-{\alpha }_1\overrightarrow{e_1}\Vert }\text{\hspace{1em}(11)}$$
$$\overrightarrow{e_1}=[1,0,0,0{]}^T\text{\hspace{1em}(12)}$$

• 带入可得：
  $$\overrightarrow{u_1}=\frac{1}{2}[-1,1,1,1{]}^T\text{\hspace{1em}(13)}$$
• 根据householder 来定义矩阵$H_1$
  $$\overrightarrow{H_1}=\vec{I}-2\overrightarrow{u_1}(\overrightarrow{u_1}{)}^T\text{\hspace{1em}(14)}$$
• 带入可得：
  $$\overrightarrow{H_1}=\frac{1}{2}\left[\begin{array}{cccc}1& 1& 1& 1\\ & & & \\ 1& 1& -1& -1\\ & & & \\ 1& -1& 1& -1\\ & & & \\ 1& -1& -1& 1\end{array}\right]\text{\hspace{1em}(15)}$$
• 第一次householder变换 ，目的是保证第一列只有存在第一个元素，其他位置为0：
  $$\overrightarrow{H_1}A=\frac{1}{2}\left[\begin{array}{cccc}1& 1& 1& 1\\ & & & \\ 1& 1& -1& -1\\ & & & \\ 1& -1& 1& -1\\ & & & \\ 1& -1& -1& 1\end{array}\right]\left[\begin{array}{cccc}1& 2& 0& 1\\ & & & \\ 1& 0& 3& 1\\ & & & \\ 1& 0& 3& 2\\ & & & \\ 1& 2& 0& 2\end{array}\right]\text{\hspace{1em}(16)}$$
  $$\overrightarrow{H_1}A=\left[\begin{array}{cccc}2& 2& 3& 3\\ & & & \\ 0& 0& 0& -1\\ & & & \\ 0& 0& 0& 0\\ & & & \\ 0& 2& -3& 0\end{array}\right]\text{\hspace{1em}(17)}$$
• 我们重新第一新的矩阵$A_2$
  $$\overrightarrow{A_2}=\left[\begin{array}{ccc}0& 0& -1\\ & & \\ 0& 0& 0\\ & & \\ 2& -3& 0\end{array}\right]\text{\hspace{1em}(18)}$$
  从矩阵$A_2$中取出第一列 $\overrightarrow{a_2}=[0,0,2{]}^T$ ，$\Vert {\alpha }_2{\Vert }_2=\sqrt{0^2+0^2+2^2}=2$
  $$\overrightarrow{u_2}=\frac{\overrightarrow{a_2}-{\alpha }_2\overrightarrow{e_1}}{\Vert \overrightarrow{a_2}-{\alpha }_2\overrightarrow{e_1}\Vert }\text{\hspace{1em}(19)}$$
  $$\overrightarrow{e_1}={\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^T\text{\hspace{1em}(20)}$$
• 带入相关参数可得:
  $$\overrightarrow{u_2}=\frac{\sqrt{2}}{2}{\left[\begin{array}{ccc}-1& 0& 1\end{array}\right]}^T\text{\hspace{1em}(21)}$$
  $$\widetilde{H_2}=\vec{I}-2\overrightarrow{u_2}(\overrightarrow{u_2}{)}^T\text{\hspace{1em}(22)}$$
• 代入可得如下：
  $$\widetilde{H_2}=\left[\begin{array}{ccc}0& 0& 1\\ & & \\ 0& 1& 0\\ & & \\ 1& 0& 0\end{array}\right]\text{\hspace{1em}(23)}$$
• 同理可得$\widetilde{H_2}\overrightarrow{A_2}$
  $$\widetilde{H_2}\overrightarrow{A_2}=\left[\begin{array}{ccc}2& -3& 0\\ & & \\ 0& 0& 0\\ & & \\ 0& 0& -1\end{array}\right]\text{\hspace{1em}(24)}$$
  $$\overrightarrow{H_2}=\left[\begin{array}{cc}1& \\ & \\ & \widetilde{H_2}\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ & & & \\ 0& 0& 0& 1\\ & & & \\ 0& 0& 1& 0\\ & & & \\ 0& 1& 0& 0\end{array}\right]\text{\hspace{1em}(25)}$$
• 依公式，我们定义R如下：
  $$\vec{R}=\overrightarrow{H_2}(\overrightarrow{H_1}\vec{A})=\left[\begin{array}{cccc}2& 2& 3& 3\\ & & & \\ 0& 2& -3& 0\\ & & & \\ 0& 0& 0& 0\\ & & & \\ 0& 0& 0& -1\end{array}\right]\text{\hspace{1em}(26)}$$
• 依公式，我们定义Q如下：
  $$\vec{Q}=\overrightarrow{H_1}\overrightarrow{H_2}=\frac{1}{2}\left[\begin{array}{cccc}1& 1& -1& -1\\ & & & \\ 1& -1& -1& 1\\ & & & \\ 1& -1& 1& -1\\ & & & \\ 1& 1& 1& 1\end{array}\right]\text{\hspace{1em}(27)}$$

3. 结论：

• householder的意义是通过househoulder变换，可以将矩阵A 变成一个正交矩阵Q和上三角矩阵R
  $$A=QR\text{\hspace{1em}(28)}$$
  $$A=\left[\begin{array}{cccc}1& 2& 0& 1\\ & & & \\ 1& 0& 3& 1\\ & & & \\ 1& 0& 3& 2\\ & & & \\ 1& 2& 0& 2\end{array}\right]\text{\hspace{1em}(29)}$$
  $$\vec{Q}=\frac{1}{2}\left[\begin{array}{cccc}1& 1& -1& -1\\ & & & \\ 1& -1& -1& 1\\ & & & \\ 1& -1& 1& -1\\ & & & \\ 1& 1& 1& 1\end{array}\right]\text{\hspace{1em}(30)}$$
  $$\vec{R}=\left[\begin{array}{cccc}2& 2& 3& 3\\ & & & \\ 0& 2& -3& 0\\ & & & \\ 0& 0& 0& 0\\ & & & \\ 0& 0& 0& -1\end{array}\right]\text{\hspace{1em}(31)}$$

4. Python code

• 代码：

import numpy as np
np.set_printoptions(suppress=True,precision=3)

def householder_reflection(a):
    """
    计算给定向量a的Householder反射矩阵。
    """
    e = np.zeros_like(a)
    e[0] = np.linalg.norm(a)

    v = a - e
    v_norm = np.linalg.norm(v)

    # 避免除以零的情况
    if v_norm == 0:
        return np.eye(len(a))

    v = v / v_norm

    H = np.eye(len(a)) - 2 * np.outer(v, v)

    return H


def qr_decomposition(A):
    """
    使用Householder反射进行QR分解。
    参数:
        A: 输入矩阵 (m x n)
    返回:
        Q: 正交矩阵 (m x m)
        R: 上三角矩阵 (m x n)
    """
    m, n = A.shape
    Q = np.eye(m)
    R = A.copy()

    for i in range(n):
        # 从R的第i列开始，获取子向量
        x = R[i:, i]

        # 计算Householder反射矩阵
        H_i = np.eye(m)
        H = householder_reflection(x)
        H_i[i:, i:] = H

        # 更新R和Q
        R = H_i @ R
        Q = Q @ H_i

    return Q, R




A = np.array([[1,2,0,1],
              [1,0,3,1],
              [1,0,3,2],
              [1,2,0,2]])
Q, R = qr_decomposition(A)

print("Q:")
print(Q)
print("\nR:")
print(R)

# 验证 A = Q * R
print("\nQR Product (should equal A):")
print(Q @ R)

• 结果：

Q:
[[ 0.5  0.5 -0.5 -0.5]
 [ 0.5 -0.5 -0.5  0.5]
 [ 0.5 -0.5  0.5 -0.5]
 [ 0.5  0.5  0.5  0.5]]

R:
[[ 2.  2.  3.  3.]
 [ 0.  2. -3. -0.]
 [ 0. -0.  0.  1.]
 [ 0. -0.  0.  0.]]

QR Product (should equal A):
[[1. 2. 0. 1.]
 [1. 0. 3. 1.]
 [1. 0. 3. 2.]
 [1. 2. 0. 2.]]

5. 创建househoulder矩阵

• 通过给定一个任意一个单位向量u，可以创建一个单位正交矩阵，为Househoulder矩阵，具体如下：

#!/usr/bin/env python
# -*- coding:utf-8 -*-
# @FileName  :househoulder_vector.py
# @Time      :2024/10/18 20:26
# @Author    :Jason Zhang
import numpy as np

np.random.seed(2024)
np.set_printoptions(suppress=True, precision=3)

if __name__ == "__main__":
    run_code = 0
    N = 4
    for i in range(N):
        u_i = np.random.randn(N)
        u_i = u_i / (np.linalg.norm(u_i, ord=2))
        H_i = np.eye(N) - 2 * np.outer(u_i, u_i)
        result = H_i @ H_i.T
        print(f"\n")
        print("*" * 20)
        print(f"u_i={u_i}")
        print(f"H_{i}=\n{H_i}")
        print(f"HiNorm={round(np.linalg.norm(H_i, ord=2))}")
        print(f"result=\n{result}")
        print("*" * 20)

********************
u_i=[ 0.906  0.401 -0.109 -0.082]
H_0=
[[-0.642 -0.726  0.198  0.149]
 [-0.726  0.679  0.088  0.066]
 [ 0.198  0.088  0.976 -0.018]
 [ 0.149  0.066 -0.018  0.987]]
HiNorm=1
result=
[[ 1. -0.  0.  0.]
 [-0.  1.  0.  0.]
 [ 0.  0.  1. -0.]
 [ 0.  0. -0.  1.]]
********************


********************
u_i=[ 0.279  0.353 -0.797 -0.403]
H_1=
[[ 0.845 -0.197  0.444  0.225]
 [-0.197  0.751  0.563  0.285]
 [ 0.444  0.563 -0.27  -0.643]
 [ 0.225  0.285 -0.643  0.675]]
HiNorm=1
result=
[[ 1. -0.  0.  0.]
 [-0.  1.  0.  0.]
 [ 0.  0.  1. -0.]
 [ 0.  0. -0.  1.]]
********************


********************
u_i=[0.231 0.051 0.529 0.815]
H_2=
[[ 0.893 -0.024 -0.244 -0.376]
 [-0.024  0.995 -0.054 -0.084]
 [-0.244 -0.054  0.441 -0.862]
 [-0.376 -0.084 -0.862 -0.329]]
HiNorm=1
result=
[[ 1. -0. -0. -0.]
 [-0.  1. -0. -0.]
 [-0. -0.  1. -0.]
 [-0. -0. -0.  1.]]
********************


********************
u_i=[-0.708 -0.131  0.563  0.406]
H_3=
[[-0.001 -0.185  0.797  0.575]
 [-0.185  0.966  0.148  0.106]
 [ 0.797  0.148  0.366 -0.458]
 [ 0.575  0.106 -0.458  0.67 ]]
HiNorm=1
result=
[[ 1. -0.  0.  0.]
 [-0.  1.  0.  0.]
 [ 0.  0.  1. -0.]
 [ 0.  0. -0.  1.]]
********************