📑 目录

1. 引言：本征值问题的物理意义
2. 预备知识：相似变换与正交矩阵
3. Jacobi 方法 —— 逐个消灭非对角元
4. Householder 变换 —— 批量消元成三对角
5. QR 算法 —— 二十世纪十大算法之一
6. Power 方法（幂法）与反幂法
7. 奇异值分解（SVD）
8. Krylov 子空间方法
9. Lanczos 迭代
10. Arnoldi 迭代
11. Schrödinger 方程的对角化求解
12. DMRG：密度矩阵重整化群
13. 算法对比与总结

本质上 主线就是 把物理问题 变成 矩阵的本征值问题，然后利用正交变换把矩阵“变得简单但保证其本征值不变”，再利用迭代法吧本征值计算出来！

1. 引言：本征值问题的物理意义

1.1 什么是本征值问题？

在量子力学中，任何可以实验测量的物理量都对应一个算符，该算符的本征值就是测量可能得到的结果。薛定谔方程本质上是一个本征值问题：

说人话：给定一个矩阵 $H$（哈密顿量），要找那些满足「矩阵作用上去等于标量相乘」的向量 $|\psi\rangle$ 和对应的数 $E$。

1.2 量子多体问题中本征值问题的来源

如果系统里不是一个粒子，而是很多粒子互相作用，比如：

• 原子里的多个电子
• 原子核里的多个核子
• 量子多体凝聚态系统

那就没法直接写成一个简单单粒子问题了。通常会走这样一条路：

1. 先做一个平均场近似，得到单粒子基
2. 用这些单粒子基构造多体基态/激发态
3. 把总波函数展开成很多个组态的线性组合
4. 把哈密顿量写成这个基下的矩阵
5. 解这个矩阵的本征值问题

也就是：

代入 $H|\Psi\rangle = E|\Psi\rangle$，得到

即矩阵形式：

这里 $H_{IJ}=\langle \Phi_I|H|\Phi_J\rangle$。

描述量子物理和计算化学中一个经典的套路：如果利用计算机去求解有多个粒子互相影响的复杂系统 （比如原子里的多个电子，或者原子核里的多个质子与中子）

“组态相互作用”（Configuration Interaction） 或者 精确对角化法

独立粒子图像    → 残余相互作用  → 组态混合    → 对角化
(Mean Field)    (Residual)    (CI)        (Diagonalization)

核心链条：

💡 物理直觉：组态空间（Fock 空间）是不同粒子数 Hilbert 空间的直和。组态混合（Configuration Mixing）——多体关联/纠缠——最终都归结为「在一个巨大矩阵中求本征值」。排列组合使得组态空间极其巨大，例如 $N$ 个粒子分布在 $M$ 个轨道上，组态数呈指数增长。

1.3 本征值问题的数学本质

对于 $n \times n$ 矩阵 $A$，求标量 $\lambda$ 和非零向量 $\mathbf{x}$ 使得：

等价于解特征方程：

数值例子（2×2 矩阵）：

得 $\lambda_1 = 5,\ \lambda_2 = 1$。

💡 但是！对于 $n > 4$，直接展开行列式求根是数值灾难——这就是我们需要各种迭代算法的原因。

特征多项式数值上非常不划算！

多项式系数会涉及大量相消，高次根对于系数变化非常敏感，直接展开很不稳定，对大矩阵完全不现实！！！

2. 预备知识：相似变换与正交矩阵

2.1 相似变换保持本征值不变

如果 $B = S^{-1}AS$，则 $A$ 和 $B$ 有相同的本征值。

说人话：用矩阵 $S$ 对 $A$ 做「换个坐标系」的操作，本征值不变，只是本征矢跟着变。

验证（数值例子）：

$A$ 的本征值：$\det(A-\lambda I) = (4-\lambda)(3-\lambda)-2=0 \Rightarrow \lambda=2,5$ $B$ 的本征值：对角元即为 $2, 5$（上三角矩阵的本征值就是对对角元）。✓

2.2 实对称矩阵 → 实正交相似对角化

