1. 先搞清楚SFA到底在解决什么问题我第一次真正动手折腾SFA慢特征分析Slow Feature Analysis的时候其实是被一个很实际的需求逼的手里有一批多变量时间序列想要提取出随时间变化最缓慢的、带有语义层级的底层信号。当初第一反应是用PCA降维但结果不太对味——PCA只按空间方差排序时间维度上的“变化快慢”它完全不管。比如一个快速抖动的噪声信号如果方差够大照样会被排在主成分里这对后续的时序建模毫无帮助。SFA恰好是反过来思考它要找的不是“方差最大的方向”而是“随时间变化最快慢的方向”。核心思想一句话就能讲清楚给定一个时间序列通过一个映射把原始高维信号变成一组标量使得这些标量随时间的变化尽可能平缓。这种技术在动作识别、表情分析、生理信号处理、甚至视频自运动估计里都有非常漂亮的应用而且理论推导干净利落代码实现量也不大。这篇文章会从数学原理讲到代码实现再到一个可以动手跑的合成实验最后把我实际踩过的一些坑也一并发出来。适合正在做时间序列特征提取、无监督表示学习或者想给分类模型加一点“时序先验”的读者。你只要有基本的线性代数底子能把numpy用顺读懂后面代码完全没问题。2. 把“变化慢”翻译成数学语言2.1 慢度怎么量化给一个时间信号 (y(t))怎么衡量它“变得快不快”最直接的做法是看它的一阶导数的平方在时间上的平均[ \Delta(y) \langle \dot{y}^2 \rangle_t ]这个值的含义非常直观如果信号是一条水平直线导数处处为0那 (\Delta(y)) 就是0它是终极的慢特征如果信号是高振幅的高频锯齿波导数平方的时间平均会很大。实际操作中我们拿到的都是离散采样点所以导数就得用差分去近似。常用的是一阶前向差分[ \dot{y}(t) \approx y(t1) - y(t) ]然后对整个时间序列取平均。这个近似虽然简单但后面会讲到它却藏着不少容易踩的坑。SFA的优化目标就是让 (\Delta(y)) 尽量小同时还要满足一些约束否则问题就没有任何意义。2.2 三个约束为什么一个都不能少如果只让 (\Delta(y)) 最小最简单的方法是把 (y(t)) 设成常数0这样导数恒为0目标函数直接归零。但这个结果没有任何信息量。所以必须加约束。第一个约束是零均值[ \langle y \rangle_t 0 ]这相当于去掉直流分量避免SFA把常数项当成最有价值的特征。第二个约束是单位方差[ \langle y^2 \rangle_t 1 ]这一步很关键。它告诉算法你要在“输出信号有一定能量”的前提下让变化最慢而不是通过把幅度压到0来投机取巧。单位方差约束让SFA提取的特征有明确尺度不同特征之间的慢度也可以直接比较。第三个约束是去相关[ \langle y_i y_j \rangle_t 0 \quad (i \neq j) ]也就是说第二个慢特征必须和第一个慢特征在时间上是正交的。否则你解出来的所有特征可能都是同一个方向的不同缩放版本。这个去相关约束让每个特征提供新的信息也保证了提取结果是按慢度从慢到快依次排列的一族正交信号。2.3 线性情形下怎么变成矩阵特征问题如果限定映射是线性的也就是[ y(t) w^T x(t) ]那么目标函数可以写成[ \Delta(y) \langle (w^T \dot{x})^2 \rangle_t w^T \langle \dot{x} \dot{x}^T \rangle_t w w^T B w ]其中 (B \langle \dot{x}(t) \dot{x}(t)^T \rangle_t)我们叫它差分协方差矩阵。零均值约束在数据已中心化时自然满足单位方差约束变成[ \langle (w^T x)^2 \rangle_t w^T \langle x x^T \rangle_t w w^T A w 1 ]其中 (A \langle x(t) x(t)^T \rangle_t) 是输入数据的协方差矩阵。于是问题转化成了一个标准的广义特征值问题[ B w \lambda A w ]特征值 (\lambda) 就是对应特征向量的慢度指标。把所有特征值从小到大排序最小的一个特征值对应的特征向量就是最慢的方向第二个小的对应次慢的方向依此类推。2.4 白化步骤的实质广义特征值问题虽然可以直接解但实现上我们通常更喜欢平凡的特征值分解。办法是先对输入做白化处理让 (A) 变成单位矩阵。具体流程是先对 (A) 做特征分解[ A U \Lambda U^T ]然后构造白化矩阵[ W U \Lambda^{-1/2} ]白化后的信号[ z(t) W^T x(t) ]它的协方差矩阵近似为单位阵。这时再对 (z(t)) 做慢度分析目标函数就变成[ \Delta(y) w_z^T B_z w_z, \quad B_z \langle \dot{z} \dot{z}^T \rangle_t ]约束条件退化为 (w_z^T w_z 1)直接对 (B_z) 做普通特征分解取最小特征值对应的特征向量即可。这个“先白化再分解差分协方差”的思路非常像PCA的白化预处理区别在于PCA只关心空间方差SFA在第二步关心的是时间差分协方差。代码实现时我强烈建议按这个流程走逻辑清晰也容易排查问题。3. 用Python手写一个SFA类3.1 整体设计思路我实现SFA时习惯把它封装成一个类似sklearn接口的类。核心方法就两个fit(X)输入形状为(T, d)的时间序列T是样本数d是维度。内部完成中心化、白化、差分协方差分解保存所有变换矩阵。transform(X)把新的时间序列映射到慢特征空间。设计类的好处是代码复用时不用每次重写整套流程。另外为了避免依赖具体的第三方库这个实现只用numpy保证任何环境都能直接跑。我还会留一个degree参数用来做非线性扩展。degree1就是纯线性SFAdegree2表示先对输入做二次多项式扩展再走线性SFA流程。这样文章后面演示非线性扩展时可以直接复用同一个类。3.2 核心代码与注释这里给出一个可直接运行的实现。注释尽量写得详细方便初学者逐行对照。import numpy as np from numpy.linalg import eigh def polynomial_expand(X, degree1): 对输入做多项式扩展。 参数: X: 形状 (T, d) 的时间序列 degree: 扩展次数1表示只保留原始特征 返回: 扩展后的特征矩阵形状 (T, D) 注意: 这里不包含常数项因为中心化会把常数项变成0。 T, d X.shape cols [X] # 原始线性项 if degree 2: for i in range(d): # 二次项x_i * x_j包含 ij避免重复 for j in range(i, d): cols.append((X[:, i] * X[:, j]).reshape(-1, 1)) return np.hstack(cols) class SFA: def __init__(self, n_components1, degree1, epsilon1e-8): self.n_components n_components self.degree degree self.epsilon epsilon # 白化时的正则化常数 def fit(self, X): # 1. 多项式扩展 X_exp polynomial_expand(X, self.degree) T, D X_exp.shape # 2. 中心化 self.mean_ X_exp.mean(axis0, keepdimsTrue) Xc X_exp - self.mean_ # 3. 计算协方差矩阵并白化 # 用 (T-1) 做分母是标准无偏估计实际差别很小 cov (Xc.T Xc) / (T - 1) eigenvalues, eigenvectors eigh(cov) # 把过小的特征值钳制到一个下限避免除零和数值爆炸 eigenvalues np.maximum(eigenvalues, self.epsilon) # 白化矩阵 W U * Lambda^{-1/2} # 这里没有乘 U^T得到的是PCA白化乘上 U^T 则是ZCA白化 self.whiten_matrix_ eigenvectors np.diag(1.0 / np.sqrt(eigenvalues)) # 白化后的信号 Z Xc self.whiten_matrix_ # 4. 对白化信号做一阶差分求差分协方差矩阵 dZ np.diff(Z, axis0) diff_cov (dZ.T dZ) / (T - 1) # 5. 对差分协方差矩阵做特征分解 slow_vals, slow_vecs eigh(diff_cov) # eigh 在不同版本里特征值排列顺序可能不同显式排序最稳妥 order np.argsort(slow_vals) self.slowness_ slow_vals[order] self.rotation_ slow_vecs[:, order][:, :self.n_components] # 6. 保存整体线性映射便于分析权重 self.components_ self.whiten_matrix_ self.rotation_ return self def transform(self, X): X_exp polynomial_expand(X, self.degree) Xc X_exp - self.mean_ Z Xc self.whiten_matrix_ Y Z self.rotation_ return Y3.3 代码里的几个关键抉择第一为什么用eigh而不是eig因为协方差矩阵是实对称矩阵eigh专门针对对称矩阵做了优化数值稳定性更好特征向量也能保证正交。如果用了eig碰到数值误差可能导致特征向量不正交后面的投影就会积累偏差。第二白化时为什么要对特征值做下限钳制训练数据里如果有高度相关的维度协方差矩阵会出现接近0的特征值。直接做1/sqrt(λ)会让该方向上的信号被放大到离谱的程度噪声也跟着被放大成“伪慢特征”。加上一个epsilon之后数值上稳定得多。这个值和数据尺度有关通常取1e-8到1e-6都可以。第三为什么用一阶差分而不是中心差分一阶差分实现简单延迟只有半个采样点对周期特征的相位影响也小。如果用中心差分(x[t1] - x[t-1]) / 2频率响应会更准但边界处理比较麻烦。对小规模实验来说一阶差分已经足够。第四特征向量的符号是不确定的。特征分解保证方向但不保证正负号这意味着取出的慢特征可能是cos(theta)也可能是-cos(theta)。这个是正常现象对比信号时不能用“完全相等”去判断应该看绝对相关度。4. 用一个旋转实验验证SFA真的能找到“慢特征”4.1 实验设计思路SFA教科书里最经典的验证场景是匀速旋转的点。假设一个点以恒定角速度在圆上运动观测信号是两个坐标 (x(t)) 和 (y(t))[ x(t) r\cos(\theta(t)) noise_x ] [ y(t) r\sin(\theta(t)) noise_y ]当角度随时间匀速增加时这两个坐标本身是快速振荡的正弦信号从频谱上看没有明显的“慢”属性。真正和角度 (\theta(t)) 相关的慢特征其实是 (\cos(\theta)) 和 (\sin(\theta)) 这一对正交分量。SFA的任务就是在只有坐标观测的情况下把这俩分量自动找出来。如果线性SFA能做到这一点就说明算法确实从时间维度上识别出了信号结构而不是靠空间方差来排序。4.2 实验代码我建议把整个实验写成脚本跑完之后画图观察。import matplotlib.pyplot as plt # 生成数据4个完整圆周均匀采样2000个点 n_points 2000 t np.linspace(0, 4 * np.pi, n_points) r 1.0 theta t x r * np.cos(theta) 0.03 * np.random.randn(n_points) y r * np.sin(theta) 0.03 * np.random.randn(n_points) X np.column_stack([x, y]) # 训练线性SFA提取前两个慢特征 sfa SFA(n_components2, degree1) sfa.fit(X) Y sfa.transform(X) slow_vals sfa.slowness_ print(特征值排序:, slow_vals[:6]) # 画图比较 plt.figure(figsize(10, 4)) plt.subplot(1, 2, 1) plt.plot(t, X[:, 0], alpha0.6, labelx) plt.plot(t, X[:, 1], alpha0.6, labely) plt.title(原始信号) plt.legend() plt.subplot(1, 2, 2) plt.plot(t, Y[:, 0], labelSFA feature 1) plt.plot(t, Y[:, 1], labelSFA feature 2) plt.plot(t, np.sqrt(2) * np.cos(theta), --, labelsqrt(2)*cos(theta)) plt.title(SFA慢特征) plt.legend() plt.tight_layout() plt.show()运行这段脚本你会看到两个慢特征曲线和 (\sqrt{2}\cos(\theta))、(\sqrt{2}\sin(\theta)) 几乎重合。也就是说SFA从高频振荡的坐标信号里成功恢复出了旋转角的正交慢变量。按照第2节的推导如果采样正好覆盖整数个周期理想慢特征的特征值应该约等于角速度平方。在这个例子里角速度就是1因为 (\theta t)所以特征值应该接近1。如果你跑出来的值在1附近就说明实现没有错。4.3 实验解读和延伸这个实验最大的意义在于SFA不依赖任何标签信息只靠“随时间变化慢”这个隐含约束就发现了数据背后真正的生成变量。这在无监督学习里是非常可贵的性质。不过有两点需要提醒。第一如果加了大噪声比如把噪声标准差调到0.5SFA提取的特征仍然会比较光滑但特征值会明显偏大因为差分运算会把噪声的高频分量放大。第二如果旋转角速度不是恒定的而是时变的SFA提取出来的特征就不再是严格的正余弦而是根据角速度变化形态调整的某种非线性信号。对于真实数据来说这种情况反而更常见。5. 从线性到非线性二次扩展、核化与分层SFA5.1 多项式扩展怎么用线性SFA只能找到原始特征的线性组合表达能力有限。解决思路非常朴素先把输入映射到非线性特征空间再做线性SFA。最简单的非线性扩展就是多项式扩展。前面代码里的polynomial_expand已经实现了二次多项式扩展。对于二维输入 (x(x_1,x_2))degree2时扩展后的特征空间包含[ x_1, x_2, x_1^2, x_1x_2, x_2^2 ]共5个特征。在这个空间里做线性SFA得到的慢特征本质上就是原输入的非线性函数。用法很简单sfa_nonlinear SFA(n_components3, degree2) sfa_nonlinear.fit(X) Y_nonlinear sfa_nonlinear.transform(X)有一点必须注意多项式扩展之后维度涨得很快。如果原始维度 (d10)二次扩展后维度就有 (10 55 65)。白化需要估计65维的协方差矩阵序列长度T至少要是维度的几倍才稳定否则特征分解出来的结果意义非常存疑。我一般要求 (T \ge 10 \times D)如果达不到就老老实实加大正则化项并且只看前几个大特征值对应的特征。另外二次扩展里往往会出现接近常数的组合方向。比如匀速旋转的圆的 (x^2 y^2) 接近常数这个方向在白化之前方差极小。处理不好时SFA会把一个“方差被正则化强撑起来、实际变化几乎为0”的方向排到最慢特征里看起来像一条直线。这不是算法错而是模型退化。解决办法是不要迷信排名最靠前的特征要结合原始信号去看这个慢特征的方差和波形是否有意义。5.2 核化SFA的思路多项式扩展只是非线性化的一种手段。另一个常见思路是核化用一个核函数隐式地把数据映射到高维再生核希尔伯特空间然后在那个空间里求解SFA问题。核化SFA的计算套路和核PCA很接近先算核矩阵 (K)再做中心化和白化最后对差分核矩阵做特征分解。好处是不需要显式构造多项式特征可以用RBF这类核捕捉更复杂的非线性关系代价是核矩阵是 (T \times T) 的当样本数T上万时非常吃内存训练时间也上去了。如果你只是做中小规模实验二次多项式扩展已经足够如果数据本身复杂比如像素级的视频帧再考虑核化不迟。5.3 分层SFA把慢特征再送进去一次单个SFA提取出的慢特征往往还停留在“信号级”。想要更高层的语义表示可以把第一轮的慢特征作为输入再做一次SFA。这样形成一种层次结构每层提取的时间尺度越来越长特征越来越抽象。这种思路最早启发自视觉皮层对信息的分级处理。低层SFA提取图像中某个局部区域的光流或亮度变化高层SFA再把这些局部慢特征整合成语义级别的整体运动模式。做动作识别时这种分层结构比单层SFA要稳定得多。实现分层SFA并不复杂第一层SFA训练完把输出拼接上可选的原始特征作为第二层SFA的输入往复循环。需要注意每一层都要做独立的中心化和白化层数多了以后数值误差会累积所以中间可以时不时检查一下特征值分布确认数据没被破坏。6. SFA在实际任务中能干什么SFA最大的特点是不需要标签纯靠时序平滑性就能学到有结构的表示所以它很适合那些“标签很贵、时序信号很丰富”的场景。视频动作识别是典型例子。把视频拆成连续帧提取每帧的姿态关键点坐标或者光流特征拼接成一条高维的时间序列。动作从开始到结束在特征空间里必然有一个渐变过程SFA提取出的慢特征正好对应动作阶段的底层状态。之后拿这些慢特征做分类不仅维度低还天然过滤掉了帧间抖动。生理信号处理也是SFA的主场。脑电、心电、肌电这类信号要么基线漂移缓慢要么伪迹突变快速。SFA可以把不同时间尺度的分量按慢度分离开慢特征对应基线漂移快特征对应噪声和伪迹。这在信号预处理里非常实用甚至可以用来做盲源分离。机器人自运动估计方面SFA也有经典的衍生应用。一个移动平台配备全景相机图像像素在运动时产生的视觉流和时间变化率之间有关系SFA可以学会提取“自身运动参数”对应的慢变化特征完全不需要标注。此外SFA还适合用在时间序列降维可视化。高维神经数据经过SFA后前两维往往能展现出清晰的轨迹结构比直接PCA可视化更能反映时间演变过程。我做过一个模拟神经群体活动数据的实验SFA的可视化结果明显比PCA多了一个“时间平滑”属性曲线连起来非常容易看。7. 踩坑记录SFA实现与使用中的高频问题7.1 采样间隔和差分方向SFA的慢度指标本质上是对时间导数的估算但离散差分的结果和采样间隔直接相关。如果采样间隔是 (dt)那么一阶差分应该除以 (dt)得到的才是导数量纲。很多人没注意到这点导致同一个物理信号在不同采样率下得到的“慢度”数值差很多。如果你的分析需要对不同采样率的数据做比较最好在计算差分时除以dt并且换算成“每单位时间的慢度”。7.2 数据长度与维度比SFA要估计两个矩阵输入协方差和白化信号的差分协方差。维度一高估计误差就会变大。经验上样本数至少要是输入维度的10倍以上取20倍更稳妥。对于非线性扩展后的高维特征这个比例更难满足所以不要盲目扩展。如果序列不够长我会先降维比如先PCA压缩到能用维度再做SFA。虽然PCA会丢掉一部分非线性信息但至少把协方差结构估计稳住了。7.3 特征符号和相位不确定特征分解出来的特征向量符号不唯一所以SFA特征可能整体翻转。数据每次重跑第一慢特征的相位也可能出现90度偏移因为两个简并特征之间可以做任意旋转。在做下游分析时不要直接比较两段实验提取的慢特征的原始值而是比较它们的包络、功率谱或者互相关绝对值。7.4 正则化项设不好会出幺蛾子白化时的epsilon设置很关键。太小协方差接近0的方向噪声爆炸太大会把低方差方向强行放大产生虚假慢特征。常见做法是先用训练集统计协方差特征值的量级把epsilon设为最大特征值的1e-6到1e-4。如果数据已经做过标准差归一化1e-8一般够用。7.5 伪慢特征怎么识别SFA给出的第一个特征不一定是“有意义的慢特征”。如果最慢特征看起来像一条几乎水平的直线而且它的原始方差极小那很可能就是常数方向或退化方向。真正值得用的慢特征应该有非零波动、形态平滑、且和原始时间序列的某些事件对得上。我习惯把前3个特征的波形都画出来看一眼再做下一步决策。还有一个常见现象数据里如果有一个明显的趋势项比如传感器温漂趋势会成为天然的最慢特征SFA会把所有注意力放在它上面反而忽略了周期性慢变量。这不是算法问题而是任务目标问题。如果趋势不是你要的信号先用高通滤波或差分把趋势去掉再跑SFA。8. 后续可以继续玩的方向SFA这套“先白化再分解差分协方差”的思路其实可以和很多方法结合。比如把慢特征当作自监督学习的一个正则项加到现有的表示学习损失函数里鼓励模型学到时序平滑的隐变量。这在处理连续控制任务时尤其有用因为稳定、平滑的状态表示往往更容易学习和泛化。另一个玩法是慢特征回归。SFA提取的特征天然和时序关联强用它们去做未来趋势预测比用原始高维输入更稳也更能抵抗噪声。我在处理一段带强噪声的传感器数据时发现SFA特征作为回归输入的误差比原始信号作为输入低了不少主要原因是慢特征已经把高频噪声分量丢掉了。最后可以沿着分层的路继续走多个尺度的SFA拼接然后再接浅层分类器就是一个很轻量的时序分类pipeline。不用深度学习框架用numpy就能完成整条链路。这也是我觉得SFA最适合入门的位置——代码短、理论清晰、效果看得见。希望你也能从中发现一点新的用法。
