【数值计算方法（黄明游）】矩阵特征值与特征向量的计算（五）：Householder方法【理论到程序】

  矩阵的特征值（eigenvalue）和特征向量（eigenvector）在很多应用中都具有重要的数学和物理意义。Householder 矩阵和变换提供了一种有效的方式，通过反射变换将一个向量映射到一个标准的方向，这对于一些数值计算问题具有重要的意义。

  本文将详细介绍Householder方法的基本原理和步骤，并给出其Python实现。

一、Jacobi 旋转法

  Jacobi 旋转法的每一次迭代中，需要选择一个非对角元素最大的位置，然后构造相应的旋转矩阵，进行相似变换，使得矩阵逐渐对角化。

二、Jacobi 过关法

  Jacobi 过关法（Jacobi’s threshold method）是 Jacobi 旋转法的一种改进版本，其主要目的是减少计算工作和提高运行速度。该方法通过动态调整阈值，并根据阈值对非对角元素进行选择性的旋转变换，以逐步对角化对称矩阵。

三、Householder 方法

  如果对任意向量 $z$，我们可以将其分解为与 $u$ 平行的分量 $au$ 和与 $u$ 正交的分量 $bv$，即 $z=au+bv$，那么 Householder 变换会将 $z$ 变换为 $-au+bv$。这个变换可以理解为镜面反射，它不改变向量在与 $u$ 正交的平面上的投影，但将向量沿着 $u$ 的方向反射。数学表达式为：

$$Hz=au+bv\to -au+bv$$

  这个性质使得 Householder 变换在一些数值计算的应用中非常有用，例如 QR 分解等。

1. 旋转变换

  在 Householder 方法中，通过一系列的正交相似变换，可以将实对称矩阵 (A) 转化为三对角矩阵。不同于 Jacobi 旋转法($a_{ij}=a_{ji}=0$)，Householder 方法的旋转矩阵选择的角度使得$c_{ik}=c_{kj}=0$。

a. 旋转变换的选择

  对于实对称矩阵 $A$ 中的元素$a_{ik}$ ，选择适当的旋转角度 $\theta$ ，可以使得$a_{ik}$ 变为零。具体而言，选择 $\theta$ 使得：
$$c_{ik}=c_{kj}=a_{ik}\cos (\theta )+a_{jk}\sin (\theta )=0$$

  通过这样的选择，我们可以构造一个旋转矩阵 $P_{i,j,k}$，该矩阵对应的正交相似变换可以将 $c_{ik}$ 变为零。这一过程实现了对实对称矩阵的正交相似变换，使得某些元素变为零，逐步实现了将矩阵转化为三对角形式。

b. 旋转变换的顺序

  在进行 Householder 变换时，旋转的顺序很重要。通常，可以选择按列进行变换。例如，对于 $i=2,\dots ,n-1$，可以依次选择 $j=i+1,i+2,\dots ,n$，然后对每一对 $(i,j)$ 进行正交相似变换，将 $a_{ik}$ 变为零。整个过程的迭代将会逐步将实对称矩阵 $A$ 转化为三对角矩阵 $C$。

2. Householder矩阵（Householder Matrix）

a. H矩阵的定义

  设 $u$为单位向量，即 $\Vert u\Vert =1$。定义 Householder 矩阵 $H=I-2u{u}^T$，其中 $I$ 为单位矩阵。这个矩阵具有以下性质：

• 对称性： $H^T=H$，即 Householder 矩阵是对称的。
• 正交性： $H^{T}H=I$，即 Householder 矩阵是正交矩阵。
• 保范性： 对于任意非零向量 $x$ 和 $y$，如果 $\Vert x{\Vert }^2=\Vert y{\Vert }^2$，则存在 Householder 矩阵 $H$，使得 $Hx=y$。
  • 考虑 Householder 矩阵对向量 $u$ 的作用：$Hu=(I-2u{u}^T)u=-u$。这说明 Householder 矩阵将向量 $u$ 反射到其负向量上。
  • 对于任何与 $u$ 正交的向量 $v$，有 $Hv=(I-2u{u}^T)v=v$，即 Householder 矩阵保持与 $u$ 正交的向量不变。
  • 因此对任意向量 $z$，设 $z=au+bv$，那么 Householder 变换会将 $z$ 变换为 $-au+bv$，数学表达式为：
    $$Hz=a(Hu)+b(Hv)\to -au+bv$$

b. H变换的几何解释

  可以将 Householder 变换视为镜面反射。考虑 $u$ 为反射面上的单位法向量。对于任意$z$，Householder 变换将 $z$ 的投影反射到 $-u$方向，同时保持投影在反射面上。

c. H变换的应用场景

