本征值问题

2026-07-22

📑 目录

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

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

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

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

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

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

也就是：

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

$$ \sum_J H_{IJ} c_J = E c_I $$ 即矩阵形式：

$$ H\mathbf{c} = E\mathbf{c} $$ 这里 $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}$ 使得：

$$A\mathbf{x} = \lambda\mathbf{x}$$

等价于解特征方程：

$$\det(A - \lambda I) = 0$$

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

$$ A = \begin{bmatrix} 3 & 2 \\ 2 & 3 \end{bmatrix} $$

$$\det\begin{bmatrix} 3-\lambda & 2 \\ 2 & 3-\lambda \end{bmatrix} = (3-\lambda)^2 - 4 = \lambda^2 - 6\lambda + 5 = 0$$

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

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

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

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

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

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

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

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

验证（数值例子）：

$$ A = \begin{bmatrix} 4 & 1 \\ 2 & 3 \end{bmatrix},\quad S = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix},\quad S^{-1} = \begin{bmatrix} 1 & -1 \\ 0 & 1 \end{bmatrix} $$

$$ B = S^{-1}AS = \begin{bmatrix} 1 & -1 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} 4 & 1 \\ 2 & 3 \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 2 & 0 \\ 2 & 5 \end{bmatrix} $$

$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 ，使得

$$ Q^{T} A Q = \Lambda $$ 其中这个 $\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 = \begin{bmatrix} 1 & & & & \\ & \ddots & & & \\ & & \cos\theta & \cdots & -\sin\theta \\ & & \vdots & 1 & \vdots \\ & & \sin\theta & \cdots & \cos\theta \\ & & & & & \ddots & \\ & & & & & & 1 \end{bmatrix} \quad \text{——} \quad P_{pp}=P_{qq}=\cos\theta,\ P_{pq}=-P_{qp}=-\sin\theta $$

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

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

$$ A' = P^{T} A P $$ 可以保持对称性和本征值不变！！！

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

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

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

$$ \tan 2\theta = \frac{2a_{pq}}{a_{pp} - a_{qq}} $$ 选择较小的根确保 $|\theta| \leq \pi/4$。

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

$$ \tan 2\theta = \frac{2 \times 1}{3-3} = \infty \quad\Rightarrow\quad 2\theta = \frac{\pi}{2},\ \theta = \frac{\pi}{4} $$

$$ \cos\theta = \sin\theta = \frac{1}{\sqrt{2}},\quad P = \begin{bmatrix} 1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{bmatrix} $$

$$ A' = P^T A P = \begin{bmatrix} 4 & 0 \\ 0 & 2 \end{bmatrix} $$

一次旋转直接对角化！本征值 $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$），构造反射矩阵：

$$ P = I - 2\mathbf{w}\mathbf{w}^T $$

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

$$ P_{ij} = \delta_{ij} - 2w_i w_j = \begin{cases} 1 - 2w_i^2, & i=j \\ -2w_i w_j, & i \neq j \end{cases} $$

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

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

$$ v=v_\parallel+v_\perp. $$ 其中 $v_\parallel$ 是沿 $w$ 方向的分量：

$$ v_\parallel=w(w^Tv), $$ 而垂直于 $w$ 的部分是：

$$ v_\perp=v-w(w^Tv). $$ Householder 作用在 $v$ 上：

$$ Pv=(I-2ww^T)v=v-2w(w^Tv). $$ 代入分解：

$$ v=v_\parallel+v_\perp, $$ 且

$$ v_\parallel=w(w^Tv), $$ 所以

$$ Pv=v_\parallel+v_\perp-2v_\parallel =v_\perp-v_\parallel. $$ 也就是说：

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

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

$$ P\mathbf{v} = \mathbf{v} - 2\mathbf{w}(\mathbf{w}^T\mathbf{v}) $$

性质：$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{v_1^2 + v_2^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$ 的公式

令

$$ u=v-\alpha e_1. $$ 然后令

$$ w=\frac{u}{\|u\|}. $$ 于是

$$ P=I-2ww^T. $$ 我们要证明：

$$ Pv=\alpha e_1. $$ 先计算

$$ u^Tv=(v-\alpha e_1)^Tv. $$ 展开：

$$ u^Tv=v^Tv-\alpha e_1^Tv. $$ 其中

$$ v^Tv=\|v\|^2=\alpha^2, $$ 因为 $|\alpha|=\|v\|$，而

$$ e_1^Tv=v_1. $$ 所以

$$ u^Tv=\alpha^2-\alpha v_1. $$ 再计算

$$ u^Tu=(v-\alpha e_1)^T(v-\alpha e_1). $$ 展开：

$$ u^Tu=v^Tv-2\alpha e_1^Tv+\alpha^2 e_1^Te_1. $$ 即

$$ u^Tu=\alpha^2-2\alpha v_1+\alpha^2. $$ 所以

$$ u^Tu=2(\alpha^2-\alpha v_1). $$ 因此

$$ \|u\|=\sqrt{2(\alpha^2-\alpha v_1)}. $$ 于是

$$ w=\frac{v-\alpha e_1}{\sqrt{2(\alpha^2-\alpha v_1)}}. $$ 可以验证 $\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$ 实对称矩阵：

$$ A=A^T. $$ 我们希望通过一系列正交相似变换：

$$ A^{(k+1)}=H_k^TA^{(k)}H_k, $$ 把它变成三对角矩阵。

因为 $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$

$$ T = \begin{bmatrix} \alpha_1 & \beta_1 & & & \\ \beta_1 & \alpha_2 & \beta_2 & & \\ & \beta_2 & \alpha_3 & \ddots & \\ & & \ddots & \ddots & \beta_{n-1} \\ & & & \beta_{n-1} & \alpha_n \end{bmatrix} $$

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

详细展开过程

设

$$ A= \begin{bmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n}\\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n}\\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} \end{bmatrix}. $$ 因为 $A$ 对称，所以第一列下面是：

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

$$ \begin{bmatrix} a_{11}\\ \beta_1\\ 0\\ \vdots\\ 0 \end{bmatrix}. $$ 取子向量：

$$ v= \begin{bmatrix} a_{21}\\a_{31}\\ \vdots\\a_{n1} \end{bmatrix}. $$ 构造一个 $(n-1)\times(n-1)$ 的 Householder 矩阵 $P_1$，使得：

$$ P_1v= \begin{bmatrix} \beta_1\\0\\ \vdots\\0 \end{bmatrix}. $$ 然后把它扩展成 $n\times n$ 矩阵：

$$ H_1= \begin{bmatrix} 1&0\\ 0&P_1 \end{bmatrix}. $$ 做相似变换： $$ A^{(2)}=H_1^TAH_1. $$ 由于 $H_1=H_1^T$，也可写成： $$ A^{(2)}=H_1AH_1. $$ 这会把第一列的第 3 到第 $n$ 个元素消掉。

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

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

4.14 第二步及以后

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

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

