量化交易中文教材

第 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 方程求平稳协方差——底层都是本章的酉变换。

学习目标

读完本章,你应当能够:

  1. 说清酉矩阵(实情形为正交矩阵)的几种等价刻画,理解"保长度 ⇔ 保内积 ⇔ 列标准正交",并会用 Givens 旋转与 Householder 反射构造酉矩阵。
  2. 掌握 QR 分解的存在性、唯一性和 Householder 构造;理解为什么用 QR 解最小二乘比正规方程稳定,以及 QR 与因子逐步正交化(Gram–Schmidt)的关系。
  3. 理解酉相似保持 Frobenius 范数,并能陈述和证明 Schur 三角化定理及其实版本(实 Schur 形),会用 Schur 不等式度量"偏离正规性"。
  4. 会用 Schur 定理推出 Cayley–Hamilton 定理、\(p(A)\) 的特征值、Sylvester 方程唯一可解条件、分块对角化和特征值的连续性。
  5. 能把这些工具用在量化场景: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^*\),相似变换就成了

\[A\mapsto U^*AU .\]

这叫酉相似(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\),则

\[0=\Big\|\sum\alpha_ix_i\Big\|_2^2=\sum_{i,j}\bar\alpha_i\alpha_jx_i^*x_j=\sum_i|\alpha_i|^2 .\]

定义 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\),则

\[\operatorname{rank}U_{12}=\operatorname{rank}U_{21},\qquad\operatorname{rank}U_{22}=\operatorname{rank}U_{11}+n-2k .\]

特别地 \(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=I-\frac{2}{w^*w}ww^* .\]

性质:\(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\),则

\[U_1A=\begin{bmatrix}r_1&\star\\0&A_2\end{bmatrix}.\]

对 \(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}\),所以

\[|\det A|\le\prod_{i=1}^n\|a_i\|_2 ,\]

等号成立当且仅当某列为零或各列两两正交。几何含义:平行多面体的体积不超过各棱长之积。用于协方差矩阵 \(\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\) 酉,可以不同),则

\[\sum_{i,j}|a_{ij}|^2=\sum_{i,j}|b_{ij}|^2 ,\]

即 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\),而最大元的平方不小于平均值,所以

\[\sum_{p\ne q}b_{pq}^2\le\Big(1-\frac{2}{n^2-n}\Big)\sum_{p\ne q}a_{pq}^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^*\),

\[\lambda_\ell=\sum_{k=0}^{n-1}a_{k+1}\omega^{k(\ell-1)},\qquad\ell=1,\dots,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\),

\[U_1^*AU_1=\begin{bmatrix}\lambda_1&\star\\0&A_1\end{bmatrix},\qquad A_1\in M_{n-1},\]

\(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\),于是

\[\sum_{i=1}^n|\lambda_i|^2\le\sum_{i,j}|a_{ij}|^2=\operatorname{tr}(A^*A),\tag{2.3.2a}\]

等号当且仅当 \(T\) 是对角阵,即 \(A\) 可酉对角化(第 02b 章将证明这等价于 \(A\) 正规)。差值

\[\Delta(A)=\|A\|_F^2-\sum|\lambda_i|^2=\sum_{i<j}|t_{ij}|^2\]

称为偏离正规性(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):

\[A=\begin{bmatrix}1&3&0\\0&2&4\\0&0&3\end{bmatrix},\qquad B=\begin{bmatrix}1&0&0\\0&2&5\\0&0&3\end{bmatrix},\]

特征值都是 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\) 块形如

\[\begin{bmatrix}a&b\\-b&a\end{bmatrix},\qquad b>0,\quad a\pm ib\text{ 是 }A\text{ 的特征值},\]
且对角块的顺序可任意指定。 (b) 存在实正交 \(Q\) 使 \(Q^TAQ\) 为实上拟三角阵,\(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\),所以

\[A^2=3A-2I,\quad A^3=7A-6I,\quad A^4=15A-14I,\quad A^{-1}=-\tfrac12A+\tfrac32I=\begin{bmatrix}0&-1/2\\1&3/2\end{bmatrix}.\]

一般地,若 \(p_A(t)=t^n+a_{n-1}t^{n-1}+\cdots+a_1t+a_0\) 且 \(a_0\ne0\),则

\[A^{-1}=-\frac1{a_0}\big(A^{n-1}+a_{n-1}A^{n-2}+\cdots+a_1I\big).\]

更系统的做法(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\) 的多项式:

\[\operatorname{adj}A=(-1)^{n-1}\big(A^{n-1}+a_{n-1}A^{n-2}+\cdots+a_2A+a_1I\big).\]

注意 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\),则

