PACE：一种基于Pareto-AICc-CV共识估计方法的自动化阈值选择工具，用于稀疏识别非线性动力系统

《Communications in Nonlinear Science and Numerical Simulation》：PACE: Pareto–AICc–CV consensus estimator for automated threshold selection in sparse identification of nonlinear dynamical systems

　　摘要：非线性动力学的稀疏识别（SINDy）可以从时间序列数据中发现控制方程，但其核心的顺序阈值最小二乘（STLSQ）算法严重依赖于一个稀疏阈值λ，该阈值通常通过手动调整或网格搜索来确定，这限制了可重复性。我们引入了PACE（Pareto–AICc–CV共识估计器），这是一个自动化框架，可以自动选择阈值参数。PACE遍历对数间隔的网格，并计算两种方法上不同的估计值：一种是从残差-稀疏权衡中得到的Pareto肘点估计λP，另一种是考虑时间自相关的h块交叉验证（CV）估计λCV。这些估计定义了一个候选区域，在该区域内，修正后的赤池信息量准则（AICc）选择最终阈值λ*；CV检验验证了预测竞争力，而AICc宽度诊断则赋予校准的置信度标签（高、中等、低）。在200次基准试验中——包括四个系统（Lorenz-63、Van der Pol、Lotka–Volterra、R?ssler）和五种噪声水平（1-20%）以及十个随机种子——PACE在112次试验中实现了精确的支持恢复（56.0%）。根据这一严格标准，PACE与最强的基线方法（交叉验证，55.0%）相匹配，并且在任何噪声水平下都不会低于它，并且恢复的模型数量大约是任何基于信息量准则的方法的两倍（最多为33.5%），而在高噪声情况下，其他方法几乎无法恢复模型，而PACE在80次试验中仍恢复了30次。此外，PACE还提供了一个基线方法所缺乏的校准置信度标签：高置信度标记的预测准确率为77%，而低置信度标记的准确率为11%，而Pareto-CV的一致性选择实现了100%的正确恢复率。

引言：直接从实验或模拟的时间序列数据中提取可解释的数学模型是现代非线性动力学和系统识别的核心目标[1][2]。经典的灰箱和黑箱方法要么需要对模型结构有很强的先验知识，要么为了预测准确性而牺牲可解释性。数据驱动的稀疏回归提供了一个原则性的中间路径：通过假设控制方程只涉及用户指定候选库中的一小部分项，可以将发现控制方程的问题转化为一个稀疏选择问题。稀疏识别非线性动力学（SINDy）框架由Brunton、Proctor和Kutz提出[1]，通过顺序阈值最小二乘（STLSQ）算法实现了这一理念。从沿测量轨迹评估的候选函数库Θ(x)开始，STLSQ迭代解决最小二乘问题，同时移除幅度低于阈值参数λ的系数。由此产生的稀疏系数矩阵Ξ编码了发现的控制方程。自SINDy问世以来，它已成功应用于广泛的领域，包括流体动力学[3]、化学反应网络[4]、生物振荡器[5]、模型预测控制[6][7]。尽管取得了这些成功，但一个根本的实际挑战仍未解决：阈值参数λ的选择。如果λ太小，噪声驱动的虚假项会通过阈值筛选，导致恢复的模型过拟合；相反，如果λ太大，则会移除物理上有意义的项，得到不完整或欠定的模型。在测量噪声存在的情况下，以及真实系数动态范围较广的系统中，这种敏感性尤为严重。在原始的SINDy工作及后续的大部分文献中，λ通常是通过手动检查系数谱或遍历值网格并选择产生物理上合理的方程来确定的[1][8]。这些方法具有主观性，对于大型库来说计算效率低下，并且难以在不同的噪声实现中重复。为了应对这一挑战，已经提出了几种原则性的方法，包括信息量准则模型选择[8]、集成和稳定性选择方法[9][10]，以及像SR3这样的放宽正则化框架[11]，它同时概括了STLSQ和LASSO[12]。这些方法在第2节中进行了详细讨论；然而，每种方法要么引入了次要的超参数，要么没有产生单一的确定性阈值，要么依赖于在强测量噪声下会退化的残差统计。尽管有这些进展，但目前没有任何方法能够同时提供：(i)自动化的、数据驱动的阈值参数λ的确定；(ii)多种互补估计器的原则性组合，以降低选择次优解的风险；(iii)基于信息量准则的主要仲裁步骤，在明确定义的候选区域内进行选择，而不是在整个网格上；(iv)校准的置信度标签，向最终用户传达每个恢复模型的可靠性。为了应对这些挑战，我们提出了PACE（Pareto–AICc共识估计器），这是一个基于标准STLSQ算法的六阶段自动化λ选择流程。PACE遍历对数间隔的候选阈值值网格，并计算两种方法上不同的λ估计值：一种来自残差-稀疏权衡曲线的Pareto肘点（λP），另一种使用h块交叉验证（λCV）获得，后者考虑了动态系统数据中的时间自相关性。PACE不是直接比较两个候选值的稀疏模式，而是定义一个候选区域[min(λP, λCV), max(λP, λCV)]，并在此区域内最小化修正后的赤池信息量准则（AICc）以获得共识估计λ*。然后通过交叉验证检验验证AICc选择的解决方案的预测竞争力；如果检验失败，算法会回到λP。随后通过AICc宽度诊断分配置信度。当Pareto和交叉验证估计器收敛到相同的稀疏模式时（一致决定），恢复的模型在所有测试的基准系统上达到100%的准确率。我们在200次基准试验中对PACE进行了验证，这些试验涵盖了四个典型的非线性动力学系统——混沌的Lorenz-63系统[14]、极限环Van der Pol振荡器[15]、保守的Lotka–Volterra方程[16]和R?ssler吸引子[17]——涵盖了五种噪声水平和每个配置的十个独立随机种子。本文的其余部分组织如下：第2节回顾了SINDy框架和相关工作；第3节阐述了问题；第4节描述了Savitzky–Golay导数估计方法；第5节详细介绍了PACE流程；第6节概述了基准设计；第7节展示了结果；第8节讨论了关键发现、设计决策和局限性；第9节总结了本文。