$$ v= \begin{bmatrix} a_{32}^{(2)}\\ a_{42}^{(2)}\\ \vdots\\ a_{n2}^{(2)} \end{bmatrix}. $$ 构造 $(n-2)\times(n-2)$ 的 Householder 矩阵 $P_2$，扩展为：

$$ H_2= \begin{bmatrix} I_2&0\\ 0&P_2 \end{bmatrix}. $$ 然后做：

$$ A^{(3)}=H_2^TA^{(2)}H_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分解 就是把一个矩阵 写成

$$ A = QR $$ Q：正交矩阵

R：上三角矩阵

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

5.1 核心迭代格式

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

基本迭代：

从

$$ A^{(0)}=A $$ 开始，每一步做 QR 分解：

$$ A^{(k)}=Q^{(k)}R^{(k)}. $$ 然后反过来相乘：

$$ A^{(k+1)}=R^{(k)}Q^{(k)}. $$ 看起来只是把乘法顺序反过来，但关键是它等价于一个正交相似变换。

因为

$$ A^{(k)}=Q^{(k)}R^{(k)}, $$ 左乘 $(Q^{(k)})^T$：

$$ (Q^{(k)})^TA^{(k)}=R^{(k)}. $$ 所以

$$ A^{(k+1)} = R^{(k)}Q^{(k)} = (Q^{(k)})^TA^{(k)}Q^{(k)}. $$ 因此

$$ A^{(k+1)} $$ 和

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

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

收敛性：

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

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

设矩阵可以对角化：

$$ A=X\Lambda X^{-1}, $$ 其中

$$ \Lambda=\operatorname{diag}(\lambda_1,\ldots,\lambda_n), $$ 并假设本征值大小满足：

$$ |\lambda_1|>|\lambda_2|>\cdots>|\lambda_n|. $$ 幂法中：

$$ A^k=X\Lambda^kX^{-1}. $$ 因为

$$ \Lambda^k=\operatorname{diag}(\lambda_1^k,\ldots,\lambda_n^k), $$ 所以高次幂会越来越强调最大本征值对应的方向。

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_{k+1} = c_{k+1},\quad \sin\theta_{k+1} = s_{k+1}$$

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

$$ \begin{bmatrix} c & s \\ -s & c \end{bmatrix} \begin{bmatrix} a_{k,k} \\ a_{k+1,k} \end{bmatrix} = \begin{bmatrix} r_{k,k} \\ 0 \end{bmatrix} $$

给定一个向量的两个分量： $$ \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}$）：

$$ A = \begin{bmatrix} 3 & 2 \\ 4 & 1 \end{bmatrix} $$

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

$$ r = \sqrt{3^2+4^2} = 5,\quad c = \frac{3}{5} = 0.6,\quad s = \frac{4}{5} = 0.8 $$

$$ P_1 = \begin{bmatrix} 0.6 & 0.8 \\ -0.8 & 0.6 \end{bmatrix},\quad R = P_1 A = \begin{bmatrix} 5 & 2 \\ 0 & 1.4 \end{bmatrix} $$

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

$$ Q = P_1^T P_2^T \cdots P_{n-1}^T $$

因为

$$ R=PA $$ 所以

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

5.3 Householder 做 QR 分解

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

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

$$ x=A_{1:n,1}. $$ 构造 Householder 矩阵 $P_1$，使得：

$$ P_1x= \begin{bmatrix} \alpha\\0\\ \vdots\\0 \end{bmatrix}. $$ 于是

$$ P_1A $$ 的第一列下面全部为零。

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

$$ \begin{bmatrix} 1&0\\ 0&P_2 \end{bmatrix}. $$ 继续做，最终得到：

$$ P_{n-1}\cdots P_2P_1A=R. $$ 因为每个 $P_i$ 都正交，所以它们乘积也正交。

令

$$ Q=P_1P_2\cdots P_{n-1}, $$ 则

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

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

$$ Q = P_1 P_2 \cdots P_{n-1},\quad R = P_{n-1} \cdots P_2 P_1 A,\quad A=QR $$ 💡 Householder 做 QR 分解比 Givens 旋转效率更高。

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

基础 QR 迭代：

$$ A_k=Q_kR_k, $$ 如果本征值很接近，收敛很慢。

于是引入位移 $\sigma_k$：

$$ A_k-\sigma_k I=Q_kR_k, $$ 然后

$$ A_{k+1}=R_kQ_k+\sigma_k I. $$ 由

$$ A_k-\sigma_k I=Q_kR_k $$ 得：

$$ R_k=Q_k^T(A_k-\sigma_k I). $$ 于是

$$ A_{k+1} = R_kQ_k+\sigma_k I = Q_k^T(A_k-\sigma_k I)Q_k+\sigma_k I. $$ 展开：

$$ A_{k+1} = Q_k^TA_kQ_k-\sigma_k Q_k^TIQ_k+\sigma_k I. $$ 因为

$$ Q_k^TIQ_k=I, $$ 所以

$$ A_{k+1} = Q_k^TA_kQ_k-\sigma_k I+\sigma_k I = Q_k^TA_kQ_k. $$ 因此带位移 QR 仍然是正交相似变换，本征值不变。

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

对

$$ A-\sigma I $$ 而言，本征值变成：

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

$$ |\lambda_j-\sigma| $$ 会很小。

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

一旦某个次对角元满足：

$$ |a_{i+1,i}|\approx 0, $$ 矩阵就近似分裂成：

$$ \begin{bmatrix} A_1&*\\ 0&A_2 \end{bmatrix}. $$ 如果矩阵是对称三对角，则更干净地分裂成两个独立块：

$$ \begin{bmatrix} T_1&0\\ 0&T_2 \end{bmatrix}. $$ 这叫 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$ 使得：

$$ A = Q U Q^{\dagger} $$ $U$ 的对角元就是 $A$ 的全部本征值。

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

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

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

证明思路用数学归纳法。

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

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

$$ Q_1=[q_1,q_2,\ldots,q_n] $$ 因为 $Q_1$ 酉：

$$ Q_1^\dagger Q_1=I $$ 考察：

$$ Q_1^\dagger A Q_1 $$ 它的第一列是：

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

2、带入特征方程

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

展开来写就是

$$ Q_1^{\dagger}q1 = \begin{bmatrix}q_1^{\dagger} \\q_1^{\dagger} \\q_2^{\dagger} \\\vdots \\q_n^{\dagger} \end{bmatrix} \cdot q_1 = \begin{bmatrix}q_1^{\dagger}q_1 \\q_1^{\dagger}q_1 \\q_2^{\dagger}q_1 \\\vdots \\q_n^{\dagger}q_1 \end{bmatrix} $$ 进一步，利用标准正交基的性质，得到上述式子其实就是 e1

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

所以矩阵形状为：

$$ Q_1^\dagger A Q_1 = \begin{bmatrix} \lambda_1 & *\\ 0 & A_2 \end{bmatrix}. $$ 然后对右下角 $A_2$ 继续做同样操作。归纳下去，最终得到上三角矩阵。

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

