1 Householder reflector

  Householder反射是这样子的(图片来自瑞典皇家理工学院)：

  图中u是长度为1的向量。x是任意向量，H是u的Householder reflector。可见无论x是什么向量，$Hx$始终除于和u正交的平面上。H和u的关系是：
$$H=I-2u{u}^{\ast }$$
  $u^{\ast }$是u的共轭转置，$u$是单位向量。在有些书上，也写成$u^H$.也就是说H这个矩阵是由$u$确定的。那么接下来的问题是如何精确控制变换方向了。如果能精确控制，把向量投影到标准基的方向上，那么QR分解就容易了。
  那么这和QR分解有什么关系呢？

2 分解步骤

  首先把矩阵A按列向量分块，为$(a_1,a_2,\cdots ,a_n)$,设单位矩阵的第一列为$e_1$,取$u$向量为以下向量：
$$u=\frac{a_1-\parallel a_1\parallel e_1}{\parallel (a_1-\parallel a_1\parallel e_1)\parallel }$$
  用这个$u$向量构造的Householder矩阵$H$会把$a_1$投影到$e_1$的方向上。也就是：
$$H{a}_1=\parallel a_1\parallel e_1$$
  那么$HA$就是这个样子：
$$HA=H(a_1,a_2,\cdots ,a_n)=\left(\begin{array}{ccc}\parallel a_1\parallel & \cdots & \ast \\ 0& \cdots & \ast \\ ⋮& \ddots & \ast \\ 0& \cdots & \ast \end{array}\right)$$
  把右边的矩阵叫做B，那么有：
$$\begin{gathered} HA=B \\ {H}^{-1}HA=H^{-1}B \\ A=H^{-1}B \end{gathered}$$
  Householder变换矩阵是自逆矩阵，也是对合矩阵，所以$H=H^{-1}$.所以又有：
$$A=HB$$
  $B$矩阵的右下角，再当成新的子矩阵，继续完成这个过程，这样得到一系列的H矩阵，编号为$H_1,H_2,\cdots ,H_n$.这写矩阵乘以$A$,得到的结果也是不断把对角线以下的元素变成$0$,所以最终结果就是上三角矩阵$R$。但是会发现这些Householder矩阵是不同阶的。为了同阶，可以在对角线上补1，凑成$n$阶矩阵，这种左上角补1的矩阵乘以任何矩阵不会改变左上角的元素。如下面这样处理：
$${\tilde{H}}_i={\left(\begin{array}{cccc}1& & & \\ & 1& & \\ & & \ddots & \\ & & & H_i\end{array}\right)}_n$$
  所以整个过程就是：
$$\begin{gathered} {\tilde{H}}_n{\tilde{H}}_{n-1}\cdots {\tilde{H}}_2{\tilde{H}}_{1}A=R \\ {\tilde{H}}_1^{-1}{\tilde{H}}_2^{-1}\cdots {\tilde{H}}_{n-1}^{-1}{\tilde{H}}_n^{-1}{\tilde{H}}_n{\tilde{H}}_{n-1}\cdots {\tilde{H}}_2{\tilde{H}}_{1}A={\tilde{H}}_1^{-1}{\tilde{H}}_2^{-1}\cdots {\tilde{H}}_{n-1}^{-1}{\tilde{H}}_n^{-1}R \\ A={\tilde{H}}_1^{-1}{\tilde{H}}_2^{-1}\cdots {\tilde{H}}_{n-1}^{-1}{\tilde{H}}_n^{-1}R \\ A={\tilde{H}}_1{\tilde{H}}_2\cdots {\tilde{H}}_{n-1}{\tilde{H}}_{n}R \end{gathered}$$
  最终这些拼凑为$n$阶的矩阵连乘起来，就是Q矩阵，最终的结果就是R矩阵。

3 举例

  以这个矩阵为例子：
$$\left(\begin{array}{ccc}0& 3& 1\\ 0& 4& -2\\ 2& 1& 2\end{array}\right)$$
  第一次分解得到：
$$\left(\begin{array}{ccc}0& 3& 1\\ 0& 4& -2\\ 2& 1& 2\end{array}\right)=\left(\begin{array}{ccc}0& 0& 1\\ 0& 1& 0\\ 1& 0& 0\end{array}\right)\left(\begin{array}{ccc}2& 1& 2\\ 0& 4& -2\\ 0& 3& 1\end{array}\right)$$
  第二次分解得到：
$$\left(\begin{array}{ccc}2& 1& 2\\ 0& 4& -2\\ 0& 3& 1\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0.8& 0.6\\ 0& 0.6& -0.8\end{array}\right)\left(\begin{array}{ccc}2& 1& 2\\ 0& 5& -1\\ 0& 0& -2\end{array}\right)$$
  第三次分解得到：
$$\left(\begin{array}{ccc}2& 1& 2\\ 0& 5& -1\\ 0& 0& -2\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& -1\end{array}\right)\left(\begin{array}{ccc}2& 1& 2\\ 0& 5& -1\\ 0& 0& 2\end{array}\right)$$

4 Python实现

  代码实现如下：

    @staticmethod
    def householder_u(x, z, inner_product_matrix=None):
        x_len = vector_len(x, inner_product_matrix)
        xz = mul_num(z, x_len)
        complement = sub(x, xz)
        complement_len = vector_len(complement, inner_product_matrix)
        return mul_num(complement, 1 / complement_len)

    # 获取householder变换的u向量
    def get_householder_u(self, inner_product_matrix=None):
        a1 = self.__vectors[0]
        n = len(a1)
        e1 = [0] * n
        e1[0] = 1
        return Matrix.householder_u(a1, e1, inner_product_matrix)

    # 获取householder矩阵
    def householder_matrix(self, inner_product_matrix=None):
        n = len(self.__vectors)
        unit = Matrix.unit_matrix(n)
        u = self.get_householder_u(inner_product_matrix)
        u_matrix = Matrix([u])
        h = Matrix(unit) - u_matrix * u_matrix.transpose_matrix() * 2
        return h

    # householder变换
    def householder(self, inner_product_matrix=None):
        n = len(self.__vectors)
        r = self
        q = Matrix(Matrix.unit_matrix(n))
        for i in range(n):
            # 子矩阵
            sub_matrix_array = [r.__vectors[j][j:] for j in range(i, n)]
            sub_matrix = Matrix(sub_matrix_array)
            sub_h = sub_matrix.householder_matrix(inner_product_matrix)
            # 左上角补1
            h_array = Matrix.unit_matrix(n)
            for j in range(i, n):
                for k in range(i,n):
                    h_array[j][k] = sub_h.__vectors[j-i][k-i]
            h = Matrix(h_array)
            r = h * r
            print((h * r).to_latex(), '=',h.to_latex(), r.to_latex())
            # sub_r 的内容
            q = q * h
        return q, r