基于Householder变换的约化对称矩阵为三对角形式的算法详解
题目描述
本题目讲解如何利用Householder变换将一个实对称矩阵 \(A\) 通过正交相似变换（即 \(Q^TAQ = T\)）化为一个实对称三对角矩阵 \(T\) 的算法过程。这个过程是许多计算对称矩阵特征值算法（如QR算法、分而治之法）的关键预处理步骤，因为它能在保持矩阵所有特征值不变的前提下，将问题简化到具有特殊稀疏结构（仅有主对角线和次对角线非零）的三对角矩阵上，从而极大提高后续计算的效率。
核心目标：给定一个 \(n \times n\) 的实对称矩阵 \(A\)，构造一系列Householder反射矩阵 \(P_1, P_2, ..., P_{n-2}\)，使得 \(T = P_{n-2} ... P_2 P_1 A P_1 P_2 ... P_{n-2}\) 是一个实对称三对角矩阵。
解题过程详解
第一步：理解Householder变换的核心思想
- Householder反射：对于一个非零向量 \(v \in \mathbb{R}^n\)，可以定义一个Householder矩阵（反射矩阵）为：
\[ P = I - 2 \frac{vv^T}{v^T v} \]
其中 $I$ 是单位矩阵。矩阵 $P$ 是**对称** ($P^T = P$) 且**正交** ($P^{-1} = P^T = P$) 的。其几何意义是将任意向量 $x$ 关于以向量 $v$ 为法向的超平面进行反射。
- 关键性质：给定一个向量 \(x\) 和一个单位向量 \(e_1 = (1,0,...,0)^T\)，我们可以构造一个特定的Householder矩阵 \(P\)，使得 \(Px = \sigma e_1\)，其中 \(\sigma = \pm ||x||_2\)。这个性质使得Householder变换能“将向量的一部分分量消为零”，这是将矩阵化为特定形式（如上Hessenberg形、三对角形）的基础。
第二步：将对称矩阵三对角化的整体策略
对于一个 \(n \times n\) 的对称矩阵 \(A\)，其对称性保证了其下三角和上三角部分是对称的。三对角化的目标是，在保持对称性的前提下，依次将矩阵的每一行和每一列中，处于“次对角线”之外的元素（即距离主对角线“距离”大于1的元素）消为零。
这个过程是迭代进行的，从第一列（和第一行）开始，一直到第 \(n-2\) 列（和 \(n-2\) 行）结束。每一步，我们处理矩阵的一个“缩小”的子块。
第三步：逐步推导算法（以4x4矩阵为例）
设初始对称矩阵 \(A^{(0)} = A\)。
第k=1步：处理第一列（和第一行）
-
目标：将 \(A^{(0)}\) 的第一列中，第3行到第n行（即从 \(a_{31}\) 到 \(a_{n1}\)）的元素消为零。由于对称性，对应的第一行这些位置也会被自动消零。
-
操作范围：我们只关注 \(A^{(0)}\) 的第一列的下方部分 \(x_1 = \begin{bmatrix} a_{21}^{(0)} \\ a_{31}^{(0)} \\ \vdots \\ a_{n1}^{(0)} \end{bmatrix}\)。我们要找一个Householder变换 \(P_1'\)，作用在 \((n-1)\) 维空间上，使得 \(P_1' x_1 = \beta_1 e_1^{(n-1)}\)，其中 \(e_1^{(n-1)} = (1,0,...,0)^T \in \mathbb{R}^{n-1}\)，\(\beta_1 = \pm ||x_1||_2\)。
-
构造：计算 \(u_1 = x_1 - \beta_1 e_1^{(n-1)}\)，则 \((n-1)\) 维的Householder矩阵为 \(P_1' = I^{(n-1)} - 2 \frac{u_1 u_1^T}{u_1^T u_1}\)。然后，我们将其“嵌入”到 \(n\) 维空间，构造 \(n \times n\) 的Householder矩阵：
\[ P_1 = \begin{bmatrix} 1 & 0 \\ 0 & P_1' \end{bmatrix} \]
- 变换矩阵：
\[ A^{(1)} = P_1 A^{(0)} P_1 \]
这个正交相似变换保持对称性和特征值。
* $P_1 A^{(0)}$ 会将 $A^{(0)}$ 第一列的 $x_1$ 部分（第2到n行）变换为 $(\beta_1, 0, ..., 0)^T$。
* 由于 $P_1$ 是块对角形式，右乘 $P_1$ 相当于用 $P_1'$ 变换乘积结果的第2到n行和第2到n列。因为左乘时已经将第一列的特定部分化零，右乘会确保对称性，并将第一行的对应部分也化零。结果是：
\[ A^{(1)} = \begin{bmatrix} a_{11}^{(1)} & \beta_1 & 0 & 0 \\ \beta_1 & a_{22}^{(1)} & * & * \\ 0 & * & * & * \\ 0 & * & * & * \end{bmatrix} \]
其中*表示该位置在本次变换后还未被“清理”。
第k=2步：处理第二列（和剩余子块的第一行）
-
目标：现在聚焦于 \(A^{(1)}\) 的右下角 \((n-1) \times (n-1)\) 的子矩阵。我们想将这个子矩阵的其第一列中，从第3个元素到最后一个元素（即相对于全局矩阵的第4到n行，第3列）消为零。
-
操作范围：考虑子矩阵 \(A^{(1)}(2:n, 2:n)\)。取这个子矩阵的第一列（对应于全局矩阵的第3到n行，第2列）的下方部分：\(x_2 = \begin{bmatrix} a_{32}^{(1)} \\ a_{42}^{(1)} \\ \vdots \\ a_{n2}^{(1)} \end{bmatrix}\)。
-
构造：类似于第一步，构造一个 \((n-2)\) 维的Householder矩阵 \(P_2'\)，使得 \(P_2' x_2 = \beta_2 e_1^{(n-2)}\)，其中 \(\beta_2 = \pm ||x_2||_2\)。然后嵌入到 \(n\) 维空间：
\[ P_2 = \begin{bmatrix} I_2 & 0 \\ 0 & P_2' \end{bmatrix} \]
其中 $I_2$ 是 $2 \times 2$ 的单位矩阵。
- 变换矩阵：
\[ A^{(2)} = P_2 A^{(1)} P_2 \]
这个变换不会破坏第一步已得到的三对角结构（前两行和前两列），只作用于右下角的 $(n-1) \times (n-1)$ 子块。结果是：
\[ A^{(2)} = \begin{bmatrix} a_{11}^{(1)} & \beta_1 & 0 & 0 \\ \beta_1 & a_{22}^{(2)} & \beta_2 & 0 \\ 0 & \beta_2 & a_{33}^{(2)} & * \\ 0 & 0 & * & * \end{bmatrix} \]
第k=3步（对于n=4的情况，这是最后一步）
对于 \(4 \times 4\) 矩阵，经过前两步，我们只剩下一个 \(2 \times 2\) 的右下角子块 \(A^{(2)}(3:4, 3:4)\)。对于一个 \(2 \times 2\) 的对称矩阵，它本身就是一个三对角矩阵（只有对角线和紧邻的上下对角线）。因此，对于 \(n=4\)，我们只需要 \(n-2=2\) 步就完成了三对角化。结果矩阵 \(T = A^{(2)}\) 已经是三对角形式。
第四步：算法的一般形式和输出
对于一般的 \(n \times n\) 对称矩阵 \(A\)，算法需要 \(n-2\) 步。
- 初始化：\(A^{(0)} = A\)。
- 循环：对于 \(k = 1, 2, ..., n-2\)：
a. 令 \(x_k = A^{(k-1)}(k+1:n, k)\) （当前列的待处理部分）。
b. 计算 \(\beta_k = \text{sign}(x_{k,1}) \cdot ||x_k||_2\)，其中 \(x_{k,1}\) 是 \(x_k\) 的第一个分量。这个符号选择（通常取与 \(x_{k,1}\) 相反）是为了数值稳定性，避免在构造 \(u_k = x_k + \beta_k e_1\) 时产生抵消。
c. 计算反射向量 \(u_k = x_k + \beta_k e_1^{(n-k)}\) （注意这里用 \(+\) 是因为我们通常定义 \(Px = -\beta e_1\) 的形式，与第二步描述等价但实现细节不同，目标一致）。
d. 计算缩放因子 \(\tau_k = 2 / (u_k^T u_k)\)。
e. 构造对应的 \(n\) 维 Householder 矩阵 \(P_k = I - \tau_k \begin{bmatrix} 0 \\ u_k \end{bmatrix} \begin{bmatrix} 0 & u_k^T \end{bmatrix}\)，其结构是左上角 \(k \times k\) 块为单位阵。
f. 更新矩阵：\(A^{(k)} = P_k A^{(k-1)} P_k\)。在实际计算中，这一步通过高效的矩阵-向量和向量-向量运算完成，无需显式形成 \(P_k\)：
* 计算 \(w = \tau_k A^{(k-1)} u_k\) （利用了矩阵对称性，只计算一半）。
* 计算 \(p = w - (0.5 * \tau_k * (w^T u_k)) * u_k\)。
* 更新右下角子矩阵：\(A^{(k)}(k+1:n, k+1:n) = A^{(k-1)}(k+1:n, k+1:n) - u_k p^T - p u_k^T\)。
* 设置三对角元素：\(A^{(k)}(k, k+1) = A^{(k)}(k+1, k) = -\beta_k\)，并将第 \(k\) 列和第 \(k\) 行的对应位置（k+2:n）置零。 - 输出：最终得到的 \(T = A^{(n-2)}\) 是一个实对称三对角矩阵，其主对角线元素为 \(T_{ii} = a_{ii}^{(n-2)}\)，次对角线元素为 \(T_{i,i+1} = T_{i+1,i} = -\beta_i\)。
总结：这个算法通过 \(n-2\) 次精心设计的Householder反射，逐步将原始对称矩阵的“边角”元素清除，最终得到一个具有相同特征值的三对角矩阵。它是数值线性代数中一个经典、稳定且高效的标准预处理过程。