如果 A 是 实对称矩阵，那么有谱定理：

• 本征值都是实数 借助定义直接去证明
• 不同本征值对应的本征矢正交 利用对称矩阵 $A=A^T$的性质，在内积运算中“转移”矩阵 A 的作用
• 存在正交矩阵 Q ，使得

其中这个 $\Lambda$ 是对角矩阵， 对角元就是本征值！

证明：使用数学归纳法，每次找一个特征向量，把空间降维，并保证降维后的子矩阵依然对称

为什么 PCA（主成分分析）、图神经网络（GNN）、量子力学、流形学习等算法都死死抱住实对称矩阵不放？

1. 物理/几何意义明确：对称矩阵代表二次型 $x^TAx$（如能量、方差）。谱定理告诉我们，总可以通过旋转坐标系（正交矩阵 Q ），消除交叉项，找到事物的“主成分”（对角元 $Λ$）。
2. 计算极其稳定：一般矩阵对角化 $A=PΛP^{−1}$ 中，P可能是病态的（条件数很大），求逆会放大误差。而对称矩阵的 Q是正交矩阵，$Q^{−1}=Q^T$，转置计算不仅没有误差，而且速度极快。
3. 降维与截断：在算法中，我们常常只取前 k 个最大的特征值（因为对称矩阵的特征值是实数，可以比大小！）。由于特征向量正交，这前 k个方向保留了数据中最大的“能量”或“信息量”，且彼此不相关。

简而言之，实对称矩阵的谱定理保证了：复杂的世界（一般矩阵）可以通过旋转（正交变换），被拆解为几个相互独立、互不干扰的简单方向（对角矩阵）。 这就是它在科学计算中拥有统治地位的原因。

关键定理：若 $A$ 是实对称矩阵，则存在实正交矩阵 $S$（即 $S^T S = I$），使得 $S^T A S$ 为对角矩阵，对角元就是本征值。

3. Jacobi 方法 —— 逐个消灭非对角元

3.1 核心思想

一次干掉一个最大的非对角元素，用 Givens 旋转矩阵做相似变换。迭代到所有非对角元都接近零。

最直观的对称矩阵对角化的方法！

每次只处理一个非对角元 $a_{pq}$ 借助一个平面旋转将其消掉，如果把矩阵看成一个房间， Jacobi 就像每次只把一个歪掉的方向修正！

人话版：就像拧魔方，每次只转一个面，转完检查还有哪里不对，再转下一面，直到六面全对齐。

3.2 Givens 旋转矩阵

Givens 旋转 $P(p,q,\theta)$ 只在第 $p,q$ 行/列上做旋转，其他位置保持单位矩阵：

$P$ 是正交矩阵（$P^T P = I$）

对称矩阵 做这样的相似变换：

可以保持对称性和本征值不变！！！

3.3 消灭非对角元 $a_{pq}$

对于对称矩阵 $A$，做相似变换 $A' = P^T A P$：

变换后对应的非对角项 $a'_{pq} = a'_{qp}$ 变为（选 $\theta$ 使得这个位置 = 0）：

选择较小的根确保 $|\theta| \leq \pi/4$。

具体数值例子： $$ A = \begin{bmatrix} 3 & 1 \ 1 & 3 \end{bmatrix},\quad a{pq}=1,\ a{pp}=a_{qq}=3 $$

一次旋转直接对角化！本征值 $4, 2$。

3.4 算法流程

1. 找一个大的非对角元 $a_{pq}$
2. 算旋转角 $\theta$
3. 构造旋转矩阵 $P$
4. 做 $A \leftarrow P^TAP$
5. 重复，直到所有非对角元都足够小

最后对角元就是本征值。

1. 找到 |a_pq| 最大的非对角元 (p≠q)
2. 计算旋转角 θ：tan(2θ) = 2a_pq / (a_pp - a_qq)
3. 构造 Givens 旋转矩阵 P
4. 做相似变换：A ← P^T A P
5. 重复直到 max|非对角元| < ε
6. 对角元即为本征值

3.5 Jacobi的优缺点

优点：思路直观 对称矩阵上作用很稳定 适合进行并行运算

