本节属“补充篇 · 数值线性代数工具箱” :2.x 是补充篇内部编号,不是课堂第 2 章的节号。课堂第 2 章是线性方程组的解法 。
前四节把矩阵当作「算子」「方程组」「投影工具」,这一节换视角:矩阵的谱(特征值集合)刻画它的本质行为 ——条件数回答「解被放大多少」(2.4),谱半径 ρ ( A ) \rho(A) ρ ( A ) 回答「迭代收不收敛、多快」,SVD(2.6)与 QR 迭代也都建立在谱上。本节重点是特征值的扰动敏感性 :1.2.6 中 conditioning 与 stability 的区分在此重演,且结论反直觉——对称矩阵温和,非正规矩阵「开方敏感」。路线:谱半径及其上界 → Gerschgorin 圆盘(不算特征值也能定位它们)→ 对称情形(Rayleigh 商与 Courant–Fischer)→ 非正规情形的定量警告。
直觉 · 特征向量 = 「只被拉伸不被转向」的方向
矩阵乘在一个向量上,一般既会把它拉长,也会把它转向。但有些特殊方向上,矩阵只拉长、不转向 ——这些方向就是特征向量,拉长的倍数就是特征值:
A v = λ v , v ≠ 0 . A\boldsymbol{v}=\lambda\boldsymbol{v},\qquad \boldsymbol{v}\neq\boldsymbol{0}. A v = λ v , v = 0 .
λ \lambda λ 可以是负数(方向反过来)、复数(带一点旋转)、零(把向量压扁到零)。全部特征值合起来就是矩阵的「指纹」 :行列式 = ∏ λ i =\prod\lambda_i = ∏ λ i (体积缩放)、迹 = ∑ λ i =\sum\lambda_i = ∑ λ i 、可逆性(是否 λ i ≠ 0 \lambda_i\neq0 λ i = 0 )、迭代收敛(是否 ∣ λ i ∣ < 1 \lvert\lambda_i\rvert\lt1 ∣ λ i ∣ < 1 )、稳定性(看特征值的位置),统统由这份指纹决定。而指纹里最常用的一个数,就是谱半径 ρ ( A ) = max i ∣ λ i ∣ \rho(A)=\max_i\lvert\lambda_i\rvert ρ ( A ) = max i ∣ λ i ∣ 。
定义 2.5.1 · 特征值、谱与谱半径
设 A ∈ R n × n A\in\mathbb{R}^{n\times n} A ∈ R n × n 。若存在 λ ∈ C \lambda\in\mathbb{C} λ ∈ C 与非零向量 v ∈ C n \boldsymbol{v}\in\mathbb{C}^n v ∈ C n 使
A v = λ v (2.25) A\boldsymbol{v}=\lambda\boldsymbol{v} \tag{2.25} A v = λ v ( 2.25 )
则称 λ \lambda λ 为 A A A 的特征值 (eigenvalue),v \boldsymbol{v} v 为对应的特征向量 (eigenvector)。A A A 的全部特征值(含代数重数)构成 A A A 的谱 (spectrum),记 σ ( A ) \sigma(A) σ ( A ) (注意与奇异值 σ i \sigma_i σ i 区分)。谱半径 (spectral radius)定义为
ρ ( A ) = max λ ∈ σ ( A ) ∣ λ ∣ . (2.26) \rho(A)=\max_{\lambda\in\sigma(A)}\lvert\lambda\rvert . \tag{2.26} ρ ( A ) = λ ∈ σ ( A ) max ∣ λ ∣ . ( 2.26 )
λ \lambda λ 是特征值 ⟺ \iff ⟺ det ( A − λ I ) = 0 \det(A-\lambda I)=0 det ( A − λ I ) = 0 (特征多项式 为 0)。
定理 2.5.1(新增) · 基本性质
设 A ∈ R n × n A\in\mathbb{R}^{n\times n} A ∈ R n × n ,λ 1 , … , λ n \lambda_1,\dots,\lambda_n λ 1 , … , λ n (含重数)为它的特征值,则
tr A = ∑ i = 1 n λ i \operatorname{tr}A=\sum_{i=1}^n\lambda_i tr A = ∑ i = 1 n λ i ,det A = ∏ i = 1 n λ i \det A=\prod_{i=1}^n\lambda_i det A = ∏ i = 1 n λ i ;
σ ( A T ) = σ ( A ) \sigma(A^{\mathrm T})=\sigma(A) σ ( A T ) = σ ( A ) (转置不改变特征值;但特征向量一般不同);
若 A A A 可逆,则 σ ( A − 1 ) = { 1 / λ : λ ∈ σ ( A ) } \sigma(A^{-1})=\{1/\lambda:\lambda\in\sigma(A)\} σ ( A − 1 ) = { 1/ λ : λ ∈ σ ( A )} ;
σ ( A m ) = { λ m : λ ∈ σ ( A ) } \sigma(A^m)=\{\lambda^m:\lambda\in\sigma(A)\} σ ( A m ) = { λ m : λ ∈ σ ( A )} ;
若 A A A 对称,则全部特征值为实数 ,且存在正交 矩阵 V V V 使 V T A V = diag ( λ 1 , … , λ n ) V^{\mathrm T}AV=\operatorname{diag}(\lambda_1,\dots,\lambda_n) V T A V = diag ( λ 1 , … , λ n ) ;
若 A A A 相似于 B B B (B = P − 1 A P B=P^{-1}AP B = P − 1 A P ),则 σ ( B ) = σ ( A ) \sigma(B)=\sigma(A) σ ( B ) = σ ( A ) 。
证明
1: A A A 的特征多项式 p ( λ ) = det ( A − λ I ) = ( − 1 ) n ( λ − λ 1 ) ⋯ ( λ − λ n ) p(\lambda)=\det(A-\lambda I)=(-1)^n(\lambda-\lambda_1)\cdots(\lambda-\lambda_n) p ( λ ) = det ( A − λ I ) = ( − 1 ) n ( λ − λ 1 ) ⋯ ( λ − λ n ) 。比较两边 λ n − 1 \lambda^{n-1} λ n − 1 的系数得 tr A = ∑ λ i \operatorname{tr}A=\sum\lambda_i tr A = ∑ λ i (因为 det ( A − λ I ) \det(A-\lambda I) det ( A − λ I ) 展开中 λ n − 1 \lambda^{n-1} λ n − 1 项来自「取对角元 a i i − λ a_{ii}-\lambda a ii − λ 、其余取常数项」),取 λ = 0 \lambda=0 λ = 0 得 det A = ∏ λ i \det A=\prod\lambda_i det A = ∏ λ i 。
2、3、4: A v = λ v ⇒ A T A\boldsymbol{v}=\lambda\boldsymbol{v}\Rightarrow A^{\mathrm T} A v = λ v ⇒ A T 的(左)特征向量满足 v T A = λ v T \boldsymbol{v}^{\mathrm T}A=\lambda\boldsymbol{v}^{\mathrm T} v T A = λ v T ,故 λ \lambda λ 也是 A T A^{\mathrm T} A T 的特征值;A v = λ v ⇒ A − 1 v = λ − 1 v A\boldsymbol{v}=\lambda\boldsymbol{v}\Rightarrow A^{-1}\boldsymbol{v}=\lambda^{-1}\boldsymbol{v} A v = λ v ⇒ A − 1 v = λ − 1 v (A A A 可逆时 λ ≠ 0 \lambda\neq0 λ = 0 );归纳得 A m v = λ m v A^m\boldsymbol{v}=\lambda^m\boldsymbol{v} A m v = λ m v 。
5: 见 1.3.4 已用到的谱定理;完整证明见习题 B2 (要点:实对称矩阵的任一特征值必为实数——在 A v = λ v A\boldsymbol{v}=\lambda\boldsymbol{v} A v = λ v 两边左乘 v ‾ T \overline{\boldsymbol{v}}^{\mathrm T} v T 得 λ ∥ v ∥ 2 = v ‾ T A v \lambda\lVert\boldsymbol{v}\rVert^2=\overline{\boldsymbol{v}}^{\mathrm T}A\boldsymbol{v} λ ∥ v ∥ 2 = v T A v ,其共轭转置给出 λ ‾ ∥ v ∥ 2 = v T A T v ‾ \overline{\lambda}\lVert\boldsymbol{v}\rVert^2=\boldsymbol{v}^{\mathrm T}A^{\mathrm T}\overline{\boldsymbol{v}} λ ∥ v ∥ 2 = v T A T v ,由 A = A T A=A^{\mathrm T} A = A T 得 λ = λ ‾ \lambda=\overline{\lambda} λ = λ )。
6: det ( B − λ I ) = det ( P − 1 A P − λ I ) = det ( P − 1 ( A − λ I ) P ) = det ( A − λ I ) \det(B-\lambda I)=\det(P^{-1}AP-\lambda I)=\det(P^{-1}(A-\lambda I)P)=\det(A-\lambda I) det ( B − λ I ) = det ( P − 1 A P − λ I ) = det ( P − 1 ( A − λ I ) P ) = det ( A − λ I ) 。
∎
三者关系:ρ ( A ) ≤ ∥ A ∥ \rho(A)\le\lVert A\rVert ρ ( A ) ≤ ∥ A ∥
由定义立刻可得:ρ ( A ) ≤ ∥ A ∥ \rho(A)\le\lVert A\rVert ρ ( A ) ≤ ∥ A ∥ 对任何 相容矩阵范数都成立。严格说明如下:取任意特征值 λ \lambda λ 与它对应的特征向量 v \boldsymbol{v} v ,
∣ λ ∣ ∥ v ∥ = ∥ λ v ∥ = ∥ A v ∥ ≤ ∥ A ∥ ∥ v ∥ ⟹ ∣ λ ∣ ≤ ∥ A ∥ , \lvert\lambda\rvert\lVert\boldsymbol{v}\rVert=\lVert\lambda\boldsymbol{v}\rVert=\lVert A\boldsymbol{v}\rVert\le\lVert A\rVert\lVert\boldsymbol{v}\rVert\ \Longrightarrow\ \lvert\lambda\rvert\le\lVert A\rVert, ∣ λ ∣ ∥ v ∥ = ∥ λ v ∥ = ∥ A v ∥ ≤ ∥ A ∥ ∥ v ∥ ⟹ ∣ λ ∣ ≤ ∥ A ∥ ,
两边同除以 ∥ v ∥ > 0 \lVert\boldsymbol{v}\rVert\gt0 ∥ v ∥ > 0 、再对 λ \lambda λ 取最大即得。这条不等式是「用范数估算收敛速度」的理论基础 (第五章 的迭代法分析就从它起步)。反过来一般不成立:Jordan 块 A = ( 1 1000 0 1 ) A=\begin{pmatrix}1&1000\\0&1\end{pmatrix} A = ( 1 0 1000 1 ) 的 ρ ( A ) = 1 \rho(A)=1 ρ ( A ) = 1 ,而 ∥ A ∥ 2 ≈ 1000 \lVert A\rVert_2\approx1000 ∥ A ∥ 2 ≈ 1000 ——谱半径可以远小于范数 ,这正是非正规矩阵「谱与范数对不上账」的根源。
例 2.5.1(新增) · ρ ≤ ∥ A ∥ \rho\le\lVert A\rVert ρ ≤ ∥ A ∥ 的数值核对
取三个随机 4 × 4 4\times4 4 × 4 矩阵,把 ρ ( A ) \rho(A) ρ ( A ) 与各个范数逐一比较。
解
import numpy as np
rng = np.random.default_rng(1)
for _ in range(3):
M = rng.normal(size=(4, 4))
print(f"rho={max(abs(np.linalg.eigvals(M))):.4f} "
f"||A||2={np.linalg.norm(M,2):.4f} ||A||inf={np.linalg.norm(M,np.inf):.4f} "
f"||A||1={np.linalg.norm(M,1):.4f} ||A||F={np.linalg.norm(M,'fro'):.4f}")
rho=1.6564 ||A||2=1.8582 ||A||inf=2.8008 ||A||1=3.0298 ||A||F=2.4343
rho=1.8239 ||A||2=3.3748 ||A||inf=5.1971 ||A||1=4.3686 ||A||F=4.5245
rho=1.3385 ||A||2=2.8698 ||A||inf=3.7301 ||A||1=4.4696 ||A||F=3.3676
三行里 ρ ( A ) \rho(A) ρ ( A ) 都严格小于每一个范数,而且 ρ ( A ) ≤ ∥ A ∥ 2 ≤ ∥ A ∥ F \rho(A)\le\lVert A\rVert_2\le\lVert A\rVert_F ρ ( A ) ≤ ∥ A ∥ 2 ≤ ∥ A ∥ F 与 ∥ A ∥ 2 ≤ ∥ A ∥ 1 ∥ A ∥ ∞ \lVert A\rVert_2\le\sqrt{\lVert A\rVert_1\lVert A\rVert_\infty} ∥ A ∥ 2 ≤ ∥ A ∥ 1 ∥ A ∥ ∞ 也都成立(习题 B6 已证)✓。
为什么拿 ρ ≤ ∥ A ∥ 2 \rho\le\lVert A\rVert_2 ρ ≤ ∥ A ∥ 2 当最快的估计? 因为 ∥ A ∥ 2 \lVert A\rVert_2 ∥ A ∥ 2 能直接算出来 (靠 SVD 或幂法迭代),而 ρ ( A ) \rho(A) ρ ( A ) 原则上要算出全部特征值;反过来,只要 ∥ A ∥ 2 < 1 \lVert A\rVert_2\lt1 ∥ A ∥ 2 < 1 ,立刻就能断定 ρ < 1 \rho\lt1 ρ < 1 、迭代收敛——根本不用算特征值 。这就是「用范数做收敛判据」最大的便利。
∎
直觉 · 每个特征值被「困在一个圆盘里」
特征值满足 det ( A − λ I ) = 0 \det(A-\lambda I)=0 det ( A − λ I ) = 0 。能不能把这个条件翻译成一句更直观的话——「λ \lambda λ 不会离对角元 a i i a_{ii} a ii 太远」?能。把方程 A v = λ v A\boldsymbol{v}=\lambda\boldsymbol{v} A v = λ v 改写成 ( A − λ I ) v = 0 (A-\lambda I)\boldsymbol{v}=\boldsymbol{0} ( A − λ I ) v = 0 ,挑出 v \boldsymbol{v} v 中绝对值最大的那个分量 v k v_k v k ,只看第 k k k 个方程,两边再同除以 v k v_k v k ,就得到
( λ − a k k ) v k = ∑ j ≠ k a k j v j . (\lambda-a_{kk})v_k=\sum_{j\neq k}a_{kj}v_j . ( λ − a kk ) v k = j = k ∑ a kj v j .
由于 ∣ v j / v k ∣ ≤ 1 \lvert v_j/v_k\rvert\le1 ∣ v j / v k ∣ ≤ 1 ,右端的绝对值不超过 ∑ j ≠ k ∣ a k j ∣ \sum_{j\neq k}\lvert a_{kj}\rvert ∑ j = k ∣ a kj ∣ 。于是
∣ λ − a k k ∣ ≤ ∑ j ≠ k ∣ a k j ∣ = : 第 k 行的「非对角绝对值和」 . \lvert\lambda-a_{kk}\rvert\le\sum_{j\neq k}\lvert a_{kj}\rvert=:\text{第 }k\text{ 行的「非对角绝对值和」}. ∣ λ − a kk ∣ ≤ j = k ∑ ∣ a kj ∣ =: 第 k 行的「非对角绝对值和」 .
于是每个特征值都被某个这样的圆盘困住。 这条定理的妙处在于:几乎不花计算量——只要把矩阵读一遍——就能得到特征值的可靠围栏 。工程上常用它快速判断「有没有零特征值」「是不是全在左半平面」。
定理 2.5.2(原书结论,此处给出证明) · Gerschgorin 圆盘定理
设 A = [ a i j ] ∈ R n × n A=[a_{ij}]\in\mathbb{R}^{n\times n} A = [ a ij ] ∈ R n × n ,定义 n n n 个Gerschgorin 圆盘
D i = { z ∈ C : ∣ z − a i i ∣ ≤ r i } , r i = ∑ j ≠ i ∣ a i j ∣ ( i = 1 , … , n ) , D_i=\Bigl\{z\in\mathbb{C}:\lvert z-a_{ii}\rvert\le r_i\Bigr\},\qquad r_i=\sum_{j\neq i}\lvert a_{ij}\rvert\qquad(i=1,\dots,n), D i = { z ∈ C : ∣ z − a ii ∣ ≤ r i } , r i = j = i ∑ ∣ a ij ∣ ( i = 1 , … , n ) ,
(r i r_i r i 是第 i i i 行的非对角绝对值和 ),则
(包含性) A A A 的全部特征值都包含在 ⋃ i = 1 n D i \bigcup_{i=1}^nD_i ⋃ i = 1 n D i 中;
(连通性) 若这 n n n 个圆盘中有 k k k 个构成一个与其他圆盘不相交的连通区域 ,则该区域中恰好含有 A A A 的 k k k 个特征值(重数计)。
(列版本) 把 r i r_i r i 换成列 的非对角绝对值和 ∑ j ≠ i ∣ a j i ∣ \sum_{j\neq i}\lvert a_{ji}\rvert ∑ j = i ∣ a ji ∣ (即对 A T A^{\mathrm T} A T 用上述结论),结论同样成立;两个版本的交集 可用来收紧围栏。
证明
(1) 设 λ \lambda λ 是 A A A 的特征值,v ≠ 0 \boldsymbol{v}\neq\boldsymbol{0} v = 0 满足 A v = λ v A\boldsymbol{v}=\lambda\boldsymbol{v} A v = λ v ,即对每个 i i i ,
∑ j = 1 n a i j v j = λ v i ⟹ ( λ − a i i ) v i = ∑ j ≠ i a i j v j . (2.27) \sum_{j=1}^{n}a_{ij}v_j=\lambda v_i\qquad\Longrightarrow\qquad (\lambda-a_{ii})v_i=\sum_{j\neq i}a_{ij}v_j . \tag{2.27} j = 1 ∑ n a ij v j = λ v i ⟹ ( λ − a ii ) v i = j = i ∑ a ij v j . ( 2.27 )
取 k k k 使 ∣ v k ∣ = max j ∣ v j ∣ > 0 \lvert v_k\rvert=\max_j\lvert v_j\rvert\gt0 ∣ v k ∣ = max j ∣ v j ∣ > 0 。看(2.27)在 i = k i=k i = k 的情形,两边取绝对值:
∣ λ − a k k ∣ ∣ v k ∣ = ∣ ∑ j ≠ k a k j v j ∣ ≤ ∑ j ≠ k ∣ a k j ∣ ∣ v j ∣ ≤ ( ∑ j ≠ k ∣ a k j ∣ ) ∣ v k ∣ . \lvert\lambda-a_{kk}\rvert\,\lvert v_k\rvert=\Bigl\lvert\sum_{j\neq k}a_{kj}v_j\Bigr\rvert\le\sum_{j\neq k}\lvert a_{kj}\rvert\,\lvert v_j\rvert\le\Bigl(\sum_{j\neq k}\lvert a_{kj}\rvert\Bigr)\lvert v_k\rvert . ∣ λ − a kk ∣ ∣ v k ∣ = j = k ∑ a kj v j ≤ j = k ∑ ∣ a kj ∣ ∣ v j ∣ ≤ ( j = k ∑ ∣ a kj ∣ ) ∣ v k ∣ .
约去 ∣ v k ∣ > 0 \lvert v_k\rvert\gt0 ∣ v k ∣ > 0 ,得 ∣ λ − a k k ∣ ≤ r k \lvert\lambda-a_{kk}\rvert\le r_k ∣ λ − a kk ∣ ≤ r k ,即 λ ∈ D k \lambda\in D_k λ ∈ D k 。
(2) 这是圆盘定理的「加强版」,标准证法用连续性。考虑 A t = D + t ( A − D ) A_t=D+t(A-D) A t = D + t ( A − D ) (D = diag ( a i i ) D=\operatorname{diag}(a_{ii}) D = diag ( a ii ) 是对角部分)。A 0 = D A_0=D A 0 = D 的特征值就是各个 a i i a_{ii} a ii ,即第 i i i 个圆盘的圆心 ;A 1 = A A_1=A A 1 = A 的特征值按(1)落在各圆盘内。t t t 从 0 0 0 连续增加到 1 1 1 时,特征值也随 t t t 连续变化(特征值是特征多项式系数的连续函数,而多项式的根对系数连续)。既然某个连通区域与其他圆盘不相交,特征值在连续移动中就没有办法跳过中间的空白地带,因此进入该区域的特征值个数始终不变——而 t = 0 t=0 t = 0 时该区域恰好含 k k k 个圆心,也就是 k k k 个特征值。所以 t = 1 t=1 t = 1 时它仍含 k k k 个。
(3) 对 A T A^{\mathrm T} A T 应用(1)。A T A^{\mathrm T} A T 与 A A A 有相同的特征值(定理 2.5.1 第 2 条),故结论对「列和」也成立。两个版本同时成立时,取交集即可。
∎
Gerschgorin 的三种用法
判断可逆性 / 是否有零特征值。 若所有圆盘都不含原点(0 ∉ ⋃ D i 0\notin\bigcup D_i 0 ∈ / ⋃ D i ),则 A A A 没有零特征值,因而非奇异。判断「不含原点」最常用的充分条件是严格对角占优 (strictly diagonally dominant):若 ∣ a i i ∣ > r i \lvert a_{ii}\rvert\gt r_i ∣ a ii ∣ > r i 对每一行都成立,每个圆盘就都以对角元为中心、半径小于到原点的距离,于是 0 ∉ ⋃ D i 0\notin\bigcup D_i 0 ∈ / ⋃ D i ——第五章 迭代法(Jacobi、Gauss–Seidel)的收敛性判据就建立在这条之上 。
判断特征值是否全在左半平面。 若所有圆盘都落在左半平面(Re a i i + r i < 0 \operatorname{Re}a_{ii}+r_i\lt0 Re a ii + r i < 0 ),则全部特征值的实部为负——这是动力系统 / 控制理论的稳定性判据 。
给出特征值的近似。 对角元占优时圆盘很小,a i i a_{ii} a ii 就是 λ i \lambda_i λ i 的良好近似(误差 ≤ r i \le r_i ≤ r i )。这是「对角元是特征值的粗糙近似」的定量版本。
例 2.5.2(新增) · 用 Gerschgorin 圆盘围出特征值
对
A = ( 10 − 1 2 0 4 1 − 2 1 5 ) A=\begin{pmatrix}10&-1&2\\0&4&1\\-2&1&5\end{pmatrix} A = 10 0 − 2 − 1 4 1 2 1 5
写出全部 Gerschgorin 圆盘,用行列两个版本收紧围栏,并判断 A A A 是否非奇异。
解
行圆盘:
i i i
中心 a i i a_{ii} a ii
半径 r i r_i r i (该行非对角绝对值和)
区间(实轴上)
1
10 10 10
∣ − 1 ∣ + ∣ 2 ∣ = 3 \lvert-1\rvert+\lvert2\rvert=3 ∣ − 1 ∣ + ∣ 2 ∣ = 3
[ 7 , 13 ] [7,\ 13] [ 7 , 13 ]
2
4 4 4
0 + ∣ 1 ∣ = 1 0+\lvert1\rvert=1 0 + ∣ 1 ∣ = 1
[ 3 , 5 ] [3,\ 5] [ 3 , 5 ]
3
5 5 5
∣ − 2 ∣ + ∣ 1 ∣ = 3 \lvert-2\rvert+\lvert1\rvert=3 ∣ − 2 ∣ + ∣ 1 ∣ = 3
[ 2 , 8 ] [2,\ 8] [ 2 , 8 ]
列圆盘: 半径换成「该列非对角绝对值和」——以第 3 列为例,整列绝对值和 ∣ 2 ∣ + ∣ 1 ∣ + ∣ 5 ∣ \lvert2\rvert+\lvert1\rvert+\lvert5\rvert ∣ 2 ∣ + ∣ 1 ∣ + ∣ 5 ∣ 中属于对角元的那一项不计入:
j j j
半径 c j c_j c j (该列非对角绝对值和)
区间
1
∣ 0 ∣ + ∣ − 2 ∣ = 2 \lvert0\rvert+\lvert-2\rvert=2 ∣ 0 ∣ + ∣ − 2 ∣ = 2
[ 8 , 12 ] [8,\ 12] [ 8 , 12 ]
2
∣ − 1 ∣ + ∣ 1 ∣ = 2 \lvert-1\rvert+\lvert1\rvert=2 ∣ − 1 ∣ + ∣ 1 ∣ = 2
[ 2 , 6 ] [2,\ 6] [ 2 , 6 ]
3
∣ 2 ∣ + ∣ 1 ∣ = 3 \lvert2\rvert+\lvert1\rvert=3 ∣ 2 ∣ + ∣ 1 ∣ = 3
[ 2 , 8 ] [2,\ 8] [ 2 , 8 ]
收紧后的围栏: 定理严格给出的是两个并集的交 :σ ( A ) ⊆ ( 行圆盘的并 ) ∩ ( 列圆盘的并 ) = [ 2 , 13 ] ∩ [ 2 , 12 ] = [ 2 , 12 ] \sigma(A)\subseteq\bigl(\text{行圆盘的并}\bigr)\cap\bigl(\text{列圆盘的并}\bigr)=[2,13]\cap[2,12]=[2,12] σ ( A ) ⊆ ( 行圆盘的并 ) ∩ ( 列圆盘的并 ) = [ 2 , 13 ] ∩ [ 2 , 12 ] = [ 2 , 12 ] 。至于逐盘配对取交([ 8 , 12 ] [8,12] [ 8 , 12 ] 、[ 3 , 5 ] [3,5] [ 3 , 5 ] 、[ 2 , 8 ] [2,8] [ 2 , 8 ] ),那只是启发式收紧 ,并不是定理的结论——用的时候要当心:[ 8 , 12 ] [8,12] [ 8 , 12 ] 与 [ 2 , 8 ] [2,8] [ 2 , 8 ] 在 8 8 8 处相切 ,按连通性条款它们其实连成了一块(连通性只对同一版本 的 n n n 个圆盘适用)。
老实说: 列版本的三个圆盘连成一条链,严格讲并不存在孤立分量;不过相切是「刀锋」情形——任何微小扰动都会把链断开,所以不妨把它当作一种极限情形下的直觉。数值验证(NumPy):特征值为 9.0687 9.0687 9.0687 、6.2778 6.2778 6.2778 、3.6535 3.6535 3.6535 ——9.0687 9.0687 9.0687 确实独占 [ 8 , 12 ] [8,12] [ 8 , 12 ] ,另两个落在 [ 2 , 8 ] [2,8] [ 2 , 8 ] 内,与直觉一致。
非奇异性: 原点 0 0 0 不在任何行圆盘内(最靠左的区间是 [ 2 , 8 ] [2,8] [ 2 , 8 ] ),故 A A A 非奇异。
真值(NumPy): σ ( A ) = { 3.6535 , 6.2778 , 9.0687 } \sigma(A)=\{3.6535,\ 6.2778,\ 9.0687\} σ ( A ) = { 3.6535 , 6.2778 , 9.0687 } 。逐条核对:9.0687 ∈ [ 8 , 12 ] 9.0687\in[8,12] 9.0687 ∈ [ 8 , 12 ] ✓;6.2778 ∈ [ 2 , 8 ] 6.2778\in[2,8] 6.2778 ∈ [ 2 , 8 ] ✓;3.6535 ∈ [ 3 , 5 ] 3.6535\in[3,5] 3.6535 ∈ [ 3 , 5 ] ✓;一个特征值落在孤立圆盘、两个落在连通区域 ✓。
顺带一提: A A A 其实还是严格对角占优 的——第 1 行 ∣ 10 ∣ > 3 \lvert10\rvert\gt3 ∣ 10 ∣ > 3 ✓,第 2 行 ∣ 4 ∣ > 1 \lvert4\rvert\gt1 ∣ 4 ∣ > 1 ✓,第 3 行 ∣ 5 ∣ > 3 \lvert5\rvert\gt3 ∣ 5 ∣ > 3 ✓,三行都满足 ∣ a i i ∣ > r i \lvert a_{ii}\rvert\gt r_i ∣ a ii ∣ > r i ,所以「A A A 非奇异」这个结论直接由此也能得到。
∎
直觉 · 反复乘,最大的特征值会「浮出水面」
把一个向量按 A A A 的特征向量展开:x = ∑ i c i v i \boldsymbol{x}=\sum_i c_i\boldsymbol{v}_i x = ∑ i c i v i 。然后反复乘 A A A :
A k x = ∑ i c i λ i k v i = λ 1 k [ c 1 v 1 + ∑ i > 1 c i ( λ i λ 1 ) k v i ] . A^k\boldsymbol{x}=\sum_ic_i\lambda_i^k\boldsymbol{v}_i=\lambda_1^k\Bigl[c_1\boldsymbol{v}_1+\sum_{i\gt1}c_i\Bigl(\frac{\lambda_i}{\lambda_1}\Bigr)^k\boldsymbol{v}_i\Bigr]. A k x = i ∑ c i λ i k v i = λ 1 k [ c 1 v 1 + i > 1 ∑ c i ( λ 1 λ i ) k v i ] .
提取公因子 λ 1 k \lambda_1^k λ 1 k 后,括号里除了第一项,其余各项都带着因子 ( λ i / λ 1 ) k → 0 (\lambda_i/\lambda_1)^k\to0 ( λ i / λ 1 ) k → 0 (只要 ∣ λ 1 ∣ > ∣ λ i ∣ \lvert\lambda_1\rvert\gt\lvert\lambda_i\rvert ∣ λ 1 ∣ > ∣ λ i ∣ )。于是这个向量的方向趋于 v 1 \boldsymbol{v}_1 v 1 ,而模长每乘一次就放大 ∣ λ 1 ∣ \lvert\lambda_1\rvert ∣ λ 1 ∣ 。这就是幂法 :不需要任何矩阵分解,只有反复的「乘一次、归一化一次」——只要矩阵–向量乘法足够便宜(稀疏矩阵!),它就能求出主特征值 。
收敛速度取决于 ∣ λ 2 / λ 1 ∣ \lvert\lambda_2/\lambda_1\rvert ∣ λ 2 / λ 1 ∣ (这个比值越小,收敛越快)。它也叫做分离度 ,是衡量「主特征值问题好不好算」的一把尺子。
算法 2.5.1 · 幂法(power method)+ Rayleigh 商
输入 :A ∈ R n × n A\in\mathbb{R}^{n\times n} A ∈ R n × n ,初始向量 x 0 \boldsymbol{x}_0 x 0 (u 1 T x 0 ≠ 0 \boldsymbol{u}_1^{\mathrm T}\boldsymbol{x}_0\neq0 u 1 T x 0 = 0 ),步数 k max k_{\max} k m a x ;输出 :λ 1 \lambda_1 λ 1 的近似 R ( x k ) R(\boldsymbol{x}_k) R ( x k ) 与特征向量近似 x k \boldsymbol{x}_k x k ;假设 :∣ λ 1 ∣ > ∣ λ 2 ∣ \lvert\lambda_1\rvert\gt\lvert\lambda_2\rvert ∣ λ 1 ∣ > ∣ λ 2 ∣ (主特征值单占);复杂度 :每步一次矩阵–向量乘 O ( n 2 ) O(n^2) O ( n 2 ) (A A A 稀疏时 O ( n n z ) O(\mathrm{nnz}) O ( nnz ) ),共 O ( k max n 2 ) O(k_{\max}n^2) O ( k m a x n 2 ) 。
x = 任意非零向量,规范化:x = x / ||x||
repeat:
y = A @ x # 一次矩阵-向量乘法
x_new = y / ||y|| # 规范化(防止溢出/下溢)
lambda = x_new' @ A @ x_new # Rayleigh 商 ← 更准的特征值估计
if ||x_new - x|| < tol: break
x = x_new
return lambda, x
为什么用 Rayleigh 商,而不是 ∥ y ∥ \lVert\boldsymbol{y}\rVert ∥ y ∥ 这类比值? 因为 ∥ y ∥ / ∥ x ∥ \lVert y\rVert/\lVert x\rVert ∥ y ∥ / ∥ x ∥ 只有 ∣ λ 1 ∣ \lvert\lambda_1\rvert ∣ λ 1 ∣ 的一阶 精度(误差 O ( ∣ λ 2 / λ 1 ∣ k ) O(\lvert\lambda_2/\lambda_1\rvert^k) O (∣ λ 2 / λ 1 ∣ k ) ),而 Rayleigh 商 R ( x ) = x T A x / x T x R(\boldsymbol{x})=\boldsymbol{x}^{\mathrm T}A\boldsymbol{x}/\boldsymbol{x}^{\mathrm T}\boldsymbol{x} R ( x ) = x T A x / x T x 的精度是二阶 (误差 O ( ∣ λ 2 / λ 1 ∣ 2 k ) O(\lvert\lambda_2/\lambda_1\rvert^{2k}) O (∣ λ 2 / λ 1 ∣ 2 k ) )——到达同样的精度,所需步数减半 。这是 Rayleigh 商最实际的用途。
定理 2.5.3(新增) · 幂法收敛速度与 Rayleigh 商精度
设 A A A 对称,特征值满足 ∣ λ 1 ∣ > ∣ λ 2 ∣ ≥ ⋯ ≥ ∣ λ n ∣ \lvert\lambda_1\rvert\gt\lvert\lambda_2\rvert\ge\cdots\ge\lvert\lambda_n\rvert ∣ λ 1 ∣ > ∣ λ 2 ∣ ≥ ⋯ ≥ ∣ λ n ∣ ,对应的标准正交特征向量为 u 1 , … , u n \boldsymbol{u}_1,\dots,\boldsymbol{u}_n u 1 , … , u n 。设初值 x 0 \boldsymbol{x}_0 x 0 满足 u 1 T x 0 ≠ 0 \boldsymbol{u}_1^{\mathrm T}\boldsymbol{x}_0\neq0 u 1 T x 0 = 0 ,x k = A k x 0 / ∥ A k x 0 ∥ \boldsymbol{x}_k=A^k\boldsymbol{x}_0/\lVert A^k\boldsymbol{x}_0\rVert x k = A k x 0 / ∥ A k x 0 ∥ 。则
∥ x k − u 1 ∥ 2 = O ( ∣ λ 2 λ 1 ∣ k ) , ∣ R ( x k ) − λ 1 ∣ = O ( ∣ λ 2 λ 1 ∣ 2 k ) , \lVert\boldsymbol{x}_k-\boldsymbol{u}_1\rVert_2=O\!\left(\left\lvert\frac{\lambda_2}{\lambda_1}\right\rvert^{k}\right),\qquad \lvert R(\boldsymbol{x}_k)-\lambda_1\rvert=O\!\left(\left\lvert\frac{\lambda_2}{\lambda_1}\right\rvert^{2k}\right), ∥ x k − u 1 ∥ 2 = O ( λ 1 λ 2 k ) , ∣ R ( x k ) − λ 1 ∣ = O ( λ 1 λ 2 2 k ) ,
其中 R ( x ) = x T A x x T x R(\boldsymbol{x})=\dfrac{\boldsymbol{x}^{\mathrm T}A\boldsymbol{x}}{\boldsymbol{x}^{\mathrm T}\boldsymbol{x}} R ( x ) = x T x x T A x 为 Rayleigh 商 。
证明(要点)
记 x 0 = ∑ i c i u i \boldsymbol{x}_0=\sum_ic_i\boldsymbol{u}_i x 0 = ∑ i c i u i ,c 1 ≠ 0 c_1\neq0 c 1 = 0 。则
A k x 0 = ∑ i c i λ i k u i = λ 1 k [ c 1 u 1 + ∑ i ≥ 2 c i ( λ i λ 1 ) k u i ] = λ 1 k ( c 1 u 1 + η k ) , A^k\boldsymbol{x}_0=\sum_ic_i\lambda_i^k\boldsymbol{u}_i=\lambda_1^k\Bigl[c_1\boldsymbol{u}_1+\sum_{i\ge2}c_i\Bigl(\frac{\lambda_i}{\lambda_1}\Bigr)^k\boldsymbol{u}_i\Bigr] =\lambda_1^k\bigl(c_1\boldsymbol{u}_1+\boldsymbol{\eta}_k\bigr), A k x 0 = i ∑ c i λ i k u i = λ 1 k [ c 1 u 1 + i ≥ 2 ∑ c i ( λ 1 λ i ) k u i ] = λ 1 k ( c 1 u 1 + η k ) ,
其中 ∥ η k ∥ = O ( ∣ λ 2 / λ 1 ∣ k ) \lVert\boldsymbol{\eta}_k\rVert=O(\lvert\lambda_2/\lambda_1\rvert^k) ∥ η k ∥ = O (∣ λ 2 / λ 1 ∣ k ) 。规范化后 x k = ( u 1 + c 1 − 1 η k ) / ∥ ⋅ ∥ \boldsymbol{x}_k=(\boldsymbol{u}_1+c_1^{-1}\boldsymbol{\eta}_k)/\lVert\cdot\rVert x k = ( u 1 + c 1 − 1 η k ) / ∥ ⋅ ∥ ,故方向误差为一阶量 O ( ∣ λ 2 / λ 1 ∣ k ) O(\lvert\lambda_2/\lambda_1\rvert^k) O (∣ λ 2 / λ 1 ∣ k ) 。
Rayleigh 商的部分 :设 x = u 1 + δ \boldsymbol{x}=\boldsymbol{u}_1+\boldsymbol{\delta} x = u 1 + δ 、δ ⊥ u 1 \boldsymbol{\delta}\perp\boldsymbol{u}_1 δ ⊥ u 1 。则
x T A x = ( u 1 + δ ) T A ( u 1 + δ ) = λ 1 + 2 λ 1 u 1 T δ + δ T A δ = λ 1 + δ T A δ \boldsymbol{x}^{\mathrm T}A\boldsymbol{x}=(\boldsymbol{u}_1+\boldsymbol{\delta})^{\mathrm T}A(\boldsymbol{u}_1+\boldsymbol{\delta})=\lambda_1+2\lambda_1\boldsymbol{u}_1^{\mathrm T}\boldsymbol{\delta}+\boldsymbol{\delta}^{\mathrm T}A\boldsymbol{\delta} =\lambda_1+\boldsymbol{\delta}^{\mathrm T}A\boldsymbol{\delta} x T A x = ( u 1 + δ ) T A ( u 1 + δ ) = λ 1 + 2 λ 1 u 1 T δ + δ T A δ = λ 1 + δ T A δ
(交叉项因 u 1 T A δ = λ 1 u 1 T δ = 0 \boldsymbol{u}_1^{\mathrm T}A\boldsymbol{\delta}=\lambda_1\boldsymbol{u}_1^{\mathrm T}\boldsymbol{\delta}=0 u 1 T A δ = λ 1 u 1 T δ = 0 消失)。又 x T x = 1 + 2 u 1 T δ + ∥ δ ∥ 2 = 1 + O ( ∥ δ ∥ 2 ) \boldsymbol{x}^{\mathrm T}\boldsymbol{x}=1+2\boldsymbol{u}_1^{\mathrm T}\boldsymbol{\delta}+\lVert\boldsymbol{\delta}\rVert^2=1+O(\lVert\boldsymbol{\delta}\rVert^2) x T x = 1 + 2 u 1 T δ + ∥ δ ∥ 2 = 1 + O (∥ δ ∥ 2 ) 。于是
R ( x ) = λ 1 + δ T A δ 1 + O ( ∥ δ ∥ 2 ) = λ 1 + O ( ∥ δ ∥ 2 ) . R(\boldsymbol{x})=\frac{\lambda_1+\boldsymbol{\delta}^{\mathrm T}A\boldsymbol{\delta}}{1+O(\lVert\boldsymbol{\delta}\rVert^2)}=\lambda_1+O(\lVert\boldsymbol{\delta}\rVert^2). R ( x ) = 1 + O (∥ δ ∥ 2 ) λ 1 + δ T A δ = λ 1 + O (∥ δ ∥ 2 ) .
关键就在这一点:特征向量的误差是一阶量,Rayleigh 商的误差是它的平方 ——因为 λ 1 \lambda_1 λ 1 正是 Rayleigh 商在 u 1 \boldsymbol{u}_1 u 1 处取到的驻点 ,一阶项自动消失(这也是 定理 2.5.5 的 Courant–Fischer 定理所刻画的极值性质)。
∎
例 2.5.3(新增) · 幂法的实测收敛速度
对
A = ( 5 2 0 2 4 1 0 1 3 ) A=\begin{pmatrix}5&2&0\\2&4&1\\0&1&3\end{pmatrix} A = 5 2 0 2 4 1 0 1 3
用幂法算主特征值,报告每一步的 Rayleigh 商误差,验证速度约为 ∣ λ 2 / λ 1 ∣ 2 \lvert\lambda_2/\lambda_1\rvert^2 ∣ λ 2 / λ 1 ∣ 2 。
解
σ ( A ) = { 6.6691 , 3.4760 , 1.8549 } \sigma(A)=\{6.6691,\ 3.4760,\ 1.8549\} σ ( A ) = { 6.6691 , 3.4760 , 1.8549 } (NumPy),故 ∣ λ 2 / λ 1 ∣ ≈ 0.52122 \lvert\lambda_2/\lambda_1\rvert\approx0.52122 ∣ λ 2 / λ 1 ∣ ≈ 0.52122 ,Rayleigh 商误差每步应乘约 0.5212 2 2 ≈ 0.2717 0.52122^2\approx0.2717 0.5212 2 2 ≈ 0.2717 。
import numpy as np
A = np.array([[5.,2,0],[2,4,1],[0,1,3]])
w = np.sort(np.abs(np.linalg.eigvals(A)))[::-1]
x = np.ones(3)/np.sqrt(3)
prev = None
for k in range(1, 13):
y = A @ x; x = y/np.linalg.norm(y); r = x @ A @ x
ratio = "" if prev is None else f" 误差比={abs(r-w[0])/prev:.4f}"
prev = abs(r - w[0])
if k in (1,2,3,5,8,12):
print(f"k={k:2d} Rayleigh={r:.10f} 误差={abs(r-w[0]):.3e}{ratio}")
k= 1 Rayleigh=6.5000000000 误差=1.691e-01
k= 2 Rayleigh=6.6248462485 误差=4.423e-02
k= 3 Rayleigh=6.6572150909 误差=1.186e-02
k= 5 Rayleigh=6.6682079209 误差=8.712e-04
k= 8 Rayleigh=6.6690616301 误差=1.746e-05
k=12 Rayleigh=6.6690789932 误差=9.509e-08
逐段核对比值: 4.423 × 1 0 − 2 / 1.691 × 1 0 − 1 = 0.2616 4.423\times10^{-2}/1.691\times10^{-1}=0.2616 4.423 × 1 0 − 2 /1.691 × 1 0 − 1 = 0.2616 ;1.186 × 1 0 − 2 / 4.423 × 1 0 − 2 = 0.2682 1.186\times10^{-2}/4.423\times10^{-2}=0.2682 1.186 × 1 0 − 2 /4.423 × 1 0 − 2 = 0.2682 ;8.712 × 1 0 − 4 / … ≈ 0.269 8.712\times10^{-4}/\ldots\approx0.269 8.712 × 1 0 − 4 / … ≈ 0.269 ;9.509 × 1 0 − 8 / 1.746 × 1 0 − 5 = 0.00545 9.509\times10^{-8}/1.746\times10^{-5}=0.00545 9.509 × 1 0 − 8 /1.746 × 1 0 − 5 = 0.00545 (这一段恰好是理论比值的四步复合 0.271 7 4 0.2717^4 0.271 7 4 ,仍走在理论轨道上)。
结论: 实测比值 0.26 ∼ 0.27 0.26\sim0.27 0.26 ∼ 0.27 ,与理论值 0.2717 0.2717 0.2717 吻合 ✓。这就是「收敛速度 = ( λ 2 / λ 1 ) 2 =(\lambda_2/\lambda_1)^2 = ( λ 2 / λ 1 ) 2 」的直接证据。
对比一下: 如果不用 Rayleigh 商,而用 ∥ A x k ∥ \lVert A\boldsymbol{x}_k\rVert ∥ A x k ∥ 当估计,比值会是 0.5212 0.5212 0.5212 (一阶),达到同样精度需要大约两倍的步数 。这正是 定理 2.5.3 的实用价值。
反幂法(预告): 把 A A A 换成 ( A − μ I ) − 1 (A-\mu I)^{-1} ( A − μ I ) − 1 ,主特征值就变成 1 / ( λ ∗ − μ ) 1/(\lambda^{*}-\mu) 1/ ( λ ∗ − μ ) ——迭代会收敛到最接近 μ \mu μ 的那个特征值 ,速度由 ∣ ( λ i − μ ) / ( λ j − μ ) ∣ \lvert(\lambda_i-\mu)/(\lambda_j-\mu)\rvert ∣( λ i − μ ) / ( λ j − μ )∣ 决定。这既是「已知近似值、反过来求精确特征值」的标准技巧(μ \mu μ 取 Rayleigh 商时,精度还能再升到三阶),也是第三章 3.1.2 的内容。
∎
对称矩阵在数值分析里享有「特权」:它的谱全是实数,可以正交对角化,条件数等于奇异值之比而且不需要求逆 ,扰动敏感性温和。也就是说,对称特征值问题属于「总能算得又准又稳」的那一类。为什么对称矩阵的世界特别美好?根源在于:它的特征值可以被一个优化问题刻画 ,而刻画的核心工具就是 Rayleigh 商 。下面先把 Rayleigh 商的变分性质讲清楚,再用 Courant–Fischer 定理把全部特征值一次性纳入同一个框架。
定义 2.5.2 · Rayleigh 商
设 A ∈ R n × n A\in\mathbb{R}^{n\times n} A ∈ R n × n 对称、x ≠ 0 \boldsymbol{x}\neq\boldsymbol{0} x = 0 ,称
R ( x ) = x T A x x T x (2.28) R(\boldsymbol{x})=\frac{\boldsymbol{x}^{\mathrm T}A\boldsymbol{x}}{\boldsymbol{x}^{\mathrm T}\boldsymbol{x}} \tag{2.28} R ( x ) = x T x x T A x ( 2.28 )
为 A A A 在 x \boldsymbol{x} x 处的 Rayleigh 商 。若 x \boldsymbol{x} x 是单位向量,则 R ( x ) = x T A x R(\boldsymbol{x})=\boldsymbol{x}^{\mathrm T}A\boldsymbol{x} R ( x ) = x T A x 。
直觉 · Rayleigh 商是「在某个方向上测出的特征值」
把 x \boldsymbol{x} x 归一化以后,R ( x ) R(\boldsymbol{x}) R ( x ) 就是 A A A 在方向 x \boldsymbol{x} x 上的「平均拉伸倍数」:它一定落在 λ min \lambda_{\min} λ m i n 与 λ max \lambda_{\max} λ m a x 之间,而当 x \boldsymbol{x} x 恰好是特征向量时,R ( x ) = λ R(\boldsymbol{x})=\lambda R ( x ) = λ ——分毫不差。Courant–Fischer 定理说的就是:把所有方向都试一遍,最小值是 λ min \lambda_{\min} λ m i n ,最大值是 λ max \lambda_{\max} λ m a x 。 更漂亮的是,夹在中间的那些特征值也有变分刻画——「先在 k k k 维子空间上取最小,再在所有 k k k 维子空间上取最大」,得到的就是 λ k \lambda_k λ k 。这条定理是把特征值问题变成优化问题 的桥梁,也是现代特征值算法(Lanczos、Arnoldi、子空间迭代)的全部理论基础。
定理 2.5.4(新增) · 对称矩阵 Rayleigh 商的变分性质
设 A ∈ R n × n A\in\mathbb{R}^{n\times n} A ∈ R n × n 对称,特征值(从小到大)λ 1 ≤ λ 2 ≤ ⋯ ≤ λ n \lambda_1\le\lambda_2\le\cdots\le\lambda_n λ 1 ≤ λ 2 ≤ ⋯ ≤ λ n ,对应的标准正交特征向量 u 1 , … , u n \boldsymbol{u}_1,\dots,\boldsymbol{u}_n u 1 , … , u n ,则
λ 1 = min x ≠ 0 R ( x ) , λ n = max x ≠ 0 R ( x ) . \lambda_1=\min_{\boldsymbol{x}\neq\boldsymbol{0}}R(\boldsymbol{x}),\qquad \lambda_n=\max_{\boldsymbol{x}\neq\boldsymbol{0}}R(\boldsymbol{x}). λ 1 = x = 0 min R ( x ) , λ n = x = 0 max R ( x ) .
且 x \boldsymbol{x} x 是极值点 ⟺ \iff ⟺ x \boldsymbol{x} x 是对应 λ 1 \lambda_1 λ 1 (或 λ n \lambda_n λ n )的特征向量。
证明
任何 x ≠ 0 \boldsymbol{x}\neq\boldsymbol{0} x = 0 可展开为 x = ∑ i k i u i \boldsymbol{x}=\sum_ik_i\boldsymbol{u}_i x = ∑ i k i u i ,由正交性 x T x = ∑ i k i 2 \boldsymbol{x}^{\mathrm T}\boldsymbol{x}=\sum_ik_i^2 x T x = ∑ i k i 2 ,且
x T A x = ∑ i , j k i k j u i T A u j = ∑ i , j k i k j λ j δ i j = ∑ i λ i k i 2 . \boldsymbol{x}^{\mathrm T}A\boldsymbol{x}=\sum_{i,j}k_ik_j\boldsymbol{u}_i^{\mathrm T}A\boldsymbol{u}_j=\sum_{i,j}k_ik_j\lambda_j\delta_{ij}=\sum_i\lambda_ik_i^2 . x T A x = i , j ∑ k i k j u i T A u j = i , j ∑ k i k j λ j δ ij = i ∑ λ i k i 2 .
于是
R ( x ) = ∑ i λ i k i 2 ∑ i k i 2 . R(\boldsymbol{x})=\frac{\sum_i\lambda_ik_i^2}{\sum_ik_i^2} . R ( x ) = ∑ i k i 2 ∑ i λ i k i 2 .
一旦看出它是加权平均,结论就自动出来了: R ( x ) R(\boldsymbol{x}) R ( x ) 是以 k i 2 / ( ∑ j k j 2 ) k_i^2/(\sum_jk_j^2) k i 2 / ( ∑ j k j 2 ) 为权的 λ i \lambda_i λ i 的加权平均,必然落在 [ λ 1 , λ n ] [\lambda_1,\lambda_n] [ λ 1 , λ n ] 内。等号出现在下界,当且仅当全部权重都集中到 λ 1 \lambda_1 λ 1 上,即 x ∈ span { u 1 } \boldsymbol{x}\in\operatorname{span}\{\boldsymbol{u}_1\} x ∈ span { u 1 } ;上界的情形同理。
∎
定理 2.5.5(新增) · Courant–Fischer 极小极大原理
设 A A A 对称,特征值从小到大为 λ 1 ≤ ⋯ ≤ λ n \lambda_1\le\cdots\le\lambda_n λ 1 ≤ ⋯ ≤ λ n ,则对 k = 1 , … , n k=1,\dots,n k = 1 , … , n ,
λ k = min dim S = k max 0 ≠ x ∈ S R ( x ) = max dim T = n − k + 1 min 0 ≠ x ∈ T R ( x ) , \lambda_k=\min_{\dim\mathcal{S}=k}\ \max_{\boldsymbol{0}\neq\boldsymbol{x}\in\mathcal{S}}R(\boldsymbol{x}) =\max_{\dim\mathcal{T}=n-k+1}\ \min_{\boldsymbol{0}\neq\boldsymbol{x}\in\mathcal{T}}R(\boldsymbol{x}) , λ k = d i m S = k min 0 = x ∈ S max R ( x ) = d i m T = n − k + 1 max 0 = x ∈ T min R ( x ) ,
其中 S , T \mathcal{S},\mathcal{T} S , T 取遍 R n \mathbb{R}^n R n 的相应维数子空间。
证明(第一等式)
先证 ≥ \ge ≥ (取特定的 k k k 维子空间)。 令 S 0 = span { u 1 , … , u k } \mathcal{S}_0=\operatorname{span}\{\boldsymbol{u}_1,\dots,\boldsymbol{u}_k\} S 0 = span { u 1 , … , u k } 。对任意 x ∈ S 0 \boldsymbol{x}\in\mathcal{S}_0 x ∈ S 0 、x = ∑ i ≤ k k i u i \boldsymbol{x}=\sum_{i\le k}k_i\boldsymbol{u}_i x = ∑ i ≤ k k i u i ,有 R ( x ) ≤ λ k R(\boldsymbol{x})\le\lambda_k R ( x ) ≤ λ k (加权平均的权只落在 λ 1 , … , λ k \lambda_1,\dots,\lambda_k λ 1 , … , λ k 上)。故
max x ∈ S 0 R ( x ) = λ k ( 取 x = u k 达到 ) ⟹ min dim S = k max x ∈ S R ( x ) ≤ λ k . \max_{\boldsymbol{x}\in\mathcal{S}_0}R(\boldsymbol{x})=\lambda_k\qquad(\text{取}\ \boldsymbol{x}=\boldsymbol{u}_k\ \text{达到})\quad\Longrightarrow\quad \min_{\dim\mathcal{S}=k}\max_{\boldsymbol{x}\in\mathcal{S}}R(\boldsymbol{x})\le\lambda_k . x ∈ S 0 max R ( x ) = λ k ( 取 x = u k 达到 ) ⟹ d i m S = k min x ∈ S max R ( x ) ≤ λ k .
再证 ≤ \le ≤ (维数计数)。 设 S \mathcal{S} S 是任意 k k k 维子空间,令 T = span { u k , … , u n } \mathcal{T}=\operatorname{span}\{\boldsymbol{u}_k,\dots,\boldsymbol{u}_n\} T = span { u k , … , u n } ,dim T = n − k + 1 \dim\mathcal{T}=n-k+1 dim T = n − k + 1 。由维数公式
dim ( S ∩ T ) ≥ dim S + dim T − n = k + ( n − k + 1 ) − n = 1 , \dim(\mathcal{S}\cap\mathcal{T})\ge\dim\mathcal{S}+\dim\mathcal{T}-n=k+(n-k+1)-n=1, dim ( S ∩ T ) ≥ dim S + dim T − n = k + ( n − k + 1 ) − n = 1 ,
故存在 0 ≠ x ∈ S ∩ T \boldsymbol{0}\neq\boldsymbol{x}\in\mathcal{S}\cap\mathcal{T} 0 = x ∈ S ∩ T 。写 x = ∑ i ≥ k k i u i \boldsymbol{x}=\sum_{i\ge k}k_i\boldsymbol{u}_i x = ∑ i ≥ k k i u i ,则 R ( x ) ≥ λ k R(\boldsymbol{x})\ge\lambda_k R ( x ) ≥ λ k (权只落在 λ k , … , λ n \lambda_k,\dots,\lambda_n λ k , … , λ n 上)。于是在这个 S \mathcal{S} S 上 max R ≥ λ k \max R\ge\lambda_k max R ≥ λ k ;S \mathcal{S} S 是任意的,对 S \mathcal{S} S 取下确界后仍 ≥ λ k \ge\lambda_k ≥ λ k 。
两个方向合起来就是第一等式。至于第二等式,对 − A -A − A 套用第一等式即可(( − A ) (-A) ( − A ) 的特征值是 − λ n ≤ ⋯ ≤ − λ 1 -\lambda_n\le\cdots\le-\lambda_1 − λ n ≤ ⋯ ≤ − λ 1 ,其中第 k k k 大 = − λ n − k + 1 =-\lambda_{n-k+1} = − λ n − k + 1 ,再用 min max \min\max min max 与 max min \max\min max min 的互换)。
∎
例 2.5.4(新增) · 数值验证 Courant–Fischer
对
A = ( 4 1 0 1 3 1 0 1 2 ) A=\begin{pmatrix}4&1&0\\1&3&1\\0&1&2\end{pmatrix} A = 4 1 0 1 3 1 0 1 2
数值验证:随机方向上的 Rayleigh 商极值趋于 λ min , λ max \lambda_{\min},\lambda_{\max} λ m i n , λ m a x ;k k k 维子空间上的极小值最大化后趋于 λ 2 \lambda_2 λ 2 。
解
σ ( A ) = { 1.2679 , 3.0000 , 4.7321 } \sigma(A)=\{1.2679,\ 3.0000,\ 4.7321\} σ ( A ) = { 1.2679 , 3.0000 , 4.7321 } (从小到大)。
import numpy as np
A = np.array([[4.,1,0],[1,3,1],[0,1,2]])
w = np.linalg.eigvalsh(A) # 升序
rng = np.random.default_rng(0)
# (a) 随机单位方向上的 Rayleigh 商极值
vals = []
for _ in range(200000):
x = rng.normal(size=3); x /= np.linalg.norm(x); vals.append(x @ A @ x)
print(f"随机方向: min R={min(vals):.8f} max R={max(vals):.8f}")
print(f"理论: λmin={w[0]:.8f} λmax={w[-1]:.8f}")
# (b) 2 维子空间上的极小值,再在所有 2 维子空间上取最大
best = -np.inf
for _ in range(200000):
Q, _ = np.linalg.qr(rng.normal(size=(3, 2))) # 随机 2 维子空间的正交基
best = max(best, np.linalg.eigvalsh(Q.T @ A @ Q)[0]) # 该子空间上的 min R
print(f"max over 2-dim subspaces of min R = {best:.8f} (理论 λ2={w[1]:.8f})")
随机方向: min R=1.26795245 max R=4.73203903
理论: λmin=1.26794919 λmax=4.73205081
max over 2-dim subspaces of min R = 3.00000000 (理论 λ2=3.00000000)
从结果里读出三件事:
随机方向上的极值只是逼近 λ min , λ max \lambda_{\min},\lambda_{\max} λ m i n , λ m a x (差 ∼ 1 0 − 5 \sim10^{-5} ∼ 1 0 − 5 ):20 万次随机采样并没有命中精确的特征方向——变分原理的本质是优化 ,随机撒点只不过是一种近似的优化方式。
第二项精确等于 λ 2 = 3 \lambda_2=3 λ 2 = 3 :取等的子空间是 S = span { u 2 , u 3 } \mathcal{S}=\operatorname{span}\{\boldsymbol{u}_2,\boldsymbol{u}_3\} S = span { u 2 , u 3 } (其上 min R = λ 2 \min R=\lambda_2 min R = λ 2 );而任何 2 维子空间都与 span { u 1 , u 2 } \operatorname{span}\{\boldsymbol{u}_1,\boldsymbol{u}_2\} span { u 1 , u 2 } 有非零交(维数公式 2 + 2 > 3 2+2\gt3 2 + 2 > 3 ),在该交集上 R ≤ λ 2 R\le\lambda_2 R ≤ λ 2 ,故 max \max max 不可能超过 λ 2 \lambda_2 λ 2 。这就是 Courant–Fischer 的取等构造在起作用。
实用含义: 求 λ 2 \lambda_2 λ 2 可以「在 2 维子空间上做 Rayleigh 商最小化」——这就是 Lanczos 算法 / 子空间迭代 的基本思想,也是 习题 C6 的动手内容。
∎
这是本节最重要、也最反直觉的一节。前面反复出现的「对称矩阵很温和」这句话,到这里要兑现成定量的结论,同时也把「非正规矩阵有多难缠」摆到台面上。
直觉 · 重特征值会被扰动「撕开」成 ε \sqrt{\varepsilon} ε 或更糟
设想 A A A 有一个重特征值 (例如 λ = 1 \lambda=1 λ = 1 ,二重)。如果往矩阵上加一个大小为 ε \varepsilon ε 的扰动,特征值会跑多远?
对称(或正规)情形: 移动量是 O ( ε ) O(\varepsilon) O ( ε ) ——温和。
非正规情形(如 Jordan 块): 移动量是 O ( ε ) O(\sqrt{\varepsilon}) O ( ε ) ——开方敏感 !ε = 1 0 − 16 \varepsilon=10^{-16} ε = 1 0 − 16 时移动 1 0 − 8 10^{-8} 1 0 − 8 :丢了 8 位有效数字 。
为什么会开方? 把 A = ( 1 1 0 1 ) A=\begin{pmatrix}1&1\\0&1\end{pmatrix} A = ( 1 0 1 1 ) 扰动成 ( 1 1 t 1 ) \begin{pmatrix}1&1\\t&1\end{pmatrix} ( 1 t 1 1 ) ,特征方程是 ( 1 − λ ) 2 = t (1-\lambda)^2=t ( 1 − λ ) 2 = t ,于是 λ = 1 ± t \lambda=1\pm\sqrt{t} λ = 1 ± t 。重根处的特征多项式是「平方」形态,解出来自然带根号。 重数越高越糟:k k k 重特征值遇上 O ( ε ) O(\varepsilon) O ( ε ) 量级的扰动,可以移动 O ( ε 1 / k ) O(\varepsilon^{1/k}) O ( ε 1/ k ) 。
工程含义: 别指望用数值方法算准重特征值——重特征值问题是固有病态的 (1.2.6 的 conditioning 概念)。反过来,对称矩阵即使有重特征值也照样温和 :它的特征值由 Rayleigh 商的极值刻画(定理 2.5.4 ),而极值对扰动的响应是 Lipschitz 的——扰动有多大,极值最多也就变多大。
定理 2.5.6(新增) · 特征值扰动界(两个极端情形)
(一)对称情形(Weyl 不等式): 设 A , E ∈ R n × n A,E\in\mathbb{R}^{n\times n} A , E ∈ R n × n 对称,特征值分别按升序为 λ 1 ≤ ⋯ ≤ λ n \lambda_1\le\cdots\le\lambda_n λ 1 ≤ ⋯ ≤ λ n 与 μ 1 ≤ ⋯ ≤ μ n \mu_1\le\cdots\le\mu_n μ 1 ≤ ⋯ ≤ μ n ,则对每个 i i i ,
∣ λ i − μ i ∣ ≤ ∥ E ∥ 2 . \lvert\lambda_i-\mu_i\rvert\le\lVert E\rVert_2 . ∣ λ i − μ i ∣ ≤ ∥ E ∥ 2 .
即对称矩阵的特征值对扰动是 1 1 1 -Lipschitz 的 ——误差与扰动同阶,一点都不放大。
(二)一般情形(Bauer–Fike 定理): 设 A A A 可对角化,A = V Λ V − 1 A=V\Lambda V^{-1} A = V Λ V − 1 ,E E E 为扰动,则 σ ( A + E ) \sigma(A+E) σ ( A + E ) 中的每个特征值 μ \mu μ 都满足
min λ ∈ σ ( A ) ∣ λ − μ ∣ ≤ ∥ V ∥ 2 ∥ V − 1 ∥ 2 ∥ E ∥ 2 = c o n d 2 ( V ) ∥ E ∥ 2 . \min_{\lambda\in\sigma(A)}\lvert\lambda-\mu\rvert\le\lVert V\rVert_2\,\lVert V^{-1}\rVert_2\,\lVert E\rVert_2=\mathrm{cond}_2(V)\,\lVert E\rVert_2 . λ ∈ σ ( A ) min ∣ λ − μ ∣ ≤ ∥ V ∥ 2 ∥ V − 1 ∥ 2 ∥ E ∥ 2 = cond 2 ( V ) ∥ E ∥ 2 .
放大因子就是特征向量矩阵 V V V 的条件数。 对称矩阵的 V V V 可以取成正交矩阵(此时 c o n d 2 ( V ) = 1 \mathrm{cond}_2(V)=1 cond 2 ( V ) = 1 ),界便退化为「扰动的一倍」——这正是情形(一)的来源。而非正规矩阵的特征向量矩阵 V V V 可以病态到极点 ,能把 O ( ε ) O(\varepsilon) O ( ε ) 的扰动放大成巨大的特征值移动。
(三)重特征值情形(无一般界!) 若 A A A 有 k k k 重特征值且非正规(Jordan 块),扰动 ∥ E ∥ = ε \lVert E\rVert=\varepsilon ∥ E ∥ = ε 可使特征值移动 O ( ε 1 / k ) O(\varepsilon^{1/k}) O ( ε 1/ k ) ——这是无界的放大因子(ε 1 / k / ε = ε 1 / k − 1 → ∞ \varepsilon^{1/k}/\varepsilon=\varepsilon^{1/k-1}\to\infty ε 1/ k / ε = ε 1/ k − 1 → ∞ ) ,且任何算法都无法改善 。
证明
(一)Weyl 不等式。 用 Courant–Fischer(定理 2.5.5 ):
λ i = min dim S = i max x ∈ S R A ( x ) , μ i = min dim S = i max x ∈ S R A + E ( x ) . \lambda_i=\min_{\dim\mathcal{S}=i}\max_{\boldsymbol{x}\in\mathcal{S}}R_A(\boldsymbol{x}),\qquad \mu_i=\min_{\dim\mathcal{S}=i}\max_{\boldsymbol{x}\in\mathcal{S}}R_{A+E}(\boldsymbol{x}) . λ i = d i m S = i min x ∈ S max R A ( x ) , μ i = d i m S = i min x ∈ S max R A + E ( x ) .
对任意固定的 S \mathcal{S} S 与单位向量 x ∈ S \boldsymbol{x}\in\mathcal{S} x ∈ S ,∣ R A ( x ) − R A + E ( x ) ∣ = ∣ x T E x ∣ ≤ ∥ E ∥ 2 \lvert R_{A}(\boldsymbol{x})-R_{A+E}(\boldsymbol{x})\rvert=\lvert\boldsymbol{x}^{\mathrm T}E\boldsymbol{x}\rvert\le\lVert E\rVert_2 ∣ R A ( x ) − R A + E ( x )∣ = ∣ x T E x ∣ ≤ ∥ E ∥ 2 。因此对每个固定的 S \mathcal{S} S ,两个「max \max max 」之差都不超过 ∥ E ∥ 2 \lVert E\rVert_2 ∥ E ∥ 2 ;再对 S \mathcal{S} S 取 min \min min ,同样得到 ∣ λ i − μ i ∣ ≤ ∥ E ∥ 2 \lvert\lambda_i-\mu_i\rvert\le\lVert E\rVert_2 ∣ λ i − μ i ∣ ≤ ∥ E ∥ 2 。
(二)Bauer–Fike。 设 ( A + E ) w = μ w (A+E)\boldsymbol{w}=\mu\boldsymbol{w} ( A + E ) w = μ w 、∥ w ∥ 2 = 1 \lVert\boldsymbol{w}\rVert_2=1 ∥ w ∥ 2 = 1 。则
( A − μ I ) w = − E w , w = − ( A − μ I ) − 1 E w ( μ ∉ σ ( A ) 时 ) . (A-\mu I)\boldsymbol{w}=-E\boldsymbol{w},\qquad \boldsymbol{w}=-(A-\mu I)^{-1}E\boldsymbol{w}\ \ (\mu\notin\sigma(A)\ \text{时}) . ( A − μ I ) w = − E w , w = − ( A − μ I ) − 1 E w ( μ ∈ / σ ( A ) 时 ) .
设 μ \mu μ 使 μ I − A \mu I-A μ I − A 可逆(否则 μ \mu μ 已是 A A A 的特征值,结论显然)。取范数:
1 = ∥ w ∥ ≤ ∥ ( A − μ I ) − 1 ∥ 2 ∥ E ∥ 2 . 1=\lVert\boldsymbol{w}\rVert\le\lVert(A-\mu I)^{-1}\rVert_2\lVert E\rVert_2 . 1 = ∥ w ∥ ≤ ∥( A − μ I ) − 1 ∥ 2 ∥ E ∥ 2 .
关键一步:估计 ∥ ( A − μ I ) − 1 ∥ 2 \lVert(A-\mu I)^{-1}\rVert_2 ∥( A − μ I ) − 1 ∥ 2 。 用 A = V Λ V − 1 A=V\Lambda V^{-1} A = V Λ V − 1 :
∥ ( A − μ I ) − 1 ∥ 2 = ∥ V ( Λ − μ I ) − 1 V − 1 ∥ 2 ≤ ∥ V ∥ 2 ∥ V − 1 ∥ 2 max i 1 ∣ λ i − μ ∣ = c o n d 2 ( V ) min i ∣ λ i − μ ∣ . \lVert(A-\mu I)^{-1}\rVert_2=\lVert V(\Lambda-\mu I)^{-1}V^{-1}\rVert_2\le\lVert V\rVert_2\lVert V^{-1}\rVert_2\max_i\frac{1}{\lvert\lambda_i-\mu\rvert} =\frac{\mathrm{cond}_2(V)}{\min_i\lvert\lambda_i-\mu\rvert} . ∥( A − μ I ) − 1 ∥ 2 = ∥ V ( Λ − μ I ) − 1 V − 1 ∥ 2 ≤ ∥ V ∥ 2 ∥ V − 1 ∥ 2 i max ∣ λ i − μ ∣ 1 = min i ∣ λ i − μ ∣ cond 2 ( V ) .
代入上式得 1 ≤ c o n d 2 ( V ) ∥ E ∥ 2 min i ∣ λ i − μ ∣ 1\le\dfrac{\mathrm{cond}_2(V)\lVert E\rVert_2}{\min_i\lvert\lambda_i-\mu\rvert} 1 ≤ min i ∣ λ i − μ ∣ cond 2 ( V ) ∥ E ∥ 2 ,整理后就是待证的结论。
(三) 用一个 Jordan 块把现象显式摆出来:A = ( 1 1 0 1 ) A=\begin{pmatrix}1&1\\0&1\end{pmatrix} A = ( 1 0 1 1 ) (k = 2 k=2 k = 2 重特征值 1 1 1 ),E = ( 0 0 t 0 ) E=\begin{pmatrix}0&0\\t&0\end{pmatrix} E = ( 0 t 0 0 ) ,则 A + E A+E A + E 的特征方程 ( 1 − λ ) 2 = t (1-\lambda)^2=t ( 1 − λ ) 2 = t 给出 λ = 1 ± t \lambda=1\pm\sqrt{t} λ = 1 ± t ,移动量就是 t = ε 1 / 2 \sqrt{t}=\varepsilon^{1/2} t = ε 1/2 (取 ε = t \varepsilon=t ε = t )。
∎
例 2.5.5(新增) · 开方敏感 vs 温和:数值对照
(1)对 Jordan 块 A = ( 1 1 0 1 ) A=\begin{pmatrix}1&1\\0&1\end{pmatrix} A = ( 1 0 1 1 ) 加扰动 E = ( 0 0 t 0 ) E=\begin{pmatrix}0&0\\t&0\end{pmatrix} E = ( 0 t 0 0 ) ,看特征值的移动怎样随 t t t 变化。
(2)对对称矩阵 A = ( 1 0 0 2 ) A=\begin{pmatrix}1&0\\0&2\end{pmatrix} A = ( 1 0 0 2 ) 加扰动 E = ( 0 t t 0 ) E=\begin{pmatrix}0&t\\t&0\end{pmatrix} E = ( 0 t t 0 ) ,同样观察一遍。
解
import numpy as np
print("== (1) Jordan 块(非正规)==")
J = np.array([[1.,1],[0,1.]])
for t in [1e-2,1e-4,1e-8,1e-12,1e-16]:
B = J + np.array([[0.,0],[t,0]])
print(f" t={t:.0e} 特征值={np.round(np.linalg.eigvals(B),10)} "
f"1±sqrt(t)={1-np.sqrt(t):.10f},{1+np.sqrt(t):.10f}")
print("== (2) 对称矩阵 ==")
S = np.array([[1.,0],[0,2.]])
for t in [1e-2,1e-8,1e-16]:
B = S + np.array([[0.,t],[t,0]])
print(f" t={t:.0e} 特征值={np.round(np.linalg.eigvalsh(B),14)} (移动量 ~t)")
== (1) Jordan 块(非正规)==
t=1e-02 特征值=[1.1 0.9] 1±sqrt(t)=0.9000000000,1.1000000000
t=1e-04 特征值=[1.01 0.99] 1±sqrt(t)=0.9900000000,1.0100000000
t=1e-08 特征值=[1.0001 0.9999] 1±sqrt(t)=0.9999000000,1.0001000000
t=1e-12 特征值=[1.000001 0.999999] 1±sqrt(t)=0.9999990000,1.0000010000
t=1e-16 特征值=[1. 1.] 1±sqrt(t)=0.9999999900,1.0000000100
== (2) 对称矩阵 ==
t=1e-02 特征值=[0.99990001 2.00009999] (移动量 ≈t²)
t=1e-08 特征值=[1. 2.]
t=1e-16 特征值=[1. 2.]
(1)的读法: 移动量精确等于 t \sqrt{t} t 。t = 1 0 − 2 t=10^{-2} t = 1 0 − 2 时移动 1 0 − 1 10^{-1} 1 0 − 1 、t = 1 0 − 12 t=10^{-12} t = 1 0 − 12 时移动 1 0 − 6 10^{-6} 1 0 − 6 ——扰动缩小了 10 个数量级,误差才从 1 0 − 1 10^{-1} 1 0 − 1 压到 1 0 − 6 10^{-6} 1 0 − 6 (只缩了 5 个数量级) ,这就是「开方」的含义。t = 1 0 − 16 t=10^{-16} t = 1 0 − 16 时(t = 1 0 − 8 \sqrt{t}=10^{-8} t = 1 0 − 8 )双精度仍能分辨 1 ± 1 0 − 8 1\pm10^{-8} 1 ± 1 0 − 8 (1 1 1 附近的机器精度是 2 × 1 0 − 16 2\times10^{-16} 2 × 1 0 − 16 ),所以这一行依然正确;若取 t = 1 0 − 20 t=10^{-20} t = 1 0 − 20 ,1 ± 1 0 − 10 1\pm10^{-10} 1 ± 1 0 − 10 就再也算不出来了——在重特征值附近,双精度只能保证约 8 位有效数字 (精度是 u \sqrt{u} u 而不是 u u u )。
(2)的读法: 特征值精确为 1.5 ± 0.25 + t 2 1.5\pm\sqrt{0.25+t^2} 1.5 ± 0.25 + t 2 ,移动量约 t 2 t^2 t 2 ——比 t t t 还小一个数量级 (t = 1 0 − 2 t=10^{-2} t = 1 0 − 2 时只移动 1 0 − 4 10^{-4} 1 0 − 4 ,见表中 0.99990001 0.99990001 0.99990001 )。t = 1 0 − 8 t=10^{-8} t = 1 0 − 8 与 1 0 − 16 10^{-16} 1 0 − 16 时,移动量(∼ 1 0 − 16 \sim10^{-16} ∼ 1 0 − 16 与 ∼ 1 0 − 32 \sim10^{-32} ∼ 1 0 − 32 )已经低于机器精度,完全测不出来 ——对称矩阵的特征值对扰动极其稳健 。
工程结论: ① 对称特征值问题(结构力学、协方差矩阵、量子力学)总能算得很准;② 非对称 / 非正规矩阵的重特征值问题(如流体稳定性、某些控制问题)在数值上极难处理 ——这不是算法不行,是问题本身病态;③ 实践中应尽量避免 让矩阵既非正规又带重特征值,或者改用 Schur 形式 / 奇异值(σ i = λ i ( A T A ) \sigma_i=\sqrt{\lambda_i(A^{\mathrm T}A)} σ i = λ i ( A T A ) 总在对称 矩阵上算,原理上更稳)。
∎
例 2.5.6(新增) · 非正规矩阵的「伪谱」直觉(可选)
取 N = ( 1 1000 0 1 ) N=\begin{pmatrix}1&1000\\0&1\end{pmatrix} N = ( 1 0 1000 1 ) ,它的 ρ ( N ) = 1 \rho(N)=1 ρ ( N ) = 1 、∥ N ∥ 2 ≈ 1000 \lVert N\rVert_2\approx1000 ∥ N ∥ 2 ≈ 1000 。请解释:为什么「ρ ( N ) < 1 \rho(N)\lt1 ρ ( N ) < 1 就意味着迭代收敛」这句话,在非正规矩阵上不能 按直觉来用。
解
ρ ( N ) = 1 \rho(N)=1 ρ ( N ) = 1 (特征值 1 , 1 1,1 1 , 1 ),∥ N ∥ 2 ≈ 1000 \lVert N\rVert_2\approx1000 ∥ N ∥ 2 ≈ 1000 。谱半径远小于范数 。
问题出在「迭代 N k N^k N k 到底衰不衰减」上。 对正规矩阵,∥ A k ∥ 2 = ρ ( A ) k = ∣ λ max ∣ k \lVert A^k\rVert_2=\rho(A)^k=\lvert\lambda_{\max}\rvert^k ∥ A k ∥ 2 = ρ ( A ) k = ∣ λ m a x ∣ k ,所以 ρ < 1 \rho\lt1 ρ < 1 立刻给出几何衰减。对非正规矩阵,只能保证 ρ ( A ) = lim k ∥ A k ∥ 1 / k \rho(A)=\lim_k\lVert A^k\rVert^{1/k} ρ ( A ) = lim k ∥ A k ∥ 1/ k (Gelfand 公式,证明与瞬态增长的定量分析见附录 A.2.7 )——这个极限式子固然成立,但有限步内的 ∥ N k ∥ \lVert N^k\rVert ∥ N k ∥ 可以要多大有多大 。本例:
N k = ( 1 k ⋅ 1000 0 1 ) , ∥ N k ∥ 2 ≈ 1000 k ( 线性增长! ) N^k=\begin{pmatrix}1&k\cdot1000\\0&1\end{pmatrix},\qquad \lVert N^k\rVert_2\approx1000k\ (\text{线性增长!}) N k = ( 1 0 k ⋅ 1000 1 ) , ∥ N k ∥ 2 ≈ 1000 k ( 线性增长! )
即使 ρ ( N ) = 1 \rho(N)=1 ρ ( N ) = 1 (注意还不是 < 1 \lt1 < 1 ),∥ N k ∥ \lVert N^k\rVert ∥ N k ∥ 也能增长到任意大 ;哪怕取 ρ ( N ) = 0.9 < 1 \rho(N)=0.9\lt1 ρ ( N ) = 0.9 < 1 (例如 ( 0.9 M 0 0.9 ) \begin{pmatrix}0.9&M\\0&0.9\end{pmatrix} ( 0.9 0 M 0.9 ) ),在 k ≲ 1 / log ( 1 / 0.9 ) k\lesssim1/\log(1/0.9) k ≲ 1/ log ( 1/0.9 ) 的步数内,∥ N k ∥ \lVert N^k\rVert ∥ N k ∥ 仍可能先上升——这就是「瞬态增长」(transient growth) 。
工程含义: 第五章 的迭代法收敛性判据如果只写「ρ ( H ) < 1 \rho(H)\lt1 ρ ( H ) < 1 」,在非正规情形下就是理论上成立、实践中可能失效 (有限步内误差先涨后落,甚至溢出)。更稳妥的判据是 ∥ H ∥ < 1 \lVert H\rVert\lt1 ∥ H ∥ < 1 (1.3.5 定理 1.3.6),它给出每一步都在收缩 的保证。这就是「用范数,不用谱半径」的实用理由。
(这类现象在流体力学与数值天气预报里尤其常见,有专门的工具对付它——伪谱(pseudospectrum) ,见 延伸阅读 。)
∎
本节公式速查
概念
公式 / 结论
特征值
A v = λ v A\boldsymbol{v}=\lambda\boldsymbol{v} A v = λ v ;det ( A − λ I ) = 0 \det(A-\lambda I)=0 det ( A − λ I ) = 0
迹 / 行列式
tr A = ∑ λ i \operatorname{tr}A=\sum\lambda_i tr A = ∑ λ i ;det A = ∏ λ i \det A=\prod\lambda_i det A = ∏ λ i
谱半径
ρ ( A ) = max ∣ λ i ∣ \rho(A)=\max\lvert\lambda_i\rvert ρ ( A ) = max ∣ λ i ∣
谱半径与范数
ρ ( A ) ≤ ∥ A ∥ \rho(A)\le\lVert A\rVert ρ ( A ) ≤ ∥ A ∥ (任一相容范数)
Gerschgorin
∣ λ − a i i ∣ ≤ ∑ j ≠ i ∣ a i j ∣ \lvert\lambda-a_{ii}\rvert\le\sum_{j\neq i}\lvert a_{ij}\rvert ∣ λ − a ii ∣ ≤ ∑ j = i ∣ a ij ∣ ;连通区域含对应个数特征值
对角占优
∣ a i i ∣ > ∑ j ≠ i ∣ a i j ∣ ∀ i ⇒ A \lvert a_{ii}\rvert\gt\sum_{j\neq i}\lvert a_{ij}\rvert\ \forall i\Rightarrow A ∣ a ii ∣ > ∑ j = i ∣ a ij ∣ ∀ i ⇒ A 非奇异
Rayleigh 商
R ( x ) = x T A x / x T x R(\boldsymbol{x})=\boldsymbol{x}^{\mathrm T}A\boldsymbol{x}/\boldsymbol{x}^{\mathrm T}\boldsymbol{x} R ( x ) = x T A x / x T x (A A A 对称)
变分刻画
λ 1 = min R \lambda_1=\min R λ 1 = min R ,λ n = max R \lambda_n=\max R λ n = max R ;Courant–Fischer 给全部中间特征值
幂法收敛
向量误差 O ( ∣ λ 2 / λ 1 ∣ k ) O(\lvert\lambda_2/\lambda_1\rvert^k) O (∣ λ 2 / λ 1 ∣ k ) ;Rayleigh 商误差 O ( ∣ λ 2 / λ 1 ∣ 2 k ) O(\lvert\lambda_2/\lambda_1\rvert^{2k}) O (∣ λ 2 / λ 1 ∣ 2 k )
对称扰动
∣ λ i ( A + E ) − λ i ( A ) ∣ ≤ ∥ E ∥ 2 \lvert\lambda_i(A+E)-\lambda_i(A)\rvert\le\lVert E\rVert_2 ∣ λ i ( A + E ) − λ i ( A )∣ ≤ ∥ E ∥ 2 (Weyl)
一般扰动
≤ c o n d 2 ( V ) ∥ E ∥ 2 \le\mathrm{cond}_2(V)\lVert E\rVert_2 ≤ cond 2 ( V ) ∥ E ∥ 2 (Bauer–Fike)
重特征值
非正规 k k k 重:移动 O ( ∥ E ∥ 1 / k ) O(\lVert E\rVert^{1/k}) O (∥ E ∥ 1/ k ) ——固有病态
常见错误清单 · 自查六条
拿 ∥ A ∥ \lVert A\rVert ∥ A ∥ 当 ρ ( A ) \rho(A) ρ ( A ) 。 ρ ≤ ∥ A ∥ \rho\le\lVert A\rVert ρ ≤ ∥ A ∥ 是不等式 ,在非正规矩阵上两者可以差几个数量级(例 2.5.6 的 1 1 1 对 1000 1000 1000 )。用 ρ \rho ρ 判收敛时还要记住:ρ < 1 \rho\lt1 ρ < 1 只是渐近结论,有限步内误差可能先涨 ——稳妥的判据是 ∥ A ∥ < 1 \lVert A\rVert\lt1 ∥ A ∥ < 1 。(这正是课堂上「迭代收敛判定选 ∥ G ∥ < 1 \lVert G\rVert\lt1 ∥ G ∥ < 1 而非 ρ ( G ) < 1 \rho(G)\lt1 ρ ( G ) < 1 」的原因,见第 3 讲精编 定理后的预告与第 4 讲 的收敛定理;∥ A ∥ 2 = σ 1 \lVert A\rVert_2=\sigma_1 ∥ A ∥ 2 = σ 1 的公式见 2.6 。)
以为重特征值也能算准。 非正规矩阵的 k k k 重特征值遇上 O ( ε ) O(\varepsilon) O ( ε ) 扰动会移动 O ( ε 1 / k ) O(\varepsilon^{1/k}) O ( ε 1/ k ) :双精度下 k = 2 k=2 k = 2 只剩约 8 位有效数字,k = 4 k=4 k = 4 只剩约 4 位。这是问题的固有属性 ,换个算法没用(1.2.6 )。
把 Gerschgorin 圆盘当精确值。 它给出的是包含区间 ;只有对角占优时圆心才接近特征值。若圆盘相交,只能断定「这个连通区里有多少个特征值」,无法把它们一一配到具体位置上。
忘记 Gerschgorin 的列版本。 行、列两个版本的交集 能给出更紧的围栏(例 2.5.2 的 [ 8 , 12 ] [8,12] [ 8 , 12 ] 就是取交的收获)。对 A T A^{\mathrm T} A T 用行版本,即得列版本。
幂法中不做规范化。 A k x A^k\boldsymbol{x} A k x 的模长按 ∣ λ 1 ∣ k \lvert\lambda_1\rvert^k ∣ λ 1 ∣ k 增长或衰减,会溢出 / 下溢。每步都要归一化 (除以范数)。同样地,提精度要用 Rayleigh 商,而不是 ∥ A x ∥ / ∥ x ∥ \lVert A\boldsymbol{x}\rVert/\lVert\boldsymbol{x}\rVert ∥ A x ∥ / ∥ x ∥ 。
以为特征值都是实的。 实矩阵也可以有复特征值(旋转矩阵就是!);只有对称矩阵才保证特征值全实。算非对称矩阵时要用 np.linalg.eigvals(复)而不是 np.linalg.eigvalsh(假定对称)——错误地假定对称会静默给出错误结果 。
课程预告 · 本节之后
SVD 就是对称谱定理在长方形矩阵上的推广 :σ i ( A ) = λ i ( A T A ) \sigma_i(A)=\sqrt{\lambda_i(A^{\mathrm T}A)} σ i ( A ) = λ i ( A T A ) 在对称半正定 矩阵上算,因此总是良态的——2.6 ;
幂法 / 反幂法 / QR 方法 / Jacobi 方法求特征值 :原书第 3 章 3.1–3.3 节(课堂已于第 6 讲 开篇,并在第 7 讲 展开规范化幂法、反幂法与位移反幂法;本节的 Rayleigh 商正是幂法提取特征值的工具);
谱半径决定迭代收敛速度 :原书第 2 章 2.4 节与第 2 章 (ρ ( H ) < 1 \rho(H)\lt1 ρ ( H ) < 1 是 Jacobi / Gauss–Seidel / SOR 收敛的判据,课堂实录见第 5 讲 、第 6 讲 );
条件数与谱的关系 :c o n d 2 ( A ) = σ 1 / σ r \mathrm{cond}_2(A)=\sigma_1/\sigma_r cond 2 ( A ) = σ 1 / σ r (2.4 ),而 σ \sigma σ 来自对称矩阵的谱——两节在这里汇合。