5.6 广义本征值问题

广义本征值问题是：

$$ Ax=\lambda Bx. $$ 假设：

$$ A=A^T,\qquad B=B^T,\qquad B>0. $$ 因为 $B$ 对称正定，可以做 Cholesky 分解：

$$ B=LL^T. $$ 原方程：

$$ Ax=\lambda LL^Tx. $$ 令

$$ y=L^Tx. $$ 则

$$ x=L^{-T}y. $$ 代回：

$$ A L^{-T}y=\lambda L y. $$ 左乘 $L^{-1}$：

$$ L^{-1}AL^{-T}y=\lambda y. $$ 定义

$$ C=L^{-1}AL^{-T}. $$ 则：

$$ Cy=\lambda y. $$ 而且 $C$ 是对称的：

$$ C^T=(L^{-1}AL^{-T})^T = L^{-1}A^TL^{-T} = L^{-1}AL^{-T} = C. $$ 因此广义对称本征值问题被化成了标准对称本征值问题。

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

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

$$ C = L^{-1} A (L^{-1})^T $$ $C$ 是对称矩阵，其本征值即原问题的本征值，但本征矢需要回代：

$$ \mathbf{y} = (L^{-1})^T \mathbf{x} $$

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

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

6.1 幂法 —— 求最大本征值

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

$$ \mathbf{x}^{(0)} = c_1\mathbf{v}_1 + c_2\mathbf{v}_2 + \cdots + c_n\mathbf{v}_n $$

反复左乘 $A$：

$$ \mathbf{x}^{(k)} = A^k \mathbf{x}^{(0)} = c_1\lambda_1^k \mathbf{v}_1 + c_2\lambda_2^k \mathbf{v}_2 + \cdots $$ 当 $k$ 很大时，$\lambda_1^k$ 主导（假设 $c_1 \neq 0$）：

$$ \mathbf{x}^{(k)} \approx c_1\lambda_1^k \mathbf{v}_1 $$ 说人话：不断被矩阵乘，向量方向会被拽到「最大本征值对应的本征矢方向」上。

6.2 算法步骤（带归一化）

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

$$ \|A^kx^{(0)}\| $$ 会指数增长，导致数值溢出。

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

所以每步都做归一化：

$$ y^{(k)}=Ax^{(k-1)}, \\ \space \\ x^{(k)} = \frac{y^{(k)}}{||y^{(k)}||} $$

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

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

$$ A = \begin{bmatrix} 2 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 2 \end{bmatrix},\quad \mathbf{x}^{(0)} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} $$

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

6.3 对称矩阵的幂法

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

$$ \rho(x)=\frac{x^TAx}{x^Tx}. $$ 如果 $x$ 已归一化，则

$$ \rho(x)=x^TAx. $$ 如果 $x=v_i$，则

$$ \rho(v_i)=\frac{v_i^TAv_i}{v_i^Tv_i}. $$ 因为

$$ Av_i=\lambda_i v_i, $$ 所以

$$ \rho(v_i)=\frac{v_i^T\lambda_i v_i}{v_i^Tv_i} = \lambda_i. $$ 所以当 $x^{(k)}$ 越接近本征向量时，Rayleigh 商越接近对应本征值。

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

$$ \lambda_{\max} \approx \frac{\mathbf{x}^T A \mathbf{x}}{\mathbf{x}^T\mathbf{x}} \quad\text{（Rayleigh 商）} $$

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

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

$$ q_1,q_2,\ldots,q_n. $$ 令当前归一化向量为：

$$ x=c_1q_1+c_2q_2+\cdots+c_nq_n. $$ 且

$$ \sum_i c_i^2=1. $$ Rayleigh 商：

$$ \rho(x)=x^TAx. $$ 代入：

$$ Ax=\sum_i c_i\lambda_i q_i. $$ 于是：

$$ x^TAx= \left(\sum_i c_iq_i\right)^T \left(\sum_j c_j\lambda_jq_j\right). $$ 由于正交归一：

$$ q_i^Tq_j=\delta_{ij}, $$ 所以：

$$ \rho(x)=\sum_i c_i^2\lambda_i. $$ 与 $\lambda_1$ 的误差：

$$ \rho(x)-\lambda_1 = \sum_i c_i^2\lambda_i-\lambda_1\sum_i c_i^2. $$ 因为

$$ \sum_i c_i^2=1, $$ 所以：

$$ \rho(x)-\lambda_1 = \sum_i c_i^2(\lambda_i-\lambda_1). $$ 其中 $i=1$ 项为 0，因此：

$$ \rho(x)-\lambda_1 = \sum_{i\geq2}c_i^2(\lambda_i-\lambda_1). $$ 注意误差由 $c_i^2$ 控制。

而幂法中：

$$ c_i^{(k)}\sim \left(\frac{\lambda_i}{\lambda_1}\right)^k. $$ 所以 Rayleigh 商误差大致是：

$$ O\left(\left|\frac{\lambda_2}{\lambda_1}\right|^{2k}\right). $$ 这就是对称矩阵下本征值估计更快的原因。

$$ \text{收敛速度} \propto \left(\frac{\lambda_2}{\lambda_1}\right)^{2k} \quad \text{（对称矩阵）} $$

vs

$$ \text{收敛速度} \propto \left(\frac{\lambda_2}{\lambda_1}\right)^{k} \quad \text{（一般矩阵）} $$

幂法失败的情形

1、如果说 c1 = 0 ， 初始向量没有最大本征方向分量！

2、如果最大本征值出现重，就是 1和2本征值相同，没有唯一的主导方向

3、如果 最大本征值是小于0的，那么向量可能是正负交替的（分奇偶）

4、如果 $\frac{\lambda_1}{\lambda_2} \sim 1$ 会出现收敛很慢的情况！

5、幂法只自然给出 最大模本征值，而不是 全部本征值！

6.4 反幂法（Inverse Power Method）—— 求指定本征值

如果想求接近某值 $q$ 的本征值：

1. 对 $(A - qI)^{-1}$ 做幂法
2. $(A - qI)^{-1}$ 的最大本征值 = $(\lambda_{\text{最接近}q} - q)^{-1}$

如果我们想找离某个数 $q$ 最近的本征值，考虑矩阵：

$$ (A-qI)^{-1}. $$ 如果

$$ Av_i=\lambda_i v_i, $$ 那么

$$ (A-qI)v_i=(\lambda_i-q)v_i. $$ 所以

$$ (A-qI)^{-1}v_i=\frac1{\lambda_i-q}v_i. $$ 也就是说，$(A-qI)^{-1}$ 的本征值是：

$$ \mu_i=\frac1{\lambda_i-q}. $$ 如果 $\lambda_i$ 离 $q$ 最近，则

$$ |\lambda_i-q| $$ 最小，于是

$$ \left|\frac1{\lambda_i-q}\right| $$ 最大。

所以对 $(A-qI)^{-1}$ 做幂法，会收敛到原矩阵中最靠近 $q$ 的本征值。

反幂法迭代：