缺点：每次只能消去一个非对角元，效率并不高！ 对大矩阵不如QR快

适合 较小规模稠密矩阵，或者需要进行并行化的场景

3.6 并行 Jacobi 方法

怎么变？ 用旋转（Rotation）。 想象你在二维平面上旋转坐标轴。比如你选中第1行和第2行交叉的那个数 a12，做一次旋转，就能把这个 a12 精准地变成0。

但是，当你把(1,2)位置变成0后，刚才已经消成0的(1,3)位置可能又变成非0了（牵一发而动全身）。所以你得不停地扫荡整个矩阵，反复消。

既然要消很多次，那我同时消 a12a12 和 a34a34 行不行？答案是：可以，但有严格条件。

并行加速，我们必须把所有的旋转操作分组。同一组里的旋转，涉及的行号列号必须完全不重叠（不然会冲突，比如想要同时消去 12和13位置，都涉及到了第一行，改第一行的时候，会影响另一个操作的数据；如果是12和34，完全不重叠，互不影响，可以同时消去）

标准结论是： 对 $n$ 个变量的完全配对，可以分成大约 $n-1$ 轮，每轮里有若干个互不冲突的旋转，能并行执行。

对于一个 n阶矩阵，总共有 n(n−1)/2个非对角元素需要消（也就是所有配对）。

• 因为每次并行旋转需要 n/2 对（比如8个变量，一次并行做4对，正好把1~8全用上）。
• 总对数 ÷ 每轮并行对数 = n-1 轮。

所以，标准结论是 n−1 轮。

对于 $n=2m$，可以分出 $2m-1$ 个选择集，每个集有 $m$ 个互不冲突的旋转（即旋转指标全不同），可以并行分配。

$n=8$ 的分组示例（7 组，每组 4 个旋转，可分配给 4 个处理器）：

💡 算法特点：简单直观，但每次只消一个元素，收敛较慢。

收敛较慢主要是 只消一个非对角，同时之前消去的还可能会被再次污染回来，相比之下，现在的算法（分而治之，QR算法）都是一次性处理大块数据！

适合小型稠密矩阵。Jacobi 天然适合并行！ 极其稳定 计算量完全可以预测

4. Householder 变换 —— 批量消元成三对角

4.1 核心思想

我们知道理论上 实对称矩阵一定是可以被正交对角化的 但直接一步求出来用于对角化的正交矩阵Q是很难，所以数值算法通常选择线借助 Householder变换将原本实对称矩阵A变成 T（三对角矩阵） 然后再利用 QR算法 算出来A

三对角：只保留了主对角线和上下两条次对角线！！！

元素数量从 $O(n^{2}) -> O(n)$ 后续的迭代算法会快很多！

用一个正交反射变换 一次消灭一整列的非对角元，把对称矩阵变成三对角形式。

对比 Jacobi：

4.2 Householder 反射矩阵

给定单位向量 $\mathbf{w}$（$\mathbf{w}^T\mathbf{w}=1$），构造反射矩阵：

展开为矩阵元形式（设 $\mathbf{w} = (w_1, w_2, \ldots, w_n)^T$）：

几何意义：$P$ 将任意向量 $\mathbf{v}$ 变成它关于超平面（法向量为 $\mathbf{w}$）的镜像反射

任意向量 $v$ 可以分解成两部分：

其中 $v_\parallel$ 是沿 $w$ 方向的分量：

而垂直于 $w$ 的部分是：

Householder 作用在 $v$ 上：

代入分解：

且

所以

也就是说：

• 垂直于 $w$ 的分量保持不变；
• 沿 $w$ 的分量反号。

这正是关于法向量为 $w$ 的超平面的镜像反射。

性质：$P$ 是对称的、正交的（$P = P^T = P^{-1}$）

4.3 构造 $\mathbf{w}$ 的方法

目标：给定向量 $\mathbf{v}$，找一个反射矩阵 $P$ 使 $P\mathbf{v} = \alpha\mathbf{e}_1$（即 $P\mathbf{v}$ 只有第一个分量非零）

