第 28 章 矩阵运算
本章对应原书第 28 章。量化研究里最常写的一行代码大概是"解一个线性方程组"：回归系数、最小方差权重、对冲比率、风险模型里的因子收益，背后都是 \(Ax=b\)。本章从算法的角度回答三个问题：怎样解得快（LUP 分解，\(\Theta(n^3)\) 分解、\(\Theta(n^2)\) 求解）；怎样解得稳（选主元，不求逆）；为什么求逆和矩阵乘法一样难。最后讲对称正定矩阵和最小二乘。矩阵理论本身见第 01 册（QR 分解见第 02a 章，奇异值分解见第 02b 章，条件数见第 05b 章，Cholesky 分解见第 07a 章，舒尔补见第 07d 章）；回归的统计性质见第 03、05 册。本章侧重算法和工程。
学习目标
读完本章，你应当能够：
- 解释 LUP 分解 \(PA=LU\) 如何把解方程组拆成前代和回代两步，写出 LU、LUP 分解的伪代码并分析 \(\Theta(n^3)\) 的复杂度。
- 用舒尔补理解高斯消元的每一步，说明部分选主元为什么既避免除以零又改善数值稳定性。
- 理解"矩阵求逆与矩阵乘法同样难"的两个定理，知道实践中为什么不该先求逆再相乘。
- 证明对称正定矩阵的顺序主子矩阵和舒尔补仍对称正定、LU 分解的主元全为正；把舒尔补与条件协方差、对冲残余风险联系起来。
- 推导最小二乘的正规方程 \(A^TAc=A^Ty\)，知道正规方程会把条件数平方，在因子高度共线时改用 QR 或 SVD。
- 用 \(O(n)\) 的三对角求解（Thomas 算法）构造自然三次样条收益率曲线。
读前导读
这一章在解决什么问题。 你在 CFA 二级做多元回归时，软件一瞬间就给出了系数；算最小方差组合时，公式里有个 \(\Sigma^{-1}\)。这一章讲的是这些"一瞬间"背后计算机实际在做什么，以及什么时候它会悄悄算错。核心是解线性方程组 \(Ax=b\)：回归系数是它的解，最小方差权重是它的解，多个期货的对冲比率也是它的解。
本章的答案可以浓缩成三句话。第一，不要先求逆矩阵再相乘，而是把 \(A\) 拆成"下三角 × 上三角"（LU 分解），然后像剥洋葱一样一个未知数一个未知数地解出来；这就是中学"消元法"的系统化版本。第二，消元时要挑绝对值最大的数做除数（选主元），否则除以一个极小的数会把舍入误差放大到面目全非。第三，协方差矩阵这类"对称正定"矩阵性质特别好，不用选主元，还有一半计算量的专用方法（Cholesky）。贯穿全章的一个对象是舒尔补：它既是消元法每一步剩下的那个小矩阵，又是统计里的条件协方差，也是你做最小方差对冲后剩下的风险。最后，本章解释了为什么"正规方程"这一教科书公式在因子高度相关时会给出完全错误的回归系数，以及该用什么代替。
需要先想起来的数学。
- 矩阵乘法、转置、逆。 \(AB\) 的第 \((i,j)\) 元是 \(A\) 第 \(i\) 行与 \(B\) 第 \(j\) 列的内积；\(A^T\) 是把行变成列；\(A^{-1}\) 满足 \(AA^{-1}=I\)（\(I\) 是对角线为 1 的单位矩阵）。例：\(\begin{pmatrix}2&0\\0&4\end{pmatrix}^{-1}=\begin{pmatrix}0.5&0\\0&0.25\end{pmatrix}\)。注意 \((AB)^T=B^TA^T\)，顺序要反过来。见 第 00 册第 06 章 线性代数速成。
- 非奇异、秩与行列式。 \(A\) 非奇异（可逆）等价于：\(Ax=b\) 有唯一解；\(A\) 的各列线性无关（没有哪一列能由其他列组合出来）；\(\det A\ne0\)。秩是线性无关的列（或行）的最大个数。回归里"完全共线"就是设计矩阵不满秩。同见第 06 章。
- 二次型与正定。 \(x^TAx=\sum_i\sum_ja_{ij}x_ix_j\)。若 \(A\) 是协方差矩阵 \(\Sigma\)、\(x\) 是组合权重 \(w\)，\(w^T\Sigma w\) 就是组合方差。正定指对任何非零 \(x\) 都有 \(x^TAx>0\)：任何非零组合都有正的方差。例：\(\Sigma=\begin{pmatrix}0.04&0.01\\0.01&0.09\end{pmatrix}\)，\(w=(1,-1)\) 时 \(w^T\Sigma w=0.04-0.02+0.09=0.11>0\)。同见第 06 章。
- 偏导数与求极值。 多变量函数对某一个变量求导、其余视为常数，叫偏导数；最小值点处所有偏导数为 0。28.4 节由此推出正规方程。例：\(f(c)=(c-3)^2\) 在 \(f'(c)=2(c-3)=0\) 即 \(c=3\) 处最小。见 第 00 册第 05 章 多元微积分与优化。
- 范数与条件数。 \(\|x\|\) 是向量的"长度"，\(\|x\|^2=\sum x_i^2\)。条件数 \(\kappa(A)\) 衡量"输入有一点误差，解会被放大多少倍"，\(\kappa=10^k\) 大约意味着损失 \(k\) 位有效数字。本册之外的系统讲解见第 01 册第 05b 章。
怎么读这一章。 必读：28.1.2–28.1.3（LUP 分解如何解方程）、28.1.4 的消元与舒尔补（跟着下面的数值例子走一遍）、28.1.5 选主元的动机、28.3（对称正定与舒尔补的统计解读）、28.4（正规方程及其数值隐患）、28.5.2 的三个应用。28.1.4–28.1.5 的分块推导和伪代码第一次可以只读懂一个数值例子；28.2.2"求逆与乘法同样难"的两个定理第一次可以只看结论和 28.13 式的金融用法；伪逆的四个 Moore-Penrose 条件可以跳过。
28.1 求解线性方程组
28.1.1 问题
\(n\) 个方程、\(n\) 个未知数：
若 \(A\) 非奇异，\(x=A^{-1}b\) 是唯一解。本节只讨论这种情形。方程数少于未知数或秩小于 \(n\) 的系统称为欠定（underdetermined）的，通常有无穷多解或无解；方程数多于未知数的称为超定（overdetermined）的，通常无精确解，28.4 节求它的最小二乘近似解。
不要先求逆：先算 \(A^{-1}\) 再乘 \(b\) 的方法在数值上不稳定；LUP 分解既更稳定，实践中也更快。
关于浮点误差：舍入误差可能在计算中被放大而导致错误结果，这叫数值不稳定（numerical instability）。原书只偶尔提及，详细讨论见 Golub 与 Van Loan 的《Matrix Computations》。
28.1.2 LUP 分解的思路
找 \(n\times n\) 矩阵 \(L,U,P\) 使
其中 \(L\) 是单位下三角矩阵（对角线全为 1），\(U\) 是上三角矩阵，\(P\) 是置换矩阵（permutation matrix）。每个非奇异矩阵都有 LUP 分解。有了它，\(Ax=b\) 两边左乘 \(P\)，得 \(LUx=Pb\)。令 \(y=Ux\)：
- 前代（forward substitution）：解下三角系统 \(Ly=Pb\)；
- 回代（back substitution）：解上三角系统 \(Ux=y\)。
验证：\(Ax=P^{-1}LUx=P^{-1}Ly=P^{-1}Pb=b\)。
白话解释：为什么要拆成三角矩阵？因为三角方程组"一眼就能解"。下三角例子：\(y_1=4\)；\(2y_1+y_2=10\)；\(y_1+3y_2+y_3=12\)。第一行直接给出 \(y_1=4\)，代入第二行得 \(y_2=2\)，再代入第三行得 \(y_3=2\)，每一步只有一个新未知数，这就是前代。上三角从最后一行往回解，叫回代。这样把一个"所有未知数纠缠在一起"的问题，变成两个"逐个剥开"的问题。 \(P\) 是置换矩阵，作用只是"换行顺序"：每行每列恰好一个 1，\(Pb\) 就是把 \(b\) 的分量重新排列。下面的例子里 \(Pb=(8,3,7)^T\)，就是把 \(b=(3,7,8)^T\) 的第 3 个分量放到最前面。
28.1.3 前代与回代
用数组 \(\pi[1..n]\) 紧凑地表示 \(P\)：\(P_{i,\pi[i]}=1\)，于是 \(Pb\) 的第 \(i\) 个元素是 \(b_{\pi[i]}\)。因为 \(L\) 的对角线全为 1：
def lup_solve(L, U, pi, b): # Θ(n^2) 时间
n = len(L)
y = [0.0] * n; x = [0.0] * n
for i in range(n): # 前代
y[i] = b[pi[i]] - sum(L[i][j] * y[j] for j in range(i))
for i in reversed(range(n)): # 回代
x[i] = (y[i] - sum(U[i][j] * x[j] for j in range(i + 1, n))) / U[i][i]
return x
两步各是两重循环，共 \(\Theta(n^2)\)。
例
前代解 \(Ly=Pb=(8,3,7)^T\)，得 \(y=(8,1.4,1.5)^T\)；回代解 \(Ux=y\)，得 \(x=(-1.4,2.2,0.6)^T\)。
28.1.4 LU 分解：高斯消元的递归表述
先看不选主元（\(P=I\)）的情形。高斯消元（Gaussian elimination）从其余方程中减去第一个方程的倍数以消去 \(x_1\)，再用第二个方程消去后面的 \(x_2\)……最后剩下上三角形式，即 \(U\)；消元所用的乘数组成 \(L\)。写成分块形式：
其中 \(v=(a_{21},\dots,a_{n1})^T\)，\(w^T=(a_{12},\dots,a_{1n})\)。\((n-1)\times(n-1)\) 矩阵
叫 \(A\) 关于 \(a_{11}\) 的舒尔补（Schur complement）。
推导拆解：用原书图 28.1 矩阵的左上 \(2\times2\) 块 \(\begin{pmatrix}2&3\\6&13\end{pmatrix}\) 走一遍。这里 \(a_{11}=2\)，\(w^T=(3)\)，\(v=(6)\)，\(A'=(13)\)。 第 1 步（乘数）：\(v/a_{11}=6/2=3\)，即"第二行要减去 3 倍的第一行"，这个 3 进入 \(L\) 的左下角。 第 2 步（舒尔补）：\(A'-vw^T/a_{11}=13-6\times3/2=4\)。这正是"第二行减 3 倍第一行"后第二行第二列剩下的数：\(13-3\times3=4\)，它成为下一个主元（图 28.1 的主元依次为 2、4、1、3）。 第 3 步（验证）：\(\begin{pmatrix}1&0\\3&1\end{pmatrix}\begin{pmatrix}2&3\\0&4\end{pmatrix}=\begin{pmatrix}2&3\\6&13\end{pmatrix}\)。 一般情形下 \(vw^T\) 是"列向量乘行向量"得到的矩阵（外积），\((vw^T)_{ij}=v_iw_j\)，所以 \(A'-vw^T/a_{11}\) 的每个元素就是"该位置减去 \(\frac{v_i}{a_{11}}\) 倍的第一行对应元素"，即一次消元。舒尔补就是"消掉 \(x_1\) 之后剩下的方程组的系数矩阵"。
- 若 \(A\) 非奇异，舒尔补也非奇异：否则 (28.8) 右边第二个矩阵的下面 \(n-1\) 行（第一列全为 0）的行秩小于 \(n-1\)，整个矩阵行秩小于 \(n\)，与 \(A\) 非奇异矛盾。
- 递归地分解舒尔补 \(A'-vw^T/a_{11}=L'U'\)，就得到
除数 \(a_{11}\) 和后续各步舒尔补的左上角叫主元（pivot），它们出现在 \(U\) 的对角线上。
def lu_decomposition(A): # Θ(n^3) 时间；尾递归改写为循环
n = len(A)
L = [[1.0 if i == j else 0.0 for j in range(n)] for i in range(n)]
U = [[0.0] * n for _ in range(n)]
for k in range(n):
U[k][k] = A[k][k] # 主元
for i in range(k + 1, n):
L[i][k] = A[i][k] / A[k][k] # v_i / 主元
U[k][i] = A[k][i] # w_i
for i in range(k + 1, n): # 舒尔补写回 A
for j in range(k + 1, n):
A[i][j] -= L[i][k] * U[k][j]
return L, U
三重循环，\(\Theta(n^3)\)。标准优化是把 \(L\)（严格下三角部分）和 \(U\)（上三角部分）原地存放在 \(A\) 中。
例（原书图 28.1）
主元依次为 2、4、1、3。
28.1.5 LUP 分解：部分选主元
若 \(a_{11}=0\)，或某一步舒尔补的左上角为 0，上面的方法会除以零而失败；即使不为零，除以一个很小的数也会放大舍入误差。解决办法是选主元（pivoting）：每一步把当前列中绝对值最大的元素所在行换到顶上，叫部分选主元（partial pivoting）。
推导拆解：一个经典的数值例子说明"小主元"的危害。解 \(\begin{pmatrix}10^{-20}&1\\1&1\end{pmatrix}x=\begin{pmatrix}1\\2\end{pmatrix}\)，真解非常接近 \(x=(1,1)\)。 不选主元：乘数 \(=1/10^{-20}=10^{20}\)，舒尔补 \(=1-10^{20}\)。双精度只有约 16 位有效数字，\(1-10^{20}\) 被舍入成 \(-10^{20}\)，那个"1"整个丢了。右端同样变成 \(2-10^{20}\approx-10^{20}\)，于是 \(x_2=1\)；回代 \(x_1=(1-x_2)/10^{-20}=0\)。\(x_1\) 完全算错。 选主元：先交换两行，用 1 作主元，乘数是 \(10^{-20}\)，舒尔补 \(=1-10^{-20}\approx1\)，解得 \(x_2\approx1\)、\(x_1=2-1=1\)，正确。 道理是：除以小数会产生巨大的乘数，巨大的数和普通的数相加时，普通的数在舍入中被"吞掉"。挑绝对值最大的主元，保证所有乘数的绝对值都不超过 1。
第一列不可能全为 0（否则 \(A\) 奇异）。设 \(|a_{k1}|\) 最大，交换第 1 行与第 \(k\) 行，相当于左乘置换矩阵 \(Q\)：
舒尔补非奇异，递归得 \(P'(A'-vw^T/a_{k1})=L'U'\)。令 \(P=\begin{pmatrix}1&0\\0&P'\end{pmatrix}Q\)，则
与 LU 不同，已算出的 \(L\) 的列 \(v/a_{k1}\) 也要乘 \(P'\)——所以实现时交换整行（包括已经存放 \(L\) 的部分）。
def lup_decomposition(A): # Θ(n^3)，原地；返回置换数组 pi
n = len(A)
pi = list(range(n))
for k in range(n):
p, kp = 0.0, None
for i in range(k, n): # 找第 k 列绝对值最大的元素
if abs(A[i][k]) > p:
p, kp = abs(A[i][k]), i
if p == 0:
raise ValueError("singular matrix")
pi[k], pi[kp] = pi[kp], pi[k]
A[k], A[kp] = A[kp], A[k] # 交换整行
for i in range(k + 1, n):
A[i][k] /= A[k][k] # L 的元素
for j in range(k + 1, n):
A[i][j] -= A[i][k] * A[k][j] # 舒尔补
return pi # 结束时 a_ij = l_ij (i>j)，u_ij (i<=j)
时间仍是 \(\Theta(n^3)\)——选主元只多花一个低阶项。
例（原书图 28.2）
第一步第一列最大元是第 3 行的 5，交换第 1、3 行；以下各步类推。
复用分解：分解只做一次（\(\Theta(n^3)\)），之后每个新的右端项只需 \(\Theta(n^2)\)。量化中常见的情形是同一个协方差矩阵要对多个收益预测向量求最优权重，或者同一个设计矩阵要对多个因变量做回归——都应该先分解一次再反复求解。
28.2 矩阵求逆
28.2.1 由 LUP 分解求逆
\(AX=I_n\) 可以看成 \(n\) 个方程组 \(AX_i=e_i\)（\(X_i\) 是 \(X\) 的第 \(i\) 列）。分解一次 \(\Theta(n^3)\)，每列求解 \(\Theta(n^2)\)，共 \(\Theta(n^3)\)。
实践中一般不用逆矩阵解方程，但有时确实需要逆本身——比如最小方差组合的解析式、精度矩阵（precision matrix）\(\Sigma^{-1}\) 的稀疏结构分析、协方差逆的增量更新。
28.2.2 求逆与乘法同样难
定理 28.1（乘法不比求逆难） 若能在 \(I(n)\) 时间内求 \(n\times n\) 矩阵的逆，\(I(n)=\Omega(n^2)\) 且满足 \(I(3n)=O(I(n))\)，则可在 \(O(I(n))\) 时间内计算两个 \(n\times n\) 矩阵的乘积。
证明：构造 \(3n\times3n\) 矩阵
\(AB\) 就是 \(D^{-1}\) 右上角的块。\(\square\)
白话解释：这类定理叫"归约"：如果你有一个求逆的黑箱，就能借它做乘法，所以乘法不会比求逆更难。验证 \(D\cdot D^{-1}=I\) 只需按块相乘，例如第一行块乘第三列块：\(I\cdot AB+A\cdot(-B)+0\cdot I=AB-AB=0\)，正是单位矩阵右上角应有的 0。条件 \(I(3n)=O(I(n))\) 是说规模放大 3 倍，求逆时间只放大常数倍，这样对 \(3n\) 阶矩阵求逆的代价仍可记作 \(O(I(n))\)。第一次读可以只记结论：求逆和乘法的计算难度同一个量级，都是 \(\Theta(n^3)\) 左右。
定理 28.2（求逆不比乘法难） 若能在 \(M(n)\) 时间内计算两个 \(n\times n\) 实矩阵的乘积，\(M(n)=\Omega(n^2)\)，且满足 \(M(n+k)=O(M(n))\)（\(0\le k\le n\)）与 \(M(n/2)\le cM(n)\)（某常数 \(c<1/2\)），则任何实非奇异矩阵可在 \(O(M(n))\) 时间内求逆。
证明思路：
- 补成 2 的幂阶：\(\begin{pmatrix}A&0\\0&I_k\end{pmatrix}^{-1}=\begin{pmatrix}A^{-1}&0\\0&I_k\end{pmatrix}\)。
- 先设 \(A\) 对称正定，分块
\(B\) 和 \(S\) 都对称正定（28.3 节），递归求它们的逆。共 2 次 \(n/2\) 阶递归求逆、4 次 \(n/2\) 阶乘法和 \(O(n^2)\) 的其他工作：\(I(n)\le2I(n/2)+4M(n/2)+O(n^2)=O(M(n))\)。 3. 一般非奇异 \(A\)：\(A^TA\) 对称正定，且 \(A^{-1}=(A^TA)^{-1}A^T\)。
所以用 Strassen 算法（\(O(n^{\lg7})\)）也能在 \(O(n^{\lg7})\) 内求逆——这正是 Strassen 原论文的动机。
定理 28.2 第 3 步提示了另一种解法：两边乘 \(A^T\) 得 \((A^TA)x=A^Tb\)，对对称正定的 \(A^TA\) 做不选主元的 LU 分解。理论上正确，但实践中 LUP 更好：运算量少一个常数因子，而且 \(A^TA\) 会把条件数平方——这一点在 28.4 节的最小二乘中至关重要。
分块求逆公式 (28.13) 在量化中很常用：右下块 \(S^{-1}\) 告诉我们，精度矩阵的一个对角块等于对应舒尔补的逆。28.5.2 节会用它来解释对冲后的残余风险。
28.3 对称正定矩阵
28.3.1 定义与基本性质
对称正定矩阵（symmetric positive-definite，SPD）：\(A=A^T\)，且对所有非零 \(x\)，\(x^TAx>0\)。协方差矩阵（满秩时）就是对称正定的。
金融直觉：取 \(A=\Sigma\)（协方差矩阵）、\(x=w\)（组合权重），\(x^TAx\) 就是组合方差 \(w^T\Sigma w\)。"正定"的意思是：任何不全为零的持仓组合，方差都严格大于 0，即不存在用这些资产拼出的无风险组合。若某个资产恰好是其他资产的线性组合（比如同时放进指数和它的全部成分股），就能拼出方差为 0 的组合，\(\Sigma\) 只是半正定、不可逆，最小方差权重的公式 \(\Sigma^{-1}\mathbf 1/(\mathbf 1^T\Sigma^{-1}\mathbf 1)\) 就失效了。
引理 28.3 正定矩阵非奇异。（若 \(Ax=0\) 有非零解，则 \(x^TAx=0\)。）
\(A\) 的第 \(k\) 个顺序主子矩阵（leading submatrix）\(A_k\) 是前 \(k\) 行与前 \(k\) 列的交。
引理 28.4 对称正定矩阵的每个顺序主子矩阵都对称正定。（若 \(x_k^TA_kx_k\le0\)，取 \(x=(x_k^T,0)^T\)，则 \(x^TAx\le0\)，矛盾。）
28.3.2 舒尔补引理
把 \(A\) 分块为 \(\begin{pmatrix}A_k&B^T\\B&C\end{pmatrix}\)，\(A\) 关于 \(A_k\) 的舒尔补为
引理 28.5（舒尔补引理） 若 \(A\) 对称正定，则 \(S\) 对称正定。
证明：把 \(x\) 按分块写成 \((y,z)\)，配方得
对任意 \(z\ne0\)，取 \(y=-A_k^{-1}B^Tz\)，第一项为 0，于是 \(z^TSz=x^TAx>0\)。\(\square\)
统计解读：若 \((X_1,X_2)\) 服从协方差为 \(\begin{pmatrix}\Sigma_{11}&\Sigma_{12}\\\Sigma_{21}&\Sigma_{22}\end{pmatrix}\) 的多元正态分布，则已知 \(X_1\) 时 \(X_2\) 的条件协方差是 \(\Sigma_{22}-\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}\)——正是舒尔补。配方式 (28.16) 中的 \(y=-A_k^{-1}B^Tz\) 对应最优线性预测（回归系数）。在对冲问题中，\(X_1\) 是对冲工具，\(X_2\) 是被对冲的组合，舒尔补就是用最优对冲比率对冲后剩下的风险。
金融直觉：用你在 CFA 里学过的单一期货对冲核对一下。组合波动 \(\sigma_p=0.22\)，期货波动 \(\sigma_f=0.18\)，相关系数 \(\rho=0.85\)。协方差矩阵按"期货在前、组合在后"排列，\(A_k=\sigma_f^2\)，\(B^T=\rho\sigma_p\sigma_f\)，\(C=\sigma_p^2\)。 舒尔补 \(S=\sigma_p^2-\frac{(\rho\sigma_p\sigma_f)^2}{\sigma_f^2}=\sigma_p^2(1-\rho^2)\)，对冲后波动 \(=0.22\sqrt{1-0.7225}\approx0.116\)。 而 \(A_k^{-1}B^T=\rho\sigma_p\sigma_f/\sigma_f^2=\rho\sigma_p/\sigma_f\approx1.04\)，正是 CFA 里的最小方差对冲比率。配方式 (28.16) 中让第一项为 0 的 \(y=-A_k^{-1}B^Tz\)，就是"每持有 \(z\) 单位组合，卖出 \(1.04z\) 单位期货"。28.5.2 节用两个期货对冲时，同一公式里的 \(A_k\) 换成 \(2\times2\) 矩阵，结论不变。
推论 28.6 对称正定矩阵的 LU 分解不会除以 0，而且每个主元都严格为正。第一个主元 \(a_{11}=e_1^TAe_1>0\)；LU 第一步得到关于 \(A_1=(a_{11})\) 的舒尔补，由引理 28.5 它仍对称正定，归纳即得。
所以对称正定矩阵不需要选主元。把每个主元开方分到两边，\(A=LDL^T=(LD^{1/2})(LD^{1/2})^T\)，就是 Cholesky 分解——它只需 LU 一半的运算量，是协方差矩阵的标准分解方式。Cholesky 能否成功，也是检验一个矩阵是否正定的最快方法（第 01 册第 07a 章）。
28.4 最小二乘逼近
28.4.1 正规方程
给定 \(m\) 个带测量误差的数据点 \((x_1,y_1),\dots,(x_m,y_m)\)，求函数 \(F\) 使逼近误差 \(\eta_i=F(x_i)-y_i\) 尽量小。设 \(F\) 是 \(n\) 个基函数的线性组合
常用 \(f_j(x)=x^{j-1}\)，即多项式拟合。
\(n=m\) 时可以精确穿过每个点，但高次 \(F\) 会"拟合噪声"，对新的 \(x\) 预测很差——这就是过拟合。通常取 \(n\ll m\)，希望抓住主要模式而不追逐噪声，于是得到超定方程组。
记 \(A=(a_{ij})\)，\(a_{ij}=f_j(x_i)\)（\(m\times n\)），则 \(\eta=Ac-y\)。最小化
对每个 \(c_k\) 求导并令其为 0，得 \(A^T(Ac-y)=0\)，即正规方程（normal equation）
推导拆解：把"对每个 \(c_k\) 求导"写出来。目标函数是 \(\sum_i(\sum_ja_{ij}c_j-y_i)^2\)，对 \(c_k\) 求偏导，用链式法则（外层平方的导数乘以内层对 \(c_k\) 的导数 \(a_{ik}\)）：\(\sum_i2(\sum_ja_{ij}c_j-y_i)a_{ik}=0\)。去掉 2，左边正是向量 \(A^T(Ac-y)\) 的第 \(k\) 个分量。\(k=1,\dots,n\) 全部为 0，合起来就是 \(A^T(Ac-y)=0\)。 和 CFA 的对应：只有截距和一个自变量时，\(A\) 的两列是 \((1,\dots,1)\) 和 \((x_1,\dots,x_m)\)，解这个 \(2\times2\) 的正规方程，得到的斜率就是熟悉的 \(\hat b=\frac{\sum(x_i-\bar x)(y_i-\bar y)}{\sum(x_i-\bar x)^2}=\frac{\mathrm{Cov}(x,y)}{\mathrm{Var}(x)}\)。几何上，\(A^T(Ac-y)=0\) 说的是"残差与每一个自变量都不相关"，这是 OLS 的基本性质。
若 \(A\) 列满秩，\(A^TA\) 对称正定、可逆，
\(A^+=(A^TA)^{-1}A^T\) 称为 \(A\) 的伪逆（pseudoinverse），它把逆矩阵推广到非方阵。伪逆满足四个 Moore-Penrose 条件（习题 28.3-7）：\(AA^+A=A\)，\(A^+AA^+=A^+\)，\((AA^+)^T=AA^+\)，\((A^+A)^T=A^+A\)。
例（原书图 28.3） 5 个点 \((-1,2),(1,1),(2,1),(3,0),(5,3)\)，拟合二次多项式 \(F(x)=c_1+c_2x+c_3x^2\)：
\(c=A^+y=(1.200,-0.757,0.214)^T\)，即 \(F(x)=1.200-0.757x+0.214x^2\)。
原书给出的实际做法是：算 \(A^Ty\)，对 \(A^TA\) 做分解后前代、回代。形成 \(A^TA\) 需 \(O(mn^2)\)，分解 \(O(n^3)\)。
28.4.2 正规方程的数值隐患
\(A^TA\) 的条件数是 \(A\) 的条件数的平方：\(\kappa(A^TA)=\kappa(A)^2\)。双精度浮点约有 16 位有效数字，若 \(\kappa(A)\approx10^7\)，\(\kappa(A^TA)\approx10^{14}\)，正规方程的解只剩一两位有效数字，甚至完全错误。原书注记引用 Golub 与 Van Loan：行列式不是衡量稳定性的好指标，应看条件数 \(\|A\|\,\|A^{-1}\|\)。
更稳定的做法：
- QR 分解：\(A=QR\)，\(Q\) 列正交、\(R\) 上三角，则 \(Rc=Q^Ty\)。只涉及 \(\kappa(A)\)，不平方。
- 奇异值分解（SVD）：\(A=Q_1\Sigma Q_2^T\)（原书注记），对极小的奇异值截断，可以稳健地处理近似秩亏。
numpy.linalg.lstsq用的就是它。
28.5.1 节用因子回归演示这一点。
白话解释：条件数可以理解为"误差放大倍数"。\(\kappa=10^7\) 时，数据或舍入中的相对误差 \(10^{-16}\) 可能被放大为 \(10^{-9}\)，结果仍有约 9 位可信数字；但正规方程先把矩阵变成 \(A^TA\)，放大倍数平方成 \(10^{14}\)，结果只剩约 2 位可信数字。在因子回归里，两个高度相关的因子（例如两种定义略有不同的动量）就会让 \(\kappa\) 变大。 原文这里说"28.5.1 节"用因子回归演示，实际演示在 28.5.2 节的"应用一：因子回归与条件数"，28.5.1 节只是核对原书例题。
28.5 量化实战
28.5.1 先核对原书例题
import numpy as np
import scipy.linalg as sla
def lup_decomposition(A):
"""原地 LUP 分解（部分选主元）。返回 (LU 合并存储, pi)，满足 A[pi] = L @ U。"""
A = np.array(A, dtype=float); n = len(A); pi = np.arange(n)
for k in range(n):
kp = k + np.argmax(np.abs(A[k:, k])) # 第 k 列绝对值最大的元素
if A[kp, k] == 0:
raise ValueError("singular matrix")
A[[k, kp]] = A[[kp, k]]; pi[[k, kp]] = pi[[kp, k]] # 整行交换
A[k+1:, k] /= A[k, k] # L 的第 k 列
A[k+1:, k+1:] -= np.outer(A[k+1:, k], A[k, k+1:]) # 舒尔补
return A, pi
def lup_solve(LU, pi, b):
n = len(LU); y = np.zeros(n); x = np.zeros(n)
for i in range(n): # 前代：L y = P b
y[i] = b[pi[i]] - LU[i, :i] @ y[:i]
for i in reversed(range(n)): # 回代：U x = y
x[i] = (y[i] - LU[i, i+1:] @ x[i+1:]) / LU[i, i]
return x
def split(LU):
return np.tril(LU, -1) + np.eye(len(LU)), np.triu(LU)
# 28.1 节例题：A x = b
A = np.array([[1, 2, 0], [3, 4, 4], [5, 6, 3]]); b = np.array([3, 7, 8])
LU, pi = lup_decomposition(A); L, U = split(LU)
print("pi =", pi, "\nL =\n", L, "\nU =\n", U)
print("x =", lup_solve(LU, pi, b))
# 图 28.2：4x4 的 LUP 分解
A2 = np.array([[2, 0, 2, 0.6], [3, 3, 4, -2], [5, 5, 4, 2], [-1, -2, 3.4, -1]])
LU2, pi2 = lup_decomposition(A2); L2, U2 = split(LU2)
print("图28.2 pi =", pi2, "\nU =\n", U2.round(3))
print("与 scipy.linalg.lu_factor 的 U 一致:", np.allclose(np.triu(sla.lu_factor(A2)[0]), U2))
# 28.2 节：用 n 次 LUP-SOLVE 求逆
Ainv = np.column_stack([lup_solve(LU2, pi2, e) for e in np.eye(4)])
print("A @ A^{-1} = I:", np.allclose(A2 @ Ainv, np.eye(4)))
# 图 28.3：最小二乘拟合二次多项式
xs = np.array([-1, 1, 2, 3, 5.]); ys = np.array([2, 1, 1, 0, 3.])
X = np.column_stack([np.ones(5), xs, xs ** 2])
Aplus = np.linalg.solve(X.T @ X, X.T) # 伪逆 (A^T A)^{-1} A^T
print("A+ =\n", Aplus.round(3))
print("c = A+ y =", (Aplus @ ys).round(3))
输出：
pi = [2 0 1]
L =
[[1. 0. 0. ]
[0.2 1. 0. ]
[0.6 0.5 1. ]]
U =
[[ 5. 6. 3. ]
[ 0. 0.8 -0.6]
[ 0. 0. 2.5]]
x = [-1.4 2.2 0.6]
图28.2 pi = [2 0 3 1]
U =
[[ 5. 5. 4. 2. ]
[ 0. -2. 0.4 -0.2]
[ 0. 0. 4. -0.5]
[ 0. 0. 0. -3. ]]
与 scipy.linalg.lu_factor 的 U 一致: True
A @ A^{-1} = I: True
A+ =
[[ 0.5 0.3 0.2 0.1 -0.1 ]
[-0.388 0.093 0.19 0.193 -0.088]
[ 0.06 -0.036 -0.048 -0.036 0.06 ]]
c = A+ y = [ 1.2 -0.757 0.214]
\(L,U,x\)、图 28.2 的 \(U\) 与置换（\(\pi=[2,0,3,1]\) 即 \(P\) 的第 1 行取 \(A\) 的第 3 行，依此类推，下标从 0 起）、伪逆和拟合系数都与原书一致。代码把原书伪代码中最内两层循环换成了 np.outer 的秩 1 更新，这正是"每一步对舒尔补做秩 1 修正"的向量化写法。
28.5.2 三个量化应用
import numpy as np
import scipy.linalg as sla
from scipy.interpolate import CubicSpline
rng = np.random.default_rng(28)
# ---------- 1. 截面回归：正规方程 vs QR，共线性如何放大误差 ----------
n, k = 3000, 6
F = rng.normal(size=(n, k)) # 6 个风格因子暴露（标准化）
beta_true = np.array([0.8, -0.5, 0.3, 0.0, 0.2, -0.1]) / 100
for eps in (1e-1, 1e-4, 1e-7):
X = np.column_stack([np.ones(n), F, F[:, 0] + eps * rng.normal(size=n)]) # 第 8 列几乎等于"因子1"
b_true = np.r_[0.0, beta_true, 0.004]
y = X @ b_true + 0.02 * rng.normal(size=n)
b_ne = sla.solve(X.T @ X, X.T @ y, assume_a="pos") # 正规方程 + Cholesky
Q, R = np.linalg.qr(X); b_qr = sla.solve_triangular(R, Q.T @ y)
b_ls = np.linalg.lstsq(X, y, rcond=None)[0] # 基于 SVD
print(f"eps={eps:.0e}: cond(X)={np.linalg.cond(X):.1e}, cond(X'X)={np.linalg.cond(X.T @ X):.1e}, "
f"正规方程与 QR 解的差 {np.abs(b_ne - b_qr).max():.1e}, QR 与 SVD 的差 {np.abs(b_qr - b_ls).max():.1e}")
# ---------- 2. 舒尔补 = 条件协方差 = 对冲后的残余风险 ----------
# 资产 0 是要对冲的股票组合，资产 1、2 是两个股指期货
vol = np.array([0.22, 0.18, 0.20]); corr = np.array([[1, 0.85, 0.80], [0.85, 1, 0.90], [0.80, 0.90, 1]])
S = corr * np.outer(vol, vol)
Spp, Sph, Shh = S[0, 0], S[0, 1:], S[1:, 1:]
h = sla.solve(Shh, Sph, assume_a="pos") # 最小方差对冲比率 Σ_hh^{-1} Σ_hp
resid_var = Spp - Sph @ h # 舒尔补 S_pp - S_ph S_hh^{-1} S_hp
print("对冲比率（每 1 元组合卖出期货）:", h.round(3))
print(f"对冲前波动 {np.sqrt(Spp):.3f}，对冲后波动 {np.sqrt(resid_var):.3f}")
# 验证：舒尔补等于 (S^{-1})_{00} 的倒数（分块求逆公式 28.13）
print("舒尔补 = 1/(Σ^{-1})_00:", np.isclose(resid_var, 1 / np.linalg.inv(S)[0, 0]))
# 模拟验证：用 Cholesky 因子生成相关收益，看对冲后的样本波动
Lc = np.linalg.cholesky(S / 252)
r = rng.normal(size=(252 * 20, 3)) @ Lc.T
print(f"模拟 20 年日收益：对冲后年化波动 {np.std(r[:, 0] - r[:, 1:] @ h) * np.sqrt(252):.3f}")
print("Cholesky 的对角元（LU 主元的平方根）全为正:", np.all(np.diag(Lc) > 0))
# ---------- 3. 思考题 28-1/28-2：三对角系统与自然三次样条收益率曲线 ----------
def thomas(a, b, c, d):
"""解三对角系统：a 为次对角线（长 n-1），b 主对角线（长 n），c 超对角线（长 n-1）。O(n)。"""
n = len(b); cp = np.zeros(n - 1); dp = np.zeros(n)
cp[0] = c[0] / b[0]; dp[0] = d[0] / b[0]
for i in range(1, n): # 前向消元 = 不选主元的 LU 分解
m = b[i] - a[i - 1] * cp[i - 1]
if i < n - 1: cp[i] = c[i] / m
dp[i] = (d[i] - a[i - 1] * dp[i - 1]) / m
x = np.zeros(n); x[-1] = dp[-1]
for i in range(n - 2, -1, -1): # 回代
x[i] = dp[i] - cp[i] * x[i + 1]
return x
T = np.array([0.25, 0.5, 1, 2, 3, 5, 7, 10, 20, 30]) # 期限（年）
Y = np.array([4.10, 4.05, 3.90, 3.70, 3.62, 3.60, 3.68, 3.80, 4.15, 4.20]) # 即期收益率（%）
hh = np.diff(T)
# 自然样条：二阶导 M_0 = M_n = 0，内部 M_i 满足对称、严格对角占优的三对角方程
rhs = 6 * (np.diff(Y[1:]) / hh[1:] - np.diff(Y[:-1]) / hh[:-1])
M = np.r_[0, thomas(hh[1:-1], 2 * (hh[:-1] + hh[1:]), hh[1:-1], rhs), 0]
def spline_eval(t):
i = np.clip(np.searchsorted(T, t) - 1, 0, len(T) - 2); h = hh[i]
A, B = (T[i + 1] - t) / h, (t - T[i]) / h
return A * Y[i] + B * Y[i + 1] + ((A ** 3 - A) * M[i] + (B ** 3 - B) * M[i + 1]) * h ** 2 / 6
grid = np.array([1.5, 4, 6, 8.5, 15, 25])
ours = spline_eval(grid); ref = CubicSpline(T, Y, bc_type="natural")(grid)
print("插值收益率(%):", dict(zip(grid.tolist(), ours.round(4).tolist())))
print("与 scipy CubicSpline(natural) 一致:", np.allclose(ours, ref))
# 三对角求解的规模：O(n) 的 Thomas vs O(n^3) 的稠密 LU
import time
for n_ in (500, 2000):
a_ = -np.ones(n_ - 1); b_ = 2.5 * np.ones(n_); d_ = rng.normal(size=n_)
t0 = time.perf_counter(); x1 = thomas(a_, b_, a_, d_); t1 = time.perf_counter()
Afull = np.diag(b_) + np.diag(a_, 1) + np.diag(a_, -1); x2 = sla.lu_solve(sla.lu_factor(Afull), d_); t2 = time.perf_counter()
print(f"n={n_}: Thomas {1e3*(t1-t0):.2f} ms，稠密 LU {1e3*(t2-t1):.2f} ms，结果一致 {np.allclose(x1, x2)}")
输出（计时随机器而异）：
eps=1e-01: cond(X)=2.1e+01, cond(X'X)=4.3e+02, 正规方程与 QR 解的差 1.3e-15, QR 与 SVD 的差 4.3e-18
eps=1e-04: cond(X)=2.0e+04, cond(X'X)=4.0e+08, 正规方程与 QR 解的差 2.0e-07, QR 与 SVD 的差 3.9e-15
eps=1e-07: cond(X)=2.0e+07, cond(X'X)=5.4e+14, 正规方程与 QR 解的差 1.9e+02, QR 与 SVD 的差 3.3e-12
对冲比率（每 1 元组合卖出期货）: [0.836 0.203]
对冲前波动 0.220，对冲后波动 0.115
舒尔补 = 1/(Σ^{-1})_00: True
模拟 20 年日收益：对冲后年化波动 0.114
Cholesky 的对角元（LU 主元的平方根）全为正: True
插值收益率(%): {1.5: 3.778, 4.0: 3.5923, 6.0: 3.6344, 8.5: 3.7426, 15.0: 3.9965, 25.0: 4.2053}
与 scipy CubicSpline(natural) 一致: True
n=500: Thomas 0.25 ms，稠密 LU 1.34 ms，结果一致 True
n=2000: Thomas 0.97 ms，稠密 LU 31.38 ms，结果一致 True
应用一：因子回归与条件数。截面回归 \(r=X\beta+\varepsilon\) 是 Barra 类风险模型和 Fama-MacBeth 回归的核心计算。实验中第 8 列因子与"因子 1"几乎重合（相当于两个高度相关的风格因子，比如两种定义略有不同的动量）。\(\kappa(X)\) 从 21 升到 \(2\times10^7\) 时，\(\kappa(X^TX)\) 恰好是它的平方，到 \(5\times10^{14}\)；此时正规方程的解与 QR 的解相差 190（系数的量级是 \(10^{-2}\)），完全不可用，而 QR 与 SVD 仍一致到 \(10^{-12}\)。
需要分清两个问题：数值问题是"算出来的是不是最小二乘解"，QR/SVD 能解决；统计问题是"这个最小二乘解有没有意义"——两列几乎相同时，它们各自的系数方差极大，即使算得准也不可信（多重共线性，第 05 册）。工程上的对策是：用 QR/SVD 保证算对；用因子正交化、合并或岭回归（第 03 册）保证有意义。
应用二：舒尔补与对冲。要用两个股指期货对冲一个股票组合，最小方差对冲比率是 \(h=\Sigma_{hh}^{-1}\Sigma_{hp}\)（这是一个 2×2 的对称正定方程组，直接用 Cholesky 求解，不求逆），对冲后的残余方差是舒尔补 \(\Sigma_{pp}-\Sigma_{ph}\Sigma_{hh}^{-1}\Sigma_{hp}\)。本例对冲把年化波动从 22% 降到 11.5%。代码还验证了两个联系：残余方差等于精度矩阵对应对角元的倒数（分块求逆公式 28.13）；用 Cholesky 因子 \(L\)（\(\Sigma=LL^T\)）把独立正态变换为相关收益，模拟得到的对冲后波动 11.4% 与理论值一致。Cholesky 对角元全为正，即推论 28.6"对称正定矩阵的主元全正"。
应用三：三对角系统与收益率曲线。思考题 28-2 的自然三次样条在节点处二阶导连续、两端二阶导为 0，求各节点二阶导 \(M_i\) 归结为一个对称、三对角、严格对角占优（因而正定）的线性方程组：
其中 \(h_i=x_{i+1}-x_i\)。原书思考题用等距节点和一阶导 \(D_i\) 写出等价的方程 \(D_{i-1}+4D_i+D_{i+1}=3(y_{i+1}-y_{i-1})\)；这里用二阶导形式是为了直接处理不等距的期限节点（思考题 28-2(f)）。三对角系统用不选主元的 LU 分解（Thomas 算法，思考题 28-1）只需 \(O(n)\)；先求逆则至少 \(\Omega(n^2)\)，因为三对角矩阵的逆一般是稠密的。实测 \(n=2000\) 时 Thomas 算法比稠密 LU 快 30 倍以上，差距随 \(n\) 线性扩大。
同一个三对角求解器也是期权定价中有限差分法的核心：隐式格式和 Crank-Nicolson 格式（第 08 册）每个时间步都要解一个三对角系统。
28.5.3 工程要点汇总
- 永远不要写
inv(A) @ b。用solve；对称正定矩阵用assume_a="pos"或cho_factor/cho_solve。 - 分解一次，求解多次。
lu_factor/cho_factor返回的分解可以对多个右端项复用。 - 回归用 QR 或
lstsq，不要手写正规方程，除非确认 \(\kappa(X)\) 很小。 - 协方差矩阵先用 Cholesky 检查正定性。样本协方差在资产数接近样本数时会接近奇异，需要收缩估计（第 03、06 册）。
- 利用结构：三对角、带状、稀疏矩阵都有远快于 \(\Theta(n^3)\) 的专用解法。
本章小结
解 \(Ax=b\) 不应先求逆，而应做 LUP 分解 \(PA=LU\)（\(\Theta(n^3)\)），再用前代和回代（各 \(\Theta(n^2)\)）求解，分解一次可复用于多个右端项。LU 分解就是高斯消元，每一步把问题化为对舒尔补的分解；部分选主元既避免除以零，又改善数值稳定性。求逆可由 \(n\) 次 LUP-SOLVE 完成；理论上求逆与矩阵乘法同样难，可借 Strassen 加速。对称正定矩阵的顺序主子矩阵和舒尔补仍对称正定，LU 主元全为正，无需选主元，并导出 Cholesky 分解。最小二乘拟合归结为正规方程 \(A^TAc=A^Ty\)，解为伪逆 \(c=A^+y\)，但正规方程会把条件数平方，共线时应改用 QR 或 SVD。量化中，舒尔补就是条件协方差和对冲后的残余风险；三对角系统让样条插值和有限差分定价在 \(O(n)\) 内完成。
| 概念 / 结论 | 公式或要点 |
|---|---|
| LUP 分解 | \(PA=LU\)，\(\Theta(n^3)\) |
| 求解 | 前代 \(Ly=Pb\)，回代 \(Ux=y\)，各 \(\Theta(n^2)\) |
| 消元一步 | \(A=\begin{pmatrix}1&0\\v/a_{11}&I\end{pmatrix}\begin{pmatrix}a_{11}&w^T\\0&A'-vw^T/a_{11}\end{pmatrix}\) |
| 部分选主元 | 每步取当前列绝对值最大的元素作主元，整行交换 |
| 求逆 | \(n\) 次 LUP-SOLVE，\(\Theta(n^3)\)；与乘法同阶难度 |
| 分块求逆 | 右下块 \(=S^{-1}\)，\(S=D-CB^{-1}C^T\) |
| 舒尔补引理 | \(A\) 对称正定 ⇒ \(S=C-BA_k^{-1}B^T\) 对称正定 |
| SPD 的 LU | 主元全正，无需选主元；\(A=LL^T\)（Cholesky） |
| 正规方程 | \(A^TAc=A^Ty\)，\(c=A^+y\)，\(A^+=(A^TA)^{-1}A^T\) |
| 条件数 | \(\kappa(A^TA)=\kappa(A)^2\)，共线时改用 QR/SVD |
| 最小方差对冲 | \(h=\Sigma_{hh}^{-1}\Sigma_{hp}\)，残余方差 = 舒尔补 |
| 三对角系统 | Thomas 算法 \(O(n)\)；自然三次样条、Crank-Nicolson |
练习
基础
- 用前代求解 \(\begin{pmatrix}1&0&0\\4&1&0\\-6&5&1\end{pmatrix}x=\begin{pmatrix}3\\14\\-7\end{pmatrix}\)。（原书 28.1-1。答：\(x=(3,2,1)^T\)。）
- 求 \(\begin{pmatrix}4&-5&6\\8&-6&7\\12&-7&12\end{pmatrix}\) 的 LU 分解。（原书 28.1-2。可用 28.5.1 的代码核对，注意那里用的是 LUP。）
- 用 LUP 分解求解 \(\begin{pmatrix}1&5&4\\2&0&3\\5&8&2\end{pmatrix}x=\begin{pmatrix}12\\9\\5\end{pmatrix}\)。（原书 28.1-3。）
- 证明对称正定矩阵的对角元都为正，且最大元素在对角线上。（原书 28.3-1、28.3-3。）
- 证明 LU 分解的第 \(k\) 个主元等于 \(\det(A_k)/\det(A_{k-1})\)（\(\det A_0=1\)）。（原书 28.3-5。）
- 用 \(F(x)=c_1+c_2x\lg x+c_3e^x\) 对点 \((1,1),(2,1),(3,3),(4,8)\) 做最小二乘拟合。（原书 28.3-6。提示：设计矩阵的三列是 \(1\)、\(x\lg x\)、\(e^x\)。）
进阶
- 证明伪逆满足四个 Moore-Penrose 条件。（原书 28.3-7。）
- 思考题 28-1：对 5×5 三对角矩阵（主对角线 \(1,2,2,2,2\)，次对角线全为 \(-1\)）求 LU 分解，解 \(Ax=(1,1,1,1,1)^T\)，并求逆。说明为什么任何先求 \(A^{-1}\) 的方法最坏情况下都比 \(O(n)\) 慢。
- 在 28.5.2 的应用一中，加入岭回归 \(\min\|Xc-y\|^2+\lambda\|c\|^2\)。证明它的正规方程是 \((X^TX+\lambda I)c=X^Ty\)，并说明 \(\lambda>0\) 如何改善条件数。（提示：\(X^TX+\lambda I\) 的特征值是 \(\sigma_i^2+\lambda\)。）
- 最小方差组合 \(\min w^T\Sigma w\)，s.t. \(\mathbf 1^Tw=1\)，解为 \(w=\Sigma^{-1}\mathbf 1/(\mathbf 1^T\Sigma^{-1}\mathbf 1)\)。写一个只做一次 Cholesky 分解、不显式求逆的实现，并与
np.linalg.inv版本在 \(\kappa(\Sigma)=10^{12}\) 时比较精度。
原书推荐习题：28.1-3、28.1-7、28.2-3、28.3-2、28.3-5、28.3-6、28.3-7，思考题 28-1（PDE 定价必备）、28-2（曲线插值必备）。
原书对照
页码换算：原书页码 = PDF 页码 − 21。
| 本章小节 | 原书章节 | PDF 页码（原书页码） |
|---|---|---|
| 28.1 求解线性方程组 | 第 28 章导言；28.1 Solving systems of linear equations | p.834–848（813–827） |
| 28.2 矩阵求逆 | 28.2 Inverting matrices | p.848–853（827–832） |
| 28.3 对称正定矩阵、28.4 最小二乘逼近 | 28.3 Symmetric positive-definite matrices and least-squares approximation | p.853–861（832–840） |
| — | Problems 28-1、28-2；Chapter notes | p.861–863（840–842） |