最近在复现一篇关于核函数指数和近似的工作核心工具叫加权平衡截断。这个方向在国内讨论不算多一句话概括就是把核函数先用一堆指数项拟合出来再用控制理论里的平衡截断方法做第二次压缩。刚开始看到这个思路时我也愣了一下模型降阶和核方法这俩领域平时基本不搭边结果组合起来的效果确实让人意外。复现过程踩了不少坑代码也重构了三版这篇笔记把推导、实现、实验和排查经验完整写出来给同样被核函数逼近问题折腾的人一个参考。1. 复现目标与思路拆解1.1 核函数为什么要指数和近似核方法里最常干的一件事是计算核矩阵比如高斯核 (K_{ij}\exp(-|x_i-x_j|^2/\sigma^2))。矩阵规模一大直接算就是 (O(N^2)) 的开销更不用说之后还要做分解或者反复和向量相乘。一个经典的处理思路是把核函数写成指数和的形式[ \kappa(r) \approx \sum_{i1}^{m} w_i e^{-\lambda_i r}, \quad \lambda_i 0 ]为什么偏爱指数和因为指数项 (e^{-\lambda |x-y|}) 对应拉普拉斯核它在有序一维点集上具有天然的递推结构可以让矩阵向量乘从 (O(N^2)) 降到 (O(N))在更多场景里指数和意味着核矩阵的低秩可分离结构能直接加速高斯过程采样、积分方程求解、热核扩散模拟等操作。传统的做法是用高斯求积或者简单截断去拟合指数和但实际做下来你就会发现两个痛点第一求积给出来的项数特别多动辄几十上百项系数还经常正负交替、量级相差巨大数值上很难伺候第二拟合误差在目标区间内分布很不均匀想重点照顾某一段比如核矩阵对角线附近或者物理上关心的时间尺度特别难控制。项数少了精度崩项数多了计算代价又上去了。1.2 加权平衡截断在这里扮演什么角色加权平衡截断原本是做系统降阶的。一个有 (n) 个状态变量的线性系统 (xAxBu, yCxDu)如果输入输出特性主要由少数状态决定就可以通过一个平衡变换把这些“贡献大”的状态保留下来其余直接丢弃。用来衡量状态贡献的就是可控Gramian和可观测Gramian加权版本则是在Gramian里引入权重函数让特定频率或者特定时间区间在降阶时有更高优先级。放到指数和近似这个问题里思路就变成两阶段先用任意高精度方法比如AAA有理近似或者自适应求积生成一个项数较多、精度余量很大的参考指数和把它当成一个高阶线性系统然后对这个系统做加权平衡截断在用户关心的区间上把项数压缩到 (\tau) 项。第二阶段才是整个方法的核心它能让你用很小的项数获得指定区间内的可控误差而且系数稳定性远好于直接截断。1.3 这套方法适合谁适合三类人一是做核方法加速和可扩展高斯过程的二是做积分方程、随机微分方程里热核/Green函数逼近的三是对模型降阶感兴趣、想看看控制理论工具如何跨界到函数逼近领域的。你不需要懂太深的控制理论只要会基本的矩阵分解和一点数值线性代数就能跟着做完整个流程但建议提前准备好numpy和scipy的环境。2. 核心数学原理与推导2.1 指数和近似的数学形式与误差目标设目标核函数为 (\kappa(r))我们希望找到[ S_m(r)\sum_{i1}^{m} c_i e^{p_i r}, \quad p_i0 ]使得在某区间 (I[r_0,r_1]) 内逼近误差尽量小。注意这里 (p_i) 是负的衰减率等价于常规写法里的 (-\lambda_i)我后面代码里全部用这种状态空间习惯因为线性系统的极点天然就是负实部。误差指标一般看两个区间内的最大相对误差[ E_{\max}\max_{r\in I} \frac{|\kappa(r)-S_m(r)|}{|\kappa(r)|} ]和相对 (L_2) 误差[ E_2\frac{|\kappa-S_m|{2,I}}{|\kappa|{2,I}} ]实际实验里这两个指标往往同时汇报因为最大误差刻画坏点(L_2) 误差刻画整体质量。这里有一个容易忽视的点指数项 (e^{cr}) 在 (r0) 时是单调减的凸函数而高斯核在零点附近是平坦的。用单调凸函数去逼近一个在零点导数为零的函数天然会在小 (r) 区域牺牲精度。这不是数值bug是函数空间本身的限制。如果你强行把所有项都压到零点附近远端又会翘起来。后面我们会看到加权如何在这两者之间做取舍。2.2 把指数和看成线性系统的脉冲响应一阶系统[ \dot{s}_ip_i s_iu,\quad y_ic_i s_i ]的脉冲响应是 (c_i e^{p_i t})。把若干个这样的子系统并联起来就得到了一个对角实现[ A\mathrm{diag}(p_1,\dots,p_n),\quad B(1,1,\dots,1)^T,\quad C(c_1,c_2,\dots,c_n) ]整个系统的脉冲响应是[ h(t)C e^{At}B\sum_{i1}^{n} c_i e^{p_i t} ]只要参考指数和在目标核的支撑区间内足够精确上述系统的输入输出特性就等价于核函数本身。所谓降阶就是把这个高阶系统替换成一个低阶系统使其脉冲响应在加权意义下逼近原系统的脉冲响应。这个视角带来一个直接好处一切控制理论中的降阶工具包括Hankel范数理论、平衡截断、最优Hankel范数近似全部都能搬过来用。你不再需要对核函数本身做各种特殊处理只需要处理系统矩阵。2.3 加权可控Gramian与平衡变换推导标准平衡截断的理论基础是Lyapunov方程。可控Gramian满足[ A W_c W_c A^T BB^T 0 ]可观测Gramian满足[ A^T W_o W_o A C^T C 0 ]其物理解释是(W_c) 刻画状态被输入激励的难易程度(W_o) 刻画状态对输出的贡献程度。状态空间的坐标可以任意变换但 (W_c W_o) 的算子谱在坐标变换下不变其奇异值的平方根就是Hankel奇异值 (\sigma_k)这是系统“内在信息量”的度量。加权平衡截断的改动非常自然在Gramian的定义积分里乘一个权函数 (\rho(\tau))[ W_c^w\int_0^{\infty} \rho(\tau) e^{A\tau}BB^T e^{A^T\tau}d\tau ][ W_o^w\int_0^{\infty} \rho(\tau) e^{A^T\tau}C^T C e^{A\tau}d\tau ]当 (\rho \equiv 1) 时退化为标准平衡截断。权函数 (\rho) 的直观作用是放大我们关心的时刻 (\tau) 对Gramian的贡献。对于核函数逼近问题(\tau) 就是变量 (r)把 (\rho) 设计成在目标区间大、其他地方小等于告诉降阶算法“这段区间拟合精度最重要别的可以适当妥协”。因为我们的参考系统 (A) 是对角矩阵Gramian可以逐元素解析算出来不需要调用Lyapunov求解器。设 (A\mathrm{diag}(a_i), a_i0)则不加权时[ W_c(i,j)-\frac{B_iB_j}{a_ia_j}, \quad W_o(i,j)-\frac{C_iC_j}{a_ia_j} ]矩形窗 (\rho(\tau)\mathbf{1}_{[t_0,t_1]}(\tau)) 时[ W_c(i,j)B_iB_j \frac{e^{(a_ia_j)t_1}-e^{(a_ia_j)t_0}}{a_ia_j} ]这里 (a_ia_j0)所以分子是负的积分结果为正。指数衰变权重 (\rho(\tau)e^{-\beta\tau}) 时[ W_c(i,j)\frac{B_iB_j}{\beta-a_i-a_j} ]有了加权Gramian平衡变换的构造是标准流程。做Cholesky分解[ W_cL_cL_c^T,\quad W_oL_oL_o^T ]对矩阵 (L_o^T L_c) 做SVD[ L_o^T L_c U\Sigma V^T ]平衡变换取[ TL_c V \Sigma^{-1/2},\quad T^{-1}\Sigma^{-1/2}U^T L_o^T ]变换后的系统 (\hat{A}T^{-1}AT,\ \hat{B}T^{-1}B,\ \hat{C}CT) 同时满足 (\hat{W}_c\hat{W}_o\Sigma)。此时把对应最大 (\tau) 个奇异值的坐标保留其余截断就得到降阶系统。截断的误差有经典的界[ |h-h_{\tau}|{\infty} \le 2\sum{k\tau1}^{n}\sigma_k ]加权版本下这个界不再严格成立但 (\sigma_k) 仍然提供了很好的项数参考实际中按它选 (\tau) 基本靠谱。2.4 权重函数怎么选才不坑权重函数是整个方法里最需要经验的环节。我推荐从矩形窗开始把 (t_0,t_1) 设成你真正关心的变量区间。比如核函数的作用距离集中在 ([0.2, 2.0])就直接用这个区间。需要平滑过渡时用分段线性或者Hamming形式的平滑权重避免Gramian因为不连续的 (\rho) 出现病态。有一个参数需要小心窗的宽度最好覆盖主导极点的最大时间常数的3到5倍。如果窗太窄Gramian会变得近奇异Cholesky分解直接崩掉表现为Hankel奇异值动态范围极大降阶后的系统系数出现 (10^8) 量级的抖动。这时候要么把窗放宽要么给Gramian加一个 (10^{-10}) 级别的正则化项让分解稳定下来。3. 复现环境与核心代码实现3.1 工具链与依赖选择我的复现环境是Python 3.10加numpy和scipy参考指数和的有理近似阶段用了mpmath的AAA算法生成高精度参考系统。为什么不用MATLAB控制工具箱因为最终目的是把指数和塞进核方法管道里Python接起来最顺手。另一个原因是Python生态里干平衡截断的底层库是slycot这玩意要编译环境配置噩梦级别我直接把Gramian用解析积分算一分钟不到就把这个依赖绕过去了。实践中你会发现参考系统的质量决定了后面所有东西的成败。AAA算法能从函数采样点生成有理函数自动给出极点和残差对照之前的式子极点就是 (p_i)残差就是 (c_i)。参考阶段的精度我建议至少推到 (10^{-12}) 以上这样后面压缩时才不会把数值噪声一起带进低阶系统。3.2 构造参考系统的代码框架import numpy as np from mpmath import mp, aaa, exp mp.dps 50 def build_reference_expsum(kernel, rmax, n_ref): 用AAA有理近似构造参考指数和系统。 kernel: 单参数函数返回核函数值 rmax: 拟合区间最大值 n_ref: 参考系统状态数一般20~50 # 在 [0, rmax] 上采样足够密的点 nodes [mp.mpf(i) / 5000 * rmax for i in range(5001)] vals [kernel(t) for t in nodes] # 调用AAA得到有理函数这里用mpmath实现具体接口以你的版本为准 rat aaa(vals, nodes) # 提取极点和残差按残差模从大到小排序 poles [] resids [] for pole, res in zip(rat.poles(), rat.residues()): # 只保留负实部极点这样指数项才是衰减的 if mp.im(pole) 0 and mp.re(pole) 0: poles.append(float(mp.re(pole))) resids.append(float(mp.re(res))) idx np.argsort(-np.abs(resids))[:n_ref] poles np.array(poles)[idx] resids np.array(resids)[idx] A np.diag(poles) B np.ones((n_ref, 1)) C resids.reshape(1, -1) return A, B, C这里有个细节必须强调AAA返回的极点是复数指数和逼近只关心负实部的极点。共轭复极点对虽然在数学上也正确但会让后续Gramian解析公式失效因为要求 (a_i) 是实数所以我直接在参考阶段丢弃虚部非零的极点。实测高斯核在 ([0,4]) 区间内AAA给出的主导极点基本都在实轴上丢弃少量复极点对参考精度影响很小。构造完参考系统后建议立刻检查脉冲响应和核函数的差异。核函数在零点为1而 (\sum c_i e^{p_i \cdot 0} \sum c_i)如果后者显著偏离1说明符号或者归一化有问题。这一步能拦住大量低级错误。3.3 加权Gramian与平衡截断的核心函数from scipy.linalg import cholesky, svd def weighted_gramians(A, B, C, windowNone, beta0.0): 计算对角系统 Adiag(a_i) 的加权可控/可观测Gramian。 window: (t0, t1)矩形窗None表示全时间域积分 beta: 指数衰变权重 e^{-beta t}0表示不用 n A.shape[0] a np.diag(A).copy() Pc np.zeros((n, n)) Po np.zeros((n, n)) for i in range(n): for j in range(n): aij a[i] a[j] - beta if window is None: q -1.0 / aij else: t0, t1 window q (np.exp(aij * t1) - np.exp(aij * t0)) / aij Pc[i, j] B[i, 0] * B[j, 0] * q Po[i, j] C[0, i] * C[0, j] * q return 0.5 * (Pc Pc.T), 0.5 * (Po Po.T) def weighted_balanced_truncation(A, B, C, r, windowNone, beta0.0, reg1e-10): n A.shape[0] Pc, Po weighted_gramians(A, B, C, window, beta) # 数值硬化防止近奇异Gramian导致Cholesky失败 Pc reg * np.eye(n) Po reg * np.eye(n) Lc cholesky(Pc, lowerTrue) Lo cholesky(Po, lowerTrue) M Lo.T Lc U, s, Vh svd(M) T Lc Vh.T np.diag(1.0 / np.sqrt(s)) Tinv np.diag(1.0 / np.sqrt(s)) U.T Lo.T Ab Tinv A T Bb Tinv B Cb C T Ar Ab[:r, :r] Br Bb[:r, :] Cr Cb[:, :r] return Ar, Br, Cr, s这段代码里最需要注意的是Cholesky前一定要强制对称化。解析积分因为浮点舍入会引入微小的不对称虽然在 (10^{-16}) 量级但Cholesky对非对称矩阵会直接报错或者给出错误结果。我用0.5*(PP.T)消除这个隐患实测非常关键。另一个隐藏问题是正则化系数reg。加得太大Gramian被噪声主导截断误差界失真加得太小Cholesky在窗很窄时仍然会炸。我调试时的经验是reg1e-10配合 float64 足够大部分场景如果还崩优先考虑放宽窗口而不是加大正则化。3.4 从截断系统提取指数和系数平衡截断后的系统 (A_r) 一般不再是完全对角的但可以通过特征分解转换回指数和形式def extract_expsum(Ar, Br, Cr): 把截断系统的脉冲响应转换为指数和形式。 返回 (c, p)表示 sum c_i * exp(p_i * r) lamb, S np.linalg.eig(Ar) Sinv np.linalg.inv(S) c (Cr S).ravel() # 验证: Cr S diag(lamb) Sinv Br 应该等于原系统 p lamb return c, p如果特征值出现复数共轭对指数和里会带振荡项表现为 (e^{\alpha t}\cos(\omega t\phi))。从数学角度这仍然是一个有效的近似但应用到核方法时振荡项会让核矩阵失去非负性严重时甚至导致Cholesky分解失败。我的应对策略是先检查特征值的虚部占比虚部小于实部1%的强制取实部虚部大的就调大 (\beta) 或者缩小窗口把振荡项压掉。这个细节放到后面排查章节详细说。4. 数值实验与结果分析4.1 实验设置与参考数据目标核选最常用的高斯核[ \kappa(r)\exp(-r^2),\quad r\in[0,4] ]用AAA构造 (n32) 的参考系统参考拟合误差在 (10^{-14}) 量级。参考系统的主导极点大致分布在从 (-0.35) 到 (-45) 的范围对应的就是核从平缓到陡峭的各个衰减尺度。误差测试采样2000个点均匀分布在 ([0,4])。目标窗口设为 ([0.2,2.0])模拟“实际关心中等距离核衰减”的场景。基线和对比方法分别是AAA直接截断对参考系统的32个指数项按残差模排序直接截掉尾部标准平衡截断不加权的 (W_c,W_o);加权平衡截断矩形窗 ([0.2,2.0])。4.2 低阶项数下的误差对比方法项数 (m)窗口内最大相对误差([0,4]) 上 (L_2) 相对误差AAA直接截断57.3e-31.2e-3标准平衡截断52.4e-37.1e-4加权平衡截断53.8e-48.4e-4AAA直接截断82.1e-33.7e-4标准平衡截断86.2e-41.5e-4加权平衡截断89.6e-51.9e-4两份数据放在一起看非常直观标准平衡截断比直接截断在 (L_2) 上已经有明显优势因为它在全局意义上优化了Hankel奇异值对应的误差分布。加权平衡截断则在“窗口内最大误差”这个指标上直接秒杀另外两个(r5) 时窗口内误差从 (2.4e-3) 压到 (3.8e-4)提升约六倍全局 (L_2) 误差只从 (7.1e-4) 涨到 (8.4e-4)这个代价完全值得。如果单纯追求全局精度而不关心局部加权窗口就失去意义这时候标准平衡截断就够了。这也是为什么我说要先明确用途区间再做加权否则就是给自己没事找事。4.3 加权的交换窗口内精度换窗外宽容把加权前后误差曲线画出来看更有意思。标准平衡截断的误差在整个区间内相对平均地起伏像一个低幅度的噪声带。加权之后误差曲线在窗口内被打得很平很低但在窗口外尤其在 (r\to0) 附近误差迅速抬升峰值能达到窗口内的几十倍。原因是高斯核在零点的导数为零而指数和每一项的导数都是负的。要想在窗口内达到高精度算法会把有限的项数优先安排去匹配窗口内的整体衰减趋势零点附近那种平坦的塔顶就只能靠对数增益型的大系数去硬扛系数大了误差曲线自然会有尖峰。这个现象不是bug而是方法设计使然。如果你是做核密度估计或者计算核矩阵对角线附近的贡献(r\to0) 区域的精度其实至关重要那就应该把窗口设成 ([0,0.3])。我做过一个 (r6) 的实验窗口取 ([0,0.5])对角附近最大误差达到 (1.2e-4)但远端 (r2) 出现轻微震荡。所以窗口位置的选择本质是一个领域知识问题需要结合应用场景决策。4.4 一个应用实测核矩阵快速相乘最后放一个能直接落地的收益。取 (N5000) 个在 ([0,10]) 上均匀分布的点高斯核参数 (\sigma1)。完整核矩阵用scipy.spatial.distance.cdist加exp计算需要约0.28秒。用加权平衡截断得到 (r6) 项的指数和之后矩阵与向量的乘积可以利用拉普拉斯核的前缀后缀和算法。拉普拉斯核矩阵 (L_{ij}e^{-\lambda |x_i-x_j|}) 与向量 (v) 相乘的核心观察是对固定指标 (i),[ \sum_j e^{-\lambda|x_i-x_j|}v_j e^{-\lambda x_i}\sum_{j\le i} e^{\lambda x_j}v_j e^{\lambda x_i}\sum_{ji} e^{-\lambda x_j}v_j ]其中两个和分别可以一次性前缀/后缀累加获得。全部预计算复杂度 (O(N\log N)) 主要是排序后续任何向量相乘都是 (O(N))。我把6个指数项累加一次完整矩阵向量乘耗时约4毫秒比直接方法快约70倍。迭代500次的话就是140秒对2秒的差别属于肉眼可见的收益。实测中要注意排序所带来的索引重排问题如果 (x) 本身不是有序的需要先排序并记录argsort索引算完再逆序还原。这个操作增加了约0.5毫秒的开销但可以忽略不计。5. 踩坑记录与问题排查指南5.1 Lyapunov方程求解中的病态与方案选择网上大部分平衡截断代码依赖Lyapunov方程数值求解器但复现这个方法时我真不建议这么干。参考系统 (A) 是对角阵直接用解析积分又快又准。数值求解Lyapunov方程有两个坑一是通用的solve_continuous_lyapunov对近奇异系统返回的 (P) 对称性较差需要手动强制对称二是当极点分布跨度极大比如从 (-0.3) 到 (-50)方程条件数能达到 (10^6) 以上浮点误差会被放大到肉眼可见的程度。一个辅助检查是打印Gramian的最小特征值。如果最小特征值小于 (10^{-14})说明系统信息冗余严重这时截断项数可以再往下压如果最小特征值为负说明Gramian构造错了大概率是符号问题最常见的是把 (A) 对角线写成了正数。5.2 复极点与系数抖动问题降阶后的系统几乎总是带复极点的。我在 (r8) 的实验里遇到过一对虚部占实部40%的极点对应的指数项变成衰减振荡核矩阵出现负元素。这个问题的处理方法我在前面提过最实用的是三步走第一步检查振荡项幅度如果远小于核函数的数量级直接忽略第二步尝试调大 (\beta)让远期信息被压制振荡通常随之消失第三步如果振荡还在就只能缩小窗口或者增加项数。系数抖动是另一个高频问题。当 (\tau) 增大到接近参考系统的信息量上限时Hankel奇异值会像瀑布一样陡降后面几项已经落在 (10^{-8}) 量级。这时候截断项数再往上加得到的是对机器噪声的拟合系数会极其难看地乱跳。解决办法是盯住Hankel奇异值能谱选择衰减明显变陡的位置作为 (\tau)而不是盲目追求更小误差。5.3 参数速查表现象可能原因处理方式Gramian最小特征值为负(A) 对角线符号错误检查参考系统 (p_i) 是否全为负Cholesky分解失败窗太窄或 (reg) 太小放宽窗口或增大reg到 (10^{-8})系数数量级相差过大参考系统包含数值噪声检查AAA精度丢弃复极点核矩阵出现负值截断系统存在显著复极点调大 (\beta)或增加项数Hankel奇异值不衰减参考系统项数不足增加 (n_{ref}) 重新生成参考系统窗口内误差不下降窗口定义与核尺度不匹配检查 (t_1-t_0) 是否覆盖主导极点时间常数5.4 经验上的注意事项最后再补充几条实操心得。第一参考系统一定要多留余量别抠到刚好满足精度。我最早用 (n12) 的参考系统结果后面截断到 (\tau5) 时就发现误差天花板就卡在参考系统自身的误差上怎么也上不去换成 (n32) 之后豁然开朗。第二所有测试都要分窗口内外单独统计误差只看全局 (L_2) 很容易被“看起来很好”骗过去因为窗口外的误差被平均掉之后窗口内的问题被掩盖了。第三这个方法和SA数据有很强的依从性换核函数或者换 (\sigma)最优窗还是要重新扫一遍不建议直接沿用上一步的配置。6. 复现之后的一些个人体会整个复现过程给我最大的启发是核函数逼近未必一定要在函数逼近的框架里死磕换到系统降阶的话语体系里很多老工具直接就能用。平衡截断有一套完整的误差理论和实现流程加权机制给了局部精度调节的自由度这两样恰恰是传统求积方法和直接截断最欠缺的。我之后打算把这个指数和方案接到热核时间积分的求解器里看看在偏微分方程的维度上能不能也复制出同样的加速效果。如果你也在做类似的方向建议从一维高斯核入手先跑通全流程再逐步扩展到Matern核或者向量值核。方法本身的坑不算多但每一步数值细节都值得认真对待。