线性注意力的线性代数前置知识

1. Householder 矩阵

1.1 定义

设 $\bm{u} \in \mathbb{R}^n$ 且 $\bm{u}^{\intercal}\bm{u} = 1$（即 $\|\bm{u}\|_2 = 1$），定义

$$\mathbf{H} = \mathbf{I} - 2\bm{u}\bm{u}^{\intercal} \in \mathbb{R}^{n \times n}$$

题外一句：Householder 矩阵也是做 QR 分解的重要工具——逐列用反射把对角线以下的元素打成零，就得到 $\mathbf{A} = \mathbf{Q}\mathbf{R}$，比 Gram-Schmidt 数值上稳得多。

1.1.1 矩阵形式是怎么被「提」出来的

定义式 $\mathbf{H} = \mathbf{I} - 2\bm{u}\bm{u}^{\intercal}$ 常被当成天降公式，其实它只是把「反射两步走」这句话提取公因子的结果。取 $n = 2$、$\bm{u} = \bm{e}_1 = [1, 0]^{\intercal}$、$\bm{x} = [3, 2]^{\intercal}$ 走一遍：

几何上：$\bm{u} = \bm{e}_1$ 时镜面就是 $x_2$ 轴。从 $\bm{x}$ 出发沿 $-\bm{u}$ 方向走一步 $-(\bm{u}^{\intercal}\bm{x})\bm{u}$ 落到镜面上（此时第一坐标归零），再走同样的一步就到 $\bm{x}'$。两步等长，这正是系数 2 的来源。

推导上唯一需要的一步是 $\bm{u}^{\intercal}\bm{x}$ 是标量，标量与向量相乘可换序，所以 $(\bm{u}^{\intercal}\bm{x})\bm{u} = \bm{u}(\bm{u}^{\intercal}\bm{x})$。换序之后 $\bm{x}$ 变成右端公因子，反着用分配律提出来：

$$\bm{x}' = \bm{x} - 2(\bm{u}^{\intercal}\bm{x})\bm{u} = \bm{x} - 2\bm{u}(\bm{u}^{\intercal}\bm{x}) = \big(\mathbf{I} - 2\bm{u}\bm{u}^{\intercal}\big)\bm{x} = \mathbf{H}\bm{x}$$

换序是全部的技巧所在：$(\bm{u}^{\intercal}\bm{x})\bm{u}$ 里 $\bm{x}$ 被夹在中间提不出来，写成 $\bm{u}(\bm{u}^{\intercal}\bm{x})$ 后 $\bm{u}\bm{u}^{\intercal}$ 自然聚成一个矩阵。

1.2 对称性证明

$$\mathbf{H}^{\intercal} = (\mathbf{I} - 2\bm{u}\bm{u}^{\intercal})^{\intercal} = \mathbf{I}^{\intercal} - 2(\bm{u}\bm{u}^{\intercal})^{\intercal} = \mathbf{I} - 2\bm{u}\bm{u}^{\intercal} = \mathbf{H} \qquad \blacksquare$$

关键一步用的是乘积转置法则 $(\mathbf{A}\mathbf{B})^{\intercal} = \mathbf{B}^{\intercal}\mathbf{A}^{\intercal}$（转置后因子次序颠倒），以及 $(\mathbf{A}^{\intercal})^{\intercal} = \mathbf{A}$：

$$(\bm{u}\bm{u}^{\intercal})^{\intercal} = (\bm{u}^{\intercal})^{\intercal}\bm{u}^{\intercal} = \bm{u}\bm{u}^{\intercal}$$

外积矩阵因此天然对称。

1.3 广义形式：把归一化收进一个系数 $\tau$

$\mathbf{H} = \mathbf{I} - 2\bm{u}\bm{u}^{\intercal}$ 要求 $\|\bm{u}\|_2 = 1$。想对任意非零 $\bm{u}$ 都能写，把归一化带在系数里就行：

$$\mathbf{H} = \mathbf{I} - \frac{2\bm{u}\bm{u}^{\intercal}}{\bm{u}^{\intercal}\bm{u}} = \mathbf{I} + \tau\,\bm{u}\bm{u}^{\intercal}, \qquad \tau := -\frac{2}{\bm{u}^{\intercal}\bm{u}}$$

把 $\tau$ 当成一个普通标量系数写在外积前面，而不是留在分母上，后面的推导会干净很多。$\|\bm{u}\|_2 = 1$ 时 $\tau = -2$，退回 §1.1 的定义。

系数取 $-2/(\bm{u}^{\intercal}\bm{u})$ 不是美学选择，而是被「$\mathbf{H}$ 是反射」逆推出来的。反射必须满足对合性 $\mathbf{H}^2 = \mathbf{I}$（照两次镜子回到原处），展开：

$$(\mathbf{I} + \tau\bm{u}\bm{u}^{\intercal})^2 = \mathbf{I} + 2\tau\bm{u}\bm{u}^{\intercal} + \tau^2\bm{u}\underbrace{(\bm{u}^{\intercal}\bm{u})}_{\text{标量}}\bm{u}^{\intercal} = \mathbf{I} + \tau\big(2 + \tau\,\bm{u}^{\intercal}\bm{u}\big)\bm{u}\bm{u}^{\intercal}$$

要它等于 $\mathbf{I}$，就要 $\tau(2 + \tau\bm{u}^{\intercal}\bm{u}) = 0$，除去平常解 $\tau = 0$ 只剩

$$\tau = -\frac{2}{\bm{u}^{\intercal}\bm{u}}$$

所以 $\tau$ 取任何值时 $\mathbf{H}$ 都对称（外积天然对称，与系数无关），但只有这一个 $\tau$ 才让它成为对合的正交反射、才有「沿法向走两步等长」的几何意义。其余 $\tau$ 给出的是一般的秩一修正 $\mathbf{I} + \tau\bm{u}\bm{u}^{\intercal}$，不再是等距变换；但 §3 的推导对 $\tau$ 不做任何限制，对这类一般秩一修正同样成立。

2. WY 表达式：把一串反射压成一个低秩修正

结论先给：$p$ 个 Householder 反射的乘积，永远能写成单位矩阵加一个秩不超过 $p$ 的修正

$$\mathbf{H}_1\mathbf{H}_2\cdots\mathbf{H}_p = \mathbf{I} + \mathbf{W}\mathbf{Y}^{\intercal}, \qquad \mathbf{W},\,\mathbf{Y} \in \mathbb{R}^{n \times p}$$