$$ \mathbf{x}^{(k)} = (A - qI)^{-1} \mathbf{x}^{(k-1)} $$ 实际计算中解线性方程组

$$ (A - qI)\mathbf{y} = \mathbf{x}^{(k-1)} $$ 然后进行归一化，本征值仍利用 Rayleigh商进行估计！！！

如果每一步把 shift 取为当前 Rayleigh 商：

$$ q_k=\frac{(x^{(k)})^TAx^{(k)}}{(x^{(k)})^Tx^{(k)}}, $$ 然后解：

$$ (A-q_kI)y^{(k+1)}=x^{(k)}, $$ 再归一化：

$$ x^{(k+1)}=\frac{y^{(k+1)}}{\|y^{(k+1)}\|}, $$ 这就是 Rayleigh 商迭代。

对实对称矩阵，它在接近本征向量后通常有三次收敛速度，非常快。

💡 加速技巧：如果 $q$ 很接近本征值，收敛极快。实际应用中常结合位移 QR 法的思想。

6.5 带位移的幂法（Shifted Power Method）

本征矢不变，改变收敛目标：

$$ (A - \sigma I)\mathbf{x}^{(k-1)} = \mathbf{y}^{(k)} $$ 收敛到离 $\sigma$ 最近的本征值。

6.6 Wielandt 收缩法 —— 求次大本征值

假设已知最大本征对：

$$ Av_1=\lambda_1v_1. $$ 对于对称矩阵，且 $v_1$ 已归一化：

$$ v_1^Tv_1=1. $$ 定义：

$$ B=A-\lambda_1v_1v_1^T. $$ 看它对 $v_1$ 的作用：

$$ Bv_1=Av_1-\lambda_1v_1v_1^Tv_1. $$ 因为

$$ Av_1=\lambda_1v_1,\qquad v_1^Tv_1=1, $$ 所以

$$ Bv_1=\lambda_1v_1-\lambda_1v_1=0. $$ 也就是说，原来 $\lambda_1$ 被压成了 0。

对于其他正交本征向量 $v_i$，$i\geq2$：

$$ v_1^Tv_i=0. $$ 于是

$$ Bv_i=Av_i-\lambda_1v_1v_1^Tv_i = \lambda_iv_i-0 = \lambda_iv_i. $$ 所以 $B$ 的本征值为：

$$ 0,\lambda_2,\lambda_3,\ldots,\lambda_n. $$ 这样可以继续对 $B$ 做幂法，求次大本征值。

$B'$ 的本征值为 $\lambda_2, \lambda_3, \ldots, \lambda_n$，可以继续用幂法求 $\lambda_2$。

💡 类似量子计算中 VQD（Variational Quantum Deflation）算法求激发态的思路。

6.7 三种 Power 方法对比

⚠️ Power 方法的局限：只求一个本征值，而且当 $|\lambda_2/\lambda_1|$ 接近 1 时收敛极慢

7. 奇异值分解（SVD）

7.1 定义

对 $m \times n$ 矩阵 $A$，存在分解：

$$A = U S V^T$$

• $U$：$m \times m$ 正交矩阵，列是「左奇异向量」
• $V$：$n \times n$ 正交矩阵，列是「右奇异向量」
• $S$：$m \times n$ 对角矩阵，对角元 $\sigma_i \geq 0$ 是奇异值

矩阵形式（$m=3, n=2$ 的例子）：

