第 02a 章 酉矩阵、QR 分解与 Schur 三角化
对应原书:Horn & Johnson《Matrix Analysis》第 2 版,第 2 章 Unitary Similarity and Unitary Equivalence 的 2.0–2.4 节(书 p.83–131,PDF p.103–151)。第 2 章的后半(正规矩阵、谱定理、SVD、CS 分解)见第 02b 章。
第 01 章用一般的可逆矩阵 \(S\) 做相似变换 \(S^{-1}AS\)。这一章把 \(S\) 限制为酉矩阵:换的是"标准正交基",长度和夹角都不变。代价是不能再要求把每个矩阵化成对角阵;收获是任何方阵都能用酉相似化成上三角阵(Schur 定理),而且计算上数值稳定。量化中几乎所有稠密矩阵算法——最小二乘的 QR 解法、协方差矩阵的特征分解、VAR 模型的特征根、Lyapunov 方程求平稳协方差——底层都是本章的酉变换。
学习目标
读完本章,你应当能够:
- 说清酉矩阵(实情形为正交矩阵)的几种等价刻画,理解"保长度 ⇔ 保内积 ⇔ 列标准正交",并会用 Givens 旋转与 Householder 反射构造酉矩阵。
- 掌握 QR 分解的存在性、唯一性和 Householder 构造;理解为什么用 QR 解最小二乘比正规方程稳定,以及 QR 与因子逐步正交化(Gram–Schmidt)的关系。
- 理解酉相似保持 Frobenius 范数,并能陈述和证明 Schur 三角化定理及其实版本(实 Schur 形),会用 Schur 不等式度量"偏离正规性"。
- 会用 Schur 定理推出 Cayley–Hamilton 定理、\(p(A)\) 的特征值、Sylvester 方程唯一可解条件、分块对角化和特征值的连续性。
- 能把这些工具用在量化场景:QR 回归、VAR(1) 稳定性分析、平稳协方差的 Lyapunov 方程、非正规系统的瞬时放大。
读前导读
这一章在解决什么问题
第 01 章的对角化有两个毛病:不是每个矩阵都能对角化,而且换基矩阵 \(S\) 可能很"歪",数值计算时会放大误差。这一章换一个思路:只允许用不变形的换基——旋转和反射。它们不改变任何向量的长度和任何两个向量的夹角,所以误差不会被放大。这种矩阵在实数里叫正交矩阵,在复数里叫酉矩阵。
你其实熟悉正交矩阵。PCA 的载荷矩阵(各主成分的权重向量排成的矩阵)就是正交的:每个主成分组合的权重平方和为 1,不同主成分两两正交。把原始收益换成主成分收益,就是乘了一个正交矩阵;总方差(各资产方差之和)在换前换后不变。
本章的两个主角:QR 分解把任意矩阵写成"正交矩阵 × 上三角矩阵",它就是因子逐步正交化(新因子对老因子回归取残差),也是软件解 OLS 的标准做法,比 CFA 课本里的 \((X^TX)^{-1}X^Ty\) 更稳定。Schur 定理说任何方阵都能用酉矩阵化成上三角阵,对角线上就是特征值。退而求其次(上三角而不是对角),换来"永远能做到、数值稳定"。由它可以推出一串有用结论,包括 Cayley–Hamilton 定理和 VAR 平稳协方差满足的 Lyapunov 方程。
需要先想起来的数学
1. 正交矩阵(实数版的酉矩阵)。 \(Q^TQ=I\),即每列长度为 1、两两垂直。于是 \(Q^{-1}=Q^T\),求逆只要转置。例:旋转 \(\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix}\),取 \(\theta=90^\circ\) 得 \(\begin{bmatrix}0&-1\\1&0\end{bmatrix}\),把 \([1,0]^T\) 转成 \([0,1]^T\),长度不变。复数版把 \(T\) 换成 \(*\):\(U^*U=I\)。见 第 00 册第 06 章 线性代数速成。
2. 复数的模与 \(e^{i\theta}\)。 \(|a+bi|=\sqrt{a^2+b^2}\);\(e^{i\theta}=\cos\theta+i\sin\theta\) 是单位圆上的点,模为 1。"酉矩阵的特征值在单位圆上"就是说它们都形如 \(e^{i\theta}\):实正交矩阵只能是 \(\pm1\) 或成对的 \(e^{\pm i\theta}\)(旋转角 \(\theta\))。
3. 最小二乘与正规方程。 OLS 求 \(\min_\beta\|y-X\beta\|_2^2\),令导数为零得正规方程 \(X^TX\beta=X^Ty\)。这是 CFA 里回归系数公式的矩阵形式。本章会说明为什么实际计算不直接用它。见 第 00 册第 05 章 多元微积分与优化。
4. 迹与 Frobenius 范数。 \(\operatorname{tr}A\) 是对角元之和;\(\|A\|_F=\sqrt{\sum_{i,j}|a_{ij}|^2}=\sqrt{\operatorname{tr}A^*A}\),即把矩阵当作一长串数字算长度。例:\(\left\|\begin{bmatrix}1&2\\3&4\end{bmatrix}\right\|_F=\sqrt{1+4+9+16}=\sqrt{30}\)。迹有"循环性":\(\operatorname{tr}(ABC)=\operatorname{tr}(BCA)\)。
5. 收敛子列与紧性(只用于 2.2.2 和 2.7.7)。 有界数列不一定收敛(如 \(1,-1,1,-1,\dots\)),但一定能挑出一个收敛的子列(只取奇数项)。"紧"的集合就是"每个序列都有收敛子列、而且极限还在集合里"。见 第 00 册第 01 章 函数极限与连续。记号 \(\perp\) 表示"垂直/正交",\(w^\perp\) 是与 \(w\) 垂直的全部向量;\(\operatorname{Re}\)、\(\operatorname{Im}\) 是取复数的实部、虚部。
怎么读这一章
核心必读:2.1(为什么要"酉")、2.2.1(酉矩阵的刻画)、2.3.2(Householder 反射的思想)、2.4.1–2.4.2(QR 分解与 Gram–Schmidt)、2.5.1(Frobenius 范数不变)、2.6.1–2.6.2(Schur 定理与偏离正规性)、2.7.2(Cayley–Hamilton)、2.7.3(Sylvester 与 Lyapunov 方程)。
第一遍可以只看结论:2.2.2 紧性、2.2.3、2.3.1 的旋转个数计算、2.4.3 的变体、2.5.2–2.5.3、2.6.3、2.6.4 的证明、2.7.4–2.7.9。实战 1(QR 回归)和实战 2(VAR 分析)与日常工作最接近,建议跑一遍。
2.1 为什么要"酉":从相似到酉相似
第 01 章的相似变换 \(A\mapsto S^{-1}AS\) 对应"换一组基"。如果新基是标准正交基,换基矩阵 \(U\) 满足 \(U^{-1}=U^*\),相似变换就成了
这叫酉相似(unitary similarity)。原书在 2.0 节给出两个理由说明它值得单独研究:
- 概念更简单:求逆变成了共轭转置,不需要解方程。
- 数值更稳定:酉矩阵不放大误差(\(\|Ux\|_2=\|x\|_2\)),所以由一串酉变换组成的算法,舍入误差不会被逐步放大。这正是 LAPACK 中 QR、特征值、SVD 算法都以酉变换为核心的原因。
两条主线贯穿第 2 章:
| 变换 | 形式 | 能达到的最简形 | 本册位置 |
|---|---|---|---|
| 酉相似 | \(A\mapsto U^*AU\)(方阵) | 上三角阵,对角元为特征值(Schur) | 本章 |
| 酉等价 | \(A\mapsto UAV\)(可为长方阵,\(U,V\) 独立) | 非负对角阵,对角元为奇异值(SVD) | 第 02b 章 |
另一个相关的变换是 \(A\mapsto S^*AS\)(\(S\) 非奇异),称为 *合同(*congruence),原书第 4 章研究;酉相似同时是相似和 *合同。
2.2 酉矩阵
2.2.1 定义与等价刻画
向量组 \(x_1,\dots,x_k\in\mathbf C^n\) 称为正交(orthogonal),若 \(x_i^*x_j=0\)(\(i\ne j\));再要求 \(x_i^*x_i=1\) 就是标准正交(orthonormal)。标准正交组一定线性无关:若 \(\sum\alpha_ix_i=0\),则
定义 2.1.3 \(U\in M_n\) 称为酉矩阵(unitary),若 \(U^*U=I\);实矩阵 \(U\) 满足 \(U^TU=I\) 时称为实正交矩阵(real orthogonal)。
定理 2.1.4(酉矩阵的等价刻画) 对 \(U\in M_n\),以下等价:
(a) \(U\) 酉;(b) \(U\) 非奇异且 \(U^{-1}=U^*\);(c) \(UU^*=I\);(d) \(U^*\) 酉;(e) \(U\) 的列标准正交;(f) \(U\) 的行标准正交;(g) 对所有 \(x\),\(\|Ux\|_2=\|x\|_2\)。
白话解释:先看实数情形。\(U^TU\) 的第 \((i,j)\) 元是第 \(i\) 列与第 \(j\) 列的内积,所以 \(U^TU=I\) 就是"每列长度为 1(对角元为 1)、不同列互相垂直(非对角元为 0)",即 (a)⇔(e)。旋转矩阵 \(\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix}\) 两列长度都是 \(\sqrt{\cos^2\theta+\sin^2\theta}=1\),内积 \(-\cos\theta\sin\theta+\sin\theta\cos\theta=0\),所以正交。 复数情形只是把内积换成带共轭的 \(x^*y\)。例:\(U=\frac1{\sqrt2}\begin{bmatrix}1&1\\i&-i\end{bmatrix}\)。第一列自身:\(\frac12(\bar1\cdot1+\bar i\cdot i)=\frac12(1+1)=1\);两列之间:\(\frac12(\bar1\cdot1+\bar i\cdot(-i))=\frac12(1-1)=0\)。所以 \(U^*U=I\),\(U\) 是酉矩阵。但若误用 \(U^TU\),第一列自身得 \(\frac12(1+i^2)=0\),结论全错。这就是复数里必须用 \(*\) 的原因。 金融含义:(g) 说酉/正交矩阵"保长度"。把资产收益 \(r\) 换成主成分收益 \(U^Tr\),\(\|U^Tr\|^2=\|r\|^2\),收益平方和不变;协方差的迹(总方差)也不变,只是重新分配到各个主成分上。
(a)–(f) 只是"方阵的左逆即右逆"加上矩阵乘法的列/行解读。值得细看的是 (g)⇒(a):保持长度为什么就能保持内积?令 \(A=U^*U\),对 \(x=z+w\) 展开 \(\|Ux\|^2=\|x\|^2\),得到 \(\operatorname{Re}z^*w=\operatorname{Re}z^*Aw\) 对所有 \(z,w\) 成立。取 \(z=e_p\)、\(w=ie_q\) 得 \(\operatorname{Im}a_{pq}=0\);取 \(z=e_p\)、\(w=e_q\) 得 \(a_{pq}=\delta_{pq}\)。所以 \(A=I\)。这就是极化恒等式的思想:内积可以由范数还原。
推导拆解:实数情形更容易看懂。展开 \(\|x+y\|^2=(x+y)^T(x+y)=\|x\|^2+2x^Ty+\|y\|^2\),于是
\[x^Ty=\tfrac12\big(\|x+y\|^2-\|x\|^2-\|y\|^2\big).\]右边只用到长度。若 \(U\) 保持所有向量的长度,就保持右边三项,从而保持左边的内积:\((Ux)^T(Uy)=x^Ty\),即 \(x^TU^TUy=x^Ty\) 对所有 \(x,y\) 成立;取 \(x=e_p\)、\(y=e_q\) 得 \(U^TU\) 的 \((p,q)\) 元等于 \(\delta_{pq}\),所以 \(U^TU=I\)。 这和组合方差公式是同一个恒等式:\(\operatorname{Var}(a+b)=\operatorname{Var}a+\operatorname{Var}b+2\operatorname{Cov}(a,b)\),所以协方差可以由三个方差还原。复数情形 \(z^*w\) 有实部和虚部两部分,所以原书要分别取 \(w=e_q\) 和 \(w=ie_q\) 两次,才能把 \(a_{pq}\) 的实部和虚部都确定下来。
满足 \(\|Tx\|_2=\|x\|_2\) 的线性映射叫欧氏等距(Euclidean isometry)。对方阵,等距与酉是一回事。\(2\times2\) 实正交矩阵只有两类:旋转 \(\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix}\) 与反射 \(\begin{bmatrix}\cos\theta&\sin\theta\\\sin\theta&-\cos\theta\end{bmatrix}\)。
几个直接性质(原书 2.1.P1–P4、P9):\(|\det U|=1\);酉矩阵的特征值都在单位圆上;\(x\) 是 \(U\) 的右特征向量当且仅当它是同一特征值的左特征向量;酉变换保持正交性;对角酉矩阵恰为 \(\operatorname{diag}(e^{i\theta_1},\dots,e^{i\theta_n})\),对角实正交矩阵的对角元为 \(\pm1\)。
2.2.2 酉群与紧性
酉矩阵之积仍是酉矩阵,逆也是,所以 \(n\times n\) 酉矩阵构成一个群,称为酉群(unitary group);实正交矩阵构成实正交群。\(n!\) 个置换矩阵是实正交群的子群(2.1.P5)。
更重要的是紧性:\(U^*U=I\) 使每列长度为 1,于是 \(|u_{ij}|\le1\),酉群有界;又若 \(U_k\to U\),则 \(U^*U=\lim U_k^*U_k=I\),酉群是闭集。有限维中有界闭集紧,得到:
引理 2.1.8(选择原理,selection principle) 任一酉矩阵序列都有逐元收敛的子序列,且极限仍是酉矩阵。
这条看似抽象的引理,是本章证明"特征值连续依赖于矩阵元素"(2.4.9)和第 02b 章"奇异值连续"的关键:先在酉群里取收敛子列,再取极限。注意极限依赖子序列:\(\begin{bmatrix}0&1\\1&0\end{bmatrix}^k\) 的偶数项收敛到 \(I\),奇数项收敛到它自己。
白话解释:这段的意思是"酉矩阵不会跑到无穷远,也不会在极限处突然失去酉性"。有界性:酉矩阵每列长度为 1,所以每个元素的模都不超过 1,所有酉矩阵都住在一个有限的"盒子"里。闭性:\(U_k^*U_k=I\) 对每个 \(k\) 成立,取极限两边仍相等。有界加闭,就保证任何酉矩阵序列都能挑出收敛的子列,极限仍是酉矩阵——这和"有界数列必有收敛子列"是同一件事的矩阵版。 为什么证明里需要它?比如要证"特征值随矩阵连续变化",做法是对一列矩阵 \(A_k\) 各取一个 Schur 分解 \(A_k=U_kT_kU_k^*\)。\(U_k\) 可能不收敛(Schur 形不唯一,\(U_k\) 会乱跳),但选择原理允许先挑一个收敛子列,再在这个子列上取极限。一般可逆矩阵没有这个性质:\(\operatorname{diag}(k,1/k)\) 都可逆,却跑向无穷。
与之对比,复正交矩阵(\(A^TA=I\),\(A\) 可以是复的)构成的群无界:取 \(S=\begin{bmatrix}0&1\\-1&0\end{bmatrix}\),\(A(t)=(\cosh t)I+(i\sinh t)S\) 对所有实 \(t\) 都复正交,但元素随 \(t\) 无界增长,且只在 \(t=0\) 时是酉矩阵(2.1.P8)。复正交矩阵既酉又复正交 ⇔ 它是实正交矩阵。
2.2.3 酉矩阵的分块结构
引理 2.1.10 设 \(U=\begin{bmatrix}U_{11}&U_{12}\\U_{21}&U_{22}\end{bmatrix}\) 酉,\(U_{11}\in M_k\),则
特别地 \(U_{12}=0\iff U_{21}=0\),此时 \(U_{11},U_{22}\) 都是酉矩阵。推论:上三角的酉矩阵一定是对角阵。证明用第 00 章的互补零化度定律作用于 \(U^{-1}=U^*\)。这条引理在 Schur 形唯一性、正规矩阵和 CS 分解(第 02b 章)中反复出现。
2.3 两把基本工具:Givens 旋转与 Householder 反射
所有实用的酉变换算法都由两类"积木"拼成。
2.3.1 平面旋转(Givens 旋转)
例 2.1.11 \(U(\theta;i,j)\)(\(i<j\))由单位阵把 \((i,i),(j,j)\) 元换成 \(\cos\theta\)、\((i,j)\) 元换成 \(-\sin\theta\)、\((j,i)\) 元换成 \(\sin\theta\) 得到。它在 \((i,j)\) 坐标平面内旋转角 \(\theta\),其余坐标不动。左乘只改第 \(i,j\) 行,右乘只改第 \(i,j\) 列,\(U(\theta;i,j)^{-1}=U(-\theta;i,j)\)。
用途是逐个消元:选 \(\theta\) 使 \(\cos\theta=x_{n-1}/\sqrt{x_{n-1}^2+x_n^2}\)、\(\sin\theta=-x_n/\sqrt{\cdot}\),旋转后向量第 \(n\) 个分量变成 0 而 2-范数不变(2.1.P27)。重复下去,任意实 \(A\in M_{n,m}\) 可经 \(m(n-\frac{m+1}2)\) 个平面旋转化成上三角(2.1.P28)。稀疏矩阵、滚动回归中"加一行删一行"的 QR 更新都用 Givens 旋转,因为它只碰两行。
2.3.2 Householder 反射
例 2.1.12 对非零 \(w\),令
性质:\(U_w\) 既酉又 Hermite,\(U_w^{-1}=U_w\);\(U_ww=-w\),而 \(x\perp w\) 时 \(U_wx=x\)。所以它是关于超平面 \(w^\perp\) 的镜面反射。实情形 \(\det U_w=-1\),特征值为 \(-1,1,\dots,1\);因而实 Householder 矩阵永远不是"真旋转"(行列式为 \(+1\) 的实正交阵)。
关键用法:把一个向量反射到另一个同长度的向量。在 \(\mathbf R^n\) 中,若 \(\|x\|_2=\|y\|_2\)、\(x\ne y\),取 \(w=x-y\) 就有 \(U_wx=y\)。复情形需要多乘一个模为 1 的数:
推导拆解:先看为什么 \(U_w\) 是镜面反射。\(U_wv=v-2\dfrac{w^*v}{w^*w}w\),其中 \(\dfrac{w^*v}{w^*w}w\) 是 \(v\) 在 \(w\) 方向上的投影(就像回归里 \(v\) 对 \(w\) 的拟合值)。从 \(v\) 减去一次投影落到镜面 \(w^\perp\) 上,减去两次就到了镜面另一侧,正是反射。所以 \(w\) 本身被翻成 \(-w\),与 \(w\) 垂直的向量不动。 数值例子:\(x=[3,4]^T\),目标 \(y=[5,0]^T\)(同为长度 5)。取 \(w=x-y=[-2,4]^T\),\(w^Tw=20\),\(w^Tx=-6+16=10\)。于是 \(U_wx=x-2\cdot\dfrac{10}{20}w=x-w=y=[5,0]^T\)。一次反射就把第二个分量消成了 0,长度保持 5。 这正是 QR 分解每一步做的事:把一列"压"到第一个坐标轴上,下方全变成零。和 Givens 旋转一次消一个元素不同,Householder 一次消掉一整列的下半部分。
定理 2.1.13 设 \(\|x\|_2=\|y\|_2>0\)。若 \(y=e^{i\theta}x\),取 \(U(y,x)=e^{i\theta}I\);否则取 \(\phi\) 使 \(x^*y=e^{i\phi}|x^*y|\),令 \(w=e^{i\phi}x-y\)、\(U(y,x)=e^{i\phi}U_w\)。则 \(U(y,x)\) 酉、\(U(y,x)x=y\),且 \(z\perp x\Rightarrow U(y,x)z\perp y\)。\(x,y\) 实时 \(U(y,x)\) 是实 Householder 矩阵 \(U_{x-y}\)。
两个立即可用的推论:\(U(\|x\|_2e_1,x)\) 把 \(x\) 变成 \(\|x\|_2e_1\)(QR 分解的一步);\(U(y,e_1)\) 的第一列是 \(y\)——这给出"把一个单位向量扩充为酉矩阵"的显式构造,Schur 定理的证明需要它。
两个 Householder 反射的乘积是一个旋转。原书 2.1.P16 的 Palais 矩阵 \(P_{x,y}=U_wU_x\)(\(w=x+y\))就是在平面 \(\operatorname{span}\{x,y\}\) 内把单位向量 \(x\) 转到 \(y\)、在正交补上为恒等的真旋转。
2.4 QR 分解
2.4.1 定理与构造
定理 2.1.14(QR 分解) 设 \(A\in M_{n,m}\)。
(a) 若 \(n\ge m\),存在列标准正交的 \(Q\in M_{n,m}\) 和对角元非负的上三角 \(R\in M_m\),使 \(A=QR\)("瘦" QR,thin QR); (b) 若 \(\operatorname{rank}A=m\)(列满秩),则 (a) 中 \(Q,R\) 唯一,且 \(R\) 对角元全为正; (c) 若 \(m=n\),\(Q\) 是酉矩阵; (d) 总存在酉 \(Q\in M_n\) 和对角元非负的上三角 \(R\in M_{n,m}\) 使 \(A=QR\)("全" QR); (e) \(A\) 实时以上因子都可取实。
构造(Householder 算法) 设 \(a_1\) 是 \(A\) 的第一列,\(r_1=\|a_1\|_2\)。用定理 2.1.13 取酉 \(U_1\) 使 \(U_1a_1=r_1e_1\),则
对 \(A_2\) 的第一列重复,用 \(U_2=I_1\oplus V_2\) 不破坏已经做好的第一行第一列。\(m\) 步后 \(U_m\cdots U_1A=\begin{bmatrix}R\\0\end{bmatrix}\)。令 \(U^*=(U_m\cdots U_1)^*=[Q\ Q_2]\) 即得 \(A=QR\)。
唯一性的论证很短,值得记住:若 \(A=QR=\tilde Q\tilde R\) 且列满秩,则 \(A^*A=R^*R=\tilde R^*\tilde R\),于是 \(\tilde R^{-*}R^*=\tilde RR^{-1}\)。左边下三角、右边上三角,只能是对角阵 \(D\),且对角元为正;再由 \(D=D^{-1}\) 得 \(D=I\)。
2.4.2 QR 就是 Gram–Schmidt
第 00 章 0.5 节的 Gram–Schmidt 过程把 \(x_1,\dots,x_m\) 依次正交化,矩阵形式正是 \(X=QR\)。原书 2.1.P18 点出了每个量的含义:
- \(q_1,\dots,q_k\) 是 \(\operatorname{span}\{a_1,\dots,a_k\}\) 的标准正交基;
- \(r_{kk}\) 等于 \(a_k\) 到 \(\operatorname{span}\{a_1,\dots,a_{k-1}\}\) 的欧氏距离,也就是 \(a_k\) 对前 \(k-1\) 列做最小二乘回归的残差范数。
用量化的话说:把因子按某个顺序排好做 QR,第 \(k\) 个 \(q_k\) 就是"第 \(k\) 个因子剔除前面所有因子之后的纯净部分"(归一化的回归残差),\(r_{kk}\) 衡量它还剩多少独立信息。\(r_{kk}\) 很小意味着第 \(k\) 个因子几乎可以由前面的因子线性表示——这是共线性诊断。实战 1 会验证这一点。
推导拆解:手算一个 \(2\times2\) 的 QR。\(A=\begin{bmatrix}3&1\\4&2\end{bmatrix}\),两列 \(a_1=[3,4]^T\),\(a_2=[1,2]^T\)。 第一步,单位化第一列:\(r_{11}=\|a_1\|=5\),\(q_1=[0.6,\,0.8]^T\)。 第二步,把 \(a_2\) 对 \(q_1\) "回归":系数 \(r_{12}=q_1^Ta_2=0.6+1.6=2.2\)。残差 \(a_2-2.2q_1=[1-1.32,\ 2-1.76]^T=[-0.32,\,0.24]^T\)。 第三步,单位化残差:\(r_{22}=\sqrt{0.32^2+0.24^2}=0.4\),\(q_2=[-0.8,\,0.6]^T\)。 结果 \(Q=\begin{bmatrix}0.6&-0.8\\0.8&0.6\end{bmatrix}\),\(R=\begin{bmatrix}5&2.2\\0&0.4\end{bmatrix}\),可验证 \(QR=A\)。\(R\) 的每一列记录了"原第 \(k\) 列 = 前 \(k\) 个正交方向的什么组合",所以 \(R\) 是上三角:第 1 列只用到 \(q_1\),第 2 列用到 \(q_1,q_2\)。\(r_{22}=0.4\) 相对 \(\|a_2\|\approx2.24\) 较小,说明 \(a_2\) 大部分已被 \(a_1\) 解释。顺带 \(|\det A|=r_{11}r_{22}=2\)。
推导拆解:为什么用 QR 解 OLS?设 \(X=QR\)(瘦 QR,\(Q^TQ=I\),\(R\) 可逆)。把它代入正规方程 \(X^TX\beta=X^Ty\): 左边 \(X^TX=R^TQ^TQR=R^TR\)(用了 \(Q^TQ=I\)),右边 \(X^Ty=R^TQ^Ty\)。两边左乘 \(R^{-T}\),得 \(R\beta=Q^Ty\)。 结论一样,但计算路线不同:\(R\) 是上三角,从最后一个方程往上逐个回代即可;整个过程从不形成 \(X^TX\)。\(X^TX\) 的条件数是 \(X\) 的平方(实战 1 前言),形成它就等于先把误差放大了一遍。这就像算组合收益时先把权重四舍五入到两位小数再相乘,误差在第一步就已经引入,后面算得再准也补不回来。
教科书式的 Gram–Schmidt 数值上不稳定(正交性会逐步丢失);Householder QR 不存在这个问题,numpy.linalg.qr 调用的就是 LAPACK 的 Householder 实现。
2.4.3 QR 的几个推论
Hadamard 不等式(2.1.P23) 由 \(A=QR\),\(|\det A|=\det R=r_{11}\cdots r_{nn}\),而 \(\|a_i\|_2=\|r_i\|_2\ge r_{ii}\),所以
等号成立当且仅当某列为零或各列两两正交。几何含义:平行多面体的体积不超过各棱长之积。用于协方差矩阵 \(\Sigma\) 时(对其 Cholesky 因子用上式),得到 \(\det\Sigma\le\prod\sigma_{ii}\):相关性只会让"广义方差"变小。
Cholesky 分解 任意 \(B=A^*A\) 都能写成 \(B=LL^*\),\(L\) 下三角、对角元非负:由 \(A=QR\) 得 \(B=R^*R\),取 \(L=R^*\)。\(A\) 非奇异时唯一。每个半正定矩阵都有这种分解(原书 7.2.9)。蒙特卡洛模拟相关正态收益时,\(x=Lz\)(\(z\sim N(0,I)\))的协方差恰为 \(LL^T=\Sigma\)。
QR 的变体 对 \(A^*\) 做 QR 再取共轭转置得 LQ 分解 \(A=LQ\);借助反序矩阵 \(K\)(\(K^2=I\),\(KRK\) 把上三角变下三角)还有 QL 和 RQ 分解。它们在原理上都一样,只是三角因子的位置不同。
定理 2.1.18 若 \(X,Y\in M_{n,k}\) 都列标准正交,则存在酉 \(U\) 使 \(Y=UX\)。做法:把两者分别扩充成酉矩阵 \(V=[X\ X_2]\)、\(W=[Y\ Y_2]\),取 \(U=WV^*\)。因子模型里,载荷矩阵右乘任意正交矩阵(因子旋转,如 varimax)不改变模型拟合,背后就是这个事实。原书 2.1.P22 进一步说:两个列标准正交矩阵的列空间相同,当且仅当 \(X=YU\)(\(U\) 酉)。
2.5 酉相似
2.5.1 Frobenius 范数不变
定义 2.2.1 若存在酉 \(U\) 使 \(A=UBU^*\),称 \(A\) 与 \(B\) 酉相似;\(U\) 可取实正交时称实正交相似。\(A\) 酉相似于对角阵时称可酉对角化。
定理 2.2.2 若 \(A=UBV\)(\(U,V\) 酉,可以不同),则
即 Frobenius 范数 \(\|A\|_F=(\operatorname{tr}A^*A)^{1/2}\) 不变。证明一行:\(\operatorname{tr}A^*A=\operatorname{tr}(V^*B^*U^*UBV)=\operatorname{tr}(B^*BVV^*)=\operatorname{tr}B^*B\)。
这给出判断"不酉相似"的最简单工具。例:\(\begin{bmatrix}3&1\\-2&0\end{bmatrix}\) 与 \(\begin{bmatrix}1&1\\0&2\end{bmatrix}\) 特征值都是 1、2,所以相似;但元素平方和分别是 14 和 6,所以不酉相似。酉相似是比相似更细的等价关系。
推导拆解:那一行证明用了三件事。(1) \((UBV)^*=V^*B^*U^*\)(反序律)。(2) 中间的 \(U^*U=I\) 直接消掉。(3) 迹的循环性 \(\operatorname{tr}(XY)=\operatorname{tr}(YX)\):把 \(V^*\) 从最前面挪到最后面,\(\operatorname{tr}(V^*\cdot B^*BV)=\operatorname{tr}(B^*BV\cdot V^*)\),再用 \(VV^*=I\)。 白话:Frobenius 范数是"矩阵所有元素的平方和",相当于把矩阵看成一个长向量后的长度。旋转坐标系不改变长度,所以也不改变它。一般相似 \(S^{-1}AS\) 会把矩阵"拉歪",元素平方和可以任意变,所以 Frobenius 范数能区分相似和酉相似。
2.5.2 两个可达形式
例 2.2.3(对角元全相等) 任意 \(A\) 都酉相似于对角元全等于 \(\frac1n\operatorname{tr}A\) 的矩阵(\(A\) 实时可用实正交相似)。\(2\times2\) 情形先把迹平移为 0,再构造单位向量 \(u\) 使 \(u^*Au=0\);一般情形每次挑对角元相差最大的一对 \((i,j)\),用只作用于第 \(i,j\) 行列的 \(2\times2\) 酉相似把两者都换成平均值,再用紧性论证收尾。推论(2.2.P9):\(\operatorname{tr}A=0\) ⇔ \(A\) 是两个幂零矩阵之和。
例 2.2.4(上 Hessenberg 形) 任意 \(A\) 酉相似于次对角元非负的上 Hessenberg 矩阵(\(a_{ij}=0\) 当 \(i>j+1\))。做法是对第一列下方的 \(n-1\) 个元素用 Householder 反射 \(V_1=I_1\oplus U_1\) 压成一个元素,再对右下子块重复;注意右乘 \(V_1^*\) 不会破坏第一列。若 \(A\) Hermite,结果是三对角 Hermite 矩阵。这是对称特征值算法(以及 numpy.linalg.eigh)的第一步:先 \(O(n^3)\) 化成三对角,再在三对角矩阵上迭代。
2.5.3 Specht 定理(了解)
既然 Frobenius 范数(即 \(\operatorname{tr}A^*A\))是酉不变量,更一般地:对两个非交换变量的任意"词" \(W(s,t)=s^{m_1}t^{n_1}\cdots s^{m_k}t^{n_k}\),\(\operatorname{tr}W(A,A^*)\) 都是酉不变量。
定理 2.2.6(Specht) \(A,B\) 酉相似 ⇔ 对每个词 \(W\),\(\operatorname{tr}W(A,A^*)=\operatorname{tr}W(B,B^*)\)。
原书 2.2.8 给出有限化:\(n=2\) 只需检查 \(s,s^2,st\) 三个词(即 \(\operatorname{tr}A\)、\(\operatorname{tr}A^2\)、\(\operatorname{tr}AA^*\)),\(n=3\) 需 7 个,\(n=4\) 需 20 个。推论:每个 \(2\times2\) 复矩阵都酉相似于其转置(2.2.P3),而 \(3\times3\) 不一定(2.2.P4)。这部分与量化没有直接关系,知道它存在即可。
2.5.4 Jacobi 方法与循环矩阵(两个可计算的例子)
Jacobi 方法(2.2.P1) 对实对称 \(A\),取绝对值最大的非对角元 \(a_{ij}\),用平面旋转 \(U(\theta;i,j)\) 做实正交相似把它消成 0。由定理 2.2.2,非对角元平方和恰好减少 \(2a_{ij}^2\),而最大元的平方不小于平均值,所以
非对角"质量"几何收敛到 0,极限对角阵给出特征值,旋转之积给出特征向量。实战 3 实现了它。
循环矩阵与 Fourier 矩阵(2.2.P10) 设 \(\omega=e^{2\pi i/n}\),\(F_n=n^{-1/2}[\omega^{(i-1)(j-1)}]\) 是对称酉矩阵(离散 Fourier 变换)。首行为 \([a_1,\dots,a_n]\) 的循环矩阵(每行是上一行右移一位)满足 \(A=F_n\Lambda F_n^*\),
即所有循环矩阵被同一个酉矩阵 \(F_n\) 对角化,特征值就是首行的离散 Fourier 变换。时间序列里的循环卷积、平稳过程协方差矩阵的近似对角化(Toeplitz 矩阵在大 \(n\) 时近似循环矩阵)、FFT 快速计算都以此为基础。
2.6 Schur 三角化
2.6.1 定理与证明
原书称下面这个定理是"初等矩阵理论中最根本有用的事实"。
定理 2.3.1(Schur 三角化) 设 \(A\in M_n\) 的特征值按任意指定的顺序为 \(\lambda_1,\dots,\lambda_n\),\(x\) 是对应 \(\lambda_1\) 的单位特征向量。
(a) 存在酉 \(U=[x\ u_2\ \cdots\ u_n]\) 使 \(U^*AU=T\) 为上三角,且 \(t_{ii}=\lambda_i\); (b) 若 \(A\) 实且特征值全实,则 \(U\) 可取实正交矩阵。
证明(逐次收缩,deflation) 取以 \(x\) 为第一列的酉矩阵 \(U_1\)(例如 2.1.13 的 \(U(x,e_1)\))。由于 \(Ax=\lambda_1x\) 且 \(u_i^*x=0\),
\(A_1\) 的特征值是 \(\lambda_2,\dots,\lambda_n\)。对 \(A_1\) 取 \(\lambda_2\) 的单位特征向量做同样的事,用 \(V_2=[1]\oplus U_2\) 嵌回去;\(n-1\) 步后得到上三角阵,\(U=U_1V_2\cdots V_{n-1}\)。∎
推导拆解:为什么 \(U_1^*AU_1\) 的第一列是 \([\lambda_1,0,\dots,0]^T\)?矩阵乘积的第一列等于"矩阵乘以第一列"。\(U_1\) 的第一列是 \(x\),所以 \(U_1^*AU_1\) 的第一列是 \(U_1^*Ax=U_1^*(\lambda_1x)=\lambda_1U_1^*x\)。而 \(U_1^*x\) 的第 \(i\) 个分量是 \(U_1\) 第 \(i\) 列与 \(x\) 的内积:第一列就是 \(x\) 本身,内积为 1;其余列与 \(x\) 正交,内积为 0。所以 \(U_1^*x=e_1\),第一列是 \(\lambda_1e_1\)。 为什么 \(A_1\) 的特征值是剩下那些?分块上三角阵的特征值是对角块特征值的并(第 01 章 1.3.2),总共还是 \(\lambda_1,\dots,\lambda_n\),去掉左上角的 \(\lambda_1\) 就是 \(A_1\) 的。 "\(\oplus\)" 是直和:\([1]\oplus U_2=\begin{bmatrix}1&0\\0&U_2\end{bmatrix}\),它不动已做好的第一行第一列,只处理右下角。整个证明就是"剥洋葱":每次剥出一个特征值放到对角线上。
白话解释:Schur 定理和对角化比,差在哪、好在哪?对角化要求换到的基恰好是特征向量,而特征向量可能不够 \(n\) 个,也可能彼此夹角很小("歪"),所以不总能做到、做到了也可能数值上不稳。Schur 定理只要求换到一组标准正交基,使得矩阵变成上三角,这永远能做到。 上三角的含义可以用 VAR 理解。在新坐标 \(z=U^*x\) 下,\(z_{t+1}=Tz_t+\cdots\):最后一个坐标 \(z_n\) 只受自己影响(\(z_{n,t+1}=\lambda_nz_{n,t}\)),倒数第二个只受自己和 \(z_n\) 影响,依此类推。系统被拆成一条"单向传导链",每个环节的自身衰减率就是对角线上的特征值。所以看平稳性只要看对角元。
两个常被忽略的细节:
- 特征值顺序可以任意指定,这在后面"把相同的特征值排在一起"(通向 Jordan 形)时很关键。
- Schur 形不唯一。原书例 2.3.2 中
\[T_1=\begin{bmatrix}1&1&4\\0&2&2\\0&0&3\end{bmatrix},\quad T_2=\begin{bmatrix}2&-1&3\sqrt2\\0&1&\sqrt2\\0&0&3\end{bmatrix}\]经 \(U=\frac1{\sqrt2}\begin{bmatrix}1&1&0\\1&-1&0\\0&0&\sqrt2\end{bmatrix}\) 酉相似,上三角部分完全不同。对转置用 Schur 定理还能得到下三角形式。
2.6.2 Schur 不等式与"偏离正规性"
设 \(A=UTU^*\)。由定理 2.2.2,\(\|A\|_F^2=\|T\|_F^2=\sum|\lambda_i|^2+\sum_{i<j}|t_{ij}|^2\),于是
等号当且仅当 \(T\) 是对角阵,即 \(A\) 可酉对角化(第 02b 章将证明这等价于 \(A\) 正规)。差值
称为偏离正规性(defect from normality),它与选哪个 Schur 形无关。\(\Delta(A)\) 大,说明 \(A\) 的特征向量彼此"挤在一起"、远离正交。
推导拆解:用 2.5.1 的两个例子算一下。\(A=\begin{bmatrix}3&1\\-2&0\end{bmatrix}\):元素平方和 \(9+1+4+0=14\),特征值 1、2 的模平方和 \(1+4=5\),所以 \(\Delta(A)=9\)。它的某个 Schur 形必然是 \(\begin{bmatrix}1&t\\0&2\end{bmatrix}\)(或对角元顺序相反),且 \(|t|^2=9\),即上三角部分的"质量"是 9。\(B=\begin{bmatrix}1&1\\0&2\end{bmatrix}\) 本身就是 Schur 形,\(\Delta(B)=1\)。对称矩阵 \(\begin{bmatrix}2&1\\1&2\end{bmatrix}\):平方和 \(4+1+1+4=10\),特征值 3、1 的平方和 10,\(\Delta=0\)。 白话:\(\|A\|_F^2\) 是矩阵的总"能量",其中 \(\sum|\lambda_i|^2\) 是装在特征值里的部分,\(\Delta(A)\) 是装在非对角耦合里的部分。对称矩阵(协方差)的能量全在特征值里,所以 PCA 能把风险完全分解到主成分上;非对称的 VAR 系数矩阵有一部分能量在耦合里,这部分会造成下文说的"先放大后衰减"。
一个反例:\(2\times2\) 时,特征值相同且元素平方和相同就足以推出酉相似(Specht 的 \(n=2\) 情形);\(3\times3\) 不行。原书 (2.3.2b):
特征值都是 1、2、3,平方和都是 39,但不酉相似(它们相似,因为特征值互异)。
量化含义:对 VAR(1) 系数矩阵 \(A\),若 \(A\) 正规,则 \(\|A^k\|_2=\rho(A)^k\) 单调衰减;若 \(A\) 强烈非正规,即使 \(\rho(A)<1\),\(\|A^k\|\) 也可能先放大很多倍再衰减——冲击在系统里先被放大,之后才消散。实战 2 给出一个 \(\rho=0.9\) 但瞬时放大 13 倍的例子。
2.6.3 交换族的同时三角化
定理 2.3.3 若 \(\mathcal F\subseteq M_n\) 是交换族(两两交换),则存在同一个酉 \(U\) 使所有 \(U^*AU\)(\(A\in\mathcal F\))都是上三角。
证明沿用 2.3.1,每一步取族中矩阵的公共特征向量(第 01 章:交换族有公共特征向量)。代价是不能再指定对角元的顺序。交换性是充分而非必要的(2.3.P4 给出不交换但可同时三角化的例子);必要条件是 \(AB-BA\) 的特征值全为 0(2.3.P6)。精确刻画见 2.4 节的 McCoy 定理。
2.6.4 实 Schur 形
实矩阵若有非实特征值,就不可能经实相似化为上三角(对角元会是复数)。退一步,可以化成上拟三角(upper quasitriangular):对角线上是 \(1\times1\) 或 \(2\times2\) 的块。
定理 2.3.4(实 Schur 形) 设 \(A\in M_n(\mathbf R)\)。
(a) 存在实非奇异 \(S\) 使 \(S^{-1}AS\) 为实上拟三角阵,\(1\times1\) 块是实特征值,每个 \(2\times2\) 块形如
(b) 由 (a) 推出:对 \(S\) 做实 QR 分解 \(S=QR\),则 \(Q^TAQ=R(S^{-1}AS)R^{-1}\) 仍是上拟三角,\(2\times2\) 块被相似变形了。scipy.linalg.schur(A, output='real') 计算的就是 (b)。
\(\begin{bmatrix}a&b\\-b&a\end{bmatrix}=r\begin{bmatrix}\cos\theta&\sin\theta\\-\sin\theta&\cos\theta\end{bmatrix}\)(\(r=\sqrt{a^2+b^2}\))是"伸缩 \(r\) 倍再旋转"。在离散动力系统 \(x_{t+1}=Ax_t\) 中,这样的块对应以 \(2\pi/\theta\) 为周期、以 \(r\) 为每期衰减率的振荡模态——商业周期、AR(2) 中的"伪周期"都来自这里。
原书还给出交换实矩阵族的同时实拟三角化(2.3.6),以及满足 \(A\bar A=\bar AA\) 的复矩阵可实正交相似于复上拟三角形(2.3.7)。这些属于理论补充。
2.7 Schur 定理的推论
这一节是 Schur 定理的"丰收"。每一条都能用"先化成上三角,再看对角线"的套路证明。
2.7.1 迹、行列式与多项式的特征值
由 \(A=UTU^*\),\(\operatorname{tr}A=\sum\lambda_i\)、\(\det A=\prod\lambda_i\) 立刻可见。更进一步,\(p(A)=Up(T)U^*\),而 \(p(T)\) 上三角、对角元为 \(p(\lambda_i)\),所以:
\(p(A)\) 的特征值(计重数)恰为 \(p(\lambda_1),\dots,p(\lambda_n)\)。特别地 \(\operatorname{tr}A^k=\sum_i\lambda_i^k\)。
第 01 章只得到了不计重数的版本,Schur 定理补上了重数。
严格上三角阵 \(T\) 的 \(T^p\) 主对角线和前 \(p-1\) 条超对角线都是 0,所以 \(T^n=0\)。于是以下等价:\(A\) 幂零;\(A^n=0\);\(A\) 的特征值全为 0。
2.7.2 Cayley–Hamilton 定理
定理 2.4.3.2(Cayley–Hamilton) 设 \(p_A(t)=\det(tI-A)\),则 \(p_A(A)=0\)。
证明 \(p_A(A)=U\big[(T-\lambda_1I)(T-\lambda_2I)\cdots(T-\lambda_nI)\big]U^*\)。\(T-\lambda_1I\) 的 \((1,1)\) 元为 0;乘上 \((2,2)\) 元为 0 的 \(T-\lambda_2I\),乘积左上 \(2\times2\) 块为 0;归纳地,乘到第 \(k\) 个因子时左上 \(k\times k\) 块为 0(引理 2.4.3.1)。乘完 \(n\) 个因子,整个矩阵为 0。∎
推导拆解:\(2\times2\) 时把乘积写出来。\(T=\begin{bmatrix}\lambda_1&t\\0&\lambda_2\end{bmatrix}\),
\[(T-\lambda_1I)(T-\lambda_2I)=\begin{bmatrix}0&t\\0&\lambda_2-\lambda_1\end{bmatrix}\begin{bmatrix}\lambda_1-\lambda_2&t\\0&0\end{bmatrix}=\begin{bmatrix}0\cdot(\lambda_1-\lambda_2)+t\cdot0&0\cdot t+t\cdot0\\0&0\end{bmatrix}=0.\]第一个因子的第一列为零,第二个因子的第二行为零,相乘时每一项都会碰到一个零。第一步 \(p_A(A)=Up_A(T)U^*\) 用的是 \(A^k=(UTU^*)^k=UT^kU^*\)(中间的 \(U^*U\) 抵消)。 再用例 2.4.3.3 核对:\(A=\begin{bmatrix}3&1\\-2&0\end{bmatrix}\),\(A^2=\begin{bmatrix}7&3\\-6&-2\end{bmatrix}\),\(3A-2I=\begin{bmatrix}7&3\\-6&-2\end{bmatrix}\),确实 \(A^2-3A+2I=0\)。
原书特意给出两个错误论证作为练习:(1) "\(p_A(A)\) 的特征值都是 \(p_A(\lambda_i)=0\),所以 \(p_A(A)=0\)"——错,特征值全为 0 只说明幂零,\(\begin{bmatrix}0&1\\0&0\end{bmatrix}\ne0\);(2) "\(p_A(A)=\det(AI-A)=\det0=0\)"——错,左边是矩阵,右边是数,\(p_A(A)\) 是先算出标量多项式再把矩阵代进去。
用途:高次幂与逆都降为 \(n-1\) 次多项式。例 2.4.3.3:\(A=\begin{bmatrix}3&1\\-2&0\end{bmatrix}\),\(p_A(t)=t^2-3t+2\),所以
一般地,若 \(p_A(t)=t^n+a_{n-1}t^{n-1}+\cdots+a_1t+a_0\) 且 \(a_0\ne0\),则
更系统的做法(2.4.P30):用多项式除法 \(t^k=h(t)p_A(t)+r(t)\),\(\deg r<n\),则 \(A^k=r(A)\)。实战 2 用它计算了 \(A^{25}\)。仿射期限结构模型中的 \(e^{At}\)、信用评级迁移矩阵的多期转移,都能这样化为低次多项式(第 03a 章还会从 Jordan 形和最小多项式的角度再看一次)。
原书 2.4.P3 给出不依赖 Schur 定理的证明(只用加法和乘法,对任意交换环成立),并顺带得到 \(\operatorname{adj}A\) 是 \(A\) 的多项式:
注意 Cayley–Hamilton 给出的零化多项式不一定次数最低。例 2.4.3.5:\(A=\begin{bmatrix}1&0&0\\0&1&1\\0&0&1\end{bmatrix}\) 的特征多项式是 \((t-1)^3\),但 \((A-I)^2=0\) 已经成立。次数最低的零化多项式(最小多项式)是第 03a 章的主题。
Newton 恒等式(2.4.P9) 设幂和 \(\mu_k=\sum\lambda_i^k=\operatorname{tr}A^k\),则
所以前 \(n\) 个幂迹 \(\operatorname{tr}A,\dots,\operatorname{tr}A^n\) 唯一决定特征多项式(2.4.P10);\(A\) 幂零 ⇔ \(\operatorname{tr}A^k=0\)(\(k=1,\dots,n\))。
2.7.3 Sylvester 方程
形如
的方程称为 Sylvester 方程。\(AX=XB\)(缠绕关系,intertwining relation)是其齐次形式,\(AX=XA\) 即交换性。
定理 2.4.4.1(Sylvester) \(AX-XB=C\) 对每个 \(C\) 都有唯一解 ⇔ \(\sigma(A)\cap\sigma(B)=\varnothing\)。\(A,B,C\) 实时解也是实的。
证明 \(X\mapsto AX-XB\) 是 \(M_{n,m}\) 上的线性映射,只需证其核为零。若 \(AX=XB\),则对任意多项式 \(g\) 有 \(g(A)X=Xg(B)\);取 \(g=p_B\)(\(B\) 的特征多项式),由 Cayley–Hamilton,\(p_B(A)X=Xp_B(B)=0\)。而 \(p_B(A)=\prod_j(A-\mu_jI)\)(\(\mu_j\) 为 \(B\) 的特征值),谱不交时每个因子可逆,所以 \(X=0\)。反向类似。∎
原书 2.4.P13 给出更直观的另一种看法:\(X\mapsto AX\) 与 \(X\mapsto XB\) 是两个交换的线性映射,它们之差的特征值恰好是所有差 \(\lambda_i-\mu_j\)。
量化中最常见的两个特例:
- 离散 Lyapunov(Stein)方程 \(\Sigma=A\Sigma A^T+Q\):VAR(1) \(x_{t+1}=Ax_t+\varepsilon_t\)(\(\operatorname{Cov}\varepsilon=Q\))的平稳协方差满足它。\(A\) 可逆时改写为 \(A^{-1}\Sigma-\Sigma A^T=A^{-1}Q\),是 Sylvester 方程,唯一可解 ⇔ \(1/\lambda_i\ne\lambda_j\),即所有 \(\lambda_i\lambda_j\ne1\)。\(\rho(A)<1\)(平稳)时显然满足。
- 连续 Lyapunov 方程 \(A\Sigma+\Sigma A^T=-Q\):多维 Ornstein–Uhlenbeck 过程 \(dx=Ax\,dt+dW\) 的平稳协方差。唯一可解 ⇔ \(\lambda_i+\lambda_j\ne0\);\(A\) 的特征值实部全为负(均值回复)时满足。
金融直觉:Stein 方程是怎么来的?对 \(x_{t+1}=Ax_t+\varepsilon_{t+1}\) 两边取协方差,\(\varepsilon_{t+1}\) 与 \(x_t\) 独立,所以 \(\operatorname{Cov}(x_{t+1})=A\operatorname{Cov}(x_t)A^T+Q\)(用了 \(\operatorname{Cov}(Ax)=A\operatorname{Cov}(x)A^T\),这是组合方差 \(w^T\Sigma w\) 的矩阵推广)。平稳意味着前后两期协方差相同,都记为 \(\Sigma\),就得到 \(\Sigma=A\Sigma A^T+Q\)。 一维时它就是熟悉的 AR(1) 结论:\(x_{t+1}=\phi x_t+\varepsilon\),\(\sigma^2=\phi^2\sigma^2+q\),解出 \(\sigma^2=q/(1-\phi^2)\)。这里需要 \(\phi^2\ne1\),矩阵版的条件"所有 \(\lambda_i\lambda_j\ne1\)"就是它的推广(\(i=j\) 时正是 \(\lambda_i^2\ne1\))。\(\phi\) 越接近 1,平稳方差越大,矩阵版同理:\(A\) 有接近单位圆的特征值时,平稳协方差在对应方向上很大。
推论 2.4.4.2(缠绕关系的分块) 若 \(B=B_1\oplus\cdots\oplus B_k\)、\(C=C_1\oplus\cdots\oplus C_k\),\(i\ne j\) 时 \(\sigma(B_i)\cap\sigma(C_j)=\varnothing\),且 \(AB=CA\),则 \(A\) 也共形块对角。特别地(2.4.4.3),与可对角化矩阵 \(S\Lambda S^{-1}\) 交换的矩阵,在同一组基下按相同特征值分块对角。原书把它总结为一条基本原则:若 \(AX=XB\) 且 \(A,B\) 结构特殊,则 \(X\) 也很可能结构特殊——先把 \(A,B\) 换成标准形再研究。
2.7.4 分块对角化与"几乎可对角化"
定理 2.4.6.1 设 \(A\) 的不同特征值为 \(\lambda_1,\dots,\lambda_d\),重数 \(n_1,\dots,n_d\)。由 Schur 定理(把相同特征值排在一起),\(A\) 酉相似于块上三角 \(T=[T_{ij}]\),\(T_{ii}=\lambda_iI_{n_i}+(\text{严格上三角})\)。进一步,\(A\) 相似于
做法:\(T=\begin{bmatrix}T_{11}&Y\\0&S_2\end{bmatrix}\),\(T_{11}\) 与 \(S_2\) 谱不交,由 Sylvester 定理解出 \(T_{11}X-XS_2=-Y\),再用 \(M=\begin{bmatrix}I&X\\0&I\end{bmatrix}\) 做相似:
这一步一般不能用酉相似完成:若某个 \(T_{ij}\ne0\),块对角化会改变元素平方和。这是通往 Jordan 形的第二步(第 03a 章)。
定理 2.4.7.1(可对角化矩阵稠密) 对任意 \(A\) 和 \(\epsilon>0\),存在特征值互异(因而可对角化)的 \(A(\epsilon)\) 使 \(\|A-A(\epsilon)\|_F^2<\epsilon\)。做法:在 Schur 形的对角线上加小扰动使特征值互异。另有(2.4.7.2):可用非奇异 \(S_\epsilon\) 把 \(A\) 化成非对角元都不超过 \(\epsilon\) 的上三角阵(用 \(\operatorname{diag}(1,\epsilon,\dots,\epsilon^{n-1})\) 压缩超对角元)。
这给了一个常用证明技巧:先对可对角化矩阵证明,再取极限。但也要警惕:\(\begin{bmatrix}0&\epsilon\\0&0\end{bmatrix}\) 与零矩阵任意接近,一个不可对角化一个可对角化——可对角化性本身对扰动不稳定。
2.7.5 Schur 形的唯一性程度
定理 2.4.5.1(要点) 若两个对角线相同、按相同特征值分块的上三角阵 \(T,T'\) 满足 \(WT=T'W\),则 \(W\) 块上三角;若 \(W\) 还是酉的,则 \(W\) 块对角;若每个对角块的第一超对角元都非零且 \(W\) 酉,则 \(W\) 是对角酉阵;若第一超对角元都是正实数,则 \(W=wI\),\(T=T'\)。证明的方法是从左下角开始逐块比较 \(WT=T'W\),每一步用 Sylvester 定理。
2.7.6 交换族与特征值的加法
定理 2.4.8.1 若 \(AB=BA\),则存在特征值的某种排序 \(\alpha_i,\beta_i\) 使 \(A+B\) 的特征值为 \(\alpha_i+\beta_i\),\(AB\) 的特征值为 \(\alpha_i\beta_i\)。因此交换时谱半径次可加、次可乘。
三个例子划出了边界:
- 例 2.4.8.2:\(A=\operatorname{diag}(1,2)\)、\(B=\operatorname{diag}(3,4)\),\(1+4=5\notin\sigma(A+B)=\{4,6\}\)——不是任意配对都行。
- 例 2.4.8.3:\(A=\begin{bmatrix}0&1\\0&0\end{bmatrix}\)、\(B=\begin{bmatrix}0&0\\1&0\end{bmatrix}\) 不交换,\(\sigma(A+B)=\{\pm1\}\) 而 \(\sigma(A)=\sigma(B)=\{0\}\)。谱半径在一般矩阵上不次可加。
- 例 2.4.8.4:存在矩阵对使 \(\alpha A+\beta B\) 的特征值对所有 \(\alpha,\beta\) 都"相加",但 \(AB\) 的特征值不相乘——逆命题不成立。
量化提醒:两个风险来源叠加时,只有当两个协方差(或两个动力学矩阵)交换,"总特征值 = 分特征值之和"才可能成立;一般情况下必须重算整体矩阵。
定理 2.4.8.6 可同时相似上三角化 ⇔ 可同时酉相似上三角化(对 \(S\) 做 QR 即可)。
定理 2.4.8.7(McCoy,了解) \(A_1,\dots,A_m\) 可同时酉三角化 ⇔ 对每个非交换多项式 \(p\) 和每对 \(k,\ell\),\(p(A_1,\dots,A_m)(A_kA_\ell-A_\ell A_k)\) 都幂零。由此可得 Laffey 定理:\(\operatorname{rank}(AB-BA)\le1\) 时 \(A,B\) 可同时三角化(2.4.P11)。
2.7.7 特征值的连续性
定理 2.4.9.2 若 \(A_k\to A\)(逐元),则 \(A_k\) 的特征值在适当配对下收敛到 \(A\) 的特征值:对任意 \(\varepsilon>0\),当 \(k\) 充分大时
证明同时用了 Schur 定理的两面:酉(酉群紧,可以取 \(U_k\) 的收敛子列)与三角(对角线就是特征值,上三角阵的极限仍是上三角阵)。推论(2.4.P1):特征值互异的矩阵构成开集。
连续性只说"会收敛",不说"收敛多快"。对非正规矩阵,特征值可以对扰动极其敏感(第 01 章实战 4 的移位矩阵例子);对 Hermite 矩阵(协方差矩阵)则有 Weyl 不等式这类定量的稳定性(原书第 4 章)。
2.7.8 秩一扰动与双正交原理
定理 2.4.10.1(Brauer) 若 \(Ax=\lambda x\),则对任意 \(v\),\(A+xv^*\) 的特征值是 \(\lambda+v^*x,\lambda_2,\dots,\lambda_n\)。用 Schur 定理证明:以 \(x/\|x\|\) 为第一列做酉相似,\(xv^*\) 只影响第一行。第 01 章已讨论其在收缩法、Google 矩阵中的用途。
定理 2.4.11.1(完整的双正交原理,要点) 设 \(Ax=\lambda x\)、\(y^*A=\mu y^*\)(\(x,y\) 单位向量)。
- \(\lambda\ne\mu\) 时 \(y^*x=0\),且以 \([x\ y\ \cdots]\) 为前两列的酉相似使 \(A\) 的第二行出现额外的零;
- \(\lambda=\mu\) 且 \(y^*x\ne0\) 时,\(A\) 相似于 \([\lambda]\oplus A_{n-1}\),可以把 \(\lambda\) "干净地"分离出来;
- 若 \(x\) 同时是右、左特征向量(\(x^*A=\lambda x^*\),称正规特征向量),这种分离可用酉相似完成。
2.7.9 其他结论(选读)
- 2.4.P2:\(\operatorname{rank}A\ge|\operatorname{tr}A|^2/\operatorname{tr}(A^*A)\)。
- 2.4.P21:以 \(\operatorname{tr}A^{i+j-2}\) 为元素的 Hankel 矩阵 \(K\) 的行列式是判别式 \(\prod_{i<j}(\lambda_j-\lambda_i)^2\);特征值互异 ⇔ \(K\) 非奇异,且 \(\operatorname{rank}K\) 等于不同特征值的个数。
- 2.4.P31–P32:\(\operatorname{tr}(AB-BA)=0\),所以有限维中不可能有 \(AB-BA=cI\)(\(c\ne0\))——量子力学的 Heisenberg 关系只能在无穷维中成立。
- 2.3.P12:复合矩阵 \(C_r(A)\) 的特征值是所有 \(r\) 个特征值的乘积。
量化实战
实战 1:用 QR 解因子回归,QR 就是因子逐步正交化
场景:截面或时间序列因子回归 \(y=X\beta+\varepsilon\)。教科书公式 \(\hat\beta=(X^TX)^{-1}X^Ty\) 要显式构造 \(X^TX\),它的条件数是 \(X\) 的平方(第 02b 章用 SVD 证明 \(\kappa(X^TX)=\kappa(X)^2\))。QR 解法 \(X=QR\)、\(R\hat\beta=Q^Ty\) 只和 \(\kappa(X)\) 打交道。
import numpy as np
rng = np.random.default_rng(42)
# ---------- 1. 两个高度共线的因子:正规方程 vs QR ----------
T = 500
mkt = rng.normal(0, 0.01, T)
size = 0.98 * mkt + 0.002 * rng.normal(0, 0.01, T) # 与市场几乎共线
X = np.column_stack([np.ones(T), mkt, size])
beta_true = np.array([0.0002, 1.1, -0.4])
y = X @ beta_true + rng.normal(0, 0.005, T)
print("cond(X) = %.3e" % np.linalg.cond(X))
print("cond(X^T X) = %.3e (= cond(X)^2 = %.3e)" % (np.linalg.cond(X.T @ X), np.linalg.cond(X) ** 2))
# 正规方程:显式构造 X^T X
b_ne = np.linalg.solve(X.T @ X, X.T @ y)
# QR:X = QR,R b = Q^T y(回代)
Q, R = np.linalg.qr(X) # 瘦 QR
b_qr = np.linalg.solve(R, Q.T @ y) # R 是上三角,这里为简洁直接 solve
b_ls = np.linalg.lstsq(X, y, rcond=None)[0]
print("正规方程 β:", b_ne.round(6))
print("QR β:", b_qr.round(6))
print("lstsq β:", b_ls.round(6))
# 用单精度放大差异:float32 下正规方程失去全部有效数字
X32, y32 = X.astype(np.float32), y.astype(np.float32)
b_ne32 = np.linalg.solve(X32.T @ X32, X32.T @ y32)
Q32, R32 = np.linalg.qr(X32)
b_qr32 = np.linalg.solve(R32, Q32.T @ y32)
print("\nfloat32 正规方程 β:", b_ne32.round(4), " 相对误差 %.2e" % (np.linalg.norm(b_ne32 - b_ls) / np.linalg.norm(b_ls)))
print("float32 QR β:", b_qr32.round(4), " 相对误差 %.2e" % (np.linalg.norm(b_qr32 - b_ls) / np.linalg.norm(b_ls)))
# ---------- 2. QR = 逐步正交化因子(Gram–Schmidt) ----------
# 列顺序:常数、市场、规模、价值;价值因子与前两者相关
value = 0.5 * mkt - 0.3 * size + rng.normal(0, 0.006, T)
F = np.column_stack([np.ones(T), mkt, size, value])
Qf, Rf = np.linalg.qr(F)
s = np.sign(np.diag(Rf)); Qf, Rf = Qf * s, (Rf.T * s).T # 统一为 R 对角元为正
# 第 4 列:价值因子对 [1, mkt, size] 回归的残差,归一化后就是 q4
G = F[:, :3]
resid = value - G @ np.linalg.lstsq(G, value, rcond=None)[0]
print("\nq4 与 '价值对前三列回归残差/范数' 的最大差:", np.abs(Qf[:, 3] - resid / np.linalg.norm(resid)).max())
print("r44 = %.6f, 残差范数 = %.6f(原书 2.1.P18:r_kk = 第 k 列到前 k-1 列张成空间的距离)"
% (Rf[3, 3], np.linalg.norm(resid)))
# ---------- 3. Householder 反射:把向量变成 ‖x‖e1 ----------
x = rng.normal(size=5)
w = x - np.linalg.norm(x) * np.eye(5)[0]
H = np.eye(5) - 2 * np.outer(w, w) / (w @ w)
print("\nHx =", (H @ x).round(10))
print("H 对称正交:", np.allclose(H, H.T), np.allclose(H @ H, np.eye(5)), " det(H) = %.1f" % np.linalg.det(H))
关键输出:
cond(X) = 6.881e+04
cond(X^T X) = 4.735e+09 (= cond(X)^2 = 4.735e+09)
正规方程 β: [ 1.7600000e-04 1.0177651e+01 -9.6350400e+00]
QR β: [ 1.7600000e-04 1.0177651e+01 -9.6350400e+00]
lstsq β: [ 1.7600000e-04 1.0177651e+01 -9.6350400e+00]
float32 正规方程 β: [ 2.0000e-04 9.8375e+00 -9.2880e+00] 相对误差 3.47e-02
float32 QR β: [ 2.00000e-04 1.01778e+01 -9.63520e+00] 相对误差 1.50e-05
q4 与 '价值对前三列回归残差/范数' 的最大差: 3.740063814205996e-14
r44 = 0.133932, 残差范数 = 0.133932(原书 2.1.P18:r_kk = 第 k 列到前 k-1 列张成空间的距离)
Hx = [ 1.69495937 0. -0. -0. 0. ]
H 对称正交: True True det(H) = -1.0
读法:
- \(\kappa(X^TX)=\kappa(X)^2\) 精确成立。双精度下 \(\kappa\approx5\times10^9\) 还不至于出事,三种方法给出相同的解;换成单精度(约 7 位有效数字),正规方程的相对误差是 3.5%,QR 只有 \(1.5\times10^{-5}\)。经验法则:正规方程大约损失 \(2\log_{10}\kappa(X)\) 位有效数字,QR 损失 \(\log_{10}\kappa(X)\) 位。在 GPU(常用单精度)上跑大规模截面回归时,这个差别是实打实的。
- 注意解本身:真实系数是 \((1.1,-0.4)\),估计却是 \((10.2,-9.6)\)。这不是数值误差,三种算法完全一致——这是共线性导致的统计不稳定(估计方差巨大)。数值稳定的算法只能保证"准确地算出一个不靠谱的解"。解决统计问题要靠截断 SVD、岭回归等正则化(第 02b 章实战 2)。
- QR 的第 4 列 \(q_4\) 与"价值因子对常数、市场、规模回归后的残差(归一化)"在机器精度内相同,\(r_{44}\) 就是残差范数。因子研究中"把新因子对已有因子正交化再检验增量信息",一次 QR 就全部完成,而且排列顺序决定了谁"吃掉"共同部分。
- Householder 反射把任意向量变成 \(\|x\|e_1\),它对称、自逆、行列式为 \(-1\)。
实战 2:VAR(1) 的 Schur 分析、平稳协方差与瞬时放大
场景:三变量 VAR(1) \(x_{t+1}=Ax_t+\varepsilon_t\)(比如利率水平、斜率、通胀预期),需要判断平稳性、求平稳协方差、计算多步预测所需的 \(A^k\)。
import numpy as np
from scipy.linalg import schur, solve_discrete_lyapunov, solve_sylvester
rng = np.random.default_rng(7)
# ---------- 1. VAR(1):x_{t+1} = A x_t + e_t,用实 Schur 形看稳定性 ----------
A = np.array([[0.50, 0.80, 0.00],
[-0.30, 0.60, 0.10],
[0.05, 0.00, 0.90]])
T_real, Z = schur(A, output="real") # A = Z T Z^T,Z 实正交
print("实 Schur 形 T(2x2 块对应共轭复特征值对):\n", T_real.round(4))
print("‖Z T Z^T - A‖ =", np.linalg.norm(Z @ T_real @ Z.T - A))
ev = np.linalg.eigvals(A)
print("特征值:", np.round(ev, 4), " 谱半径 ρ(A) = %.4f" % max(abs(ev)))
# ---------- 2. 平稳协方差:Σ = A Σ A^T + Q(Stein/离散 Lyapunov 方程) ----------
Q = np.diag([1.0, 0.5, 0.2])
Sigma = solve_discrete_lyapunov(A, Q)
print("\n平稳协方差 Σ:\n", Sigma.round(4))
print("残差 ‖AΣA^T + Q - Σ‖ =", np.linalg.norm(A @ Sigma @ A.T + Q - Sigma))
# 等价的 vec 形式:(I - A⊗A) vec Σ = vec Q,唯一可解 ⟺ 所有 λ_i λ_j ≠ 1
K = np.eye(9) - np.kron(A, A)
print("I - A⊗A 的最小奇异值 = %.4f(>0 即唯一可解)" % np.linalg.svd(K, compute_uv=False).min())
# 模拟核对
x = np.zeros(3); X = []
L = np.linalg.cholesky(Q)
for t in range(200_000):
x = A @ x + L @ rng.normal(size=3); X.append(x)
print("模拟样本协方差:\n", np.cov(np.array(X[1000:]).T).round(3))
# Sylvester 方程 AX - XB = C:谱不交时唯一解
B = np.diag([2.0, -1.5])
C = rng.normal(size=(3, 2))
Xs = solve_sylvester(A, -B, C) # scipy 解 AX + XB' = C,这里 B' = -B
print("\nSylvester 残差:", np.linalg.norm(A @ Xs - Xs @ B - C))
# ---------- 3. 偏离正规性与瞬时放大 ----------
def dep_normal(M):
return np.sqrt(max(np.linalg.norm(M, "fro") ** 2 - np.sum(abs(np.linalg.eigvals(M)) ** 2), 0))
N1 = np.diag([0.9, 0.8]) # 正规
N2 = np.array([[0.9, 5.0], [0.0, 0.8]]) # 同特征值,强非正规
for name, M in [("正规 diag(0.9,0.8)", N1), ("非正规 [[0.9,5],[0,0.8]]", N2), ("VAR 系数 A", A)]:
norms = [np.linalg.norm(np.linalg.matrix_power(M, k), 2) for k in range(0, 40)]
print("%-24s ρ=%.2f 偏离正规性=%.3f max_k‖M^k‖=%.2f(k=%d)"
% (name, max(abs(np.linalg.eigvals(M))), dep_normal(M), max(norms), int(np.argmax(norms))))
# ---------- 4. Cayley–Hamilton:高次幂降为 n-1 次多项式 ----------
p = np.poly(A) # 特征多项式系数 [1, c2, c1, c0]
k = 25
r = np.polydiv(np.r_[1.0, np.zeros(k)], p)[1] # t^k = h(t) p(t) + r(t)
Ak_ch = sum(c * np.linalg.matrix_power(A, i) for i, c in enumerate(r[::-1]))
print("\nA^25 用 C–H 余式计算的误差:", np.abs(Ak_ch - np.linalg.matrix_power(A, k)).max())
关键输出:
实 Schur 形 T(2x2 块对应共轭复特征值对):
[[ 0.5446 -0.8007 -0.0741]
[ 0.3015 0.5446 0.0468]
[ 0. 0. 0.9109]]
‖Z T Z^T - A‖ = 9.94379369323482e-16
特征值: [0.5446+0.4913j 0.5446-0.4913j 0.9109+0.j ] 谱半径 ρ(A) = 0.9109
平稳协方差 Σ:
[[2.3911 0.1624 0.1898]
[0.1624 1.0366 0.053 ]
[0.1898 0.053 1.174 ]]
残差 ‖AΣA^T + Q - Σ‖ = 1.32021460364232e-16
I - A⊗A 的最小奇异值 = 0.1664(>0 即唯一可解)
模拟样本协方差:
[[2.379 0.156 0.194]
[0.156 1.031 0.052]
[0.194 0.052 1.189]]
Sylvester 残差: 1.246531298897714e-15
正规 diag(0.9,0.8) ρ=0.90 偏离正规性=0.000 max_k‖M^k‖=1.00(k=0)
非正规 [[0.9,5],[0,0.8]] ρ=0.90 偏离正规性=5.000 max_k‖M^k‖=13.48(k=6)
VAR 系数 A ρ=0.91 偏离正规性=0.507 max_k‖M^k‖=1.04(k=1)
A^25 用 C–H 余式计算的误差: 6.245004513516506e-16
读法:
- 实 Schur 形左上角是一个 \(2\times2\) 块,对应共轭复特征值 \(0.5446\pm0.4913i\)(模约 0.733,辐角约 0.734 弧度,伪周期约 \(2\pi/0.734\approx8.6\) 期);它对角元相等但非对角元不互为相反数——这正是定理 2.3.4(b) 说的"实正交相似下 \(2\times2\) 块不一定是特殊形式"。右下角是实特征值 0.9109,决定最慢的衰减速度(半衰期 \(\ln0.5/\ln0.9109\approx7.4\) 期)。
- 平稳协方差由 Stein 方程直接解出,与 20 万步模拟的样本协方差吻合。\(I-A\otimes A\) 的最小奇异值 0.166 远离 0,说明方程良态;当 \(\rho(A)\to1\) 时它趋于 0,平稳协方差爆炸,这就是"近单位根"在矩阵层面的表现。
- 两个矩阵谱半径都是 0.9,正规的那个 \(\|M^k\|\) 单调下降;非正规的那个(偏离正规性为 5)在第 6 期把冲击放大 13.5 倍才开始衰减。估计 VAR 时只看特征根判断"系统稳定、冲击很快消散"可能严重低估中期风险,应同时看 \(\|A^k\|\) 或脉冲响应。
- 用 Cayley–Hamilton 余式计算 \(A^{25}\),只需 \(I,A,A^2\) 的线性组合,结果与直接乘幂一致。
实战 3:Jacobi 方法求协方差矩阵的特征分解
import numpy as np
rng = np.random.default_rng(1)
n = 8
R = rng.normal(size=(250, n)) @ rng.normal(size=(n, n)) # 相关的收益
S = np.cov(R.T) # 样本协方差(实对称)
def off(M):
return np.sqrt(np.sum(M ** 2) - np.sum(np.diag(M) ** 2))
def jacobi_eig(A, tol=1e-12, max_rot=10_000):
A = A.copy(); n = A.shape[0]; V = np.eye(n); hist = [off(A)]
for _ in range(max_rot):
Ao = np.abs(A - np.diag(np.diag(A)))
i, j = np.unravel_index(np.argmax(Ao), Ao.shape) # 经典 Jacobi:选最大非对角元
if Ao[i, j] < tol:
break
theta = 0.5 * np.arctan2(2 * A[i, j], A[j, j] - A[i, i]) # 使旋转后 b_ij = 0
c, s = np.cos(theta), np.sin(theta)
G = np.eye(n); G[i, i] = G[j, j] = c; G[i, j] = s; G[j, i] = -s
A = G.T @ A @ G; V = V @ G; hist.append(off(A))
return np.diag(A), V, hist
lam, V, hist = jacobi_eig(S)
print("旋转次数:", len(hist) - 1)
print("非对角 Frobenius 质量(每 10 次旋转):", np.array(hist[::10]).round(6)[:8])
bound = 1 - 2 / (n * n - n)
print("理论每步压缩因子 ≤ %.4f,实测前 20 步平均压缩比 %.4f"
% (bound, np.mean(np.array(hist[1:21]) ** 2 / np.array(hist[:20]) ** 2)))
print("Jacobi 特征值:", np.sort(lam).round(6))
print("eigh 特征值:", np.linalg.eigvalsh(S).round(6))
print("V 正交:", np.allclose(V.T @ V, np.eye(n)), " ‖S V - V Λ‖ =", np.linalg.norm(S @ V - V * lam))
关键输出:
旋转次数: 99
非对角 Frobenius 质量(每 10 次旋转): [2.3385926e+01 8.2197490e+00 3.0412990e+00 1.1247410e+00 2.3916600e-01
6.2879000e-02 5.8230000e-03 1.6900000e-04]
理论每步压缩因子 ≤ 0.9643,实测前 20 步平均压缩比 0.8160
Jacobi 特征值: [ 0.030384 0.528662 1.084879 3.389409 7.101252 12.289573 21.935015
25.527836]
eigh 特征值: [ 0.030384 0.528662 1.084879 3.389409 7.101252 12.289573 21.935015
25.527836]
V 正交: True ‖S V - V Λ‖ = 1.046536950214474e-12
读法:每次旋转后非对角平方和的压缩比(实测平均 0.82)优于理论上界 \(1-2/(n^2-n)\approx0.964\);后期进入二次收敛,99 次旋转后达到 \(10^{-12}\)。特征值与 LAPACK 一致,旋转之积 \(V\) 是正交的特征向量矩阵。Jacobi 方法比三对角化 + QR 迭代慢,但它天然并行、对小特征值的相对精度高,在需要精确估计最小特征值(最小方差组合、协方差矩阵条件数)时仍有用武之地。
本章小结
酉矩阵是保持长度与内积的线性变换,等价于列(行)标准正交的方阵;酉群是紧群,这使"取收敛子列再取极限"成为可用的证明工具。Givens 旋转(一次消一个元素)和 Householder 反射(一次压缩一整列)是构造一切酉变换算法的积木;用 Householder 反射逐列消元得到 QR 分解,列满秩时 \(R\) 对角元为正的 QR 分解唯一。QR 等价于 Gram–Schmidt:\(q_k\) 是第 \(k\) 列对前面各列回归的归一化残差,\(r_{kk}\) 是残差范数;用 QR 解最小二乘只受 \(\kappa(X)\) 影响,正规方程受 \(\kappa(X)^2\) 影响。酉相似保持 Frobenius 范数,由此可快速判定不酉相似。Schur 定理说任意方阵都酉相似于上三角阵,特征值按任意顺序排在对角线上;实矩阵则实正交相似于上拟三角阵,\(2\times2\) 块对应共轭复特征值(振荡模态)。Schur 不等式 \(\sum|\lambda_i|^2\le\|A\|_F^2\) 的差值度量偏离正规性,它决定了 \(\|A^k\|\) 可能出现的瞬时放大。Schur 定理推出了一连串基本结果:\(p(A)\) 的特征值计重数为 \(p(\lambda_i)\);Cayley–Hamilton 定理;Sylvester 方程 \(AX-XB=C\) 唯一可解 ⇔ 谱不交(包含离散、连续 Lyapunov 方程);任意方阵相似于"每块只有一个特征值"的块对角阵;可对角化矩阵稠密;交换矩阵的特征值可配对相加相乘;特征值连续依赖于矩阵元素。
| 概念/公式 | 表达式 | 用途 |
|---|---|---|
| 酉矩阵 | \(U^*U=I\iff|Ux|_2=|x|_2\) | 换标准正交基,数值稳定 |
| Givens 旋转 | \(U(\theta;i,j)\) | 逐个消元、QR 更新 |
| Householder 反射 | \(U_w=I-2ww^*/w^*w\) | 整列消元、QR、三对角化 |
| QR 分解 | \(A=QR\),列满秩时唯一 | 稳定的最小二乘 |
| QR 的含义 | \(r_{kk}=\operatorname{dist}(a_k,\operatorname{span}\{a_1..a_{k-1}\})\) | 因子逐步正交化、共线性诊断 |
| Hadamard 不等式 | \(\vert \det A\vert \le\prod|a_i|_2\) | 广义方差上界 |
| Cholesky | \(A^*A=R^*R=LL^*\) | 模拟相关随机数 |
| Frobenius 不变 | \(|UAV|_F=|A|_F\) | 判定不酉相似 |
| Schur 定理 | \(A=UTU^*\),\(t_{ii}=\lambda_i\) | 几乎所有推论的出发点 |
| 实 Schur 形 | \(Q^TAQ\) 上拟三角 | VAR 特征根、振荡模态 |
| Schur 不等式 | \(\sum\vert \lambda_i\vert ^2\le|A|_F^2\) | 偏离正规性、瞬时放大 |
| 谱映射(计重数) | \(\sigma(p(A))=\{p(\lambda_i)\}\) | \(\operatorname{tr}A^k=\sum\lambda_i^k\) |
| Cayley–Hamilton | \(p_A(A)=0\) | \(A^k\)、\(A^{-1}\) 降为低次多项式 |
| Sylvester 方程 | \(AX-XB=C\) 唯一可解 \(\iff\sigma(A)\cap\sigma(B)=\varnothing\) | Lyapunov 方程、块对角化 |
| 离散 Lyapunov | \(\Sigma=A\Sigma A^T+Q\),需 \(\lambda_i\lambda_j\ne1\) | VAR 平稳协方差 |
| 交换矩阵的特征值 | \(\sigma(A+B)=\{\alpha_i+\beta_i\}\) | 风险叠加的前提 |
| 特征值连续性 | \(A_k\to A\Rightarrow\lambda(A_k)\to\lambda(A)\) | 估计误差的定性保证 |
练习
基础
- 证明 \(2\times2\) 实正交矩阵只有旋转和反射两类,并说明 Householder 矩阵 \(U_w\)(\(w\in\mathbf R^2\))属于哪一类。 提示:列为单位向量且正交,第一列写成 \((\cos\theta,\sin\theta)^T\),第二列只有两种选择;\(\det U_w=-1\),属于反射。
- 用 Frobenius 范数判断 \(\begin{bmatrix}1&2\\0&3\end{bmatrix}\) 与 \(\begin{bmatrix}1&0\\0&3\end{bmatrix}\) 是否酉相似,是否相似。 答案要点:平方和 14 与 10 不等,不酉相似;特征值互异且相同,相似。
- 对 \(A=\begin{bmatrix}3&1\\-2&0\end{bmatrix}\),用 Cayley–Hamilton 求 \(A^5\) 与 \(A^{-2}\)。 答案要点:\(A^5=31A-30I\);\(A^{-2}=-\frac34A+\frac74I\)。
- 设 \(X=[\mathbf 1,\ f]\)(常数列与一个因子列),对 \(X\) 做 QR。说明 \(q_2\) 是"去均值后再归一化为单位长度的 \(f\)",且 \(r_{22}=\sqrt{T}\cdot\operatorname{sd}(f)\)(标准差取除以 \(T\) 的口径)。 提示:\(f\) 对常数回归的残差就是 \(f-\bar f\)。
- 证明 Sylvester 方程 \(AX-XB=C\) 在 \(A\) 的特征值全在单位圆内、\(B\) 的特征值全在单位圆外时有唯一解。 提示:谱不交。
- 设 VAR(1) 的 \(A=\begin{bmatrix}0.5&0.4\\0&0.5\end{bmatrix}\)。求 \(\rho(A)\)、偏离正规性 \(\Delta(A)\),并用 Cayley–Hamilton 写出 \(A^k\) 的闭式。 答案要点:\(\rho=0.5\),\(\Delta=0.16\)(\(\sqrt\Delta=0.4\));\(p_A(t)=(t-0.5)^2\),\(A^k=0.5^kI+k\,0.5^{k-1}(A-0.5I)\),右上元 \(0.4k\,0.5^{k-1}\)。
进阶
- 证明 Schur 定理的实版本 2.3.4(b) 可以由 (a) 加实 QR 分解推出,并解释为什么 \(2\times2\) 块失去了特殊形式。
- 设 \(A\) 可逆。把离散 Lyapunov 方程 \(\Sigma-A\Sigma A^T=Q\) 改写成 Sylvester 方程,推出唯一可解条件 \(\lambda_i\lambda_j\ne1\);再用 \(\operatorname{vec}\) 与 Kronecker 积写出 \((I-A\otimes A)\operatorname{vec}\Sigma=\operatorname{vec}Q\),说明 \(I-A\otimes A\) 的特征值是 \(1-\lambda_i\lambda_j\)。 提示:\(A\otimes A\) 的特征值是 \(\lambda_i\lambda_j\)(对 \(A\) 做 Schur 分解,\(A\otimes A\) 酉相似于上三角的 \(T\otimes T\))。
- 证明:若 \(A\) 正规(\(A=U\Lambda U^*\)),则 \(\|A^k\|_2=\rho(A)^k\);构造一个 \(2\times2\) 上三角矩阵,谱半径为 0.95,但 \(\max_k\|A^k\|_2>10\)。 提示:取 \(\begin{bmatrix}0.95&c\\0&0.95\end{bmatrix}\),\(A^k\) 的右上元为 \(kc\,0.95^{k-1}\)。
- 用 Newton 恒等式由 \(\operatorname{tr}A=3\)、\(\operatorname{tr}A^2=5\)、\(\operatorname{tr}A^3=9\) 求 \(3\times3\) 矩阵 \(A\) 的特征多项式。 答案要点:\(a_2=-\mu_1=-3\);\(2a_1+a_2\mu_1+\mu_2=0\Rightarrow a_1=2\);\(3a_0+a_1\mu_1+a_2\mu_2+\mu_3=0\Rightarrow a_0=0\)。\(p(t)=t^3-3t^2+2t\),特征值 0、1、2。
原书推荐习题
- 2.1.P18、P23:QR 中 \(r_{kk}\) 的几何意义;Hadamard 不等式。
- 2.1.P27–P29:用平面旋转做 QR(可编程)。
- 2.2.P1:Jacobi 方法的收敛性;2.2.P10:Fourier 矩阵与循环矩阵的对角化。
- 2.3.P1:显式构造以给定单位向量为首列的酉矩阵;2.3.P10:\(|\det A|\le c^nn^{n/2}\)。
- 2.4.P3:无除法的 Cayley–Hamilton 证明与 \(\operatorname{adj}A\) 的多项式表示;2.4.P9–P10:Newton 恒等式;2.4.P13:Sylvester 方程的谱刻画;2.4.P30:用余式计算高次幂。
原书对照
| 本章小节 | 原书小节 | 书页 | PDF 页 |
|---|---|---|---|
| 2.1 为什么要"酉" | 2.0 Introduction | 83 | 103 |
| 2.2–2.4 酉矩阵、Givens/Householder、QR | 2.1 Unitary matrices and the QR factorization(含习题 2.1) | 83–94 | 103–114 |
| 2.5 酉相似 | 2.2 Unitary similarity(含习题 2.2) | 94–101 | 114–121 |
| 2.6 Schur 三角化 | 2.3 Unitary and real orthogonal triangularizations(含习题 2.3) | 101–108 | 121–128 |
| 2.7 Schur 定理的推论 | 2.4 Consequences of Schur's triangularization theorem(正文 p.108–124,习题 p.124–131) | 108–131 | 128–151 |
下一章(第 02b 章)把 Schur 定理推到最"漂亮"的一类矩阵——正规矩阵,得到谱定理;再放开到酉等价,得到对任意长方阵都成立的奇异值分解。