数学推导与算法

本节推导正文中“直接给出”的几个结论，给出度保持随机化零模型的算法，并补充格兰杰因果、加权聚集系数与带符号模块度、扩散映射嵌入和指纹识别的数学细节。符号约定：$N$ 为脑区（节点）数，$T$ 为时间点数，$x_{it}$ 为脑区 $i$ 在时刻 $t$ 的 BOLD 信号，$\Sigma$ 为 $N\times N$ 协方差矩阵，$\Theta=\Sigma^{-1}$ 为精度矩阵。

A. 偏相关为什么等于精度矩阵的归一化负值

偏相关的“定义式”是残差相关：把 $x_i$ 和 $x_j$ 分别对其余 $N-2$ 个脑区做线性回归，取两组残差 $e_i, e_j$，再求它们的相关。下面说明这与 $-\Theta_{ij}/\sqrt{\Theta_{ii}\Theta_{jj}}$ 相同。

记 $A=\{i,j\}$，$B$ 为其余脑区。把协方差矩阵分块为 $\Sigma=\begin{pmatrix}\Sigma_{AA}&\Sigma_{AB}\\ \Sigma_{BA}&\Sigma_{BB}\end{pmatrix}$。$x_A$ 对 $x_B$ 回归后的残差协方差是舒尔补

$$\operatorname{Cov}(e_A)=\Sigma_{AA}-\Sigma_{AB}\Sigma_{BB}^{-1}\Sigma_{BA}.$$

分块矩阵求逆公式告诉我们，精度矩阵中对应的 $2\times 2$ 块恰好是这个舒尔补的逆：$\Theta_{AA}=\big(\operatorname{Cov}(e_A)\big)^{-1}$。于是

$$\operatorname{Cov}(e_A)=\Theta_{AA}^{-1}=\frac{1}{\Theta_{ii}\Theta_{jj}-\Theta_{ij}^2}\begin{pmatrix}\Theta_{jj}&-\Theta_{ij}\\-\Theta_{ij}&\Theta_{ii}\end{pmatrix},$$

对这个 $2\times2$ 协方差求相关，公共因子约去，得到

$$\rho_{ij\cdot\text{rest}}=\frac{-\Theta_{ij}}{\sqrt{\Theta_{ii}\Theta_{jj}}}.$$

这个推导只用到线性回归，不需要正态假设；正态假设的作用是让“偏相关为 0”进一步等价于“条件独立”（详见症状网络的数学推导）。实际问题是 $\Theta$ 要靠 $\hat\Sigma^{-1}$ 估计：当 $T$ 与 $N$ 同量级时 $\hat\Sigma$ 接近奇异，逆矩阵噪声极大；$T<N$ 时根本不可逆，这就是需要收缩估计（如 Ledoit–Wolf）或图形 LASSO 的原因。

B. Fisher $z$ 变换：方差稳定化

若 $(x_i,x_j)$ 为二元正态且 $T$ 个样本独立，样本相关 $r$ 的渐近方差为 $\operatorname{Var}(r)\approx(1-\rho^2)^2/T$——它依赖未知的真值 $\rho$，$\rho$ 越接近 $\pm1$ 方差越小，分布越偏。寻找变换 $g$ 使 $\operatorname{Var}(g(r))$ 不依赖 $\rho$。由 delta 方法 $\operatorname{Var}(g(r))\approx g'(\rho)^2\operatorname{Var}(r)$，令 $g'(\rho)=1/(1-\rho^2)$，积分得

$$g(\rho)=\tfrac12\ln\frac{1+\rho}{1-\rho}=\operatorname{artanh}(\rho),\qquad \operatorname{Var}(z)\approx\frac{1}{T-3}.$$

（$T-3$ 是小样本修正。）由此得到两种常用组间检验：