证明用归纳构造，一列一列把反射「吸收」进 $\mathbf{W},\mathbf{Y}$。

2.1 基例 $p = 1$

$\mathbf{H}_1 = \mathbf{I} - 2\bm{u}_1\bm{u}_1^{\intercal} = \mathbf{I} + \mathbf{W}_1\mathbf{Y}_1^{\intercal}$，取

$$\mathbf{W}_1 = \big[-2\bm{u}_1\big], \qquad \mathbf{Y}_1 = \big[\bm{u}_1\big]$$

均为 $n \times 1$。命题对 $p = 1$ 成立。

2.2 归纳步：一列一列「吸收」反射

归纳假设：已存在 $\mathbf{W}_k,\,\mathbf{Y}_k \in \mathbb{R}^{n \times k}$ 使 $\mathbf{H}_1\cdots\mathbf{H}_k = \mathbf{I} + \mathbf{W}_k\mathbf{Y}_k^{\intercal}$。要证 $\mathbf{H}_1\cdots\mathbf{H}_k\mathbf{H}_{k+1}$ 也能写成 $\mathbf{I} + \mathbf{W}_{k+1}\mathbf{Y}_{k+1}^{\intercal}$，其中 $\mathbf{W}_{k+1},\mathbf{Y}_{k+1} \in \mathbb{R}^{n \times (k+1)}$。

记新反射方向 $\bm{u} = \bm{u}_{k+1}$，并记

$$\bm{w} := -2\,(\mathbf{I} + \mathbf{W}_k\mathbf{Y}_k^{\intercal})\,\bm{u}$$

则把 $(\mathbf{I} + \mathbf{W}_k\mathbf{Y}_k^{\intercal})(\mathbf{I} - 2\bm{u}\bm{u}^{\intercal})$ 乘开后的四项可以合并（用 $(\mathbf{I} + \mathbf{W}_k\mathbf{Y}_k^{\intercal})\bm{u}\bm{u}^{\intercal} = \big[(\mathbf{I} + \mathbf{W}_k\mathbf{Y}_k^{\intercal})\bm{u}\big]\bm{u}^{\intercal}$，$\bm{u}^{\intercal}$ 原样留在最右）

$$\mathbf{I} + \mathbf{W}_k\mathbf{Y}_k^{\intercal} - 2\bm{u}\bm{u}^{\intercal} - 2\mathbf{W}_k\mathbf{Y}_k^{\intercal}\bm{u}\bm{u}^{\intercal} = \mathbf{I} + \mathbf{W}_k\mathbf{Y}_k^{\intercal} + \bm{w}\bm{u}^{\intercal}$$

于是

$$\mathbf{H}_1\cdots\mathbf{H}_{k+1} = \mathbf{I} + \big[\,\mathbf{W}_k \;\; \bm{w}\,\big]\big[\,\mathbf{Y}_k \;\; \bm{u}\,\big]^{\intercal}$$

即 $\mathbf{W}_{k+1} = [\mathbf{W}_k,\ \bm{w}]$，$\mathbf{Y}_{k+1} = [\mathbf{Y}_k,\ \bm{u}_{k+1}]$。定理得证。$\blacksquare$

3. compact WY：把纠缠收进一个上三角因子 $\mathbf{T}$

WY 有个不便：$\mathbf{W}$ 的列互相纠缠。由 §2.2 的递推回代可得 $\bm{w}_t = -2\mathbf{H}_1\cdots\mathbf{H}_{t-1}\bm{u}_t$，每根列都带着一个前缀乘积，$\mathbf{W}$ 里存的已经不是原始法向了。

1989 年 Schreiber–Van Loan 的 compact WY 把全部纠缠塞进一个 $p \times p$ 小矩阵，让两个 $n \times p$ 因子里只剩下原始法向：

$$\mathbf{H}_1\cdots\mathbf{H}_p = \mathbf{I} + \mathbf{Y}\mathbf{T}\mathbf{Y}^{\intercal}, \qquad \mathbf{Y} = [\,\bm{u}_1, \dots, \bm{u}_p\,], \quad \mathbf{T} \in \mathbb{R}^{p \times p}\ \text{上三角}$$

下面用 §1.3 的广义形式 $\mathbf{H} = \mathbf{I} + \tau\bm{u}\bm{u}^{\intercal}$ 写，$\bm{u}_t$ 不需归一，$\tau_t$ 也不限于 $-2/(\bm{u}_t^{\intercal}\bm{u}_t)$，当成一串给定系数即可。

递推构造：

$$\mathbf{T}_1 = [\tau_1], \qquad \mathbf{T}_{k+1} = \begin{bmatrix} \mathbf{T}_k & \tau_{k+1}\,\mathbf{T}_k\,(\mathbf{Y}_k^{\intercal}\bm{u}_{k+1}) \\ \bm{0} & \tau_{k+1} \end{bmatrix}$$

证明与 §2.2 是同一个机制，只是新增的一列不往 $n$ 维因子里塞，而往 $\mathbf{T}$ 里塞。记 $\bm{u} = \bm{u}_{k+1}$、$\tau = \tau_{k+1}$，乘开：

$$(\mathbf{I} + \mathbf{Y}_k\mathbf{T}_k\mathbf{Y}_k^{\intercal})(\mathbf{I} + \tau\bm{u}\bm{u}^{\intercal}) = \mathbf{I} + \mathbf{Y}_k\mathbf{T}_k\mathbf{Y}_k^{\intercal} + \tau\bm{u}\bm{u}^{\intercal} + \tau\,\mathbf{Y}_k\mathbf{T}_k(\mathbf{Y}_k^{\intercal}\bm{u})\,\bm{u}^{\intercal}$$

三个修正项正好是 $\big[\mathbf{Y}_k,\, \bm{u}\big]$ 夹一个 $(k+1) \times (k+1)$ 矩阵再夹 $\big[\mathbf{Y}_k,\, \bm{u}\big]^{\intercal}$ 的展开结果：$\mathbf{Y}_k\mathbf{T}_k\mathbf{Y}_k^{\intercal}$ 对应左上块，$\tau\bm{u}\bm{u}^{\intercal}$ 对应右下角，交叉项 $\tau\mathbf{Y}_k\mathbf{T}_k(\mathbf{Y}_k^{\intercal}\bm{u})\bm{u}^{\intercal}$ 对应右上列，左下为零：