$$ A = \underbrace{\begin{bmatrix} | & | & | \\ \mathbf{u}_1 & \mathbf{u}_2 & \mathbf{u}_3 \\ | & | & | \end{bmatrix}}_{U} \underbrace{\begin{bmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \\ 0 & 0 \end{bmatrix}}_{S} \underbrace{\begin{bmatrix} — & \mathbf{v}_1^T & — \\ — & \mathbf{v}_2^T & — \end{bmatrix}}_{V^T} $$

Proof

考察：

$$ A^TA. $$ 它是 $n\times n$ 对称矩阵：

$$ (A^TA)^T=A^TA. $$ 并且半正定，因为对任意 $x$：

$$ x^TA^TAx=(Ax)^T(Ax)=\|Ax\|^2\geq0. $$ 所以 $A^TA$ 可以正交对角化：

$$ A^TA=V\Lambda V^T. $$ 其中

$$ \Lambda=\operatorname{diag}(\lambda_1,\ldots,\lambda_n), $$ 且

$$ \lambda_i\geq0. $$ 定义奇异值：

$$ \sigma_i=\sqrt{\lambda_i}. $$ 于是：

$$ A^TAv_i=\sigma_i^2v_i. $$ 这里 $v_i$ 是 $V$ 的第 $i$ 列，称为右奇异向量。

构造左奇异向量

当 $\sigma_i>0$ 时，定义：

$$ u_i=\frac{1}{\sigma_i}Av_i. $$ 现在证明 $u_i$ 正交归一。

计算：

$$ u_i^Tu_j = \frac1{\sigma_i\sigma_j}(Av_i)^T(Av_j). $$ 即：

$$ u_i^Tu_j = \frac1{\sigma_i\sigma_j}v_i^TA^TAv_j. $$ 因为

$$ A^TAv_j=\sigma_j^2v_j, $$

所以：

$$ u_i^Tu_j = \frac{\sigma_j^2}{\sigma_i\sigma_j}v_i^Tv_j = \frac{\sigma_j}{\sigma_i}\delta_{ij}. $$ 当 $i=j$ 时：

$$ u_i^Tu_i=1. $$ 当 $i\neq j$ 时：

$$ u_i^Tu_j=0. $$ 因此这些 $u_i$ 是标准正交的。

得到 $Av_i=\sigma_i u_i$

由定义：

$$ u_i=\frac1{\sigma_i}Av_i, $$ 所以

$$ Av_i=\sigma_i u_i. $$ 把所有 $v_i$ 合起来：

$$ AV=U\Sigma. $$ 右乘 $V^T$：

$$ A=U\Sigma V^T. $$ 这就是 SVD。

如果某些 $\sigma_i=0$，则 $Av_i=0$，这些方向属于零空间。对应的 $u_i$ 可以通过补全正交基得到。

进一步的

同样可以证明：

$$ AA^Tu_i=\sigma_i^2u_i. $$ 因为：

$$ Av_i=\sigma_i u_i. $$ 左乘 $A^T$ 得：

$$ A^Tu_i=\sigma_i v_i. $$ 再左乘 $A$：

$$ AA^Tu_i=\sigma_i Av_i=\sigma_i^2u_i. $$ 所以：

• $A^TA$ 的本征向量是右奇异向量；
• $AA^T$ 的本征向量是左奇异向量；
• 二者的非零本征值相同，都是 $\sigma_i^2$。

7.2 SVD 与 $A^T A$ 的关系

$A^T A$ 是 $n \times n$ 对称半正定矩阵，其本征值为 $\sigma_i^2$，本征矢为 $\mathbf{v}_i$：

$$ A^T A \mathbf{v}_i = \sigma_i^2 \mathbf{v}_i $$ 计算 SVD 的步骤：

1. 求 A^T A 的本征值 σ_i^2 和本征矢 v_i
2. 奇异值 σ_i = sqrt(σ_i^2)
3. 构造 V = [v_1, v_2, ..., v_n]
4. U 的前 k 列：u_i = (1/σ_i) A v_i  (i=1,...,k, k=rank(A))
5. U 的剩余 m-k 列：补充线性无关向量，用 Gram-Schmidt 正交化

7.3 SVD与本征值

💡 如果 $A$ 是对称正定矩阵，则奇异值 = 本征值。一般情况下 $\sigma_i \neq |\lambda_i|$。

如果 $A$ 是对称正定矩阵，则： $$ A=Q\Lambda Q^T, $$ 且 $\lambda_i>0$。这时： $$ A^TA=A^2=Q\Lambda^2Q^T, $$ 所以： $$ \sigma_i=\lambda_i. $$ 如果 $A$ 是对称但不正定，则： $$ \sigma_i=|\lambda_i|. $$ 如果 $A$ 是一般矩阵，奇异值通常不等于本征值绝对值。

7.4 SVD 的图像压缩应用

低秩近似

SVD 可以写成外积和：

$$ A=\sum_{i=1}^r\sigma_i u_iv_i^T, $$ 其中 $r=\operatorname{rank}(A)$。

如果只保留前 $k$ 个最大奇异值：

$$ A_k=\sum_{i=1}^k\sigma_i u_iv_i^T. $$ 这就是最好的秩 $k$ 近似，称为 Eckart–Young 定理：

$$ \|A-A_k\|_F = \min_{\operatorname{rank}(B)\leq k}\|A-B\|_F. $$ 误差为：

$$ \|A-A_k\|_F^2 = \sum_{i=k+1}^r\sigma_i^2. $$ 这就是图像压缩的数学基础。

每个像素对应矩阵的一个元素。取前 $k$ 个最大奇异值（$k \ll N$）：

$$A \approx \sum_{i=1}^{k} \sigma_i \mathbf{u}_i \mathbf{v}_i^T$$

存储量对比：

8. Krylov 子空间方法

8.1 Krylov 子空间的定义

给定矩阵 $A$ 和初始向量 $\mathbf{q}_1$，$m$ 阶 Krylov 子空间为：

$$ \mathcal{K}_m(A, \mathbf{q}_1) = \text{span}\{\mathbf{q}_1, A\mathbf{q}_1, A^2\mathbf{q}_1, \ldots, A^{m-1}\mathbf{q}_1\} $$ 说人话：不断用矩阵 $A$ 去乘初始向量，每次乘出来的新向量构成一个子空间的基。这个子空间自动包含了对解最有用的方向。

8.2 为什么 Krylov 子空间有效？

幂法只看：

$$ q_1,Aq_1,A^2q_1,\ldots $$ 的最后一个方向。

Krylov 方法则把这些向量张成的整个空间都利用起来。

也就是说，它不是只取一个向量，而是在空间：

$$ \mathcal K_m(A,q_1) $$ 里面寻找最佳近似本征向量。

• 幂法本质上就是在 $\mathcal{K}_1$ 中找最优近似
• Krylov 子空间方法用 $\mathcal{K}_m$（$m$ 维），信息量大得多
• $m$ 维 Krylov 子空间能捕获前 $m$ 个主导本征矢的信息

8.3 Krylov 基的性质

8.3.1 多项式滤波视角

Krylov 子空间中的任意向量都可以写成：

$$ y=p_{m-1}(A)q_1, $$ 其中 $p_{m-1}$ 是次数不超过 $m-1$ 的多项式。

如果因为 $\mathbf{v}_i$ 是本征矢，满足 $A\mathbf{v}_i = \lambda_i \mathbf{v}_i$，所以 $A^2\mathbf{v}_i = \lambda_i^2 \mathbf{v}_i$。同理，多项式作用上去就会变成：

$$ q_1=\sum_i c_iv_i, $$ 则

$$ p(A)q_1=\sum_i c_ip(\lambda_i)v_i. $$ 你想得到系统最低能量（基态本征值 $\lambda_1$）对应的本征矢 $\mathbf{v}_1$。 你可以设计（或者让算法自动寻找）一个多项式 $p(x)$。这个多项式在 $x = \lambda_1$ 处的值非常大，而在其他本征值处的值都接近 $0$（就像一个滤波器）。 这样，在求和项 $\sum c_i p(\lambda_i) \mathbf{v}_i$ 中，非目标项都被压制到了几乎为 0，只有 $c_1 p(\lambda_1)\mathbf{v}_1$ 这一项脱颖而出！

所以 Krylov 方法本质上是在构造多项式 $p$，使得：

• 目标本征值对应的 $p(\lambda_i)$ 大；
• 非目标本征值对应的 $p(\lambda_i)$ 小。

这就是“自动筛选重要方向”的原因。

8.3.2 Krylov 投影思想

直接解超级大矩阵 $A\mathbf{x} = \lambda \mathbf{x}$ 是不可能的。 既然如此，我们能不能只在 $m$ 维的 Krylov 子空间里找近似解？（通常 $m$ 很小，比如 $m=30$ 或 $50$）

我们把这个小空间里的一组正交基按列排好，拼成一个大矩阵 $Q_m = [\mathbf{q}_1, \mathbf{q}_2, \dots, \mathbf{q}_m]$（尺寸是 $n \times m$）。 因为基向量是正交归一的，所以满足 $Q_m^T Q_m = I_m$（$m \times m$ 的单位阵）

我们不直接在 $n$ 维大空间里解：

$$ Ax=\lambda x. $$ 而是在 $m$ 维子空间中近似（我们假设真实的本征矢 $\mathbf{x}$ 可以由这个小空间的基线性组合得到），这里y是一个只有m维的小向量：

$$ x\approx Q_my, $$ 其中 $Q_m$ 的列是 Krylov 子空间的一组正交基。

要求残差（近似解和精确解的差距）：

$$ r=AQ_my-\theta Q_my $$ 利用伽辽金条件（Golerkin Condition）：我们不能让残差等于 0（因为空间不够大，不可能完全精确），但我们可以要求残差 $\mathbf{r}$ 与我们这个小空间 $Q_m$ 垂直。 意思是，误差不能在我们关心的空间里露头。即要求与子空间正交：

$$ Q_m^Tr=0. $$ 代入：

$$ Q_m^T(AQ_my-\theta Q_my)=0. $$ 得到：

$$ Q_m^TAQ_my-\theta Q_m^TQ_my=0. $$ 因为

$$ Q_m^TQ_m=I, $$ 所以：

$$ (Q_m^TAQ_m)y=\theta y. $$ 定义小矩阵 m ✖️ m的小矩阵：

$$ H_m=Q_m^TAQ_m. $$ 于是只需解小规模本征值问题：

$$ H_my=\theta y. $$ $\theta$ 称为 Ritz 值，$Q_my$ 称为 Ritz 向量。

我们把一个 $10^{10} \times 10^{10}$ 的巨无霸本征值问题，通过投影，变成了一个 $50 \times 50$ 的微型本征值问题。 我们用前几章介绍的常规方法（如 QR 算法）轻松对角化 $H_m$，求出本征值 $\theta$（称为 Ritz 值）和本征矢 $\mathbf{y}$。 然后通过 $\mathbf{x} \approx Q_m \mathbf{y}$，就能非常精准地还原出原大矩阵的本征矢（称为 Ritz 向量）。

在 krylov子空间中， A的投影矩阵是 上 Hessenberg形式 （即除了主对角线、紧挨着的下一条次对角线以外，下三角其余元素全是0）

Krylov基生成过程

这个基的生成过程（也就是著名的 Arnoldi 迭代算法）一步一步彻底拆解。

首先，我们要明白为什么不直接用 $\{ \mathbf{b}, A\mathbf{b}, A^2\mathbf{b}, \ldots \}$ 做基？ 因为当 $k$ 稍微变大一点，根据幂法的原理，$A^k\mathbf{b}$ 会迅速向最大本征值对应的本征矢靠拢。这就导致这几个向量几乎指向同一个方向（也就是几乎共线），数学上叫“高度线性相关”。用它们做基，矩阵极其病态，计算机稍微一算就全是舍入误差。

所以，我们需要一边用 $A$ 乘，一边做正交化。这个完美的结合，就是 Arnoldi 过程。

第一步：准备启动（Step 0）

我们有一个初始向量 $\mathbf{b}$。

1. 计算它的模长：$\beta = \|\mathbf{b}\|$。
2. 把它归一化，得到我们第一根正交基向量 $\mathbf{q}_1$：

$$\mathbf{q}_1 = \frac{\mathbf{b}}{\beta}$$ 此时，我们的基里面只有一根向量：$\{\mathbf{q}_1\}$。

第二步：生成第二根基向量 $\mathbf{q}_2$（Step 1）

现在我们想通过 $A$ 作用到已有的基上来扩展空间。

1. 用 $A$ 乘已有的最新基向量：

$$\mathbf{v} = A\mathbf{q}_1$$

1. 正交化（剔除在 $\mathbf{q}_1$ 方向上的投影）：

计算 $\mathbf{v}$ 与 $\mathbf{q}_1$ 的内积，记为 $h_{11}$： $$h_{11} = \mathbf{q}_1^T \mathbf{v}$$ 然后从 $\mathbf{v}$ 中减去 $\mathbf{q}_1$ 方向的分量，得到一个纯净的、与 $\mathbf{q}_1$ 垂直的新方向向量 $\mathbf{r}$： $$\mathbf{r} = \mathbf{v} - h_{11}\mathbf{q}_1$$

1. 归一化：

计算这根新向量的模长，记为 $h_{21}$： $$h_{21} = \|\mathbf{r}\|$$ 如果 $h_{21} \neq 0$，归一化得到第二根基向量 $\mathbf{q}_2$： $$\mathbf{q}_2 = \frac{\mathbf{r}}{h_{21}}$$ 此时，我们的基向量组变成了：$\{\mathbf{q}_1, \mathbf{q}_2\}$。它们彼此正交且长度为 1。

第三步：生成第三根基向量 $\mathbf{q}_3$（Step 2）

重复这个思路：

1. 用 $A$ 乘上一步新生成的基向量 $\mathbf{q}_2$：

$$\mathbf{v} = A\mathbf{q}_2$$

1. 正交化（剔除在已有的 $\mathbf{q}_1$ 和 $\mathbf{q}_2$ 上的所有投影）：

• 剔除 $\mathbf{q}_1$ 方向：$h_{12} = \mathbf{q}_1^T \mathbf{v}$
• 剔除 $\mathbf{q}_2$ 方向：$h_{22} = \mathbf{q}_2^T \mathbf{v}$

从 $\mathbf{v}$ 中把这两份投影全部扣掉： $$\mathbf{r} = \mathbf{v} - h_{12}\mathbf{q}_1 - h_{22}\mathbf{q}_2$$

1. 归一化：

计算模长：$h_{32} = \|\mathbf{r}\|$，得到第三根基向量： $$\mathbf{q}_3 = \frac{\mathbf{r}}{h_{32}}$$ 此时，基向量组变成了：$\{\mathbf{q}_1, \mathbf{q}_2, \mathbf{q}_3\}$。

归纳：第 $j$ 步的通用公式（Step $j$）

已知前 $j$ 个已经相互正交归一的基 $\{\mathbf{q}_1, \mathbf{q}_2, \ldots, \mathbf{q}_j\}$。我们要生出第 $j+1$ 个：

1. 作用算符：

$$\mathbf{v} = A\mathbf{q}_j$$

1. 扣除前 $j$ 个方向的投影（施密特正交化）：

对于每一个 $i = 1, 2, \ldots, j$： $$h_{ij} = \mathbf{q}_i^T \mathbf{v}$$ 扣除投影后的残差向量： $$\mathbf{r} = \mathbf{v} - \sum_{i=1}^{j} h_{ij}\mathbf{q}_i$$

1. 归一化：

$$h_{j+1, j} = \|\mathbf{r}\|$$ $$\mathbf{q}_{j+1} = \frac{\mathbf{r}}{h_{j+1, j}}$$

这就是完整的 Krylov 基生成过程（Arnoldi 迭代）。

💡 为什么这样会自然产生“上 Hessenberg 矩阵”？

我们把第 $j$ 步归一化公式 $\mathbf{q}_{j+1} = \frac{\mathbf{r}}{h_{j+1, j}}$ 变形一下（乘以分母）： $$\mathbf{r} = h_{j+1, j} \mathbf{q}_{j+1}$$

再把第二步正交化公式代入这个式子： $$\mathbf{v} - \sum_{i=1}^{j} h_{ij}\mathbf{q}_i = h_{j+1, j} \mathbf{q}_{j+1}$$

因为 $\mathbf{v} = A\mathbf{q}_j$，我们移项整理得到： $$A\mathbf{q}_j = h_{1j}\mathbf{q}_1 + h_{2j}\mathbf{q}_2 + \ldots + h_{jj}\mathbf{q}_j + h_{j+1, j}\mathbf{q}_{j+1}$$

这个式子的物理意义非常震撼： 当你用大矩阵 $A$ 作用到第 $j$ 个基向量上时，产生的新向量最多只能投影到第 $j+1$ 个基向量上！它跟再后面的基向量（例如 $\mathbf{q}_{j+2}, \mathbf{q}_{j+3}$）绝对是垂直的（内积为0）。

这就意味着，如果把所有系数 $h_{ij}$ 排成一个矩阵： $$ H = \begin{bmatrix} h_{11} & h_{12} & h_{13} & \cdots \\ h_{21} & h_{22} & h_{23} & \cdots \\ 0 & h_{32} & h_{33} & \cdots \\ 0 & 0 & h_{43} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} $$ 由于对于所有 $i > j+1$ 的位置，其系数 $h_{ij} = 0$，所以主对角线下方的第二条线开始，全部都是 0！ 这正是上 Hessenberg 矩阵！

