ch08+矩阵特征计算+本硕博讲义

8.1 特征值性质和估计 详细讲解与证明

我将以多年数值分析教学与研究的经验，从定义出发，逐步推导每一个核心结论，所有证明依据均以粗体标注，确保逻辑严谨、脉络清晰。

一、特征值问题的基本定义

1. 特征值与特征向量

设矩阵 \(A \in \mathbb{R}^{n \times n}\)，特征值问题是求复数 \(\lambda \in \mathbb{C}\) 和非零向量 \(x \in \mathbb{C}^n\)，使得：

\[Ax = \lambda x \tag{8.1} \]

其中 \(\lambda\) 称为矩阵 \(A\) 的特征值，非零向量 \(x\) 称为 \(A\) 属于特征值 \(\lambda\) 的特征向量。

关键说明：必须强调 \(x \neq 0\)。若允许 \(x=0\)，则对任意 \(\lambda\) 等式都成立，失去了特征值问题的意义。

2. 特征多项式与特征方程

将式(8.1)移项得：

\[(\lambda I - A)x = 0 \]

这是一个齐次线性方程组。依据：齐次线性方程组有非零解的充要条件是系数矩阵的行列式为零，因此：

\[\det(\lambda I - A) = 0 \]

定义矩阵 \(A\) 的特征多项式为：

\[p(\lambda) = \det(\lambda I - A) = \begin{vmatrix} \lambda - a_{11} & -a_{12} & \dots & -a_{1n} \\ -a_{21} & \lambda - a_{22} & \dots & -a_{2n} \\ \vdots & \vdots & & \vdots \\ -a_{n1} & -a_{n2} & \dots & \lambda - a_{nn} \end{vmatrix}\]

展开后得到 \(n\) 次多项式：

\[p(\lambda) = \lambda^n + c_1\lambda^{n-1} + \dots + c_{n-1}\lambda + c_n \tag{8.2} \]

方程 \(p(\lambda)=0\) 称为矩阵 \(A\) 的特征方程。依据：代数基本定理，\(n\) 次代数方程在复数域中有且仅有 \(n\) 个根（重根按重数计算），记为 \(\lambda_1, \lambda_2, \dots, \lambda_n\)，因此特征多项式可因式分解为：

\[p(\lambda) = (\lambda - \lambda_1)(\lambda - \lambda_2)\dots(\lambda - \lambda_n) \]

二、迹与行列式的特征值表示

1. 核心结论

• 矩阵的迹（主对角线元素之和）等于所有特征值之和：
  \[\operatorname{tr}A = \sum_{i=1}^n a_{ii} = \sum_{i=1}^n \lambda_i \]

• 矩阵的行列式等于所有特征值之积：
  \[\det A = \lambda_1\lambda_2\dots\lambda_n \]

2. 详细推导

将特征多项式的两种展开式进行比较：

1. 行列式展开式：

   \[p(\lambda) = \lambda^n - (a_{11}+a_{22}+\dots+a_{nn})\lambda^{n-1} + \dots + (-1)^n\det A \]

   （主对角线乘积贡献最高次项和次高次项，常数项由 \(\lambda=0\) 代入得 \(\det(-A)=(-1)^n\det A\)）

2. 因式分解式：

   \[p(\lambda) = \lambda^n - (\lambda_1+\lambda_2+\dots+\lambda_n)\lambda^{n-1} + \dots + (-1)^n\lambda_1\lambda_2\dots\lambda_n \]

依据：多项式恒等的充要条件是对应次数的系数相等，因此：

• 比较 \(\lambda^{n-1}\) 项系数：\(-(a_{11}+\dots+a_{nn}) = -(\lambda_1+\dots+\lambda_n)\)，即得迹的公式。
• 比较常数项：\((-1)^n\det A = (-1)^n\lambda_1\lambda_2\dots\lambda_n\)，约去 \((-1)^n\) 即得行列式的公式。

三、特征值与特征向量的基本性质

性质1：转置矩阵与原矩阵有相同的特征值

证明依据：行列式的性质 \(\det(B^T) = \det(B)\)

\[\det(\lambda I - A^T) = \det\left((\lambda I - A)^T\right) = \det(\lambda I - A) \]

因此 \(A^T\) 与 \(A\) 有相同的特征多项式，故有相同的特征值。

备注：特征值相同，但特征向量不一定相同。\(A^T\) 的特征向量满足 \(A^T y = \lambda y\)，与 \(A\) 的特征向量 \(x\) 一般不同。

性质2：非奇异矩阵的逆矩阵的特征值

若 \(A\) 非奇异，则 \(A^{-1}\) 的特征值为 \(\lambda^{-1}\)，特征向量仍为 \(x\)。

证明依据：非奇异矩阵的定义（\(\det A \neq 0\) 等价于所有特征值非零）
已知 \(Ax = \lambda x\)，因 \(A\) 非奇异，故 \(\lambda \neq 0\)。两边左乘 \(A^{-1}\) 得：

\[x = \lambda A^{-1}x \implies A^{-1}x = \frac{1}{\lambda}x \]

证毕。

性质3：相似矩阵有相同的特征多项式

若 \(B = S^{-1}AS\)（\(S\) 非奇异），则 \(A\) 与 \(B\) 有相同的特征多项式，从而有相同的特征值。

证明依据：行列式的乘法性质 \(\det(CD) = \det C \det D\) 及 \(\det(S^{-1}) = 1/\det S\)

\[\begin{align*} \det(\lambda I - B) &= \det(\lambda I - S^{-1}AS) \\ &= \det(S^{-1}\lambda I S - S^{-1}AS) \\ &= \det\left(S^{-1}(\lambda I - A)S\right) \\ &= \det(S^{-1}) \det(\lambda I - A) \det S \\ &= \frac{1}{\det S} \cdot \det(\lambda I - A) \cdot \det S \\ &= \det(\lambda I - A) \end{align*} \]

备注：相似矩阵特征值相同，但特征向量不同。若 \(x\) 是 \(A\) 的特征向量，则 \(S^{-1}x\) 是 \(B\) 的特征向量。

性质4：实矩阵的复特征值成对共轭出现

实矩阵的复特征值一定相互共轭成对出现，且对应的复特征向量的实部和虚部线性无关。

证明依据：实矩阵的共轭等于自身（\(\bar{A} = A\)）及不同特征值对应的特征向量线性无关

1. 设 \(\lambda = \alpha + i\beta\) 是实矩阵 \(A\) 的复特征值，\(x = u + iv\) 是对应的复特征向量，则 \(Ax = \lambda x\)。
2. 两边取共轭得 \(\bar{A}\bar{x} = \bar{\lambda}\bar{x}\)，因 \(\bar{A}=A\)，故 \(A\bar{x} = \bar{\lambda}\bar{x}\)，即 \(\bar{\lambda} = \alpha - i\beta\) 也是 \(A\) 的特征值，对应特征向量 \(\bar{x} = u - iv\)。
3. 假设 \(u\) 和 \(v\) 线性相关，则存在不全为零的实数 \(k_1,k_2\) 使 \(k_1u + k_2v = 0\)，从而 \(k_1x + k_2\bar{x} = 0\)。但 \(x\) 和 \(\bar{x}\) 是不同特征值对应的特征向量，依据：不同特征值对应的特征向量线性无关，故 \(k_1=k_2=0\)，矛盾。因此 \(u\) 和 \(v\) 线性无关。

四、特征值的运算性质（定理8.1）

设 \(\lambda\) 是 \(A \in \mathbb{R}^{n \times n}\) 的特征值，\(x\) 是对应的特征向量（\(Ax = \lambda x, x \neq 0\)），则：

1. \(c\lambda\) 是 \(cA\) 的特征值（\(c\) 为非零常数）
2. \(\lambda - \mu\) 是 \(A - \mu I\) 的特征值（\(\mu\) 为任意常数）
3. \(\lambda^k\) 是 \(A^k\) 的特征值（\(k\) 为正整数）

证明：

1. 两边乘以 \(c\) 得 \(cAx = c\lambda x\)，即 \((cA)x = (c\lambda)x\)。
2. \((A - \mu I)x = Ax - \mu x = \lambda x - \mu x = (\lambda - \mu)x\)。
3. 依据：数学归纳法原理
   • 基例：\(k=1\) 时显然成立。
   • 归纳假设：设 \(k=m\) 时 \(A^m x = \lambda^m x\) 成立。
   • 归纳步骤：\(A^{m+1}x = A(A^m x) = A(\lambda^m x) = \lambda^m (Ax) = \lambda^m \cdot \lambda x = \lambda^{m+1}x\)。

推论（多项式矩阵的特征值）：若 \(f(t) = a_0 + a_1t + \dots + a_k t^k\) 是多项式，则 \(f(\lambda)\) 是矩阵多项式 \(f(A) = a_0I + a_1A + \dots + a_k A^k\) 的特征值，特征向量仍为 \(x\)。

五、矩阵可对角化的条件（定理8.2）

定理8.2(1)：可对角化的充要条件

\(A \in \mathbb{R}^{n \times n}\) 可对角化（即存在非奇异矩阵 \(P\) 使 \(P^{-1}AP = \operatorname{diag}(\lambda_1, \lambda_2, \dots, \lambda_n)\)）的充分必要条件是 \(A\) 具有 \(n\) 个线性无关的特征向量。

证明：

• 必要性：设 \(P^{-1}AP = \Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_n)\)，则 \(AP = P\Lambda\)。
  将 \(P\) 按列分块为 \(P = [p_1, p_2, \dots, p_n]\)，则：

  \[AP = [Ap_1, Ap_2, \dots, Ap_n], \quad P\Lambda = [\lambda_1 p_1, \lambda_2 p_2, \dots, \lambda_n p_n] \]

  故 \(Ap_i = \lambda_i p_i\)（\(i=1,\dots,n\)）。依据：非奇异矩阵的列向量线性无关，因此 \(p_1, \dots, p_n\) 是 \(A\) 的 \(n\) 个线性无关的特征向量。

• 充分性：设 \(A\) 有 \(n\) 个线性无关的特征向量 \(p_1, \dots, p_n\)，对应特征值 \(\lambda_1, \dots, \lambda_n\)。
  令 \(P = [p_1, \dots, p_n]\)，则 \(P\) 非奇异，且：

  \[AP = [Ap_1, \dots, Ap_n] = [\lambda_1 p_1, \dots, \lambda_n p_n] = P\Lambda \]

  两边左乘 \(P^{-1}\) 得 \(P^{-1}AP = \Lambda\)，即 \(A\) 可对角化。

定理8.2(2)：不同特征值对应的特征向量线性无关

若 \(A\) 有 \(m\) 个（\(m \leq n\)）不同的特征值 \(\lambda_1, \dots, \lambda_m\)，则对应的特征向量 \(x_1, \dots, x_m\) 线性无关。

证明依据：数学归纳法原理

• 基例：\(m=1\) 时，\(x_1 \neq 0\)，显然线性无关。
• 归纳假设：设 \(m=k\) 时结论成立，即 \(k\) 个不同特征值对应的特征向量线性无关。
• 归纳步骤：设有常数 \(c_1, \dots, c_{k+1}\) 使：
  \[c_1 x_1 + \dots + c_k x_k + c_{k+1} x_{k+1} = 0 \tag{1} \]
  两边左乘 \(A\) 得：
  \[c_1 \lambda_1 x_1 + \dots + c_k \lambda_k x_k + c_{k+1} \lambda_{k+1} x_{k+1} = 0 \tag{2} \]
  用(2)式减去 \(\lambda_{k+1} \times (1)\) 式得：
  \[c_1(\lambda_1 - \lambda_{k+1})x_1 + \dots + c_k(\lambda_k - \lambda_{k+1})x_k = 0 \]
  由归纳假设，\(x_1, \dots, x_k\) 线性无关，故 \(c_i(\lambda_i - \lambda_{k+1}) = 0\)（\(i=1,\dots,k\)）。
  因 \(\lambda_i \neq \lambda_{k+1}\)，故 \(c_i = 0\)（\(i=1,\dots,k\)）。代入(1)式得 \(c_{k+1}x_{k+1}=0\)，而 \(x_{k+1} \neq 0\)，故 \(c_{k+1}=0\)。
  因此 \(x_1, \dots, x_{k+1}\) 线性无关。

推论：若 \(A\) 有 \(n\) 个不同的特征值，则 \(A\) 一定可对角化。

六、对称矩阵的瑞利商性质（定理8.3）

1. 瑞利商定义

设 \(A \in \mathbb{R}^{n \times n}\) 为对称矩阵，对任意非零向量 \(x \in \mathbb{R}^n\)，定义：

\[R(x) = \frac{(Ax, x)}{(x, x)} \]

称为矩阵 \(A\) 的瑞利(Rayleigh)商。

2. 定理8.3内容

设对称矩阵 \(A\) 的特征值按降序排列为 \(\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_n\)，则：

1. 对任意非零向量 \(x \in \mathbb{R}^n\)，有 \(\lambda_n \leq R(x) \leq \lambda_1\)
2. \(\lambda_1 = \max_{\substack{x \in \mathbb{R}^n \\ x \neq 0}} R(x)\)，\(\lambda_n = \min_{\substack{x \in \mathbb{R}^n \\ x \neq 0}} R(x)\)

3. 详细证明

证明依据：实对称矩阵的正交对角化定理（实对称矩阵一定存在正交矩阵 \(Q\)，使 \(Q^T A Q = \operatorname{diag}(\lambda_1, \dots, \lambda_n)\)，且 \(Q\) 的列向量是标准正交特征向量 \(q_1, \dots, q_n\)，满足 \(A q_i = \lambda_i q_i\)，\((q_i, q_j) = \delta_{ij}\)）

对任意非零向量 \(x \in \mathbb{R}^n\)，可表示为标准正交基的线性组合：

\[x = \alpha_1 q_1 + \alpha_2 q_2 + \dots + \alpha_n q_n \]

其中 \(\alpha_i = (x, q_i)\) 为实数，且 \(\|x\|^2 = (x, x) = \sum_{i=1}^n \alpha_i^2 \neq 0\)。

计算内积：

\[Ax = A\left(\sum_{i=1}^n \alpha_i q_i\right) = \sum_{i=1}^n \alpha_i A q_i = \sum_{i=1}^n \alpha_i \lambda_i q_i \]

\[(Ax, x) = \left(\sum_{i=1}^n \alpha_i \lambda_i q_i, \sum_{j=1}^n \alpha_j q_j\right) = \sum_{i=1}^n \sum_{j=1}^n \alpha_i \alpha_j \lambda_i (q_i, q_j) = \sum_{i=1}^n \alpha_i^2 \lambda_i \]

因此瑞利商为：

\[R(x) = \frac{\sum_{i=1}^n \alpha_i^2 \lambda_i}{\sum_{i=1}^n \alpha_i^2} \]

证明结论(1)

因 \(\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_n\)，故对所有 \(i\) 有 \(\lambda_n \leq \lambda_i \leq \lambda_1\)。两边乘以 \(\alpha_i^2 \geq 0\) 并求和得：

\[\lambda_n \sum_{i=1}^n \alpha_i^2 \leq \sum_{i=1}^n \alpha_i^2 \lambda_i \leq \lambda_1 \sum_{i=1}^n \alpha_i^2 \]

两边除以正数 \(\sum_{i=1}^n \alpha_i^2\)，即得：

\[\lambda_n \leq R(x) \leq \lambda_1 \]

证明结论(2)

由结论(1)知 \(R(x) \leq \lambda_1\) 对所有非零 \(x\) 成立，故 \(\max R(x) \leq \lambda_1\)。
另一方面，取 \(x = q_1\)（对应最大特征值的特征向量），则：

\[R(q_1) = \frac{(A q_1, q_1)}{(q_1, q_1)} = \frac{\lambda_1 (q_1, q_1)}{(q_1, q_1)} = \lambda_1 \]

因此 \(\max_{\substack{x \neq 0}} R(x) = \lambda_1\)。

同理可证 \(\min_{\substack{x \neq 0}} R(x) = \lambda_n\)（取 \(x = q_n\) 即可）。

重要意义：瑞利商为对称矩阵的特征值估计提供了理论基础，无需求解特征方程即可得到特征值的上下界。

七、知识点系统总结表

知识点类别	核心内容	证明依据	关键备注与应用
基本定义	特征值 \(\lambda\) 和特征向量 \(x\) 满足 \(Ax = \lambda x\)（\(x \neq 0\)）	齐次线性方程组有非零解的充要条件	特征向量必须非零；特征值可以是复数
特征多项式	\(p(\lambda) = \det(\lambda I - A)\)，特征方程 \(p(\lambda)=0\)	行列式定义、代数基本定理	\(n\) 阶矩阵有 \(n\) 个特征值（重根按重数计）
迹与行列式	\(\operatorname{tr}A = \sum a_{ii} = \sum \lambda_i\) \(\det A = \prod \lambda_i\)	多项式恒等系数相等	快速计算迹和行列式；判断矩阵奇异性（\(\det A=0\) 等价于有零特征值）
转置矩阵性质	\(A^T\) 与 \(A\) 有相同特征值	行列式转置不变性	特征向量一般不同
逆矩阵性质	非奇异矩阵 \(A^{-1}\) 的特征值为 \(\lambda^{-1}\)	非奇异矩阵定义、矩阵乘法	特征向量与原矩阵相同
相似矩阵性质	相似矩阵有相同特征多项式和特征值	行列式乘法性质	特征向量不同；相似变换不改变矩阵的迹和行列式
实矩阵复特征值	复特征值共轭成对出现，对应特征向量的实部虚部线性无关	实矩阵共轭不变性、不同特征值特征向量线性无关	实矩阵的特征值集合关于实轴对称
特征值运算性质	\(cA\) 特征值为 \(c\lambda\) \(A-\mu I\) 特征值为 \(\lambda-\mu\) \(A^k\) 特征值为 \(\lambda^k\)	特征值定义、数学归纳法	可推广到任意矩阵多项式 \(f(A)\)，其特征值为 \(f(\lambda)\)
可对角化条件	充要条件：有 \(n\) 个线性无关的特征向量 充分条件：有 \(n\) 个不同特征值	矩阵分块乘法、数学归纳法	可对角化矩阵可通过相似变换化为对角矩阵，简化计算
瑞利商性质	对称矩阵的瑞利商 \(R(x)=(Ax,x)/(x,x)\) 满足 \(\lambda_n \leq R(x) \leq \lambda_1\) 最大最小特征值分别是瑞利商的最大值和最小值	实对称矩阵正交对角化定理	用于对称矩阵特征值的估计和迭代算法（如瑞利商迭代法）

---

8.1.2 特征值估计与扰动 详细讲解与证明

我将继续以资深数值分析研究员的身份，系统讲解格什戈林圆盘定理这一特征值估计的核心工具，所有关键证明依据均以粗体标注，并补充实用技巧与应用场景。

一、格什戈林圆盘的定义（定义8.1）

设 \(A=(a_{ij}) \in \mathbb{C}^{n \times n}\)，定义：

1. 行和：\(r_i = \sum_{\substack{j=1 \\ j \neq i}}^n |a_{ij}| \quad (i=1,2,\dots,n)\)
   即矩阵第 \(i\) 行所有非对角元的绝对值之和。
2. 格什戈林圆盘：复平面上以对角元 \(a_{ii}\) 为圆心，以行和 \(r_i\) 为半径的圆盘
   \[D_i = \{ z \in \mathbb{C} \mid |z - a_{ii}| \leq r_i \} \]

3. 格什戈林圆盘集合：所有 \(n\) 个圆盘的并集 \(\bigcup_{i=1}^n D_i\)。

几何意义：格什戈林圆盘定理告诉我们，矩阵的所有特征值都落在这些圆盘的并集中，无需计算特征多项式即可快速估计特征值范围。

二、格什戈林圆盘定理（定理8.4）

定理8.4(1)：特征值的包含性

矩阵 \(A\) 的每一个特征值都必属于它的某一个格什戈林圆盘之中，即：

\[\forall \lambda \in \sigma(A), \exists i \in \{1,2,\dots,n\}, \quad |\lambda - a_{ii}| \leq r_i \]

其中 \(\sigma(A)\) 表示矩阵 \(A\) 的特征值集合。

详细证明

设 \(\lambda\) 是 \(A\) 的任意一个特征值，\(x=(x_1,x_2,\dots,x_n)^T \neq 0\) 是对应的特征向量，则：

\[Ax = \lambda x \tag{1} \]

依据：特征值与特征向量的定义

由于 \(x \neq 0\)，其分量中至少有一个非零。设第 \(k\) 个分量的绝对值最大，即：

\[|x_k| = \max_{1 \leq i \leq n} |x_i| = \|x\|_\infty \neq 0 \]

依据：无穷范数的定义，非零向量的无穷范数必大于零

考虑方程(1)的第 \(k\) 个分量：

\[\sum_{j=1}^n a_{kj} x_j = \lambda x_k \]

将对角元项移到右边：

\[(\lambda - a_{kk})x_k = \sum_{\substack{j=1 \\ j \neq k}}^n a_{kj} x_j \tag{2} \]

对等式(2)两边取绝对值：

\[|\lambda - a_{kk}| \cdot |x_k| = \left| \sum_{\substack{j \neq k}} a_{kj} x_j \right| \]

依据：复数绝对值的乘法性质 \(|ab|=|a||b|\)

对右边应用三角不等式 \(|\sum z_j| \leq \sum |z_j|\)：

\[|\lambda - a_{kk}| \cdot |x_k| \leq \sum_{\substack{j \neq k}} |a_{kj}| \cdot |x_j| \tag{3} \]

由于 \(|x_k|\) 是最大分量，故对所有 \(j \neq k\) 有 \(|x_j| \leq |x_k|\)，代入(3)式得：

\[|\lambda - a_{kk}| \cdot |x_k| \leq \sum_{\substack{j \neq k}} |a_{kj}| \cdot |x_k| = |x_k| \cdot r_k \]

因为 \(|x_k| \neq 0\)，两边同时除以 \(|x_k|\)，得到：

\[|\lambda - a_{kk}| \leq r_k \]

这就证明了特征值 \(\lambda\) 必属于第 \(k\) 个格什戈林圆盘 \(D_k\)。证毕。

关键洞察：特征值 \(\lambda\) 所在的圆盘下标 \(k\)，恰好是对应特征向量 \(x\) 中绝对值最大分量的下标。

定理8.4(2)：分离圆盘的特征值计数

如果矩阵 \(A\) 有 \(m\) 个格什戈林圆盘组成一个连通的并集 \(S\)，且 \(S\) 与其余 \(n-m\) 个格什戈林圆盘是分离的（即不相交），则 \(S\) 内恰好包含 \(A\) 的 \(m\) 个特征值（按代数重数计算）。

特别地：如果某一个格什戈林圆盘 \(D_i\) 与其他所有圆盘都分离（称为孤立圆盘），则 \(D_i\) 中精确地包含 \(A\) 的一个特征值。

证明思路（教材未给出，补充核心思想）

依据：特征值的连续性原理
考虑带参数的矩阵族 \(A(t) = D + t(A-D)\)，其中 \(D = \operatorname{diag}(a_{11},a_{22},\dots,a_{nn})\) 是 \(A\) 的对角部分，\(t \in [0,1]\)。

• 当 \(t=0\) 时，\(A(0)=D\)，其特征值就是对角元 \(a_{11},a_{22},\dots,a_{nn}\)，每个圆盘 \(D_i\) 中恰好有一个特征值。
• 当 \(t\) 从 \(0\) 连续变化到 \(1\) 时，\(A(t)\) 的特征值也连续变化，且每个特征值始终在某个格什戈林圆盘中。
• 如果 \(m\) 个圆盘组成连通区域 \(S\) 且与其他圆盘分离，则这 \(m\) 个特征值不会跑出 \(S\)，也不会有其他特征值进入 \(S\)，因此 \(S\) 中恰好有 \(m\) 个特征值。

三、利用相似变换改进特征值估计

格什戈林圆盘定理给出的特征值范围有时会比较宽松，特别是当圆盘重叠严重时。我们可以利用相似矩阵有相同特征值的性质，通过适当的相似变换来缩小某些圆盘的半径，从而得到更精确的估计。

基本方法

选取非奇异对角矩阵：

\[D = \operatorname{diag}(\alpha_1, \alpha_2, \dots, \alpha_n), \quad \alpha_i > 0 \]

则矩阵 \(B = D^{-1}AD\) 与 \(A\) 相似，有相同的特征值。

计算 \(B\) 的元素：

\[B = (b_{ij}) = \left( \frac{\alpha_j}{\alpha_i} a_{ij} \right) \]

\(B\) 的第 \(i\) 个格什戈林圆盘的半径为：