构造步骤：

首先就是 P 是正交矩阵，所以它保持了向量长度！ $$ \alpha = \pm|\mathbf{v}| = \pm\sqrt{v1^2 + v2^2 + \cdots + v_n^2} $$

符号选取：取 $\alpha$ 与 $v_1$ 反号，避免相近数相减导致精度损失。

构造 Householder 向量时，会用到：

$$u=v-\alpha e_1.$$

如果 $\alpha$ 与 $v_1$ 同号，且 $|\alpha|\approx |v_1|$，那么第一个分量

$$v_1-\alpha$$

会发生相近数相减，造成严重精度损失。

所以数值上通常选：

$$\alpha=-\operatorname{sign}(v_1)\|v\|.$$

如果 $v_1>0$，取

$$\alpha=-\|v\|.$$

如果 $v_1<0$，取

$$\alpha=+\|v\|.$$

这样

$$v_1-\alpha$$

的绝对值会变大，数值更稳定。

推导 $\omega$ 的公式

令

然后令

于是

我们要证明：

先计算

展开：

其中

因为 $|\alpha|=\|v\|$，而

所以

再计算

展开：

即

所以

因此

于是

可以验证 $\mathbf{w}^T\mathbf{w} = 1$，且 $P = I - 2\mathbf{w}\mathbf{w}^T$ 满足 $P\mathbf{v} = \alpha\mathbf{e}_1$。

数值例子：设 $\mathbf{v} = (2, 2, 1)^T$

1. $\|\mathbf{v}\| = \sqrt{4+4+1} = 3$
2. $v_1 = 2 > 0$，取 $\alpha = -3$（与 $v_1$ 反号）
3. $\mathbf{v} - \alpha\mathbf{e}_1 = (2, 2, 1)^T - (-3, 0, 0)^T = (5, 2, 1)^T$
4. $\mathbf{w}^T\mathbf{v} = \frac{v_1 - \alpha}{2} = \frac{2 - (-3)}{2} = 2.5$
5. $\mathbf{w} = \frac{(5, 2, 1)^T}{\sqrt{2 \times (-3)^2 - \cdots}}$（具体计算略）
6. 最终 $P\mathbf{v} = (-3, 0, 0)^T$ ✓

4.4 Householder 三对角化算法

现在考虑一个 $n\times n$ 实对称矩阵：

我们希望通过一系列正交相似变换：

把它变成三对角矩阵。

因为 $H_k$ 正交，所以本征值不变。

对 $n \times n$ 对称矩阵 $A$，经过 $n-2$ 次 Householder 变换：

for k = 1 to n-2:
    取 A 的第 k 列第 k+1 到 n 行作为向量 v
    构造 Householder 矩阵 P_k（大小为 n-k）
    用 P_k 对 A 做相似变换（扩充为 n×n 的块对角形式）
end

效果：对称 $A$ → 三对角矩阵 $T$

💡 如果 $A$ 不对称，则变换成 上 Hessenberg 矩阵（下三角除次对角线外全为零）

详细展开过程

设

因为 $A$ 对称，所以第一列下面是：

我们希望保留 $a_{21}$，把 $a_{31},\ldots,a_{n1}$ 全部消掉，使第一列变成：

取子向量：

构造一个 $(n-1)\times(n-1)$ 的 Householder 矩阵 $P_1$，使得：

然后把它扩展成 $n\times n$ 矩阵：

做相似变换： $$ A^{(2)}=H1^TAH1. $$ 由于 $H_1=H_1^T$，也可写成： $$ A^{(2)}=H1AH1. $$ 这会把第一列的第 3 到第 $n$ 个元素消掉。

由于变换是对称相似变换，矩阵仍然对称，因此第一行对应的第 3 到第 $n$ 个元素也同时为 0。

于是第一行/列已经满足三对角结构。

4.14 第二步及以后

第二步不再动第一行第一列。

我们取第 2 列从第 3 行到第 $n$ 行的子向量：

构造 $(n-2)\times(n-2)$ 的 Householder 矩阵 $P_2$，扩展为：

