目录
Paper Reading 是从个人角度进行的一些总结分享，受到个人关注点的侧重和实力所限，可能有理解不到位的地方。具体的细节还需要以原文的内容为准，博客中的图表若未另外说明则均来自原文。
| 论文概况 | 详细 |
|---|---|
| 标题 | 《PS-Tree: A piecewise symbolic regression tree》 |
| 作者 | Hengzhe Zhanga, Aimin Zhoua, Hong Qiana, Hu Zhang |
| 发表期刊 | 《Swarm and Evolutionary Computation》 |
| 发表年份 | 2022 |
| 期刊等级 | 中科院 SCI 期刊分区(2022年12月最新升级版）1区 |
| 论文代码 | https://github.com/zhenlingcn/PS-Tree |
作者单位：
- School of Computer Science and Technology, Shanghai Institute of AI for Education, East China Normal University;
- Science and Technology on Complex System Control and Intelligent Agent Cooperation Laboratory, Beijing Electro-Mechanical Engineering Institute
研究动机
目前已经有各种各样的回归方法被提出，从可解释性的角度可以将这些方法分为黑盒方法和白盒方法。黑盒方法的一个典型例子是深度神经网络，虽然它们被广泛用于各种问题，但是也因其低解释性而被诟病。基于树的方法具有较高的准确性和可解释性，其中回归树的基本思想是将特征空间划分为几个子区域，然后使用几个简单的模型在这些子区域中拟合数据，例如 CART。与神经网络相比，回归树更容易实现和理解，但是它们本质上是分段常数函数。所以传统的回归树可能很难描绘出自变量和因变量之间的平滑关系，因此开发一种可解释的算法来解决非线性分段回归问题仍然是一个关键问题。
为了改进 CART 的缺点，一些新方法使用复杂的模型来取代 DT 的叶子节点中的简单模型，例如分段线性回归树(PL-Tree)。PL-Tree 采用线性模型而不是常数作为回归树的叶节点模型，但当训练数据集中存在复杂的特征时，PL-Tree 在学习非线性关系时可能会遇到一些困难。下图展示了 Regression Tree、PL-Tree、GP 和本文算法拟合同一条曲线的例子：
以 PREM 模型为例，PREM 是一个复杂的非线性分段模型，描述了地球的性质。直接将这四个模型应用于采样数据，并将预测结果绘制在下图中。可以看到 CART 可以描述数据的突变，可以构建分段模型，但是 CART 在每一段都产生了一条崎岖的线。GP 只能产生一个模拟 PREM 突变的斜曲线，但更容易产生平滑的模型。PL-Tree 不仅能分段，还能在每个子区域构建线性模型，但是可能因为构建的都是线性模型导致性能比 CART 还差。
根据上面的分析以上实验结果，本文作者计划结合 CART 和 GP 的优点，构造出一个分段的非线性模型来解决这些问题。
文章贡献
综合 CART 和 GP 方法的优点，本文提出了一种新的回归方法——分段符号回归树(PS-Tree)。基本思想是通过CART将特征空间划分为若干个子区域，然后使用 GP 和岭回归为每个子区域构建一个简单的回归模型。在模型训练算法方面，使用分类树动态学习每个分区的最合适的数据分配方案，并进化出一组 GP 个体来表达非线性特征，在所有区域构建局部岭回归模型。GP 的特征构建被转化为一个多目标优化问题，从而获得了用于所有子区域的一组重要的特征。由于初始空间分区可能不正确，因此算法还部署了一种动态调整分区方案的自适应方法。
分段符号回归树
PS-Tree 算法包括两个步骤，第一步是通过决策树将特征空间划分为子区域，算法采用传统的决策树归纳法构造。第二步是为每个子区域构造一个符号回归函数，主要使用多目标优化算法来进化多个新特征，并构造几个符号回归量。
以房价预测问题为例，首先根据犯罪率将数据集分为低风险组和高风险组。接着利用多目标优化技术得到多个新的非线性特征，并用这些特征在每个子区域建立线性模型。预测阶段时可以根据新数据的犯罪率，选择最合适的决策树局部模型，然后基于非线性模型预测新的房价。
个体表示
PS-Tree 树由两部分组成，上面的部分是一棵 DT，用于将特征空间划分为几个子区域，下面的树包含许多非线性模型，用于学习特定区域的模式。对于上面的部分，在算法开始时会随机或基于某些策略生成一些 DT，在进化算法的每代结束时，算法将基于合适的分区方案重建决策树。对于下半部分，每个分支都建立了一个基于 GP 个体的非线性模型，本文的 GP 通过最小二乘(或梯度下降方法)对每个 GP 模型的系数进行拟合。
此处的 GP 模型使用进化特征合成(EFS)的思想，也就是在种群中的每个 GP 个体代表一个特征，种群中的所有特征共同作用形成一个回归模型。不同的子区域使用相同的新的构造特征的原因是，虽然不同的子区域的数据模式不同，但重用相同的高阶特征集不见得是不可行的。为了适配这个设计，GP 个体的适应度评估使用的是特征重要度而不是训练损失。
特征分区初始化
特征空间分区过程主要包含两个阶段，首先在进化过程之前进行分区初始化，接着在进化过程中进行分区重构。在进化开始前需要先构造一个决策树，将特征空间(数据集)划分为 m(超参数) 个子区域。本文使用 CART 来构造决策树，并使用贪婪分割策略根据 MSE 将特征空间划分为几个互不关联的子区域。
也可以使用 K-means 或 Bayesian Gaussian mixture 模型来生成初始 DT，基本思想是将训练点聚类为几个簇再分别构建 DT。
自适应决策树重建
提供给 PS-Tree 的分区数可能是错误的，或者上文的初始化方法不能对空间进行正确划分，因此需要一种自适应决策树重构方法进行动态调整。算法在每轮进化后都得到了一组局部模型 {f1, …fi…, fm}，此时将每个数据项 xi 对应的分类器 fi 作为标签，可以得到一个用于选择最适合的局部模型的“分类数据集”。通过对此数据集训练一棵分类树，就可以通过它将每个新数据分配到其最合适的局部模型。
例如将数据分为 4 个特征分区，通过训练得到了 (1, 2, 3, 4) 4 个局部模型。如下图所示，此时将原始数据集和 4 个局部模型(1, 2, 3, 4) 作为标签，可以训练出如下图所示的分类树，用于选择最优的局部模型。本文使用 CART 算法，以 Gini 系数作为分裂准则来生成这棵分类树，当然可以选择其他分类算法。
为了降低计算开销和提升可解释性，需要限制分类树的高度来避免生成过多的局部模型，这可能会导致不是所有的数据都能被分配到最适合的局部模型。同时也有可能有某个数据符合多个叶节点的条件，适合使用多个局部模型，因此设置一种软分配策略比较合适。本文将每个数据点 x 的分区分配结果定义为一个 m 维向量，每个元素表示对每个子区域的隶属度，向量的和为 1。在每个子区域构建局部模型时，也将该隶属度向量作为损失函数中的权重向量，对局部模型进行优化。
例如对于下面的图中，第二个分区中的一部分数据点可能在第一个分区中的局部模型中产生更好的预测结果。因此算法将更新标签然后重建一个新的分类树，例如此时 CART 可能会将第一个分区决策边界向右移动一点。
如果不断地合并这两个分区，也就是只可以使用一个模型来拟合这两个区域中的所有数据点，两个分区最终可能合并成一个分区。
多目标训练
为了节省计算资源，PS-Tree 的所有局部模型使用相同的特征集来构建所有局部模型。每个局部模型可以提供给每个 GP 个体一个特征重要性分数，由于很难确定一个权重向量来将所有这些分数汇总为一个标量，因此算法保持分数向量不变，并将 GP 优化问题表示为一个多目标优化问题。
假设当前特征 Φ 在第 k 个子区域得到局部线性模型 fk，如果使用 MSE 作为损失函数，则可以将优化目标表述为下面的公式。其中 pk 是一个权重向量，表示每个数据点属于分区 k 的概率。这个目标表明进化的目标是得到一组特征 Φ，它们在所有子区域 1, …, m 上都能带来好的效果。
如果使用 Wrapper 特征构建方法，虽然可以直接根据损失值来优化这个目标函数，但是这样会带来很大的计算开销。为了避免这个问题，此处采用 Embedded 特征构建方法，使用所有 GP 个体构建单个模型。使用 Embedded 方法无法为每个 GP 个体计算上述的目标函数，因此改为使用每个 GP 个体的特征重要性来代替训练损失作为优化目标。
新的优化目标来源于特征的重要性排序，也就是如果一个特征在模型中很重要，对该特征的微小变化也会导致预测值的显著变化，反之亦然。因此将特征 Φi 的特征重要性定义为特征 Φi 对预测值 fk(Φ) 的贡献 |∂fk/∂Φi|，就可以得到如下面公式所示的真正的目标函数。
此处可以使用多任务优化方法找到一组在所有分区上表现良好的特征，也可以找到一组在不同分区上表现不同的特征。更具体地说，通过采用多目标进化算法(MOEA)，可以允许在除一个分区外的所有分区上都表现不佳的特征被保存。
PS-Tree 算法流程
PS-Tree 的算法流程如下伪代码所示，本文的多目标遗传算法与一般遗传算法一样，主要包括初始化、生成、评价和选择四个步骤。在算法的最后还有一个可选的剪枝步骤，以此来进一步简化得到的模型。
伪代码中的一些符号及其含义如下表所示：
| 符号 | 说明 |
|---|---|
| rand() | 来自统一分布 U(0,1) 的随机数 |
| pc | 交叉概率 |
| pm | 突变概率 |
| max_gen | 最大代数 |
算法中的主要步骤如下表所示：
| 步骤 | 伪代码行数 | 说明 |
|---|---|---|
| 初始化 | 2-6 | 随机初始化一组 s 特征 Φ={Φ1, Φ2, …, Φs}，并将随机初始化的模型 f 设置为最好的模型(ffinal← f)。 |
| 生成 | 11-15 | 通过对父特征 Φold 进行交叉和突变生成一组新的候选特征 Φ。 |
| 评估 | 4–5、17–18 | 特征 Φ 在训练数据 Nk 的每个子集上构建新的局部模型，获得局部模型的系数后，从系数向量中提取值来组成每个特征的适应度向量。 |
| 选择 | 22 | 从当前种群 Φ 和后代种群 Φold 中选择 s 个有价值的特征作为下一代的父特征。 |
| 修剪 | 24 | 使用 Lasso 算法对构建的特征进行修剪，以提高最终模型最终的可解释性，即自动将一些无需使用的特征的系数设置为零。 |
注意即使将优化问题表述为 MOEA，算法也应该只返回一个最终模型来预测那些看不见的数据。因此在评估步骤中需要计算每个局部模型的交叉验证得分，然后将所有子区域的所有交叉验证得分相加为 CV(f)，权重值等于该子区域中的数据数量。最后如果 f 的总分优于目前最佳模型的总分，也就是 CV(f) < CV(ffinal)，则我们将最佳模型 ffinal 替换为当前的模型 f。
生成算子
与其他 GP 算法一样，交叉算子通过随机交换两个旧个体的两个子树生成两个新的个体，突变算子通过随机生成的子树替换旧个体的一部分来生成一个新的个体。
由于本文的方法使用最小二乘法来确定每个特征的系数，多重共线性(指线性回归模型中的解释变量之间由于存在精确相关关系或高度相关关系，使模型估计失真或难以估计准确)将影响特征重要性的正确分配，也就是说生成的特征最好是不相关的。
为了解决这一问题，本文设计了一种预选策略，拒绝那些与其任一父个体的相关性绝对值大于 0.95 的新生成个体，并强制其重新生成。为了节省计算资源，预选只使用训练数据中的前二十个数据点 {x1，…，x20} 来计算每个个体的表示向量 {Φi(x1)，…，(x20)}，和父个体的表示向量来计算相关系数。为了生成更多样化的个体，算法会记录所有曾出现过的的表示向量。
适应度评估
在评估阶段将根据每个特征在多个局部回归模型中的重要性，为其分配几个适应度值。因为岭回归模型易于解释，在 PS-Tree 中选择岭回归模型作为局部回归模型。由于建立这个模型是基于许多非线性特征，因此即使在每个局部区域都构建了岭回归模型，但它仍然是一个非线性模型。
令 Φ∗= Φ∪Φold 表示一个包含新旧特征的特征集，则可以将岭回归模型的损失函数表示为下面的式子。
式子中的符号和含义如下表所示，式子使用了 Tikhonov 正则化来鼓励简约的系数。由于很难为给定的问题确定最佳的 α 值，PS-Tree 将从一组预定义值 α∈[1,10 ^4] 中选择最佳的 α 值。将有 50 个 α 值在对数尺度上均匀分布在 [1,10^4]，这样就可以通过交叉验证法来选择最佳的 α 值。
| 符号 | 说明 |
|---|---|
| w | 线性模型的系数向量 |
| α | 正则化项的系数 |
因为在线性模型 fk 中，fk 关于 Φi 的偏导数等价于该特征的系数(|∂fk/∂Φi| = |wi|)，因此可以直接将系数 wi 的绝对值作为各个特征 Φi 的适应度。PS-Tree 有 m 个局部模型，每个模型可以为每个特征 Φi 赋予一个特征重要性值，最后可以得到长度为 m 的特征重要性向量。为了得到简约的模型，PS-Tree 将个体大小附加到适应度向量作为另一个目标，因此拟合向量的长度将变成 m + 1。
PS-Tree 使用 NSGA-II 作为默认选择操作符，由于目标的大小会影响 NSGA-II 中的拥挤距离计算，因此需要用平均大小来归一化个体大小，并将其缩放 0.01。缩放个体大小的公式如下：
式子中的符号和含义如下表所示：
| 符号 | 说明 |
|---|---|
| size(Φi) | 个体 Φi 的大小 |
| avg_size(Φ) | 所有个体的平均大小 |
根据归一化后的个体大小，可以将目标函数定义为如下的式子：
选择算子
选择算子从父个体 Φold 和新个体 Φ 的并集中选择 |Φ| 个体作为新的父个体。本文使用 NSGA-II 作为选择算子，对于个体 Φi 来说，如果在所有目标上 Φi 的表现至少和 Φj 一样好，并且在至少一个目标上 Φi 的表现严格优于 Φj，那么它就优于 Φj。根据这一定义，NSGA-II 根据其优势关系对所有个体进行排名，然后对于那些具有相同级别的个体使用拥挤距离(crowding distance)进行排名。选择了一组个体 Φ 之后，此处使用基于优势的锦标赛选择算子通过交叉变异算子生成后代。
实验分析
数据集
使用 PMLB 基准集中的 122 个回归数据集，最小的数据集有 47 个数据项，而最大的数据集有 1025010 个数据项，特征数量为 2 到 124 不等。训练数据集先进行标准化，然后使用 80% 的数据项来训练回归函数，剩余的 20% 的数据项作为测试集。在消融实验选择的是 200 到 10,000 个数据项之间的共 83 个数据集，然后运行 50 次。算法比较部分，使用 PMLB 中所有 122 个数据集。
PS-Tree 实验设置
PS-Tree 的参数和实验设置如下表所示，其中为了避免零除法误差，用分析商算子(AQ)代替除法算子。AQ(a,b) = a/√(1+b^2)，其中 a 和 b 是除法运算的两个参数。在结果描述方面，论文使用符号 “+”、“∼” 和 “-” 来表示算法显著优于、类似于或劣于其他算法或。
| 参数 | 设置 |
|---|---|
| 评价指标 | R^2 |
| 最大分区数 | 4 |
| 种群大小 | 25 |
| 最大代数 | 50 |
| 特征树的最大高度 | 6 |
| 常数的随机常数范围 | [−2, 2] |
| 新生成的子树的高度 | [1, 3] |
| primitive_function | +、−、∗、AQ |
| 交叉概率 | 0.9 |
| 突变概率 | 0.1 |
消融实验
消融实验主要用来验证以下五点：
- 使用 GP 增强 PL-Tree 的优势；
- 空间分区的优势和各种空间分区技术的影响；
- 自适应空间划分的优势；
- 使用 EA 的优势和各种选择操作的影响；
- 使用 Ridge 而不是 Lasso 的理由；
总体而言，消融实验说明使用 GP 进行特征构造是 PS-Tree 的必要组成部分，空间划分也起到了重要的作用。自适应分区(adaptive partitioning, AP)和 MOEA 等策略也可以稍微提高性能，使用 Ridge 而不是 Lasso 也是合理的。
特征构建
特征构建的消融实验比较了 PS-Tree 和 PL-Tree，实验用相同的方式划分特征空间，PL-Tree 基于每个分区的原始特征构建线性模型，PS-Tree 则进化几个非线性特征来构建局部模型。实验结果如下表所示，从训练数据的预测性能来看，特征构建策略可以有效地提高模型的拟合能力。
空间分区
空间分区的消融实验主要包括 3 个内容：
- 空间分区的优势；
- 自适应分区的优势；
- 不同分区数量的影响。
对于初始空间分区策略，首先选择了四种初始划分策略：决策树、k-means、Bayesian Gaussian mixture、不划分。下图为四个数据集的训练和测试曲线，可以看出空间分区有助于提高性能，但很难确定最佳的初始空间分区策略。
下表所示的实验结果支持这一结论，即空间分区策略可以显著提高 32 个数据集上的性能，不同初始空间分区策略之间没有显著差异。
对于自适应分区的效果，下表的结果表明，自适应分区是缓解过度分割带来的过拟合问题的有效方法。但是因为这样过大的小分区可能会导致局部模型不准确，所以它也增加了过拟合的风险。
对于分区数量的影响，从下表可以看出 PS-Tree 的预测性能随着分区数的增加而提高。但即使使用自适应分区策略，更大数量的分区可能会导致过拟合训。根据实验结果，本文选择 4 作为 PS-Tree 的默认分区数。
多目标优化
为了证明多目标优化的好处，此处选用了四种不同的选择算子，分别是 NSGA-II、SPEA-II、锦标赛选择、随机选择。根据数据集的大小将所有数据集分为四组，然后在下图呈现训练和测试 loss 分布。结论是三种基于进化的选择算子在训练 loss 和测试 loss 没有显著差异，随机选择的表现明显更差。
下表中给出了所有实验结果，结论是与随机搜索方案相比，多目标训练方案可以在四分之三的数据集上显著提高测试损失性能，多目标训练方案略优于单目标。
适应度评估
GP 的相关研究通常使用 Lasso 回归的结果来计算每个个体的适应度值，但根据下表所示的实验结果，使用 Lasso 回归代替 Ridge 回归没有优势。
但是因为作者希望 PS-Tree 模型是可解释的，所以在进化结束时应用 Lasso 回归来获得一组简约特征。下图展示了特征选择过程之前和之后的特征，可见 Lasso 回归将特征的数量从 26 个减少到 15 个，而性能只降低了很小的程度，进一步提高了可解释性。
与其他算法的比较
对比实验中，PS-Tree 与其他 21 种算法在 122 个数据集上进行比较。所有算法的超参数，如 XGBoost 中的树数和 Operon 中的种群大小都通过对半网格搜索方法进行了微调。本文中作者直接使用了 SRBench 提供的实验脚本和实验结果，确保提出的算法和最先进的算法之间的公平比较。
为了提高性能，PS-Tree 的代数增加到 500，并在 primitive_function 中包含 sin 和 cos 函数。对所有数据集使用默认超参数，不进行任何超参数调优。下图显示了测试精度、模型大小和训练时间方面的实验结果，可以看出 PS-Tree 在精度上优于其他 22 种算法。在模型大小方面，PS-Tree 的规模类似于 SBP-GP，比经典遗传规划方法大一个数量级。这是因为 PS-Tree 使用整个种群作为最终模型，而不是只使用单个最优个体。这种方式的好处是 PS-Tree 的训练时间比其他符号回归算法短，比其他基于 python 的符号回归包性能好一个数量级，与 C++ 实现一些算法也可相提并论。
下图给出了不同算法之间模型大小和测试精度的成对显著性检验，PS-Tree 在测试精度方面与最先进的符号回归方法 Operon 有显著差异。
下图中给出了测试精度的临界差值图，图中数字表示特定算法在所有数据集中的平均排名。如果任意两个算法没有显著性差异，将它们用一条线连接起来。排名第一和第二的算法分别是 PS-Tree 和 Operon，显著差异结果与上一张图相似。PS-Tree 与其他算法之间没有界限，表明 PS-Tree 在测试精度方面与所有其他算法有显著不同。
优点和创新点
个人认为，本文有如下一些优点和创新点可供参考学习：
- PS-Tree 结合了分类树、非线性回归和 GP，这几种不同的算法构成了强烈的耦合，缺一不可，都为问题的求解提供了正面的效益，非常有创意；
- 对于特征空间分区部分使用了自适应的策略，在训练过程中可以动态调整，不仅能修正一定的分区错误，更能得到简洁的模型；
- 为了避免在适应度评估阶段带来的巨大时空开销，本文将评估问题转换成了一个多目标优化问题，从而能得到个体的适应度；
- 本文通过实验为每个部分选择了合适的组件，从而使整个算法的实现带来了可行性；
- 模型的设计本身就具有良好的可解释性，同时通过合适的设置使模型更加简洁，且不失可用性；
- 实验数据丰富有力，在上百个数据集上进行了实验，并选择了大量的对比模型和做了大量的消融实验证明算法的优势，说服力强。