寻找控制复杂网络的全部可能输入节点的高效算法

一句话看懂

一个网络的所有最小输入节点集（MIS）的并集，就是全部可能输入节点。这篇论文证明：只要算出一个最大匹配，再从它的输入节点出发沿交替路径搜索，路径能到达的节点就是全部可能输入节点。这一步恰好可以并入 Hopcroft–Karp 算法的最后一轮，因此总复杂度与求一次最大匹配相同，为 $O(N^{1/2}L)$；在大型真实网络上比已有方法快几个数量级。

研究背景与问题

结构可控性与最大匹配中讲过：驱动节点（本文称输入节点）由最大匹配中的非匹配节点给出，最大匹配不唯一，所以 MIS 不唯一。知道”哪些节点有资格当输入节点”在应用中很有用：理解节点在控制中的角色、在约束条件下设计 MIS、在信号通路中识别关键基因。

已有方法（Jia 等 2013，Nature Communications）先求一个最大匹配，再逐个检验匹配节点：对每个节点判断它能否在某个最大匹配中成为非匹配节点。求最大匹配需 $O(N^{1/2}L)$，逐点检验需 $O(NL)$，总复杂度 $O(NL)$。对几十万、上百万节点的网络，这意味着数小时甚至十几小时的计算。

Zhang、Lv 与 Pu（2016）的输入图论文已经证明：从输入节点出发沿控制路径可达的节点都是可能输入节点，反之不可达者为冗余节点。本文的问题是：能否不构建完整输入图、不逐点检验，以求一次最大匹配的代价拿到全部可能输入节点？

核心思想与方法

二分图与交替路径

把有向网络 $G(V,E)$ 拆成二分图：每个节点 $n$ 拆成入端 $n_{in}$ 和出端 $n_{out}$，有向边 $u\to v$ 对应 $u_{out}$–$v_{in}$ 的无向边。最大匹配下，$V_{in}$ 中的非匹配节点就是输入节点。交替路径是匹配边与非匹配边交替出现的路径；两端都是非匹配节点的交替路径称为增广路径，Berge 定理说最大匹配等价于不存在增广路径。

充分性的证明很像”翻转”：$m_{in}$ 和 $n_{in}$ 都在入端，所以连接它们的交替路径长度为偶数，以 $m_{in}$ 处的非匹配边开始、以 $n_{in}$ 处的匹配边结束。把路径上所有边的匹配状态翻转，匹配边数不变，得到新的最大匹配 $M'$；此时 $m_{in}$ 被匹配，$n_{in}$ 失去匹配，成为输入节点。必要性用反证：若 $n$ 在另一个最大匹配 $M'$ 中未被匹配，可以构造一条偶数长度的路径，回到 $M$ 下会发现它恰好从 $M$ 的某个输入节点出发到达 $n_{in}$，与”不可达”矛盾。

小例子：边为 $1\to2$、$1\to3$、$4\to3$、$4\to5$，取最大匹配 $\{1\to2, 4\to3\}$，输入节点为 1、4、5（1 和 4 入度为 0）。从 $5_{in}$ 出发：非匹配边到 $4_{out}$，匹配边到 $3_{in}$，非匹配边到 $1_{out}$，匹配边到 $2_{in}$。这条交替路径经过的入端是 $3_{in}$、$2_{in}$，所以 3 和 2 都是可能输入节点，也都是 5 的替代节点：MIS $\{1,4,5\}$ 可以换成 $\{1,4,3\}$ 或 $\{1,4,2\}$。全部可能输入节点为 $\{1,2,3,4,5\}$。

把搜索并入 Hopcroft–Karp 算法

Hopcroft–Karp 算法反复做两件事：从所有非匹配节点出发搜索交替路径，若找到增广路径就沿其扩大匹配。当某一轮再也找不到增广路径时，算法结束——而这最后一轮恰好已经搜索了从所有输入节点出发的交替路径。根据定理 1，这一轮经过的入端节点就是答案。算法 All_Input(G)：

1. 由有向网络构建二分图，初始匹配为空；
2. 从 $V_{in}$ 中所有非匹配节点出发找交替路径，把路径上的入端节点记为候选；
3. 若存在增广路径，扩大匹配、清空候选，回到第 2 步；
4. 若不存在增广路径，候选节点即全部可能输入节点，当前非匹配节点构成一个 MIS。

由于只是在原算法中”顺手记录”候选节点，复杂度保持 $O(N^{1/2}L)$。同时，每个输入节点的交替路径可达集就是它的替代节点集。

主要发现