所以，Arnoldi 算法妙就妙在：它在用施密特正交化建立空间基底的同时，顺手就把在这个空间里投影的大矩阵 $A$，记录成了一个极其简单、几乎是对角阵的 Hessenberg 小矩阵 $H$！接下来对这个 $H$ 做对角化，就能轻松拿到我们想要的本征值。

直观解释：

回忆一下 Krylov 基的生成过程：

• $\mathbf{q}_2$ 是由 $A\mathbf{q}_1$ 正交化得到的；
• $\mathbf{q}_3$ 是由 $A\mathbf{q}_2$ 正交化得到的；
• 一般地， $A \mathbf{q}_j$ 只能跟前面的 $\mathbf{q}_1, \dots, \mathbf{q}_j, \mathbf{q}_{j+1}$ 发生关联（投影不为 0）。 对于比 $j+1$ 更靠后的基向量（如 $\mathbf{q}_{j+2}$ 甚至更远），由于生成顺序的关系，它们和 $A \mathbf{q}_j$ 必然是正交的：

$$ \mathbf{q}_i^T A \mathbf{q}_j = 0 \quad (\text{当 } i > j+1) $$

而 $H_{ij} = \mathbf{q}_i^T A \mathbf{q}_j$ 正是矩阵 $H_m$ 的矩阵元！ 这说明，只要行数 $i$ 大于列数 $j+1$，对应的矩阵元就全部为 0。

