# 稀疏求解单元终点：从填充账本到多层误差证书

## 任务与验收

本任务使用精确有理数核对小型算法，结构填充另用图计算。交付不只是一份最终近似解，而是让别人能逐步复算的状态记录。

1. 对五节点链消去节点 1、3，求解保留系统并恢复原变量
2. 对六节点图及七乘七网格执行完整符号消元，记录新增边、列模式与工作量
3. 检查 IC(0) 的成功例与正定输入失败例，区分原系统和移位预条件器
4. 对九节点链计算重叠 Schwarz 校正及一个全局粗方向的效果
5. 在七节点 Poisson 系统上完整执行一次 V-cycle，逐层核对矩阵、残差和回程校正
6. 对八节点聚合计算平滑延拓、粗矩阵及候选再现缺陷

成功条件：最终所有恢复解满足原方程；结构账本与数值相消分开；每个粗右端来自残差；Galerkin 尺度一致；真实残差和能量误差均能核验；若使用 PCG，应用算子固定、线性、对称正定。

## 一、静态凝聚

本题所有节点编号从零开始。令 $T_n=\operatorname{tridiag}(-1,2,-1)$，求解 $T_5x=\mathbf1$。消去 $I=(1,3)$，保留 $B=(0,2,4)$。内部方程为

$$
x_1=(1+x_0+x_2)/2,\qquad x_3=(1+x_2+x_4)/2.
$$

代回得

$$
S=\begin{pmatrix}3/2&-1/2&0\\-1/2&1&-1/2\\0&-1/2&3/2\end{pmatrix},
\qquad g=(3/2,2,3/2)^T.
$$

求解 $Sx_B=g$ 得 $(5/2,9/2,5/2)^T$，恢复内部后得到

$$
x=(5/2,4,9/2,4,5/2)^T.
$$

原残差严格为零。直接删行列则得到保留变量均为 $1/2$，已经删除了内部传递的耦合和载荷。

正定性证书也能独立给出：对任意界面值 $y$，令内部值 $z=-A_{II}^{-1}A_{IB}y$，则 $y^TSy=(z,y)^TA(z,y)>0$。逆矩阵在算法中以内部系统求解实现。

## 二、符号消元与排列

### 六节点的两份账本

原边为 $01,12,23,34,45,50,14$。自然序每步新增边依次为

$$
\{15\},\quad\{24,25\},\quad\{35\},\quad\varnothing,\quad\varnothing,\quad\varnothing.
$$

对应列行号为

$$
(0,1,5),(1,2,4,5),(2,3,4,5),(3,4,5),(4,5),(5).
$$

于是因子结构项数为 17，CSC 列指针为 $(0,3,7,11,14,16,17)$。各列后继数 $d=(2,3,3,2,1,0)$，故

$$
\sum(d_j+1)^2=55,\qquad \sum d_j(d_j+1)/2=19.
$$

改用原顶点序 $(0,2,5,3,1,4)$，只新增 $15,13$ 两条边。排列后的列行号为

$$
(0,2,4),(1,3,4),(2,4,5),(3,4,5),(4,5),(5).
$$

结构项数为 15，平方和为 41，对称乘减更新数为 13。原顶点号和排列后矩阵位置号要同时保存，不能混用。

路径判据认证所有位置：$i>j$ 时，$L_{ij}$ 出现，当且仅当原图有从 $j$ 到 $i$ 的路径，内部节点均早于 $j$。补团归纳分别证明两个方向。脚本穷举全部 1,024 个五顶点无向图，每图测试自然序与反序，共 2,048 份账本，逐位置交叉检查路径判据与消元树祖先性质。

### 七乘七嵌套剖分

坐标 $(r,c)$ 的原编号为 $7r+c$。每个方块先递归左上、右上、左下、右下四个角区，再按自然编号追加中间行列。最外十字有 13 点，各角区是三乘三。

自然序有 349 个因子结构位置，十字嵌套剖分有 288 个。原矩阵下三角含对角为 $49+84=133$ 项，因此新增填充分别为 216 与 155。统一平方和工作指标为 2643 与 1926。

深度 $\ell$ 的四次分区有 $4^\ell$ 个边长约 $q/2^\ell$ 的子域，每个前沿含 $O(q/2^\ell)$ 个变量。该层存储 $O(q^2)$，算术 $O(q^3/2^\ell)$；总计二维 $O(N\log N)$ 存储、$O(N^{3/2})$ 算术。这个证明依赖规则局部二维结构，不能只说“递归”就推广到任意图。

完整 49 节点顺序及逐列模式见结果文件的 nested_dissection 字段；同一脚本另算三、十五、三十一阶边长，验证统计口径一致。

