最近在整理动力学系统的资料时又重新把“从数据里找方程”这个方向完整过了一遍。很多人第一次听到“动力学重构”会觉得玄乎其实说白了就是手里只有一条或几条时间序列观测数据想反推出系统背后那条微分方程。这不是纯理论问题而是很多真实工程场景的刚需——化学反应中间变量测不全、传感器只有某些状态量的读数、仿真模型太复杂想简化……诸如此类。这篇博文想聊的是论文泛读里一个非常经典的技术路线神经网络自动编码器Autoencoder配合 SINDySparse Identification of Nonlinear Dynamics稀疏非线性动力学辨识一起完成动力学重构并附上可运行的 Python 源码思路。这套方法的论文原型是 2019 年前后 Champion、Lusch、Kutz、Brunton 那篇关于“dynamic autoencoder”的工作核心思想一句话就能概括先用神经网络找到一个隐藏坐标系让系统在这个坐标系下的动力学表达式尽量稀疏、尽量简单再用 SINDy 把这个方程“读”出来。非常适合做时间序列分析、系统辨识、物理仿真、甚至神经网络科学计算方向的朋友参考复现。1. 先搞懂一个问题动力学重构究竟在解决什么1.1 一条轨迹可以有无数个方程所以这是个不适定问题先给你一个直觉层面的例子。假设我在一个二维平面上画了一条圆的轨迹问你这个轨迹对应哪个动力学方程你大概率会回答[ \dot{x} -y,\quad \dot{y} x ]这是最简单的匀速圆周运动。但如果我告诉你其实真实系统是[ \dot{x} -y \epsilon x^3 ]在 (\epsilon) 非常小、并且你只观测了很短一段轨迹时两条方程产生的曲线几乎叠在一起。你根本没法从有限数据里区分到底哪个是“真”的。这就是动力学重构的第一个难点它不是简单地“拟合一条曲线”而是要在无数个能产生相似轨迹的候选方程里挑出一个最合理、最可解释的。靠纯黑箱拟合解决不了这个问题因为拟合出来的参数再多也很难告诉你说“这个系统就是这三个非线性项在起主导作用”。1.2 从“拟合”到“发现”稀疏辨识带来的视角转换传统系统辨识的做法是先假设模型结构比如“我猜这是一个二阶线性系统”然后去估计里面的系数。问题是模型结构一旦猜错后面所有结果都白搭。SINDy 换了个思路先别急着定结构而是给系统准备一个内容足够丰富的候选函数库把常数项、一次项、二次项、三次项、甚至三角函数项都摆进去然后让回归算法在里面挑出真正有用的那一小部分。绝大多数项的系数是 0只有少量非零系数这个“稀疏解”本身就是模型结构。你可以把它理解成侦探破案先列出一堆可疑人员候选函数经过排查之后锁定三五个关键人物非零系数。其他的要么排除要么不太重要。这种思路的好处是你不用再人工猜测“它是不是多项式系统”“是不是含三角函数”只要候选库里包含真实项稀疏回归就有机会把它找出来。1.3 坐标问题方程复杂不复杂取决于你站在哪里看直接 SINDy 有个隐性前提你手上的观测坐标差不多是“好坐标”。什么叫好坐标对经典 Lorenz 系统来说直角坐标系 (x,y,z) 就是好坐标因为它的方程在三次多项式库里只有几个非零项非常稀疏。但现实里传感器量测的坐标往往不是“上帝视角”的好坐标而是各种状态量的线性或非线性混合。一旦坐标系取得不好原本很简单的动力学方程会变得极其臃肿。举个例子还是 Lorenz 系统如果我把观测变量改成 (p x^2 y^2)、(q e^z) 这种非线性组合那么 (p,q) 满足的微分方程会包含大量的高次交叉项SINDy 根本筛不干净。这正是引入自动编码器的原因让网络自己学着把观测数据映射到一个隐空间里并在这个隐空间上要求动力学方程尽量稀疏。换句话说自动编码器负责换坐标系SINDy 负责在换好坐标系后找方程两者是分工协作的关系。2. SINDy 的数学内核稀疏回归怎么“发现”方程2.1 候选函数库与线性化结构SINDy 的精髓在于把非线性动力学问题转化为线性回归问题做法是这样的。设系统状态为 (z(t) \in \mathbb{R}^d)可以是观测变量也可以是隐变量假设动力学可以写成[ \frac{dz}{dt} f(z) ]我们不知道 (f) 的形式但可以假设 (f) 能表示成一组候选基函数的线性组合。把候选基函数堆成一排向量[ \Theta(z) [1,\ z_1,\ z_2,\ \dots,\ z_1^2,\ z_1z_2,\ \dots,\ \sin(z_1),\ \dots] ]于是[ \frac{dz}{dt} \Theta(z),\Xi ]其中 (\Xi) 是一个系数矩阵每一列对应一个状态变量的方程。只要 (\Theta(z)) 里包含了真实的非线性项那么 (\Xi) 就会很稀疏——绝大多数位置是 0只有少数位置非零。这一步就是 SINDy 的“稀疏辨识”雏形。这里特别提醒一点为了能捕捉二次非线性动力学(z) 的候选函数库至少要把多项式阶数开到 2 或 3。Lorenz、Rossler、Lotka-Volterra 这类经典系统在三次多项式库内都能被非常简洁地表示出来。反过来如果你只给线性函数库SINDy 再厉害也识别不出非线性项所以候选库的完整性直接决定了结果上限。2.2 数值求导实操里最容易被忽略的误差源SINDy 的输入听起来很简单左边是导数右边是函数库矩阵。但“导数”从哪来真实传感器数据不会直接给你导数你需要从离散时间序列里数值估计。最常用的中心差分[ \dot{x}(t_k) \approx \frac{x(t_{k1}) - x(t_{k-1})}{2\Delta t} ]这个公式问题非常大如果信号里有高频噪声差分运算会把噪声放大到离谱的程度放大幅度大约是 (1/\Delta t)。假设采样时间 (\Delta t 0.01)噪声标准差是 0.01那么导数的噪声标准差可能直接飙到 1 左右比真实导数本身还大。实操中我一般先做两步预处理采样率要足够高至少保证每个振荡周期有 50 个点以上时序数据先做平滑比如 Savitzky-Golay 滤波或简单卷积平滑再计算导数。如果你用的是干净仿真数据中心差分通常够了。但在真实实验数据上SINDy 的识别效果受导数噪声影响非常大这一点后面遇到的坑多半都出在这里。2.3 稀疏回归与阈值迭代知道了数据矩阵之后剩下的就是解一个稀疏线性回归[ \min_{\Xi} \ \frac{1}{2}|\dot{Z} - \Theta(Z)\Xi|_F^2 \lambda|\Xi|_1 ]这看起来就是一个 Lasso 问题但 SINDy 原论文里更常用的一种做法叫STLSQSequentially Thresholded Least Squares顺序阈值最小二乘流程非常直观先用最小二乘求一个初始系数矩阵(\Xi \Theta^ \dot{Z})给定一个阈值 (\lambda)把绝对值小于 (\lambda) 的系数全部置零用剩下的非零列重新做一次最小二乘重复 2、3 步直到系数矩阵不再变化。阈值 (\lambda) 怎么选没有一个万能值。我的经验是先设一个相对较大的阈值比如 0.5把明显多余的项全部干掉然后逐步调小到 0.1、0.05观察哪些项稳定保持非零。你还可以用“系数稀疏程度 vs 拟合误差”的 Pareto 曲线来判断曲线拐点附近的阈值通常是比较好的选择。2.4 一个 30 行代码的 SINDy 最小实现我先给一个不依赖任何深度学习库的 SINDy 最小实现用 Lorenz 数据直接跑这样后面的自动编码器路线会更容易理解。import numpy as np from scipy.integrate import solve_ivp def lorenz(t, state, sigma10.0, rho28.0, beta8/3): x, y, z state return [sigma*(y-x), x*(rho-z)-y, x*y-beta*z] # 生成 Lorentz 轨迹 dt 0.001 t_span (0, 20) t_eval np.arange(0, 20, dt) sol solve_ivp(lorenz, t_span, [1.0, 1.0, 1.0], t_evalt_eval, rtol1e-9, atol1e-9) X sol.y.T # (N,3) # 中心差分求导 dX np.gradient(X, dt, axis0) def poly_library(Z, order3): n, d Z.shape cols [np.ones(n)] for deg in range(1, order1): # 生成所有次数为 deg 的单项式 if d 1: idx_list [(i,) for i in range(d)] # 用递归或 itertools 生成 d 维单项式索引这里简化为 d 较小场景 from itertools import combinations_with_replacement for comb in combinations_with_replacement(range(d), deg): col np.ones(n) for idx in comb: col * Z[:, idx] cols.append(col) return np.column_stack(cols) Theta poly_library(X) def stlsq(Theta, dX, threshold0.05): Xi np.linalg.lstsq(Theta, dX, rcondNone)[0] for _ in range(5): small np.abs(Xi) threshold Xi[small] 0 for j in range(Xi.shape[1]): if np.any(Xi[:, j] ! 0): idx Xi[:, j] ! 0 Xi[idx, j] np.linalg.lstsq(Theta[:, idx], dX[:, j], rcondNone)[0] return Xi Xi stlsq(Theta, dX, threshold0.05) np.set_printoptions(precision2, suppressTrue) print(Xi)这个代码跑完你会看到系数矩阵在对应位置出现 (10)、(-10)、(28) 等接近真实参数的数值而其他位置基本是 0。这就是 SINDy “发现方程”的过程——非常朴素但非常有效。3. 自动编码器补上关键一环找到合适的坐标系3.1 为什么直接 SINDy 常常不够如果你手上的观测本身就已经是简洁坐标SINDy 确实够用。但现实问题常常是观测维度高比如高速相机拍一个摆动物体一帧就是几百维像素观测值是多个物理量的混合真实“状态”不直接可见变量之间存在强耦合直接在每个观测维度上写方程候选库会爆炸。举个例子如果观测维度是 100 维哪怕只用二次多项式库候选函数数量也有大约 5000 个。每个变量的方程都要在 5000 个基函数里做稀疏回归计算量大不说还极容易把噪声项也学进去。这个时候需要降维。传统 PCA 可以做线性降维但真实动力学系统的低维流形往往是非线性的。自动编码器的优势就在这里——它能学习一个非线性映射把高维观测压缩到低维隐变量 (z)再从 (z) 解码回观测空间。3.2 关键结合点不仅要“重构得好”还要“动力学简单”但是如果只用一个普通自动编码器去做降维得到的隐变量不一定会让动力学方程变简单。普通 AE 的损失函数是[ \mathcal{L}_{recon} |x - D(E(x))|_2^2 ]它只要求解码后能还原原数据至于隐空间里的时序演化长什么样它完全不管。于是 AE 很可能学习到一个扭曲的隐空间重构误差很小但 (z(t)) 的时间演化极其复杂SINDy 在隐空间上直接失效。所以要把 AE 和 SINDy 联合训练总的损失函数变成[ \mathcal{L} |x - D(z)|_2^2 \lambda \left| \frac{dz}{dt} - \Theta(z)\Xi \right|_2^2 \lambda_1 |\Xi|_1 ]第一项是重构损失保证 (z) 能完整表示观测信息第二项是动力学损失要求隐变量满足 SINDy 给出的稀疏方程第三项是稀疏正则迫使 (\Xi) 尽量稀疏。这个联合训练的精妙之处在于网络可以自由变换坐标但变换后的坐标必须满足两个条件——既能重构回原来的观测数据又让隐空间里的方程变得极其稀疏。这在数学上等价于去搜索一个“能稀疏化动力学的坐标变换”。我自己的理解是这个过程有点像把一团乱麻的毛线整理好AE 负责把缠绕的线重新卷到线轴上SINDy 则负责告诉你这个线轴每转一圈线长变化遵循什么规律。两者缺一不可。3.3 完整前向流程整套方法跑起来时数据流是这样的拿到观测轨迹 (x(t))把长轨迹切成固定长度的小片段每个片段输入编码器 (E)得到隐变量轨迹 (z(t))用中心差分从 (z(t)) 算出隐变量的导数 (\dot{z}(t))构造候选函数库 (\Theta(z))用 SINDy 或者直接用梯度优化求解稀疏系数 (\Xi)用解码器 (D) 把 (z(t)) 映射回观测空间计算重构误差把重构误差和动力学误差合并反向传播更新编码器、解码器和 (\Xi)。训练完成之后你不只是得到一个能重构数据的自编码器更重要的是你拿到了一组明确的常微分方程(\dot{z} \Theta(z)\Xi)。以后你想预测系统未来行为不需要再跑神经网络直接积分这组方程就行——这才是真正的“发现方程”。4. PyTorch 源码拆解从数据生成到联合训练4.1 准备训练数据多轨迹、切片、归一化直接拿一整条长轨迹训练AE 很容易过拟合到初始条件的特殊路径上。我一般会准备 20 条以上不同初始条件的轨迹然后按固定窗口切成小样本。以 Lorenz 为例import torch import numpy as np from scipy.integrate import solve_ivp def generate_lorenz_data(n_traj20, traj_len5000, dt0.002, window128): data [] for _ in range(n_traj): y0 np.random.uniform(-10, 10, size3) t_eval np.arange(0, traj_len*dt, dt) sol solve_ivp(lorenz, (0, traj_len*dt), y0, t_evalt_eval, rtol1e-9, atol1e-9) data.append(sol.y.T) data np.stack(data, axis0) # (n_traj, traj_len, 3) # 按窗口切段 windows [] for seq in data: for i in range(0, traj_len - window, window//2): windows.append(seq[i:iwindow]) windows np.stack(windows, axis0) # 归一化 mean windows.reshape(-1, 3).mean(axis0) std windows.reshape(-1, 3).std(axis0) windows (windows - mean) / std X torch.tensor(windows, dtypetorch.float32) # (N, window, 3) return X注意归一化非常关键。Lorenz 的三个变量尺度差异不大但如果是真实系统变量间量纲差异可能很大不归一化的话高次多项式项会数值溢出或者某个变量主导整个损失。4.2 搭自编码器MLP 编码器/解码器网络结构不需要很深我用的是三到四层全连接import torch.nn as nn class Encoder(nn.Module): def __init__(self, obs_dim3, hidden_dim64, latent_dim3): super().__init__() self.net nn.Sequential( nn.Linear(obs_dim, hidden_dim), nn.Tanh(), nn.Linear(hidden_dim, hidden_dim), nn.Tanh(), nn.Linear(hidden_dim, latent_dim), ) def forward(self, x): return self.net(x) class Decoder(nn.Module): def __init__(self, latent_dim3, hidden_dim64, obs_dim3): super().__init__() self.net nn.Sequential( nn.Linear(latent_dim, hidden_dim), nn.Tanh(), nn.Linear(hidden_dim, hidden_dim), nn.Tanh(), nn.Linear(hidden_dim, obs_dim), ) def forward(self, z): return self.net(z)这里一个实用小技巧数据是按时间窗口过来的所以前向传播时给每个时刻的状态单独过 MLP然后把整个窗口拼接起来。在 PyTorch 里直接用Encoder(x)会默认作用在最后一维上所以形状(N, window, obs_dim)可以直接传进 MLP因为nn.Linear会作用在最后一维。4.3 构造 SINDy 候选库与导数计算隐变量 (z) 的候选库构造方式和前面原生 SINDy 一样只是这里要支持 PyTorch 的自动求导所以用张量操作实现from itertools import combinations_with_replacement def poly_library_torch(z, order3): # z: (N, T, d) n, T, d z.shape z z.reshape(-1, d) # 合并 N,T columns [torch.ones(n*T, 1, devicez.device)] for deg in range(1, order1): for comb in combinations_with_replacement(range(d), deg): col torch.ones(n*T, 1, devicez.device) for idx in comb: col col * z[:, idx:idx1] columns.append(col) return torch.cat(columns, dim1) # (n*T, n_terms)导数计算可以直接用torch.gradientdef compute_dz_dt(z, dt): # z: (N, T, d) dz torch.gradient(z, spacingdt, dim1)[0] return dz如果你的数据噪声不大torch.gradient就行如果噪声明显我建议先对 z 轨迹做平滑再算导数。4.4 联合训练损失函数与稀疏系数更新这里有个实现细节要特别提醒(z) 的导数一旦通过torch.gradient算出来它是依赖于z的计算图的可以正常反向传播。但 (\Xi) 的更新有两种常见方式直接把 (\Xi) 当成可训练参数用 Adam 优化器更新每隔若干轮对 (\Xi) 做一次硬阈值收缩其实也就是 STLSQ 的思路。我比较推荐混合方案Adam 负责更新编码器和解码器(\Xi) 则每隔 10 轮做一次阈值收缩。因为如果让 Adam 直接优化 (\Xi) 并加 L1 正则系数很容易全部缩到 0或者某个稀疏结构很难被“钉”住。class DynamicAutoencoder(nn.Module): def __init__(self, obs_dim3, latent_dim3, hidden_dim64, order3): super().__init__() self.encoder Encoder(obs_dim, hidden_dim, latent_dim) self.decoder Decoder(latent_dim, hidden_dim, obs_dim) self.order order # Thêta(z) 的项数 n_terms 0 for deg in range(order1): n_terms len(list(combinations_with_replacement(range(latent_dim), deg))) self.Xi nn.Parameter(0.1 * torch.randn(n_terms, latent_dim)) def forward(self, x): z self.encoder(x) x_hat self.decoder(z) return z, x_hat def train(model, X, dt0.002, epochs300, lr1e-3, lambda_dyn10.0, threshold0.05): optimizer torch.optim.Adam(model.parameters(), lrlr) for epoch in range(epochs): optimizer.zero_grad() z, x_hat model(X) # X: (N, T, obs_dim) recon_loss torch.mean((x_hat - X)**2) dz torch.gradient(z, spacingdt, dim1)[0] Theta poly_library_torch(z, model.order) # (N*T, n_terms) Xi model.Xi # 动力学损失: || dz - Theta Xi ||^2 dyn_loss torch.mean((dz.reshape(-1, z.shape[-1]) - Theta Xi)**2) sparse_reg torch.mean(torch.abs(Xi)) loss recon_loss lambda_dyn * dyn_loss 1e-3 * sparse_reg loss.backward() optimizer.step() # 每 10 轮做一次硬阈值收缩 if epoch % 10 0: with torch.no_grad(): small torch.abs(model.Xi) threshold model.Xi[small] 0.0 if epoch % 50 0: print(fepoch {epoch}, recon{recon_loss.item():.4f}, dyn{dyn_loss.item():.4f})这里有几个参数调起来很讲究。lambda_dyn是动力学损失在总损失里的权重我试过 0.1、1、10、100结论是太小的话 AE 学到一个重构很好但动力学极其混乱的隐空间太大的话 AE 的重构能力会被牺牲隐变量和原始观测之间的信息无法完全保持。一般我会先在 PCA 或普通 AE 预训练基础上从lambda_dyn 1开始试再看重构损失和动力学损失的平衡。4.5 训练后验证隐空间积分 解码回观测空间当你把模型训练完拿到了稀疏的 (\Xi)怎么确认它真的发现了系统方程最直观的验证方式是从某个初始隐状态出发用识别出的 ODE 做数值积分然后把整条轨迹解码回观测空间再和真实观测轨迹做对比。from scipy.integrate import solve_ivp def identify_rhs(z, Xi, order3): # 返回 Theta(z) Xi Theta_z poly_library_torch(torch.tensor(z.reshape(1, 1, -1), dtypetorch.float32), order).detach().numpy() return (Theta_z Xi.detach().numpy()).flatten() def predict_trajectory(model, z0, t_eval): Xi_np model.Xi.detach().numpy() sol solve_ivp(identify_rhs, (t_eval[0], t_eval[-1]), z0, t_evalt_eval, args(Xi_np, model.order), rtol1e-8, atol1e-8) z_pred sol.y.T z_tensor torch.tensor(z_pred, dtypetorch.float32) x_pred model.decoder(z_tensor).detach().numpy() return x_pred, z_pred如果识别正确预测出来的相位图会和真实数据的吸引子形状高度重合。如果只是重构损失很小、但预测轨迹很快漂移发散说明动力学损失没学好隐空间的微分方程是错的。5. 调参与踩坑记录复现过程中最值得注意的细节5.1 AE 发现的方程不一定是原始方程这是很多人第一次复现时最容易困惑的地方。你用 Lorenz 数据喂给 AESINDy最后得到的隐变量方程往往不是标准形式[ \dot{x} 10(y-x),\quad \dot{y}x(28-z)-y,\quad \dot{z}xy-\frac{8}{3}z ]而是某个可逆线性变换或者非线性变换之后的形式。比如隐变量可能是 (z_1 x y)、(z_2 x - y)、(z_3 z) 这样的组合方程系数完全不同但动力学性质等价。所以验证结果的时候不要执着于“系数必须等于真实系统”而要看两点隐空间 ODE 预测出的轨迹解码回观测空间后是否逼近真实轨迹识别出的方程是否稀疏非零项数量是否远小于候选库规模。只要轨迹对得上、方程又稀疏那这套方法就算成功了。坐标本身并不唯一允许任意可逆变换是这一类“表示学习动力学辨识”方案的理论基础。5.2 损失权重的“跷跷板”重构损失和动力学损失之间存在天然博弈。如果你把lambda_dyn调很大AE 会倾向于找到一条“最好拟合成稀疏方程”的隐空间轨迹但代价是 z 可能丢失一些重构所需的信息导致从隐变量解码回观测的时候面目全非。反之lambda_dyn太小动力学方程会变得特别复杂SINDy 很多非零系数都撂在那里根本无法解释。我建议的顺序是先单独用重构损失预训练 AE得到一组正常的编码器/解码器权重再加载预训练权重从lambda_dyn 1开始联合训练观察两个损失的比例如果重构损失开始飙升就降低lambda_dyn如果动力学损失不降就增大它。还有一种比较省心的做法动力学损失不直接对隐变量的导数做拟合而是对“积分之后的轨迹”做拟合——也就是说把整条 z 轨迹用识别出的 ODE 积分出来再与编码器给出的 z 轨迹做对比。这样能避免数值微分带来的高频噪声影响但计算量更大实现也更复杂。新手阶段先用导数拟合就够了。5.3 数据长度与状态空间覆盖SINDy 和 AE 都很吃数据质量。如果你的训练数据只是 Lorenz 吸引子上很窄的一小段系统从来没有经历过大幅振荡那么识别出的方程很可能只在小范围适用。换个初始条件预测几步误差瞬间爆炸。我一般会刻意让训练轨迹覆盖更大的状态空间初始条件在 ([-15, 15]^3) 区间里随机采样然后把多条轨迹混合训练。这比单纯把一条轨迹拉长更有效因为 SINDy 需要看到系统在不同区域的“反应”才能真正判断哪些非线性项是必需的。5.4 阈值、候选库尺度和数值稳定性候选函数库如果包含三阶甚至四阶多项式项隐变量稍微超出训练范围时库函数数值会迅速膨胀。比如 (z_1^3)当 (|z_1|10) 时是 1000 量级跟常数项 1 完全不在一个量级。所以训练时归一化很重要最好让隐变量基本落在 ([-2, 2]) 区间内。另外阈值不是越大越好阈值设置效果问题0.5 以上方程非常稀疏可能把真实非零项也删掉识别失效0.1 ~ 0.3多数情况下效果不错需要联合损失配合避免遗漏小系数项0.01 以下非零项很多几乎不稀疏难以解释我是从 0.1 起步看识别出的方程里非零项数量如果太多就升高阈值如果出现动力学损失完全不收敛就降低阈值。这个调参过程不复杂但必须结合轨迹验证一起看单看损失曲线很容易被骗。6. 常见问题速查表与排查思路6.1 常见问题速查表下面是我在调这套方法时经常遇到的问题整理成一张速查表方便你对照排查问题现象可能原因解决思路重构损失很低但预测轨迹发散发散动力学损失权重不够隐空间方程没学准增大lambda_dyn增加联合训练轮数动力学损失很低但重构损失很高编码器把信息丢了AE 退化成恒等映射附近的变形降低lambda_dyn增大隐变量维度增加编码器容量系数矩阵几乎全零稀疏阈值太大把真实项砍掉了调低阈值或者先增大预训练轮数系数矩阵一堆非零小系数阈值太小稀疏性不足调高阈值或增加 L1 正则权重隐变量轨迹高频抖动数值微分放大了编码器输出的微小波动对编码器输出做平滑或减小采样时间步dt训练到最后损失还在震荡学习率太高或 AE 结构和 SINDy 联合训练不稳定降低学习率加入学习率衰减增大 batch 或序列长度换了真实数据后完全失效真实数据噪声大候选库没有覆盖系统真实结构先做数据降噪换更高阶库或增加三角函数项6.2 排查顺序从损失到系数再到轨迹如果你复现时出了诡异问题别急着改代码。我一般是按下面这个顺序拆解先看重构损失。如果重构损失不降说明 AE 本身就没训练好后面的动力学都白搭再看动力学损失。如果重构很好但动力学损失下不去大概率是隐空间维度不够或者候选函数库缺了关键项打印 (\Xi) 矩阵看非零项是否稳定。如果每次训练非零项位置都不一样说明数据覆盖不足或者阈值选择不当最后做轨迹积分验证。如果预测轨迹和真实轨迹在短时间内能对齐说明整套流程是通的剩下的只是调参精度问题。这个顺序能帮你快速定位问题出在表示学习端、稀疏回归端还是数据端。6.3 复现完整流程后的三条心得体会第一别把这个方案当做一个“一键找方程”的黑盒。它更适合做“寻找可解释的候选模型”或者“低维动力学发现”的辅助工具。真实数据里的噪声、缺失点和尺度差异都会显著影响结果预处理往往比网络结构更重要。第二自动编码器和 SINDy 的结合最大价值不是把神经网络替换成了方程而是提供了一个“可解释的中间表示”。也就是说你获得的不只是一堆参数而是“系统在某个坐标系下的动力学规则”。这对于后续控制、预测和物理分析都很有意义。第三如果想让方法更鲁棒可以考虑后续加上符号回归Symbolic Regression、贝叶斯推断或者更精细的验证集选择策略。不过作为入门和论文泛读复现先跑通 AESINDy 这一条通路你对“数据驱动发现方程”这件事的理解会比看十篇论文都深刻。最后分享一个小技巧训练结束之后不要只盯着最终结果看把中间步骤的可视化全部留好——隐变量轨迹、重构轨迹、系数矩阵热力图、稀疏结构变化过程。至少在调试的时候这些图片比任何一张 loss 曲线都更能告诉你模型到底学会了什么。这套可视化思路放到其他神经网络科学计算结合的项目里也一样好用。 动力学重构实战：自动编码器与SINDy联合发现系统方程
最近在整理动力学系统的资料时又重新把“从数据里找方程”这个方向完整过了一遍。很多人第一次听到“动力学重构”会觉得玄乎其实说白了就是手里只有一条或几条时间序列观测数据想反推出系统背后那条微分方程。这不是纯理论问题而是很多真实工程场景的刚需——化学反应中间变量测不全、传感器只有某些状态量的读数、仿真模型太复杂想简化……诸如此类。这篇博文想聊的是论文泛读里一个非常经典的技术路线神经网络自动编码器Autoencoder配合 SINDySparse Identification of Nonlinear Dynamics稀疏非线性动力学辨识一起完成动力学重构并附上可运行的 Python 源码思路。这套方法的论文原型是 2019 年前后 Champion、Lusch、Kutz、Brunton 那篇关于“dynamic autoencoder”的工作核心思想一句话就能概括先用神经网络找到一个隐藏坐标系让系统在这个坐标系下的动力学表达式尽量稀疏、尽量简单再用 SINDy 把这个方程“读”出来。非常适合做时间序列分析、系统辨识、物理仿真、甚至神经网络科学计算方向的朋友参考复现。1. 先搞懂一个问题动力学重构究竟在解决什么1.1 一条轨迹可以有无数个方程所以这是个不适定问题先给你一个直觉层面的例子。假设我在一个二维平面上画了一条圆的轨迹问你这个轨迹对应哪个动力学方程你大概率会回答[ \dot{x} -y,\quad \dot{y} x ]这是最简单的匀速圆周运动。但如果我告诉你其实真实系统是[ \dot{x} -y \epsilon x^3 ]在 (\epsilon) 非常小、并且你只观测了很短一段轨迹时两条方程产生的曲线几乎叠在一起。你根本没法从有限数据里区分到底哪个是“真”的。这就是动力学重构的第一个难点它不是简单地“拟合一条曲线”而是要在无数个能产生相似轨迹的候选方程里挑出一个最合理、最可解释的。靠纯黑箱拟合解决不了这个问题因为拟合出来的参数再多也很难告诉你说“这个系统就是这三个非线性项在起主导作用”。1.2 从“拟合”到“发现”稀疏辨识带来的视角转换传统系统辨识的做法是先假设模型结构比如“我猜这是一个二阶线性系统”然后去估计里面的系数。问题是模型结构一旦猜错后面所有结果都白搭。SINDy 换了个思路先别急着定结构而是给系统准备一个内容足够丰富的候选函数库把常数项、一次项、二次项、三次项、甚至三角函数项都摆进去然后让回归算法在里面挑出真正有用的那一小部分。绝大多数项的系数是 0只有少量非零系数这个“稀疏解”本身就是模型结构。你可以把它理解成侦探破案先列出一堆可疑人员候选函数经过排查之后锁定三五个关键人物非零系数。其他的要么排除要么不太重要。这种思路的好处是你不用再人工猜测“它是不是多项式系统”“是不是含三角函数”只要候选库里包含真实项稀疏回归就有机会把它找出来。1.3 坐标问题方程复杂不复杂取决于你站在哪里看直接 SINDy 有个隐性前提你手上的观测坐标差不多是“好坐标”。什么叫好坐标对经典 Lorenz 系统来说直角坐标系 (x,y,z) 就是好坐标因为它的方程在三次多项式库里只有几个非零项非常稀疏。但现实里传感器量测的坐标往往不是“上帝视角”的好坐标而是各种状态量的线性或非线性混合。一旦坐标系取得不好原本很简单的动力学方程会变得极其臃肿。举个例子还是 Lorenz 系统如果我把观测变量改成 (p x^2 y^2)、(q e^z) 这种非线性组合那么 (p,q) 满足的微分方程会包含大量的高次交叉项SINDy 根本筛不干净。这正是引入自动编码器的原因让网络自己学着把观测数据映射到一个隐空间里并在这个隐空间上要求动力学方程尽量稀疏。换句话说自动编码器负责换坐标系SINDy 负责在换好坐标系后找方程两者是分工协作的关系。2. SINDy 的数学内核稀疏回归怎么“发现”方程2.1 候选函数库与线性化结构SINDy 的精髓在于把非线性动力学问题转化为线性回归问题做法是这样的。设系统状态为 (z(t) \in \mathbb{R}^d)可以是观测变量也可以是隐变量假设动力学可以写成[ \frac{dz}{dt} f(z) ]我们不知道 (f) 的形式但可以假设 (f) 能表示成一组候选基函数的线性组合。把候选基函数堆成一排向量[ \Theta(z) [1,\ z_1,\ z_2,\ \dots,\ z_1^2,\ z_1z_2,\ \dots,\ \sin(z_1),\ \dots] ]于是[ \frac{dz}{dt} \Theta(z),\Xi ]其中 (\Xi) 是一个系数矩阵每一列对应一个状态变量的方程。只要 (\Theta(z)) 里包含了真实的非线性项那么 (\Xi) 就会很稀疏——绝大多数位置是 0只有少数位置非零。这一步就是 SINDy 的“稀疏辨识”雏形。这里特别提醒一点为了能捕捉二次非线性动力学(z) 的候选函数库至少要把多项式阶数开到 2 或 3。Lorenz、Rossler、Lotka-Volterra 这类经典系统在三次多项式库内都能被非常简洁地表示出来。反过来如果你只给线性函数库SINDy 再厉害也识别不出非线性项所以候选库的完整性直接决定了结果上限。2.2 数值求导实操里最容易被忽略的误差源SINDy 的输入听起来很简单左边是导数右边是函数库矩阵。但“导数”从哪来真实传感器数据不会直接给你导数你需要从离散时间序列里数值估计。最常用的中心差分[ \dot{x}(t_k) \approx \frac{x(t_{k1}) - x(t_{k-1})}{2\Delta t} ]这个公式问题非常大如果信号里有高频噪声差分运算会把噪声放大到离谱的程度放大幅度大约是 (1/\Delta t)。假设采样时间 (\Delta t 0.01)噪声标准差是 0.01那么导数的噪声标准差可能直接飙到 1 左右比真实导数本身还大。实操中我一般先做两步预处理采样率要足够高至少保证每个振荡周期有 50 个点以上时序数据先做平滑比如 Savitzky-Golay 滤波或简单卷积平滑再计算导数。如果你用的是干净仿真数据中心差分通常够了。但在真实实验数据上SINDy 的识别效果受导数噪声影响非常大这一点后面遇到的坑多半都出在这里。2.3 稀疏回归与阈值迭代知道了数据矩阵之后剩下的就是解一个稀疏线性回归[ \min_{\Xi} \ \frac{1}{2}|\dot{Z} - \Theta(Z)\Xi|_F^2 \lambda|\Xi|_1 ]这看起来就是一个 Lasso 问题但 SINDy 原论文里更常用的一种做法叫STLSQSequentially Thresholded Least Squares顺序阈值最小二乘流程非常直观先用最小二乘求一个初始系数矩阵(\Xi \Theta^ \dot{Z})给定一个阈值 (\lambda)把绝对值小于 (\lambda) 的系数全部置零用剩下的非零列重新做一次最小二乘重复 2、3 步直到系数矩阵不再变化。阈值 (\lambda) 怎么选没有一个万能值。我的经验是先设一个相对较大的阈值比如 0.5把明显多余的项全部干掉然后逐步调小到 0.1、0.05观察哪些项稳定保持非零。你还可以用“系数稀疏程度 vs 拟合误差”的 Pareto 曲线来判断曲线拐点附近的阈值通常是比较好的选择。2.4 一个 30 行代码的 SINDy 最小实现我先给一个不依赖任何深度学习库的 SINDy 最小实现用 Lorenz 数据直接跑这样后面的自动编码器路线会更容易理解。import numpy as np from scipy.integrate import solve_ivp def lorenz(t, state, sigma10.0, rho28.0, beta8/3): x, y, z state return [sigma*(y-x), x*(rho-z)-y, x*y-beta*z] # 生成 Lorentz 轨迹 dt 0.001 t_span (0, 20) t_eval np.arange(0, 20, dt) sol solve_ivp(lorenz, t_span, [1.0, 1.0, 1.0], t_evalt_eval, rtol1e-9, atol1e-9) X sol.y.T # (N,3) # 中心差分求导 dX np.gradient(X, dt, axis0) def poly_library(Z, order3): n, d Z.shape cols [np.ones(n)] for deg in range(1, order1): # 生成所有次数为 deg 的单项式 if d 1: idx_list [(i,) for i in range(d)] # 用递归或 itertools 生成 d 维单项式索引这里简化为 d 较小场景 from itertools import combinations_with_replacement for comb in combinations_with_replacement(range(d), deg): col np.ones(n) for idx in comb: col * Z[:, idx] cols.append(col) return np.column_stack(cols) Theta poly_library(X) def stlsq(Theta, dX, threshold0.05): Xi np.linalg.lstsq(Theta, dX, rcondNone)[0] for _ in range(5): small np.abs(Xi) threshold Xi[small] 0 for j in range(Xi.shape[1]): if np.any(Xi[:, j] ! 0): idx Xi[:, j] ! 0 Xi[idx, j] np.linalg.lstsq(Theta[:, idx], dX[:, j], rcondNone)[0] return Xi Xi stlsq(Theta, dX, threshold0.05) np.set_printoptions(precision2, suppressTrue) print(Xi)这个代码跑完你会看到系数矩阵在对应位置出现 (10)、(-10)、(28) 等接近真实参数的数值而其他位置基本是 0。这就是 SINDy “发现方程”的过程——非常朴素但非常有效。3. 自动编码器补上关键一环找到合适的坐标系3.1 为什么直接 SINDy 常常不够如果你手上的观测本身就已经是简洁坐标SINDy 确实够用。但现实问题常常是观测维度高比如高速相机拍一个摆动物体一帧就是几百维像素观测值是多个物理量的混合真实“状态”不直接可见变量之间存在强耦合直接在每个观测维度上写方程候选库会爆炸。举个例子如果观测维度是 100 维哪怕只用二次多项式库候选函数数量也有大约 5000 个。每个变量的方程都要在 5000 个基函数里做稀疏回归计算量大不说还极容易把噪声项也学进去。这个时候需要降维。传统 PCA 可以做线性降维但真实动力学系统的低维流形往往是非线性的。自动编码器的优势就在这里——它能学习一个非线性映射把高维观测压缩到低维隐变量 (z)再从 (z) 解码回观测空间。3.2 关键结合点不仅要“重构得好”还要“动力学简单”但是如果只用一个普通自动编码器去做降维得到的隐变量不一定会让动力学方程变简单。普通 AE 的损失函数是[ \mathcal{L}_{recon} |x - D(E(x))|_2^2 ]它只要求解码后能还原原数据至于隐空间里的时序演化长什么样它完全不管。于是 AE 很可能学习到一个扭曲的隐空间重构误差很小但 (z(t)) 的时间演化极其复杂SINDy 在隐空间上直接失效。所以要把 AE 和 SINDy 联合训练总的损失函数变成[ \mathcal{L} |x - D(z)|_2^2 \lambda \left| \frac{dz}{dt} - \Theta(z)\Xi \right|_2^2 \lambda_1 |\Xi|_1 ]第一项是重构损失保证 (z) 能完整表示观测信息第二项是动力学损失要求隐变量满足 SINDy 给出的稀疏方程第三项是稀疏正则迫使 (\Xi) 尽量稀疏。这个联合训练的精妙之处在于网络可以自由变换坐标但变换后的坐标必须满足两个条件——既能重构回原来的观测数据又让隐空间里的方程变得极其稀疏。这在数学上等价于去搜索一个“能稀疏化动力学的坐标变换”。我自己的理解是这个过程有点像把一团乱麻的毛线整理好AE 负责把缠绕的线重新卷到线轴上SINDy 则负责告诉你这个线轴每转一圈线长变化遵循什么规律。两者缺一不可。3.3 完整前向流程整套方法跑起来时数据流是这样的拿到观测轨迹 (x(t))把长轨迹切成固定长度的小片段每个片段输入编码器 (E)得到隐变量轨迹 (z(t))用中心差分从 (z(t)) 算出隐变量的导数 (\dot{z}(t))构造候选函数库 (\Theta(z))用 SINDy 或者直接用梯度优化求解稀疏系数 (\Xi)用解码器 (D) 把 (z(t)) 映射回观测空间计算重构误差把重构误差和动力学误差合并反向传播更新编码器、解码器和 (\Xi)。训练完成之后你不只是得到一个能重构数据的自编码器更重要的是你拿到了一组明确的常微分方程(\dot{z} \Theta(z)\Xi)。以后你想预测系统未来行为不需要再跑神经网络直接积分这组方程就行——这才是真正的“发现方程”。4. PyTorch 源码拆解从数据生成到联合训练4.1 准备训练数据多轨迹、切片、归一化直接拿一整条长轨迹训练AE 很容易过拟合到初始条件的特殊路径上。我一般会准备 20 条以上不同初始条件的轨迹然后按固定窗口切成小样本。以 Lorenz 为例import torch import numpy as np from scipy.integrate import solve_ivp def generate_lorenz_data(n_traj20, traj_len5000, dt0.002, window128): data [] for _ in range(n_traj): y0 np.random.uniform(-10, 10, size3) t_eval np.arange(0, traj_len*dt, dt) sol solve_ivp(lorenz, (0, traj_len*dt), y0, t_evalt_eval, rtol1e-9, atol1e-9) data.append(sol.y.T) data np.stack(data, axis0) # (n_traj, traj_len, 3) # 按窗口切段 windows [] for seq in data: for i in range(0, traj_len - window, window//2): windows.append(seq[i:iwindow]) windows np.stack(windows, axis0) # 归一化 mean windows.reshape(-1, 3).mean(axis0) std windows.reshape(-1, 3).std(axis0) windows (windows - mean) / std X torch.tensor(windows, dtypetorch.float32) # (N, window, 3) return X注意归一化非常关键。Lorenz 的三个变量尺度差异不大但如果是真实系统变量间量纲差异可能很大不归一化的话高次多项式项会数值溢出或者某个变量主导整个损失。4.2 搭自编码器MLP 编码器/解码器网络结构不需要很深我用的是三到四层全连接import torch.nn as nn class Encoder(nn.Module): def __init__(self, obs_dim3, hidden_dim64, latent_dim3): super().__init__() self.net nn.Sequential( nn.Linear(obs_dim, hidden_dim), nn.Tanh(), nn.Linear(hidden_dim, hidden_dim), nn.Tanh(), nn.Linear(hidden_dim, latent_dim), ) def forward(self, x): return self.net(x) class Decoder(nn.Module): def __init__(self, latent_dim3, hidden_dim64, obs_dim3): super().__init__() self.net nn.Sequential( nn.Linear(latent_dim, hidden_dim), nn.Tanh(), nn.Linear(hidden_dim, hidden_dim), nn.Tanh(), nn.Linear(hidden_dim, obs_dim), ) def forward(self, z): return self.net(z)这里一个实用小技巧数据是按时间窗口过来的所以前向传播时给每个时刻的状态单独过 MLP然后把整个窗口拼接起来。在 PyTorch 里直接用Encoder(x)会默认作用在最后一维上所以形状(N, window, obs_dim)可以直接传进 MLP因为nn.Linear会作用在最后一维。4.3 构造 SINDy 候选库与导数计算隐变量 (z) 的候选库构造方式和前面原生 SINDy 一样只是这里要支持 PyTorch 的自动求导所以用张量操作实现from itertools import combinations_with_replacement def poly_library_torch(z, order3): # z: (N, T, d) n, T, d z.shape z z.reshape(-1, d) # 合并 N,T columns [torch.ones(n*T, 1, devicez.device)] for deg in range(1, order1): for comb in combinations_with_replacement(range(d), deg): col torch.ones(n*T, 1, devicez.device) for idx in comb: col col * z[:, idx:idx1] columns.append(col) return torch.cat(columns, dim1) # (n*T, n_terms)导数计算可以直接用torch.gradientdef compute_dz_dt(z, dt): # z: (N, T, d) dz torch.gradient(z, spacingdt, dim1)[0] return dz如果你的数据噪声不大torch.gradient就行如果噪声明显我建议先对 z 轨迹做平滑再算导数。4.4 联合训练损失函数与稀疏系数更新这里有个实现细节要特别提醒(z) 的导数一旦通过torch.gradient算出来它是依赖于z的计算图的可以正常反向传播。但 (\Xi) 的更新有两种常见方式直接把 (\Xi) 当成可训练参数用 Adam 优化器更新每隔若干轮对 (\Xi) 做一次硬阈值收缩其实也就是 STLSQ 的思路。我比较推荐混合方案Adam 负责更新编码器和解码器(\Xi) 则每隔 10 轮做一次阈值收缩。因为如果让 Adam 直接优化 (\Xi) 并加 L1 正则系数很容易全部缩到 0或者某个稀疏结构很难被“钉”住。class DynamicAutoencoder(nn.Module): def __init__(self, obs_dim3, latent_dim3, hidden_dim64, order3): super().__init__() self.encoder Encoder(obs_dim, hidden_dim, latent_dim) self.decoder Decoder(latent_dim, hidden_dim, obs_dim) self.order order # Thêta(z) 的项数 n_terms 0 for deg in range(order1): n_terms len(list(combinations_with_replacement(range(latent_dim), deg))) self.Xi nn.Parameter(0.1 * torch.randn(n_terms, latent_dim)) def forward(self, x): z self.encoder(x) x_hat self.decoder(z) return z, x_hat def train(model, X, dt0.002, epochs300, lr1e-3, lambda_dyn10.0, threshold0.05): optimizer torch.optim.Adam(model.parameters(), lrlr) for epoch in range(epochs): optimizer.zero_grad() z, x_hat model(X) # X: (N, T, obs_dim) recon_loss torch.mean((x_hat - X)**2) dz torch.gradient(z, spacingdt, dim1)[0] Theta poly_library_torch(z, model.order) # (N*T, n_terms) Xi model.Xi # 动力学损失: || dz - Theta Xi ||^2 dyn_loss torch.mean((dz.reshape(-1, z.shape[-1]) - Theta Xi)**2) sparse_reg torch.mean(torch.abs(Xi)) loss recon_loss lambda_dyn * dyn_loss 1e-3 * sparse_reg loss.backward() optimizer.step() # 每 10 轮做一次硬阈值收缩 if epoch % 10 0: with torch.no_grad(): small torch.abs(model.Xi) threshold model.Xi[small] 0.0 if epoch % 50 0: print(fepoch {epoch}, recon{recon_loss.item():.4f}, dyn{dyn_loss.item():.4f})这里有几个参数调起来很讲究。lambda_dyn是动力学损失在总损失里的权重我试过 0.1、1、10、100结论是太小的话 AE 学到一个重构很好但动力学极其混乱的隐空间太大的话 AE 的重构能力会被牺牲隐变量和原始观测之间的信息无法完全保持。一般我会先在 PCA 或普通 AE 预训练基础上从lambda_dyn 1开始试再看重构损失和动力学损失的平衡。4.5 训练后验证隐空间积分 解码回观测空间当你把模型训练完拿到了稀疏的 (\Xi)怎么确认它真的发现了系统方程最直观的验证方式是从某个初始隐状态出发用识别出的 ODE 做数值积分然后把整条轨迹解码回观测空间再和真实观测轨迹做对比。from scipy.integrate import solve_ivp def identify_rhs(z, Xi, order3): # 返回 Theta(z) Xi Theta_z poly_library_torch(torch.tensor(z.reshape(1, 1, -1), dtypetorch.float32), order).detach().numpy() return (Theta_z Xi.detach().numpy()).flatten() def predict_trajectory(model, z0, t_eval): Xi_np model.Xi.detach().numpy() sol solve_ivp(identify_rhs, (t_eval[0], t_eval[-1]), z0, t_evalt_eval, args(Xi_np, model.order), rtol1e-8, atol1e-8) z_pred sol.y.T z_tensor torch.tensor(z_pred, dtypetorch.float32) x_pred model.decoder(z_tensor).detach().numpy() return x_pred, z_pred如果识别正确预测出来的相位图会和真实数据的吸引子形状高度重合。如果只是重构损失很小、但预测轨迹很快漂移发散说明动力学损失没学好隐空间的微分方程是错的。5. 调参与踩坑记录复现过程中最值得注意的细节5.1 AE 发现的方程不一定是原始方程这是很多人第一次复现时最容易困惑的地方。你用 Lorenz 数据喂给 AESINDy最后得到的隐变量方程往往不是标准形式[ \dot{x} 10(y-x),\quad \dot{y}x(28-z)-y,\quad \dot{z}xy-\frac{8}{3}z ]而是某个可逆线性变换或者非线性变换之后的形式。比如隐变量可能是 (z_1 x y)、(z_2 x - y)、(z_3 z) 这样的组合方程系数完全不同但动力学性质等价。所以验证结果的时候不要执着于“系数必须等于真实系统”而要看两点隐空间 ODE 预测出的轨迹解码回观测空间后是否逼近真实轨迹识别出的方程是否稀疏非零项数量是否远小于候选库规模。只要轨迹对得上、方程又稀疏那这套方法就算成功了。坐标本身并不唯一允许任意可逆变换是这一类“表示学习动力学辨识”方案的理论基础。5.2 损失权重的“跷跷板”重构损失和动力学损失之间存在天然博弈。如果你把lambda_dyn调很大AE 会倾向于找到一条“最好拟合成稀疏方程”的隐空间轨迹但代价是 z 可能丢失一些重构所需的信息导致从隐变量解码回观测的时候面目全非。反之lambda_dyn太小动力学方程会变得特别复杂SINDy 很多非零系数都撂在那里根本无法解释。我建议的顺序是先单独用重构损失预训练 AE得到一组正常的编码器/解码器权重再加载预训练权重从lambda_dyn 1开始联合训练观察两个损失的比例如果重构损失开始飙升就降低lambda_dyn如果动力学损失不降就增大它。还有一种比较省心的做法动力学损失不直接对隐变量的导数做拟合而是对“积分之后的轨迹”做拟合——也就是说把整条 z 轨迹用识别出的 ODE 积分出来再与编码器给出的 z 轨迹做对比。这样能避免数值微分带来的高频噪声影响但计算量更大实现也更复杂。新手阶段先用导数拟合就够了。5.3 数据长度与状态空间覆盖SINDy 和 AE 都很吃数据质量。如果你的训练数据只是 Lorenz 吸引子上很窄的一小段系统从来没有经历过大幅振荡那么识别出的方程很可能只在小范围适用。换个初始条件预测几步误差瞬间爆炸。我一般会刻意让训练轨迹覆盖更大的状态空间初始条件在 ([-15, 15]^3) 区间里随机采样然后把多条轨迹混合训练。这比单纯把一条轨迹拉长更有效因为 SINDy 需要看到系统在不同区域的“反应”才能真正判断哪些非线性项是必需的。5.4 阈值、候选库尺度和数值稳定性候选函数库如果包含三阶甚至四阶多项式项隐变量稍微超出训练范围时库函数数值会迅速膨胀。比如 (z_1^3)当 (|z_1|10) 时是 1000 量级跟常数项 1 完全不在一个量级。所以训练时归一化很重要最好让隐变量基本落在 ([-2, 2]) 区间内。另外阈值不是越大越好阈值设置效果问题0.5 以上方程非常稀疏可能把真实非零项也删掉识别失效0.1 ~ 0.3多数情况下效果不错需要联合损失配合避免遗漏小系数项0.01 以下非零项很多几乎不稀疏难以解释我是从 0.1 起步看识别出的方程里非零项数量如果太多就升高阈值如果出现动力学损失完全不收敛就降低阈值。这个调参过程不复杂但必须结合轨迹验证一起看单看损失曲线很容易被骗。6. 常见问题速查表与排查思路6.1 常见问题速查表下面是我在调这套方法时经常遇到的问题整理成一张速查表方便你对照排查问题现象可能原因解决思路重构损失很低但预测轨迹发散发散动力学损失权重不够隐空间方程没学准增大lambda_dyn增加联合训练轮数动力学损失很低但重构损失很高编码器把信息丢了AE 退化成恒等映射附近的变形降低lambda_dyn增大隐变量维度增加编码器容量系数矩阵几乎全零稀疏阈值太大把真实项砍掉了调低阈值或者先增大预训练轮数系数矩阵一堆非零小系数阈值太小稀疏性不足调高阈值或增加 L1 正则权重隐变量轨迹高频抖动数值微分放大了编码器输出的微小波动对编码器输出做平滑或减小采样时间步dt训练到最后损失还在震荡学习率太高或 AE 结构和 SINDy 联合训练不稳定降低学习率加入学习率衰减增大 batch 或序列长度换了真实数据后完全失效真实数据噪声大候选库没有覆盖系统真实结构先做数据降噪换更高阶库或增加三角函数项6.2 排查顺序从损失到系数再到轨迹如果你复现时出了诡异问题别急着改代码。我一般是按下面这个顺序拆解先看重构损失。如果重构损失不降说明 AE 本身就没训练好后面的动力学都白搭再看动力学损失。如果重构很好但动力学损失下不去大概率是隐空间维度不够或者候选函数库缺了关键项打印 (\Xi) 矩阵看非零项是否稳定。如果每次训练非零项位置都不一样说明数据覆盖不足或者阈值选择不当最后做轨迹积分验证。如果预测轨迹和真实轨迹在短时间内能对齐说明整套流程是通的剩下的只是调参精度问题。这个顺序能帮你快速定位问题出在表示学习端、稀疏回归端还是数据端。6.3 复现完整流程后的三条心得体会第一别把这个方案当做一个“一键找方程”的黑盒。它更适合做“寻找可解释的候选模型”或者“低维动力学发现”的辅助工具。真实数据里的噪声、缺失点和尺度差异都会显著影响结果预处理往往比网络结构更重要。第二自动编码器和 SINDy 的结合最大价值不是把神经网络替换成了方程而是提供了一个“可解释的中间表示”。也就是说你获得的不只是一堆参数而是“系统在某个坐标系下的动力学规则”。这对于后续控制、预测和物理分析都很有意义。第三如果想让方法更鲁棒可以考虑后续加上符号回归Symbolic Regression、贝叶斯推断或者更精细的验证集选择策略。不过作为入门和论文泛读复现先跑通 AESINDy 这一条通路你对“数据驱动发现方程”这件事的理解会比看十篇论文都深刻。最后分享一个小技巧训练结束之后不要只盯着最终结果看把中间步骤的可视化全部留好——隐变量轨迹、重构轨迹、系数矩阵热力图、稀疏结构变化过程。至少在调试的时候这些图片比任何一张 loss 曲线都更能告诉你模型到底学会了什么。这套可视化思路放到其他神经网络科学计算结合的项目里也一样好用。