# 矩阵函数单元题解：把函数值、导数与分支分别验清

## 任务与运行方式

研究矩阵族

$$
T_\varepsilon=\begin{pmatrix}1&2&1\\0&1+\varepsilon&3\\0&0&4\end{pmatrix},
\qquad \varepsilon\in\{0.1,10^{-6},10^{-12},0\}.
$$

要求计算指数、主平方根和主对数，给出缺陷极限的精确答案，计算非交换方向 $E=e_2e_1^T$ 的指数 Fréchet 导数，并比较完整指数与只求 $e^{T_\varepsilon}e_3$ 的路线。最后将相同的验错原则用于 Sylvester 方程、四节点扩散和极分解。

配套的 `foundation-matrix-functions-capstone.py` 是完整独立脚本，不读取本地题库或研究资料。它依赖 Python 3、NumPy、SciPy、mpmath；本次记录环境为 Python 3.12.14、NumPy 2.3.5、SciPy 1.17.0、mpmath 1.3.0。运行：

```sh
python foundation-matrix-functions-capstone.py --output matrix-functions-results.json
```

脚本内有独立写出的三角 Sylvester、固定十三阶 Padé、局部 Taylor 加块 Parlett、三角平方根、逆缩放对数、Denman–Beavers、Newton 极因子与 Arnoldi 指数作用核心。Schur 分解、线性求解和稠密乘法使用库例程；SciPy 专用函数及 mpmath 八十位计算作为交叉参照。参照计算不计入教学算法的操作计数。末位舍入可能随平台变化，因此接受标准使用范数容差，而不要求浮点结果的文本逐字一致。

## 一、先得到不能被数值 Jordan 分解替代的精确极限

当 $\varepsilon=0$，极小多项式为 $(t-1)^2(t-4)$。记 $a=f(1)$、$b=f'(1)$、$c=f(4)$，Hermite 多项式为

$$
p(t)=a+b(t-1)+\frac{c-a-3b}{9}(t-1)^2.
$$

置 $M=T_0-I$，直接乘法给出

$$
M=\begin{pmatrix}0&2&1\\0&0&3\\0&0&3\end{pmatrix},\qquad
M^2=\begin{pmatrix}0&0&9\\0&0&9\\0&0&9\end{pmatrix}.
$$

所以

$$
f(T_0)=\begin{pmatrix}a&2b&c-a-2b\\0&a&c-a\\0&0&c\end{pmatrix}.
$$

分别代入指数、主平方根和主对数：

$$
F_0=e^{T_0}=\begin{pmatrix}e&2e&e^4-3e\\0&e&e^4-e\\0&0&e^4\end{pmatrix},
$$

$$
S_0=T_0^{1/2}=\begin{pmatrix}1&1&0\\0&1&1\\0&0&2\end{pmatrix},
\qquad
L_0=\operatorname{Log}T_0=
\begin{pmatrix}0&2&\log4-2\\0&0&\log4\\0&0&\log4\end{pmatrix}.
$$

乘法可立即核对 $S_0^2=T_0$。主根的三个特征值是 $1,1,2$，都在开右半平面；主对数的三个特征值是 $0,0,\log4$，都在主值条带内。这些条件连同残差分别验证“算对方程”和“选对分支”。

这里的 Jordan 极限用于推导参照，不要求数值求 Jordan 形。事实上，接近缺陷矩阵时强行提取数值 Jordan 块会把一个连续求值问题改造成不稳定的结构识别问题。

## 二、接近而不相等的特征值怎样计算

### 指数：跨块递推，无小分母

把指标 $\{1,2\}$ 合成一个二阶对角块，第三个指标单独成块。令 $a=e$、$b=e^{1+\varepsilon}$、$c=e^4$。当 $\varepsilon\ne0$ 时，精确表达为

$$
F_{12}=2e\frac{\operatorname{expm1}(\varepsilon)}{\varepsilon},\qquad
F_{23}=\frac{3(c-b)}{3-\varepsilon},
$$

$$
F_{13}=\frac{c-a-3F_{12}+2F_{23}}3.
$$

