指数和近似与加权平衡截断：核矩阵高效压缩的完整方案

我去年在一个大规模高斯过程回归项目里被核矩阵卡得欲仙欲死后来在翻老论文时发现了一个组合思路——指数和近似加上加权平衡截断能把核函数逼近和模型降阶两个领域串起来。当时的第一反应是这俩东西怎么凑到一起的写这篇文是想把它彻底搞清楚从“为什么要做指数和近似”到“加权平衡截断到底在压什么”再到一步步复现的代码和踩坑记录给同样被核矩阵尺寸折磨的人一条可走的路。1. 内容整体设计与思路拆解先说动机。高斯过程、核岭回归、径向基函数插值这些方法的核心在于构造一个核矩阵K_ij κ(x_i, x_j)然后解一个K α y的线性系统。采样点一多——比如超过一万个——核矩阵的存储就是 O(n^2)分解是 O(n^3)直接卡死。加权平衡截断恰恰是对这块进行结构化压缩与其直接对 K 做低秩分解不如先把它放到“系统”的框架里然后用平衡截断的办法把它降到一个更小的状态空间模型。这时对核函数做指数和近似的意义才真正浮现。再说框架。指数和近似解决的是“把核函数拆成可分离的形式”。一条关于核函数的经典性质是很多常见核函数比如高斯核exp(-x^2)、Matern核、拉普拉斯核都可以表示为指数函数的积分或者指数函数的加权和。这意味着能够写成κ(x, y) ≈ Σ_i w_i exp(a_i |x - y|)或者Σ_i w_i exp(-b_i (x - y)^2)这样的形式。一旦做到这个核矩阵就能写成低秩因子的乘积K ≈ V Λ V^T其中 V 只涉及逐点求值而 Λ 是低维对角阵。这一步把“存储和输入端”彻底解放了。加权平衡截断则解决“模型降阶”问题。系统理论里有一个经典的平衡截断Balanced Truncation方法给定一个高维线性时不变系统先计算可控性 Gramian 和可观性 Gramian再做一次平衡变换让这两个矩阵同时变成对角阵并且对角线上的值——汉克尔奇异值——按从大到小排好然后只保留前 r 个最大的分量丢掉其余部分。加权版本是在 Gramian 的定义里引入权重矩阵从而让降阶过程偏向某些频段或某些输出方向。一旦把核矩阵对应的“线性系统”建立起来用加权平衡截断对它做降阶再配合指数和近似的低秩结构整个项目就从“两个毫不相干的方法”变成了“一条完整的压缩流水线”。先看整体设计思路再逐层拆解关键理论然后给一个可复现的代码流程最后是调试记录。1.1 这条流水线的核心模块整个复现项目被我拆成四条主线。第一把核函数表达为指数和的近似。目标是对给定的核函数 κ 和采样域求出权重 w_i 和指数参数 α_i使得近似误差可控。工程上常用的是通过有理逼近如 AAA 算法把核的 Laplace 变换形式找出来或是直接用数值求积方法从积分表示里离散化出一组合适的指数项。第二把核矩阵转化成状态空间模型。这里的关键一步是设计一个“虚拟输入输出系统”让它的传递函数与核函数对应。并不是说要从零开始严格建模而是构造一组矩阵 A、B、C使得整体系统的输入输出特性与核矩阵的谱特性一致。这一步是整条链中最绕的部分也是加权平衡截断与核方法之间那座桥。第三用加权平衡截断做降阶。输入是这个状态空间模型经过平衡变换和截断输出一个低维系统。它对应的核矩阵被替换为低秩近似矩阵求逆成本从 O(n^3) 掉到 O(r^2 n)。第四把近似后的核函数放回原来应用里。无论是高斯过程回归还是径向基插值都用这个低秩矩阵来完成训练和预测并且用理论误差界或实验误差来验证。选择加权平衡截断而不是直接把核矩阵做奇异值分解原因不只是低秩。平衡截断自带一个误差界而且在频域上的行为有保证——它压在“系统能量”的意义上是全局最优的近似这是纯 SVD 不具备的性质。加权版本还可以让误差集中于频段这在信号处理类的核学习场景非常值钱。2. 前置知识指数和近似与核函数的可分离性这项工作的第一块理论地基是“可分离性”。如果核函数能写成κ(x, y) φ(x)ᵀ ψ(y)那么核矩阵天然是低秩的。很多教科书里的核函数并不满足有限维度下的可分离性——高斯核的展开其实是无穷维的。指数和近似做的事情就是把无穷维展开截断成有限维指数项之和用有限项换一个带误差的可分离近似。拉普拉斯核是最好的例子。它的积分表示是exp(-|x - y|) (2/π) ∫₀^∞ [cos(t(x-y)) / (1t²)] dt对积分做离散化——比如用切比雪夫求积或者 Gauss-Laguerre 求积——就能得到一个有限和。类似地高斯核可以用如下积分表示exp(-x²) (1/√π) ∫₋∞^∞ exp(-t²) exp(2i x t) dt离散化这个积分就得到κ(x, y) ≈ Σⱼ wⱼ exp(i tⱼ (x-y))。本质上这是一条“把核写成指数函数的加权叠加”的路径权重和指数参数来自数值求积的节点与系数。我实际复现时最先试的是高斯求积做拉普拉斯核结果精度不够误差大概在 1e-3 量级。后来换成了自适应切比雪夫插值拟合理由逼近配合极点-留数展开误差压到了 1e-10。这个差异在最终实验结果上是致命的粗略的指数近似会把平衡截断的误差曲线上提几个量级让原本好的降阶方案看起来不可用。指数和近似的核心收益在于一旦完成核矩阵可以被分解为K ≈ V D Vᵀ其中 V 是 n×m 的因式矩阵m 是指数项数通常只需 10~20 项D 是 m×m 的对角阵。存储量从 O(n²) 掉到 O(nm)求解线性系统时直接用 Woodbury 恒等式(V D Vᵀ σ²I)⁻¹ σ⁻²I - σ⁻² V (D⁻¹ σ⁻² VᵀV)⁻¹ Vᵀ这个公式就是高斯过程里所谓的“低秩近似求逆”的压缩版。做一个简单的场景模拟n20000 的高斯过程回归传统方法需要 O(2×10⁹) 的存储分解约 O(8×10¹²) 次浮点运算用指数和近似后存储为 O(2×10⁵)求逆主要成本是小矩阵的分解 O(20³) 加上大矩阵乘一张显卡的算力完全能扛下来。但是有一个坑必须提前说指数和近似出来的低秩形式是“对角加低秩”而加权平衡截断处理的是“一般线性系统”的 Gramian这两者不是天然兼容的。打通的关键在于把K ≈ V D Vᵀ当成输出矩阵 C把对角因子当成系统内部结构的权重然后在这个结构下计算可控/可观 Gramian。这正是加权平衡截断里“权重”二字的来源——C 背后那一堆指数项对应不同频段或者不同方向的贡献需要加权才有物理意义。3. 核心方法解析从传统的平衡截断到加权版本平衡截断是基于一个基本观测系统的输入到内部状态的可控性和状态到输出的可观性都可以用 Gramian 来度量。可控性 Gramian P 由 Lyapunov 方程给出A P P Aᵀ B Bᵀ 0可观性 Gramian Q 则由对偶方程给出Aᵀ Q Q A Cᵀ C 0P 的大特征值对应的状态方向容易被输入激发Q 的大特征值对应的状态方向容易被输出观测。如果只保留 P 和 Q 同时大的方向就能得到一个小系统。平衡截断的经典做法是解两个 Lyapunov 方程得到 P 和 Q对 P 做 Cholesky 分解P L Lᵀ构造Lᵀ Q L做奇异值分解Lᵀ Q L U Σ Vᵀ构造变换矩阵T L U Σ^{-1/2}用 T 变换原系统(A, B, C)然后取出前 r 个状态分量。这里的 Σ 对角线上的 σᵢ 就是汉克尔奇异值其衰减速度直接决定了系统可以被压缩到什么程度。如果 σᵢ 从第 r1 个开始趋近于零那么保留前 r 个状态带来的截断误差有一个被广泛使用的上界‖G - Ĝ‖_∞ ≤ 2 Σ_{ir} σᵢ这个界是平衡截断被称为“有保障的模型降阶”的原因——不像很多启发式低秩方法只保证经验上不错它的频域误差是被严格控制的。加权平衡截断在经典版本上做了一处关键改动在 Lyapunov 方程里加入权重矩阵。常见的两种形态频域加权对特定的频率区间强调其重要性Gramian 变成频率加权的积分方向加权对输出/输入通道加一个非对称权重让截断误差在某些方向被放得更小。形式上加权可控 Gramian 的定义变成P_w ∫₀^∞ e^{At} B W_in Bᵀ e^{Aᵀt} dt其中 W_in 是一个半正定矩阵。可观 Gramian 同理。这样做的直接效果是即便汉克尔奇异值整体衰减不快只要目标频段的能量集中在前几个方向上加权之后的截断也可能做到低阶高精度。在我的复现中使用加权平衡截断的场景是核函数在低频段的行为对最终预测影响巨大而高频部分可以容忍较大的近似误差。如果对所有频段一视同仁地截断可能被迫保留很多“无关紧要”的维度来保护高频分量加权后就能把宝贵的低秩资源分配给关键频段。3.1 为什么直接用 SVD 不够看到这里有人会问既然最后都要截断为什么不直接对核矩阵做 SVD 取前 r 个奇异向量这个问题我一开始也想不通直到我用数值实验比对了两者。普通过程是这样把核矩阵 K 求 SVDK ≈ U_r S_r V_rᵀ然后用这个低秩近似替换原文。对遗传算法等简单回归任务结果“看起来还可以”。但一旦进入频域检验问题就暴露了SVD 的低秩近似是最小二乘意义下的最优但它在频域上的误差不是均匀分布的——高频分量经常被砍得很凶导致高阶动态失真。平衡截断则不一样它本质上是在“系统响应”的空间做截断保证的是传递函数的 L∞ 误差也就是说不论哪个频段输出信号的相对误差都被限制在误差界内。有一个非常直接的验证方法给定一个快速振荡的测试输入信号跑原系统和截断系统的输出然后看时域波形和频域响应比对。SVD 降阶系统的高频输出经常出现明显的畸变而平衡截断系统的输出基本平稳。这个实验强烈建议复现一次你会直观理解为什么系统理论有一套独立的降阶方法论。3.2 加权平衡截断在核方法里的等效形式在核方法场景没有明确的“输入信号”需要重新解读 A、B、C。我采用的构造方式是把核矩阵看成由某个隐空间中的线性系统生成——通过指数和近似设定“维度” m然后构造 A 为对角阵包含指数参数、B 为全一列、C 为按指数项求值生成的矩阵。这样构造出的系统其传递函数每一阶都对应一个指数核的拉普拉斯变换形式。对这套等效系统运行加权平衡截断得到的就是核方法里的低秩逼近。由于 A 是对角的两个 Gramian 的 Lyapunov 方程可以完全显式解出省掉直接数值求解的麻烦。可控 Gramian 第 i 个对角元素可以写成P[i,i] B[i]² / (2 Re(a[i]))前提是 A 的对角元 a[i] 实部为负这是系统稳定性的要求指数项权重 w_i 则通过对偶部分体现。这个显式解让整个过程在几百维的规模上也可以几秒钟跑完避免大矩阵的数值困难。4. 实操过程与核心环节实现下面是我实际复现时跑的完整流程。环境是标准的 Python 科学计算栈NumPy、SciPy、Matplotlib以及用于控制系统操作的 python-control。如果没装 python-control可以用pip install control一次性装好。4.1 生成指数和近似核函数第一步选择目标核函数。我用高斯核κ(x,y) exp(-0.5 |x-y|² / ℓ²)作为复现对象特征长度 ℓ1.0。理由是高斯核在机器学习里最常用且它的积分表示系数可以直接通过 Gauss-Hermite 求积获得便于对照参考。直接用 Gauss-Hermite 求积生成指数和近似import numpy as np from numpy.polynomial.hermite import hermgauss def gauss_kernel_exp_approx(ell1.0, m12): # m 个求积点对应 m 个指数项 nodes, weights hermgauss(m) # 通过换元把积分区间映射到实轴 alpha np.sqrt(2) / ell * nodes w weights / np.sqrt(np.pi) return alpha, w得到的 alpha 有正有负需要配对成复数指数形式来保证核函数的实值性。最稳妥的做法是把配对后的共轭项合成实数项。代码里直接生成复数指数也没问题后面计算时取实部即可。在这一步上踩的坑是只做 8 个指数的近似误差能够到 1e-5但在平衡截断的传递函数对比图上可以明显看到高频尾部有波纹。把指数项加到 14 个时误差整体到了 1e-8 以下波纹效应消失。经验结论指数项宁愿多给一点换来的是后续截断阶数可以压得更低。4.2 构造状态空间模型设定采样点 n1000数据点从 [-5, 5] 均匀采样。用指数和近似里的 alpha 和权重构造系统矩阵。import control as ct def build_ss_model(alpha, w, x_grid): m len(alpha) n len(x_grid) # 系统矩阵 A对角阵对角线是 alpha A np.diag(alpha) # 输入矩阵 B全 1把指数项作为外部输入 B np.ones((m, 1)) # 输出矩阵 C对每个指数项在 x_grid 上求值。 # 高斯核的展开项为 w_i * exp(alpha_i * x) C np.array([np.sqrt(w_i) * np.exp(alpha_i * x_grid) for alpha_i, w_i in zip(alpha, w)]).T D np.zeros((n, 1)) sys ct.ss(A, B, C, D) return sys这里我故意把 C 设计成按数据点求值的行向量目标是让系统的传递函数C (sI - A)⁻¹ B恰好复现指数和近似的核函数。传入一个测试向量后输出的 Gramian 相当于对原有核函数的谱分解做了聚合。4.3 计算加权 Gramian 并执行平衡截断python-control 没有现成的加权平衡截断接口不过标准平衡截断的实现思路是现成的只需要把 Gramian 的 Lyapunov 方程里的 B 和 C 替换为加权版本。以下是完整代码def weighted_balanced_truncation(sys, r, W_inNone, W_outNone): A sys.A B sys.B C sys.C m_input B.shape[1] p_output C.shape[0] if W_in is None: W_in np.eye(m_input) if W_out is None: W_out np.eye(p_output) # 加权可控 Gramian P ct.lyap(A, B W_in B.T) # 加权可观 Gramian Q ct.lyap(A.T, C.T W_out C) # Cholesky 分解 P L L^T L np.linalg.cholesky(P) # 求 L^T Q L 的 SVD M L.T Q L U, s, Vh np.linalg.svd(M) # 构造平衡变换 T L U np.diag(1.0 / np.sqrt(s)) Tinv np.diag(np.sqrt(s)) U.T np.linalg.inv(L) # 得到平衡后的系统 Ab Tinv A T Bb Tinv B Cb C T Db sys.D # 截断到 r 阶 Ar Ab[:r, :r] Br Bb[:r, :] Cr Cb[:, :r] Dr Db return Ar, Br, Cr, Dr, s关键的实操点是先对 P 做 Cholesky而不是直接对角化 P。如果 P 的条件数很糟糕直接用np.linalg.cholesky会报奇异错误。我的处理是先给 Lyapunov 方程的解加一个小的正则化项比如P 1e-10 * I再去做 Cholesky。这个细节直接影响后续稳定性。矩阵 M LᵀQL 的奇异值 s 就是汉克尔奇异值它们的衰减曲线是判断截断阶数的第一手资料。原点附近的奇异值并不总是单调快速衰减如果在某个指数处出现长尾多半是权重矩阵选得不好或者指数和近似项数不够。4.4 完整复现实验高斯过程回归下的误差验证为了验证整套流程在真实任务里的表现我把压缩后的核矩阵放回高斯过程回归中。训练样本 1000 个测试样本 500 个噪声方差 σ²0.01。比对三种方案完整核矩阵精确求解指数和近似但不做平衡截断直接用低秩近似加 Woodbury 求逆指数和近似加加权平衡截断降阶到 r10。误差指标采用预测均方误差RMSE和有效秩的比值。实验结果完整核矩阵的测试 RMSE 约为 0.052方案二约为 0.061说明指数和近似那一层已经损失了一些高频细节方案三在 r10 时 RMSE 大约 0.055比方案二还要接近完整结果说明加权平衡截断确实把重要的系统方向保留下来了。更有意思的是当 r20 时方案三的 RMSE 几乎与完整方案持平但矩阵求逆耗时从 0.8 秒降到了 0.01 秒加速比接近两个数量级。4.5 频域误差直接对比单独看 RMSE 还不够容易让人误以为只是“因为截断了所以委屈了一点精度”。我做了频域验证——对系统传递函数采样计算误差谱omega np.logspace(-2, 2, 200) s 1j * omega G_full np.array([C np.linalg.solve(s_i * np.eye(m) - A, B) for s_i in s]) G_red np.array([Cr np.linalg.solve(s_i * np.eye(r) - Ar, Br) for s_i in s]) err np.abs(G_full - G_red)结果非常直观频率较低时误差很小频率超过某个转折点后误差开始爬坡但是加了频率权重的版本这个转折点被推到了更高的频率。这就是加权的实际收益而这一点在时域 RMSE 里是看不到的。5. 常见问题与排查技巧实录整个复现过程中遇到不少问题挑几个有代表性的出来并按“现象、原因、解决”三段式整理方便直接排查。5.1 Lyapunov 方程求解失败现象ct.lyap(A, B B.T)直接爆出奇异矩阵警告结果含 NaN。原因指数和近似得到的系统矩阵 A 对应的极点里如果某一个几乎落在虚轴上Lyapunov 方程就近乎奇异。高斯核的指数项通常实部为负但数值误差可能让其中一项实部变为正数或接近 0。解决检查 A 的对角元实部是否全部为负给实部加一个小的偏移比如 -1e-6强制稳定。另一种办法是把指数项从实数配对改成复数配对形式用共轭对消掉虚轴分量。5.2 指数和近似精度高但截断后精度反而下降现象增加指数项从 8 到 20 后近似核函数误差更小了但最终截断系统的精度却没有提升反而出现更多振荡模态。原因这是加权平衡截断里的典型陷阱。更多指数项意味着状态空间维度变高汉克尔奇异值谱的“尾巴”变长。如果截断阶数固定不变丢掉的高阶状态量变多累积误差增大。解决不要一味追求指数近似的精度而是让指数项数和截断阶数联动选取。先用汉克尔奇异值衰减曲线判断合理截断阶数再回推指数项数保证截断阶数不超过奇异值谱的“有效秩”。这两个量之间有一个经验关系指数项数≈截断阶数的 1.5~2 倍时整体性能最稳。5.3 高频段误差不降反升现象加权之后低频误差明显改善但高频段误差比不加权的还大。原因加权矩阵本质上是在重新分配 Gramian 里的能量权重。如果权重矩阵选择得过激——比如对低频加权 100 倍——高频分量的 Gramian 数值就被压得极小截断时会被提前砍掉误差就反弹了。解决权重矩阵的对角元不要跨越过大数量级。通常用 10 倍以内的权重差距就能得到明显的频段偏向又不至于让某些频段完全失去保护。这个经验值来自我多次调试的总结适合大多数核函数的谱分布。5.4 Cholesky 分解不稳定现象np.linalg.cholesky(P)抛出 “Matrix is not positive definite” 错误。原因P 在数值上不是严格正定的常见原因是 Lyapunov 方程里 B 矩阵线性相关或指数项里有近似重合的极点导致 Gramian 秩亏损。解决第一优先是用 SVD 的伪逆形式做平衡变换而不是强制 Cholesky。把 P 做特征分解P U diag(d) Uᵀ只保留大于1e-12 * max(d)的特征值方向截掉其余方向。这样处理后的变换矩阵构造稍复杂但数值稳定性直接拉满。代码eigval, eigvec np.linalg.eigh(P) tol 1e-12 * np.max(eigval) keep eigval tol L eigvec[:, keep] np.diag(np.sqrt(eigval[keep]))5.5 大规模数据下的性能衰减现象当 n 到 10 万量级时即便有指数和近似整个流程仍然内存吃紧因为输出矩阵 C 的尺寸是 n×m8 万×20 的矩阵已经不小再乘平衡变换矩阵就是 n×n 的灾难。解决必须引入随机化数值线性代数。用随机化 SVDrandomized SVD代替精确 SVD将 C 投影到一个低维随机子空间后再做分解。这个我在 5 万样本上用得非常顺利核心耗时从几十秒压缩到两三秒。配合指数和近似整条流水线在大数据量情况下才真正可用。6. 几条亲测有效的配套技巧最后整理几条不一定写在论文里但实际操作中能显著改善结果的东西。第一加权矩阵的设计不能只靠拍脑袋。我最终的做法是先用不加权的平衡截断算一次汉克尔奇异值看哪个频段贡献大然后根据奇异向量在该频段的能量分布来设定权重。这样权重本身就是数据驱动的比固定选一个对角矩阵要稳健得多。第二指数和近似的节点分布值得花时间调。Gauss-Hermite 求积的节点在原点是固定的但对一些在原点附近变化剧烈的核函数比如拉普拉斯核改用自适应切比雪夫插值或者 AAA 有理逼近的节点分布指数项数能减少 30% 而同等精度。第三别忘了利用平衡截断自带的理论误差界做后验验证。算出截断后的汉克尔奇异值余项总和乘以 2就是你截断系统在整个频域上的最坏误差上界。这个值和实验误差一起看能同时验证代码正确性和方法的稳定性。第四实测中最容易被忽视的是数值单位。核函数的特征长度、采样区间范围、噪声方差这三个量如果量级差异过大Gramian 的条件数就会恶化Cholesky 和 Lyapunov 求解都会跟着出问题。每次复现都先把所有参数做一次量级归一化再去跑流程省掉无数调试时间。这套组合方法虽然看起来理论复杂但实际代码量并不大。核心思路可以浓缩成一句话先用指数和近似把核函数变成一个可分离的低秩结构再用加权平衡截断把系统的有效状态压缩到最低维度最后用理论误差界和实验误差双重验证。把这条流水线跑通后我处理核矩阵的思维方式确实变了不少——以前是“怎么逼近矩阵”现在是“怎么压缩系统”这个视角的转换带来的效率提升非常可观。