写出来就是： $$ H_m = \begin{bmatrix} h_{11} & h_{12} & h_{13} & \cdots \\ h_{21} & h_{22} & h_{23} & \cdots \\ 0 & h_{32} & h_{33} & \cdots \\ 0 & 0 & h_{43} & \cdots \\ \end{bmatrix} $$ 这正是上 Hessenberg 矩阵！这种矩阵不仅容易存储，而且极易使用 QR 算法进行对角化。Krylov 方法之所以如此高效，秘密全在于此！

💡 这正是 Krylov 子空间方法的核心：在 Krylov 子空间中，$A$ 的投影矩阵是上 Hessenberg 形式——比原矩阵更容易对角化。

9. Lanczos 迭代

9.1 核心思想

对对称矩阵 $A$，在 Krylov 子空间中构造一组正交基（Lanczos 基），使得 $A$ 在这组基下的投影是三对角矩阵。

9.2 Lanczos 三对角化

初始化：

$$ q_0=0,\qquad \beta_0=0. $$ 给定归一化 $q_1$。

第 $i$ 步：

$$ u=Aq_i-\beta_{i-1}q_{i-1}. $$ 计算对角元：

$$ \alpha_i=q_i^Tu. $$ 去掉当前方向：

$$ u=u-\alpha_iq_i. $$ 计算次对角元：

$$ \beta_i=\|u\|. $$ 如果

$$ \beta_i=0, $$ 说明 Krylov 子空间已经不再扩展，算法终止。

否则：

$$ q_{i+1}=\frac{u}{\beta_i}. $$ 于是得到三项关系：

$$ Aq_i=\beta_{i-1}q_{i-1}+\alpha_iq_i+\beta_iq_{i+1}. $$ 矩阵形式

把列向量合成： $$ Q_m=[q_1,q_2,\ldots,q_m]. $$ 则： $$ AQ_m=Q_mT_m+\beta_mq_{m+1}e_m^T. $$ 其中： $$ T_m= \begin{bmatrix} \alpha_1&\beta_1\\ \beta_1&\alpha_2&\beta_2\\ &\beta_2&\alpha_3&\ddots\\ &&\ddots&\ddots&\beta_{m-1}\\ &&&\beta_{m-1}&\alpha_m \end{bmatrix}. $$ 如果忽略最后残差项，就是： $$ Q_m^TAQ_m=T_m. $$ 所以 Lanczos 把大矩阵 $A$ 投影成小三对角矩阵 $T_m$。

给定对称矩阵 $A$ 和初始归一化向量 $\mathbf{q}_1$（$\mathbf{q}_0 = \mathbf{0}$, $\beta_0 = 0$）：

Lanczos 迭代：

for i = 1, 2, ..., m:
    u = A q_i - β_{i-1} q_{i-1}       ← 矩阵乘 + 减去上一步贡献
    α_i = q_i^T u                      ← 对角元
    u = u - α_i q_i                    ← 再减当前方向
    β_i = ||u||                        ← 次对角元
    if β_i ≈ 0: break
    q_{i+1} = u / β_i                  ← 归一化得到下一个基向量
end

每次迭代只需要3 个向量（$\mathbf{q}_{i-1}, \mathbf{q}_i, \mathbf{q}_{i+1}$），内存效率极高！

9.3 得到的结构

$m$ 步后得到三对角矩阵 $T_m$：

$$T_m = \begin{bmatrix} \alpha_1 & \beta_1 & & \\ \beta_1 & \alpha_2 & \beta_2 & \\ & \beta_2 & \ddots & \ddots \\ & & \ddots & \ddots & \beta_{m-1} \\ & & & \beta_{m-1} & \alpha_m \end{bmatrix}$$

且 $A Q_m = Q_m T_m + \beta_m \mathbf{q}_{m+1} \mathbf{e}_m^T$

9.4 数值例子

$$ A = \begin{bmatrix} 3 & 2 & 5 \\ 2 & 3 & 5 \\ 5 & 5 & 0 \end{bmatrix},\quad \mathbf{q}_1 = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} $$

精确本征值：$10,\ -5,\ 1$

$m=2$ 步 Lanczos：

$$ \mathbf{q}_2 = \frac{A\mathbf{q}_1 - \alpha_1\mathbf{q}_1}{\beta_1} = \text{计算得} \begin{bmatrix} 0 \\ 0.371 \\ 0.928 \end{bmatrix} $$

得到 $2\times2$ 三对角矩阵：

$$T_2 = \begin{bmatrix} 3 & \sqrt{29} \\ \sqrt{29} & 29/29 \approx 1 \end{bmatrix}$$

对角化 $T_2$ 得到近似本征值 $8.83$ 和 $-1.97$。仅 2 步就找到了"极端本征值"（最大和最小）的近似！

完整 $m=3$（$m=n$，全矩阵） 则严格恢复本征值 $[10, 1, -5]$。

9.5 Lanczos 的收敛特性

设

$$ T_my=\theta y. $$ 则近似本征向量为：

$$ x=Q_my. $$ 残差：

$$ r=Ax-\theta x. $$ 代入：

$$ r=AQ_my-\theta Q_my. $$ 利用 Lanczos 关系：

$$ AQ_m=Q_mT_m+\beta_mq_{m+1}e_m^T. $$ 所以：

$$ r=(Q_mT_m+\beta_mq_{m+1}e_m^T)y-\theta Q_my. $$ 因为