然后做：

这样第 2 列中第 4 行以下的元素都被消掉。

继续这个过程，直到第 $n-2$ 步。

n-2是因为 每k列只需要消掉 第k+2行 以下的元素

当 k = n - 2的时候，最后两列已经自然满足三对角

还有核心就是 后边的变换并不会破坏前面的0，就是每次构造矩阵的时候，我们把左上角构造成 单元矩阵 得以 保证已经做好的三对角零结构不会被破坏！

4.5 三对角矩阵的求解

• 如果某个 $\beta_i = 0$，矩阵可以分裂成两个更小的独立块分别对角化
• 然后用 Jacobi 或 QR 算法对角化三对角矩阵
• 计算量：Householder 三对角化 ≈ $2n^3/3$（不含本征矢）或 $4n^3/3$（含本征矢）

5. QR 算法 —— 二十世纪十大算法之一

核心：QR分解

QR分解 就是把一个矩阵 写成

Q：正交矩阵

R：上三角矩阵

如果 A 是 n ✖️ n 非奇异矩阵 则通常可以得到完整的 QR分解

5.1 核心迭代格式

QR 算法是计算矩阵全部本征值的最重要算法。

基本迭代：

从

开始，每一步做 QR 分解：

然后反过来相乘：

看起来只是把乘法顺序反过来，但关键是它等价于一个正交相似变换。

因为

左乘 $(Q^{(k)})^T$：

所以

因此

和

相似，且由于 $Q^{(k)}$ 正交，所以是正交相似。

因此每一步都不改变本征值。

收敛性：

$A^{(k)}$ 趋于上三角矩阵（Schur 形式），对角元就是本征值。

严格证明比较长，但核心直觉可以从幂迭代理解。

设矩阵可以对角化：

其中

并假设本征值大小满足：

幂法中：

因为

所以高次幂会越来越强调最大本征值对应的方向。

QR 迭代可以理解成一种“同时幂法”：

• 幂法追踪一个向量；
• QR 迭代追踪一组正交基；
• 每一步 QR 分解相当于对这组基重新正交化。

在不断迭代中，正交基逐渐对齐到不变子空间，因此矩阵在这组基下逐渐变成上三角形式，也就是 Schur 形式。

对于实对称矩阵，上三角 Schur 形式进一步退化为对角矩阵，因为实对称矩阵被正交相似变换后仍然对称，而一个既对称又上三角的矩阵只能是对角矩阵。

人话版：

1. 把矩阵拆成「正交×上三角」
2. 把乘积顺序反过来再乘一次
3. 重复，矩阵就自动趋近对角形

5.2 QR 分解 —— Givens 旋转法

用 $n-1$ 次 Givens 旋转把矩阵消成上三角：

第 $k$ 次旋转：选 $P_k$ 使 $(P_k A_{k-1}^{(1)})$ 的第 $(k, k-1)$ 位置为 0。

取 $\cos\theta, \sin\theta$ 满足：

给定一个向量的两个分量： $$ \begin{bmatrix} a\b \end{bmatrix}, $$ 我们想构造旋转矩阵： $$ G= \begin{bmatrix} c&s\ -s&c \end{bmatrix}, $$ 使得： $$ G \begin{bmatrix} a\b \end{bmatrix} = \begin{bmatrix} r\0 \end{bmatrix}. $$

推导 $c,s$

计算： $$ \begin{bmatrix} c&s\ -s&c \end{bmatrix} \begin{bmatrix} a\b \end{bmatrix} = \begin{bmatrix} ca+sb\ -sa+cb \end{bmatrix}. $$ 我们希望第二个分量为 0： $$ -sa+cb=0. $$ 即 $$ cb=sa. $$ 一个自然选择是： $$ c=\frac{a}{r},\qquad s=\frac{b}{r}, $$ 其中 $$ r=\sqrt{a^2+b^2}. $$ 验证第二个分量： $$ -sa+cb = -\frac{b}{r}a+\frac{a}{r}b = 0. $$ 第一个分量： $$ ca+sb = \frac{a}{r}a+\frac{b}{r}b = \frac{a^2+b^2}{r} = r. $$ 所以 $$ G \begin{bmatrix} a\b \end{bmatrix} = \begin{bmatrix} r\0 \end{bmatrix}. $$

