1. 拓扑信号处理从图到高阶网络的范式演进在传统的图信号处理领域我们习惯于将数据点视为网络中的节点将关系视为边并利用图拉普拉斯算子对节点上的信号进行滤波和分析。然而现实世界的数据结构远比简单的“点-边”关系复杂。想象一下大脑中的神经元集群、社交网络中的群体互动或者海洋中洋流的路径这些现象不仅涉及成对的关系还涉及三个、四个甚至更多元素之间的高阶交互。为了刻画这种复杂性单纯复形和高阶网络应运而生它们允许我们定义在节点、边、三角形等不同维度“单形”上的信号。拓扑信号处理正是为了处理这些定义在复杂拓扑结构上的高阶数据而生。其核心挑战在于如何联合分析不同维度如节点和边上的信号并捕捉它们之间由拓扑结构所决定的深层耦合关系。传统的基于Hodge拉普拉斯算子的方法虽然强大但通常将不同维度的信号分开处理或者通过升维/降维操作进行间接耦合这可能会损失掉维度间直接的、物理意义明确的关联信息。这时狄拉克算子及其对应的狄拉克方程为我们打开了一扇新的大门。狄拉克方程最初是描述相对论性量子粒子如电子行为的核心方程其离散版本——拓扑狄拉克方程天然地将节点和边信号耦合在一个统一的“旋量”框架中。基于此发展出的狄拉克方程信号处理算法特别是DESP及其迭代版本IDESP提供了一种物理启发的、自适应的滤波新思路能够同时利用节点和边的信息并引入一个可学习的“质量”参数来调节信号尺度在处理真实世界的高阶拓扑信号时展现出了显著优势。2. 核心基石拓扑狄拉克算子与方程要理解DESP和IDESP我们必须先深入其数学基础拓扑狄拉克算子及其对应的特征值问题。2.1 从图到单纯复形高阶信号的舞台首先我们明确操作对象。一个单纯复形是由点0-单形、线段1-单形、三角形2-单形等基本几何单元按照一定规则如下面的单形其所有面也必须存在于复形中组合而成的结构。一个网络可以看作是一个1维单纯复形只有节点和边。在单纯复形上我们可以定义k-链群其元素是k-单形的形式线性组合。边界算子 ∂_k 将一个k-单形映射到其边界构成的(k-1)-链。上边界算子则是边界算子的对偶。对于信号处理我们关注的是定义在这些单形上的实值函数即k-上链可以看作k-维信号。例如节点上的电位、边上的流量、三角形上的涡度等。2.2 狄拉克算子耦合节点与边的桥梁Hodge拉普拉斯算子 L_k ∂_{k1} ∂_{k1}^T ∂_k^T ∂_k 是处理k-维信号的标准工具但它作用在单一维度上。狄拉克算子 D 的精妙之处在于它提供了一个统一的框架。在1维单纯复形即图的语境下拓扑狄拉克算子可以构造为一个分块矩阵[ D \begin{pmatrix} 0 B_1^T \ B_1 0 \end{pmatrix} ]其中B_1 是关联矩阵或其某种加权/归一化变体它编码了节点与边之间的关联关系。B_1^T 是其转置。这个算子直接作用在一个由节点信号和边信号拼接而成的“拓扑旋量” ψ [χ; φ] 上其中 χ 是节点信号向量φ 是边信号向量。狄拉克算子的平方 D^2 会产生一个分块对角矩阵其对角线上的块正是Hodge拉普拉斯算子L_0 B_1^T B_1 节点上的图拉普拉斯和 L_1 B_1 B_1^T 边上的上拉普拉斯。这揭示了狄拉克算子与经典工具的内在联系。2.3 拓扑狄拉克方程与质量参数受量子力学中狄拉克方程的启发我们可以定义拓扑空间上的狄拉克方程[ D \psi E \psi ]其中 E 是特征值类比能量ψ 是特征旋量。这个方程要求旋量 ψ 同时是狄拉克算子的特征向量。然而一个更一般、更强大的形式是引入了“质量”项 mγ[ (D m \gamma) \psi E \psi ]这里γ 是一个满足 γ^2 I 且与 D 反对易的矩阵在典型构造中γ diag(I, -I)。质量参数 m 的引入具有深刻的物理和数学意义物理类比在物理学中质量项破坏了手征对称性使得粒子获得质量。在信号处理中它打破了节点信号和边信号在能量尺度上的对称性允许算法自适应地调整对不同维度信号的“重视程度”。数学功能方程 (D mγ) ψ E ψ 可以改写为 (D - E I) ψ -m γ ψ。这表明质量项 m 实质上在节点和边信号分量之间引入了一种耦合或“混合”其强度由 m 控制。当 m0 时方程退化为 D ψ E ψ节点和边信号通过 B_1 和 B_1^T 耦合当 m 不为零时额外的 γ 项引入了直接的尺度变换耦合。谱特性可以证明满足此方程的旋量 ψ其分量 χ 和 φ 分别满足 (L_0 - E^2 I m^2 I) χ 0 和 (L_1 - E^2 I m^2 I) φ 0 等关系。这意味着特征值 λ 满足相对论性色散关系E^2 λ^2 m^2其中 λ^2 是旋量在 L_0 或 L_1 上的某种期望值。这个关系是DESP算法进行参数优化的核心依据。注意质量参数 m 并非事先已知的它是DESP算法需要从带噪数据中学习的关键参数。正是这个可学习的 m使得DESP能够灵活地适应信号在不同维度上的真实尺度关系从而超越了固定 m0 的DSP算法。3. DESP算法原理深度拆解狄拉克方程信号处理算法的目标是从一个被噪声污染的观测拓扑旋量 ˜ψ ψ ϵ 中恢复出真实的信号 ψ。其核心假设是真实信号 ψ 是拓扑狄拉克方程 (D mγ) ψ E ψ 的一个特征旋量。3.1 问题建模与损失函数设计给定观测信号 ˜ψDESP试图寻找一个估计信号 ˆψ_m 和一个质量参数 m使得 ˆψ_m 尽可能接近真实信号 ψ同时尽可能满足狄拉克方程的特征结构。算法通过最小化以下损失函数来实现[ L_m(\hat{\psi}_m) | \tilde{\psi} - \hat{\psi}_m |^2 \tau \cdot R(\hat{\psi}_m) ]这是一个经典的“数据保真项正则化项”的范式。数据保真项‖ ˜ψ - ˆψ_m ‖^2迫使重建信号 ˆψ_m 接近观测数据 ˜ψ。正则化项R(ˆψ_m)注入先验知识即信号应具有狄拉克方程特征向量的结构。DESP采用的正则化项是‖ (D mγ) ˆψ_m ‖^2其目的是让 ˆψ_m 在算子 (D mγ) 的作用下尽可能接近于零向量即对应特征值 E0 的情况或者说让 ˆψ_m 尽可能位于 (D mγ) 的零空间或低能量子空间附近。参数 τ 控制着正则化的强度。因此完整的损失函数为[ L_m(\hat{\psi}_m) | \tilde{\psi} - \hat{\psi}_m |^2 \tau | (D m\gamma) \hat{\psi}_m |^2 ]对于固定的 m这是一个关于 ˆψ_m 的凸优化问题。通过求导并令导数为零我们可以得到闭合形式的解[ \hat{\psi}_m (I \tau (D m\gamma)^T (D m\gamma))^{-1} \tilde{\psi} ]这个解可以看作是对观测信号 ˜ψ 应用了一个依赖于 m 和 τ 的线性滤波器。3.2 质量参数 m 的优化两种策略算法的关键步骤是确定最优的质量参数 m。DESP提供了两种基于不同原理的优化策略。策略一最小化损失函数 L_m最直接的方法是遍历一个合理的 m 值范围对于每个 m计算其对应的最优重建信号 ˆψ_m 以及损失值 L_m(ˆψ_m)然后选择使 L_m 最小的 m 作为最优估计。这个过程可以直观地理解为寻找一个 m使得在该 m 定义的正则化项下得到的总损失数据误差结构误差最小。策略二最小化相对论性色散关系误差这是一种更具物理洞察力的方法。根据拓扑狄拉克方程一个真正的特征旋量 ψ 必须满足相对论性色散关系E^2 λ^2 m^2。对于任意一个重建信号 ˆψ_m 和假定的 m我们可以计算其“表现出的”能量 E_m 和拉普拉斯期望值 λ_m^2[ E_m \frac{\hat{\psi}_m^T H \hat{\psi}_m}{|\hat{\psi}_m|^2}, \quad \lambda_m^2 \frac{\hat{\psi}_m^T D^2 \hat{\psi}_m}{|\hat{\psi}_m|^2} ]其中 H 是与 (Dmγ) 相关的哈密顿量矩阵。然后定义色散关系误差[ S_m | E_m^2 - (\lambda_m^2 m^2) | ]S_m 衡量了重建信号 ˆψ_m 在多大程度上偏离了一个理想狄拉克特征旋量应满足的数学关系。最优的 m 就是使 S_m 最小的那个值。当 ˆψ_m 恰好是一个特征旋量时S_m 0。实操心得在噪声水平不高的情况下两种策略得到的 m 估计值和最终重建误差通常非常接近。策略一最小化L_m更通用直接与重建目标挂钩策略二最小化S_m物理意义更清晰有时在噪声结构特殊时可能更稳健。在实际应用中可以两种方法都尝试对比结果。如果计算资源允许将 S_m 作为另一个正则化项加入 L_m 进行联合优化也是一个值得探索的方向。3.3 与LSP和DSP的统合关系DESP的强大之处在于它的泛化能力当设置质量参数 m 0 且能量 E 0 时DESP的正则化项退化为τ ‖ D ˆψ ‖^2。这恰好是狄拉克信号处理算法所使用的正则化器。因此DSP是DESP在 m0 时的一个特例。当进一步地我们只考虑节点信号或者忽略狄拉克算子的耦合部分DSP可以退化为经典的拉普拉斯信号处理算法其正则化项是τ χ^T L_0 χ对节点信号。因此LSP也可以视为DESP的一个退化特例。这种关系意味着DESP提供了一个统一的框架。当真实信号确实是一个 m0 的狄拉克特征态时DESP能够通过学习发现 m≈0从而自动退化为性能相当的DSP。但当真实信号具有非零质量即节点和边信号满足更一般的尺度关系时DESP通过学习到非零的 m能够获得比DSP和LSP更好的重建效果。论文中的数值实验清晰地表明随着 |m| 的增大DESP相对于DSP的改进也越来越显著。3.4 联合处理的优势噪声不均场景下的信息互补DESP天然地联合处理节点和边信号。一个有趣且实用的优势体现在节点和边噪声水平不同的场景。假设节点信号噪声大α1高而边信号相对干净α2低。传统的、单独处理节点信号的方法会受限于高噪声。但DESP在联合优化时干净的边信号可以通过狄拉克算子 D 中的 B_1 和 B_1^T 矩阵为重建节点信号提供额外的约束和信息流从而改善节点信号的重建精度。反之亦然干净的节点信号也有助于重建噪声较大的边信号。这揭示了拓扑信号处理的一个核心价值不同维度信号之间通过拓扑算子相互约束可以提升整体去噪和重建的鲁棒性。这在高阶数据中尤为宝贵因为我们可能在某些维度拥有更高质量的测量数据。4. IDESP算法从单特征态到一般信号的迭代扩展DESP有一个基本假设真实信号是拓扑狄拉克方程的单个特征旋量。然而现实数据更为复杂通常是多个特征旋量的线性组合。迭代狄拉克方程信号处理算法正是为突破这一限制而设计。4.1 算法流程与迭代思想IDESP的核心思想是残差分解。算法流程如下输入带噪观测信号 ˜ψ估计的或真实的信噪比系数变异 c_V^{true}。初始化设残差信号 φ ˜ψ迭代次数 J1重建信号 ˆΨ 0。主循环 a. 对当前残差 φ 运行DESP算法得到当前轮次估计的主特征旋量 ˆψ_J。 b. 将 ˆψ_J 累加到总重建信号中ˆΨ ˆΨ ˆψ_J。 c. 更新残差φ ˜ψ - ˆΨ。这移除了已重建的信号成分。 d. 计算当前总重建信号 ˆΨ 与原始观测信号 ˜ψ 之间的系数变异 c_V(J) ‖ ˆΨ - ˜ψ ‖ / ‖ ˆΨ ‖。 e. 判断是否停止如果 |c_V(J) - c_V^{true}| 开始增大即本次迭代后重建信号的估计噪声水平反而离真实噪声水平更远了则停止迭代。否则JJ1回到步骤a。输出最终的重建信号 ˆΨ。这个过程的直观理解是第一轮DESP从带噪数据中提取出能量最强或最显著的狄拉克特征模式 ˆψ_1。将其从观测数据中减去后残差中可能包含能量次强的特征模式以及噪声。对残差再次应用DESP可以提取出下一个主要模式 ˆψ_2。如此反复逐步剥离出信号中的多个主要成分。4.2 停止准则的重要性防止过拟合噪声迭代算法最大的风险是“过拟合”——将噪声也当作信号成分提取出来。IDESP巧妙地利用系数变异噪声与信号的能量比作为停止准则。c_V^{true}这是观测数据中真实的噪声水平通常需要事先估计例如通过信号平稳区域的统计或先验知识。c_V(J)这是经过J轮迭代后当前重建信号 ˆΨ 与原始观测 ˜ψ 之间的差异视为残余“噪声”与重建信号能量的比值。如果重建是完美的c_V(J) 应等于 c_V^{true}。因此迭代的目标是让 c_V(J) 逼近 c_V^{true}。当 J 过小时c_V(J) c_V^{true}说明信号成分提取不足当 J 过大时c_V(J) c_V^{true}说明已经开始拟合噪声。IDESP选择在 |c_V(J) - c_V^{true}| 达到最小时停止此时重建信号包含了绝大部分真实信号成分而引入了最少的噪声成分。注意事项准确估计 c_V^{true} 是关键。如果估计值偏差太大会导致迭代过早停止或过度进行。在实际应用中可以采用多种方法交叉验证1) 如果数据有重复测量或已知的干净片段可直接计算2) 利用高频部分估计噪声功率3) 将IDESP本身用于参数搜索观察重建信号随迭代次数变化的平台区。4.3 性能验证与真实数据应用论文在海洋浮标轨迹数据上验证了IDESP。该数据将太平洋马达加斯加岛附近的海域划分为网格节点浮标在不同网格间的运动形成边上的流量信号。从边信号出发他们构造了一个合成的节点-边联合拓扑旋量作为真实信号并添加噪声。实验表明误差下降随着迭代次数 J 增加真实误差 ∆(J) ‖ ˆΨ - ψ ‖ 首先快速下降在最优迭代次数 J_opt 附近达到最低点之后由于开始拟合噪声误差缓慢上升。这验证了停止准则的有效性。超越IDSP对比同样迭代框架但使用DSP即固定m0的IDSP算法IDESP无论是用L_m还是S_m优化m能够以更少的迭代次数达到相同或更低的误差水平。这说明学习到的质量参数 m 有效地捕捉了信号的内在结构提高了每次迭代提取特征的效率。可视化改善从噪声严重的原始观测信号到第一次迭代J1提取出主要环流模式再到第二次迭代J2进一步细化结构IDESP逐步重建出了与真实流场高度相似的模式。5. 算法实现关键与实战指南5.1 数据结构与算子构建实现DESP/IDESP的一步是构建拓扑狄拉克算子 D 和 γ 矩阵。import numpy as np import scipy.sparse as sp def build_topological_dirac_operator(B1): 构建1维单纯复形图的拓扑狄拉克算子。 参数: B1: scipy.sparse.csr_matrix, 形状为 (n_nodes, n_edges)。 节点-边关联矩阵。B1[i, j] 1 如果边j以节点i为终点-1如果为起点有向图 对于无向图通常取绝对值或特定方向。 返回: D: scipy.sparse.csr_matrix, 形状为 (n_nodesn_edges, n_nodesn_edges)。 狄拉克算子 [[0, B1^T], [B1, 0]]。 gamma: scipy.sparse.csr_matrix, 相同形状。 γ 矩阵通常为 diag(I_{nodes}, -I_{edges})。 n_nodes, n_edges B1.shape # 构建狄拉克算子 D zero_nn sp.csr_matrix((n_nodes, n_nodes)) zero_ne sp.csr_matrix((n_nodes, n_edges)) zero_en sp.csr_matrix((n_edges, n_nodes)) zero_ee sp.csr_matrix((n_edges, n_edges)) # D [[0, B1^T], [B1, 0]] top_block sp.hstack([zero_nn, B1.T], formatcsr) bottom_block sp.hstack([B1, zero_ee], formatcsr) D sp.vstack([top_block, bottom_block], formatcsr) # 构建 gamma 矩阵: diag(I, -I) gamma_n sp.identity(n_nodes, formatcsr) gamma_e -sp.identity(n_edges, formatcsr) gamma sp.block_diag([gamma_n, gamma_e], formatcsr) return D, gamma5.2 DESP核心求解与参数扫描对于固定的 τ 和 m求解 ˆψ_m 需要解一个线性系统。由于矩阵 (I τ (Dmγ)^T (Dmγ)) 是大规模稀疏正定矩阵应使用迭代法如共轭梯度法CG而非直接求逆。from scipy.sparse.linalg import LinearOperator, cg def desp_single_step(psi_tilde, D, gamma, m, tau, tol1e-8): 执行单次DESP重建固定m。 参数: psi_tilde: np.ndarray, 观测拓扑旋量 (n_nodes n_edges)。 D, gamma: 狄拉克算子和gamma矩阵。 m: float, 质量参数。 tau: float, 正则化参数。 tol: float, CG求解器的容忍误差。 返回: psi_hat: np.ndarray, 重建的拓扑旋量。 n_total D.shape[0] # 构造算子 A I tau * (D m*gamma)^T (D m*gamma) H D m * gamma # 哈密顿量 H D mγ A_operator LinearOperator((n_total, n_total), matveclambda x: x tau * H.T (H x)) # 使用共轭梯度法求解 A * psi_hat psi_tilde psi_hat, info cg(A_operator, psi_tilde, toltol) if info ! 0: print(fCG solver did not converge, info: {info}) return psi_hat def optimize_m_by_loss(psi_tilde, D, gamma, tau, m_range): 通过扫描m范围并最小化损失函数L_m来优化m。 参数: m_range: np.ndarray, 待搜索的m值数组。 返回: best_m: float, 最优m值。 best_psi: np.ndarray, 对应最优m的重建信号。 losses: list, 各m对应的损失值。 best_loss float(inf) best_m None best_psi None losses [] for m in m_range: psi_hat desp_single_step(psi_tilde, D, gamma, m, tau) # 计算损失 L_m ||psi_tilde - psi_hat||^2 tau * ||(Dm*gamma)psi_hat||^2 residual psi_tilde - psi_hat data_fidelity np.linalg.norm(residual)**2 H_psi (D m*gamma) psi_hat regularization tau * np.linalg.norm(H_psi)**2 loss data_fidelity regularization losses.append(loss) if loss best_loss: best_loss loss best_m m best_psi psi_hat.copy() return best_m, best_psi, losses5.3 IDESP迭代实现IDESP的实现需要仔细管理残差和停止准则。def idesp(psi_tilde, D, gamma, tau, c_v_true, max_iter20, m_optimizerloss): IDESP算法主函数。 参数: psi_tilde: 观测信号。 c_v_true: float, 真实的系数变异噪声水平估计。 max_iter: 最大迭代次数。 m_optimizer: loss 或 rdre优化m的策略。 返回: Psi_hat: 最终重建信号。 m_list: 每次迭代学习到的m值列表。 c_v_list: 每次迭代的c_V(J)值列表。 errors: 每次迭代的真实误差如果有真实信号psi_gt。 n_total len(psi_tilde) residual psi_tilde.copy() Psi_hat np.zeros(n_total) m_list [] c_v_list [] prev_diff float(inf) for j in range(1, max_iter1): # 1. 对当前残差应用DESP优化m if m_optimizer loss: m_opt, psi_hat_j, _ optimize_m_by_loss(residual, D, gamma, tau, m_rangenp.linspace(-2, 2, 41)) else: # rdre m_opt, psi_hat_j, _ optimize_m_by_rdre(residual, D, gamma, tau, m_rangenp.linspace(-2, 2, 41)) m_list.append(m_opt) # 2. 更新总重建信号和残差 Psi_hat psi_hat_j residual psi_tilde - Psi_hat # 3. 计算当前系数变异 c_V(J) norm_Psi np.linalg.norm(Psi_hat) if norm_Psi 1e-12: c_v 1.0 else: c_v np.linalg.norm(residual) / norm_Psi c_v_list.append(c_v) # 4. 检查停止准则 (简单版本当c_V开始上升或变化很小时停止) diff abs(c_v - c_v_true) if diff prev_diff: # 本次迭代后离目标更远了 print(fStopping at iteration {j}, c_V starts to diverge.) # 回退到上一次迭代的结果 Psi_hat - psi_hat_j m_list.pop() c_v_list.pop() break prev_diff diff # 可选如果c_V已经非常接近真实值也可停止 if diff 0.01 * c_v_true: print(fStopping at iteration {j}, c_V is close enough to target.) break return Psi_hat, m_list, c_v_list5.4 参数选择与调优经验正则化参数 τ控制平滑度与数据保真度的权衡。τ 过大重建信号过于平滑可能丢失细节τ 过小去噪效果弱。建议从 τ 10 开始尝试根据信噪比调整。高噪声下可适当增大 τ。质量参数 m 的搜索范围需要根据先验知识或数据尺度设定。对于标准化后的信号m 通常在 [-3, 3] 区间内。可以先进行粗搜索如步长0.5在最优值附近再进行细搜索。真实噪声水平 c_V^{true} 的估计这是IDESP的关键。如果无法直接获得可以使用高通滤波或小波变换估计噪声方差。在数据中寻找理论上应为“平坦”或“平滑”的区域计算该区域的方差作为噪声估计。将 c_V^{true} 作为一个超参数在验证集上调整选择使最终应用指标如分类准确率、预测误差最优的值。最大迭代次数 max_iter不宜设置过大通常5-10次迭代足以提取主要模式。设置一个上限防止无限循环并通过监控 c_V(J) 曲线来观察收敛情况。6. 应用场景与未来展望DESP/IDESP算法为处理具有内在拓扑结构的高维数据提供了新工具。计算神经科学脑功能网络中的信号不仅存在于脑区节点也存在于连接边上。DESP可以联合处理节点上的BOLD信号时间序列和边上的功能连接强度动态可能更精准地识别与认知任务相关的脑网络模式。交通流与物流分析在交通网络中节点代表路口或城市边代表道路。节点信号可以是车流量边信号可以是行程时间或拥堵指数。IDESP可以分解出影响全局交通的主要传播模式如早晚高峰潮汐流。物理场重构在传感器网络中节点测量标量场如温度边可以测量梯度或通量。DESP利用狄拉克方程耦合标量和梯度信息能在部分传感器失效或噪声较大时实现更鲁棒的物理场重构。分子动力学与生物信息学蛋白质残基节点和它们之间的相互作用边构成网络。节点信号可以是氨基酸特性边信号可以是结合能。拓扑狄拉克分析可能帮助识别稳定的构象或功能位点。未来方向与挑战扩展到更高阶当前框架主要针对节点和边1-复形。如何将其推广到包含三角形、四面体等更高阶单形的复形处理更丰富的高阶相互作用是一个前沿问题。与图神经网络结合狄拉克算子可以作为GNN中的消息传递算子其固有的处理节点和边特征的能力可能缓解传统GNN的过度平滑和过度挤压问题。DESP的滤波思想可以融入图卷积层的设计。处理时变拓扑许多现实网络的拓扑结构是动态变化的。发展适用于时变单纯复形的在线DESP算法是一个实际需求。理论收敛性分析目前IDESP的停止准则虽直观有效但缺乏严格的理论收敛性保证。在什么条件下算法能稳定地恢复信号的主要成分仍需深入研究。在我自己的尝试中将IDESP应用于社交网络用户-群组联合分析时发现初始的关联矩阵 B1 的构造方式对结果影响巨大。例如在定义“用户-边”关系时是采用简单的隶属关系还是考虑交互频率的加权会导致学到的“质量参数” m 具有完全不同的解释。一个负的 m 值可能暗示着节点和边信号之间存在某种抑制性关联这为数据解读提供了新的视角。这提醒我们算法背后的拓扑模型必须紧密结合领域知识来构建否则“物理启发”就可能沦为“数学游戏”。
拓扑信号处理进阶：狄拉克方程与IDESP算法解析
发布时间：2026/5/24 6:47:48
1. 拓扑信号处理从图到高阶网络的范式演进在传统的图信号处理领域我们习惯于将数据点视为网络中的节点将关系视为边并利用图拉普拉斯算子对节点上的信号进行滤波和分析。然而现实世界的数据结构远比简单的“点-边”关系复杂。想象一下大脑中的神经元集群、社交网络中的群体互动或者海洋中洋流的路径这些现象不仅涉及成对的关系还涉及三个、四个甚至更多元素之间的高阶交互。为了刻画这种复杂性单纯复形和高阶网络应运而生它们允许我们定义在节点、边、三角形等不同维度“单形”上的信号。拓扑信号处理正是为了处理这些定义在复杂拓扑结构上的高阶数据而生。其核心挑战在于如何联合分析不同维度如节点和边上的信号并捕捉它们之间由拓扑结构所决定的深层耦合关系。传统的基于Hodge拉普拉斯算子的方法虽然强大但通常将不同维度的信号分开处理或者通过升维/降维操作进行间接耦合这可能会损失掉维度间直接的、物理意义明确的关联信息。这时狄拉克算子及其对应的狄拉克方程为我们打开了一扇新的大门。狄拉克方程最初是描述相对论性量子粒子如电子行为的核心方程其离散版本——拓扑狄拉克方程天然地将节点和边信号耦合在一个统一的“旋量”框架中。基于此发展出的狄拉克方程信号处理算法特别是DESP及其迭代版本IDESP提供了一种物理启发的、自适应的滤波新思路能够同时利用节点和边的信息并引入一个可学习的“质量”参数来调节信号尺度在处理真实世界的高阶拓扑信号时展现出了显著优势。2. 核心基石拓扑狄拉克算子与方程要理解DESP和IDESP我们必须先深入其数学基础拓扑狄拉克算子及其对应的特征值问题。2.1 从图到单纯复形高阶信号的舞台首先我们明确操作对象。一个单纯复形是由点0-单形、线段1-单形、三角形2-单形等基本几何单元按照一定规则如下面的单形其所有面也必须存在于复形中组合而成的结构。一个网络可以看作是一个1维单纯复形只有节点和边。在单纯复形上我们可以定义k-链群其元素是k-单形的形式线性组合。边界算子 ∂_k 将一个k-单形映射到其边界构成的(k-1)-链。上边界算子则是边界算子的对偶。对于信号处理我们关注的是定义在这些单形上的实值函数即k-上链可以看作k-维信号。例如节点上的电位、边上的流量、三角形上的涡度等。2.2 狄拉克算子耦合节点与边的桥梁Hodge拉普拉斯算子 L_k ∂_{k1} ∂_{k1}^T ∂_k^T ∂_k 是处理k-维信号的标准工具但它作用在单一维度上。狄拉克算子 D 的精妙之处在于它提供了一个统一的框架。在1维单纯复形即图的语境下拓扑狄拉克算子可以构造为一个分块矩阵[ D \begin{pmatrix} 0 B_1^T \ B_1 0 \end{pmatrix} ]其中B_1 是关联矩阵或其某种加权/归一化变体它编码了节点与边之间的关联关系。B_1^T 是其转置。这个算子直接作用在一个由节点信号和边信号拼接而成的“拓扑旋量” ψ [χ; φ] 上其中 χ 是节点信号向量φ 是边信号向量。狄拉克算子的平方 D^2 会产生一个分块对角矩阵其对角线上的块正是Hodge拉普拉斯算子L_0 B_1^T B_1 节点上的图拉普拉斯和 L_1 B_1 B_1^T 边上的上拉普拉斯。这揭示了狄拉克算子与经典工具的内在联系。2.3 拓扑狄拉克方程与质量参数受量子力学中狄拉克方程的启发我们可以定义拓扑空间上的狄拉克方程[ D \psi E \psi ]其中 E 是特征值类比能量ψ 是特征旋量。这个方程要求旋量 ψ 同时是狄拉克算子的特征向量。然而一个更一般、更强大的形式是引入了“质量”项 mγ[ (D m \gamma) \psi E \psi ]这里γ 是一个满足 γ^2 I 且与 D 反对易的矩阵在典型构造中γ diag(I, -I)。质量参数 m 的引入具有深刻的物理和数学意义物理类比在物理学中质量项破坏了手征对称性使得粒子获得质量。在信号处理中它打破了节点信号和边信号在能量尺度上的对称性允许算法自适应地调整对不同维度信号的“重视程度”。数学功能方程 (D mγ) ψ E ψ 可以改写为 (D - E I) ψ -m γ ψ。这表明质量项 m 实质上在节点和边信号分量之间引入了一种耦合或“混合”其强度由 m 控制。当 m0 时方程退化为 D ψ E ψ节点和边信号通过 B_1 和 B_1^T 耦合当 m 不为零时额外的 γ 项引入了直接的尺度变换耦合。谱特性可以证明满足此方程的旋量 ψ其分量 χ 和 φ 分别满足 (L_0 - E^2 I m^2 I) χ 0 和 (L_1 - E^2 I m^2 I) φ 0 等关系。这意味着特征值 λ 满足相对论性色散关系E^2 λ^2 m^2其中 λ^2 是旋量在 L_0 或 L_1 上的某种期望值。这个关系是DESP算法进行参数优化的核心依据。注意质量参数 m 并非事先已知的它是DESP算法需要从带噪数据中学习的关键参数。正是这个可学习的 m使得DESP能够灵活地适应信号在不同维度上的真实尺度关系从而超越了固定 m0 的DSP算法。3. DESP算法原理深度拆解狄拉克方程信号处理算法的目标是从一个被噪声污染的观测拓扑旋量 ˜ψ ψ ϵ 中恢复出真实的信号 ψ。其核心假设是真实信号 ψ 是拓扑狄拉克方程 (D mγ) ψ E ψ 的一个特征旋量。3.1 问题建模与损失函数设计给定观测信号 ˜ψDESP试图寻找一个估计信号 ˆψ_m 和一个质量参数 m使得 ˆψ_m 尽可能接近真实信号 ψ同时尽可能满足狄拉克方程的特征结构。算法通过最小化以下损失函数来实现[ L_m(\hat{\psi}_m) | \tilde{\psi} - \hat{\psi}_m |^2 \tau \cdot R(\hat{\psi}_m) ]这是一个经典的“数据保真项正则化项”的范式。数据保真项‖ ˜ψ - ˆψ_m ‖^2迫使重建信号 ˆψ_m 接近观测数据 ˜ψ。正则化项R(ˆψ_m)注入先验知识即信号应具有狄拉克方程特征向量的结构。DESP采用的正则化项是‖ (D mγ) ˆψ_m ‖^2其目的是让 ˆψ_m 在算子 (D mγ) 的作用下尽可能接近于零向量即对应特征值 E0 的情况或者说让 ˆψ_m 尽可能位于 (D mγ) 的零空间或低能量子空间附近。参数 τ 控制着正则化的强度。因此完整的损失函数为[ L_m(\hat{\psi}_m) | \tilde{\psi} - \hat{\psi}_m |^2 \tau | (D m\gamma) \hat{\psi}_m |^2 ]对于固定的 m这是一个关于 ˆψ_m 的凸优化问题。通过求导并令导数为零我们可以得到闭合形式的解[ \hat{\psi}_m (I \tau (D m\gamma)^T (D m\gamma))^{-1} \tilde{\psi} ]这个解可以看作是对观测信号 ˜ψ 应用了一个依赖于 m 和 τ 的线性滤波器。3.2 质量参数 m 的优化两种策略算法的关键步骤是确定最优的质量参数 m。DESP提供了两种基于不同原理的优化策略。策略一最小化损失函数 L_m最直接的方法是遍历一个合理的 m 值范围对于每个 m计算其对应的最优重建信号 ˆψ_m 以及损失值 L_m(ˆψ_m)然后选择使 L_m 最小的 m 作为最优估计。这个过程可以直观地理解为寻找一个 m使得在该 m 定义的正则化项下得到的总损失数据误差结构误差最小。策略二最小化相对论性色散关系误差这是一种更具物理洞察力的方法。根据拓扑狄拉克方程一个真正的特征旋量 ψ 必须满足相对论性色散关系E^2 λ^2 m^2。对于任意一个重建信号 ˆψ_m 和假定的 m我们可以计算其“表现出的”能量 E_m 和拉普拉斯期望值 λ_m^2[ E_m \frac{\hat{\psi}_m^T H \hat{\psi}_m}{|\hat{\psi}_m|^2}, \quad \lambda_m^2 \frac{\hat{\psi}_m^T D^2 \hat{\psi}_m}{|\hat{\psi}_m|^2} ]其中 H 是与 (Dmγ) 相关的哈密顿量矩阵。然后定义色散关系误差[ S_m | E_m^2 - (\lambda_m^2 m^2) | ]S_m 衡量了重建信号 ˆψ_m 在多大程度上偏离了一个理想狄拉克特征旋量应满足的数学关系。最优的 m 就是使 S_m 最小的那个值。当 ˆψ_m 恰好是一个特征旋量时S_m 0。实操心得在噪声水平不高的情况下两种策略得到的 m 估计值和最终重建误差通常非常接近。策略一最小化L_m更通用直接与重建目标挂钩策略二最小化S_m物理意义更清晰有时在噪声结构特殊时可能更稳健。在实际应用中可以两种方法都尝试对比结果。如果计算资源允许将 S_m 作为另一个正则化项加入 L_m 进行联合优化也是一个值得探索的方向。3.3 与LSP和DSP的统合关系DESP的强大之处在于它的泛化能力当设置质量参数 m 0 且能量 E 0 时DESP的正则化项退化为τ ‖ D ˆψ ‖^2。这恰好是狄拉克信号处理算法所使用的正则化器。因此DSP是DESP在 m0 时的一个特例。当进一步地我们只考虑节点信号或者忽略狄拉克算子的耦合部分DSP可以退化为经典的拉普拉斯信号处理算法其正则化项是τ χ^T L_0 χ对节点信号。因此LSP也可以视为DESP的一个退化特例。这种关系意味着DESP提供了一个统一的框架。当真实信号确实是一个 m0 的狄拉克特征态时DESP能够通过学习发现 m≈0从而自动退化为性能相当的DSP。但当真实信号具有非零质量即节点和边信号满足更一般的尺度关系时DESP通过学习到非零的 m能够获得比DSP和LSP更好的重建效果。论文中的数值实验清晰地表明随着 |m| 的增大DESP相对于DSP的改进也越来越显著。3.4 联合处理的优势噪声不均场景下的信息互补DESP天然地联合处理节点和边信号。一个有趣且实用的优势体现在节点和边噪声水平不同的场景。假设节点信号噪声大α1高而边信号相对干净α2低。传统的、单独处理节点信号的方法会受限于高噪声。但DESP在联合优化时干净的边信号可以通过狄拉克算子 D 中的 B_1 和 B_1^T 矩阵为重建节点信号提供额外的约束和信息流从而改善节点信号的重建精度。反之亦然干净的节点信号也有助于重建噪声较大的边信号。这揭示了拓扑信号处理的一个核心价值不同维度信号之间通过拓扑算子相互约束可以提升整体去噪和重建的鲁棒性。这在高阶数据中尤为宝贵因为我们可能在某些维度拥有更高质量的测量数据。4. IDESP算法从单特征态到一般信号的迭代扩展DESP有一个基本假设真实信号是拓扑狄拉克方程的单个特征旋量。然而现实数据更为复杂通常是多个特征旋量的线性组合。迭代狄拉克方程信号处理算法正是为突破这一限制而设计。4.1 算法流程与迭代思想IDESP的核心思想是残差分解。算法流程如下输入带噪观测信号 ˜ψ估计的或真实的信噪比系数变异 c_V^{true}。初始化设残差信号 φ ˜ψ迭代次数 J1重建信号 ˆΨ 0。主循环 a. 对当前残差 φ 运行DESP算法得到当前轮次估计的主特征旋量 ˆψ_J。 b. 将 ˆψ_J 累加到总重建信号中ˆΨ ˆΨ ˆψ_J。 c. 更新残差φ ˜ψ - ˆΨ。这移除了已重建的信号成分。 d. 计算当前总重建信号 ˆΨ 与原始观测信号 ˜ψ 之间的系数变异 c_V(J) ‖ ˆΨ - ˜ψ ‖ / ‖ ˆΨ ‖。 e. 判断是否停止如果 |c_V(J) - c_V^{true}| 开始增大即本次迭代后重建信号的估计噪声水平反而离真实噪声水平更远了则停止迭代。否则JJ1回到步骤a。输出最终的重建信号 ˆΨ。这个过程的直观理解是第一轮DESP从带噪数据中提取出能量最强或最显著的狄拉克特征模式 ˆψ_1。将其从观测数据中减去后残差中可能包含能量次强的特征模式以及噪声。对残差再次应用DESP可以提取出下一个主要模式 ˆψ_2。如此反复逐步剥离出信号中的多个主要成分。4.2 停止准则的重要性防止过拟合噪声迭代算法最大的风险是“过拟合”——将噪声也当作信号成分提取出来。IDESP巧妙地利用系数变异噪声与信号的能量比作为停止准则。c_V^{true}这是观测数据中真实的噪声水平通常需要事先估计例如通过信号平稳区域的统计或先验知识。c_V(J)这是经过J轮迭代后当前重建信号 ˆΨ 与原始观测 ˜ψ 之间的差异视为残余“噪声”与重建信号能量的比值。如果重建是完美的c_V(J) 应等于 c_V^{true}。因此迭代的目标是让 c_V(J) 逼近 c_V^{true}。当 J 过小时c_V(J) c_V^{true}说明信号成分提取不足当 J 过大时c_V(J) c_V^{true}说明已经开始拟合噪声。IDESP选择在 |c_V(J) - c_V^{true}| 达到最小时停止此时重建信号包含了绝大部分真实信号成分而引入了最少的噪声成分。注意事项准确估计 c_V^{true} 是关键。如果估计值偏差太大会导致迭代过早停止或过度进行。在实际应用中可以采用多种方法交叉验证1) 如果数据有重复测量或已知的干净片段可直接计算2) 利用高频部分估计噪声功率3) 将IDESP本身用于参数搜索观察重建信号随迭代次数变化的平台区。4.3 性能验证与真实数据应用论文在海洋浮标轨迹数据上验证了IDESP。该数据将太平洋马达加斯加岛附近的海域划分为网格节点浮标在不同网格间的运动形成边上的流量信号。从边信号出发他们构造了一个合成的节点-边联合拓扑旋量作为真实信号并添加噪声。实验表明误差下降随着迭代次数 J 增加真实误差 ∆(J) ‖ ˆΨ - ψ ‖ 首先快速下降在最优迭代次数 J_opt 附近达到最低点之后由于开始拟合噪声误差缓慢上升。这验证了停止准则的有效性。超越IDSP对比同样迭代框架但使用DSP即固定m0的IDSP算法IDESP无论是用L_m还是S_m优化m能够以更少的迭代次数达到相同或更低的误差水平。这说明学习到的质量参数 m 有效地捕捉了信号的内在结构提高了每次迭代提取特征的效率。可视化改善从噪声严重的原始观测信号到第一次迭代J1提取出主要环流模式再到第二次迭代J2进一步细化结构IDESP逐步重建出了与真实流场高度相似的模式。5. 算法实现关键与实战指南5.1 数据结构与算子构建实现DESP/IDESP的一步是构建拓扑狄拉克算子 D 和 γ 矩阵。import numpy as np import scipy.sparse as sp def build_topological_dirac_operator(B1): 构建1维单纯复形图的拓扑狄拉克算子。 参数: B1: scipy.sparse.csr_matrix, 形状为 (n_nodes, n_edges)。 节点-边关联矩阵。B1[i, j] 1 如果边j以节点i为终点-1如果为起点有向图 对于无向图通常取绝对值或特定方向。 返回: D: scipy.sparse.csr_matrix, 形状为 (n_nodesn_edges, n_nodesn_edges)。 狄拉克算子 [[0, B1^T], [B1, 0]]。 gamma: scipy.sparse.csr_matrix, 相同形状。 γ 矩阵通常为 diag(I_{nodes}, -I_{edges})。 n_nodes, n_edges B1.shape # 构建狄拉克算子 D zero_nn sp.csr_matrix((n_nodes, n_nodes)) zero_ne sp.csr_matrix((n_nodes, n_edges)) zero_en sp.csr_matrix((n_edges, n_nodes)) zero_ee sp.csr_matrix((n_edges, n_edges)) # D [[0, B1^T], [B1, 0]] top_block sp.hstack([zero_nn, B1.T], formatcsr) bottom_block sp.hstack([B1, zero_ee], formatcsr) D sp.vstack([top_block, bottom_block], formatcsr) # 构建 gamma 矩阵: diag(I, -I) gamma_n sp.identity(n_nodes, formatcsr) gamma_e -sp.identity(n_edges, formatcsr) gamma sp.block_diag([gamma_n, gamma_e], formatcsr) return D, gamma5.2 DESP核心求解与参数扫描对于固定的 τ 和 m求解 ˆψ_m 需要解一个线性系统。由于矩阵 (I τ (Dmγ)^T (Dmγ)) 是大规模稀疏正定矩阵应使用迭代法如共轭梯度法CG而非直接求逆。from scipy.sparse.linalg import LinearOperator, cg def desp_single_step(psi_tilde, D, gamma, m, tau, tol1e-8): 执行单次DESP重建固定m。 参数: psi_tilde: np.ndarray, 观测拓扑旋量 (n_nodes n_edges)。 D, gamma: 狄拉克算子和gamma矩阵。 m: float, 质量参数。 tau: float, 正则化参数。 tol: float, CG求解器的容忍误差。 返回: psi_hat: np.ndarray, 重建的拓扑旋量。 n_total D.shape[0] # 构造算子 A I tau * (D m*gamma)^T (D m*gamma) H D m * gamma # 哈密顿量 H D mγ A_operator LinearOperator((n_total, n_total), matveclambda x: x tau * H.T (H x)) # 使用共轭梯度法求解 A * psi_hat psi_tilde psi_hat, info cg(A_operator, psi_tilde, toltol) if info ! 0: print(fCG solver did not converge, info: {info}) return psi_hat def optimize_m_by_loss(psi_tilde, D, gamma, tau, m_range): 通过扫描m范围并最小化损失函数L_m来优化m。 参数: m_range: np.ndarray, 待搜索的m值数组。 返回: best_m: float, 最优m值。 best_psi: np.ndarray, 对应最优m的重建信号。 losses: list, 各m对应的损失值。 best_loss float(inf) best_m None best_psi None losses [] for m in m_range: psi_hat desp_single_step(psi_tilde, D, gamma, m, tau) # 计算损失 L_m ||psi_tilde - psi_hat||^2 tau * ||(Dm*gamma)psi_hat||^2 residual psi_tilde - psi_hat data_fidelity np.linalg.norm(residual)**2 H_psi (D m*gamma) psi_hat regularization tau * np.linalg.norm(H_psi)**2 loss data_fidelity regularization losses.append(loss) if loss best_loss: best_loss loss best_m m best_psi psi_hat.copy() return best_m, best_psi, losses5.3 IDESP迭代实现IDESP的实现需要仔细管理残差和停止准则。def idesp(psi_tilde, D, gamma, tau, c_v_true, max_iter20, m_optimizerloss): IDESP算法主函数。 参数: psi_tilde: 观测信号。 c_v_true: float, 真实的系数变异噪声水平估计。 max_iter: 最大迭代次数。 m_optimizer: loss 或 rdre优化m的策略。 返回: Psi_hat: 最终重建信号。 m_list: 每次迭代学习到的m值列表。 c_v_list: 每次迭代的c_V(J)值列表。 errors: 每次迭代的真实误差如果有真实信号psi_gt。 n_total len(psi_tilde) residual psi_tilde.copy() Psi_hat np.zeros(n_total) m_list [] c_v_list [] prev_diff float(inf) for j in range(1, max_iter1): # 1. 对当前残差应用DESP优化m if m_optimizer loss: m_opt, psi_hat_j, _ optimize_m_by_loss(residual, D, gamma, tau, m_rangenp.linspace(-2, 2, 41)) else: # rdre m_opt, psi_hat_j, _ optimize_m_by_rdre(residual, D, gamma, tau, m_rangenp.linspace(-2, 2, 41)) m_list.append(m_opt) # 2. 更新总重建信号和残差 Psi_hat psi_hat_j residual psi_tilde - Psi_hat # 3. 计算当前系数变异 c_V(J) norm_Psi np.linalg.norm(Psi_hat) if norm_Psi 1e-12: c_v 1.0 else: c_v np.linalg.norm(residual) / norm_Psi c_v_list.append(c_v) # 4. 检查停止准则 (简单版本当c_V开始上升或变化很小时停止) diff abs(c_v - c_v_true) if diff prev_diff: # 本次迭代后离目标更远了 print(fStopping at iteration {j}, c_V starts to diverge.) # 回退到上一次迭代的结果 Psi_hat - psi_hat_j m_list.pop() c_v_list.pop() break prev_diff diff # 可选如果c_V已经非常接近真实值也可停止 if diff 0.01 * c_v_true: print(fStopping at iteration {j}, c_V is close enough to target.) break return Psi_hat, m_list, c_v_list5.4 参数选择与调优经验正则化参数 τ控制平滑度与数据保真度的权衡。τ 过大重建信号过于平滑可能丢失细节τ 过小去噪效果弱。建议从 τ 10 开始尝试根据信噪比调整。高噪声下可适当增大 τ。质量参数 m 的搜索范围需要根据先验知识或数据尺度设定。对于标准化后的信号m 通常在 [-3, 3] 区间内。可以先进行粗搜索如步长0.5在最优值附近再进行细搜索。真实噪声水平 c_V^{true} 的估计这是IDESP的关键。如果无法直接获得可以使用高通滤波或小波变换估计噪声方差。在数据中寻找理论上应为“平坦”或“平滑”的区域计算该区域的方差作为噪声估计。将 c_V^{true} 作为一个超参数在验证集上调整选择使最终应用指标如分类准确率、预测误差最优的值。最大迭代次数 max_iter不宜设置过大通常5-10次迭代足以提取主要模式。设置一个上限防止无限循环并通过监控 c_V(J) 曲线来观察收敛情况。6. 应用场景与未来展望DESP/IDESP算法为处理具有内在拓扑结构的高维数据提供了新工具。计算神经科学脑功能网络中的信号不仅存在于脑区节点也存在于连接边上。DESP可以联合处理节点上的BOLD信号时间序列和边上的功能连接强度动态可能更精准地识别与认知任务相关的脑网络模式。交通流与物流分析在交通网络中节点代表路口或城市边代表道路。节点信号可以是车流量边信号可以是行程时间或拥堵指数。IDESP可以分解出影响全局交通的主要传播模式如早晚高峰潮汐流。物理场重构在传感器网络中节点测量标量场如温度边可以测量梯度或通量。DESP利用狄拉克方程耦合标量和梯度信息能在部分传感器失效或噪声较大时实现更鲁棒的物理场重构。分子动力学与生物信息学蛋白质残基节点和它们之间的相互作用边构成网络。节点信号可以是氨基酸特性边信号可以是结合能。拓扑狄拉克分析可能帮助识别稳定的构象或功能位点。未来方向与挑战扩展到更高阶当前框架主要针对节点和边1-复形。如何将其推广到包含三角形、四面体等更高阶单形的复形处理更丰富的高阶相互作用是一个前沿问题。与图神经网络结合狄拉克算子可以作为GNN中的消息传递算子其固有的处理节点和边特征的能力可能缓解传统GNN的过度平滑和过度挤压问题。DESP的滤波思想可以融入图卷积层的设计。处理时变拓扑许多现实网络的拓扑结构是动态变化的。发展适用于时变单纯复形的在线DESP算法是一个实际需求。理论收敛性分析目前IDESP的停止准则虽直观有效但缺乏严格的理论收敛性保证。在什么条件下算法能稳定地恢复信号的主要成分仍需深入研究。在我自己的尝试中将IDESP应用于社交网络用户-群组联合分析时发现初始的关联矩阵 B1 的构造方式对结果影响巨大。例如在定义“用户-边”关系时是采用简单的隶属关系还是考虑交互频率的加权会导致学到的“质量参数” m 具有完全不同的解释。一个负的 m 值可能暗示着节点和边信号之间存在某种抑制性关联这为数据解读提供了新的视角。这提醒我们算法背后的拓扑模型必须紧密结合领域知识来构建否则“物理启发”就可能沦为“数学游戏”。