$$ T_my=\theta y, $$ 所以前两项抵消：

$$ Q_mT_my-\theta Q_my=0. $$ 剩下：

$$ r=\beta_mq_{m+1}e_m^Ty. $$ 因此残差范数：

$$ \|r\|=|\beta_m|\cdot |e_m^Ty|. $$ 这提供了一个非常方便的收敛判据。 $$ \text{极端本征值（最大/最小）收敛快} \quad > \quad \text{中间本征值收敛慢} $$

数值验证（$A = B^T B$，$B \in \mathbb{R}^{10\times10}$ 正定对称，$m$ 取 2~10）：

💡 Lanczos = 在 Krylov 子空间中对对称矩阵做截断对角化。截断类似于 SVD 图像压缩——保留最重要的信息、丢弃细节。极端本征值对应"最重要成分"。

9.6 Lanczos 的优势

理论上 Lanczos 向量彼此正交：

$$ q_i^Tq_j=0,\qquad i\neq j. $$ 但浮点运算中，由于舍入误差，正交性会逐渐丢失。

后果是：

• 已经收敛的本征值可能重复出现；
• 出现 ghost eigenvalues；
• Ritz 值稳定性变差。

解决办法包括：

• full reorthogonalization；
• selective reorthogonalization；
• partial reorthogonalization；
• restarted Lanczos。

10. Arnoldi 迭代

10.1 动机

Lanczos 要求矩阵对称（或 Hermite）。对于不对称矩阵，需要 Arnoldi 迭代。

10.2 Arnoldi 算法

给定不对称矩阵 $A$ 和初始归一化向量 $\mathbf{q}_1$：

for k = 1, 2, ..., m:
    v = A q_k
    for j = 1 to k:                    ← 与前面所有 q_j 正交化
        h_{j,k} = q_j^T v              ← 内积（投影系数）
        v = v - h_{j,k} q_j
    end
    h_{k+1,k} = ||v||                   ← 次对角元
    if h_{k+1,k} ≈ 0: break
    q_{k+1} = v / h_{k+1,k}
end

关键区别：Lanczos 只需与 $\mathbf{q}_{k-1}$ 和 $\mathbf{q}_k$ 正交化（利用对称性），Arnoldi 需要与所有已生成的 $\mathbf{q}_1, \ldots, \mathbf{q}_k$ 正交化。

10.3 结果：上 Hessenberg 矩阵

$m$ 步 Arnoldi 得到 $m \times m$ 上 Hessenberg 矩阵 $H_m$：

$$H_m = \begin{bmatrix} h_{1,1} & h_{1,2} & h_{1,3} & \cdots & h_{1,m} \\ h_{2,1} & h_{2,2} & h_{2,3} & \cdots & h_{2,m} \\ 0 & h_{3,2} & h_{3,3} & \cdots & h_{3,m} \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & h_{m,m-1} & h_{m,m} \end{bmatrix}$$

（如果 $A$ 对称，则 $H_m$ 退化为三对角 = 回到 Lanczos）

10.4 Arnoldi vs Lanczos 对比

💡 ARPACK 是 Arnoldi Package 的缩写，是求解大型稀疏矩阵本征值问题的主流软件包。当 $\mathbf{q}_k$ 接近零向量时算法自然终止（精确不变子空间已找到）。

11. Schrödinger 方程的对角化求解

11.1 问题提法

径向 Schrödinger 方程是一个本征值问题：

$$H|\psi\rangle = E|\psi\rangle$$

其中 $H$ 是实对称矩阵，$\psi$ 是波函数，$E$ 是能级。

11.2 求解策略流程

                    ┌─────────────────────────────┐
                    │   哈密顿量矩阵 H (实对称)     │
                    └────────────┬────────────────┘
                                 │
                ┌────────────────┼────────────────┐
                ▼                ▼                ▼
        ┌───────────────┐ ┌───────────┐  ┌─────────────────┐
        │ 直接法 (稠密)   │ │ 幂法       │  │ Lanczos (稀疏)   │
        │ Householder   │ │ (最大E)    │  │ (极端本征值)     │
        │ → QR对角化     │ │ 反幂法     ｜  │ 3向量内存        │
        │ 全部本征值     │ │ (指定E)   │   │ 10^9~10^10基     │
        └───────────────┘ └───────────┘ └─────────────────┘

11.3 具体方案对比

12. DMRG：密度矩阵重整化群

12.1 问题背景

严格对角化计算量随系统尺寸指数增长（$2^L$），无法处理大系统。

12.2 NRG（数值重整化群, Wilson 1974）

K. Wilson 因重整化群在临界现象和 Kondo 效应中的应用获 1982 年诺贝尔物理学奖。

NRG 步骤：

① 初始块 A (长度 L)，Hilbert 空间维数 M
     ↓
② 形成复合块 AA (长度 2L)，哈密顿量 H_AA (维数 M²)
     ↓
③ 用 Lanczos 对角化 H_AA，取 M 个最低本征态
     ↓
④ 将 H_AA 投影到截断空间 → H_AA-tr
     ↓
⑤ 令 2L→L, AA→A, H_AA-tr→H_A，回到步骤②
     ↓  重复直到达到目标系统尺寸

12.3 NRG 的局限 → DMRG 的革新

💡 DMRG 的关键创新：不按能量截断，而按约化密度矩阵的奇异值（SVD）截断——本质是在量子纠缠熵中做截断。

12.4 DMRG 算法流程

Warm-up 阶段：

1. Lanczos 求基态
2. 构造 A' 和 B' 块
3. 构造约化密度矩阵
4. 对角化约化密度矩阵
5. 将 H_A 变换到密度矩阵本征基，截断
6. 重复

Sweeping 阶段： 通过多轮左右扫描优化波函数，直到收敛。

12.5 DMRG 的应用范围

自旋链、量子点、原子核、量子化学、量子信息……DMRG 已成为处理一维/准一维多体系统的黄金标准方法。

13. 算法对比与总结

13.1 本征值问题与量子多体问题

量子多体系统 → 组态混合 → 大型矩阵对角化 → 本征值问题的数值方法

13.2 所有算法一览

13.3 选择算法的决策树

需要求全部本征值？
├── 是 → 矩阵是否稠密？
│        ├── 是 → Householder + QR 算法
│        └── 否 → 能否做三对角化？
│                 └── 是 → Householder + QR（三对角矩阵）
│
└── 否 → 需要几个本征值？
         ├── 1 个（最大）→ 幂法
         ├── 1 个（指定位置）→ 反幂法
         └── 多个极端值 → 矩阵是否对称？
              ├── 是 → Lanczos
              └── 否 → Arnoldi

13.4 算法背后的数学思想

13.5 二十世纪十大算法（与本课程相关的两个）

1. QR 算法 —— 计算矩阵全部本征值的标准方法
2. Krylov 子空间方法 —— 大型稀疏矩阵本征值问题的革命性方法

13.6 关键公式速查

Akuiro 整理补充· 2026-07-08