数值例子（消去 $a_{21}$）：

第一列 $(3, 4)^T$，构造 Givens 旋转消去第 2 个分量：

$n-1$ 次旋转后得到 $R$，所有旋转矩阵的转置之积为 $Q$：

因为

所以

因此 $$ Q=P^T $$ 与 Jacobi 旋转的区别：QR 分解中的 Givens 旋转换不是相似变换（$Pk A$ 而非 $Pk^T A P_k$），目的是消成上三角而非对角化。

5.3 Householder 做 QR 分解

Givens 每次只消一个元素。Householder 每次可以消掉一整列下方的元素。

给定矩阵 $A$，第一步取第一列：

构造 Householder 矩阵 $P_1$，使得：

于是

的第一列下面全部为零。

第二步只处理右下角子矩阵，构造：

继续做，最终得到：

因为每个 $P_i$ 都正交，所以它们乘积也正交。

令

则

Householder QR 对稠密矩阵通常比 Givens 更高效，因为它每一步消一整列，而不是一个元素一个元素消。

用 $n-1$ 次 Householder 反射做 QR 分解：

💡 Householder 做 QR 分解比 Givens 旋转效率更高。

5.4 加速收敛：Shifted QR（带位移的 QR）

基础 QR 迭代：

如果本征值很接近，收敛很慢。

于是引入位移 $\sigma_k$：

然后

由

得：

于是

展开：

因为

所以

因此带位移 QR 仍然是正交相似变换，本征值不变。

设某个本征值 $\lambda_j$ 是我们想让它快速收敛出来的对象。

对

而言，本征值变成：

如果 $\sigma$ 接近某个本征值 $\lambda_j$，那么

会很小。

QR 迭代里，某些次对角元的衰减速度与本征值间隔有关。选择好的 $\sigma$ 可以把某个方向迅速分离出来，使矩阵右下角的次对角元素快速趋近 0。

一旦某个次对角元满足：

矩阵就近似分裂成：

如果矩阵是对称三对角，则更干净地分裂成两个独立块：

这叫 deflation。

Wilkinson shift

对于对称三对角矩阵，常取右下角 $2\times2$ 子块： $$ \begin{bmatrix} a&b\ b&d \end{bmatrix}. $$ 它的本征值为： $$ \mu_{\pm} = \frac{a+d}{2} \pm \sqrt{\left(\frac{a-d}{2}\right)^2+b^2}. $$ Wilkinson shift 取其中更接近 $d$ 的那个。

一种数值稳定写法是： $$ \delta=\frac{a-d}{2}, $$ 这个公式避免了相近数相减。

选 $\sigma$ 的策略：取 $A^{(k)}$ 右下角 $2 \times 2$ 子矩阵的本征值中靠近 $A_{nn}$ 的那个。

数值例子：

右下角 $2\times2$ 子矩阵： $$\begin{bmatrix} a{n-1,n-1}^{(k)} & a{n-1,n}^{(k)} \ a{n,n-1}^{(k)} & a{n,n}^{(k)} \end{bmatrix}$$

设其本征值为 $\mu_1, \mu_2$（例如 $1/2 \pm \sqrt{3}/2$），选 $\sigma = \mu$ 中更接近 $a_{nn}$ 的那个。

5.5 Schur 分解

对任意 $n \times n$ 方阵复 $A$，存在酉矩阵 $Q$ 和上三角矩阵 $U$ 使得：

$U$ 的对角元就是 $A$ 的全部本征值。

$Q^\dagger Q = Q Q^\dagger = I$。

💡 QR 算法本质上就是在计算 Schur 分解。对于不对称矩阵也可以使用。

任何复矩阵A，都可以通过“旋转坐标轴” （酉变换），变成一个“上三角矩阵”（对角线是特征值，对角线下面全是0）

证明思路用数学归纳法。