\[r_i' = \sum_{\substack{j \neq i}} |b_{ij}| = \sum_{\substack{j \neq i}} \frac{\alpha_j}{\alpha_i} |a_{ij}| \]

通过选择合适的 \(\alpha_i\)，可以有针对性地缩小某些圆盘的半径，甚至使重叠的圆盘分离。

实用技巧

• 若要缩小第 \(k\) 个圆盘的半径，可取 \(\alpha_k > 1\)，其余 \(\alpha_i=1\)。此时第 \(k\) 个圆盘的半径会缩小，而其他圆盘的半径会略有增大。
• 若要放大第 \(k\) 个圆盘的半径，可取 \(\alpha_k < 1\)，其余 \(\alpha_i=1\)。

示例

考虑矩阵：

\[A = \begin{pmatrix} 4 & 1 \\ 1 & 2 \end{pmatrix} \]

• 原格什戈林圆盘：\(D_1: |z-4| \leq 1\)，\(D_2: |z-2| \leq 1\)，两个圆盘不相交，故 \(D_1\) 中有一个特征值，\(D_2\) 中有一个特征值。
• 实际特征值：\(\lambda_1=3+\sqrt{2} \approx 4.414\)，\(\lambda_2=3-\sqrt{2} \approx 1.586\)，符合定理结论。

再考虑矩阵：

\[A = \begin{pmatrix} 3 & 2 \\ 2 & 3 \end{pmatrix} \]

• 原格什戈林圆盘：\(D_1: |z-3| \leq 2\)，\(D_2: |z-3| \leq 2\)，两个圆盘完全重合，只能得到特征值在 \([1,5]\) 范围内。
• 取 \(D = \operatorname{diag}(2,1)\)，则：
  \[B = D^{-1}AD = \begin{pmatrix} 3 & 1 \\ 4 & 3 \end{pmatrix} \]

• \(B\) 的格什戈林圆盘：\(D_1': |z-3| \leq 1\)，\(D_2': |z-3| \leq 4\)。
• 此时 \(D_1'\) 是孤立圆盘，故其中必有一个特征值，另一个特征值在 \(D_2'\) 中。
• 实际特征值：\(\lambda_1=5\)，\(\lambda_2=1\)，\(\lambda_2=1\) 确实在 \(D_1'\) 中，\(\lambda_1=5\) 在 \(D_2'\) 中，估计精度显著提高。

四、知识点系统总结表

知识点类别	核心内容	证明依据	关键应用与备注
格什戈林圆盘定义	第 \(i\) 个圆盘 \(D_i: |z-a_{ii}| \leq \sum_{j \neq i} |a_{ij}|\)	复数模的几何意义	复平面上的圆盘区域，圆心为对角元，半径为非对角元绝对值之和
圆盘定理(1)	所有特征值都在格什戈林圆盘的并集中	特征值定义、三角不等式、无穷范数性质	无需计算特征多项式即可快速估计特征值范围
圆盘定理(2)	\(m\) 个连通且分离的圆盘恰好包含 \(m\) 个特征值	特征值的连续性原理	孤立圆盘必含一个特征值，可用于确定特征值的精确位置
相似变换改进估计	用对角矩阵 \(D\) 做相似变换 \(D^{-1}AD\)，调整圆盘半径	相似矩阵有相同特征值	可针对性缩小某些圆盘半径，分离重叠圆盘，提高估计精度
列圆盘定理	所有特征值也在列格什戈林圆盘的并集中 列圆盘：$	z-a_	\leq \sum_

补充说明：列圆盘定理是行圆盘定理的直接推论，因为 \(A^T\) 与 \(A\) 有相同的特征值，而 \(A^T\) 的行圆盘就是 \(A\) 的列圆盘。因此，矩阵的特征值必属于行圆盘并集与列圆盘并集的交集，这通常能给出更精确的估计。

---

例8.1 格什戈林圆盘定理应用 详细解析

我将以资深数值分析研究员的视角，逐步骤拆解这道经典例题，所有关键步骤的理论依据均以粗体标注，并深入解释相似变换的设计思路与技巧。

一、问题重述

估计矩阵

\[A = \begin{pmatrix} 4 & 1 & 0 \\ 1 & 0 & -1 \\ 1 & 1 & -4 \end{pmatrix} \]

的特征值范围。

二、第一步：原矩阵的格什戈林圆盘估计

1. 计算格什戈林圆盘

依据：格什戈林圆盘定义8.1

• 第1行：对角元 \(a_{11}=4\)，非对角元绝对值和 \(r_1=|1|+|0|=1\)
  圆盘 \(D_1: |\lambda - 4| \leq 1\)，即实轴上的区间 \([3, 5]\)
• 第2行：对角元 \(a_{22}=0\)，非对角元绝对值和 \(r_2=|1|+|-1|=2\)
  圆盘 \(D_2: |\lambda| \leq 2\)，即实轴上的区间 \([-2, 2]\)
• 第3行：对角元 \(a_{33}=-4\)，非对角元绝对值和 \(r_3=|1|+|1|=2\)
  圆盘 \(D_3: |\lambda + 4| \leq 2\)，即实轴上的区间 \([-6, -2]\)

2. 分析圆盘连通性与特征值分布

依据：格什戈林圆盘定理8.4(2)

• \(D_1=[3,5]\) 与 \(D_2\)、\(D_3\) 均不相交，是孤立圆盘，因此 \(D_1\) 中恰好包含 \(A\) 的一个实特征值 \(\lambda_1\)，即 \(3 \leq \lambda_1 \leq 5\)。
• \(D_2=[-2,2]\) 与 \(D_3=[-6,-2]\) 在点 \(\lambda=-2\) 处相交，组成一个连通的并集 \([-6,2]\)，因此这个连通区域内包含 \(A\) 的另外两个特征值 \(\lambda_2, \lambda_3\)。

此时的估计结果比较粗糙，无法区分 \(\lambda_2\) 和 \(\lambda_3\) 的具体范围，需要通过相似变换进一步改进。

三、第二步：利用相似变换改进估计

1. 相似变换的设计思路

依据：相似矩阵有相同的特征值
我们的目标是让原来连通的 \(D_2\) 和 \(D_3\) 分离。观察到：

• \(D_3\) 的圆心在 \(-4\)，半径为 \(2\)
• \(D_2\) 的圆心在 \(0\)，半径为 \(2\)
• 两者仅在 \(\lambda=-2\) 处接触，只要稍微缩小 \(D_3\) 的半径，同时略微增大 \(D_2\) 的半径，就能让它们分离。

选择非奇异对角矩阵的逆为：

\[D^{-1} = \begin{pmatrix} 1 & & \\ & 1 & \\ & & 0.9 \end{pmatrix} \]

即 \(D = \begin{pmatrix} 1 & & \\ & 1 & \\ & & 1/0.9 \end{pmatrix}\)

2. 计算相似变换后的矩阵 \(A_1\)

\[A_1 = D^{-1}AD \]

依据：对角矩阵乘法规则 \((D^{-1}AD)_{ij} = (D^{-1})_{ii} \cdot a_{ij} \cdot D_{jj}\)

计算各元素：

• \(a_{11}' = 1 \cdot 4 \cdot 1 = 4\)，\(a_{12}' = 1 \cdot 1 \cdot 1 = 1\)，\(a_{13}' = 1 \cdot 0 \cdot (1/0.9) = 0\)
• \(a_{21}' = 1 \cdot 1 \cdot 1 = 1\)，\(a_{22}' = 1 \cdot 0 \cdot 1 = 0\)，\(a_{23}' = 1 \cdot (-1) \cdot (1/0.9) = -10/9\)
• \(a_{31}' = 0.9 \cdot 1 \cdot 1 = 0.9\)，\(a_{32}' = 0.9 \cdot 1 \cdot 1 = 0.9\)，\(a_{33}' = 0.9 \cdot (-4) \cdot (1/0.9) = -4\)

因此：

\[A_1 = \begin{pmatrix} 4 & 1 & 0 \\ 1 & 0 & -10/9 \\ 0.9 & 0.9 & -4 \end{pmatrix} \]

3. 计算 \(A_1\) 的格什戈林圆盘

依据：格什戈林圆盘定义8.1

• 第1行：\(r_1' = |1| + |0| = 1\)，圆盘 \(E_1: |\lambda - 4| \leq 1\)，即 \([3, 5]\)（与原 \(D_1\) 相同）
• 第2行：\(r_2' = |1| + |-10/9| = 19/9 \approx 2.111\)，圆盘 \(E_2: |\lambda| \leq 19/9\)，即 \([-19/9, 19/9] \approx [-2.111, 2.111]\)
• 第3行：\(r_3' = |0.9| + |0.9| = 1.8\)，圆盘 \(E_3: |\lambda + 4| \leq 1.8\)，即 \([-5.8, -2.2]\)

4. 分析变换后的圆盘连通性

依据：格什戈林圆盘定理8.4(2)

• \(E_1=[3,5]\) 仍然是孤立圆盘，包含一个特征值。
• \(E_2=[-2.111, 2.111]\) 与 \(E_3=[-5.8, -2.2]\) 之间有一个间隙 \((-2.2, -2.111)\)，不再连通，各自成为孤立圆盘。
• 因此，三个圆盘都是孤立的，每个圆盘内恰好包含 \(A\) 的一个实特征值。

5. 改进后的特征值估计

\[\begin{cases} 3 \leq \lambda_1 \leq 5 \\ -\dfrac{19}{9} \leq \lambda_2 \leq \dfrac{19}{9} \quad (\approx -2.111 \leq \lambda_2 \leq 2.111) \\ -5.8 \leq \lambda_3 \leq -2.2 \end{cases} \]

四、精确验证

1. 计算特征方程

\[\begin{align*} \det(\lambda I - A) &= \begin{vmatrix} \lambda - 4 & -1 & 0 \\ -1 & \lambda & 1 \\ -1 & -1 & \lambda + 4 \end{vmatrix} \\ &= (\lambda - 4) \begin{vmatrix} \lambda & 1 \\ -1 & \lambda + 4 \end{vmatrix} - (-1) \begin{vmatrix} -1 & 1 \\ -1 & \lambda + 4 \end{vmatrix} + 0 \\ &= (\lambda - 4)(\lambda(\lambda + 4) + 1) + (-(\lambda + 4) + 1) \\ &= (\lambda - 4)(\lambda^2 + 4\lambda + 1) - (\lambda + 3) \\ &= \lambda^3 + 4\lambda^2 + \lambda - 4\lambda^2 - 16\lambda - 4 - \lambda - 3 \\ &= \lambda^3 - 16\lambda - 7 \end{align*} \]

2. 求解特征值

解三次方程 \(\lambda^3 - 16\lambda - 7 = 0\)，得到三个实根：

• \(\lambda_1 \approx 4.2030\)（在 \([3,5]\) 内）
• \(\lambda_2 \approx -0.4429\)（在 \([-19/9, 19/9]\) 内）
• \(\lambda_3 \approx -3.7601\)（在 \([-5.8, -2.2]\) 内）

所有特征值都落在我们估计的区间内，验证了格什戈林圆盘定理的正确性。

五、例题核心要点总结

1. 基础估计：直接应用格什戈林圆盘定理可以快速得到特征值的大致范围，但当圆盘连通时，无法确定每个特征值的具体位置。
2. 相似变换技巧：通过选择适当的对角矩阵进行相似变换，可以有针对性地缩小某些圆盘的半径，使连通的圆盘分离，从而得到更精确的估计。
3. 变换原则：若要缩小第 \(k\) 个圆盘的半径，可令 \(D^{-1}\) 的第 \(k\) 个对角元小于1，这会使第 \(k\) 行非对角元的绝对值减小，同时使第 \(k\) 列非对角元的绝对值增大。
4. 实对称矩阵特性：本题中的矩阵虽然不是对称矩阵，但所有特征值都是实数，因此格什戈林圆盘在实轴上的区间就是特征值的范围。

---

8.1.2 特征值扰动与Bauer-Fike定理 详细讲解

我将继续以资深数值分析研究员的身份，系统讲解特征值扰动理论的核心——Bauer-Fike定理，所有关键证明依据均以粗体标注，并深入阐释特征值条件数的物理意义与应用价值。

一、特征值扰动问题的背景

在实际工程与科学计算中，矩阵 \(A\) 的元素往往存在误差（如测量误差、舍入误差等），这相当于给矩阵 \(A\) 加上了一个微小的扰动矩阵 \(E\)，得到扰动矩阵 \(A+E\)。我们关心的核心问题是：当 \(E\) 很小时，\(A+E\) 的特征值与 \(A\) 的特征值相差多少？

这个问题的答案决定了特征值计算结果的可靠性：如果微小的扰动会导致特征值发生巨大变化，我们称该矩阵的特征值是病态的；反之则是良态的。

二、Bauer-Fike定理（定理8.5）

Bauer-Fike定理是可对角化矩阵特征值扰动分析的基础定理，它给出了扰动后特征值与原特征值之间距离的统一上界。

1. 定理内容

设矩阵 \(A \in \mathbb{R}^{n \times n}\) 可对角化，即存在非奇异矩阵 \(P\) 使得 \(P^{-1}AP = D = \operatorname{diag}(\lambda_1, \lambda_2, \dots, \lambda_n)\)，其中 \(\lambda_i\) 是 \(A\) 的特征值。设 \(\mu\) 是扰动矩阵 \(A+E\) 的任意一个特征值，则有：

\[\min_{\lambda \in \sigma(A)} |\lambda - \mu| \leq \|P^{-1}\|_p \|P\|_p \|E\|_p \tag{8.3} \]

其中 \(\sigma(A)\) 表示 \(A\) 的特征值集合，\(\|\cdot\|_p\) 为矩阵的 \(p\) 范数（\(p=1,2,\infty\)）。

2. 详细证明

证明思路：从扰动后的特征方程出发，通过相似变换将 \(A\) 对角化，再利用范数的性质推导不等式。

• 步骤1：处理特殊情况
  若 \(\mu \in \sigma(A)\)，则 \(\min_{\lambda \in \sigma(A)} |\lambda - \mu| = 0\)，不等式(8.3)显然成立。因此只需证明 \(\mu \notin \sigma(A)\) 的情况。

• 步骤2：从特征方程出发
  设 \(x \neq 0\) 是 \(A+E\) 对应于特征值 \(\mu\) 的特征向量，则：

  \[(A+E - \mu I)x = 0 \]

  依据：特征值与特征向量的定义

• 步骤3：引入相似变换
  已知 \(A = PDP^{-1}\)，代入上式得：

  \[(PDP^{-1} + E - \mu I)x = 0 \]

  两边左乘 \(P^{-1}\)：

  \[(D - \mu I)(P^{-1}x) + (P^{-1}EP)(P^{-1}x) = 0 \]

  依据：矩阵乘法的结合律与分配律

• 步骤4：移项并利用可逆性
  因 \(\mu \notin \sigma(A)\)，故 \(D - \mu I\) 是对角矩阵且对角元均非零，因此可逆。移项得：

  \[P^{-1}x = -(D - \mu I)^{-1}(P^{-1}EP)(P^{-1}x) \]

  注意到 \(x \neq 0\) 且 \(P\) 非奇异，故 \(P^{-1}x \neq 0\)。

• 步骤5：两边取范数
  对上式两边取 \(p\) 范数：

  \[\|P^{-1}x\|_p = \|(D - \mu I)^{-1}(P^{-1}EP)(P^{-1}x)\|_p \]

  依据：范数的相容性 \(\|ABC\| \leq \|A\| \|B\| \|C\|\)，得：

  \[\|P^{-1}x\|_p \leq \|(D - \mu I)^{-1}\|_p \|P^{-1}EP\|_p \|P^{-1}x\|_p \]

  两边除以正数 \(\|P^{-1}x\|_p\)，得：

  \[1 \leq \|(D - \mu I)^{-1}\|_p \|P^{-1}EP\|_p \tag{1} \]

• 步骤6：计算对角矩阵的范数
  依据：对角矩阵的 \(p\) 范数性质，对角矩阵的 \(p\) 范数（\(p=1,2,\infty\)）等于其对角元绝对值的最大值。因此：

  \[\|(D - \mu I)^{-1}\|_p = \max_{1 \leq i \leq n} \left| \frac{1}{\lambda_i - \mu} \right| = \frac{1}{\min_{1 \leq i \leq n} |\lambda_i - \mu|} \tag{2} \]

• 步骤7：估计乘积项的范数
  再次利用范数的相容性：

  \[\|P^{-1}EP\|_p \leq \|P^{-1}\|_p \|E\|_p \|P\|_p \tag{3} \]

• 步骤8：整理得到结论
  将(2)和(3)代入(1)式：

  \[1 \leq \frac{1}{\min |\lambda_i - \mu|} \cdot \|P^{-1}\|_p \|P\|_p \|E\|_p \]

  两边乘以 \(\min |\lambda_i - \mu|\)，即得：

  \[\min_{\lambda \in \sigma(A)} |\lambda - \mu| \leq \|P^{-1}\|_p \|P\|_p \|E\|_p \]

  证毕。

3. 定理的关键说明

1. 适用条件：定理仅适用于可对角化矩阵。对于不可对角化矩阵（亏损矩阵），特征值扰动可能会大得多，需要更复杂的扰动理论。
2. 范数无关性：定理对 \(p=1,2,\infty\) 范数均成立，这是其重要优点。
3. 上界性质：定理给出的是特征值扰动的上界，实际扰动可能远小于这个上界。

三、特征值问题的条件数

1. 定义与意义

由Bauer-Fike定理可知，\(\|P^{-1}\| \|P\| = \operatorname{cond}(P)\) 是特征值扰动的放大系数。但将 \(A\) 对角化的相似变换矩阵 \(P\) 不是唯一的，因此我们定义：

\[\nu(A) = \inf \left\{ \operatorname{cond}(P) \mid P^{-1}AP = \operatorname{diag}(\lambda_1, \lambda_2, \dots, \lambda_n) \right\} \]

称为矩阵 \(A\) 的特征值问题的条件数。

• 若 \(\nu(A)\) 很小（接近1），则矩阵 \(A\) 的特征值是良态的，微小的扰动只会带来特征值的微小变化。
• 若 \(\nu(A)\) 很大，则矩阵 \(A\) 的特征值是病态的，微小的扰动可能导致特征值的巨大变化。

2. 重要特例：正规矩阵

当 \(A\) 是正规矩阵（满足 \(A^T A = A A^T\)，包括对称矩阵、反对称矩阵、正交矩阵等）时，存在正交矩阵 \(Q\) 使得 \(Q^T A Q = D\)。此时 \(\|Q\|_2 = \|Q^T\|_2 = 1\)，故 \(\operatorname{cond}_2(Q) = 1\)，因此 \(\nu_2(A) = 1\)。

结论：正规矩阵的特征值对扰动总是良态的，这是数值计算中非常重要的性质。

3. 与线性方程组条件数的区别

特征值问题的条件数 \(\nu(A)\) 与解线性方程组时矩阵 \(A\) 的条件数 \(\operatorname{cond}(A) = \|A\| \|A^{-1}\|\) 是两个完全不同的概念，两者可能相差极大。

教材示例：二阶对角矩阵 \(A = \operatorname{diag}(1, 10^{-10})\)

• 特征值条件数：取 \(P=I\)，则 \(\operatorname{cond}(P)=1\)，故 \(\nu(A) \leq 1\)，特征值非常良态。
• 线性方程组条件数：\(\operatorname{cond}_\infty(A) = \|A\|_\infty \|A^{-1}\|_\infty = 1 \times 10^{10} = 10^{10}\)，解线性方程组时非常病态。

四、矩阵特征值数值方法的引入

1. 解析方法的局限性

当矩阵阶数 \(n=2,3\) 时，我们可以通过行列式展开求出特征多项式 \(p(\lambda) = \det(\lambda I - A)\)，然后求解代数方程 \(p(\lambda)=0\) 得到特征值。但当 \(n\) 较大时，这种方法存在严重缺陷：

1. 计算量巨大：计算 \(n\) 阶行列式的工作量为 \(O(n!)\)，当 \(n \geq 5\) 时已难以承受。
2. 数值稳定性差：高次多项式的根对系数的微小变化极其敏感，即使特征多项式的系数只有很小的误差，也可能导致根的巨大偏差。

因此，对于高阶矩阵，解析方法完全不适用，必须研究专门的数值方法。

2. 数值方法的分类

本章将介绍计算机上常用的两类特征值数值方法：

1. 幂法及反幂法（迭代法）：用于求解按模最大或最小的特征值及其对应的特征向量。这类方法计算量小，适合求解大型稀疏矩阵的部分特征值。
2. 正交相似变换方法：用于求解矩阵的全部特征值和特征向量。其中最著名的是QR方法，它是目前求解中小型稠密矩阵全部特征值的最有效方法。

五、知识点系统总结表

知识点类别	核心内容	证明依据	关键备注与应用
特征值扰动问题	研究矩阵微小扰动对特征值的影响，判断特征值的良态/病态	特征值定义、矩阵范数性质	是数值计算结果可靠性分析的基础
Bauer-Fike定理	可对角化矩阵的特征值扰动满足 \(\min |\lambda - \mu| \leq \operatorname{cond}(P) |E|\)	相似变换、范数相容性、对角矩阵范数性质	给出了特征值扰动的统一上界，适用于p=1,2,∞范数
特征值条件数	\(\nu(A) = \inf \{\operatorname{cond}(P) \mid P^{-1}AP=D\}\)，衡量特征值对扰动的敏感程度	Bauer-Fike定理	\(\nu(A)\) 越小，特征值越良态；正规矩阵的 \(\nu_2(A)=1\)
条件数对比	特征值条件数与线性方程组条件数是不同概念，数值可能相差极大	条件数定义	不能用线性方程组的条件数判断特征值的敏感性
数值方法引入	高阶矩阵特征值求解不能用解析方法，必须用数值方法	计算复杂度分析、数值稳定性分析	主要分为迭代法（幂法/反幂法）和正交变换法（QR方法）

---

8.2 幂法及反幂法（一）：幂法 详细讲解与证明

我将继续以资深数值分析研究员的身份，系统讲解求解矩阵主特征值的经典迭代法——幂法，所有关键证明依据均以粗体标注，深入剖析其收敛原理、数值稳定性问题及实用算法设计。

一、幂法的基本概念与适用场景

1. 定义

幂法是一种计算矩阵主特征值（按模最大的特征值）及其对应特征向量的迭代方法。

2. 核心优势

• 特别适用于大型稀疏矩阵（如有限元、有限差分法产生的矩阵）
• 算法简单，易于实现，计算量小（每步迭代仅需一次矩阵-向量乘法，复杂度为 \(O(n^2)\)）
• 存储需求低，无需存储整个矩阵的全部元素，只需存储非零元素

3. 基本假设

设矩阵 \(A \in \mathbb{R}^{n \times n}\) 可对角化，即存在 \(n\) 个线性无关的特征向量 \(x_1, x_2, \dots, x_n\)，对应的特征值满足：

\[|\lambda_1| > |\lambda_2| \geq |\lambda_3| \geq \dots \geq |\lambda_n| \]

其中 \(\lambda_1\) 称为主特征值，对应的特征向量 \(x_1\) 称为主特征向量。

二、幂法的基本原理（理论迭代）

1. 迭代思想

任取非零初始向量 \(v_0 \in \mathbb{R}^n\)，由于 \(x_1, \dots, x_n\) 是一组基，故 \(v_0\) 可唯一表示为：

\[v_0 = \alpha_1 x_1 + \alpha_2 x_2 + \dots + \alpha_n x_n \]

依据：线性空间基的定义

用矩阵 \(A\) 反复左乘 \(v_0\)，得到向量序列：

\[\begin{align*} v_1 &= A v_0 = \alpha_1 \lambda_1 x_1 + \alpha_2 \lambda_2 x_2 + \dots + \alpha_n \lambda_n x_n \\ v_2 &= A v_1 = A^2 v_0 = \alpha_1 \lambda_1^2 x_1 + \alpha_2 \lambda_2^2 x_2 + \dots + \alpha_n \lambda_n^2 x_n \\ &\vdots \\ v_k &= A^k v_0 = \alpha_1 \lambda_1^k x_1 + \alpha_2 \lambda_2^k x_2 + \dots + \alpha_n \lambda_n^k x_n \end{align*} \]

将上式提取公因子 \(\lambda_1^k\)：

\[v_k = \lambda_1^k \left( \alpha_1 x_1 + \alpha_2 \left( \frac{\lambda_2}{\lambda_1} \right)^k x_2 + \dots + \alpha_n \left( \frac{\lambda_n}{\lambda_1} \right)^k x_n \right) \tag{8.4} \]

2. 收敛性分析

由假设 \(|\lambda_1| > |\lambda_i|\)（\(i=2,3,\dots,n\)），故：

\[\lim_{k \to \infty} \left( \frac{\lambda_i}{\lambda_1} \right)^k = 0, \quad i=2,3,\dots,n \]

依据：等比数列极限性质，公比绝对值小于1时极限为0

因此，当 \(k\) 充分大时，式(8.4)中除第一项外，其余各项均可忽略不计，得到近似式：

\[v_k \approx \alpha_1 \lambda_1^k x_1 \]

这表明：当 \(k \to \infty\) 时，向量序列 \(\{v_k\}\) 将越来越接近主特征向量 \(x_1\) 的方向，两者仅相差一个常数因子。

3. 主特征值的计算

考虑相邻两个迭代向量的对应分量之比：

\[\frac{(v_{k+1})_i}{(v_k)_i} \approx \frac{\alpha_1 \lambda_1^{k+1} (x_1)_i}{\alpha_1 \lambda_1^k (x_1)_i} = \lambda_1 \]

依据：向量分量的线性性质

因此，当 \(k\) 充分大时，相邻迭代向量对应分量的比值将收敛到主特征值 \(\lambda_1\)。

4. 收敛速度

幂法的收敛速度由比值：

\[r = \left| \frac{\lambda_2}{\lambda_1} \right| \]

决定。\(r\) 越小，收敛越快；当 \(r\) 接近1时，收敛会非常缓慢。

三、实用幂法：规范化迭代

1. 理论迭代的数值问题

上述理论迭代在计算机上无法直接实现，因为存在两个严重的数值问题：

• 溢出问题：当 \(|\lambda_1| > 1\) 时，\(\lambda_1^k\) 会随着 \(k\) 的增大而指数增长，导致 \(v_k\) 的分量趋于无穷大，发生数值溢出。
• 有效数字损失：当 \(|\lambda_1| < 1\) 时，\(\lambda_1^k\) 会随着 \(k\) 的增大而指数衰减，导致 \(v_k\) 的分量趋于零，丢失有效数字。

2. 规范化方法

为了解决上述问题，我们在每一步迭代中都对向量进行规范化处理，使向量的无穷范数（最大分量的绝对值）等于1。

无穷范数定义：对于向量 \(v=(v_1,v_2,\dots,v_n)^T\)，其无穷范数为：

\[\|v\|_\infty = \max_{1 \leq i \leq n} |v_i| \]

若有多个分量同时达到最大值，取下标最小的那个分量。

3. 规范化幂法的迭代公式

\[\begin{cases} u_0 = v_0 \neq 0 \quad (\text{通常取 } u_0=(1,1,\dots,1)^T) \\ v_k = A u_{k-1} \\ \mu_k = \max(v_k) = \|v_k\|_\infty \\ u_k = \frac{v_k}{\mu_k} \end{cases} \quad k=1,2,\dots \tag{8.5} \]

4. 规范化迭代的收敛性证明

证明思路：通过数学归纳法建立 \(u_k, v_k\) 与 \(A^k v_0\) 的关系，再代入式(8.4)分析极限。

• 步骤1：数学归纳法证明表达式

  • 基例：\(k=1\) 时，\(v_1=A u_0=A v_0\)，\(u_1=v_1/\max(v_1)=A v_0/\max(A v_0)\)，成立。
  • 归纳假设：设 \(k=m\) 时，有：
    \[v_m = \frac{A^m v_0}{\max(A^{m-1} v_0)}, \quad u_m = \frac{A^m v_0}{\max(A^m v_0)} \]

  • 归纳步骤：当 \(k=m+1\) 时：
    \[v_{m+1} = A u_m = A \cdot \frac{A^m v_0}{\max(A^m v_0)} = \frac{A^{m+1} v_0}{\max(A^m v_0)} \]
    \[u_{m+1} = \frac{v_{m+1}}{\max(v_{m+1})} = \frac{A^{m+1} v_0 / \max(A^m v_0)}{\max(A^{m+1} v_0) / \max(A^m v_0)} = \frac{A^{m+1} v_0}{\max(A^{m+1} v_0)} \]

  由数学归纳法，对所有 \(k \geq 1\)，表达式成立。

• 步骤2：分析 \(u_k\) 的极限
  将式(8.4)代入 \(u_k\) 的表达式：

  \[\begin{align*} u_k &= \frac{A^k v_0}{\max(A^k v_0)} \\ &= \frac{\lambda_1^k \left( \alpha_1 x_1 + \sum_{i=2}^n \alpha_i \left( \frac{\lambda_i}{\lambda_1} \right)^k x_i \right)}{\max\left( \lambda_1^k \left( \alpha_1 x_1 + \sum_{i=2}^n \alpha_i \left( \frac{\lambda_i}{\lambda_1} \right)^k x_i \right) \right)} \\ &= \frac{\alpha_1 x_1 + \sum_{i=2}^n \alpha_i \left( \frac{\lambda_i}{\lambda_1} \right)^k x_i}{\max\left( \alpha_1 x_1 + \sum_{i=2}^n \alpha_i \left( \frac{\lambda_i}{\lambda_1} \right)^k x_i \right)} \end{align*} \]

  当 \(k \to \infty\) 时，分子分母中的求和项均趋于零，故：

  \[\lim_{k \to \infty} u_k = \frac{\alpha_1 x_1}{\max(\alpha_1 x_1)} = \frac{x_1}{\max(x_1)} \]

  即 \(u_k\) 收敛到主特征向量 \(x_1\) 的规范化形式。

• 步骤3：分析 \(\mu_k\) 的极限

  \[\begin{align*} \mu_k &= \max(v_k) = \max\left( \frac{A^k v_0}{\max(A^{k-1} v_0)} \right) \\ &= \frac{\max\left( \lambda_1^k \left( \alpha_1 x_1 + \sum_{i=2}^n \alpha_i \left( \frac{\lambda_i}{\lambda_1} \right)^k x_i \right) \right)}{\max\left( \lambda_1^{k-1} \left( \alpha_1 x_1 + \sum_{i=2}^n \alpha_i \left( \frac{\lambda_i}{\lambda_1} \right)^{k-1} x_i \right) \right)} \\ &= \lambda_1 \cdot \frac{\max\left( \alpha_1 x_1 + \sum_{i=2}^n \alpha_i \left( \frac{\lambda_i}{\lambda_1} \right)^k x_i \right)}{\max\left( \alpha_1 x_1 + \sum_{i=2}^n \alpha_i \left( \frac{\lambda_i}{\lambda_1} \right)^{k-1} x_i \right)} \end{align*} \]

  当 \(k \to \infty\) 时，分子分母的比值趋于1，故：

  \[\lim_{k \to \infty} \mu_k = \lambda_1 \]

  即 \(\mu_k\) 收敛到主特征值 \(\lambda_1\)。

四、幂法的收敛定理（定理8.6）

1. 定理内容

设 \(A \in \mathbb{R}^{n \times n}\) 有 \(n\) 个线性无关的特征向量 \(x_1, x_2, \dots, x_n\)，对应的特征值满足：

\[|\lambda_1| > |\lambda_2| \geq |\lambda_3| \geq \dots \geq |\lambda_n| \]

若非零初始向量 \(v_0 = \sum_{i=1}^n \alpha_i x_i\) 中的系数 \(\alpha_1 \neq 0\)，则幂法(8.5)产生的序列满足：

1. \(\displaystyle \lim_{k \to \infty} u_k = \frac{x_1}{\max(x_1)}\)
2. \(\displaystyle \lim_{k \to \infty} \mu_k = \lambda_1\)

且收敛速度由比值 \(r = \left| \frac{\lambda_2}{\lambda_1} \right|\) 确定。

2. 主特征值为重根的情况

若主特征值是 \(r\) 重根，即：

\[\lambda_1 = \lambda_2 = \dots = \lambda_r, \quad |\lambda_r| > |\lambda_{r+1}| \geq \dots \geq |\lambda_n| \]

且初始向量 \(v_0\) 中前 \(r\) 个特征向量的系数 \(\alpha_1, \alpha_2, \dots, \alpha_r\) 不全为零，则幂法仍然收敛，且：

\[\lim_{k \to \infty} u_k = \frac{\sum_{i=1}^r \alpha_i x_i}{\max\left( \sum_{i=1}^r \alpha_i x_i \right)}, \quad \lim_{k \to \infty} \mu_k = \lambda_1 \]

说明：此时迭代得到的是主特征子空间 \(\operatorname{span}\{x_1, \dots, x_r\}\) 中的一个向量，而不是某个特定的特征向量，但特征值的估计仍然是准确的。

五、实用注意事项

1. 初始向量的选择

定理中要求 \(\alpha_1 \neq 0\)，这个条件在实际计算中几乎总是满足的：

• 即使初始向量的 \(\alpha_1 = 0\)，计算过程中的舍入误差也会不可避免地引入 \(x_1\) 方向的分量，使得迭代最终收敛。
• 为了加快收敛速度，通常选择初始向量为全1向量 \(u_0 = (1,1,\dots,1)^T\)。

2. 迭代终止条件

通常采用以下两种终止准则之一：

• 特征向量差：\(\|u_k - u_{k-1}\|_\infty < \varepsilon\)
• 特征值差：\(|\mu_k - \mu_{k-1}| < \varepsilon\)

其中 \(\varepsilon\) 是预设的精度要求（如 \(10^{-6}\)）。

3. 收敛加速

当 \(r = |\lambda_2/\lambda_1|\) 接近1时，幂法收敛很慢，此时可以采用原点平移法或瑞利商加速法来提高收敛速度，这将在后续内容中详细讲解。

六、知识点系统总结表

知识点类别	核心内容	证明依据	关键备注与应用
幂法定义	计算矩阵主特征值（按模最大）及对应特征向量的迭代法	特征值与特征向量定义	特别适合大型稀疏矩阵，计算量小，存储需求低
基本原理	反复用矩阵左乘初始向量，主特征值对应的分量将占主导	线性空间基的定义、等比数列极限性质	收敛速度由 \(r=|\lambda_2/\lambda_1|\) 决定，\(r\) 越小收敛越快
数值问题	理论迭代存在溢出和有效数字损失问题	指数函数的增长/衰减性质	必须引入规范化步骤
规范化方法	每步迭代用无穷范数将向量规范化，使最大分量为1	无穷范数定义	解决数值稳定性问题，不改变收敛性
收敛定理	若矩阵可对角化且主特征值按模严格最大，初始向量α₁≠0，则幂法收敛	数学归纳法、极限运算性质	规范化向量序列收敛到主特征向量，最大分量序列收敛到主特征值
重根情况	主特征值为重根时，幂法仍收敛到主特征值，但特征向量是子空间中的一个向量	线性组合的极限性质	特征值估计准确，特征向量不唯一
实用技巧	初始向量取全1向量，终止条件用相邻迭代差的范数	数值计算经验	收敛慢时可采用原点平移或瑞利商加速

---

例8.2 幂法计算主特征值与特征向量 详细解析

我将逐步骤拆解这道经典幂法例题，所有计算依据均以粗体标注，并深入分析收敛过程、误差来源与理论验证。

一、问题重述

用幂法计算对称矩阵

\[A = \begin{pmatrix} 1.0 & 1.0 & 0.5 \\ 1.0 & 1.0 & 0.25 \\ 0.5 & 0.25 & 2.0 \end{pmatrix} \]

的主特征值及其对应的特征向量。

二、幂法迭代过程详解

依据：规范化幂法迭代公式(8.5)

\[\begin{cases} u_0 = (1, 1, 1)^T \quad (\text{初始向量}) \\ v_k = A u_{k-1} \\ \mu_k = \max(v_k) = \|v_k\|_\infty \\ u_k = \frac{v_k}{\mu_k} \end{cases} \]

1. 第1次迭代（k=1）

• 计算 \(v_1 = A u_0\)：
  \[v_1 = \begin{pmatrix} 1.0 & 1.0 & 0.5 \\ 1.0 & 1.0 & 0.25 \\ 0.5 & 0.25 & 2.0 \end{pmatrix} \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix} = \begin{pmatrix} 1+1+0.5 \\ 1+1+0.25 \\ 0.5+0.25+2 \end{pmatrix} = \begin{pmatrix} 2.5 \\ 2.25 \\ 2.75 \end{pmatrix} \]

• 计算无穷范数 \(\mu_1 = \max(2.5, 2.25, 2.75) = 2.75\)
• 规范化得到 \(u_1\)：
  \[u_1 = \frac{v_1}{2.75} = \begin{pmatrix} 2.5/2.75 \\ 2.25/2.75 \\ 2.75/2.75 \end{pmatrix} \approx \begin{pmatrix} 0.9090909 \\ 0.8181818 \\ 1.0 \end{pmatrix} \]
  与教材中舍入值 \((0.9091, 0.8182, 1)\) 一致。

2. 第2次迭代（k=2，教材未给出，补充计算）

• 计算 \(v_2 = A u_1\)：
  \[v_2 = \begin{pmatrix} 1.0 & 1.0 & 0.5 \\ 1.0 & 1.0 & 0.25 \\ 0.5 & 0.25 & 2.0 \end{pmatrix} \begin{pmatrix} 0.9090909 \\ 0.8181818 \\ 1.0 \end{pmatrix} = \begin{pmatrix} 0.9090909+0.8181818+0.5 \\ 0.9090909+0.8181818+0.25 \\ 0.4545455+0.2045455+2.0 \end{pmatrix} \approx \begin{pmatrix} 2.2272727 \\ 1.9772727 \\ 2.6590909 \end{pmatrix} \]

• 计算无穷范数 \(\mu_2 = \max(2.2272727, 1.9772727, 2.6590909) \approx 2.6590909\)
• 规范化得到 \(u_2\)：
  \[u_2 \approx \begin{pmatrix} 2.2272727/2.6590909 \\ 1.9772727/2.6590909 \\ 1.0 \end{pmatrix} \approx \begin{pmatrix} 0.8376 \\ 0.7436 \\ 1.0 \end{pmatrix} \]

3. 完整迭代结果表

教材中给出的8位浮点运算结果如下：

k	规范化向量 \(u_k^T\)	\(\mu_k = \max\{v_k\}\)
0	(1, 1, 1)	-
1	(0.9091, 0.8182, 1)	2.7500000
5	(0.7651, 0.6674, 1)	2.5587919
10	(0.7494, 0.6508, 1)	2.5380032
15	(0.7483, 0.6497, 1)	2.5366257
16	(0.7483, 0.6497, 1)	2.5365841
17	(0.7482, 0.6497, 1)	2.5365598
18	(0.7482, 0.6497, 1)	2.5365457
19	(0.7482, 0.6497, 1)	2.5365374
20	(0.7482, 0.6497, 1)	2.5365326

三、收敛性分析

1. 收敛趋势观察

• 特征值收敛：\(\mu_k\) 从初始的2.75逐渐下降，前5步下降较快，第10步后收敛速度明显变慢，第18步后基本稳定在2.5365附近。
• 特征向量收敛：\(u_k\) 的分量逐渐稳定，第15步后前两个分量的变化已小于 \(10^{-4}\)，第18步后完全稳定。

2. 理论收敛速度验证

该矩阵的精确特征值为：

• \(\lambda_1 \approx 2.5365259\)（主特征值）
• \(\lambda_2 \approx 1.4801215\)（次大特征值）
• \(\lambda_3 \approx -0.0166474\)（最小特征值）

依据：幂法收敛速度公式，收敛速度由比值

\[r = \left| \frac{\lambda_2}{\lambda_1} \right| \approx \frac{1.4801215}{2.5365259} \approx 0.5835 \]

决定。

这个比值意味着每迭代一次，误差大约缩小到原来的58%。经过20次迭代后，误差约为初始误差的 \(0.5835^{20} \approx 2.7 \times 10^{-5}\)，与教材中计算结果的精度一致。

3. 特殊现象解释

在整个迭代过程中，\(u_k\) 的第三个分量始终为1，这是因为：

• 主特征向量 \(x_1 = (0.74822115, 0.64966114, 1)^T\) 的第三个分量是最大分量
• 我们采用无穷范数规范化，每次都除以最大分量，因此规范化后第三个分量恒为1
• 这一现象简化了计算和收敛判断，是幂法应用中的常见情况

四、结果验证与误差分析

1. 最终结果

• 计算得到的主特征值：\(\lambda_1 \approx 2.5365326\)
• 计算得到的主特征向量：\(x_1 \approx (0.7482, 0.6497, 1)^T\)

2. 与精确值对比

• 主特征值精确值：\(\lambda_1 = 2.5365259\)
• 主特征向量精确值：\(x_1 = (0.74822115, 0.64966114, 1)^T\)

3. 误差来源

1. 舍入误差：教材中使用8位浮点数字进行运算，每一步计算都会引入舍入误差，最终结果的误差约为 \(10^{-6}\) 量级。
2. 截断误差：迭代在第20步终止，没有达到理论上的完全收敛，剩余误差约为 \(|\mu_{20} - \lambda_1| \approx 6.7 \times 10^{-6}\)。

五、例题核心要点总结

1. 适用性验证：该矩阵是对称矩阵，所有特征值均为实数，且主特征值按模严格最大（\(|\lambda_1| > |\lambda_2| > |\lambda_3|\)），完全满足幂法的收敛条件。
2. 初始向量选择：选择全1向量作为初始向量，简单且有效，避免了初始向量与主特征向量正交的极端情况。
3. 收敛速度：收敛速度由次大特征值与主特征值的比值决定，本例中比值约为0.58，属于中等收敛速度，需要约20次迭代达到 \(10^{-6}\) 精度。
4. 规范化的重要性：通过无穷范数规范化，有效避免了数值溢出和有效数字损失，保证了迭代过程的数值稳定性。
5. 终止条件：当相邻两次迭代的特征值差小于预设精度（如 \(10^{-6}\)）时，即可终止迭代，本例中第18步后已满足该条件。

---

8.2.2 幂法的加速方法 详细讲解与证明

我将继续以资深数值分析研究员的身份，系统讲解两种最常用的幂法加速技术——原点平移法和瑞利商加速法，所有关键证明依据均以粗体标注，并深入分析其加速原理、参数选择技巧与适用场景。

一、原点平移法

1. 加速动机

幂法的收敛速度由比值 \(r = \left| \frac{\lambda_2}{\lambda_1} \right|\) 决定。当 \(r\) 接近1时，收敛会变得极其缓慢。例如，若 \(r=0.99\)，则需要约460次迭代才能将误差缩小到原来的1%，这在实际计算中是不可接受的。

原点平移法是一种简单而有效的加速技术，其核心思想是通过对矩阵进行平移变换，改变特征值的分布，从而显著提高收敛速度。

2. 基本原理

引入平移参数 \(p\)，构造新矩阵：

\[B = A - pI \]

依据：特征值的运算性质（定理8.1(2)）
若 \(\lambda\) 是 \(A\) 的特征值，\(x\) 是对应的特征向量，则 \(\lambda - p\) 是 \(B\) 的特征值，且特征向量仍为 \(x\)。

这意味着：

• \(A\) 和 \(B\) 有完全相同的特征向量
• \(B\) 的特征值是 \(A\) 的特征值减去平移量 \(p\)

如果我们能选择合适的 \(p\)，使得 \(B\) 的主特征值对应的收敛比值：

\[r' = \left| \frac{\lambda_2 - p}{\lambda_1 - p} \right| < \left| \frac{\lambda_2}{\lambda_1} \right| = r \]

那么对 \(B\) 应用幂法的收敛速度将远快于对 \(A\) 直接应用幂法。求出 \(B\) 的主特征值 \(\mu_1 = \lambda_1 - p\) 后，即可得到 \(A\) 的主特征值 \(\lambda_1 = \mu_1 + p\)。

3. 最优平移参数的选择

当 \(A\) 的所有特征值都是实数，且按降序排列为：

\[\lambda_1 > \lambda_2 \geq \dots \geq \lambda_{n-1} > \lambda_n \tag{8.6} \]

时，\(B = A - pI\) 的主特征值只能是以下两者之一：

• \(\mu_1 = \lambda_1 - p\)（当 \(p < \frac{\lambda_1 + \lambda_n}{2}\) 时）
• \(\mu_n = \lambda_n - p\)（当 \(p > \frac{\lambda_1 + \lambda_n}{2}\) 时）

计算主特征值 \(\lambda_1\) 的最优 \(p\)

若我们的目标是计算 \(\lambda_1\)，则应选择 \(p\) 使得：

1. \(|\lambda_1 - p| > |\lambda_n - p|\)（保证 \(\lambda_1 - p\) 是 \(B\) 的主特征值）
2. 收敛速度的比值 \(\omega = \max\left( \left| \frac{\lambda_2 - p}{\lambda_1 - p} \right|, \left| \frac{\lambda_n - p}{\lambda_1 - p} \right| \right)\) 最小

最优参数推导：
要使 \(\omega\) 最小，需让两个比值的绝对值相等：

\[\frac{\lambda_2 - p}{\lambda_1 - p} = -\frac{\lambda_n - p}{\lambda_1 - p} \]

两边同乘 \(\lambda_1 - p\)（正数）得：

\[\lambda_2 - p = -(\lambda_n - p) \]

解得最优平移参数：

\[p^* = \frac{\lambda_2 + \lambda_n}{2} \]

此时最小收敛比值为：

\[\omega_{\text{min}} = \frac{\lambda_2 - \lambda_n}{2\lambda_1 - \lambda_2 - \lambda_n} \]

计算最小特征值 \(\lambda_n\) 的最优 \(p\)

若目标是计算按模最小的特征值 \(\lambda_n\)，则最优平移参数为：

\[p^* = \frac{\lambda_1 + \lambda_{n-1}}{2} \]

4. 例题8.3 加速效果演示

设 \(A \in \mathbb{R}^{4 \times 4}\) 的特征值为：

\[\lambda_j = 15 - j, \quad j=1,2,3,4 \]

即 \(\lambda_1=14, \lambda_2=13, \lambda_3=12, \lambda_4=11\)。

• 原幂法收敛速度：\(r = \frac{\lambda_2}{\lambda_1} = \frac{13}{14} \approx 0.9\)，收敛很慢。
• 选择平移参数：\(p=12\)
• 变换后矩阵 \(B=A-12I\) 的特征值：\(\mu_1=2, \mu_2=1, \mu_3=0, \mu_4=-1\)
• 加速后收敛速度：\(r' = \left| \frac{\mu_2}{\mu_1} \right| = \frac{1}{2} = 0.5 < 0.9\)

加速效果对比：

• 原幂法：误差缩小到原来的1%需要约 \(\log_{0.9} 0.01 \approx 44\) 次迭代
• 加速后：误差缩小到原来的1%需要约 \(\log_{0.5} 0.01 \approx 7\) 次迭代
• 收敛速度提高了约6倍！

5. 例题8.4 实际应用

用原点平移法加速计算例8.2中矩阵

\[A = \begin{pmatrix} 1.0 & 1.0 & 0.5 \\ 1.0 & 1.0 & 0.25 \\ 0.5 & 0.25 & 2.0 \end{pmatrix} \]

的主特征值。

步骤1：估计平移参数 \(p\)

已知矩阵的迹 \(\operatorname{tr}A = 1.0 + 1.0 + 2.0 = 4\)，且由例8.2已得到主特征值的近似值 \(\lambda_1 \approx 2.5365259\)。

依据：矩阵迹的性质，所有特征值之和等于迹，故：

\[\lambda_2 + \lambda_3 = \operatorname{tr}A - \lambda_1 \approx 4 - 2.5365259 = 1.4634741 \]

根据最优参数公式：

\[p^* = \frac{\lambda_2 + \lambda_3}{2} \approx \frac{1.4634741}{2} \approx 0.7317 \]

教材中取近似值 \(p=0.75\)。

步骤2：构造变换矩阵 \(B\)

\[B = A - 0.75I = \begin{pmatrix} 0.25 & 1.0 & 0.5 \\ 1.0 & 0.25 & 0.25 \\ 0.5 & 0.25 & 1.25 \end{pmatrix} \]

步骤3：对 \(B\) 应用幂法

迭代结果如下表：

k	规范化向量 \(u_k^T\)	\(\mu_k = \max\{v_k\}\)
0	(1, 1, 1)	-
5	(0.7516, 0.6522, 1)	1.7914013
6	(0.7491, 0.6511, 1)	1.7888444
7	(0.7488, 0.6501, 1)	1.7873301
8	(0.7484, 0.6499, 1)	1.7869153
9	(0.7483, 0.6497, 1)	1.7866589
10	(0.7482, 0.6497, 1)	1.7865914

步骤4：还原得到 \(A\) 的主特征值

\[\lambda_1 \approx \mu_{10} + p = 1.7865914 + 0.75 = 2.5365914 \]

结果对比：

• 原幂法：20次迭代得到 \(\lambda_1 \approx 2.5365326\)，误差约 \(6.7 \times 10^{-6}\)
• 加速后：10次迭代得到 \(\lambda_1 \approx 2.5365914\)，误差约 \(6.5 \times 10^{-5}\)
• 若迭代15次，加速后得到 \(\lambda_1 \approx 2.5365277\)，误差仅约 \(1.8 \times 10^{-6}\)，精度超过原幂法20次迭代的结果

6. 原点平移法的优缺点

• 优点：
  1. 算法简单，仅需对矩阵做一次平移变换
  2. 不破坏矩阵的稀疏性，特别适合大型稀疏矩阵
  3. 加速效果显著，可将收敛速度提高数倍甚至数十倍
• 缺点：
  1. 需要对矩阵的特征值分布有大致了解，才能选择合适的平移参数 \(p\)
  2. 自动选择最优参数 \(p\) 的过程比较困难

二、瑞利商加速法

瑞利商加速法是专门针对对称矩阵设计的加速技术，它利用对称矩阵的特殊性质，在不增加太多计算量的情况下，将幂法的收敛速度提高一倍。

1. 瑞利商的定义与性质

回顾定理8.3，对称矩阵 \(A\) 的瑞利商定义为：

\[R(x) = \frac{(Ax, x)}{(x, x)}, \quad x \neq 0 \]

瑞利商具有以下重要性质：

• 对任意非零向量 \(x\)，有 \(\lambda_n \leq R(x) \leq \lambda_1\)
• 当 \(x\) 是 \(A\) 的特征向量时，\(R(x)\) 等于对应的特征值

2. 瑞利商加速定理（定理8.7）

设 \(A \in \mathbb{R}^{n \times n}\) 为对称矩阵，特征值满足：

\[|\lambda_1| > |\lambda_2| \geq |\lambda_3| \geq \dots \geq |\lambda_n| \]

对应的特征向量满足 \((x_i, x_j) = \delta_{ij}\)（标准正交）。应用幂法得到规范化向量序列 \(\{u_k\}\)，则 \(u_k\) 的瑞利商给出 \(\lambda_1\) 的更好近似：

\[R(u_k) = \frac{(A u_k, u_k)}{(u_k, u_k)} = \lambda_1 + O\left( \left( \frac{\lambda_2}{\lambda_1} \right)^{2k} \right) \]

对比：原幂法中 \(\mu_k = \max(v_k)\) 的收敛阶为 \(O\left( \left( \frac{\lambda_2}{\lambda_1} \right)^k \right)\)，瑞利商加速将收敛速度提高了一倍！

3. 详细证明

依据：幂法的迭代公式

\[u_k = \frac{A^k u_0}{\max(A^k u_0)} \]

将初始向量 \(u_0\) 表示为标准正交特征向量的线性组合：

\[u_0 = \sum_{i=1}^n \alpha_i x_i, \quad \alpha_1 \neq 0 \]

则：

\[A^k u_0 = \sum_{i=1}^n \alpha_i \lambda_i^k x_i \]

\[A^{k+1} u_0 = \sum_{i=1}^n \alpha_i \lambda_i^{k+1} x_i \]

计算瑞利商：

\[\begin{align*} R(u_k) &= \frac{(A u_k, u_k)}{(u_k, u_k)} \\ &= \frac{(A \cdot \frac{A^k u_0}{\max(A^k u_0)}, \frac{A^k u_0}{\max(A^k u_0)})}{(\frac{A^k u_0}{\max(A^k u_0)}, \frac{A^k u_0}{\max(A^k u_0)})} \\ &= \frac{(A^{k+1} u_0, A^k u_0)}{(A^k u_0, A^k u_0)} \end{align*} \]

依据：对称矩阵的特征向量正交性，内积中交叉项均为零，故：

\[(A^{k+1} u_0, A^k u_0) = \left( \sum_{i=1}^n \alpha_i \lambda_i^{k+1} x_i, \sum_{j=1}^n \alpha_j \lambda_j^k x_j \right) = \sum_{i=1}^n \alpha_i^2 \lambda_i^{2k+1} \]

\[(A^k u_0, A^k u_0) = \sum_{i=1}^n \alpha_i^2 \lambda_i^{2k} \]

因此：

\[R(u_k) = \frac{\sum_{i=1}^n \alpha_i^2 \lambda_i^{2k+1}}{\sum_{i=1}^n \alpha_i^2 \lambda_i^{2k}} \]

提取分子分母的公因子 \(\lambda_1^{2k}\)：

\[R(u_k) = \lambda_1 \cdot \frac{\alpha_1^2 + \sum_{i=2}^n \alpha_i^2 \left( \frac{\lambda_i}{\lambda_1} \right)^{2k+1}}{\alpha_1^2 + \sum_{i=2}^n \alpha_i^2 \left( \frac{\lambda_i}{\lambda_1} \right)^{2k}} \]

当 \(k \to \infty\) 时，由于 \(|\lambda_i/\lambda_1| < 1\)（\(i \geq 2\)），高阶小项趋于零，展开得：

\[R(u_k) = \lambda_1 + O\left( \left( \frac{\lambda_2}{\lambda_1} \right)^{2k} \right) \]

证毕。

4. 瑞利商加速的实现

在幂法的每一步迭代中，只需增加一次内积计算即可得到更精确的特征值估计：

\[\begin{cases} u_0 = v_0 \neq 0 \\ v_k = A u_{k-1} \\ \mu_k = \frac{(v_k, u_{k-1})}{(u_{k-1}, u_{k-1})} \quad (\text{瑞利商}) \\ u_k = \frac{v_k}{\max(v_k)} \end{cases} \]

计算量分析：每步迭代增加了两次内积计算（\(O(n)\) 操作），与矩阵-向量乘法（\(O(n^2)\) 操作）相比，计算量增加可以忽略不计，但收敛速度提高了一倍。

三、知识点系统总结表

加速方法	核心原理	适用条件	收敛速度	优点	缺点
原点平移法	构造 \(B=A-pI\)，改变特征值分布，减小收敛比值	所有可对角化矩阵	\(O\left( \left| \frac{\lambda_2-p}{\lambda_1-p} \right|^k \right)\)	算法简单，不破坏稀疏性，加速效果显著	需要对特征值分布有大致了解，自动选参困难
瑞利商加速法	利用对称矩阵瑞利商的性质，将收敛阶提高一倍	仅适用于对称矩阵	\(O\left( \left( \frac{\lambda_2}{\lambda_1} \right)^{2k} \right)\)	无需额外参数，收敛速度加倍，计算量增加极少	仅适用于对称矩阵

重要说明：两种加速方法可以结合使用。对于对称矩阵，先应用原点平移法选择合适的 \(p\)，再对平移后的矩阵应用瑞利商加速，可以获得极其优异的收敛效果。

---

8.2.3 反幂法 详细讲解与证明

我将继续以资深数值分析研究员的身份，系统讲解求解单个高精度特征值的"黄金标准"——反幂法，所有关键证明依据均以粗体标注，深入剖析其收敛原理、实用算法设计与工业级应用技巧。

一、反幂法的基本概念与核心思想

1. 定义与适用场景

反幂法是一种计算矩阵按模最小的特征值及其对应特征向量的迭代方法，更重要的是，它可以用来计算对应于一个给定近似特征值的精确特征向量。

反幂法是目前已知求解单个高精度特征值和特征向量最有效的方法，广泛应用于结构振动分析、量子力学、控制系统设计等领域。

2. 核心原理

依据：逆矩阵的特征值性质（定理8.1(2)）
设 \(A \in \mathbb{R}^{n \times n}\) 为非奇异矩阵，其特征值按模降序排列为：

\[|\lambda_1| \geq |\lambda_2| \geq \dots \geq |\lambda_n| > 0 \]

对应的特征向量为 \(x_1, x_2, \dots, x_n\)。则逆矩阵 \(A^{-1}\) 的特征值为：

\[\frac{1}{\lambda_n}, \frac{1}{\lambda_{n-1}}, \dots, \frac{1}{\lambda_1} \]

且满足：

\[\left| \frac{1}{\lambda_n} \right| \geq \left| \frac{1}{\lambda_{n-1}} \right| \geq \dots \geq \left| \frac{1}{\lambda_1} \right| \]

对应的特征向量与 \(A\) 完全相同，仍为 \(x_n, x_{n-1}, \dots, x_1\)。

这意味着：矩阵 \(A\) 的按模最小特征值 \(\lambda_n\)，恰好对应于逆矩阵 \(A^{-1}\) 的按模最大特征值 \(1/\lambda_n\)。因此，我们可以对 \(A^{-1}\) 应用幂法，求出 \(1/\lambda_n\)，再取倒数得到 \(\lambda_n\)，这就是反幂法的基本思想。

二、反幂法的迭代公式与收敛性

1. 基本迭代公式

对 \(A^{-1}\) 应用规范化幂法，得到反幂法的迭代公式：

\[\begin{cases} u_0 = v_0 \neq 0 \quad (\text{初始向量}) \\ v_k = A^{-1} u_{k-1} \\ \mu_k = \max(v_k) = \|v_k\|_\infty \\ u_k = \frac{v_k}{\mu_k} \end{cases} \quad k=1,2,\dots \]

关键说明：在实际计算中，我们从不直接计算逆矩阵 \(A^{-1}\)，因为这既昂贵又不稳定。迭代式 \(v_k = A^{-1} u_{k-1}\) 等价于解线性方程组：

\[A v_k = u_{k-1} \]

这是反幂法实现的核心技巧。

2. 收敛性定理

定理8.8 设 \(A \in \mathbb{R}^{n \times n}\) 有 \(n\) 个线性无关的特征向量，其特征值满足：

\[|\lambda_1| \geq |\lambda_2| \geq \dots \geq |\lambda_n| > 0 \]

若非零初始向量 \(u_0 = \sum_{i=1}^n \alpha_i x_i\) 中的系数 \(\alpha_n \neq 0\)，则反幂法产生的序列满足：

1. \(\displaystyle \lim_{k \to \infty} u_k = \frac{x_n}{\max(x_n)}\)（收敛到按模最小特征值对应的规范化特征向量）
2. \(\displaystyle \lim_{k \to \infty} \mu_k = \frac{1}{\lambda_n}\)，即 \(\displaystyle \lim_{k \to \infty} \frac{1}{\mu_k} = \lambda_n\)（收敛到按模最小特征值）

且收敛速度由比值：

\[r = \left| \frac{\lambda_n}{\lambda_{n-1}} \right| \]

确定。\(r\) 越小，收敛越快。

证明：与幂法收敛定理（定理8.6）的证明完全类似，只需将 \(A\) 替换为 \(A^{-1}\)，特征值 \(\lambda_i\) 替换为 \(1/\lambda_i\) 即可。

三、带原点平移的反幂法（工业标准版本）

基本反幂法只能求解按模最小的特征值，这大大限制了其应用范围。将原点平移法与反幂法结合，得到的带原点平移的反幂法可以求解矩阵任意指定位置的单个特征值，这是目前工程应用中最常用的版本。

1. 基本原理

设 \(p\) 是目标特征值 \(\lambda_j\) 的一个近似值，且矩阵 \(A-pI\) 非奇异。

依据：特征值的运算性质
\(A-pI\) 的特征值为 \(\lambda_1-p, \lambda_2-p, \dots, \lambda_n-p\)，对应的特征向量仍为 \(x_1, x_2, \dots, x_n\)。

因此，\((A-pI)^{-1}\) 的特征值为：

\[\frac{1}{\lambda_1-p}, \frac{1}{\lambda_2-p}, \dots, \frac{1}{\lambda_n-p} \]

如果 \(p\) 非常接近 \(\lambda_j\)，且 \(\lambda_j\) 与其他特征值分离，即：

\[|\lambda_j - p| \ll |\lambda_i - p|, \quad \forall i \neq j \]

那么：

\[\left| \frac{1}{\lambda_j - p} \right| \gg \left| \frac{1}{\lambda_i - p} \right|, \quad \forall i \neq j \]

即 \(1/(\lambda_j-p)\) 是 \((A-pI)^{-1}\) 的绝对占优的主特征值。

此时对 \((A-pI)^{-1}\) 应用反幂法，将以极快的速度收敛到特征向量 \(x_j\)，并得到特征值的精确估计。

2. 迭代公式

带原点平移的反幂法迭代公式为：

\[\begin{cases} u_0 = v_0 \neq 0 \\ v_k = (A-pI)^{-1} u_{k-1} \\ \mu_k = \max(v_k) \\ u_k = \frac{v_k}{\mu_k} \end{cases} \tag{8.7} \]

同样，迭代式 \(v_k = (A-pI)^{-1} u_{k-1}\) 等价于解线性方程组：

\[(A-pI) v_k = u_{k-1} \]

3. 收敛性与收敛速度

定理8.8（推广） 设 \(p\) 是 \(\lambda_j\) 的近似值，\(A-pI\) 非奇异，且：

\[|\lambda_j - p| \ll |\lambda_i - p|, \quad \forall i \neq j \]

初始向量 \(u_0\) 中 \(\alpha_j \neq 0\)，则：

1. \(\displaystyle \lim_{k \to \infty} u_k = \frac{x_j}{\max(x_j)}\)
2. \(\displaystyle \lim_{k \to \infty} \left( p + \frac{1}{\mu_k} \right) = \lambda_j\)

收敛速度由比值：

\[r = \frac{|\lambda_j - p|}{\min_{i \neq j} |\lambda_i - p|} \]

确定。

惊人的收敛速度：当 \(p\) 与 \(\lambda_j\) 的误差为 \(\varepsilon\)，且最近的其他特征值距离为 \(\delta\) 时，\(r \approx \varepsilon/\delta\)。如果 \(\varepsilon = 10^{-3}\)，\(\delta = 1\)，则 \(r = 10^{-3}\)，每次迭代误差缩小1000倍，通常只需2-3次迭代即可达到机器精度！这是所有特征值迭代方法中最快的收敛速度。

四、反幂法的实用算法实现

为了提高计算效率，我们在迭代开始前先对矩阵 \(A-pI\) 进行LU分解，这样每步迭代只需解两个三角形方程组，复杂度从 \(O(n^3)\) 降为 \(O(n^2)\)。

完整算法步骤

1. 预处理：LU分解
   对矩阵 \(A-pI\) 进行带部分选主元的LU分解：

   \[P(A-pI) = LU \]

   其中 \(P\) 是排列矩阵，\(L\) 是单位下三角矩阵，\(U\) 是上三角矩阵。保存 \(L, U, P\) 供后续迭代使用。

2. 初始向量选择（教材推荐最优方法）
   解方程组 \(U v_1 = L^{-1} P u_0 = (1,1,\dots,1)^T\)，得到初始迭代向量 \(v_1\)。这种方法避免了初始向量与目标特征向量正交的问题，数值稳定性最好。

3. 迭代过程
   对于 \(k=2,3,\dots\)：

   1. 解下三角方程组 \(L y_k = P u_{k-1}\)，得到 \(y_k\)
   2. 解上三角方程组 \(U v_k = y_k\)，得到 \(v_k\)
   3. 计算 \(\mu_k = \max(v_k) = \|v_k\|_\infty\)
   4. 规范化得到 \(u_k = v_k / \mu_k\)
   5. 检查收敛性：若 \(|\mu_k - \mu_{k-1}| < \varepsilon\) 或 \(\|u_k - u_{k-1}\|_\infty < \varepsilon\)，则终止迭代

4. 结果计算

   • 特征向量：\(x_j \approx u_k\)
   • 特征值：\(\lambda_j \approx p + 1/\mu_k\)

五、例题8.5 详细解析

用反幂法求矩阵

\[A = \begin{pmatrix} 2 & 1 & 0 \\ 1 & 3 & 1 \\ 0 & 1 & 4 \end{pmatrix} \]

对应于近似特征值 \(\lambda \approx 1.2679\)（精确值为 \(\lambda_3 = 3-\sqrt{3} \approx 1.26794919\)）的特征向量。

步骤1：构造平移矩阵并进行LU分解

取平移参数 \(p=1.2679\)，构造矩阵：

\[A-pI = \begin{pmatrix} 2-1.2679 & 1 & 0 \\ 1 & 3-1.2679 & 1 \\ 0 & 1 & 4-1.2679 \end{pmatrix} = \begin{pmatrix} 0.7321 & 1 & 0 \\ 1 & 1.7321 & 1 \\ 0 & 1 & 2.7321 \end{pmatrix} \]

对其进行带部分选主元的LU分解，得到：

\[P = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix}, \quad L = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0.7321 & -0.26807 & 1 \end{pmatrix}, \quad U = \begin{pmatrix} 1 & 1.7321 & 1 \\ 0 & 1 & 2.7321 \\ 0 & 0 & 0.29517 \times 10^{-3} \end{pmatrix} \]

关键观察：\(U\) 的最后一个对角元非常小（约 \(3 \times 10^{-4}\)），这说明 \(A-pI\) 非常接近奇异，即 \(p\) 非常接近真实特征值，预示着反幂法将以极快的速度收敛。

步骤2：第一次迭代（k=1）

解方程组 \(U v_1 = (1,1,1)^T\)：

\[\begin{cases} v_{11} + 1.7321 v_{12} + v_{13} = 1 \\ v_{12} + 2.7321 v_{13} = 1 \\ 0.29517 \times 10^{-3} v_{13} = 1 \end{cases} \]

从下往上回代求解：

• \(v_{13} = 1 / (0.29517 \times 10^{-3}) \approx 3387.9\)
• \(v_{12} = 1 - 2.7321 \times 3387.9 \approx -9250.3\)
• \(v_{11} = 1 - 1.7321 \times (-9250.3) - 3387.9 \approx 12632\)

因此：

\[v_1 \approx (12632, -9250.3, 3387.9)^T \]

\[\mu_1 = \max(v_1) \approx 12632 \]

\[u_1 = v_1 / \mu_1 \approx (1, -0.7323, 0.2682)^T \]

步骤3：第二次迭代（k=2）

1. 解下三角方程组 \(L y_2 = P u_1\)：

   \[P u_1 = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} \begin{pmatrix} 1 \\ -0.7323 \\ 0.2682 \end{pmatrix} = \begin{pmatrix} -0.7323 \\ 0.2682 \\ 1 \end{pmatrix} \]
   \[\begin{cases} y_{21} = -0.7323 \\ y_{22} = 0.2682 \\ 0.7321 y_{21} - 0.26807 y_{22} + y_{23} = 1 \end{cases} \]

   解得 \(y_2 \approx (-0.7323, 0.2682, 1.608)^T\)

2. 解上三角方程组 \(U v_2 = y_2\)：

   \[\begin{cases} v_{21} + 1.7321 v_{22} + v_{23} = -0.7323 \\ v_{22} + 2.7321 v_{23} = 0.2682 \\ 0.29517 \times 10^{-3} v_{23} = 1.608 \end{cases} \]

   解得 \(v_2 \approx (20404, -14937, 5467.4)^T\)

3. 规范化：

   \[\mu_2 = \max(v_2) \approx 20404 \]
   \[u_2 = v_2 / \mu_2 \approx (1, -0.73206, 0.26796)^T \]

步骤4：结果计算与验证

• 计算得到的特征向量：\(x_3 \approx (1, -0.73206, 0.26796)^T\)
• 精确特征向量：\(x_3 = (1, 1-\sqrt{3}, 2-\sqrt{3})^T \approx (1, -0.73205, 0.26795)^T\)
• 计算得到的特征值：\(\lambda_3 \approx p + 1/\mu_2 = 1.2679 + 1/20404 \approx 1.26794901\)
• 精确特征值：\(\lambda_3 = 3-\sqrt{3} \approx 1.26794919\)

惊人的精度：仅经过2次迭代，特征值的误差就小于 \(2 \times 10^{-7}\)，特征向量的误差小于 \(1 \times 10^{-5}\)，充分展示了反幂法的超强收敛能力。

六、知识点系统总结表

知识点类别	核心内容	理论依据	关键备注与应用
基本反幂法	求解按模最小的特征值及对应特征向量	逆矩阵的特征值性质、幂法收敛定理	从不直接求逆，而是解线性方程组
带原点平移的反幂法	求解任意指定近似值附近的特征值及对应特征向量	特征值的平移性质、反幂法收敛定理	工业标准版本，应用最广泛
收敛速度	收敛速度由比值 \(r=|\lambda_j-p|/\min_{i≠j}|\lambda_i-p|\) 决定	幂法收敛速度分析	当p接近λⱼ时，r趋近于0，通常只需2-3次迭代
实用算法	先对A-pI做LU分解，每步迭代解两个三角形方程组	LU分解的性质、三角形方程组解法	预处理复杂度O(n³)，每步迭代复杂度O(n²)
初始向量选择	推荐解Uv₁=(1,1,…,1)ᵀ得到初始向量	数值稳定性分析	避免初始向量与目标特征向量正交
优势	精度最高，收敛最快，可求解任意位置的单个特征值	理论收敛性证明	是目前求解单个高精度特征值的首选方法
局限性	每次只能求解一个特征值，需要特征值的近似值	算法本质	通常与其他方法（如格什戈林定理、QR方法）结合使用

---

8.3 矩阵分解与正交变换（一）：基本QR方法 详细讲解与证明

我将继续以资深数值分析研究员的身份，系统讲解求解矩阵全部特征值的工业标准算法——QR方法，所有关键证明依据均以粗体标注，深入剖析其理论基础、迭代原理与收敛特性。

一、QR分解的背景与定义

1. 从LR算法到QR算法

求解矩阵全部特征值的早期方法是1958年Rutishauser提出的LR算法，其基本思想是：

• 对矩阵 \(A\) 进行LU分解：\(A = LU\)
• 交换因子顺序得到新矩阵：\(B = UL = L^{-1}AL\)
• \(B\) 与 \(A\) 相似，具有相同的特征值
• 重复上述过程，矩阵序列将逐渐收敛到上三角矩阵，对角元即为特征值

LR算法虽然原理简单，但存在严重的数值稳定性问题，因为LU分解中的选主元操作会破坏相似性，导致算法失效。为了解决这个问题，1961年Francis提出了用正交变换代替LU分解的QR方法，它具有极佳的数值稳定性，成为目前求解中小型稠密矩阵全部特征值的最有效方法。

2. 格拉姆-施密特正交化过程

QR分解的基础是格拉姆-施密特(Gram-Schmidt)正交化，它可以将一组线性无关的向量转化为正交向量组，再单位化得到标准正交向量组。

设矩阵 \(A \in \mathbb{R}^{n \times n}\) 非奇异，其列向量 \(a_1, a_2, \dots, a_n\) 线性无关。格拉姆-施密特正交化过程如下：

1. 取第一个正交向量：\(\eta_1 = a_1\)
2. 对 \(k=2,3,\dots,n\)，将 \(a_k\) 减去它在前面所有正交向量上的投影：
   \[\eta_k = a_k - \sum_{i=1}^{k-1} \frac{(a_k, \eta_i)}{(\eta_i, \eta_i)} \eta_i \]

3. 单位化得到标准正交向量：
   \[q_i = \frac{\eta_i}{\|\eta_i\|}, \quad i=1,2,\dots,n \]

依据：正交投影定理，向量 \(a_k\) 减去它在子空间 \(\operatorname{span}\{\eta_1,\dots,\eta_{k-1}\}\) 上的正交投影后，得到的向量 \(\eta_k\) 与该子空间正交。

3. QR分解的矩阵形式

将上述正交化过程表示为矩阵相乘的形式：

\[\begin{align*} (a_1, a_2, \dots, a_n) &= (\eta_1, \eta_2, \dots, \eta_n) \begin{pmatrix} 1 & \frac{(a_2, \eta_1)}{(\eta_1, \eta_1)} & \dots & \frac{(a_n, \eta_1)}{(\eta_1, \eta_1)} \\ 0 & 1 & \dots & \frac{(a_n, \eta_2)}{(\eta_2, \eta_2)} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 1 \end{pmatrix} \\ &= \left( \frac{\eta_1}{\|\eta_1\|}, \frac{\eta_2}{\|\eta_2\|}, \dots, \frac{\eta_n}{\|\eta_n\|} \right) \begin{pmatrix} \|\eta_1\| & (a_2, q_1) & \dots & (a_n, q_1) \\ 0 & \|\eta_2\| & \dots & (a_n, q_2) \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \|\eta_n\| \end{pmatrix} \end{align*} \]

令 \(Q = (q_1, q_2, \dots, q_n)\)（正交矩阵，满足 \(Q^T Q = I\)），\(R\) 为上式右边的上三角矩阵，则得到：

\[A = QR \tag{8.9} \]

4. QR分解定理（定理8.9）

设 \(A \in \mathbb{R}^{n \times n}\) 非奇异，则存在正交矩阵 \(Q\) 和上三角矩阵 \(R\)，使得 \(A = QR\)。若要求 \(R\) 的对角元为正，则分解是唯一的。

证明：存在性由上述格拉姆-施密特正交化过程直接得到。唯一性证明略。

二、基本QR方法的迭代原理

1. 基本迭代公式

对于非奇异矩阵 \(A\)，基本QR方法的迭代过程如下：

1. 初始化：\(A_1 = A\)
2. 对 \(k=1,2,\dots\)：
   • 对 \(A_k\) 进行QR分解：\(A_k = Q_k R_k\)
   • 交换因子顺序构造新矩阵：\(A_{k+1} = R_k Q_k\)

2. 核心性质：相似不变性

定理：QR迭代产生的矩阵序列 \(\{A_k\}\) 中的所有矩阵都与原矩阵 \(A\) 相似，因此具有完全相同的特征值。

证明：

\[A_{k+1} = R_k Q_k = Q_k^T Q_k R_k Q_k = Q_k^T A_k Q_k \]

依据：相似矩阵的定义，若存在非奇异矩阵 \(P\) 使得 \(B = P^{-1}AP\)，则 \(B\) 与 \(A\) 相似。这里 \(Q_k\) 是正交矩阵，\(Q_k^{-1} = Q_k^T\)，因此 \(A_{k+1}\) 与 \(A_k\) 相似。由数学归纳法，所有 \(A_k\) 都与 \(A_1 = A\) 相似。

这是QR方法最关键的性质，它保证了迭代过程中特征值保持不变，我们的目标是通过迭代让 \(A_k\) 逐渐收敛到上三角矩阵，此时对角元就是原矩阵的特征值。

3. QR迭代的基本性质（定理8.10）

记 \(\bar{Q}_k = Q_1 Q_2 \dots Q_k\)（前 \(k\) 个正交矩阵的乘积，仍为正交矩阵），\(\bar{R}_k = R_k R_{k-1} \dots R_1\)（前 \(k\) 个上三角矩阵的乘积，仍为上三角矩阵），则有：

1. \(A_{k+1} = \bar{Q}_k^T A \bar{Q}_k\)
2. \(A^k = \bar{Q}_k \bar{R}_k\)

证明：

1. 由相似不变性递推可得：

   \[A_{k+1} = Q_k^T A_k Q_k = Q_k^T Q_{k-1}^T A_{k-1} Q_{k-1} Q_k = \dots = (Q_1 Q_2 \dots Q_k)^T A (Q_1 Q_2 \dots Q_k) = \bar{Q}_k^T A \bar{Q}_k \]

2. 依据：数学归纳法原理

   • 基例：\(k=1\) 时，\(A^1 = A = Q_1 R_1 = \bar{Q}_1 \bar{R}_1\)，成立。
   • 归纳假设：设 \(k=m\) 时，\(A^m = \bar{Q}_m \bar{R}_m\) 成立。
   • 归纳步骤：当 \(k=m+1\) 时：
     \[\begin{align*} A^{m+1} &= A \cdot A^m = A \cdot \bar{Q}_m \bar{R}_m \\ &= \bar{Q}_m (\bar{Q}_m^T A \bar{Q}_m) \bar{R}_m \\ &= \bar{Q}_m A_{m+1} \bar{R}_m \\ &= \bar{Q}_m Q_{m+1} R_{m+1} \bar{R}_m \\ &= (\bar{Q}_m Q_{m+1}) (R_{m+1} \bar{R}_m) \\ &= \bar{Q}_{m+1} \bar{R}_{m+1} \end{align*} \]

   由数学归纳法，对所有正整数 \(k\)，结论成立。

重要意义：性质2表明，\(\bar{Q}_k\) 的列向量是 \(A^k\) 的标准正交基，这揭示了QR方法与幂法之间的深刻联系——QR方法本质上是同时对多个初始向量进行幂法迭代。

三、基本QR方法的收敛性

1. 一般矩阵的收敛性（定理8.11）

设 \(A \in \mathbb{R}^{n \times n}\) 满足以下条件：

1. \(A\) 的特征值按模严格分离：\(|\lambda_1| > |\lambda_2| > \dots > |\lambda_n| > 0\)
2. \(A\) 有标准形 \(A = X D X^{-1}\)，其中 \(D = \operatorname{diag}(\lambda_1, \lambda_2, \dots, \lambda_n)\)，且 \(X^{-1}\) 有三角分解 \(X^{-1} = LU\)（\(L\) 为单位下三角矩阵，\(U\) 为上三角矩阵）

则由基本QR方法产生的矩阵序列 \(\{A_k\}\) 本质上收敛于上三角矩阵：

\[A_k \xrightarrow{\text{本质上}} R = \begin{pmatrix} \lambda_1 & * & \dots & * \\ & \lambda_2 & \dots & * \\ & & \ddots & \vdots \\ & & & \lambda_n \end{pmatrix}, \quad k \to \infty \]

即：

1. 对角元收敛：\(\displaystyle \lim_{k \to \infty} a_{ii}^{(k)} = \lambda_i\)，\(i=1,2,\dots,n\)
2. 下三角元收敛到零：\(\displaystyle \lim_{k \to \infty} a_{ij}^{(k)} = 0\)，\(i > j\)
3. 上三角元的极限不一定存在

收敛速度：下三角元 \(a_{ij}^{(k)}\)（\(i>j\)）的收敛速度由比值 \(|\lambda_i/\lambda_j|\) 决定，比值越小，收敛越快。

2. 对称矩阵的收敛性（定理8.12）

如果对称矩阵 \(A\) 满足定理8.11的条件，则由QR方法产生的矩阵序列 \(\{A_k\}\) 收敛于对角矩阵：

\[A_k \to D = \operatorname{diag}(\lambda_1, \lambda_2, \dots, \lambda_n), \quad k \to \infty \]

证明：对称矩阵经过正交相似变换后仍然是对称矩阵。如果 \(A_k\) 本质上收敛于上三角矩阵，同时又是对称矩阵，那么它必须收敛于对角矩阵。

这是对称矩阵QR方法的一个巨大优势，不仅收敛到对角矩阵，而且数值稳定性极佳。

3. 重特征值与复特征值的情况

当矩阵有重特征值或复特征值时，基本QR方法仍然收敛，但收敛到分块上三角矩阵：

\[A_k \to \begin{pmatrix} \lambda_1 & \dots & * & * & \dots & * \\ & \ddots & \vdots & \vdots & & \vdots \\ & & \lambda_m & * & \dots & * \\ & & & B_1 & \dots & * \\ & & & & \ddots & \vdots \\ & & & & & B_l \end{pmatrix} \]

其中 \(m + 2l = n\)，每个 \(B_i\)（\(i=1,2,\dots,l\)）是一个 \(2 \times 2\) 的子块，它的两个特征值恰好是原矩阵 \(A\) 的一对共轭复特征值。

这意味着QR方法可以自然地处理实矩阵的复特征值问题，无需引入复数运算。

四、格拉姆-施密特正交化的数值问题

教材中提到："可惜这个算法的数值效果不太好，所计算的 \(q_i\) 之间的正交性常常会严重损失。"

原因分析：经典格拉姆-施密特正交化在计算过程中，由于舍入误差的积累，会导致计算得到的向量 \(q_i\) 之间逐渐失去正交性。当矩阵的列向量接近线性相关时，这种正交性损失会特别严重，甚至导致QR分解完全失效。

解决方案：

1. 改进的格拉姆-施密特正交化：改变计算顺序，将向量逐个正交化，数值稳定性有显著提高。
2. Householder变换法：通过一系列正交反射变换将矩阵化为上三角矩阵，数值稳定性极佳，是目前QR分解的标准方法。
3. Givens旋转变换法：通过一系列正交旋转变换将矩阵化为上三角矩阵，特别适合处理稀疏矩阵。

五、知识点系统总结表

知识点类别	核心内容	理论依据	关键备注与应用
QR分解定义	非奇异矩阵可分解为正交矩阵Q与上三角矩阵R的乘积：A=QR	格拉姆-施密特正交化定理	若要求R对角元为正，则分解唯一
基本QR迭代	A_k = Q_k R_k，A_{k+1} = R_k Q_k	相似矩阵的定义	迭代过程保持特征值不变
迭代性质	A_{k+1} = \bar{Q}_k^T A \bar{Q}_k，A^k = \bar{Q}_k \bar{R}_k	数学归纳法、矩阵乘法结合律	揭示了QR方法与幂法的内在联系
一般收敛性	若特征值按模严格分离且X^{-1}有三角分解，则本质上收敛到上三角矩阵	幂法收敛性、正交变换性质	对角元收敛到特征值，下三角元收敛到零
对称矩阵收敛性	对称矩阵收敛到对角矩阵	对称矩阵的正交相似不变性	收敛更快更稳定，是对称矩阵特征值求解的首选
复特征值处理	收敛到分块上三角矩阵，2×2子块对应共轭复特征值	实矩阵复特征值成对出现	无需复数运算即可求解复特征值
数值问题	经典格拉姆-施密特正交化会导致正交性损失	舍入误差积累	实际应用中使用Householder变换或改进的格拉姆-施密特

---

例8.6 格拉姆-施密特正交化的数值稳定性分析 详细讲解

我将以资深数值分析研究员的身份，深入剖析这道揭示经典正交化方法致命缺陷的经典例题，所有关键理论依据均以粗体标注，并系统对比不同正交化方法的数值稳定性差异。

一、例题背景与测试矩阵选择

1. 问题提出

QR分解是QR方法的核心步骤，其数值稳定性直接决定了整个特征值求解过程的可靠性。格拉姆-施密特正交化是QR分解的基础方法，但教材中明确指出其"数值效果不太好"。本例题通过对希尔伯特矩阵进行QR分解，定量展示了经典格拉姆-施密特正交化的正交性损失问题。

2. 希尔伯特矩阵的特殊性

希尔伯特矩阵是数值分析中最著名的病态矩阵，其定义为：

\[H_n = \left( \frac{1}{i+j-1} \right)_{n \times n} \]

它具有以下关键性质：

• 对称正定
• 条件数随阶数 \(n\) 呈指数增长：\(\operatorname{cond}_2(H_n) \approx O(10^{1.5n})\)
• 列向量高度线性相关，是检验数值方法稳定性的"试金石"

依据：矩阵条件数的定义，条件数越大，矩阵越病态，数值计算中舍入误差的影响越严重。

3. 正交性检验方法

理论上，正交矩阵 \(Q\) 满足 \(Q^T Q = I\)，即非对角元应全为零。由于 \(Q^T Q\) 是对称矩阵，我们只需检查其严格下三角元素的绝对值最大值：

\[e = \max_{i>j} |e_{ij}| \]

作为正交性的检验指标。\(e\) 越小，说明计算得到的 \(Q\) 的正交性越好；\(e\) 越大，说明正交性损失越严重。

二、经典格拉姆-施密特正交化的数值表现

1. 实验结果（表8-3）

用经典格拉姆-施密特正交化方法对4~9阶希尔伯特矩阵进行QR分解，得到的正交性误差如下：

矩阵阶数 \(n\)	4	5	6	7	8	9
正交性误差 \(e\)	\(4.3964 \times 10^{-11}\)	\(6.6685 \times 10^{-8}\)	\(3.0277 \times 10^{-4}\)	\(8.9561 \times 10^{-2}\)	\(9.9997 \times 10^{-1}\)	\(9.9999 \times 10^{-1}\)

2. 结果分析

• 低阶情况（n≤4）：误差很小（\(10^{-11}\) 量级），正交性保持良好。
• 中阶情况（n=5~6）：误差急剧增大，从 \(10^{-8}\) 增长到 \(10^{-4}\)，正交性开始明显损失。
• 高阶情况（n≥7）：误差接近1，说明计算得到的"正交矩阵"实际上已经完全不正交，QR分解完全失效。

3. 正交性损失的根本原因

依据：经典格拉姆-施密特正交化的串行计算特性
经典方法的计算顺序是：

\[\eta_k = a_k - \sum_{i=1}^{k-1} \frac{(a_k, \eta_i)}{(\eta_i, \eta_i)} \eta_i \]

这种"一次性减去所有投影"的计算方式存在致命缺陷：

1. 误差累积效应：每一步计算产生的舍入误差会全部传递到后续所有步骤
2. 误差放大效应：对于病态矩阵，微小的舍入误差会被矩阵的高条件数急剧放大
3. 正交性丧失：当计算到第 \(k\) 个向量时，前面的向量 \(\eta_1, \dots, \eta_{k-1}\) 已经由于舍入误差失去了正交性，因此减去它们的投影并不能保证 \(\eta_k\) 与它们正交

对于希尔伯特矩阵，当 \(n=9\) 时，其条件数约为 \(10^{13}\)，这意味着即使初始误差只有 \(10^{-16}\)（双精度浮点数的机器精度），经过放大后也会达到 \(10^{-3}\) 量级，与实验结果完全一致。

三、改进的格拉姆-施密特正交化

1. 算法改进思路

为了解决经典方法的正交性损失问题，人们提出了改进的格拉姆-施密特正交化（Modified Gram-Schmidt, MGS）。它改变了计算顺序，将"一次性减去所有投影"改为"逐个向量正交化"。

改进算法步骤：

1. 初始化：\(\eta_k = a_k\)
2. 对 \(i=1,2,\dots,k-1\)：
   \[\eta_k = \eta_k - \frac{(\eta_k, \eta_i)}{(\eta_i, \eta_i)} \eta_i \]

3. 单位化：\(q_k = \eta_k / \|\eta_k\|\)

关键区别：经典方法是用原始向量 \(a_k\) 计算所有投影；而改进方法是用不断更新的向量 \(\eta_k\) 依次与每个已正交向量正交化。这种"边正交化边更新"的方式，大大减少了误差的累积和传播。

2. 实验结果（表8-4）

用改进的格拉姆-施密特正交化方法对同样的希尔伯特矩阵进行QR分解，得到的正交性误差如下：

矩阵阶数 \(n\)	4	5	6	7	8	9
正交性误差 \(e\)	\(2.8047 \times 10^{-13}\)	\(1.2004 \times 10^{-11}\)	\(5.7325 \times 10^{-10}\)	\(1.6075 \times 10^{-8}\)	\(2.5853 \times 10^{-7}\)	\(3.0337 \times 10^{-6}\)

3. 结果对比与分析

• 精度提升显著：对于 \(n=9\) 的希尔伯特矩阵，改进方法的误差为 \(3 \times 10^{-6}\)，而经典方法的误差接近1，精度提高了约6个数量级。
• 误差增长平缓：改进方法的误差随阶数 \(n\) 呈线性增长（约每增加1阶，误差增大一个数量级），而经典方法的误差呈指数增长。
• 局限性：即使是改进方法，对于 \(n=9\) 的希尔伯特矩阵，误差仍然达到了 \(10^{-6}\) 量级，这是由希尔伯特矩阵极高的条件数和双精度浮点数的有限精度决定的，无法完全避免。

四、三种正交化方法的数值稳定性对比

目前常用的QR分解方法有三种，它们的数值稳定性和计算复杂度对比如下：

正交化方法	数值稳定性	每步计算量	适用场景
经典格拉姆-施密特	不稳定	\(2n^3/3\)	理论教学，不推荐实际使用
改进的格拉姆-施密特	条件稳定	\(2n^3/3\)	对精度要求不高的场景，或作为其他算法的子步骤
Householder变换	数值稳定	\(4n^3/3\)	工业界标准方法，推荐用于QR分解
Givens旋转	数值稳定	\(3n^3/2\)	特别适合处理稀疏矩阵和带状矩阵

重要说明：

• 数值稳定：指计算过程中舍入误差不会被无限制放大，计算结果的误差与机器精度和问题的条件数成正比。
• 条件稳定：指对于良态问题是稳定的，但对于极端病态问题仍然可能出现较大误差。
• 虽然Householder变换的计算量比改进的格拉姆-施密特大，但由于其极佳的数值稳定性，它是目前实现QR方法的标准选择。

五、例题核心结论总结

1. 经典格拉姆-施密特正交化存在严重的数值稳定性问题，对于病态矩阵，正交性会迅速丧失，导致QR分解完全失效。
2. 改进的格拉姆-施密特正交化通过改变计算顺序，显著提高了数值稳定性，但对于极端病态矩阵仍然有局限性。
3. 希尔伯特矩阵是检验数值方法稳定性的极佳工具，它能放大数值方法的任何微小缺陷。
4. 在实际工程应用中，应优先使用Householder变换或Givens旋转进行QR分解，它们具有严格的数值稳定性，能保证计算结果的可靠性。

---

8.3.2 豪斯霍尔德变换 详细讲解与证明

我将继续以资深数值分析研究员的身份，系统讲解数值稳定的正交变换——豪斯霍尔德(Householder)反射变换，所有关键证明依据均以粗体标注，深入剖析其几何意义、核心性质与QR分解应用。

一、豪斯霍尔德变换的定义与基本性质

1. 定义（定义8.2）

设向量 \(w \in \mathbb{R}^n\) 满足 \(w^T w = 1\)（单位向量），称矩阵

\[H(w) = I - 2 w w^T \]

为初等反射矩阵，也称为豪斯霍尔德变换矩阵。

其元素形式为：

\[H(w) = \begin{pmatrix} 1 - 2w_1^2 & -2w_1w_2 & \dots & -2w_1w_n \\ -2w_2w_1 & 1 - 2w_2^2 & \dots & -2w_2w_n \\ \vdots & \vdots & & \vdots \\ -2w_nw_1 & -2w_nw_2 & \dots & 1 - 2w_n^2 \end{pmatrix}\]

一般形式：对于任意非零向量 \(u \in \mathbb{R}^n\)，令 \(w = u/\|u\|_2\)，则对应的初等反射矩阵为：

\[H = I - 2 \frac{u u^T}{\|u\|_2^2} \]

2. 基本性质（定理8.13）

设初等反射矩阵 \(H = I - 2 w w^T\)，其中 \(w^T w = 1\)，则：

1. 对称性：\(H^T = H\)
2. 正交性：\(H^{-1} = H\)（即 \(H\) 是自逆矩阵）
3. 保对称性：若 \(A\) 是对称矩阵，则 \(A_1 = H^{-1} A H = H A H\) 也是对称矩阵

证明：

• 对称性：\(H^T = (I - 2 w w^T)^T = I^T - 2 (w w^T)^T = I - 2 w w^T = H\)
• 正交性：依据：矩阵乘法分配律与 \(w^T w = 1\)
  \[\begin{align*} H^T H &= H^2 = (I - 2 w w^T)(I - 2 w w^T) \\ &= I - 4 w w^T + 4 w (w^T w) w^T \\ &= I - 4 w w^T + 4 w \cdot 1 \cdot w^T \\ &= I \end{align*} \]

• 保对称性：\((H A H)^T = H^T A^T H^T = H A H\)，故 \(A_1\) 对称。

二、豪斯霍尔德变换的几何意义

豪斯霍尔德变换是一种镜面反射变换，其几何意义非常直观：

考虑以单位向量 \(w\) 为法向量且过原点的超平面 \(S: w^T x = 0\)。对任意向量 \(v \in \mathbb{R}^n\)，可唯一分解为：

\[v = x + y \]

其中 \(x \in S\)（在超平面内），\(y \in S^\perp\)（垂直于超平面），即 \(y = \mu w\)（\(\mu\) 为常数）。

变换效果：

• 对超平面内的向量 \(x\)：\(H x = (I - 2 w w^T) x = x - 2 w (w^T x) = x\)（保持不变）
• 对垂直于超平面的向量 \(y\)：\(H y = (I - 2 w w^T)(\mu w) = \mu w - 2 \mu w (w^T w) = -\mu w = -y\)（方向反转）

因此，对任意向量 \(v\)：

\[H v = H(x + y) = x - y = v' \]

其中 \(v'\) 是 \(v\) 关于超平面 \(S\) 的镜面反射。

三、向量约化定理

豪斯霍尔德变换在数值计算中的核心价值在于：它能将任意非零向量变换为与第一个单位向量 \(e_1\) 同方向的向量，这是实现矩阵QR分解的基础。

1. 一般向量约化定理（定理8.14）

设 \(x, y\) 为两个不相等的 \(n\) 维向量，但 \(\|x\|_2 = \|y\|_2\)，则存在一个初等反射矩阵 \(H\)，使得 \(H x = y\)。

证明：
构造单位向量：

\[w = \frac{x - y}{\|x - y\|_2} \]

对应的初等反射矩阵为：

\[H = I - 2 w w^T = I - 2 \frac{(x - y)(x^T - y^T)}{\|x - y\|_2^2} \]

计算 \(H x\)：

\[\begin{align*} H x &= x - 2 \frac{(x - y)(x^T - y^T) x}{\|x - y\|_2^2} \\ &= x - 2 \frac{(x - y)(x^T x - y^T x)}{\|x - y\|_2^2} \end{align*} \]

依据：向量2-范数的定义，\(\|x\|_2^2 = x^T x\)，\(\|y\|_2^2 = y^T y\)，且已知 \(\|x\|_2 = \|y\|_2\)，故：

\[\begin{align*} \|x - y\|_2^2 &= (x - y)^T (x - y) \\ &= x^T x - 2 x^T y + y^T y \\ &= 2(x^T x - x^T y) \\ &= 2(x^T x - y^T x) \end{align*} \]

代入上式得：

\[H x = x - 2 \frac{(x - y)(x^T x - y^T x)}{2(x^T x - y^T x)} = x - (x - y) = y \]

证毕。

2. 标准约化定理（定理8.15）

设 \(x = (x_1, x_2, \dots, x_n)^T \neq 0\)，则存在初等反射矩阵 \(H\)，使得

\[H x = -\sigma e_1 \]

其中 \(e_1 = (1, 0, \dots, 0)^T\) 是第一个单位向量，且：

\[\begin{cases} \sigma = \operatorname{sgn}(x_1) \|x\|_2 \\ u = x + \sigma e_1 \\ \beta = \frac{1}{2} \|u\|_2^2 = \sigma(\sigma + x_1) \\ H = I - \beta^{-1} u u^T \end{cases} \tag{8.11} \]

证明：
令 \(y = -\sigma e_1\)，则 \(\|y\|_2 = |\sigma| = \|x\|_2\)。由定理8.14，存在初等反射矩阵 \(H\) 使得 \(H x = y\)。

构造向量 \(u = x - y = x + \sigma e_1\)，则：

\[\begin{align*} \|u\|_2^2 &= (x_1 + \sigma)^2 + x_2^2 + \dots + x_n^2 \\ &= x_1^2 + 2 \sigma x_1 + \sigma^2 + x_2^2 + \dots + x_n^2 \\ &= (x_1^2 + \dots + x_n^2) + \sigma^2 + 2 \sigma x_1 \\ &= \sigma^2 + \sigma^2 + 2 \sigma x_1 \\ &= 2 \sigma(\sigma + x_1) \end{align*} \]

因此 \(\beta = \frac{1}{2} \|u\|_2^2 = \sigma(\sigma + x_1)\)，对应的初等反射矩阵为：

\[H = I - 2 \frac{u u^T}{\|u\|_2^2} = I - \beta^{-1} u u^T \]

证毕。

3. 数值稳定性考虑

• 符号选择：取 \(\sigma\) 与 \(x_1\) 同号，即 \(\sigma = \operatorname{sgn}(x_1) \|x\|_2\)。这样可以保证 \(x_1 + \sigma\) 的绝对值不会太小，避免有效数字损失。
• 溢出处理：在计算 \(\|x\|_2\) 时，可能出现上溢或下溢。为避免此问题，可先将向量 \(x\) 规范化为 \(x' = x / \|x\|_\infty\)，再对 \(x'\) 进行变换。

四、例题8.7 详细计算

设向量 \(x = (3, 5, 1, 1)^T\)，构造豪斯霍尔德矩阵 \(H\)，使得 \(H x = -\sigma e_1\)。

步骤1：计算 \(\sigma\)

\[\|x\|_2 = \sqrt{3^2 + 5^2 + 1^2 + 1^2} = \sqrt{9 + 25 + 1 + 1} = \sqrt{36} = 6 \]

因 \(x_1 = 3 > 0\)，故 \(\sigma = \operatorname{sgn}(3) \times 6 = 6\)。

步骤2：计算向量 \(u\)

\[u = x + \sigma e_1 = (3, 5, 1, 1)^T + (6, 0, 0, 0)^T = (9, 5, 1, 1)^T \]

步骤3：计算 \(\beta\)

\[\|u\|_2^2 = 9^2 + 5^2 + 1^2 + 1^2 = 81 + 25 + 1 + 1 = 108 \]

\[\beta = \frac{1}{2} \|u\|_2^2 = 54 \]

步骤4：构造豪斯霍尔德矩阵 \(H\)

\[\begin{align*} H &= I - \beta^{-1} u u^T \\ &= \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} - \frac{1}{54} \begin{pmatrix} 9 \\ 5 \\ 1 \\ 1 \end{pmatrix} \begin{pmatrix} 9 & 5 & 1 & 1 \end{pmatrix} \\ &= \frac{1}{54} \begin{pmatrix} 54 - 81 & -45 & -9 & -9 \\ -45 & 54 - 25 & -5 & -5 \\ -9 & -5 & 54 - 1 & -1 \\ -9 & -5 & -1 & 54 - 1 \end{pmatrix} \\ &= \frac{1}{54} \begin{pmatrix} -27 & -45 & -9 & -9 \\ -45 & 29 & -5 & -5 \\ -9 & -5 & 53 & -1 \\ -9 & -5 & -1 & 53 \end{pmatrix} \end{align*} \]

步骤5：验证结果

\[H x = \frac{1}{54} \begin{pmatrix} -27 & -45 & -9 & -9 \\ -45 & 29 & -5 & -5 \\ -9 & -5 & 53 & -1 \\ -9 & -5 & -1 & 53 \end{pmatrix} \begin{pmatrix} 3 \\ 5 \\ 1 \\ 1 \end{pmatrix} = \frac{1}{54} \begin{pmatrix} -324 \\ 0 \\ 0 \\ 0 \end{pmatrix} = \begin{pmatrix} -6 \\ 0 \\ 0 \\ 0 \end{pmatrix} = -\sigma e_1 \]

验证正确。

五、用豪斯霍尔德变换实现QR分解

豪斯霍尔德变换是目前实现矩阵QR分解的标准方法，它具有极佳的数值稳定性。

1. 基本思想

通过一系列初等反射变换，逐步将矩阵 \(A\) 化为上三角矩阵 \(R\)。每一步变换都保持矩阵的正交性，最终所有变换矩阵的乘积就是正交矩阵 \(Q\)。

2. 详细算法步骤

设 \(A \in \mathbb{R}^{n \times n}\) 非奇异，记 \(A^{(0)} = A\)。

第1步：约化第一列

• 取 \(A^{(0)}\) 的第一列 \(a_1^{(0)}\)
• 构造 \(n\) 阶豪斯霍尔德矩阵 \(H_1\)，使得 \(H_1 a_1^{(0)} = -\sigma_1 e_1\)
• 计算 \(A^{(1)} = H_1 A^{(0)}\)，则 \(A^{(1)}\) 具有形式：
  \[A^{(1)} = \begin{pmatrix} -\sigma_1 & b^{(1)} \\ 0 & \bar{A}^{(1)} \end{pmatrix} \]
  其中 \(\bar{A}^{(1)} \in \mathbb{R}^{(n-1) \times (n-1)}\) 是右下角的子矩阵。

第j步：约化第j列（j=2,3,…,n-1）

• 取 \(\bar{A}^{(j-1)}\) 的第一列 \(a_1^{(j-1)}\)
• 构造 \(n-j+1\) 阶豪斯霍尔德矩阵 \(\bar{H}_j\)，使得 \(\bar{H}_j a_1^{(j-1)} = -\sigma_j e_1\)
• 构造 \(n\) 阶分块对角矩阵：
  \[H_j = \begin{pmatrix} I_{j-1} & 0 \\ 0 & \bar{H}_j \end{pmatrix} \]

• 计算 \(A^{(j)} = H_j A^{(j-1)}\)，则 \(A^{(j)}\) 的前 \(j\) 列已化为上三角形式。

最终结果
经过 \(n-1\) 步变换后，得到：

\[A^{(n-1)} = H_{n-1} \dots H_2 H_1 A = R \]

其中 \(R\) 是上三角矩阵。令 \(P = H_{n-1} \dots H_2 H_1\)，则 \(P\) 是正交矩阵（正交矩阵的乘积仍为正交矩阵），且 \(P A = R\)，因此：

\[A = P^T R = Q R \]

其中 \(Q = P^T\) 是正交矩阵，\(R\) 是上三角矩阵。

3. 算法优势

• 数值稳定：豪斯霍尔德变换是正交变换，不会放大舍入误差
• 计算效率高：总计算量约为 \(4n^3/3\)，比改进的格拉姆-施密特方法略高，但稳定性好得多
• 存储需求低：无需额外存储正交矩阵 \(Q\)，可将变换信息存储在原矩阵的下三角部分

六、知识点系统总结表

知识点类别	核心内容	理论依据	关键备注与应用
豪斯霍尔德矩阵定义	\(H = I - 2ww^T\)（\(w\) 为单位向量）	矩阵乘法定义	是对称正交矩阵，且自逆（\(H^2=I\)）
几何意义	关于以 \(w\) 为法向量的超平面的镜面反射	向量分解定理	保持向量长度不变
向量约化定理	可将任意非零向量变换为与 \(e_1\) 同方向的向量	定理8.14	是QR分解的基础
数值稳定性技巧	\(\sigma\) 与 \(x_1\) 同号，避免有效数字损失	浮点数运算特性	保证算法在极端情况下仍能稳定运行
QR分解算法	通过 \(n-1\) 次豪斯霍尔德变换将矩阵化为上三角矩阵	正交变换的乘积仍为正交矩阵	工业界标准QR分解方法
算法复杂度	总计算量约为 \(4n^3/3\)	矩阵乘法复杂度分析	比格拉姆-施密特方法稳定，比Givens旋转高效

---

8.3.3 吉文斯变换与QR分解唯一性 详细讲解与证明

我将继续以资深数值分析研究员的身份，系统讲解第二种数值稳定的正交变换——吉文斯(Givens)旋转变换，所有关键证明依据均以粗体标注，并深入对比三种QR分解方法的优劣，最后证明QR分解的唯一性定理。

一、吉文斯变换的定义与基本性质

1. 二维平面旋转变换

在二维平面 \(\mathbb{R}^2\) 中，旋转变换将向量绕原点逆时针旋转 \(\theta\) 角，其变换矩阵为：

\[P(\theta) = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix} \]

性质：\(P(\theta)\) 是正交矩阵，满足 \(P(\theta)^T P(\theta) = I\)，且旋转变换保持向量的长度不变。

2. n维吉文斯变换

将二维旋转变换推广到n维空间，得到吉文斯变换矩阵（也称为平面旋转矩阵）\(P(i,j,\theta)\)，它表示在 \(\mathbb{R}^n\) 中绕 \((x_i, x_j)\) 坐标平面旋转 \(\theta\) 角的变换。

其矩阵形式为：

\[P(i,j,\theta) = \begin{pmatrix} 1 & & & & & & \\ & \ddots & & & & & \\ & & 1 & & & & \\ & & & \cos\theta & \dots & \sin\theta & \\ & & & 0 & 1 & 0 & \\ & & & \vdots & \ddots & \vdots & \\ & & & 0 & & 1 & 0 \\ & & & -\sin\theta & \dots & \cos\theta & \\ & & & & & & 1 \\ & & & & & & & \ddots \\ & & & & & & & & 1 \end{pmatrix} \begin{matrix} \\ \\ \\ \leftarrow i \text{行} \\ \\ \\ \\ \leftarrow j \text{行} \\ \\ \\ \\ \end{matrix} \]

即：

• 对角元：\(p_{kk} = 1\)（\(k \neq i,j\)），\(p_{ii} = p_{jj} = \cos\theta\)
• 非对角元：\(p_{ij} = \sin\theta\)，\(p_{ji} = -\sin\theta\)，其余非对角元均为0

3. 基本性质

吉文斯变换矩阵具有以下重要性质：

1. 稀疏性：与单位矩阵 \(I\) 仅在 \((i,i), (i,j), (j,i), (j,j)\) 四个位置元素不同，其余元素完全相同。
2. 正交性：\(P(i,j,\theta)\) 是正交矩阵，即 \(P(i,j,\theta)^{-1} = P(i,j,\theta)^T = P(i,j,-\theta)\)。
3. 左乘效果：\(P(i,j)A\) 仅改变矩阵 \(A\) 的第 \(i\) 行和第 \(j\) 行，其余行不变：
   \[\begin{pmatrix} a_{il}' \\ a_{jl}' \end{pmatrix} = \begin{pmatrix} c & s \\ -s & c \end{pmatrix} \begin{pmatrix} a_{il} \\ a_{jl} \end{pmatrix}, \quad l=1,2,\dots,n \]
   其中 \(c = \cos\theta\)，\(s = \sin\theta\)。
4. 右乘效果：\(A P(i,j)\) 仅改变矩阵 \(A\) 的第 \(i\) 列和第 \(j\) 列，其余列不变：
   \[(a_{li}', a_{lj}') = (a_{li}, a_{lj}) \begin{pmatrix} c & s \\ -s & c \end{pmatrix}, \quad l=1,2,\dots,m \]

二、向量约化定理（定理8.16）

吉文斯变换的核心作用是：可以将向量中任意指定的一个元素化为零，而不改变向量的长度。

定理8.16 设向量 \(x = (x_1, \dots, x_i, \dots, x_j, \dots, x_n)^T\)，其中 \(x_i, x_j\) 不全为零，则存在平面旋转矩阵 \(P(i,j,\theta)\)，使得：

\[P x = (x_1, \dots, x_i', \dots, 0, \dots, x_n)^T \]

其中 \(x_i' = \sqrt{x_i^2 + x_j^2}\)，\(\theta = \arctan(x_j / x_i)\)。

证明：
取 \(c = \cos\theta = x_i / x_i'\)，\(s = \sin\theta = x_j / x_i'\)，其中 \(x_i' = \sqrt{x_i^2 + x_j^2}\)。

根据吉文斯变换的左乘性质，变换后向量的第 \(i\) 和第 \(j\) 个分量为：

\[\begin{cases} x_i' = c x_i + s x_j \\ x_j' = -s x_i + c x_j \end{cases} \]

代入 \(c\) 和 \(s\) 的表达式：

\[x_i' = \frac{x_i^2}{\sqrt{x_i^2 + x_j^2}} + \frac{x_j^2}{\sqrt{x_i^2 + x_j^2}} = \sqrt{x_i^2 + x_j^2} \]

\[x_j' = -\frac{x_i x_j}{\sqrt{x_i^2 + x_j^2}} + \frac{x_i x_j}{\sqrt{x_i^2 + x_j^2}} = 0 \]

其余分量保持不变，证毕。

数值稳定的计算方法

为了避免计算过程中出现溢出或有效数字损失，实际应用中采用以下算法计算 \(c\) 和 \(s\)：

1. 如果 \(x_j = 0\)，则取 \(c=1\)，\(s=0\)（无需旋转）。
2. 否则：
   • 若 \(|x_j| \geq |x_i|\)，令 \(t = x_i / x_j\)，则 \(s = 1 / \sqrt{1 + t^2}\)，\(c = s t\)。
   • 若 \(|x_j| < |x_i|\)，令 \(t = x_j / x_i\)，则 \(c = 1 / \sqrt{1 + t^2}\)，\(s = c t\)。

这种方法通过选择较大的分量作为分母，避免了除法运算中出现大数除以小数的情况，保证了数值稳定性。

三、用吉文斯变换实现QR分解

吉文斯变换也可以用来实现矩阵的QR分解，其基本思想是通过一系列平面旋转变换，逐步将矩阵的下三角元素化为零，最终得到上三角矩阵。

1. 算法步骤

设 \(A \in \mathbb{R}^{n \times n}\) 非奇异，记 \(A^{(1)} = A\)。

第1步：约化第一列
对于 \(j=2,3,\dots,n\)：

• 构造吉文斯变换矩阵 \(P(1,j)\)，将 \(A^{(1)}\) 的 \((j,1)\) 位置元素化为零
• 计算 \(A^{(2)} = P(1,n) \dots P(1,2) A^{(1)}\)，则 \(A^{(2)}\) 的第一列除第一个元素外均为零。

第k步：约化第k列（k=2,3,…,n-1）
假设前 \(k-1\) 步已完成，矩阵 \(A^{(k)}\) 的前 \(k-1\) 列已化为上三角形式。对于 \(j=k+1,k+2,\dots,n\)：

• 构造吉文斯变换矩阵 \(P(k,j)\)，将 \(A^{(k)}\) 的 \((j,k)\) 位置元素化为零
• 计算 \(A^{(k+1)} = P(k,n) \dots P(k,k+1) A^{(k)}\)，则 \(A^{(k+1)}\) 的前 \(k\) 列已化为上三角形式。

最终结果
经过 \(n-1\) 步变换后，得到：

\[A^{(n)} = P_{n-1} \dots P_2 P_1 A = R \]

其中 \(R\) 是上三角矩阵，\(P_i\) 是第 \(i\) 步所有吉文斯变换矩阵的乘积。令 \(P = P_{n-1} \dots P_2 P_1\)（正交矩阵），则 \(P A = R\)，因此：

\[A = P^T R = Q R \]

其中 \(Q = P^T\) 是正交矩阵。

2. 与豪斯霍尔德变换的对比

对比维度	豪斯霍尔德变换	吉文斯变换
每次变换消零数	一次消去一列的多个元素	一次只能消去一个元素
稠密矩阵计算量	\(4n^3/3\)	\(3n^3/2\)（比豪斯霍尔德多50%）
稀疏矩阵适用性	差（会引入大量非零元素）	好（仅改变两行两列，保持稀疏性）
并行计算	较难	容易（不同列的变换可并行执行）
数值稳定性	极佳	极佳

结论：

• 对于稠密矩阵，优先使用豪斯霍尔德变换，计算量更小。
• 对于稀疏矩阵或带状矩阵，优先使用吉文斯变换，能有效保持矩阵的稀疏结构。

四、QR分解的唯一性定理（定理8.17）

1. 定理内容

设 \(A \in \mathbb{R}^{n \times n}\) 为非奇异矩阵，则存在正交矩阵 \(Q\) 与上三角矩阵 \(R\)，使得

\[A = Q R \]

且当 \(R\) 的对角元素为正时，分解是唯一的。

2. 详细证明

• 存在性：前面介绍的格拉姆-施密特正交化、豪斯霍尔德变换、吉文斯变换三种方法都已经证明了QR分解的存在性。

• 唯一性：假设 \(A\) 有两种QR分解：

  \[A = Q_1 R_1 = Q_2 R_2 \]

  其中 \(Q_1, Q_2\) 是正交矩阵，\(R_1, R_2\) 是对角元为正的上三角矩阵。

  则有：

  \[Q_2^T Q_1 = R_2 R_1^{-1} \]

  依据：正交矩阵的性质，左边 \(Q_2^T Q_1\) 是正交矩阵（正交矩阵的乘积和转置仍为正交矩阵）。
  依据：上三角矩阵的性质，右边 \(R_2 R_1^{-1}\) 是上三角矩阵（上三角矩阵的逆和乘积仍为上三角矩阵）。

  因此，\(Q_2^T Q_1\) 是一个既是正交矩阵又是上三角矩阵的矩阵。既是正交矩阵又是上三角矩阵的矩阵必为对角矩阵，且对角元的绝对值为1。

  又因为 \(R_1\) 和 \(R_2\) 的对角元均为正，所以 \(R_2 R_1^{-1}\) 的对角元也为正。因此，\(Q_2^T Q_1\) 是对角元为1的对角矩阵，即单位矩阵：

  \[Q_2^T Q_1 = I \implies Q_1 = Q_2 \]

  代入原式得 \(R_1 = R_2\)，故分解唯一。证毕。

3. 非唯一情况

如果不要求 \(R\) 的对角元为正，则QR分解不唯一。例如，若 \(A = Q R\) 是一个QR分解，令 \(D = \operatorname{diag}(\pm 1, \pm 1, \dots, \pm 1)\)（对角元为±1的对角矩阵），则 \(A = (Q D)(D R)\) 也是一个QR分解，其中 \(Q D\) 是正交矩阵，\(D R\) 是上三角矩阵。

五、三种QR分解方法的系统对比

方法	数值稳定性	计算量（n阶稠密矩阵）	适用场景	主要缺点
经典格拉姆-施密特	不稳定	\(2n^3/3\)	理论教学	正交性损失严重，不推荐实际使用
改进的格拉姆-施密特	条件稳定	\(2n^3/3\)	对精度要求不高的场景	极端病态矩阵仍会失效
豪斯霍尔德变换	数值稳定	\(4n^3/3\)	中小型稠密矩阵	会破坏矩阵的稀疏性
吉文斯变换	数值稳定	\(3n^3/2\)	稀疏矩阵、带状矩阵、并行计算	稠密矩阵计算量较大

教材验证：对于例8.6中的9阶希尔伯特矩阵，用正交变换（豪斯霍尔德或吉文斯）进行QR分解，得到的正交矩阵 \(Q\) 的正交性误差 \(e\) 不超过 \(10^{-15}\)，远优于格拉姆-施密特方法的 \(10^{-1}\) 量级，充分验证了正交变换的数值稳定性优势。

六、知识点系统总结表

知识点类别	核心内容	理论依据	关键备注与应用
吉文斯变换定义	n维平面旋转变换矩阵，仅改变第i,j行/列元素	二维旋转变换推广	正交矩阵，保持向量长度不变
向量约化定理	可将向量中任意指定元素化为零	旋转变换的坐标变换性质	是吉文斯QR分解的基础
吉文斯QR分解	通过一系列平面旋转变换将矩阵化为上三角	正交变换的乘积仍为正交矩阵	特别适合稀疏矩阵和并行计算
QR分解唯一性	非奇异矩阵的QR分解在R对角元为正时唯一	正交矩阵和上三角矩阵的性质	保证了QR方法的理论基础
方法对比	豪斯霍尔德适合稠密矩阵，吉文斯适合稀疏矩阵	计算复杂度和稀疏性分析	实际应用中根据矩阵类型选择合适方法

---

8.3.4 舒尔分解与上黑森伯格矩阵约化 详细讲解与证明

我将继续以资深数值分析研究员的身份，系统讲解实用QR方法的核心预处理技术——上黑森伯格矩阵约化，以及矩阵理论中的重要结果——舒尔分解，所有关键证明依据均以粗体标注，深入剖析其理论意义与工程价值。

一、舒尔分解与实舒尔分解

1. 复舒尔分解

舒尔分解是矩阵理论中最基本的分解之一，它回答了"矩阵在相似变换下能约化到什么程度"的问题。

复舒尔分解定理：对任意复矩阵 \(A \in \mathbb{C}^{n \times n}\)，存在酉矩阵 \(U\)（满足 \(U^H U = I\)），使得

\[U^H A U = R \]

其中 \(R\) 是上三角矩阵，其对角元恰好是 \(A\) 的全部特征值。

这意味着：任何复矩阵都酉相似于一个上三角矩阵。这是一个非常强的结果，因为上三角矩阵的特征值就是其对角元，一目了然。

2. 实舒尔分解定理（定理8.18）

对于实矩阵 \(A \in \mathbb{R}^{n \times n}\)，由于其特征值可能是复数，我们不能用实正交相似变换将其约化为上三角矩阵（否则复特征值会出现在实矩阵的对角线上，矛盾）。但我们可以将其约化为分块上三角矩阵，即实舒尔型。

定理8.18（实舒尔分解） 设 \(A \in \mathbb{R}^{n \times n}\)，则存在正交矩阵 \(Q\)，使得

\[Q^T A Q = \begin{pmatrix} R_{11} & R_{12} & \dots & R_{1m} \\ & R_{22} & \dots & R_{2m} \\ & & \ddots & \vdots \\ & & & R_{mm} \end{pmatrix} \tag{8.12} \]

其中每个对角块 \(R_{ii}\) 是一阶或二阶方阵：

• 一阶块 \(R_{ii}\) 是 \(A\) 的实特征值
• 二阶块 \(R_{ii}\) 的两个特征值是 \(A\) 的一对共轭复特征值

证明思路：对矩阵阶数 \(n\) 用数学归纳法。当 \(n=1\) 时显然成立。假设对 \(n-1\) 阶矩阵成立，对于 \(n\) 阶矩阵，取其一个特征值 \(\lambda\) 和对应的特征向量 \(x\)，构造正交矩阵 \(Q_1\) 使得 \(Q_1^T x = \|x\| e_1\)，则 \(Q_1^T A Q_1\) 具有分块上三角形式，再对右下角的 \(n-1\) 阶子矩阵应用归纳假设即可。

3. 舒尔分解的意义

• 理论意义：实舒尔分解是实矩阵在实正交相似变换下的标准形，它完全揭示了实矩阵的特征值结构。
• 实用意义：虽然实舒尔分解本身无法直接数值实现，但它为QR方法提供了理论基础——QR方法产生的矩阵序列本质上收敛于实舒尔型。
• 计算效率：基本QR方法每步迭代的计算量为 \(O(n^3)\)，如果先将矩阵约化为上黑森伯格矩阵，再进行QR迭代，每步计算量可降为 \(O(n^2)\)，这是QR方法实用化的关键。

二、上黑森伯格矩阵的定义与约化动机

1. 上黑森伯格矩阵定义

如果矩阵 \(H \in \mathbb{R}^{n \times n}\) 的次对角线以下的元素均为零，即

\[h_{ij} = 0, \quad \forall i > j+1 \]

则称 \(H\) 为上黑森伯格(Hessenberg)矩阵。

其形式为：

\[H = \begin{pmatrix} h_{11} & h_{12} & h_{13} & \dots & h_{1n} \\ h_{21} & h_{22} & h_{23} & \dots & h_{2n} \\ & h_{32} & h_{33} & \dots & h_{3n} \\ & & \ddots & \ddots & \vdots \\ & & & h_{n,n-1} & h_{nn} \end{pmatrix} \]

2. 约化动机

基本QR方法直接对原矩阵进行迭代，每步需要对满矩阵做QR分解，计算量为 \(O(n^3)\)。而如果我们先通过正交相似变换将原矩阵约化为上黑森伯格矩阵，那么：

• 上黑森伯格矩阵的QR分解可以高效实现，每步仅需 \(O(n^2)\) 次运算
• 正交相似变换保持特征值不变，因此约化后的矩阵与原矩阵有完全相同的特征值
• 约化过程本身的计算量为 \(O(n^3)\)，但只需做一次，后续每步迭代的计算量大幅降低

对于 \(n\) 较大的矩阵，这种预处理带来的效率提升是极其显著的。

三、豪斯霍尔德变换约化上黑森伯格矩阵（定理8.19）

1. 基本思想

通过一系列分块豪斯霍尔德变换，逐步将矩阵次对角线以下的元素消为零。每一步变换只作用于矩阵的右下角子矩阵，保持已经约化好的上黑森伯格结构不变。

2. 第一步约化（k=1）

设原矩阵为：

\[A = \begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & & \vdots \\ a_{n1} & a_{n2} & \dots & a_{nn} \end{pmatrix} = \begin{pmatrix} a_{11} & A_{12}^{(1)} \\ c_1 & A_{22}^{(1)} \end{pmatrix} \]

其中 \(c_1 = (a_{21}, a_{31}, \dots, a_{n1})^T \in \mathbb{R}^{n-1}\) 是第一列除第一个元素外的部分。

依据：豪斯霍尔德向量约化定理（定理8.15），存在 \(n-1\) 阶初等反射矩阵 \(H_1\)，使得

\[H_1 c_1 = -\sigma_1 e_1 \]

其中 \(\sigma_1 = \operatorname{sgn}(a_{21}) \|c_1\|_2\)，\(e_1 = (1, 0, \dots, 0)^T \in \mathbb{R}^{n-1}\)。

构造 \(n\) 阶分块正交矩阵：

\[U_1 = \begin{pmatrix} 1 & 0 \\ 0 & H_1 \end{pmatrix} \]

对 \(A\) 做正交相似变换：

\[\begin{align*} A_2 &= U_1 A U_1 \\ &= \begin{pmatrix} 1 & 0 \\ 0 & H_1 \end{pmatrix} \begin{pmatrix} a_{11} & A_{12}^{(1)} \\ c_1 & A_{22}^{(1)} \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & H_1 \end{pmatrix} \\ &= \begin{pmatrix} a_{11} & A_{12}^{(1)} H_1 \\ H_1 c_1 & H_1 A_{22}^{(1)} H_1 \end{pmatrix} \\ &= \begin{pmatrix} a_{11} & * & * & \dots & * \\ -\sigma_1 & * & * & \dots & * \\ 0 & * & * & \dots & * \\ \vdots & \vdots & \vdots & & \vdots \\ 0 & * & * & \dots & * \end{pmatrix} \end{align*} \]

可以看到，\(A_2\) 的第一列次对角线以下的元素已全部化为零，第一列已满足上黑森伯格结构。

3. 第k步约化（k=2,3,…,n-2）

假设经过 \(k-1\) 步变换后，矩阵 \(A_k\) 已具有如下形式：

\[A_k = \begin{pmatrix} A_{11}^{(k)} & A_{12}^{(k)} \\ 0 & c_k & A_{22}^{(k)} \end{pmatrix} \]

其中 \(A_{11}^{(k)}\) 是 \(k\) 阶上黑森伯格矩阵，\(c_k = (a_{k+1,k}^{(k)}, \dots, a_{nk}^{(k)})^T \in \mathbb{R}^{n-k}\) 是第 \(k\) 列次对角线以下的部分。

同样，构造 \(n-k\) 阶初等反射矩阵 \(H_k\) 使得 \(H_k c_k = -\sigma_k e_1\)，再构造 \(n\) 阶分块正交矩阵：

\[U_k = \begin{pmatrix} I_k & 0 \\ 0 & H_k \end{pmatrix} \]

对 \(A_k\) 做正交相似变换：

\[A_{k+1} = U_k A_k U_k = \begin{pmatrix} A_{11}^{(k)} & A_{12}^{(k)} H_k \\ 0 & H_k c_k & H_k A_{22}^{(k)} H_k \end{pmatrix} \]

变换后，第 \(k\) 列次对角线以下的元素全部化为零，而前 \(k\) 列的上黑森伯格结构保持不变。

4. 最终结果

经过 \(n-2\) 步变换后，得到：

\[H = U_{n-2} \dots U_2 U_1 A U_1 U_2 \dots U_{n-2} = U_0^T A U_0 \]

其中 \(U_0 = U_1 U_2 \dots U_{n-2}\) 是正交矩阵，\(H\) 是上黑森伯格矩阵。

定理8.19 设 \(A \in \mathbb{R}^{n \times n}\)，则存在初等反射矩阵 \(U_1, U_2, \dots, U_{n-2}\)，使得 \(U_0^T A U_0 = H\) 为上黑森伯格矩阵。

5. 计算量分析

• 仅约化矩阵 \(A\) 为上黑森伯格矩阵：约需要 \(\frac{5}{3} n^3\) 次乘法运算
• 如果还需要计算正交变换矩阵 \(U_0\)：额外增加 \(\frac{2}{3} n^3\) 次乘法运算，总共约 \(\frac{7}{3} n^3\) 次

优化技巧：利用豪斯霍尔德矩阵的特殊结构 \(H = I - \beta^{-1} u u^T\)，可以将矩阵乘法 \(H B H\) 转化为向量运算，大幅减少计算量：

\[H B H = B - u w^T - w u^T \]

其中 \(w = \beta^{-1} B u - \frac{\beta^{-2}}{2} (u^T B u) u\)。

四、对称矩阵的特殊情况：约化为对称三对角矩阵（定理8.20）

1. 定理内容

定理8.20 设 \(A \in \mathbb{R}^{n \times n}\) 为对称矩阵，则存在初等反射矩阵 \(U_1, U_2, \dots, U_{n-2}\)，使得

\[U_0^T A U_0 = C = \begin{pmatrix} b_1 & c_1 & & & \\ c_1 & b_2 & c_2 & & \\ & c_2 & b_3 & \ddots & \\ & & \ddots & \ddots & c_{n-1} \\ & & & c_{n-1} & b_n \end{pmatrix} \]

即对称矩阵正交相似于一个对称三对角矩阵。

2. 证明

依据：正交相似变换保持对称性，若 \(A\) 对称，则 \(A_{k+1} = U_k A_k U_k\) 也对称。
依据：上黑森伯格矩阵的定义，约化后的矩阵是上黑森伯格矩阵。

一个既是对称矩阵又是上黑森伯格矩阵的矩阵，必然是对称三对角矩阵。证毕。

3. 计算量优势

• 仅约化矩阵为对称三对角矩阵：约需要 \(\frac{2}{3} n^3\) 次乘法运算
• 如果还需要计算正交变换矩阵 \(U_0\)：额外增加 \(\frac{2}{3} n^3\) 次乘法运算，总共约 \(\frac{4}{3} n^3\) 次

这比一般矩阵的约化计算量减少了约60%，是对称矩阵特征值求解效率极高的根本原因。

五、例题8.9 详细解析

用豪斯霍尔德方法将矩阵

\[A = \begin{pmatrix} -4 & -3 & -7 \\ 2 & 3 & 2 \\ 4 & 2 & 7 \end{pmatrix} \]

约化为上黑森伯格矩阵。

步骤1：提取第一列待约化部分

\[c_1 = (a_{21}, a_{31})^T = (2, 4)^T \]

步骤2：构造豪斯霍尔德矩阵 \(H_1\)

• 计算 \(\sigma_1\)：\(\|c_1\|_2 = \sqrt{2^2 + 4^2} = \sqrt{20} = 2\sqrt{5} \approx 4.472136\)，因 \(c_1\) 的第一个元素 \(2 > 0\)，故 \(\sigma_1 = 2\sqrt{5} \approx 4.472136\)
• 计算向量 \(u_1\)：\(u_1 = c_1 + \sigma_1 e_1 = (2 + 2\sqrt{5}, 4)^T \approx (6.472136, 4)^T\)
• 计算 \(\beta_1\)：\(\beta_1 = \sigma_1(\sigma_1 + c_{11}) = 2\sqrt{5}(2\sqrt{5} + 2) = 20 + 4\sqrt{5} \approx 28.944272\)
• 构造 \(H_1\)：
  \[H_1 = I - \beta_1^{-1} u_1 u_1^T \approx \begin{pmatrix} -0.447214 & -0.894427 \\ -0.894427 & 0.447214 \end{pmatrix} \]

步骤3：构造分块正交矩阵 \(U_1\)

\[U_1 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & -0.447214 & -0.894427 \\ 0 & -0.894427 & 0.447214 \end{pmatrix} \]

步骤4：计算正交相似变换

\[\begin{align*} U_1 A &= \begin{pmatrix} -4 & -3 & -7 \\ -4.472136 & -3.130495 & -7.155418 \\ 2.683282 & 1.788854 & 5.366563 \end{pmatrix} \\ H = U_1 A U_1 &= \begin{pmatrix} -4 & 7.602631 & -0.447214 \\ -4.472136 & 7.8 & 0.2 \\ 0 & 0.4 & 2.2 \end{pmatrix} \end{align*} \]

结果验证

\(H\) 是上黑森伯格矩阵，且与原矩阵 \(A\) 有相同的特征值：\(\lambda_1=2, \lambda_2=3, \lambda_3=1\)，验证正确。

六、知识点系统总结表

知识点类别	核心内容	理论依据	关键备注与应用
复舒尔分解	任意复矩阵酉相似于上三角矩阵	数学归纳法、酉矩阵性质	揭示了复矩阵的本质结构
实舒尔分解	任意实矩阵正交相似于分块上三角矩阵，一阶块为实特征值，二阶块为共轭复特征值	复舒尔分解、实矩阵特征值性质	是实矩阵QR方法的理论基础
上黑森伯格矩阵	次对角线以下元素全为零的矩阵	矩阵结构定义	是实用QR方法的预处理目标
一般矩阵约化	通过n-2次豪斯霍尔德变换将矩阵约化为上黑森伯格矩阵	豪斯霍尔德约化定理、正交相似变换性质	计算量约5n³/3，将后续QR迭代计算量从O(n³)降为O(n²)
对称矩阵约化	对称矩阵正交相似于对称三对角矩阵	正交变换保持对称性、上黑森伯格结构	计算量约2n³/3，是对称矩阵特征值求解的关键步骤
计算量对比	一般矩阵约化：5n³/3；对称矩阵约化：2n³/3	矩阵乘法复杂度分析	对称矩阵特征值求解效率远高于一般矩阵

---

8.4 QR方法 详细讲解与证明

我将继续以资深数值分析研究员的身份，系统讲解求解矩阵全部特征值的工业标准算法——实用QR方法，所有关键证明依据均以粗体标注，深入剖析其迭代原理、加速技术与工程实现细节。

一、上黑森伯格矩阵的基本QR迭代

1. 预处理的必要性

基本QR方法直接对满矩阵进行迭代，每步需要对 \(n\) 阶满矩阵做QR分解，计算量为 \(O(n^3)\)。对于 \(n=1000\) 的矩阵，这需要数十亿次运算，完全不实用。

依据：上黑森伯格矩阵的结构特性，如果先将原矩阵约化为上黑森伯格矩阵，再进行QR迭代，则每步迭代的计算量可降为 \(O(n^2)\)，这是QR方法实用化的关键突破。

2. 上黑森伯格矩阵的定义与不可约性

设 \(B\) 为上黑森伯格矩阵，即

\[B = \begin{pmatrix} b_{11} & b_{12} & \dots & b_{1n} \\ b_{21} & b_{22} & \dots & b_{2n} \\ & \ddots & \ddots & \vdots \\ & & b_{n,n-1} & b_{nn} \end{pmatrix} \]

如果 \(B\) 的所有次对角元都不为零，即 \(b_{i+1,i} \neq 0\)（\(i=1,2,\dots,n-1\)），则称 \(B\) 为不可约上黑森伯格矩阵。

重要性质：不可约上黑森伯格矩阵的QR迭代具有结构保持性，即迭代产生的矩阵序列 \(\{H_k\}\) 仍然是上黑森伯格矩阵。

证明：设 \(H_k\) 是上黑森伯格矩阵，对其进行QR分解 \(H_k = Q_k R_k\)。由于上黑森伯格矩阵的次对角线以下元素全为零，用吉文斯旋转进行QR分解时，只需消去 \(n-1\) 个次对角元，因此正交矩阵 \(Q_k\) 是一个下黑森伯格矩阵（上对角线以上元素全为零）。

下黑森伯格矩阵与上黑森伯格矩阵的乘积 \(R_k Q_k\) 仍然是上黑森伯格矩阵，因此 \(H_{k+1} = R_k Q_k\) 保持上黑森伯格结构。证毕。

3. 基本QR迭代公式

对于上黑森伯格矩阵 \(H\)，基本QR迭代公式为：

\[\begin{cases} H_1 = H \\ H_k = Q_k R_k \quad (\text{QR分解}) \\ H_{k+1} = R_k Q_k \end{cases} \quad k=1,2,\dots \tag{8.13} \]

关键性质：迭代产生的所有矩阵都与原矩阵相似，因此具有完全相同的特征值。
证明：\(H_{k+1} = R_k Q_k = Q_k^T Q_k R_k Q_k = Q_k^T H_k Q_k\)，依据：相似矩阵的定义，\(H_{k+1}\) 与 \(H_k\) 相似。由数学归纳法，所有 \(H_k\) 都与原矩阵 \(H\) 相似。

4. 可约性与分治策略

在迭代过程中，如果某个次对角元 \(h_{p+1,p}\) 变得足够小，我们就可以将其视为零，此时矩阵 \(H_{k+1}\) 可约化为分块上三角矩阵：

\[H_{k+1} = \begin{pmatrix} H_{11} & H_{12} \\ 0 & H_{22} \end{pmatrix} \]

其中 \(H_{11}\) 是 \(p\) 阶上黑森伯格矩阵，\(H_{22}\) 是 \(n-p\) 阶上黑森伯格矩阵。

此时，原问题被分解为两个独立的子问题：分别求解 \(H_{11}\) 和 \(H_{22}\) 的特征值。这就是QR方法的分治策略，它可以大幅减少计算量。

特别地，当 \(p=n-1\) 时，\(H_{22}\) 是1阶矩阵，其唯一元素就是原矩阵的一个特征值；当 \(p=n-2\) 时，\(H_{22}\) 是2阶矩阵，其特征值可以通过解二次方程直接求得。

5. 收敛判据

实际计算中，当次对角元满足以下条件时，就可以将其视为零：

\[|h_{p+1,p}| \leq \varepsilon (|h_{pp}| + |h_{p+1,p+1}|) \]

其中 \(\varepsilon = 10^{-t}\)，\(t\) 是计算中有效数字的位数（通常取 \(t=6 \sim 15\)）。

这个判据考虑了对角元的大小，是相对误差判据，比绝对误差判据更合理。

二、带原点位移的QR方法

1. 基本QR方法的收敛性问题

基本QR方法的收敛速度由比值 \(r_n = |\lambda_n / \lambda_{n-1}|\) 决定。当特征值比较接近时，\(r_n\) 接近1，收敛会变得极其缓慢。例如，若 \(r_n=0.99\)，则需要约460次迭代才能将误差缩小到原来的1%。

为了解决这个问题，我们引入原点位移技术，它可以将收敛速度从线性收敛提高到二次收敛，是实用QR方法的核心。

2. 带原点位移的QR迭代公式

设 \(s_k\) 是第 \(k\) 步的位移量，带原点位移的QR迭代公式为：

\[\begin{cases} A_1 = A \\ A_k - s_k I = Q_k R_k \quad (\text{QR分解}) \\ A_{k+1} = R_k Q_k + s_k I \end{cases} \quad k=1,2,\dots \]

关键性质：迭代仍然保持相似性，特征值不变。
证明：

\[A_{k+1} = R_k Q_k + s_k I = Q_k^T Q_k R_k Q_k + s_k I = Q_k^T (A_k - s_k I) Q_k + s_k I = Q_k^T A_k Q_k \]

依据：相似矩阵的定义，\(A_{k+1}\) 与 \(A_k\) 相似，因此具有相同的特征值。证毕。

3. 加速原理

设 \(A\) 的特征值为 \(\lambda_1, \lambda_2, \dots, \lambda_n\)，则 \(A - s I\) 的特征值为 \(\lambda_1 - s, \lambda_2 - s, \dots, \lambda_n - s\)。

基本QR方法对 \(A - s I\) 的收敛速度由比值

\[r' = \left| \frac{\lambda_n - s}{\lambda_{n-1} - s} \right| \]

决定。如果我们能选择 \(s\) 非常接近 \(\lambda_n\)，则 \(r'\) 会变得非常小，收敛速度将极快。

4. 位移量的选择策略

最简单且最有效的位移策略是瑞利位移，即取当前迭代矩阵的右下角元素作为位移量：

\[s_k = h_{nn}^{(k)} \]

为什么这是一个好的选择？
当迭代接近收敛时，矩阵的右下角元素 \(h_{nn}^{(k)}\) 会趋近于最小的特征值 \(\lambda_n\)。此时，位移量 \(s_k \approx \lambda_n\)，收敛比值

\[r' = \left| \frac{\lambda_n - s_k}{\lambda_{n-1} - s_k} \right| \approx \left| \frac{\varepsilon}{\lambda_{n-1} - \lambda_n} \right| \]

其中 \(\varepsilon\) 是 \(h_{nn}^{(k)}\) 与 \(\lambda_n\) 的误差。这个比值通常非常小，使得收敛速度达到二次收敛，即每次迭代误差的位数翻倍。

对于对称矩阵，瑞利位移可以达到三次收敛，即每次迭代误差的位数翻三倍，这是极其惊人的收敛速度。

5. 实用实现技巧

在实际计算中，我们不需要显式生成和保存正交矩阵 \(Q_k\)。根据QR分解的定义，\(Q_k^T (A_k - s_k I) = R_k\)，因此：

\[A_{k+1} = Q_k^T (A_k - s_k I) Q_k + s_k I = P_{n-1} \dots P_2 P_1 (A_k - s_k I) P_1^T P_2^T \dots P_{n-1}^T + s_k I \]

其中 \(P_i\) 是吉文斯旋转矩阵。

这意味着我们可以直接通过一系列吉文斯旋转的左右乘来得到下一个迭代矩阵 \(A_{k+1}\)，无需存储 \(Q_k\)，大幅节省了存储和计算量。

三、实用QR方法的完整流程

工业级QR方法的完整流程分为以下三个阶段：

阶段1：预处理——约化为上黑森伯格矩阵

1. 输入原矩阵 \(A \in \mathbb{R}^{n \times n}\)
2. 用豪斯霍尔德变换将 \(A\) 正交相似约化为上黑森伯格矩阵 \(H\)
   • 一般矩阵：计算量约 \(\frac{5}{3} n^3\)
   • 对称矩阵：约化为对称三对角矩阵，计算量约 \(\frac{2}{3} n^3\)

阶段2：带原点位移的QR迭代

对 \(H\) 进行带瑞利位移的QR迭代：

1. 初始化：\(H_1 = H\)
2. 对于 \(k=1,2,\dots\)：
   1. 检查所有次对角元，若满足收敛判据则进行分块，将问题分解为更小的子问题
   2. 若所有子问题都已解决（都是1阶或2阶矩阵），则终止迭代
   3. 选择位移量 \(s_k = h_{nn}^{(k)}\)（瑞利位移）
   4. 用吉文斯旋转对 \(H_k - s_k I\) 进行QR分解
   5. 计算下一个迭代矩阵 \(H_{k+1} = R_k Q_k + s_k I\)

阶段3：提取特征值

1. 对于1阶子块，其唯一元素就是一个实特征值
2. 对于2阶子块 \(\begin{pmatrix} a & b \\ c & d \end{pmatrix}\)，解二次特征方程
   \[\lambda^2 - (a+d)\lambda + (ad - bc) = 0 \]
   得到两个特征值，它们可能是一对共轭复特征值

四、QR方法的收敛性与优势

1. 收敛性

• 基本QR方法：线性收敛，收敛速度由特征值的比值决定
• 带瑞利位移的QR方法：对于一般矩阵，通常是二次收敛；对于对称矩阵，是三次收敛

这意味着对于大多数矩阵，带位移的QR方法只需10~20次迭代即可达到机器精度，是目前已知收敛最快的特征值算法。

2. 核心优势

1. 数值稳定：整个过程仅使用正交变换，不会放大舍入误差
2. 收敛速度快：二次/三次收敛，迭代次数极少
3. 计算效率高：预处理后每步迭代仅需 \(O(n^2)\) 次运算
4. 功能全面：可以求解实矩阵的所有实特征值和复特征值，无需引入复数运算
5. 适用范围广：适用于所有非亏损矩阵，是目前求解中小型稠密矩阵全部特征值的工业标准

五、知识点系统总结表

知识点类别	核心内容	理论依据	关键备注与应用
预处理必要性	先将矩阵约化为上黑森伯格矩阵，再进行QR迭代	上黑森伯格矩阵的结构保持性	将每步迭代计算量从O(n³)降为O(n²)
基本QR迭代	\(H_k=Q_kR_k, H_{k+1}=R_kQ_k\)	相似矩阵的定义、正交变换性质	保持上黑森伯格结构，特征值不变
分治策略	当次对角元足够小时，将矩阵分块为两个独立子问题	分块上三角矩阵的特征值性质	大幅减少计算量，是收敛判断的基础
收敛判据	$	h_	\leq \varepsilon(
带原点位移的QR迭代	\(H_k-s_kI=Q_kR_k, H_{k+1}=R_kQ_k+s_kI\)	相似矩阵的定义、特征值平移性质	保持相似性，大幅加速收敛
瑞利位移策略	取当前矩阵右下角元素作为位移量 \(s_k=h_{nn}^{(k)}\)	迭代收敛时右下角元素趋近于最小特征值	实现简单，收敛速度达到二次/三次
实用流程	预处理→带位移QR迭代→提取特征值	正交变换的数值稳定性	工业级标准算法，广泛应用于科学计算

---

8.4.2 上黑森伯格矩阵的单步QR方法 详细讲解与证明

我将继续以资深数值分析研究员的身份，深入剖析实用QR方法的核心实现细节——上黑森伯格矩阵的单步QR迭代，所有关键证明依据均以粗体标注，并通过完整的算法步骤和例题解析，展示工业级QR方法的实际运行过程。

一、上黑森伯格矩阵单步QR方法的计算过程

对于不可约上黑森伯格矩阵 \(H\)，带原点位移的单步QR迭代公式为：

\[\begin{cases} H_k - s_k I = Q_k R_k \quad (\text{QR分解}) \\ H_{k+1} = R_k Q_k + s_k I \end{cases} \]

由于上黑森伯格矩阵的特殊结构，我们无需显式生成和存储正交矩阵 \(Q_k\)，而是通过一系列吉文斯旋转变换的左右乘来直接计算 \(H_{k+1}\)，这将每步迭代的计算量从 \(O(n^3)\) 降为 \(O(n^2)\)。

1. 左变换：将 \(H_1 - s_1 I\) 化为上三角矩阵 \(R_1\)

左变换的目标是用一系列吉文斯旋转将位移后的矩阵 \(H_1 - s_1 I\) 化为上三角矩阵。对于上黑森伯格矩阵，只有次对角线上有非零元素，因此只需 \(n-1\) 次吉文斯旋转即可完成QR分解。

第1次左变换（k=1）：

• 取矩阵的前两行前两列元素 \(\begin{pmatrix} h_{11}-s_1 & h_{12} \\ h_{21} & h_{22}-s_1 \end{pmatrix}\)
• 构造吉文斯旋转矩阵 \(P_{12} = P(1,2,\theta)\)，使得
  \[P_{12} \begin{pmatrix} h_{11}-s_1 \\ h_{21} \end{pmatrix} = \begin{pmatrix} r_{11} \\ 0 \end{pmatrix} \]

• 左乘 \(P_{12}\) 到整个矩阵上，消去 \((2,1)\) 位置的元素。

第k次左变换（k=2,3,…,n-1）：

• 经过前 \(k-1\) 次变换后，矩阵的前 \(k-1\) 列已化为上三角形式
• 取第 \(k\) 行和第 \(k+1\) 行的第 \(k\) 列元素 \(h_{kk}^{(k)}\) 和 \(h_{k+1,k}^{(k)}\)
• 构造吉文斯旋转矩阵 \(P_{k,k+1} = P(k,k+1,\theta)\)，使得
  \[P_{k,k+1} \begin{pmatrix} h_{kk}^{(k)} \\ h_{k+1,k}^{(k)} \end{pmatrix} = \begin{pmatrix} r_{kk} \\ 0 \end{pmatrix} \]

• 左乘 \(P_{k,k+1}\) 到整个矩阵上，消去 \((k+1,k)\) 位置的元素。

关键优化：每次吉文斯旋转仅改变矩阵的第 \(k\) 行和第 \(k+1\) 行，且前 \(k-1\) 个元素已经是零，无需计算。因此，每次左变换仅需计算 \(2(n-k+1)\) 次乘法运算。

经过 \(n-1\) 次左变换后，得到上三角矩阵 \(R_1\)：

\[P_{n-1,n} \dots P_{23} P_{12} (H_1 - s_1 I) = R_1 \]

2. 右变换：计算 \(H_2 = R_1 P_{12}^T P_{23}^T \dots P_{n-1,n}^T + s_1 I\)

右变换的目标是将左变换中使用的所有吉文斯旋转依次右乘到上三角矩阵 \(R_1\) 上，再加上位移矩阵 \(s_1 I\)，得到下一个迭代矩阵 \(H_2\)。

第k次右变换（k=1,2,…,n-1）：

• 依次将 \(P_{k,k+1}^T\) 右乘到当前矩阵上
• 每次右乘仅改变矩阵的第 \(k\) 列和第 \(k+1\) 列
• 由于 \(R_1\) 是上三角矩阵，右乘后矩阵的 \((j,k)\) 和 \((j,k+1)\) 位置元素（\(j \geq k+2\)）仍然为零

关键性质：右变换后得到的矩阵仍然是上黑森伯格矩阵。这是因为每次右乘吉文斯旋转矩阵只会在次对角线上引入非零元素，不会在次对角线以下产生新的非零元素。

最后，将所有对角元加上位移量 \(s_1\)，得到：

\[H_2 = R_1 P_{12}^T P_{23}^T \dots P_{n-1,n}^T + s_1 I \]

3. 计算量分析

• 左变换：共需 \(\sum_{k=1}^{n-1} 2(n-k+1) = n^2 + n - 2\) 次乘法运算
• 右变换：共需 \(\sum_{k=1}^{n-1} 2(k+1) = n^2 + 3n - 4\) 次乘法运算
• 总计：约 \(2n^2\) 次乘法运算

这就是为什么上黑森伯格矩阵的单步QR迭代计算量为 \(O(n^2)\)，比满矩阵的QR迭代快了一个数量级。

二、瑞利位移的理论基础（定理8.21）

1. 定理内容

定理8.21 设：

1. \(H \in \mathbb{R}^{n \times n}\) 为不可约上黑森伯格矩阵
2. \(\mu\) 是 \(H\) 的一个特征值

则对 \(H - \mu I\) 进行QR分解 \(H - \mu I = QR\)，并构造 \(H_2 = RQ + \mu I\)，则有：

\[h_{n,n-1}^{(2)} = 0, \quad h_{nn}^{(2)} = \mu \]

即一步迭代就精确收敛到特征值 \(\mu\)。

2. 详细证明

证明：
记 \(R\) 为上三角矩阵：

\[R = \begin{pmatrix} r_{11} & r_{12} & \dots & r_{1n} \\ & r_{22} & \dots & r_{2n} \\ & & \ddots & \vdots \\ & & & r_{nn} \end{pmatrix} \]

依据：不可约上黑森伯格矩阵的性质，若 \(H\) 不可约，则 \(H - \mu I\) 也不可约，且其所有次对角元均不为零。

依据：上黑森伯格矩阵QR分解的性质，用平面旋转变换将不可约上黑森伯格矩阵化为上三角矩阵时，得到的上三角矩阵 \(R\) 的所有对角元均不为零，即 \(r_{ii} \neq 0\)（\(i=1,2,\dots,n-1\)）。

由于 \(\mu\) 是 \(H\) 的特征值，故 \(H - \mu I\) 是奇异矩阵，因此 \(R = Q^T (H - \mu I)\) 也是奇异矩阵。而上三角矩阵奇异当且仅当至少有一个对角元为零。结合上面的结论，只能是 \(r_{nn} = 0\)。

现在计算 \(H_2 = RQ + \mu I\) 的最后一行。由于 \(R\) 的最后一行是 \((0,0,\dots,0)\)，因此 \(RQ\) 的最后一行也是 \((0,0,\dots,0)\)。因此：

\[H_2 \text{ 的最后一行} = (0,0,\dots,0) + (0,0,\dots,\mu) = (0,0,\dots,\mu) \]

即 \(h_{n,n-1}^{(2)} = 0\)，\(h_{nn}^{(2)} = \mu\)。证毕。

3. 重要推论

这个定理为瑞利位移策略提供了坚实的理论基础：如果位移量 \(s_k\) 正好是特征值 \(\lambda_n\)，那么一步迭代就能精确收敛。在实际计算中，当迭代接近收敛时，\(h_{nn}^{(k)}\) 会非常接近 \(\lambda_n\)，因此选择 \(s_k = h_{nn}^{(k)}\) 作为位移量，会使收敛速度达到二次收敛甚至三次收敛。

三、上黑森伯格矩阵QR方法的完整算法（算法8.1）

算法8.1（上黑森伯格矩阵的QR方法）

输入：上黑森伯格矩阵 \(H \in \mathbb{R}^{n \times n}\)，收敛精度 \(\varepsilon\)
输出：矩阵 \(H\) 的全部特征值

1. 初始化：令当前矩阵为 \(H\)，当前阶数为 \(m = n\)
2. 迭代循环：当 \(m > 2\) 时：
   1. 收敛检查：从下往上检查所有次对角元 \(h_{i+1,i}\)（\(i=1,2,\dots,m-1\)）
      • 若 \(|h_{i+1,i}| \leq \varepsilon (|h_{ii}| + |h_{i+1,i+1}|)\)，则将 \(h_{i+1,i}\) 置为零，将矩阵分块为两个子矩阵，对阶数较小的子矩阵递归调用本算法，令 \(m = i\)，继续处理左上角的子矩阵
   2. 选择位移量：取瑞利位移 \(s = h_{mm}\)
   3. 单步QR迭代：
      1. 对 \(H - sI\) 进行左变换，用 \(n-1\) 次吉文斯旋转化为上三角矩阵 \(R\)
      2. 对 \(R\) 进行右变换，依次右乘所有吉文斯旋转的转置
      3. 所有对角元加上位移量 \(s\)，得到新的迭代矩阵 \(H\)
3. 处理小矩阵：
   • 若 \(m = 1\)，则唯一元素 \(h_{11}\) 就是一个实特征值
   • 若 \(m = 2\)，则解二次特征方程 \(\lambda^2 - (h_{11}+h_{22})\lambda + (h_{11}h_{22} - h_{12}h_{21}) = 0\)，得到两个特征值
4. 返回所有特征值

算法说明

• 分治策略：每当一个次对角元收敛到零，就将矩阵分解为两个独立的子问题，分别求解，这大幅减少了计算量
• 瑞利位移：每次迭代都取当前矩阵的右下角元素作为位移量，保证了快速收敛
• 数值稳定性：整个过程仅使用吉文斯旋转变换，具有极佳的数值稳定性

四、例题8.10 详细解析

用QR方法计算对称三对角矩阵

\[A = \begin{pmatrix} 2 & 1 & 0 \\ 1 & 3 & 1 \\ 0 & 1 & 4 \end{pmatrix} \]

的全部特征值。

精确特征值：\(\lambda_1 = 3 + \sqrt{3} \approx 4.7320508\)，\(\lambda_2 = 3\)，\(\lambda_3 = 3 - \sqrt{3} \approx 1.2679492\)

第1次迭代（k=1）

1. 选择位移量：\(s_1 = h_{33} = 4\)
2. 构造位移矩阵：
   \[A_1 - s_1 I = \begin{pmatrix} 2-4 & 1 & 0 \\ 1 & 3-4 & 1 \\ 0 & 1 & 4-4 \end{pmatrix} = \begin{pmatrix} -2 & 1 & 0 \\ 1 & -1 & 1 \\ 0 & 1 & 0 \end{pmatrix} \]

3. 左变换（吉文斯旋转）：
   • 第1次旋转 \(P_{12}\)：消去 \((2,1)\) 元素，\(c=-2/\sqrt{5} \approx -0.8944\)，\(s=1/\sqrt{5} \approx 0.4472\)
   • 第2次旋转 \(P_{23}\)：消去 \((3,2)\) 元素，\(c=-1/\sqrt{1.8} \approx -0.7454\)，\(s=\sqrt{0.8}/\sqrt{1.8} \approx 0.6667\)
   • 得到上三角矩阵：
     \[R = \begin{pmatrix} 2.2361 & -1.3416 & 0.4472 \\ 0 & 1.0954 & -0.3651 \\ 0 & 0 & 0.8165 \end{pmatrix} \]

4. 右变换：依次右乘 \(P_{12}^T\) 和 \(P_{23}^T\)，再加上位移量 \(s_1=4\)，得到：
   \[A_2 = \begin{pmatrix} 1.4000 & -0.4899 & 0 \\ -0.4899 & 3.2667 & -0.7454 \\ 0 & -0.7454 & 4.3333 \end{pmatrix} \]

第2次迭代（k=2）

1. 选择位移量：\(s_2 = h_{33} = 4.3333\)
2. 重复上述过程，得到：
   \[A_3 = \begin{pmatrix} 1.2915 & -0.2017 & 0 \\ -0.2017 & 3.0225 & 0.2725 \\ 0 & 0.2725 & 4.6860 \end{pmatrix} \]

第3次迭代（k=3）

1. 选择位移量：\(s_3 = h_{33} = 4.6860\)
2. 迭代后得到：
   \[A_4 = \begin{pmatrix} 1.2737 & -0.0993 & 0 \\ -0.0993 & 2.9942 & -0.0732 \\ 0 & -0.0732 & 4.7320 \end{pmatrix} \]

第4次迭代（k=4）

1. 选择位移量：\(s_4 = h_{33} = 4.7320\)
2. 迭代后得到：
   \[A_5 = \begin{pmatrix} 1.2694 & -0.0498 & 0 \\ -0.0498 & 2.9986 & 0 \\ 0 & 0 & 4.7321 \end{pmatrix} \]

收敛与降阶

此时，次对角元 \(h_{32} = 1.235 \times 10^{-7}\) 已经充分小，可以视为零。因此，我们得到一个特征值 \(\lambda_3 \approx 4.7321\)，并将矩阵降阶为2阶子矩阵：

\[\tilde{A}_5 = \begin{pmatrix} 1.2694 & -0.0498 \\ -0.0498 & 2.9986 \end{pmatrix} \]

求解2阶子矩阵

对 \(\tilde{A}_5\) 再进行一次迭代，选择位移量 \(s_5 = 2.9986\)，得到：

\[\tilde{A}_6 = \begin{pmatrix} 1.2680 & 4 \times 10^{-5} \\ 4 \times 10^{-5} & 3.0000 \end{pmatrix} \]

因此，得到另外两个特征值：\(\lambda_2 \approx 3.0000\)，\(\lambda_1 \approx 1.2680\)。

结果对比

• 计算值：\(\lambda_1 \approx 1.2680\)，\(\lambda_2 \approx 3.0000\)，\(\lambda_3 \approx 4.7321\)
• 精确值：\(\lambda_1 = 1.2679492\)，\(\lambda_2 = 3.0\)，\(\lambda_3 = 4.7320508\)
• 最大误差：约 \(5 \times 10^{-5}\)，仅经过5次迭代就达到了很高的精度，充分展示了带瑞利位移的QR方法的快速收敛特性。

五、算法局限性与扩展

算法8.1主要用于对称矩阵的特征值计算。对于非对称矩阵，可能会出现复特征值，而瑞利位移是实数，无法逼近复特征值，此时需要使用双步QR方法（也称为隐式QR方法）。双步QR方法通过同时进行两次带复共轭位移的迭代，避免了复数运算，仍然保持在实数域内计算，是求解一般实矩阵全部特征值的工业标准方法。

六、知识点系统总结表

知识点类别	核心内容	理论依据	关键备注与应用
单步QR迭代计算	通过n-1次吉文斯旋转的左右乘实现上黑森伯格矩阵的QR迭代	吉文斯旋转的性质、上黑森伯格矩阵的结构保持性	每步迭代计算量为O(n²)，比满矩阵快一个数量级
左变换	用吉文斯旋转将H-sI化为上三角矩阵R	平面旋转变换的消元能力	每次旋转仅改变两行，计算量小
右变换	将吉文斯旋转依次右乘到R上，再加上位移sI	正交变换的转置性质	保持上黑森伯格结构不变
定理8.21	若位移量正好是特征值，则一步迭代精确收敛	不可约上黑森伯格矩阵的性质、QR分解的唯一性	为瑞利位移策略提供了理论基础
瑞利位移策略	取当前矩阵右下角元素作为位移量s_k=h_{nn}^	定理8.21、迭代收敛性	实现简单，收敛速度达到二次/三次
完整算法流程	预处理→带位移QR迭代→收敛检查→分治降阶→提取特征值	正交变换的数值稳定性、分治策略	工业级标准算法，广泛应用于科学计算
算法局限性	主要用于对称矩阵，不能直接处理复特征值	瑞利位移是实数	非对称矩阵需使用双步QR方法

---

8.4.3 双步QR方法（隐式QR方法）详细讲解与证明

我将继续以资深数值分析研究员的身份，系统讲解求解一般实矩阵全部特征值的工业标准算法——双步QR方法（也称为隐式QR方法），所有关键证明依据均以粗体标注，深入剖析其如何巧妙避免复数运算，同时保持极高的计算效率。

一、单步QR方法的局限性与双步方法的引入

1. 单步QR方法的致命缺陷

上一节介绍的带瑞利位移的单步QR方法对于对称矩阵非常有效，因为对称矩阵的所有特征值都是实数。但对于非对称实矩阵，其特征值可能是共轭复数对。此时，瑞利位移（取右下角元素，实数）无法逼近复特征值，导致单步QR方法收敛极慢甚至不收敛。

2. 问题的根源

当迭代接近收敛时，矩阵右下角会出现一个2×2的子块：

\[C = \begin{pmatrix} h_{n-1,n-1} & h_{n-1,n} \\ h_{n,n-1} & h_{nn} \end{pmatrix} \]

如果这个子块的特征值是一对共轭复数 \(s_1\) 和 \(s_2 = \overline{s_1}\)，那么单步QR方法的实数位移无法同时逼近这两个复特征值，导致次对角元 \(h_{n,n-1}\) 无法收敛到零。

3. 解决方案：双步QR方法

为了解决这个问题，我们引入双步QR方法，其核心思想是：

• 同时进行两次带复共轭位移 \(s_1\) 和 \(s_2\) 的单步QR迭代
• 通过巧妙的数学变换，将两次复迭代合并为一次实运算，完全避免复数运算
• 整个过程保持在实数域内进行，计算量仍为 \(O(n^2)\)

二、隐式Q定理（定理8.22）

隐式Q定理是双步QR方法的理论基石，它保证了我们可以通过构造第一列相同的正交矩阵，得到与两次单步迭代完全相同的结果。

1. 定理内容

定理8.22（隐式Q定理） 设 \(A \in \mathbb{R}^{n \times n}\)，且：

1. \(Q = (q_1, q_2, \dots, q_n)\) 和 \(V = (v_1, v_2, \dots, v_n)\) 都是正交矩阵
2. \(Q^T A Q = H\) 和 \(V^T A V = G\) 都是不可约上黑森伯格矩阵
3. \(q_1 = v_1\)（即 \(Q\) 与 \(V\) 的第一列相同）

则：

1. \(v_i = \pm q_i\)，且 \(|h_{i,i-1}| = |g_{i,i-1}|\)，\(i=2,3,\dots,n\)
2. \(G = D^{-1} H D\)，其中 \(D = \operatorname{diag}(1, \pm 1, \dots, \pm 1)\)

即 \(H\) 和 \(G\) 本质上相等，仅差一个对角元为±1的对角相似变换。

2. 定理意义

这个定理告诉我们：不可约上黑森伯格矩阵由正交相似变换的第一列唯一确定（差±1因子）。

这是一个极其强大的结论。它意味着，如果我们能构造一个正交矩阵 \(Q'\)，使得：

1. \(Q'\) 的第一列与两次单步迭代得到的正交矩阵 \(Q = Q_1 Q_2\) 的第一列相同
2. \(Q'^T H Q'\) 是不可约上黑森伯格矩阵

那么 \(Q'^T H Q'\) 就与两次单步迭代的结果 \(H_3\) 本质上相等。这正是隐式双步QR方法的核心思想。

三、双步QR方法的基本原理

1. 两次单步迭代的显式形式

设 \(s_1\) 和 \(s_2\) 是右下角2×2子块 \(C\) 的一对共轭复特征值，两次带位移的单步迭代公式为：

\[\begin{cases} H_1 - s_1 I = Q_1 R_1 \\ H_2 = R_1 Q_1 + s_1 I = Q_1^T H_1 Q_1 \\ H_2 - s_2 I = Q_2 R_2 \\ H_3 = R_2 Q_2 + s_2 I = Q_2^T H_2 Q_2 = Q^T H_1 Q \end{cases} \tag{8.15} \]

其中 \(Q = Q_1 Q_2\) 是两次正交变换的乘积。

将两次迭代合并，得到：

\[H_3 = (Q_1 Q_2)^T H_1 (Q_1 Q_2) = Q^T H_1 Q \]

2. 实矩阵M的构造

定义矩阵：

\[M = (H_1 - s_1 I)(H_1 - s_2 I) \]

关键性质：\(M\) 是实矩阵。
证明：因为 \(s_1\) 和 \(s_2\) 是共轭复数，所以它们的和 \(s = s_1 + s_2\) 与积 \(t = s_1 s_2\) 都是实数。展开得：

\[M = H_1^2 - (s_1 + s_2) H_1 + s_1 s_2 I = H_1^2 - s H_1 + t I \]

显然，\(M\) 是实矩阵。证毕。

3. M的QR分解与Q的第一列

由式(8.15)可得：

\[M = (H_1 - s_1 I)(H_1 - s_2 I) = Q_1 R_1 Q_2 R_2 = Q_1 Q_2 R_2 R_1 = Q R \]

其中 \(R = R_2 R_1\) 是上三角矩阵（上三角矩阵的乘积仍为上三角矩阵）。

这说明 \(M = QR\) 是 \(M\) 的一个QR分解。根据QR分解的性质，正交矩阵 \(Q\) 的第一列是 \(M\) 的第一列的单位化：

\[Q e_1 = \frac{M e_1}{\|M e_1\|_2} \]

依据：隐式Q定理，只要我们构造一个正交矩阵 \(Q'\)，其第一列与 \(Q\) 的第一列相同，且 \(Q'^T H_1 Q'\) 是不可约上黑森伯格矩阵，那么 \(Q'^T H_1 Q'\) 就与 \(H_3\) 本质上相等。

四、显式双步QR方法及其缺陷

1. 显式算法步骤

1. 计算实矩阵 \(M = H_1^2 - s H_1 + t I\)
2. 对 \(M\) 进行QR分解 \(M = Q R\)
3. 计算 \(H_3 = Q^T H_1 Q\)

2. 致命缺陷

计算矩阵 \(M\) 需要两次矩阵乘法，计算量为 \(O(n^3)\)，这与直接对满矩阵做QR分解的计算量相同，完全失去了上黑森伯格矩阵预处理带来的效率优势。因此，显式双步QR方法在实际中完全不实用。

五、隐式双步QR方法（工业标准实现）

隐式双步QR方法利用隐式Q定理，巧妙地避免了计算整个矩阵 \(M\)，仅需计算 \(M\) 的第一列，然后通过一系列豪斯霍尔德变换将矩阵重新约化为上黑森伯格矩阵，计算量仍为 \(O(n^2)\)。

1. 关键观察：M的第一列只有前三个元素非零

由于 \(H\) 是上黑森伯格矩阵，其结构为：

\[H = \begin{pmatrix} h_{11} & h_{12} & h_{13} & \dots & h_{1n} \\ h_{21} & h_{22} & h_{23} & \dots & h_{2n} \\ & h_{32} & h_{33} & \dots & h_{3n} \\ & & \ddots & \ddots & \vdots \\ & & & h_{n,n-1} & h_{nn} \end{pmatrix} \]

计算 \(M e_1 = (H^2 - s H + t I) e_1\)：

• \(H e_1 = (h_{11}, h_{21}, 0, \dots, 0)^T\)（只有前两个元素非零）
• \(H^2 e_1 = H (H e_1) = (h_{11}^2 + h_{12} h_{21}, h_{21} h_{11} + h_{22} h_{21}, h_{32} h_{21}, 0, \dots, 0)^T\)（只有前三个元素非零）

因此，\(M e_1\) 只有前三个元素非零：

\[M e_1 = (x, y, z, 0, \dots, 0)^T \]

其中：

\[\begin{cases} x = h_{11}^2 + h_{12} h_{21} - s h_{11} + t \\ y = h_{21}(h_{11} + h_{22} - s) \\ z = h_{21} h_{32} \end{cases} \tag{8.17} \]

这是隐式双步QR方法最关键的发现！我们不需要计算整个 \(n \times n\) 矩阵 \(M\)，只需要计算三个数 \(x, y, z\)。

2. 隐式双步QR方法的完整步骤

设 \(H\) 是不可约上黑森伯格矩阵，带弗朗西斯位移（Francis shift）的隐式双步QR迭代步骤如下：

步骤1：计算位移参数

取矩阵右下角的2×2子块：

\[C = \begin{pmatrix} h_{n-1,n-1} & h_{n-1,n} \\ h_{n,n-1} & h_{nn} \end{pmatrix} \]

计算其迹 \(s = \operatorname{tr}(C) = h_{n-1,n-1} + h_{nn}\) 和行列式 \(t = \det(C) = h_{n-1,n-1} h_{nn} - h_{n-1,n} h_{n,n-1}\)。这两个数就是 \(C\) 的两个特征值的和与积。

步骤2：计算M的第一列

根据式(8.17)计算 \(x, y, z\)，得到 \(M e_1 = (x, y, z, 0, \dots, 0)^T\)。

步骤3：构造第一个豪斯霍尔德矩阵P₀

构造3阶豪斯霍尔德矩阵 \(H_0\)，使得：

\[H_0 \begin{pmatrix} x \\ y \\ z \end{pmatrix} = -\sigma e_1 \]

其中 \(\sigma = \operatorname{sgn}(x) \sqrt{x^2 + y^2 + z^2}\)。

构造 \(n\) 阶分块正交矩阵：

\[P_0 = \begin{pmatrix} H_0 & 0 \\ 0 & I_{n-3} \end{pmatrix} \]

步骤4：第一次相似变换

对 \(H\) 做相似变换：

\[H' = P_0^T H P_0 \]

重要说明：这次变换会在矩阵的 \((3,1)\) 和 \((4,1)\) 位置引入非零元素，破坏上黑森伯格结构，但非零元素仅出现在前四行四列的子矩阵中。

步骤5：逐次恢复上黑森伯格结构

对于 \(k=1,2,\dots,n-2\)：

1. 取矩阵第 \(k\) 列中从第 \(k+2\) 行开始的元素，构造 \(n-k-1\) 阶豪斯霍尔德矩阵 \(\bar{H}_k\)，将这些元素全部消为零
2. 构造 \(n\) 阶分块正交矩阵：
   \[P_k = \begin{pmatrix} I_{k+1} & 0 \\ 0 & \bar{H}_k \end{pmatrix} \]

3. 做相似变换 \(H = P_k^T H P_k\)，恢复第 \(k\) 列的上黑森伯格结构

每一步变换只会在矩阵的下一列引入一个新的非零元素，经过 \(n-2\) 次变换后，矩阵将重新恢复为上黑森伯格矩阵。

步骤6：结果

最终得到的上黑森伯格矩阵 \(H'\) 就是两次单步迭代后的结果 \(H_3\)，与显式两次单步迭代的结果本质上相等。

3. 计算量分析

• 步骤1-3：常数计算量
• 步骤4：约 \(6n\) 次乘法运算
• 步骤5：每步约 \(6(n-k)\) 次乘法运算，总计约 \(3n^2\) 次乘法运算

总计算量：约 \(3n^2\) 次乘法运算，仅比单步QR方法多50%，仍然是 \(O(n^2)\) 量级，效率极高。

六、双步QR方法的收敛性与优势

1. 收敛性

带弗朗西斯位移的双步QR方法具有以下收敛特性：

• 对于一般实矩阵，通常是二次收敛
• 对于对称矩阵，退化为单步QR方法，是三次收敛
• 对于几乎所有矩阵，都能在10~20次迭代内收敛到机器精度

2. 核心优势

1. 完全避免复数运算：通过将两次复共轭位移的迭代合并为一次实运算，整个过程都在实数域内进行
2. 计算效率高：保持了上黑森伯格矩阵QR迭代的 \(O(n^2)\) 计算量
3. 数值稳定：仅使用豪斯霍尔德正交变换，具有极佳的数值稳定性
4. 功能全面：可以求解一般实矩阵的所有实特征值和复特征值
5. 收敛速度快：二次收敛，迭代次数极少

七、实用QR方法的完整工业流程

综合前面所有内容，求解一般实矩阵全部特征值的工业标准流程如下：

1. 预处理：用豪斯霍尔德变换将原矩阵正交相似约化为上黑森伯格矩阵（计算量约 \(\frac{5}{3} n^3\)）
2. 迭代：对约化后的上黑森伯格矩阵应用带弗朗西斯位移的隐式双步QR方法（每步约 \(3n^2\) 次运算）
3. 分治：每当次对角元收敛到零时，将矩阵分块为两个独立的子问题，分别求解
4. 提取特征值：
   • 1阶子块：元素即为实特征值
   • 2阶子块：解二次方程得到两个特征值（可能是共轭复特征值）

八、知识点系统总结表

知识点类别	核心内容	理论依据	关键备注与应用
单步方法局限性	无法处理非对称矩阵的复特征值，收敛极慢	瑞利位移是实数，无法逼近复特征值	必须引入双步方法
隐式Q定理	不可约上黑森伯格矩阵由正交相似变换的第一列唯一确定（差±1因子）	正交变换性质、上黑森伯格矩阵结构	是隐式双步QR方法的理论基石
双步方法思想	合并两次带复共轭位移的单步迭代为一次实运算	共轭复数的和与积都是实数	完全避免复数运算
显式方法缺陷	直接计算M=H²-sH+tI需要O(n³)运算量	矩阵乘法复杂度	不实用，必须使用隐式方法
隐式方法关键	M的第一列只有前三个元素非零，仅需计算三个数	上黑森伯格矩阵的乘法性质	将计算量从O(n³)降为O(n²)
隐式算法步骤	计算位移→计算M的第一列→豪斯霍尔德变换→恢复上黑森伯格结构	豪斯霍尔德约化定理、隐式Q定理	工业标准实现
收敛性	一般矩阵二次收敛，对称矩阵三次收敛	弗朗西斯位移的最优性	迭代次数极少，效率极高
工业流程	预处理→双步QR迭代→分治→提取特征值	正交变换的数值稳定性	求解一般实矩阵全部特征值的标准方法

---

格什戈林圆盘定理特征值界估计 详细解答

定理回顾

格什戈林圆盘定理（行圆盘）：设 \(A=(a_{ij}) \in \mathbb{C}^{n \times n}\)，则 \(A\) 的所有特征值都位于复平面上 \(n\) 个圆盘

\[D_i = \left\{ z \in \mathbb{C} \mid |z - a_{ii}| \leq R_i \right\}, \quad i=1,2,\dots,n \]

的并集之中，其中 \(R_i = \sum_{\substack{j=1 \\ j \neq i}}^n |a_{ij}|\) 是第 \(i\) 行非对角元的绝对值之和（圆盘半径）。

连通分支性质：若 \(k\) 个圆盘构成一个连通区域，且与其余 \(n-k\) 个圆盘不相交，则该区域中恰好包含 \(A\) 的 \(k\) 个特征值（重特征值按重数计算）。

---

(1) 3阶矩阵特征值界估计

\[A = \begin{pmatrix} -1 & 0 & 0 \\ -1 & 0 & 1 \\ -1 & -1 & 2 \end{pmatrix} \]

步骤1：计算行圆盘

圆盘序号	中心 \(a_{ii}\)	半径 \(R_i\)	圆盘表达式	实轴对应区间
1	\(-1\)	\(|0|+|0|=0\)	\(|z+1| \leq 0\)	\(\{-1\}\)（孤立点）
2	\(0\)	\(|-1|+|1|=2\)	\(|z| \leq 2\)	\([-2, 2]\)
3	\(2\)	\(|-1|+|-1|=2\)	\(|z-2| \leq 2\)	\([0, 4]\)

步骤2：连通性分析与特征值分布

1. 孤立圆盘：\(D_1 = \{-1\}\) 与其他圆盘不相交，因此矩阵必有一个特征值等于 \(\boldsymbol{-1}\)。
2. 连通区域：\(D_2\) 和 \(D_3\) 相交于 \([0,2]\)，合并为连通区域 \([-2, 4]\)，因此另外两个特征值位于区间 \([-2, 4]\) 内。

步骤3：列圆盘定理缩小范围（可选）

计算列圆盘：

• 第1列：中心 \(-1\)，半径 \(1+1=2\)，区间 \([-3, 1]\)
• 第2列：中心 \(0\)，半径 \(0+1=1\)，区间 \([-1, 1]\)
• 第3列：中心 \(2\)，半径 \(0+1=1\)，区间 \([1, 3]\)

行圆盘与列圆盘的交集为 \(\{-1\} \cup [-2, 3]\)，因此可将另外两个特征值的范围进一步缩小到 \(\boldsymbol{[-2, 3]}\)。

---

(2) n阶三对角矩阵特征值界估计

\[A = \begin{pmatrix} 4 & -1 & & & \\ -1 & 4 & -1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 4 & -1 \\ & & & -1 & 4 \end{pmatrix}\]

步骤1：计算行圆盘

• 边界行（第1行和第n行）：中心 \(4\)，半径 \(|-1|=1\)，圆盘 \(|z-4| \leq 1\)，对应区间 \([3, 5]\)
• 中间行（第2~n-1行）：中心 \(4\)，半径 \(|-1|+|-1|=2\)，圆盘 \(|z-4| \leq 2\)，对应区间 \([2, 6]\)

步骤2：连通性分析与特征值分布

所有圆盘相互连通：

• \(D_1=[3,5]\) 与 \(D_2=[2,6]\) 相交于 \([3,5]\)
• 中间圆盘完全重叠
• \(D_{n-1}=[2,6]\) 与 \(D_n=[3,5]\) 相交于 \([3,5]\)

因此所有圆盘的并集为一个连通区域 \(\boldsymbol{[2, 6]}\)，根据连通分支性质，矩阵的所有n个特征值都位于区间 \([2, 6]\) 内。

补充说明

1. 对称性验证：该矩阵是对称矩阵，所有特征值均为实数，无需考虑复平面区域。
2. 精确特征值参考：该矩阵是二阶中心差分矩阵，其精确特征值为
   \[\lambda_k = 4 + 2\cos\left( \frac{k\pi}{n+1} \right), \quad k=1,2,\dots,n \]
   由于 \(\cos\theta \in (-1,1)\)，因此特征值实际范围为 \((2,6)\)，与格什戈林定理给出的闭区间界一致。

---

矩阵特征值、特征向量计算与对角化判断 详细解答

核心理论依据

• 特征值定义：满足 \(\det(A-\lambda I)=0\) 的 \(\lambda\) 为矩阵 \(A\) 的特征值
• 特征向量定义：满足 \((A-\lambda I)x=0\) 的非零向量 \(x\) 为对应于 \(\lambda\) 的特征向量
• 可对角化判定：\(n\) 阶矩阵相似于对角矩阵的充要条件是每个特征值的代数重数等于其几何重数（即有 \(n\) 个线性无关的特征向量）

---

(1) 矩阵 \(A_1 = \begin{pmatrix} 2 & -3 & 6 \\ 0 & 3 & -4 \\ 0 & 2 & -3 \end{pmatrix}\)

步骤1：计算特征值

特征多项式按第一列展开（第一列仅第一个元素非零）：

\[\begin{align*} \det(A_1-\lambda I) &= \begin{vmatrix} 2-\lambda & -3 & 6 \\ 0 & 3-\lambda & -4 \\ 0 & 2 & -3-\lambda \end{vmatrix} \\ &= (2-\lambda) \cdot \begin{vmatrix} 3-\lambda & -4 \\ 2 & -3-\lambda \end{vmatrix} \\ &= (2-\lambda) \left[ (3-\lambda)(-3-\lambda) + 8 \right] \\ &= (2-\lambda)(\lambda^2 - 1) \\ &= -(λ-2)(λ-1)(λ+1) \end{align*} \]

令特征多项式为零，得特征值：\(\boldsymbol{\lambda_1=2,\ \lambda_2=1,\ \lambda_3=-1}\)（均为单根，代数重数均为1）。

步骤2：计算对应特征向量

• 对应 \(\lambda_1=2\)：解 \((A_1-2I)x=0\)

  \[A_1-2I = \begin{pmatrix} 0 & -3 & 6 \\ 0 & 1 & -4 \\ 0 & 2 & -5 \end{pmatrix} \xrightarrow{\text{行变换}} \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{pmatrix} \]

  得基础解系：\(\boldsymbol{\xi_1 = (1,0,0)^T}\)，特征向量为 \(k_1\xi_1\)（\(k_1 \neq 0\)）。

• 对应 \(\lambda_2=1\)：解 \((A_1-I)x=0\)

  \[A_1-I = \begin{pmatrix} 1 & -3 & 6 \\ 0 & 2 & -4 \\ 0 & 2 & -4 \end{pmatrix} \xrightarrow{\text{行变换}} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & -2 \\ 0 & 0 & 0 \end{pmatrix} \]

  得基础解系：\(\boldsymbol{\xi_2 = (0,2,1)^T}\)，特征向量为 \(k_2\xi_2\)（\(k_2 \neq 0\)）。

• 对应 \(\lambda_3=-1\)：解 \((A_1+I)x=0\)

  \[A_1+I = \begin{pmatrix} 3 & -3 & 6 \\ 0 & 4 & -4 \\ 0 & 2 & -2 \end{pmatrix} \xrightarrow{\text{行变换}} \begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & -1 \\ 0 & 0 & 0 \end{pmatrix} \]

  得基础解系：\(\boldsymbol{\xi_3 = (-1,1,1)^T}\)，特征向量为 \(k_3\xi_3\)（\(k_3 \neq 0\)）。

步骤3：对角化判断

三个特征值均为单根，对应三个线性无关的特征向量，因此矩阵 \(A_1\) 相似于对角矩阵：

\[A_1 \sim \begin{pmatrix} 2 & & \\ & 1 & \\ & & -1 \end{pmatrix} \]

---

(2) 矩阵 \(A_2 = \begin{pmatrix} 2 & 0 & 1 \\ 0 & 2 & 0 \\ 1 & 0 & 2 \end{pmatrix}\)

步骤1：计算特征值

特征多项式按第二行展开（第二行仅中间元素非零）：

\[\begin{align*} \det(A_2-\lambda I) &= \begin{vmatrix} 2-\lambda & 0 & 1 \\ 0 & 2-\lambda & 0 \\ 1 & 0 & 2-\lambda \end{vmatrix} \\ &= (2-\lambda) \cdot \begin{vmatrix} 2-\lambda & 1 \\ 1 & 2-\lambda \end{vmatrix} \\ &= (2-\lambda) \left[ (2-\lambda)^2 - 1 \right] \\ &= (2-\lambda)(\lambda-1)(\lambda-3) \end{align*} \]

令特征多项式为零，得特征值：\(\boldsymbol{\lambda_1=2,\ \lambda_2=1,\ \lambda_3=3}\)（均为单根，代数重数均为1）。

步骤2：计算对应特征向量

• 对应 \(\lambda_1=2\)：解 \((A_2-2I)x=0\)

  \[A_2-2I = \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix} \xrightarrow{\text{行变换}} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{pmatrix} \]

  得基础解系：\(\boldsymbol{\xi_1 = (0,1,0)^T}\)，特征向量为 \(k_1\xi_1\)（\(k_1 \neq 0\)）。

• 对应 \(\lambda_2=1\)：解 \((A_2-I)x=0\)

  \[A_2-I = \begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \end{pmatrix} \xrightarrow{\text{行变换}} \begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix} \]

  得基础解系：\(\boldsymbol{\xi_2 = (-1,0,1)^T}\)，特征向量为 \(k_2\xi_2\)（\(k_2 \neq 0\)）。

• 对应 \(\lambda_3=3\)：解 \((A_2-3I)x=0\)

  \[A_2-3I = \begin{pmatrix} -1 & 0 & 1 \\ 0 & -1 & 0 \\ 1 & 0 & -1 \end{pmatrix} \xrightarrow{\text{行变换}} \begin{pmatrix} 1 & 0 & -1 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix} \]

  得基础解系：\(\boldsymbol{\xi_3 = (1,0,1)^T}\)，特征向量为 \(k_3\xi_3\)（\(k_3 \neq 0\)）。

步骤3：对角化判断

三个特征值均为单根，对应三个线性无关的特征向量，且 \(A_2\) 是对称矩阵（对称矩阵必可对角化），因此矩阵 \(A_2\) 相似于对角矩阵：

\[A_2 \sim \begin{pmatrix} 2 & & \\ & 1 & \\ & & 3 \end{pmatrix} \]

---

(3) 矩阵 \(A_3 = \begin{pmatrix} 1 & 0 & 0 \\ -1 & 0 & 1 \\ -1 & -1 & 2 \end{pmatrix}\)

步骤1：计算特征值

特征多项式按第一行展开（第一行仅第一个元素非零）：

\[\begin{align*} \det(A_3-\lambda I) &= \begin{vmatrix} 1-\lambda & 0 & 0 \\ -1 & -\lambda & 1 \\ -1 & -1 & 2-\lambda \end{vmatrix} \\ &= (1-\lambda) \cdot \begin{vmatrix} -\lambda & 1 \\ -1 & 2-\lambda \end{vmatrix} \\ &= (1-\lambda) \left[ -\lambda(2-\lambda) + 1 \right] \\ &= (1-\lambda)(\lambda^2 - 2\lambda + 1) \\ &= -(λ-1)^3 \end{align*} \]

令特征多项式为零，得特征值：\(\boldsymbol{\lambda=1}\)（三重根，代数重数为3）。

步骤2：计算对应特征向量

解 \((A_3-I)x=0\)：

\[A_3-I = \begin{pmatrix} 0 & 0 & 0 \\ -1 & -1 & 1 \\ -1 & -1 & 1 \end{pmatrix} \xrightarrow{\text{行变换}} \begin{pmatrix} 1 & 1 & -1 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \]

方程为 \(x_1 + x_2 - x_3 = 0\)，自由变量为 \(x_2, x_3\)，得基础解系：

\[\boldsymbol{\xi_1 = (-1,1,0)^T,\ \xi_2 = (1,0,1)^T} \]

特征向量为 \(k_1\xi_1 + k_2\xi_2\)（\(k_1,k_2\) 不全为零）。

步骤3：对角化判断

特征值 \(\lambda=1\) 的代数重数为3，但几何重数为2（基础解系仅含2个线性无关的向量），几何重数小于代数重数，因此矩阵 \(A_3\) 不能相似于对角矩阵。

---

结果汇总表

矩阵	特征值	对应特征向量（基础解系）	是否可对角化	相似对角矩阵
\(A_1\)	2, 1, -1	\((1,0,0)^T,\ (0,2,1)^T,\ (-1,1,1)^T\)	是	\(\operatorname{diag}(2,1,-1)\)
\(A_2\)	2, 1, 3	\((0,1,0)^T,\ (-1,0,1)^T,\ (1,0,1)^T\)	是	\(\operatorname{diag}(2,1,3)\)
\(A_3\)	1（三重）	\((-1,1,0)^T,\ (1,0,1)^T\)	否	不存在

---

幂法计算主特征值与特征向量 详细解答

规范化幂法迭代公式

依据：无穷范数规范化幂法

\[\begin{cases} u_0 = (1,1,1)^T \quad (\text{初始向量}) \\ v_k = A u_{k-1} \\ \mu_k = \max(v_k) = \|v_k\|_\infty \quad (\text{主特征值近似}) \\ u_k = \frac{v_k}{\mu_k} \quad (\text{规范化特征向量近似}) \end{cases} \]

终止条件：当相邻两次迭代的 \(\mu_k\) 有3位小数稳定时，终止迭代。

---

(1) 矩阵 \(A_1 = \begin{pmatrix} 7 & 3 & -2 \\ 3 & 4 & -1 \\ -2 & -1 & 3 \end{pmatrix}\)

迭代过程

k	\(v_k^T\)	\(\mu_k\)	\(u_k^T\)（规范化）
0	-	-	(1, 1, 1)
1	(8, 6, 0)	8.000	(1.000, 0.750, 0.000)
2	(9.250, 6.000, -2.750)	9.250	(1.000, 0.649, -0.297)
3	(9.540, 5.892, -3.541)	9.540	(1.000, 0.618, -0.371)
4	(9.595, 5.842, -3.731)	9.595	(1.000, 0.609, -0.389)
5	(9.604, 5.824, -3.775)	9.604	(1.000, 0.606, -0.393)
6	(9.605, 5.819, -3.786)	9.605	(1.000, 0.606, -0.394)
7	(9.606, 5.817, -3.788)	9.606	(1.000, 0.606, -0.394)
8	(9.606, 5.817, -3.789)	9.606	(1.000, 0.606, -0.394)

最终结果

• 主特征值：\(\boldsymbol{\lambda_1 \approx 9.606}\)（精确值为 \(6+\sqrt{13} \approx 9.60555\)，误差约 \(4.5 \times 10^{-5}\)）
• 对应特征向量：\(\boldsymbol{x_1 \approx (1, 0.606, -0.394)^T}\)（精确值为 \((1, \frac{6-\sqrt{13}}{3}, \frac{\sqrt{13}-4}{3})^T\)）

---

(2) 矩阵 \(A_2 = \begin{pmatrix} 3 & -4 & 3 \\ -4 & 6 & 3 \\ 3 & 3 & 1 \end{pmatrix}\)

迭代过程

k	\(v_k^T\)	\(\mu_k\)	\(u_k^T\)（规范化）
0	-	-	(1, 1, 1)
1	(2, 5, 7)	7.000	(0.286, 0.714, 1.000)
2	(1.000, 6.143, 4.000)	6.143	(0.163, 1.000, 0.651)
3	(-1.558, 7.302, 4.140)	7.302	(-0.213, 1.000, 0.567)
4	(-2.940, 8.554, 2.927)	8.554	(-0.344, 1.000, 0.342)
5	(-4.005, 8.401, 2.311)	8.401	(-0.477, 1.000, 0.275)
10	(-5.297, 8.856, 1.381)	8.856	(-0.598, 1.000, 0.156)
15	(-5.358, 8.870, 1.340)	8.870	(-0.604, 1.000, 0.151)
16	(-5.359, 8.869, 1.339)	8.869	(-0.604, 1.000, 0.151)
17	(-5.360, 8.870, 1.338)	8.870	(-0.604, 1.000, 0.151)
18	(-5.360, 8.870, 1.338)	8.870	(-0.604, 1.000, 0.151)

最终结果

• 主特征值：\(\boldsymbol{\lambda_1 \approx 8.870}\)（精确值约为8.8699，误差约 \(1 \times 10^{-4}\)）
• 对应特征向量：\(\boldsymbol{x_1 \approx (-0.604, 1, 0.151)^T}\)（或规范化为 \((1, -1.656, -0.250)^T\)）

---

结果验证与说明

1. 收敛性分析：

   • 矩阵 \(A_1\) 的次大特征值为 \(6-\sqrt{13} \approx 2.394\)，收敛比值 \(r = |2.394/9.606| \approx 0.249\)，收敛较快，仅需8次迭代达到3位小数精度。
   • 矩阵 \(A_2\) 的次大特征值约为2.476，收敛比值 \(r \approx |2.476/8.870| \approx 0.279\)，收敛稍慢，需17次迭代达到3位小数精度。

2. 精确值验证：

   • \(A_1\) 的特征多项式为 \(-(\lambda-2)(\lambda^2-12\lambda+23)=0\)，主特征值 \(6+\sqrt{13} \approx 9.60555\)，与迭代结果一致。
   • \(A_2\) 的特征多项式为 \(-\lambda^3+10\lambda^2+7\lambda-151=0\)，主特征值约为8.8699，与迭代结果一致。

---

带原点平移的反幂法求解最接近6的特征值 详细解答

方法原理

依据：带原点平移的反幂法
要求矩阵 \(A\) 最接近给定值 \(p\) 的特征值，构造平移矩阵 \(B = A - pI\)，对 \(B^{-1}\) 应用幂法。若 \(\lambda\) 是 \(A\) 最接近 \(p\) 的特征值，则 \(1/(\lambda-p)\) 是 \(B^{-1}\) 的主特征值。迭代公式为：

\[\begin{cases} \text{解线性方程组 } (A-pI)v_k = u_{k-1} \\ \mu_k = \max(v_k) = \|v_k\|_\infty \\ u_k = \frac{v_k}{\mu_k} \\ \lambda \approx p + \frac{1}{\mu_k} \end{cases} \]

终止条件：特征值近似值有3位小数稳定。

---

问题求解

给定矩阵：

\[A = \begin{pmatrix} 6 & 2 & 1 \\ 2 & 3 & 1 \\ 1 & 1 & 1 \end{pmatrix} \]

取平移量 \(p=6\)，构造平移矩阵：

\[B = A - 6I = \begin{pmatrix} 0 & 2 & 1 \\ 2 & -3 & 1 \\ 1 & 1 & -5 \end{pmatrix} \]

步骤1：B的带部分选主元LU分解

为高效求解线性方程组，先对 \(B\) 进行LU分解：

• 置换矩阵：\(P = \begin{pmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix}\)
• 单位下三角矩阵：\(L = \begin{pmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ 0 & 0.8 & 1 \end{pmatrix}\)
• 上三角矩阵：\(U = \begin{pmatrix} 2 & -3 & 1 \\ 0 & 2.5 & -5.5 \\ 0 & 0 & 5.4 \end{pmatrix}\)

验证：\(PB = LU\)，分解正确。

步骤2：迭代过程

初始向量 \(u_0 = (1,1,1)^T\)，每次迭代通过解 \(Ly=Pu_{k-1}\) 和 \(Uv_k=y\) 得到 \(v_k\)。

k	\(v_k^T\)	\(\mu_k\)	\(u_k^T\)（规范化）	特征值近似 \(\lambda \approx 6+1/\mu_k\)
0	-	-	(1, 1, 1)	-
1	(1.111, 0.444, 0.111)	1.111	(1, 0.400, 0.100)	6.900
2	(0.700, 0.400, 0.200)	0.700	(1, 0.571, 0.286)	7.429
3	(0.804, 0.407, 0.185)	0.804	(1, 0.507, 0.230)	7.244
4	(0.768, 0.406, 0.189)	0.768	(1, 0.529, 0.246)	7.303
5	(0.779, 0.406, 0.188)	0.779	(1, 0.521, 0.241)	7.283
6	(0.775, 0.406, 0.188)	0.775	(1, 0.524, 0.243)	7.290
7	(0.777, 0.406, 0.188)	0.777	(1, 0.523, 0.242)	7.287
8	(0.776, 0.406, 0.188)	0.776	(1, 0.523, 0.242)	7.288
9	(0.776, 0.406, 0.188)	0.776	(1, 0.523, 0.242)	7.288
10	(0.776, 0.406, 0.188)	0.776	(1, 0.523, 0.242)	7.288

步骤3：收敛判断

从第8次迭代开始，特征值近似值的前3位小数（7.288）保持稳定，满足终止条件。

---

最终结果

• 最接近6的特征值：\(\boldsymbol{\lambda \approx 7.288}\)（精确值约为7.2879，误差约 \(1 \times 10^{-4}\)）
• 对应特征向量：\(\boldsymbol{x \approx (1, 0.523, 0.242)^T}\)（或其非零常数倍）

精确值验证

矩阵 \(A\) 的特征多项式为 \(\lambda^3 - 10\lambda^2 + 21\lambda - 9 = 0\)，代入 \(\lambda=7.2879\) 验证：

\[7.2879^3 - 10 \times 7.2879^2 + 21 \times 7.2879 - 9 \approx 0 \]

结果正确。

---

矩阵特征值4对应的特征向量求解

方法一：定义法（直接求解齐次线性方程组）

步骤1：构造特征矩阵

已知矩阵

\[A = \begin{pmatrix} 4 & 0 & 0 \\ 0 & 3 & 1 \\ 0 & 1 & 3 \end{pmatrix} \]

对应特征值 \(\lambda=4\)，构造特征矩阵：

\[A - 4I = \begin{pmatrix} 4-4 & 0 & 0 \\ 0 & 3-4 & 1 \\ 0 & 1 & 3-4 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & -1 & 1 \\ 0 & 1 & -1 \end{pmatrix} \]

步骤2：求解齐次线性方程组 \((A-4I)x=0\)

对系数矩阵进行初等行变换：

\[\begin{pmatrix} 0 & 0 & 0 \\ 0 & -1 & 1 \\ 0 & 1 & -1 \end{pmatrix} \xrightarrow{r_3 + r_2} \begin{pmatrix} 0 & 0 & 0 \\ 0 & -1 & 1 \\ 0 & 0 & 0 \end{pmatrix} \]

得到同解方程组：

\[-x_2 + x_3 = 0 \implies x_2 = x_3 \]

自由变量为 \(x_1\) 和 \(x_3\)，分别取：

• 令 \(x_1=1, x_3=0\)，得基础解系 \(\boldsymbol{\xi_1 = (1, 0, 0)^T}\)
• 令 \(x_1=0, x_3=1\)，得基础解系 \(\boldsymbol{\xi_2 = (0, 1, 1)^T}\)

步骤3：写出所有特征向量

矩阵 \(A\) 对应于特征值4的所有特征向量为：

\[\boldsymbol{k_1 \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + k_2 \begin{pmatrix} 0 \\ 1 \\ 1 \end{pmatrix}} \]

其中 \(k_1, k_2\) 为不全为零的任意常数。

---

方法二：带原点平移的反幂法验证

原理说明

当平移量 \(p\) 等于特征值 \(\lambda\) 时，矩阵 \(A-pI\) 奇异，理论上反幂法一步即可收敛到特征向量。实际计算中，由于舍入误差，矩阵不会严格奇异，迭代会快速收敛。

迭代过程

取初始向量 \(u_0=(1,1,1)^T\)，平移量 \(p=4.0001\)（避免严格奇异）：

1. 解线性方程组 \((A-4.0001I)v_1 = u_0\)，得 \(v_1 \approx (10000, 0.5, 0.5)^T\)
2. 规范化得 \(u_1 \approx (1, 0.00005, 0.00005)^T\)
3. 第二次迭代后，\(u_2 \approx (1, 0, 0)^T\)，已收敛到特征向量 \(\xi_1\)

若取初始向量 \(u_0=(0,1,1)^T\)，则一次迭代即可得到特征向量 \(\xi_2\)。

---

结果验证

将特征向量代入特征方程 \(Ax = \lambda x\) 验证：

• 对 \(\xi_1=(1,0,0)^T\)：\(A\xi_1 = (4,0,0)^T = 4\xi_1\)，成立
• 对 \(\xi_2=(0,1,1)^T\)：\(A\xi_2 = (0,4,4)^T = 4\xi_2\)，成立

补充说明：特征值4是矩阵 \(A\) 的二重特征值，其代数重数和几何重数均为2，因此有两个线性无关的特征向量，矩阵 \(A\) 可对角化。

---

带原点平移的反幂法求解特征值4对应的特征向量

方法原理

依据：带原点平移的反幂法
当已知特征值的近似值 \(p\) 时，构造平移矩阵 \(B = A - pI\)，对 \(B^{-1}\) 应用幂法。若 \(p\) 非常接近真实特征值 \(\lambda\)，则 \(1/(\lambda-p)\) 是 \(B^{-1}\) 的绝对占优主特征值，迭代将以极快的速度收敛到对应特征向量。

对于本题，已知特征值 \(\lambda=4\)，取平移量 \(p=4.0001\)（微小偏移避免矩阵严格奇异，同时保证极快收敛）。

---

迭代求解过程

给定矩阵：

\[A = \begin{pmatrix} 4 & 0 & 0 \\ 0 & 3 & 1 \\ 0 & 1 & 3 \end{pmatrix} \]

构造平移矩阵：

\[B = A - 4.0001I = \begin{pmatrix} -0.0001 & 0 & 0 \\ 0 & -1.0001 & 1 \\ 0 & 1 & -1.0001 \end{pmatrix} \]

步骤1：B的LU分解

对 \(B\) 进行带部分选主元的LU分解：

• 单位下三角矩阵：\(L = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -0.9999 & 1 \end{pmatrix}\)
• 上三角矩阵：\(U = \begin{pmatrix} -0.0001 & 0 & 0 \\ 0 & -1.0001 & 1 \\ 0 & 0 & -0.0002 \end{pmatrix}\)

验证：\(LU = B\)，分解正确。

步骤2：不同初始向量的迭代结果

由于特征值4是二重特征值，对应的特征子空间是二维的，反幂法的收敛结果取决于初始向量，会收敛到该子空间中的某个向量。

情况1：初始向量 \(u_0=(1,0,0)^T\)

1. 解下三角方程组 \(Ly = u_0\)，得 \(y=(1,0,0)^T\)
2. 解上三角方程组 \(Uv_1 = y\)，得 \(v_1=(-10000, 0, 0)^T\)
3. 规范化：\(u_1 = \frac{v_1}{\max(v_1)} = (1, 0, 0)^T\)

一次迭代即收敛，得到特征向量 \(\boldsymbol{\xi_1=(1,0,0)^T}\)。

情况2：初始向量 \(u_0=(0,1,0)^T\)

1. 解下三角方程组 \(Ly = u_0\)，得 \(y=(0,1,0.9999)^T\)
2. 解上三角方程组 \(Uv_1 = y\)，得 \(v_1=(0, -4999, -4999.5)^T\)
3. 规范化：\(u_1 = \frac{v_1}{\max(v_1)} \approx (0, 1, 1.0001)^T \approx (0,1,1)^T\)

一次迭代即收敛，得到特征向量 \(\boldsymbol{\xi_2=(0,1,1)^T}\)。

情况3：初始向量 \(u_0=(1,1,1)^T\)

1. 解下三角方程组 \(Ly = u_0\)，得 \(y=(1,1,1.9999)^T\)
2. 解上三角方程组 \(Uv_1 = y\)，得 \(v_1=(-10000, -9999.5, -9999.5)^T\)
3. 规范化：\(u_1 = \frac{v_1}{\max(v_1)} \approx (1, 0.99995, 0.99995)^T \approx (1,1,1)^T\)

一次迭代即收敛，得到特征向量 \(\boldsymbol{\xi_3=(1,1,1)^T}\)，它是 \(\xi_1\) 和 \(\xi_2\) 的线性组合。

---

最终结果

矩阵 \(A\) 对应于特征值4的所有特征向量为：

\[\boldsymbol{k_1 \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + k_2 \begin{pmatrix} 0 \\ 1 \\ 1 \end{pmatrix}} \]

其中 \(k_1, k_2\) 为不全为零的任意常数。

结果验证

将特征向量代入特征方程 \(Ax = 4x\) 验证：

• \(A\xi_1 = (4,0,0)^T = 4\xi_1\)，成立
• \(A\xi_2 = (0,4,4)^T = 4\xi_2\)，成立
• \(A\xi_3 = (4,4,4)^T = 4\xi_3\)，成立