$$= \mathbf{I} + \big[\mathbf{Y}_k,\, \bm{u}\big] \begin{bmatrix} \mathbf{T}_k & \tau\mathbf{T}_k(\mathbf{Y}_k^{\intercal}\bm{u}) \\ \bm{0} & \tau \end{bmatrix} \big[\mathbf{Y}_k,\, \bm{u}\big]^{\intercal}$$

即 $\mathbf{Y}_{k+1} = [\mathbf{Y}_k,\ \bm{u}_{k+1}]$、$\mathbf{T}_{k+1}$ 如上，定理得证。$\blacksquare$

三个副产品直接从递推读出来：

• $\mathbf{T}$ 上三角。左下块永远是 $\bm{0}$，上三角性逐层遗传。
• 对角线就是系数。$\mathrm{diag}(\mathbf{T}) = (\tau_1, \dots, \tau_p)$，单位向量情形下全是 $-2$。
• 与 §2 的关系是 $\mathbf{W} = \mathbf{Y}\mathbf{T}$。两个表达式并不独立：compact WY 把 WY 里那些前缀乘积分解成「原始法向 $\times$ 上三角系数」。纠缠没有消失，只是从 $n \times p$ 的 $\mathbf{W}$ 里搬到了 $p \times p$ 的 $\mathbf{T}$ 里，而 $p \ll n$。

4. 另一条路：用正交条件直接解出 $\mathbf{T}$

§3 的递推是串行的：要算 $\mathbf{T}_{k+1}$ 必须先有 $\mathbf{T}_k$。Joffrain 等人（2006）指出 $\mathbf{T}$ 有闭式，不靠归纳、不靠前缀，一步得出。

为与文献记号一致，本节写成减号形式：

$$\mathbf{Q} = \mathbf{H}_1\cdots\mathbf{H}_p = \mathbf{I} - \mathbf{U}\mathbf{S}\mathbf{U}^{\intercal}, \qquad \mathbf{U} = [\,\bm{u}_1, \dots, \bm{u}_p\,]$$

对照 §3：把那里的上三角因子记作 $\mathbf{T}_{\text{SVL}}$，则 $\mathbf{U} = \mathbf{Y}$、$\mathbf{S} = -\mathbf{T}_{\text{SVL}}$，同为上三角；由 $\mathrm{diag}(\mathbf{S}) = (2/\bm{u}_1^{\intercal}\bm{u}_1, \dots, 2/\bm{u}_p^{\intercal}\bm{u}_p)$ 全非零，且三角阵行列式等于对角元之积，知 $\mathbf{S}$ 可逆。本节的 $\mathbf{T} := \mathbf{S}^{-1}$ 与 $\mathbf{T}_{\text{SVL}}$ 差一个求逆和符号，不是同一个矩阵。以下设 $\mathbf{U}$ 列满秩。

4.1 可用的条件是正交，不是对合

先排除一个很容易误用的条件。每个 $\mathbf{H}_i$ 都对称（§1.2），但对称矩阵的乘积一般不对称：

$$(\mathbf{H}_1\mathbf{H}_2)^{\intercal} = \mathbf{H}_2^{\intercal}\mathbf{H}_1^{\intercal} = \mathbf{H}_2\mathbf{H}_1 \ne \mathbf{H}_1\mathbf{H}_2$$

转置把因子次序翻了过来，而反射一般不交换。所以 $\mathbf{Q}^{\intercal} = \mathbf{Q}$ 不成立，$\mathbf{Q}^2 = \mathbf{I}$ 也不成立——两面镜子依次反射是一个旋转，转两次不回原处。真正能用的是正交性：正交矩阵的乘积仍正交，所以

$$\mathbf{Q}^{\intercal}\mathbf{Q} = \mathbf{I}$$

好在下面的代数只需要这一条，结论一字不改。

4.2 展开正交条件

把 $\mathbf{Q} = \mathbf{I} - \mathbf{U}\mathbf{S}\mathbf{U}^{\intercal}$ 代入 $\mathbf{Q}^{\intercal}\mathbf{Q} = \mathbf{I}$：

$$(\mathbf{I} - \mathbf{U}\mathbf{S}^{\intercal}\mathbf{U}^{\intercal})(\mathbf{I} - \mathbf{U}\mathbf{S}\mathbf{U}^{\intercal}) = \mathbf{I} - \mathbf{U}\mathbf{S}^{\intercal}\mathbf{U}^{\intercal} - \mathbf{U}\mathbf{S}\mathbf{U}^{\intercal} + \mathbf{U}\mathbf{S}^{\intercal}(\mathbf{U}^{\intercal}\mathbf{U})\mathbf{S}\mathbf{U}^{\intercal} = \mathbf{I}$$

消去 $\mathbf{I}$，左提 $\mathbf{U}$、右提 $\mathbf{U}^{\intercal}$：

$$\mathbf{U}\big(-\mathbf{S}^{\intercal} - \mathbf{S} + \mathbf{S}^{\intercal}(\mathbf{U}^{\intercal}\mathbf{U})\mathbf{S}\big)\mathbf{U}^{\intercal} = \bm{0}$$

$\mathbf{U}$ 列满秩（$\mathbf{U}\bm{x} = \bm{0} \Rightarrow \bm{x} = \bm{0}$），因此括号内的 $p \times p$ 矩阵必须为零：

$$\mathbf{S}^{\intercal} + \mathbf{S} = \mathbf{S}^{\intercal}(\mathbf{U}^{\intercal}\mathbf{U})\mathbf{S}$$

$\mathbf{S}$ 可逆，左乘 $\mathbf{S}^{-\intercal}$、右乘 $\mathbf{S}^{-1}$，右端夹着的 $\mathbf{S}^{\intercal},\mathbf{S}$ 被抵消：

$$\mathbf{S}^{-\intercal}(\mathbf{S}^{\intercal} + \mathbf{S})\mathbf{S}^{-1} = \mathbf{U}^{\intercal}\mathbf{U} \qquad\Longrightarrow\qquad \mathbf{S}^{-1} + \mathbf{S}^{-\intercal} = \mathbf{U}^{\intercal}\mathbf{U}$$

令 $\mathbf{T} = \mathbf{S}^{-1}$（上三角阵的逆仍上三角），得到全部机关：