对于任意复矩阵 $A$，至少存在一个本征值 $\lambda_1$ 和对应单位本征向量 $q_1$：

把 $q_1$ 扩充成一组标准正交基，需要n个相互垂直的单位向量，把这一堆列向量拼成矩阵 Q1。因为每一列都是单位长且互相垂直，所以 Q1 是酉矩阵（Unitary），也就是 $Q_1^{†}Q_1=I$（旋转坐标轴不改变长度）：

因为 $Q_1$ 酉：

考察：

它的第一列是：

1、矩阵乘法看列：$Q1^{\dagger}A Q1$ 这个矩阵，乘上 第一列标准基 $e_1$ 就等于取出其第一列，而Q1e1其实就是 Q1矩阵的第一列，就是q1！

2、带入特征方程

3、核心就是 $Q_1^{\dagger}q_1$ 等于多少？

展开来写就是

进一步，利用标准正交基的性质，得到上述式子其实就是 e1

于是新矩阵 $Q_1^\dagger A Q_1$ 的第一列，就是 $\lambda_1$ 开头，下边都是0！

所以矩阵形状为：

然后对右下角 $A_2$ 继续做同样操作。归纳下去，最终得到上三角矩阵。

QR 算法本质上就是在数值逼近这个 Schur 分解。

5.6 广义本征值问题

广义本征值问题是：

假设：

因为 $B$ 对称正定，可以做 Cholesky 分解：

原方程：

令

则

代回：

左乘 $L^{-1}$：

定义

则：

而且 $C$ 是对称的：

因此广义对称本征值问题被化成了标准对称本征值问题。

求出 $y$ 后，原问题的本征向量是： $$ x=L^{-T}y. $$ 求解 $A\mathbf{x} = \lambda B\mathbf{x}$，其中 $A, B$ 对称且 $B$ 正定。

处理方法：做 Cholesky 分解 $B = LL^T$，化为标准本征值问题：

$C$ 是对称矩阵，其本征值即原问题的本征值，但本征矢需要回代：

💡 $B^{-1}A$ 不对称，不能直接做。通过 Cholesky 变换保对称性。

6. Power 方法（幂法）与反幂法

6.1 幂法 —— 求最大本征值

核心思想：设本征值排列 $|\lambda_1| > |\lambda_2| \geq \cdots \geq |\lambda_n|$，任意初始向量 $\mathbf{x}^{(0)}$ 在本征矢基下展开：

反复左乘 $A$：

当 $k$ 很大时，$\lambda1^k$ 主导（假设 $c1 \neq 0$）：

说人话：不断被矩阵乘，向量方向会被拽到「最大本征值对应的本征矢方向」上。

6.2 算法步骤（带归一化）

如果 $|\lambda_1|>1$，那么

会指数增长，导致数值溢出。

如果 $|\lambda_1|<1$，则会指数衰减，导致下溢。

所以每步都做归一化：

x_0 = 任意非零向量
for k = 1, 2, ...:
    y_k = A x_{k-1}           ← 矩阵乘向量
    λ_k = max(|y_k|) 的分量   ← 近似最大本征值
    x_k = y_k / λ_k           ← 归一化，防止溢出
end

数值例子（3×3 矩阵）：

精确最大本征值：$\lambda_{\max} = 2 + \sqrt{2} \approx 3.414$。

6.3 对称矩阵的幂法

对于归一化向量 $x$，常用 Rayleigh 商：

如果 $x$ 已归一化，则

如果 $x=v_i$，则

因为

所以

所以当 $x^{(k)}$ 越接近本征向量时，Rayleigh 商越接近对应本征值。

对实对称矩阵，幂法的收敛速度快得多——类似于变分法求基态能量：

当 $\lambda_2/\lambda_1$ 很小时收敛极快。

设 $A$ 对称，有标准正交本征向量：

令当前归一化向量为：

且

Rayleigh 商：

代入：

于是：

由于正交归一：

所以：

与 $\lambda_1$ 的误差：

因为

所以：

其中 $i=1$ 项为 0，因此：

注意误差由 $c_i^2$ 控制。

而幂法中：