\[ka_{n-k}+a_{n-k+1}\mu_1+\cdots+a_{n-1}\mu_{k-1}+\mu_k=0,\qquad k=1,\dots,n-1 .\]

所以前 \(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 方程

形如

\[AX-XB=C\qquad(A\in M_n,\ B\in M_m,\ X,C\in M_{n,m})\]

的方程称为 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_{11}\oplus T_{22}\oplus\cdots\oplus T_{dd}.\]

做法:\(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}\) 做相似:

\[M^{-1}TM=\begin{bmatrix}T_{11}&T_{11}X-XS_2+Y\\0&S_2\end{bmatrix}=\begin{bmatrix}T_{11}&0\\0&S_2\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\) 充分大时

\[\min_{\pi\in S_n}\max_i|\lambda_{\pi(i)}(A_k)-\lambda_i(A)|\le\varepsilon .\]

证明同时用了 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

读法:

  1. \(\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(常用单精度)上跑大规模截面回归时,这个差别是实打实的。
  2. 注意解本身:真实系数是 \((1.1,-0.4)\),估计却是 \((10.2,-9.6)\)。这不是数值误差,三种算法完全一致——这是共线性导致的统计不稳定(估计方差巨大)。数值稳定的算法只能保证"准确地算出一个不靠谱的解"。解决统计问题要靠截断 SVD、岭回归等正则化(第 02b 章实战 2)。
  3. QR 的第 4 列 \(q_4\) 与"价值因子对常数、市场、规模回归后的残差(归一化)"在机器精度内相同,\(r_{44}\) 就是残差范数。因子研究中"把新因子对已有因子正交化再检验增量信息",一次 QR 就全部完成,而且排列顺序决定了谁"吃掉"共同部分。
  4. 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

读法:

  1. 实 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\) 期)。
  2. 平稳协方差由 Stein 方程直接解出,与 20 万步模拟的样本协方差吻合。\(I-A\otimes A\) 的最小奇异值 0.166 远离 0,说明方程良态;当 \(\rho(A)\to1\) 时它趋于 0,平稳协方差爆炸,这就是"近单位根"在矩阵层面的表现。
  3. 两个矩阵谱半径都是 0.9,正规的那个 \(\|M^k\|\) 单调下降;非正规的那个(偏离正规性为 5)在第 6 期把冲击放大 13.5 倍才开始衰减。估计 VAR 时只看特征根判断"系统稳定、冲击很快消散"可能严重低估中期风险,应同时看 \(\|A^k\|\) 或脉冲响应。
  4. 用 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)\) 估计误差的定性保证

练习

基础

  1. 证明 \(2\times2\) 实正交矩阵只有旋转和反射两类,并说明 Householder 矩阵 \(U_w\)(\(w\in\mathbf R^2\))属于哪一类。 提示:列为单位向量且正交,第一列写成 \((\cos\theta,\sin\theta)^T\),第二列只有两种选择;\(\det U_w=-1\),属于反射。
  2. 用 Frobenius 范数判断 \(\begin{bmatrix}1&2\\0&3\end{bmatrix}\) 与 \(\begin{bmatrix}1&0\\0&3\end{bmatrix}\) 是否酉相似,是否相似。 答案要点:平方和 14 与 10 不等,不酉相似;特征值互异且相同,相似。
  3. 对 \(A=\begin{bmatrix}3&1\\-2&0\end{bmatrix}\),用 Cayley–Hamilton 求 \(A^5\) 与 \(A^{-2}\)。 答案要点:\(A^5=31A-30I\);\(A^{-2}=-\frac34A+\frac74I\)。
  4. 设 \(X=[\mathbf 1,\ f]\)(常数列与一个因子列),对 \(X\) 做 QR。说明 \(q_2\) 是"去均值后再归一化为单位长度的 \(f\)",且 \(r_{22}=\sqrt{T}\cdot\operatorname{sd}(f)\)(标准差取除以 \(T\) 的口径)。 提示:\(f\) 对常数回归的残差就是 \(f-\bar f\)。
  5. 证明 Sylvester 方程 \(AX-XB=C\) 在 \(A\) 的特征值全在单位圆内、\(B\) 的特征值全在单位圆外时有唯一解。 提示:谱不交。
  6. 设 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}\)。

进阶

  1. 证明 Schur 定理的实版本 2.3.4(b) 可以由 (a) 加实 QR 分解推出,并解释为什么 \(2\times2\) 块失去了特殊形式。
  2. 设 \(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\))。
  3. 证明:若 \(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}\)。
  4. 用 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 定理推到最"漂亮"的一类矩阵——正规矩阵,得到谱定理;再放开到酉等价,得到对任意长方阵都成立的奇异值分解。