## 三、不完全分解的正确接口

在三乘三网格的五点矩阵上，IC(0) 第一步消去 0，保留对角更新，丢弃新交叉位置 $(3,1)$。精确更新应为 $-1/4$，不完全因子固定该位置为零；重构后 $(A-M)_{31}=-1/4$。

使用单位下三角因子 $\widehat L$ 与对角 $D$，九个主元为

$$
4,15/4,56/15,15/4,52/15,2507/728,56/15,2507/728,8572/2507.
$$

全部为正，因此应用 $\widehat L^{-T}D^{-1}\widehat L^{-1}$ 是固定正定映射。

对失败例

$$
A=\begin{pmatrix}1&3/5&0&3/5\\3/5&1&3/5&0\\0&3/5&1&-3/5\\3/5&0&-3/5&1\end{pmatrix},
$$

写 $A=I+(3/5)C$，验证 $C^2=2I$，得到全部特征值 $1\pm3\sqrt2/5>0$。精确分解主元为 $(1,16/25,7/16,7/25)$，IC(0) 却在最后得到 $-32/175$，必须报告失败。对 $A+I/5$ 设置，主元为 $(6/5,9/10,4/5,9/20)$；这次成功，但外层依旧求解原系统。

结果文件还记录普通 CG 与 IC-PCG 的每步原残差。对选定右端，IC-PCG 在四步达到 $10^{-4}$ 绝对阈值，普通 CG 需五步；精确零残差终止步数反而分别是七和五。验收应采用真正需要的精度与总工作，不能只报一种有利的迭代数。

## 四、加性 Schwarz 与粗方向

取 $T_9$，子域为 $0\ldots5$ 与 $3\ldots8$。对全一残差，两局部解均为 $(3,5,6,6,5,3)^T$。延伸相加得到

$$
Br=(3,5,6,9,10,9,6,5,3)^T.
$$

覆盖保证 $r^TBr=\sum(R_ir)^TA_i^{-1}(R_ir)>0$，故 $B$ 正定。每个子域用的是同一个残差，不能在循环中悄悄换成更新后的残差。

令 $p=(1,2,3,4,5,4,3,2,1)^T$。有 $Ap=2e_4$、$p^TAp=10$，加入粗项 $B_0=pp^T/10$。从误差 $p$ 出发，同样以步长 $1/2$ 应用单层与两层算子：

$$
e_1=(5,10,15,15,15,15,15,10,5)^T/7,
$$

$$
e_2=(3,6,9,2,-5,2,9,6,3)^T/14.
$$

能量分别从 10 变为 $150/49$ 和 $125/98$。三个 $A$-正交投影之和的最大特征值不超过 3，故步长 $1/2$ 在 $(0,2/3)$ 内，有收敛依据。若用步长一，局部项和粗项可能对同一方向重复校正；不能因为 $B_0Ap=p$ 就声称完整加性方法一步精确。

## 五、七节点完整 V-cycle

### 层级与输入

细矩阵 $A_0=T_7$，右端 $b_0=(1,0,0,0,0,0,1)^T$，初值零。精确解 $x_*=\mathbf1_7$。延拓为

$$
P_0=\begin{pmatrix}1/2&0&0\\1&0&0\\1/2&1/2&0\\0&1&0\\0&1/2&1/2\\0&0&1\\0&0&1/2\end{pmatrix},
\qquad P_1=(1/2,1,1/2)^T.
$$

每次限制取转置，故

$$
A_1=P_0^TA_0P_0=\tfrac12T_3,
\qquad A_2=P_1^TA_1P_1=(1/2).
$$

每层前后各做一次加权 Jacobi，$\omega=2/3$；最粗层精确求解。

### 下降过程

| 层 | 本层右端 | 预平滑结果 | 预平滑后残差 | 限制右端 |
|---|---|---|---|---|
| 7 节点 | $(1,0,0,0,0,0,1)$ | $(1/3,0,0,0,0,0,1/3)$ | $(1/3,1/3,0,0,0,1/3,1/3)$ | $(1/2,0,1/2)$ |
| 3 节点 | $(1/2,0,1/2)$ | $(1/3,0,1/3)$ | $(1/6,1/3,1/6)$ | $(1/2)$ |
| 1 节点 | $(1/2)$ | 精确解 $z=1$ | 零 | 无 |

注意三节点层对角是 1，而最细层对角是 2。两次预平滑都出现 $1/3$，来自不同矩阵与右端的共同缩放，不是把一份模板输出机械复制两次。

### 回升过程