第一项在 $\varepsilon=0$ 连续延拓为 $2e$。脚本不依赖这个专门公式，而在二阶块上围绕平均对角值作 Taylor 展开，然后解一次 $2\times2$ 三角 Sylvester 耦合。截断由显式余项界控制。所选保守界使这四个输入都使用二十三次二阶块乘法；在 $\varepsilon=0$，知道幂零性的人可以直接手算两项，这与通用保守停止器的计数并不矛盾。

对 $\varepsilon=10^{-12}$，直接相减的双精度差商给出约 $5.43694494$，而高精度的 $F_{12}$ 约为 $5.436563656920809$。脚本的参照使用实际存入的浮点矩阵作为八十位输入，因此没有把十进制输入转换误差混入算法误差。

完整矩阵也可用固定十三阶 Padé 计算。所有四个矩阵的 $1$-范数都是 $8$，2005 范数阈值规则给出 $s=1$，所以每次需七次 $3\times3$ 矩阵乘法和一次多右端求解。分块法与 Padé 法的结果都同八十位指数比较。

### 平方根：递推的分母是根之和

记 $r=\sqrt{1+\varepsilon}$。三角主根的对角是 $(1,r,2)$，非对角为

$$
S_{12}=\frac2{1+r},\qquad
S_{23}=\frac3{r+2},\qquad
S_{13}=\frac{1-S_{12}S_{23}}3.
$$

每个分母均远离零。$S_{13}$ 在极限恰好为零，因此应以整个矩阵的范数误差和平方残差接受结果，不能给这个零元素规定通常的相对误差。

### 对数：主根使输入靠近单位阵

精确对数的对角是 $(0,\log(1+\varepsilon),\log4)$，非对角为

$$
(L_\varepsilon)_{12}=2\frac{\log(1+\varepsilon)}{\varepsilon},\qquad
(L_\varepsilon)_{23}=\frac{3[\log4-\log(1+\varepsilon)]}{3-\varepsilon},
$$

$$
(L_\varepsilon)_{13}=\frac{\log4-3(L_\varepsilon)_{12}+2(L_\varepsilon)_{23}}3.
$$

第一项应使用 `log1p` 型公式，并在零处取极限 $2$。脚本的算法路线则为四次主平方根，使 $X=T_\varepsilon^{1/16}-I$ 的 $1$-范数不超过 $1/4$，然后用十点 Gauss–Legendre 部分分式求值。四个输入的精确算术截断上界均小于 $4.11\times10^{-13}$。这不包含开方和求解舍入，所以还需与八十位参照及 $e^{L_\varepsilon}-T_\varepsilon$ 交叉核验。

## 三、非交换方向的导数也有精确答案

取 $E=e_2e_1^T$，即只有 $(2,1)$ 元为一。它与 $T_0$ 不交换。块指数恒等式给出

$$
\exp\begin{pmatrix}T_0&E\\0&T_0\end{pmatrix}
=\begin{pmatrix}F_0&\mathcal L\\0&F_0\end{pmatrix},
\qquad \mathcal L=L_{\exp}(T_0,E).
$$

不依赖数值块指数，也能手推 $\mathcal L$。写 $T_0=\begin{pmatrix}J&u\\0&4\end{pmatrix}$，其中

$$
J=I+N,\quad N=\begin{pmatrix}0&2\\0&0\end{pmatrix},\quad
u=\begin{pmatrix}1\\3\end{pmatrix},\quad
E_2=\begin{pmatrix}0&0\\1&0\end{pmatrix}.
$$

$N^2=0$，但 $NE_2N$ 不为零。指数导数的左上块为

$$
e\left[E_2+\frac{NE_2+E_2N}{2}+\frac{NE_2N}{6}\right]
=e\begin{pmatrix}1&2/3\\1&1\end{pmatrix}.
$$

右上耦合可写为 $h(J)u$，其中 $h(z)=(e^4-e^z)/(4-z)$。在 $z=1$ 处，