SINDy和STLSQ：考虑一个连续时间的自治动力系统˙??=???(??)，其状态???(??)=[??1?(??),…,?????(??)]?∈???，假设??:???→???是平滑且在已知函数基中稀疏的；目标是从采样的、带噪声的x(t)观测数据中识别f。SINDy[1]将这个问题视为一个稀疏回归问题，将时间导数表示为从候选库???(??)∈???×??中提取的列的稀疏线性组合：˙??=???(??)???，其中??∈???×??是轨迹矩阵，˙??是其数值时间导数，??∈???×??是系数矩阵。

问题表述：表1总结了整个文中使用的主要符号。基于第2.1节中的SINDy公式，我们观察到N个均匀采样的状态??∈???×??，这些状态受到加性高斯噪声的干扰，并从候选基函数θj(x)（例如最高为r次的单项式）构建库???(??)∈???×??。将数值估计的时间导数矩阵写作˙??∈???×??，每个状态按列建模为˙????≈???????，其中????∈???是状态i的稀疏系数向量。

Savitzky–Golay导数估计：对噪声时间序列数据进行数值微分是一个病态问题[35]，这直接限制了SINDY的恢复质量[34][36]。Savitzky–Golay（SG）滤波器[37]提供了一种实用的方法，它通过在一个滑动窗口内的w个样本上拟合一个q次的多项式来最小二乘地获得导数，从而可以分析地得到导数。窗口宽度控制了偏差-方差权衡：较大的窗口可以抑制噪声，但在快速变化区域增加了平滑偏差。

PACE算法：PACE包括六个连续阶段：(1) lambda网格构建和遍历；(2) Pareto肘点候选估计；(3) h块交叉验证候选估计；(4) 候选区域构建和区域内AICc仲裁；(5) 交叉验证检验；(6) 通过AICc宽度分配置信度。整个方法如图1所示。各个阶段在以下小节中详细描述。在介绍这些阶段之前，图1给出了概念性的说明。

动态系统：选择了四个典型的基准系统来代表不同的动态制度：一个混沌吸引子、一个稳定的极限环、一个保守的振荡器和一个螺旋吸引子。

Lorenz-63（混沌吸引子）[14]：˙??=???(?????)，˙??=???(?????)???，˙??=???????????，参数为??=10，??=28，??=8/3。真实的稀疏性为??★=[2,3,2]（在每个状态有20列的三次多项式库中有7个非零项）。混沌吸引子对λ的选择是一个严格的测试，因为它对参数的敏感依赖性很强。

整体性能：表4总结了PACE在所有200次基准试验中的性能。PACE在112次试验中实现了精确的支持恢复（56.0%）。Van der Pol振荡器在所有噪声水平和随机种子下都获得了完美的50/50（100%）成绩，Lotka–Volterra系统获得了40/50（80.0%）。Lorenz-63和R?ssler系统分别获得了12/50（24.0%）和10/50（20.0%），这反映了第7.7节中讨论的噪声下限敏感性。

一致选择和双重估计器共识：在54次一致决定中100%的通过率是本工作中最有力的发现之一。这表明，当两种方法上不同的标准——一个是几何的（Pareto残差-稀疏权衡）和一个是预测的（h块CV）——在没有AICc仲裁的情况下达成一致时，恢复的模型在测试的基准系统上始终是正确的。这一结果验证了双重估计器架构的合理性，再加上校准的置信度标签（第7.4节），鼓励将一致选择视为一种有效的模型识别方法。

结论：我们提出了PACE（Pareto–AICc–CV共识估计器），这是一个用于STLSQ算法中阈值选择的六阶段自动化流程，用于非线性动力系统的稀疏识别。PACE遍历对数间隔的候选阈值值网格，并通过结合两种方法上不同的估计器——Pareto肘点结构估计器（λP）和h块交叉验证预测估计器（λCV）——来定义一个候选区域，从而解决了手动选择λ的关键问题。