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

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

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

一、Jacobi 旋转法

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

二、Jacobi 过关法

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

三、Householder 方法

  如果对任意向量

，我们可以将其分解为与

平行的分量

和与

正交的分量

，即

，那么 Householder 变换会将

变换为

。这个变换可以理解为镜面反射，它不改变向量在与

正交的平面上的投影，但将向量沿着

的方向反射。数学表达式为：

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

1. 旋转变换

  在 Householder 方法中，通过一系列的正交相似变换，可以将实对称矩阵 (A) 转化为三对角矩阵。不同于 Jacobi 旋转法(

)，Householder 方法的旋转矩阵选择的角度使得

。

a. 旋转变换的选择

  对于实对称矩阵

中的元素

，选择适当的旋转角度

，可以使得

变为零。具体而言，选择

使得：

  通过这样的选择，我们可以构造一个旋转矩阵

，该矩阵对应的正交相似变换可以将

变为零。这一过程实现了对实对称矩阵的正交相似变换，使得某些元素变为零，逐步实现了将矩阵转化为三对角形式。

b. 旋转变换的顺序

  在进行 Householder 变换时，旋转的顺序很重要。通常，可以选择按列进行变换。例如，对于

，可以依次选择

，然后对每一对

进行正交相似变换，将

变为零。整个过程的迭代将会逐步将实对称矩阵

转化为三对角矩阵

。

2. Householder矩阵（Householder Matrix）

a. H矩阵的定义

  设

为单位向量，即

。定义 Householder 矩阵

，其中

为单位矩阵。这个矩阵具有以下性质：

• 对称性：

，即 Householder 矩阵是对称的。

• 正交性：

，即 Householder 矩阵是正交矩阵。

• 保范性： 对于任意非零向量

和

，如果

，则存在 Householder 矩阵

，使得

。

• 考虑 Householder 矩阵对向量

的作用：

。这说明 Householder 矩阵将向量

反射到其负向量上。

• 对于任何与

正交的向量

，有

，即 Householder 矩阵保持与

正交的向量不变。

• 因此对任意向量

，设

，那么 Householder 变换会将

变换为

，数学表达式为：

b. H变换的几何解释

  可以将 Householder 变换视为镜面反射。考虑

为反射面上的单位法向量。对于任意

，Householder 变换将

的投影反射到

方向，同时保持投影在反射面上。

c. H变换的应用场景

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

3. H变换过程详解

a. 过程介绍

  对于矩阵

的某一列向量

，如果我们想将向量的后

个分量化为零，即将

变为

，其中

（从

计算到

）且

，则可以引入 Householder 矩阵

，使得

。Householder 矩阵的计算方式如下：

  其中，

是单位向量

，具体位置在第

个。

b. 细节解析

• Householder 矩阵的构造：
  • 通过 Householder 变换，构造 Householder 矩阵

  $$H$$

  ，将某一列

  $$a_j$$

  的

  $$r+1$$

  到

  $$n$$

  个分量化为零。

• 计算过程的稳定性：
  • 将

  $$a$$

  的

  $$r+1$$

  到

  $$n$$

  个分量的符号设定为

  $$-\text{sign}(a_{r+1})$$

  ，以增强计算的稳定性。

• 计算相似三对角矩阵：
  • 将

  $$A$$

  逐列进行正交相似变换，得到

  $$A_1, A_2, \ldots, A_{n-1}$$

  。

  • 最终得到相似三对角矩阵

  $$G = A_n = H_{n-2} \cdot \ldots \cdot H_2 \cdot H_1 \cdot A$$

  。

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

  

4. H变换例题解析

  给定矩阵：

最终的三对角矩阵

：

的形式为

  这样，通过选择

和

进行两次 Householder 变换，矩阵

被成功地化为三对角形式

。

四、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：