$$\mathbf{T} + \mathbf{T}^{\intercal} = \mathbf{U}^{\intercal}\mathbf{U}$$

4.3 一个方程定住整个 $\mathbf{T}$

$\mathbf{T}$ 上三角、$\mathbf{T}^{\intercal}$ 下三角，两者相加时不同位置互不干扰，所以逐位置比对就能定住每一项（记 $\mathbf{G} = \mathbf{U}^{\intercal}\mathbf{U}$，它对称）：

• 严格上三角（$i < j$）：$\mathbf{T}^{\intercal}$ 在此为零，所以 $\mathbf{T}_{ij} = \mathbf{G}_{ij} = \bm{u}_i^{\intercal}\bm{u}_j$。
• 对角线：两项相等，$2\mathbf{T}_{ii} = \mathbf{G}_{ii}$，所以 $\mathbf{T}_{ii} = \tfrac12\bm{u}_i^{\intercal}\bm{u}_i$。
• 严格下三角：$\mathbf{T}_{ij} = 0$（上三角的定义）。

写成一行：

$$\mathbf{T} = \mathrm{striu}(\mathbf{U}^{\intercal}\mathbf{U}) + \tfrac12\,\mathrm{diag}(\mathbf{U}^{\intercal}\mathbf{U}), \qquad \mathbf{Q} = \mathbf{I} - \mathbf{U}\mathbf{T}^{-1}\mathbf{U}^{\intercal}$$

其中 $\mathrm{striu}$ 取严格上三角部分。换句话说，$\mathbf{T}$ 就是 Gram 矩阵 $\mathbf{U}^{\intercal}\mathbf{U}$ 的「上三角一半」。单位向量时对角线全为 $\tfrac12$，$\mathbf{S} = \mathbf{T}^{-1}$ 对角线全为 $2$，正对应 §1.1 的系数。$\mathbf{T}$ 由方程唯一确定，因此使 $\mathbf{Q} = \mathbf{I} - \mathbf{U}\mathbf{S}\mathbf{U}^{\intercal}$ 成立的上三角可逆 $\mathbf{S}$ 也唯一。

4.4 为什么这个闭式适合并行

两条路得到的是同一个 $\mathbf{T}$，但计算结构完全不同：

关键差异在于 $\mathbf{U}^{\intercal}\mathbf{U}$ 是一个普通的 Gram 矩阵乘法：$p^2$ 个内积互不依赖，能一次 GEMM 算完，而递推的第 $k+1$ 列必须等第 $k$ 列。代价是多了一个 $p \times p$ 上三角求逆（或三角求解），但 $p \ll n$，这部分开销可忽略。

更有用的是分块可合并。把 $\mathbf{U} = [\mathbf{U}_1, \mathbf{U}_2]$ 分两块，则

$$\mathbf{U}^{\intercal}\mathbf{U} = \begin{bmatrix} \mathbf{U}_1^{\intercal}\mathbf{U}_1 & \mathbf{U}_1^{\intercal}\mathbf{U}_2 \\ \mathbf{U}_2^{\intercal}\mathbf{U}_1 & \mathbf{U}_2^{\intercal}\mathbf{U}_2 \end{bmatrix} \qquad\Longrightarrow\qquad \mathbf{T} = \begin{bmatrix} \mathbf{T}_1 & \mathbf{U}_1^{\intercal}\mathbf{U}_2 \\ \bm{0} & \mathbf{T}_2 \end{bmatrix}$$

两个子块的 $\mathbf{T}_1, \mathbf{T}_2$ 可以并行算，合并时只需补上跨块耦合块 $\mathbf{U}_1^{\intercal}\mathbf{U}_2$——这正是 LAPACK 分块 QR 里 dlarft 合并相邻 panel 的做法，也是把反射串映到 Tensor Core 上时想要的结构。

最后注意一个方向约定：上三角对应的是 $\mathbf{H}_1\mathbf{H}_2\cdots\mathbf{H}_p$ 这个顺序。反序乘积 $\mathbf{H}_p\cdots\mathbf{H}_1$ 换成 $\mathbf{T}^{-\intercal}$，即下三角那一半；弄反会得到错的矩阵，实现时这是最容易踩的坑。

TileLang 实战：KDA 从零到一–Kimi Delta Attention

前三篇分别实现了无遗忘的 chunked 线性注意力、标量衰减、以及带删除项的 Gated DeltaNet。本篇是系列收尾：把 GDN 的标量门 $\alpha_t$ 换成逐通道向量门 $\bm{a}_t \in \mathbb{R}^{d_k}$。

$$\mathbf{S}_t = \mathbf{S}_{t-1}\operatorname{Diag}(\bm{a}_t)\big(\mathbf{I} - \beta_t\bm{k}_t\bm{k}_t^{\intercal}\big) + \beta_t\bm{v}_t\bm{k}_t^{\intercal}$$

改动只有一处：$\alpha_t$ 变成了 $\operatorname{Diag}(\bm{a}_t)$。但这一处把前三篇积累的所有便利拆掉了–衰减不再能从矩阵里外提，$\Gamma$ 从 $C\times C$ 变成 $C \times C \times d_k$，而累积衰减积的通道离散度会直接把 fp16 打穿。

本篇聚焦实现。KDA 的数学推导（递推式、逐通道 WY 表示、UT 变换、下界衰减与满秩门控的动机）见《KDA 的来龙去脉》§3–§4，这里不重复；本文只做一件事：把那些公式落成能跑的 kernel，并量化每一处数值边界。

前三篇见 Chunked 线性注意力、标量衰减、Gated DeltaNet。

0. 符号与三处变化

沿用前三篇：$\beta_t$ 写入强度、$C$ 块长、$[t]$ 块序号、$r \in [1,C]$ 块内位置、$\mathbf{S}^\intercal \in \mathbb{R}^{d_k \times d_v}$ 为 kernel 存储布局。

本篇的核心量改为累积 log 衰减（沿用《来龙去脉》的记号）：

$$\bm{\gamma}_i = \sum_{s=1}^{i}\bm{g}_s \in \mathbb{R}^{d_k}, \qquad \bm{g}_s = \log\bm{a}_s < \bm{0}\ \text{（逐分量）}$$

注意 $\bm{\gamma}_i$ 是向量–这是与前两篇最本质的区别。前三篇的 $\gamma^r$ 是标量，$\Gamma_{ij} = \gamma_i/\gamma_j$ 是一张 $C \times C$ 的表；本篇 $\bm{\gamma}_i - \bm{\gamma}_j$ 是 $d_k$ 维向量，"衰减矩阵"概念上是 $C \times C \times d_k$，不可能物化。