1. 矩阵三对角化： 在计算线性代数中，Householder 变换常用于将矩阵化为三对角形式，以便更容易进行特征值计算等操作。
2. QR 分解： Householder 变换是计算 QR 分解的基本工具，用于将矩阵分解为一个正交矩阵和一个上三角矩阵的乘积。

3. H变换过程详解

a. 过程介绍

  对于矩阵$A$ 的某一列向量 $\mathbf{a}=(a_1,a_2,\dots ,a_n{)}^T$，如果我们想将向量的后 $n-(r+1)$个分量化为零，即将 $\mathbf{a}$ 变为 $\mathbf{c}=(a_1,a_2,\dots ,a_r,-\text{sign}(a_{r+1})s,0,\dots ,0{)}^T$，其中 $s=\Vert a{\Vert }_2$ （从$r+1$计算到$n$）且 $a_{r+1}\ne 0$，则可以引入 Householder 矩阵 $H$，使得$Ha=c$。Householder 矩阵的计算方式如下：
$$\mathbf{w}=\mathbf{a}-\text{sign}(a_{r+1})s{\mathbf{e}}_{r+1}$$
  其中，${\mathbf{e}}_{r+1}$ 是单位向量 $(0,0,\dots ,0,1,0,\dots ,0{)}^T$，具体位置在第 $r+1$ 个。

b. 细节解析

• Householder 矩阵的构造：
  • 通过 Householder 变换，构造 Householder 矩阵 $H$，将某一列 $a_j$的 $r+1$ 到 $n$ 个分量化为零。

• 计算过程的稳定性：

  • 将 $a$ 的 $r+1$ 到 $n$ 个分量的符号设定为 $-\text{sign}(a_{r+1})$，以增强计算的稳定性。

• 计算相似三对角矩阵：

  • 将 $A$ 逐列进行正交相似变换，得到 $A_1,A_2,\dots ,A_{n-1}$。
  • 最终得到相似三对角矩阵 $G=A_n=H_{n-2}\cdot \dots \cdot H_2\cdot H_1\cdot A$。

• 实际计算中的优化：

  • 实际计算中，无需形成所有的 Householder 矩阵，也无需进行矩阵乘法运算，可以直接在原矩阵上进行计算。
    

4. H变换例题解析

  给定矩阵：
$$A=\left[\begin{array}{cccc}1& 2& 1& 2\\ 2& 2& -1& 1\\ 1& -1& 1& 1\\ 2& 1& 1& 1\end{array}\right]$$

最终的三对角矩阵 $A_2$：

• $A_2$的形式为
  $$A_2=\left[\begin{array}{cccc}1& -3& 0& 0\\ -3& 2.3333& -0.4714& 0\\ 0& -0.4714& 1.1667& -1.5003\\ 0& 0& -1.5003& 0.5002\end{array}\right]$$

  这样，通过选择 $r=1$ 和 $r=2$ 进行两次 Householder 变换，矩阵 $A$ 被成功地化为三对角形式 $A_2$。

四、Python实现

import numpy as np


def householder_matrix(v):
    """
    给定向量 v，返回 Householder 变换矩阵 H
    """
    v = np.array(v, dtype=float)
    if np.linalg.norm(v) == 0:
        raise ValueError("无效的输入向量，它应该是非零的。")

    v = v / np.linalg.norm(v)
    H = np.eye(len(v)) - 2 * np.outer(v, v)
    return H


def householder_reduction(A):
    """
    对矩阵 A 执行 Householder 变换，将其化为三对角形式。
    """
    m, n = A.shape
    if m != n:
        raise ValueError("输入矩阵 A 必须是方阵。")

    Q = np.eye(m)  # 初始化正交矩阵 Q

    for k in range(n - 2):
        x = A[k + 1:, k]
        e1 = np.zeros_like(x)
        e1[0] = 1.0
        v = np.sign(x[0]) * np.linalg.norm(x) * e1 + x
        v = v / np.linalg.norm(v)

        H = np.eye(m)
        H[k + 1:, k + 1:] -= 2.0 * np.outer(v, v)
        Q = np.dot(Q, H)
        A = np.dot(H, np.dot(A, H))

    return Q, A


# 示例矩阵
A = np.array([[1, 2, 1, 2],
              [2, 2, -1, 1],
              [1, -1, 1, 1],
              [2, 1, 1, 1]], dtype=float)

# Householder 变换
Q, tridiagonal_A = householder_reduction(A)
np.set_printoptions(precision=4, suppress=True)
print("正交矩阵 Q:")
print(Q)
print("\n三对角矩阵:")
print(tridiagonal_A)

调试过程

• 第一次：
  

• 第二次：
  

• final：