$$
h'(1)=\frac{e^4-4e}{9},\quad
h''(1)=\frac{2e^4-17e}{27},\quad
h'''(1)=\frac{2e^4-26e}{27}.
$$

于是 $L_h(J,E_2)u$ 给出右上两项，最终

$$
\mathcal L=
\begin{pmatrix}
e&2e/3&(2e^4-23e)/9\\
e&e&(e^4-7e)/3\\
0&0&0
\end{pmatrix}
\approx
\begin{pmatrix}
2.71828183&1.81218789&5.18620200\\
2.71828183&2.71828183&11.85672574\\
0&0&0
\end{pmatrix}.
$$

脚本用独立十三阶核心计算 $6\times6$ 块指数，再与上述公式及八十位块指数比较。中心差分提供第三种检查。以 $\varepsilon=0$ 为例，Frobenius 误差记录为：

| $h$ | 中心差分误差 |
| --- | --- |
| $10^{-2}$ | $1.10836\times10^{-5}$ |
| $10^{-3}$ | $1.10834\times10^{-7}$ |
| $10^{-4}$ | $1.15291\times10^{-9}$ |
| $10^{-5}$ | $2.58818\times10^{-10}$ |
| $10^{-6}$ | $8.62762\times10^{-9}$ |
| $10^{-7}$ | $1.02361\times10^{-7}$ |

前几行体现二阶截断误差，最后几行体现相减和函数求值舍入。把 $h$ 一直缩小，并不会一直提高导数精度。

## 四、残差需要怎样的敏感性解释

对 Sylvester 方程 $AX-XB=C$，解误差至多为残差范数除以 $\operatorname{sep}(A,B)$。三角测试系统使用

$$
T=\begin{pmatrix}1&2&-1\\0&2&3\\0&0&4\end{pmatrix},\quad
S=\begin{pmatrix}-1&1&2\\0&-2&-1\\0&0&-3\end{pmatrix},\quad
D=\begin{pmatrix}-2&7&5\\3&5&10\\10&-2&-11\end{pmatrix}.
$$

从左到右的三个有效右端分别为 $(-2,3,10)^T$、$(8,4,0)^T$、$(5,7,-7)^T$；每列从下到上求解后得到

$$
Y=\begin{pmatrix}1&2&0\\-1&1&2\\2&0&-1\end{pmatrix}.
$$

脚本还用一个由 $3/5,4/5$ 构成的正交旋转和一个置换矩阵把这个问题换回非三角原坐标，执行 Schur 求解并核对原残差。实二阶块另用 $A=\begin{pmatrix}0&-1\\1&0\end{pmatrix}$、$B=(2)$、$C=(1,0)^T$，得到 $X=(-2/5,-1/5)^T$；实际 LAPACK 返回的 scale 与 info 一起记录。

谱间距不能替代 sep。对 $A_M=\begin{pmatrix}1&M\\0&2\end{pmatrix}$、$B=0$、$M\in\mathbb R$，最近特征值距离始终为一，但 $M=100$ 时 sep 约为 $0.019995$，单位右端 $(0,1)^T$ 的解为 $(-50,1/2)^T$。

平方根的局部算子是 $E\mapsto S_0E+ES_0$。本例 $\operatorname{sep}(S_0,-S_0)\approx1.02004342$。若专门构造 $\widehat S=S_0+10^{-6}E$，由于 $E^2=0$，其平方残差恰为 $10^{-6}(S_0E+ES_0)$。误差范数为 $10^{-6}$，残差约为 $2.44948974\times10^{-6}$，除以 sep 给出约 $2.40135830\times10^{-6}$，覆盖真实误差。这里没有遗漏二次项，是因为选定方向满足 $E^2=0$；对一般近似根要保留二次项或使用局部估计的适用条件。

即使平方残差为零，$-S_0$ 仍不是主根。即使指数重构残差为零，$\operatorname{Log}(e^{2\pi J})=0$ 仍不等于 $2\pi J$。这两项反例检查分支，不能被残差测试替代。

## 五、只输出向量时改变成本，也改变验算方法

对本题三阶矩阵及 $b=e_3$，最多三步 Arnoldi 便得到完整不变空间，精确运算中的 Krylov 指数作用就是 $F_\varepsilon e_3$。脚本报告三次矩阵向量乘，并记录投影矩阵和基。小型指数仍有自身的成本；在这个很小的例子中，Krylov 并不必然比完整指数更快。它的优势面向 $m\ll n$ 且 $A$ 只适合做稀疏矩阵向量乘的情形。

用四节点扩散矩阵 $A=\operatorname{tridiag}(1,-2,1)$、$b=e_1$、$t=1$ 检查非平凡截断。对 $m=1,2,3$，实际误差分别为 $0.22194570,0.09434187,0.02697276$，终点缺陷分别为 $0.13533528,0.15904619,0.07972490$。

这里 $A$ 对称负定，$\|e^{\tau A}\|_2\le1$ 对 $\tau\ge0$ 成立。由于投影指数逐元素非负，缺陷积分可计算为

$$
B_m=e_m^TH_m^{-1}(e^{H_m}-I)e_1,
$$

得到 $B_1\approx0.43233236$、$B_2\approx0.15769146$、$B_3\approx0.04385172$。每项都覆盖真实误差。对本题的 $T_\varepsilon$，没有扩散收缩性，不能直接套这三个界；一般情况应从

$$
y(t)-y_m(t)=-\int_0^t e^{(t-s)A}d_m(s)\,ds
$$

出发，另行控制传播因子。终点缺陷不构成通用误差证书。

## 六、两种迭代各自要检查什么

对 $A=\begin{pmatrix}4&6\\0&9\end{pmatrix}$，Denman–Beavers 两轮得到

$$
Y_2=\begin{pmatrix}41/20&81/50\\0&17/5\end{pmatrix},\qquad
Z_2=\begin{pmatrix}41/80&-97/600\\0&17/45\end{pmatrix}.
$$

检查 $Y_k=AZ_k$、$Y_k^2-A$、$Y_kZ_k-I$，并最终检查 $Y_k$ 的谱实部为正。脚本固定运行七轮，计入十四次多右端求解。每轮两个更新都使用旧状态，不能把新 $Y$ 提前用于更新 $Z$。

对 $A=\begin{pmatrix}2&1\\0&1\end{pmatrix}$，Newton 极因子迭代最终得到

$$
U=\frac1{\sqrt{10}}\begin{pmatrix}3&1\\-1&3\end{pmatrix},\qquad
H=\frac1{\sqrt{10}}\begin{pmatrix}6&2\\2&4\end{pmatrix}.
$$

检查 $U^TU=I$、$UH=A$、$H=H^T\succ0$，还可比较到 $A$ 的 Frobenius 距离：极因子约为 $1.29439$，QR 的 $Q=I$ 则为 $\sqrt2$。七轮迭代计入七次多右端求解。

反射例 $A=\operatorname{diag}(-2,1)$ 的普通极因子行列式为 $-1$。若题目限制 $\det U=1$，最近旋转是 $-I$，距离 $\sqrt5$；不能把普通极分解答案直接交给受限旋转问题。

## 七、提交与接受标准

1. 给出 $F_0,S_0,L_0$ 的完整矩阵及 Hermite 推导，证明缺陷极限没有丢掉导数信息。
2. 对四个输入，指数、平方根、对数与八十位参照的矩阵范数误差满足脚本的相对尺度容差；对接近零的单个元素不使用无意义的相对误差。
3. 分别报告 $S^2-T$、$e^L-T$ 和主分支条件。对数截断界必须注明只控制精确算术逼近，不冒充全部浮点误差界。
4. 导数同时通过精确极限公式、块指数及有限差分的合理步长区间；解释为何 $e^TE$ 不是一般答案。
5. Sylvester 测试通过原坐标残差，正确处理实二阶块、减号约定、scale 和 info；敏感性用 sep 而非仅用谱距离。
6. 四节点扩散的缺陷积分界覆盖真实误差，且说明半群界依赖 $t,\tau\ge0$ 及对称负定条件。
7. 逐项列出实际操作计数。十三阶 Padé 的七次乘法与一次求解、块 Taylor 的二十三次小块乘法和一次耦合求解、对数的四次根与十次求解，属于不同规模的运算，不能把次数直接相加比较速度。
8. 报告 nilpotent 大范数例的过度缩放现象；其代数精确性是跨平台结论，实测舍入误差只是所记录实现的结果。
9. Denman–Beavers 和极因子迭代都通过自身的不变量与残差检查；额外的主分支或行列式约束须单独验证。

达到这些标准，说明读者能够同时回答三个问题：这个矩阵函数被如何定义，算法实际计算了什么，以及哪些证据足以支持结果的准确性。