三处结构性变化：

第三行是全篇的技术核心。

1. 衰减长在内部：问题与出路

1.1 朴素做法的代价

$M_{ci} = (\bm{k}_c \odot e^{\bm{\gamma}_c - \bm{\gamma}_i})^\intercal\bm{k}_i$ 里的指数依赖 $(c, i, d)$ 三个下标。直接算就是 $C^2/2$ 次长度 $d_k$ 的加权内积，每次都要现算 $d_k$ 个 exp。

这不是常数因子问题–它把一次 GEMM 变成了 $C^2/2$ 个独立的向量运算，完全用不上 Tensor Core。$C = 64$、$d_k = 128$ 时是 2016 次加权内积、约 26 万次 exp2。

1.2 出路：指数可分离

关键观察是指数可以按通道拆开：

$$e^{\bm{\gamma}_c - \bm{\gamma}_i} = e^{\bm{\gamma}_c} \odot e^{-\bm{\gamma}_i}$$

于是：

$$M_{ci} = \sum_{d}\underbrace{\big(k_c[d]\,e^{\gamma_c[d]}\big)}_{\widetilde{K}^{+}[c,d]}\underbrace{\big(k_i[d]\,e^{-\gamma_i[d]}\big)}_{\widetilde{K}^{-}[i,d]} = \big(\widetilde{\mathbf{K}}^{+}\widetilde{\mathbf{K}}^{-\intercal}\big)_{ci}$$

一次 GEMM 解决，前置两个 $C \times d_k$ 的逐元素加权。实测拆分与朴素计算的差异 $1.94 \times 10^{-16}$，代数上完全等价。

1.3 代价：$e^{-\bm{\gamma}_i}$ 是大于 1 的量

第二篇讲过一个教训：把比值 $\gamma_i/\gamma_j$ 因式分解成 $\gamma_i \cdot (1/\gamma_j)$ 会物化一个指数增长的量，fp16 下 $\alpha \le 0.8$ 即 NaN。这里是同一个陷阱的逐通道版本–$e^{-\gamma_i[d]} \ge 1$，且随 $i$ 与通道衰减强度指数增长。

但本篇的处境和第二篇不同：那里有替代方案（保留比值形式），这里没有。逐通道衰减无法外提，不拆就用不上 Tensor Core。所以问题从"要不要拆"变成了"怎样让拆分在数值上安全"。

实测 $e^{-\bm{\gamma}}$ 的最大值（$C = 64$、$d_k = 128$）：

门控下界在 0.9 以上时，$e^{-\bm{\gamma}}$ 连 fp16 都装得下。 这就把 KDA 论文里"下界衰减（lower-bounded decay）"这个设计从模型层面的技巧，变成了 kernel 能否用 Tensor Core 的前提条件。

2. 下界衰减：不是精度调优，是可行性前提

2.1 无下界时 fp16 全线归零

$e^{-\bm{\gamma}}$ 会溢出，另一头 $e^{\bm{\gamma}}$ 会下溢。而逐通道门让后者严重得多–总有一些通道学到很小的 $a$，它们的 $\gamma$ 累积得最快。

实测 $\min e^{\bm{\gamma}_C}$ 与 fp16 下归零的通道数（$d_k = 128$）：

$C = 128$、下界 0.5 时全部 128 个通道在 fp16 下归零–状态被彻底清空，kernel 输出恒为块内项，跨块信息完全丢失。

加下界后：

下界 0.9 是分界线。 这一个数字同时解决了两头：$e^{\bm{\gamma}}$ 不下溢、$e^{-\bm{\gamma}}$ 不溢出。

2.2 通道离散度：逐通道门特有的问题

标量门下 $\gamma^C$ 是一个数；逐通道门下它是 $d_k$ 个数，而它们的极差决定了同一个 chunk 内不同通道的数值尺度差异：

极差 $47$ 意味着同一个 fragment 里最强和最弱通道的数值相差 20 个数量级–任何浮点格式都无法同时表示。这是逐通道门独有的病：标量门只需要担心整体的尺度漂移，逐通道门还要担心通道之间的尺度撕裂。

下界 0.9 把极差压到 1.74（$e^{1.74} = 5.7$），完全可控。

2.3 sub-chunk：第二道保险

即使有下界，$C = 128$ 时 $\max e^{-\bm{\gamma}} = 1.79\times10^3$（下界 0.9）–fp16 装得下但余量不多。把块内再切成 sub-chunk，每个 sub-chunk 内部重新起算 $\bm{\gamma}$：

指数跨度只取决于 sub-chunk 长度，与总块长无关。这与第二篇"累积积按 chunk 重置"是同一个道理，只是又下降了一层：chunk 重置控制 $\bm{\gamma}$ 本身，sub-chunk 重置控制 $e^{-\bm{\gamma}}$ 的动态范围。

代价是 sub-chunk 之间需要额外的状态传递，块内变成两层循环。$C = 64$ + 下界 0.9 时 $\max e^{-\bm{\gamma}} = 55$，不需要 sub-chunk；$C = 128$ 时建议切 32 或 64。

3. 四层参考实现

沿用前三篇框架。

3.1 参考 A：逐 token 递归

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

def ref_A(Q, K, V, g, beta):
    """g: (B,N,H,DK) 逐通道 log 门控，g < 0"""
    B, N, H, DK = Q.shape
    DV = V.shape[-1]
    O = np.zeros((B, N, H, DV))
    for b in range(B):
        for h in range(H):
            S = np.zeros((DV, DK))                      # 论文方向
            for t in range(N):
                k, v, q = K[b, t, h], V[b, t, h], Q[b, t, h]
                a, bt = np.exp(g[b, t, h]), beta[b, t, h]
                # Diag(a) 在 S 与 Householder 之间
                S = S @ np.diag(a) @ (np.eye(DK) - bt * np.outer(k, k)) \
                    + bt * np.outer(v, k)
                O[b, t, h] = S @ q
    return O

$\operatorname{Diag}(\bm{a}_t)$ 的位置很关键：它夹在 $\mathbf{S}_{t-1}$ 与 Householder 之间。写成 $\operatorname{Diag}(\bm{a})\mathbf{S}$ 或挪到括号外都会算错。