• 两独立相关的比较（如两组的组平均连接）：$Z=(z_1-z_2)\big/\sqrt{1/(n_1-3)+1/(n_2-3)}$，与标准正态比较；
• 被试水平检验（更常用）：每名被试每条边一个 $z$ 值，把它当作因变量，做两样本 $t$ 检验或纳入头动、年龄等协变量的一般线性模型。

C. 阈值方式为何会“制造”组间差异

用 ER 随机网络做一个粗略的量级估计：密度为 $D$ 时，聚集系数 $C\approx D$，平均路径长度 $L\approx\ln N/\ln\big((N-1)D\big)$。即使拓扑组织完全相同，只要密度不同，$C$ 和 $L$ 就会不同。

• 绝对阈值 $\tau$：密度 $D=\Pr(|r|>\tau)$ 由该被试整体 FC 强度决定。若患者组整体 FC 偏弱（或头动更大），则 $D$ 更低，$C$ 降、$L$ 升——这一“差异”完全来自密度。
• 比例阈值 $p$：强制 $D=p$。若某被试最强的前 $p$ 条边中有一部分接近噪声，这些边的位置近似随机，会使其网络 $C$ 降低、$L$ 也降低，向随机网络靠拢。

所以两种阈值各自引入不同的偏差，这正是阈值与加权中推荐“密度扫描 + 加权指标 + 把整体 FC 强度作为协变量”的理由。

D. 度保持随机化（Maslov–Sneppen）与小世界归一化

$C_{\text{rand}}$、$L_{\text{rand}}$ 和 $\phi_{\text{rand}}(k)$ 都需要一组与真实网络度序列相同的随机网络。最常用的是 Maslov & Sneppen（2002）的边交换算法：

算法 1  度保持边交换（无向、无权、简单图）
输入：边集 E（|E| = m），每边平均交换次数 Q（常用 10–100）
输出：度序列与原图相同的随机图 E'
1. E' ← E；建立哈希集合 H 以 O(1) 判断边是否存在
2. 重复 Q·m 次：
3.     从 E' 中均匀随机抽两条边 (a,b) 与 (c,d)
4.     以 1/2 概率把第二条边写成 (d,c)，以覆盖两种重连方式
5.     若 a,b,c,d 中有重复节点，或 (a,d)、(c,b) 已在 H 中：跳过
6.     删除 (a,b)、(c,d)；加入 (a,d)、(c,b)；同步更新 H
7. （可选）若要求连通：交换后检查连通性，不连通则撤销
8. 返回 E'
复杂度：时间 O(Q·m)（不检查连通性时），空间 O(m)

每次交换不改变 $a,b,c,d$ 四个节点的度，因此度序列严格保持。$Q$ 太小时随机图仍带有原图结构的“残留”，可检查 $C$ 随交换次数的变化是否已趋于平稳。加权网络可用保持强度近似不变的变体（BCT 中的 null_model_und_sign 等），对相关矩阵本身，还可以生成与真实相关矩阵元素均值、方差相同的随机相关矩阵作为零模型（Zalesky 等, 2012），它比边交换更能体现“相关网络天然高聚集”这一基线。

有了 $M$ 个（如 100 个）随机图，取其平均得到 $C_{\text{rand}}$、$L_{\text{rand}}$，计算 $\sigma=\dfrac{C/C_{\text{rand}}}{L/L_{\text{rand}}}$。$\sigma$ 的一个弱点是：稀疏随机网络的 $C_{\text{rand}}$ 很小，分母一小，$\sigma$ 就会随网络规模膨胀。Telesford 等（2011）因此提出

$$\omega=\frac{L_{\text{rand}}}{L}-\frac{C}{C_{\text{latt}}},$$

其中 $C_{\text{latt}}$ 是同度序列“格子化”网络（把边尽量重排到矩阵对角线附近，Sporns & Zwi, 2004）的聚集系数。$\omega\in[-1,1]$：接近 0 为小世界，接近 1 偏随机，接近 $-1$ 偏规则格子。它把 $C$ 与“同度序列下能达到的最大聚集”比较，不再依赖很小的 $C_{\text{rand}}$。