1. 合成网络上提速数百到数千倍。13 个无标度网络（节点数 $10^5$ 到 $5\times10^6$，平均度 6 到 16），新算法与旧方法结果一致。$n=10^5$ 时加速 52 倍；$n=5\times10^6$ 时新算法仅用 3.276 秒，旧方法约 24 013 秒，加速约 7330 倍。
2. 网络越稠密，加速越明显。$n=10^6$ 时，平均度从 6 增至 16，加速比从 315 倍升至约 4097 倍。
3. 大型真实网络上提速可达 $10^4$ 量级。在 20 个真实网络上，Amazon0601（403 394 个节点）用时 2.886 秒，旧方法 49 293 秒，加速约 17 080 倍；Higgs_Twitter（456 626 个节点、约 1486 万条边）用时 12.589 秒，旧方法约 66 445 秒。
4. 生物网络中的可能输入节点比例差异很大。人类有向蛋白相互作用网络（6339 个节点、34 814 条边）的可能输入节点比例为 0.585，TRN-Yeast-1 转录网络为 0.999，E. coli 转录网络为 0.730。
5. 不同输入节点的可替代性不同。St. Marks 食物网的一个 MIS 有 13 个输入节点，各自的替代节点集大小不同，说明有些输入节点在控制中更”稳健”、更容易被替换；Kohonen、TRN-Yeast-1、Ythan 等网络结果类似。部分输入节点替代集大小相同，是因为它们连到同一组替代节点。

怎样读这篇论文

• 全文很短，核心是”定理 1 + 推论 1 + 算法 All_Input”。建议先看图 2（二分图与交替路径）和图 3（算法在小网络上的完整过程），再读证明。
• 读证明时自己画一个二分图，标出入端和出端，确认”为什么连接两个入端节点的交替路径长度一定是偶数”。
• 表 1、表 2 是运行时间对比，抓住”规模越大、越稠密，优势越大”即可；图 5、图 6 关于替代节点，体现了算法的第二个用途。
• 若已读过 2016 年输入图论文，可以把本文的交替路径与那里的”控制路径”一一对应起来理解。

局限与后续工作

• 理论框架仍是线性时不变系统的结构可控性和最大匹配框架，只适用于有向网络；无向网络常用的最小支配集框架不在本文讨论范围内。
• 推论 1 只保证单个输入节点可被替换；多个输入节点同时替换时可能冲突（2016 年论文已指出共享邻居的输入节点不能同时替换）。
• 运行时间比较依赖具体实现（JAVA 实现、单台工作站），加速倍数应理解为量级而非精确值。
• 与系列工作的关系：
  • 本文是 2016 年输入图的算法化：输入图回答了”为什么”（控制邻接、交换定理、控制分支），本文回答了”如何快速算”，把找出全部可能输入节点的复杂度从 $O(NL)$ 降到 $O(N^{1/2}L)$。
  • 快速判定”谁是可能输入节点、谁是冗余节点”，是后续改变控制模式工作的计算基础：Zhang 与 Li（2019）通过删边、Zhang 等（2021）通过反转边，切断驱动节点与巨输入分支之间的交替路径，使大量可能输入节点变为冗余节点（2021 年论文在参考文献中引用了本文）。
  • 后续研究基于控制路径给节点分类、寻找控制核心点（见驱动节点的枚举、分类与自适应控制），其中”所有可能的头节点即所有可能的驱动节点”这一步正可用本文方法高效计算。

关联知识

• 结构可控性与最大匹配中的二分图拆点与 Hopcroft–Karp 算法是本文的直接基础；驱动节点的枚举、分类与自适应控制提到的”在最大匹配算法最后一步的交替路径搜索中输出全部可能的驱动节点，复杂度 $O(\sqrt N L)$“即本文结果。
• 可在网络可控性与稳定性的代码示例中用 NetworkX 求出最大匹配，再从非匹配入端做交替路径的广度优先搜索，即可实现本文算法。
• 在药物靶点研究中，替代节点集给出了”某靶点不可用时的备选靶点”，而且只改动一条匹配链，对原控制方案影响最小。

延伸思考

1. 在小例子中，若再加一条边 $2\to5$，重新求最大匹配和全部可能输入节点。哪些节点的类型发生了变化？
2. 为什么旧方法需要”逐点检验”，而本文只需一轮交替路径搜索？请用”翻转交替路径得到新最大匹配”的思想解释。
3. 某输入节点的替代节点集很大，另一个很小。在设计药物干预方案时，你会如何利用这一信息？
4. 人类蛋白相互作用网络的可能输入节点比例为 0.585。这个数字本身能否说明网络处于分布式控制模式？还需要哪些信息？