3.2 参考 B：chunkwise + 逐通道 UT 变换

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

def ref_B(Q, K, V, g, beta, C):
    B, N, H, DK = Q.shape
    DV = V.shape[-1]
    O = np.zeros((B, N, H, DV))
    for b in range(B):
        for h in range(H):
            S = np.zeros((DK, DV))                      # 转置布局
            for c in range(N // C):
                sl = slice(c * C, (c + 1) * C)
                Qc, Kc, Vc = Q[b, sl, h], K[b, sl, h], V[b, sl, h]
                gc, bt = g[b, sl, h], beta[b, sl, h]

                gam = np.cumsum(gc, axis=0)             # γ_i ∈ R^{DK}
                gC  = gam[-1]

                # M[c,i] = (k_c ⊙ e^{γ_c-γ_i})·k_i,  i < c
                M = np.zeros((C, C))
                for cc in range(C):
                    for i in range(cc):
                        M[cc, i] = np.dot(Kc[cc] * np.exp(gam[cc] - gam[i]), Kc[i])

                L  = np.tril(bt[None, :] * M, -1)       # L[c,i] = β_i M[c,i]
                Tm = np.linalg.inv(np.eye(C) + L)
                Ah = bt[:, None] * Tm                   # diag(β) T

                W  = Ah @ (Kc * np.exp(gam))            # C × DK
                U  = Ah @ Vc                            # C × DV
                Vt = U - W @ S                          # 伪值 Ṽ

                # A^qk[c,j] = (q_c ⊙ e^{γ_c-γ_j})·k_j,  j ≤ c
                Aqk = np.zeros((C, C))
                for cc in range(C):
                    for j in range(cc + 1):
                        Aqk[cc, j] = np.dot(Qc[cc] * np.exp(gam[cc] - gam[j]), Kc[j])

                O[b, sl, h] = (Qc * np.exp(gam)) @ S + Aqk @ Vt
                S = np.diag(np.exp(gC)) @ S + (Kc * np.exp(gC - gam)).T @ Vt
    return O

这份参考刻意用朴素双循环算 $\mathbf{M}$ 与 $\mathbf{A}^{qk}$–慢，但与公式逐字对应，作为基准可信。§3.3 的 kernel 镜像才用 GEMM 化写法。

$L_{ci} = \beta_i M_{ci}$ 与 $\hat{\mathbf{A}} = \operatorname{diag}(\beta)\mathbf{T}$ 这一对下标必须配套。实测另一种等价写法是 $L_{ci} = \beta_c M_{ci}$ 配 $\hat{\mathbf{A}} = \mathbf{T}\operatorname{diag}(\beta)$，两者都对，混用则错。

3.3 参考 C：kernel 控制流镜像（GEMM 化）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

def ref_C(Q, K, V, g, beta, C, block_DV):
    B, N, H, DK = Q.shape
    DV = V.shape[-1]
    O = np.zeros((B, N, H, DV))
    for bv in range(DV // block_DV):                    # ← grid.x
        dv = slice(bv * block_DV, (bv + 1) * block_DV)
        for bbh in range(B * H):                        # ← grid.y
            b, h = bbh // H, bbh % H
            S = np.zeros((DK, block_DV))
            for c in range(N // C):                     # ← T.Pipelined
                sl = slice(c * C, (c + 1) * C)
                Qc, Kc = Q[b, sl, h], K[b, sl, h]
                Vc = V[b, sl, h, dv]
                gc, bt = g[b, sl, h], beta[b, sl, h]

                gam = np.cumsum(gc, axis=0)
                gC  = gam[-1]
                Ep, Em = np.exp(gam), np.exp(-gam)      # e^{γ}, e^{-γ}

                # ── 指数可分离：两次 GEMM 代替 C²/2 次加权内积 ──
                Kp, Km = Kc * Ep, Kc * Em
                M   = np.tril(Kp @ Km.T, -1)            # 严格下三角
                Aqk = np.tril((Qc * Ep) @ Km.T)         # 含对角

                L  = np.tril(bt[None, :] * M, -1)
                # 前向替换解 T @ [Kp | Vc]，再逐行乘 β（因为 Â = diag(β)T）
                RHS = np.concatenate([Kc * Ep, Vc], axis=1)
                X = np.zeros_like(RHS)
                for r in range(C):
                    X[r] = RHS[r]
                    for j in range(r):
                        X[r] -= L[r, j] * X[j]
                X = bt[:, None] * X                     # ← β 在解完之后乘
                W, U = X[:, :DK], X[:, DK:]

                Vt = U - W @ S
                O[b, sl, h, dv] = (Qc * Ep) @ S + Aqk @ Vt
                S = Ep[-1][:, None] * S + (Kc * np.exp(gC - gam)).T @ Vt
    return O

三处与参考 B 不同的实现选择：

1. $\mathbf{M}$ 与 $\mathbf{A}^{qk}$ 走 GEMM：$\widetilde{\mathbf{K}}^{+}\widetilde{\mathbf{K}}^{-\intercal}$。$\mathbf{A}^{qk}$ 含对角（$j \le c$），$\mathbf{M}$ 不含（$i < c$）。
2. $\mathbf{W}$ 与 $\mathbf{U}$ 拼成一次前向替换：两者共用同一个 $(\mathbf{I}+\mathbf{L})$，拼接后只解一遍，省一半串行开销。
3. $\beta$ 必须在解完三角系统之后乘。$\hat{\mathbf{A}} = \operatorname{diag}(\beta)\mathbf{T}$ 展开是「先 $\mathbf{T}$ 作用、再逐行乘 $\beta$」；若把 $\beta$ 提前乘进右端项，算的就是 $\mathbf{T}\operatorname{diag}(\beta)$–那是另一种配对（需搭配 $L_{ci} = \beta_c M_{ci}$），混用则错。这个错误在 $\beta$ 全部相等时完全看不出来，我第一次写就踩了：$\beta$ 随机时相对 L2 达 $1.3\times10^{-1}$。
4. 状态更新用逐行乘代替 $\operatorname{Diag}$：Ep[-1][:, None] * S 就是 $\operatorname{Diag}(e^{\bm{\gamma}_C})\mathbf{S}$，不物化对角矩阵。

3.4 一致性验证

$B=2, H=2, N=12, d_k=d_v=4, C=4$，$\bm{a} \sim \mathcal{U}(0.90, 0.999)^{d_k}$、$\beta \sim \mathcal{U}(0.1, 0.9)$，$\bm{q},\bm{k}$ 已 L2 归一化，fp64：

三个退化检验（本篇比 GDN 多一个）：

三个 $\text{block}_{DV}$（1 / 2 / 4）的相对 L2 分别为 $2.48$ / $2.57$ / $2.57 \times 10^{-16}$–切 DV 零依赖，与第一篇的结论一致。

第一个是本篇独有且最重要的：把逐通道门退化成标量门，必须精确回到第三篇的结果。这一步能抓住"指数分离时把通道维和位置维搞混"这类错误–而那类错误在通道值本来就相同时会隐身。

4. TileLang kernel

4.0 寄存器账

$d_k = d_v = 128$、$C = 64$、$\text{block}_{DV} = 32$、128 线程：

320 超限了。 逐通道门带来的 $[C, d_k]$ 量级 fragment（$e^{\pm\bm{\gamma}}$、$\mathbf{W}$）比 $C^2$ 表更吃寄存器–$C \times d_k = 64\times128$ 是 $C^2 = 64^2$ 的两倍。

三条出路：

1. $e^{-\bm{\gamma}}$ 不常驻：它只在构造 $\mathbf{M}$、$\mathbf{A}^{qk}$ 时用，用完即弃，可以放 shared memory。省 64 reg。
2. $\mathbf{W}$ 走 shared：它是 $\mathbf{U} - \mathbf{W}\mathbf{S}$ 的中间量，不参与后续 GEMM 的累加器。省 64 reg。
3. $C$ 降到 32：所有 $C$ 相关量减半。

组合 1+2 后约 192 reg，可行。下面的 kernel 采用这个方案。

4.1 声明与累积衰减

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

@tilelang.jit(out_idx=[5])
def kda_chunk(B, H, S, DK, DV, blk=64, block_DV=32,
              num_stages=2, threads=128,
              dtype=T.float16, accum_dtype=T.float32):
    C  = blk
    NS = T.ceildiv(S, C)

    @T.prim_func
    def main(
        Q:    T.Tensor([B, S, H, DK], dtype),
        K:    T.Tensor([B, S, H, DK], dtype),
        V:    T.Tensor([B, S, H, DV], dtype),
        G:    T.Tensor([B, S, H, DK], accum_dtype),   # 逐通道 log 门控，< 0
        Beta: T.Tensor([B, S, H], accum_dtype),
        O:    T.Tensor([B, S, H, DV], dtype),
    ):
        with T.Kernel(T.ceildiv(DV, block_DV), B * H, threads=threads) as (bv, bbh):
            bb, bh = bbh // H, bbh % H
            dv0 = bv * block_DV

            Q_s  = T.alloc_shared([C, DK], dtype)
            K_s  = T.alloc_shared([C, DK], dtype)
            V_s  = T.alloc_shared([C, block_DV], dtype)
            S_s  = T.alloc_shared([DK, block_DV], dtype)
            O_s  = T.alloc_shared([C, block_DV], dtype)
            # 逐通道特有：e^{±γ} 与 W 放 shared，避免寄存器超限
            Ep_s = T.alloc_shared([C, DK], accum_dtype)   # e^{γ_r}
            Em_s = T.alloc_shared([C, DK], accum_dtype)   # e^{-γ_r}
            Kp_s = T.alloc_shared([C, DK], dtype)         # k ⊙ e^{γ}
            Km_s = T.alloc_shared([C, DK], dtype)         # k ⊙ e^{-γ}
            W_s  = T.alloc_shared([C, DK], accum_dtype)
            Vt_s = T.alloc_shared([C, block_DV], dtype)

            S_f    = T.alloc_fragment([DK, block_DV], accum_dtype)  # loop-carried
            acc_o  = T.alloc_fragment([C, block_DV], accum_dtype)
            Aqk    = T.alloc_fragment([C, C], accum_dtype)
            Aq_c   = T.alloc_fragment([C, C], dtype)
            Lmat   = T.alloc_fragment([C, C], accum_dtype)          # 三角系数
            Vt     = T.alloc_fragment([C, block_DV], accum_dtype)
            RHS_v  = T.alloc_fragment([C, block_DV], accum_dtype)
            gam    = T.alloc_fragment([C, DK], accum_dtype)         # 累积 log
            be_f   = T.alloc_fragment([C], accum_dtype)

累积衰减是逐通道的前缀和–$d_k$ 个通道各自独立 cumsum，通道之间完全并行：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

T.clear(S_f)

for i_s in T.Pipelined(NS, num_stages=num_stages):
    s0 = i_s * C
    T.copy(Q[bb, s0:s0+C, bh, :], Q_s)
    T.copy(K[bb, s0:s0+C, bh, :], K_s)
    T.copy(V[bb, s0:s0+C, bh, dv0:dv0+block_DV], V_s)

    # ── 逐通道累积 log 衰减：通道维并行，位置维串行 ──
    for d in T.Parallel(DK):
        gam[0, d] = G[bb, s0, bh, d]
    for r in T.serial(1, C):
        for d in T.Parallel(DK):
            gam[r, d] = gam[r - 1, d] + G[bb, s0 + r, bh, d]

    for r, d in T.Parallel(C, DK):
        Ep_s[r, d] = T.exp2(gam[r, d] * 1.4426950408889634)   # e^{γ}
        Em_s[r, d] = T.exp2(-gam[r, d] * 1.4426950408889634)  # e^{-γ}
        Kp_s[r, d] = T.Cast(dtype, K_s[r, d] * Ep_s[r, d])
        Km_s[r, d] = T.Cast(dtype, K_s[r, d] * Em_s[r, d])
    for r in T.Parallel(C):
        be_f[r] = Beta[bb, s0 + r, bh]

对比 GDN：那里前缀和是 $C$ 步串行、每步 1 个数；这里是 $C$ 步串行、每步 $d_k$ 个通道并行。串行长度不变，但并行度从 1 涨到 128–逐通道门在这一处反而更适合 GPU。

1.4426950408889634 是 $\log_2 e$，把 $e^x$ 转成硬件 exp2（与第三篇一致）。

4.2 指数分离：两次 GEMM 取代双重循环

1
2
3
4
5
6
7
8
9
10
11
12

# ── M = strictLower(Kp @ Km^T)：衰减 KK 矩阵 ──
T.gemm(Kp_s, Km_s, Lmat, transpose_B=True, clear_accum=True)
for i, j in T.Parallel(C, C):
    Lmat[i, j] = T.if_then_else(
        j < i, be_f[j] * Lmat[i, j], 0.0)      # L[c,i] = β_i M[c,i]

# ── A^qk = tril(Qp @ Km^T)，含对角 ──
for r, d in T.Parallel(C, DK):
    Q_s[r, d] = T.Cast(dtype, Q_s[r, d] * Ep_s[r, d])   # q ⊙ e^{γ}
T.gemm(Q_s, Km_s, Aqk, transpose_B=True, clear_accum=True)
for i, j in T.Parallel(C, C):
    Aqk[i, j] = T.if_then_else(j <= i, Aqk[i, j], 0.0)

这是全篇最关键的两行 GEMM。 朴素写法需要 $C^2/2$ 次长度 $d_k$ 的加权内积（$C=64, d_k=128$ 时约 26 万次 exp）；分离后 $e^{\pm\bm{\gamma}}$ 各算一次（$C \times d_k = 8192$ 次 exp2），剩下交给 Tensor Core。

注意两个掩码的边界不同：Lmat 取 $j < i$（严格下三角，对角会重复计入 $\beta_r$），Aqk 取 $j \le i$（含对角，自己看自己零衰减）。

4.3 前向替换：W 与 U 合并求解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

# ── 右端项：[Kp | V]，注意 β 不在这里乘 ──
for r, d in T.Parallel(C, DK):
    W_s[r, d] = K_s[r, d] * Ep_s[r, d]
for r, d in T.Parallel(C, block_DV):
    RHS_v[r, d] = V_s[r, d]

# ── 前向替换，一次解出 T@[Kp | V] ──
for r in T.serial(C):
    for j in T.serial(r):
        for d in T.Parallel(DK):
            W_s[r, d] -= Lmat[r, j] * W_s[j, d]
        for d in T.Parallel(block_DV):
            RHS_v[r, d] -= Lmat[r, j] * RHS_v[j, d]

# ── 解完之后才逐行乘 β（Â = diag(β)T）──
for r, d in T.Parallel(C, DK):
    W_s[r, d] *= be_f[r]
for r, d in T.Parallel(C, block_DV):
    RHS_v[r, d] *= be_f[r]

# ── 伪值 Ṽ = U - W S ──
T.copy(S_f, S_s)
T.gemm(W_s, S_s, Vt, clear_accum=True)      # W S
for r, d in T.Parallel(C, block_DV):
    Vt[r, d] = RHS_v[r, d] - Vt[r, d]
T.copy(Vt, Vt_s)

$\mathbf{W}$ 与 $\mathbf{U}$ 共用同一个 $(\mathbf{I}+\mathbf{L})$，拼在一次前向替换里只付一遍串行代价。这一步的并行度是 $d_k + \text{block}_{DV} = 160$，比 GDN 的 32 好得多–逐通道门在这里又一次因为多了通道维而受益。

4.4 输出与状态更新

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# ② 块间 + ③ 块内
T.gemm(Q_s, S_s, acc_o, clear_accum=True)   # (Q ⊙ e^γ) S
T.copy(Aqk, Aq_c)
T.gemm(Aq_c, Vt_s, acc_o)                   # A^qk Ṽ

T.copy(acc_o, O_s)
T.copy(O_s, O[bb, s0:s0+C, bh, dv0:dv0+block_DV])

# ① 状态更新：Diag(e^{γ_C}) S + (K ⊙ e^{γ_C-γ_r})^T Ṽ
for d, j in T.Parallel(DK, block_DV):
    S_f[d, j] *= Ep_s[C - 1, d]             # 逐通道！不是标量
for r, d in T.Parallel(C, DK):
    Kp_s[r, d] = T.Cast(dtype,
        K_s[r, d] * Ep_s[C - 1, d] * Em_s[r, d])   # e^{γ_C-γ_r}
T.gemm(Kp_s, Vt_s, S_f, transpose_A=True)

最后一段体现了逐通道门与前三篇最直观的差别：

• GDN：S_f[i,j] *= alpha_pow_C，一个标量乘整个状态
• KDA：S_f[d,j] *= Ep_s[C-1,d]，每个 $d_k$ 行乘各自的衰减

状态的不同行按不同速率淡出–这就是"逐通道"在 kernel 层面的全部含义。

e^{\bm{\gamma}_C - \bm{\gamma}_r} 用 Ep_s[C-1,d] * Em_s[r,d] 算，即 $e^{\gamma_C}\cdot e^{-\gamma_r}$。这里复用了已有的两张表，不必重算 exp2；但它也是 $e^{-\bm{\gamma}}$ 参与的第三处，进一步说明为什么下界不可省。

4.5 七处次序约束

比 GDN 多两处：

1. 状态更新排在输出写回之后（右移语义，四篇一致）。
2. S_f *= Ep_s[C-1,:] 在 T.gemm 之前。
3. Q_s 被原地乘 $e^{\bm{\gamma}}$ 后不可再用于块内项–但本篇的块内项恰好也用 $\bm{q}\odot e^{\bm{\gamma}}$（$\mathbf{A}^{qk}$ 的定义里就带），所以不需要重载 Q。这是与 GDN 的一处反差，容易照抄出错。
4. Lmat 只取严格下三角（$j < i$），Aqk 含对角（$j \le i$）。
5. K_s 必须保持原始值：$\mathbf{Kp}$、$\mathbf{Km}$、右端项、状态更新四处都从 K_s 派生，任何一处原地修改都会污染后续。本篇的做法是始终写入独立的 Kp_s/Km_s，K_s 只读–比 GDN 的"污染两次再重载"更清晰。
6. $\beta$ 在前向替换之后乘，不能提前混进右端项。$\hat{\mathbf{A}} = \operatorname{diag}(\beta)\mathbf{T}$ 与 $\mathbf{T}\operatorname{diag}(\beta)$ 是两种不同配对，各自要搭配 $L_{ci} = \beta_i M_{ci}$ 或 $\beta_c M_{ci}$。$\beta$ 全相等时这个错误完全隐身。
7. T.clear(S_f) 在循环外，clear_accum=True 在循环内。

第 3、5、6 条都是"照抄上一篇会错"的地方，其中第 6 条我实际踩过。

5. 四篇对照