E. 格兰杰因果的向量自回归形式

对两个零均值时间序列 $x_t, y_t$，$p$ 阶受限模型与完整模型分别为

$$y_t=\sum_{k=1}^{p}a_k y_{t-k}+\varepsilon_t,\qquad y_t=\sum_{k=1}^{p}a'_k y_{t-k}+\sum_{k=1}^{p}b_k x_{t-k}+\varepsilon'_t .$$

用最小二乘拟合后，Geweke 形式的格兰杰因果量为

$$F_{x\to y}=\ln\frac{\operatorname{Var}(\hat\varepsilon_t)}{\operatorname{Var}(\hat\varepsilon'_t)}\ \ge 0 .$$

$F_{x\to y}=0$ 当且仅当 $b_1=\dots=b_p=0$（在总体意义上）。对 $T$ 个样本，原假设下 $T\cdot F_{x\to y}$ 渐近服从自由度为 $p$ 的 $\chi^2$ 分布，也可用嵌套模型的 $F$ 检验。阶数 $p$ 用 AIC/BIC 选择。多于两个脑区时，在两个模型中都加入其余脑区的过去值，就得到条件格兰杰因果，可以排除部分由第三方引起的虚假方向。

为什么 HRF 会反转方向？设 BOLD 为 $\tilde x=h_x*x$、$\tilde y=h_y*y$，若被驱动区 $y$ 的响应 $h_y$ 比 $h_x$ 提前 $\Delta$ 秒，而神经延迟 $\delta\ll\Delta$，则 $\tilde y$ 的变化在时间上反而领先 $\tilde x$ 约 $\Delta-\delta$，VAR 模型看到的是 $\tilde y$ 的过去帮助预测 $\tilde x$。GC 是关于观测信号的陈述，只有当观测过程对所有区域相同（或已被显式建模）时，才能把它解释为神经层面的方向。

F. 加权聚集系数与带符号模块度

Onnela 加权聚集系数的分子 $\sum_{j,h}(\hat w_{ij}\hat w_{ih}\hat w_{jh})^{1/3}$ 对节点 $i$ 参与的每个三角形取三边归一化权重的几何平均。等价的矩阵写法是：令 $\hat W^{1/3}$ 为逐元素开立方的矩阵，则分子为 $\big[(\hat W^{1/3})^3\big]_{ii}$，与二值情形 $C_i=(A^3)_{ii}/[k_i(k_i-1)]$ 结构相同，整网计算复杂度为一次矩阵乘法 $O(N^3)$（稀疏时更低）。

带符号模块度：把 $W$ 拆成正部 $W^{+}=\max(W,0)$ 与负部 $W^{-}=\max(-W,0)$，记 $v^{\pm}=\sum_{ij}w^{\pm}_{ij}$，$s_i^{\pm}=\sum_j w^{\pm}_{ij}$。分别定义

$$Q^{\pm}=\frac{1}{v^{\pm}}\sum_{ij}\Big(w^{\pm}_{ij}-\frac{s^{\pm}_i s^{\pm}_j}{v^{\pm}}\Big)\delta(c_i,c_j).$$

好的划分应使模块内正权多于期望（$Q^{+}$ 大）、负权少于期望（$Q^{-}$ 小）。Rubinov & Sporns（2011）的非对称版本为

$$Q^{*}=Q^{+}-\frac{v^{-}}{v^{+}+v^{-}}\,Q^{-},$$

负权部分按其占总权重的比例降权。当 $W$ 没有负权时 $Q^{*}$ 退化为普通加权模块度。

G. 扩散映射嵌入

给定对称非负亲和矩阵 $K$（如连接轮廓的余弦相似度），度为 $d_i=\sum_j K_{ij}$。先用参数 $\alpha\in[0,1]$ 去除采样密度的影响：

$$K^{(\alpha)}_{ij}=\frac{K_{ij}}{d_i^{\alpha}d_j^{\alpha}},\qquad d^{(\alpha)}_i=\sum_j K^{(\alpha)}_{ij},\qquad P_{ij}=\frac{K^{(\alpha)}_{ij}}{d^{(\alpha)}_i}.$$

$P$ 是图上随机游走的转移矩阵。其特征分解 $P\psi_\ell=\lambda_\ell\psi_\ell$ 中，$\lambda_0=1$ 对应常数向量（平凡解），其余特征值 $1>\lambda_1\ge\lambda_2\ge\dots$。扩散映射坐标为 $\Psi_t(i)=\big(\lambda_1^t\psi_1(i),\lambda_2^t\psi_2(i),\dots\big)$，$t$ 为扩散时间；BrainSpace 等工具常用对所有 $t$ 求和的形式，相当于把第 $\ell$ 个坐标乘以 $\lambda_\ell/(1-\lambda_\ell)$。在这一坐标中，两点的欧氏距离近似等于它们之间的“扩散距离”：若从 $i$ 和 $j$ 出发的随机游走在 $t$ 步后的分布相似，两者就靠近。

数值上，$P$ 不对称，但它与对称矩阵 $S=D_\alpha^{-1/2}K^{(\alpha)}D_\alpha^{-1/2}$ 相似（$D_\alpha=\operatorname{diag}(d^{(\alpha)})$）：$S$ 的特征向量 $\phi_\ell$ 与 $\psi_\ell=D_\alpha^{-1/2}\phi_\ell$ 对应，特征值相同。因此可以用对称特征分解（如 numpy.linalg.eigh）稳定求解，复杂度 $O(N^3)$。$\alpha=0$ 对应经典的归一化图拉普拉斯谱嵌入，$\alpha=1$ 则在流形假设下近似拉普拉斯–贝尔特拉米算子，与采样密度无关；$\alpha=0.5$ 介于两者之间。

H. 指纹识别的准确率与可识别性

设 $\mathbf a_i,\mathbf b_i\in\mathbb R^{m}$ 为被试 $i$ 两次扫描的边向量（各自标准化），识别矩阵 $C_{ij}=\operatorname{corr}(\mathbf b_i,\mathbf a_j)$。识别准确率 $\text{Acc}=\frac{1}{n}\sum_i \mathbb 1\big[\arg\max_j C_{ij}=i\big]$，$I_{\text{self}}=\frac1n\sum_i C_{ii}$，$I_{\text{others}}=\frac{1}{n(n-1)}\sum_{i\ne j}C_{ij}$，$I_{\text{diff}}=I_{\text{self}}-I_{\text{others}}$。

一个简单的加性模型可以说明各因素的作用：设 $\mathbf a_i=\mathbf g+\mathbf u_i+\mathbf e_i$，$\mathbf b_i=\mathbf g+\mathbf u_i+\mathbf e'_i$，其中 $\mathbf g$ 是群体共同模式，$\mathbf u_i$ 是方差 $\sigma_u^2$ 的个体特异成分，$\mathbf e$ 是方差 $\sigma_e^2$ 的扫描噪声，各成分独立。则同一人与不同人的协方差之差为 $\sigma_u^2$，而相关的分母随 $\sigma_e^2$ 增大，因此 $I_{\text{diff}}\propto\sigma_u^2/(\sigma_g^2+\sigma_u^2+\sigma_e^2)$。另一方面，$C_{ii}$ 与 $C_{ij}$ 的估计误差随边数 $m$ 以 $1/\sqrt m$ 缩小，所以识别准确率同时取决于个体信号比例和使用的边数。扫描时间越长，$\sigma_e^2$ 越小，这解释了为什么更多的个体内数据能提高可识别性。随机置换被试标签可以得到准确率的零分布（期望值为 $1/n$）。