一节点解延拓为 $(1/2,1,1/2)$，加到三节点预平滑解后得到 $(5/6,1,5/6)$。其残差为 $(1/6,-1/6,1/6)$；乘以本层 Jacobi 权重 $2/3$，得到后平滑返回值

$$
z_1=(17/18,8/9,17/18)^T.
$$

其剩余粗残差是 $(0,1/18,0)^T$，并非精确粗解。将 $z_1$ 延拓到七节点：

$$
P_0z_1=(17/36,17/18,11/12,8/9,11/12,17/18,17/36)^T.
$$

加回最细层预平滑值：

$$
x_{\mathrm{corr}}=(29/36,17/18,11/12,8/9,11/12,17/18,29/36)^T.
$$

再用原细右端做一次后平滑，得到

$$
x_{\mathrm{out}}=(11/12,8/9,11/12,49/54,11/12,8/9,11/12)^T.
$$

原残差与误差分别为

$$
r=(1/18,1/18,-1/27,1/54,-1/27,1/18,1/18)^T,
$$

$$
e=(1/12,1/9,1/12,5/54,1/12,1/9,1/12)^T.
$$

满足 $A_0e=r$，且

$$
\|e\|_{A_0}^2=25/1458,\qquad \|\mathbf1_7\|_{A_0}^2=2.
$$

### 误差传播与 PCG 证书

令 $W_\ell=(2/3)D_\ell^{-1}$、$S_\ell=I-W_\ell A_\ell$。若下层循环应用是 $B_{\ell+1}$，本层误差矩阵为

$$
E_\ell=S_\ell\left(I-P_\ell B_{\ell+1}P_\ell^TA_\ell\right)S_\ell.
$$

脚本独立构造这个递归矩阵，并检查它恰等于 $I-B_\ell A_\ell$，且作用于初始全一误差得到上述 $e$。七节点和十五节点均通过。

对称配对平滑还给出

$$
B_\ell=2W_\ell-W_\ell A_\ell W_\ell
+(I-W_\ell A_\ell)P_\ell B_{\ell+1}P_\ell^T(I-A_\ell W_\ell).
$$

第一项在 $0<\omega<2/\lambda_{\max}(D^{-1}A_\ell)$ 时正定，第二项半正定。最粗层为正定逆，逐层归纳得到固定 SPD 应用。脚本还用标准基调用构造实际小矩阵，检查精确对称和所有 LDL 主元为正。

该证书允许标准 PCG 调用，并不单凭 SPD 就证明任意 V-cycle 作为独立迭代都具有网格无关收敛率。

## 六、八节点平滑聚合

对 $T_8$ 按相邻两点聚合，$P_0$ 为四个聚合指示向量，$P_0\mathbf1_4=\mathbf1_8$。平滑 $P=(I-A/3)P_0$ 得

$$
P=\frac13\begin{pmatrix}2&0&0&0\\2&1&0&0\\1&2&0&0\\0&2&1&0\\0&1&2&0\\0&0&2&1\\0&0&1&2\\0&0&0&2\end{pmatrix}.
$$

由 $A_c=P^TAP$ 得

$$
A_c=\frac19\begin{pmatrix}6&-1&-1&0\\-1&4&-1&-1\\-1&-1&4&-1\\0&-1&-1&6\end{pmatrix}.
$$

四个平滑粗基的能量为 $2/3,4/9,4/9,2/3$，小于未平滑时各自的 2。细、粗算子非零数分别 22、14，延拓从 8 项增加到 14 项，两层算子复杂度 $18/11$、网格复杂度 $3/2$。

候选缺陷从恒等式而不是外观判断：

$$
\mathbf1_8-P\mathbf1_4=\frac13A\mathbf1_8
=(1/3,0,0,0,0,0,0,1/3)^T.
$$

Dirichlet 常量并非真零模。未平滑空间能精确消除常量，而平滑空间对同一误差剩下 $(1,-1,0,1,1,0,-1,1)^T/7$，能量为 $2/7$。这说明降低基函数能量与改善每一个特定向量的粗逼近是不同命题。

## 七、运行与结果解释

在任意目录放入下载脚本，执行：

```sh
python3 foundation-sparse-solvers-capstone.py result.json
```

只需 Python 标准库。脚本使用 Fraction 对线性系统、因子、Schwarz、两网格和 V-cycle 进行精确检查；平滑正弦曲线单独使用双精度三角函数，不把其近似小数记成精确分数。

结果包含逐轮新增边、列模式、完整嵌套剖分次序、IC 与精确因子、CG/PCG 原残差历史、每层循环状态和所有粗矩阵。小例子的稠密有理矩阵运算只承担核验职责，不冒称是达到正文稀疏复杂度的生产实现。