使用拉普拉斯近似贝叶斯优化探索新的超参数维度

当我发现我的模型过度拟合时，我常常想，——是时候正规化了—。但是我如何决定使用哪种正则化方法（L1、L2）以及选择哪些参数？通常，我通过网格搜索来选择设置来执行超参数优化。但是，如果自变量具有不同的规模或不同的影响水平，会发生什么情况？我可以为每个变量设计一个具有不同正则化系数的超参数网格吗？这种类型的优化在高维空间中可行吗？还有其他设计正则化的方法吗？让我们通过一个假设的例子来探讨这一点。

我的虚构示例是一个具有 3 个解释变量的二元分类用例。每个变量都是分类变量，有 7 个不同的类别。我的可重现用例是这样的笔记本。生成数据集的函数如下：

将 numpy 导入为 np
将 pandas 导入为 pd
def get_classification_dataset():
n_样本 = 200
猫= [“a”，“b”，“c”，“d”，“e”，“f”]
X = pd.DataFrame(
数据={
“col1”：np.random.choice（猫，大小= n_samples），
“col2”：np.random.choice（猫，大小= n_samples），
“col3”：np.random.choice（猫，大小= n_samples），
}
）
X_preprocessed = pd.get_dummies(X)
theta = np.random.multivariate_normal(
np.zeros(len(猫) * X.shape[1]),
np.diag(np.array([1e-1] * len(猫) + [1] * len(猫) + [1e1] * len(猫))),）
y = pd.系列(
数据 = np.random.binomial(1, expit(np.dot(X_preprocessed.to_numpy(), theta))),
索引=X_preprocessed.index,）
返回X_预处理后的y
作为参考，我特意为 theta 协方差矩阵选择了 3 个不同的值，以展示拉普拉斯近似贝叶斯优化方法的优势。
如果价值观在某种程度上相似，那么兴趣就会很小。

连同

预测平均观测值的简单基线模型在训练数据集（用于比较目的）上，我选择设计一个稍微复杂的模型。我决定对三个自变量进行一次性编码并应用逻辑回归模型在此基本预处理之上。对于正则化，我选择了 L2 设计，旨在使用两种技术找到最佳正则化系数：网格搜索和拉普拉斯近似贝叶斯优化，正如您现在可能已经预料到的那样。最后，我使用两个指标（任意选择）在测试数据集上评估模型：对数损失和 AUC ROC。

在展示结果之前，我们首先仔细看看贝叶斯模型以及我们如何优化它。

在贝叶斯框架中，参数不再是固定常数，而是随机变量。我们现在不是最大化估计这些未知参数的可能性，而是根据观察到的数据优化随机参数的后验分布。这要求我们通常有些任意地选择先验的设计和参数。然而，也可以将先验参数视为随机变量本身 - 就像盗梦空间，不确定性层层叠加……

在本研究中，我选择了以下模型：

我逻辑上为 Y_i 选择了伯努利模型 |α，对应于 α 的 L2 正则化的正态中心先验 |Σ 最后对于 Σ_i^{-1}，我选择了 Gamma 模型。我选择对精度矩阵而不是协方差矩阵进行建模，因为它在文献中是传统的，例如贝叶斯线性回归的 scikit learn 用户指南 [2]。

除了这个书面模型之外，我还假设 Y_i 和 Y_j 以及 Y_i 和 Σ 是有条件（在 α 上）独立的。

可能性

根据该模型，可能性可以写成：

为了优化，我们需要评估几乎所有项，但 P(Y=y) 除外。分子中的项可以使用所选模型进行评估。但是，分母中的其余项则不能。这就是拉普拉斯近似发挥作用。

拉普拉斯近似

为了计算分母的第一项，我们可以利用拉普拉斯近似。我们近似估计 α | 的分布Y，Σ：

其中 α* 是众数的众数，则 α | 的密度分布Y，Σ。

即使我们不知道密度函数，我们也可以通过以下分解来评估 Hessian 部分：

我们只需要知道分子的前两项就可以计算我们所做的 Hessian 矩阵。

对于那些有兴趣进一步解释的人，我建议阅读 Christopher M. Bishop 的《模式识别和机器学习》中的第 4.4 部分“拉普拉斯近似”[1]。它对我理解近似值有很大帮助。

拉普拉斯近似似然

最后，拉普拉斯近似优化可能性为：

一旦我们近似了 δ 的密度函数 |Y，Σ，如果近似值在任何地方都是准确的，我们最终可以评估我们想要的任何 θ 的可能性。为了简单起见，并且由于近似仅在接近众数时才准确，因此我们在 α* 处评估近似似然。

下面是一个函数，用于评估给定（标量）Δ=1/p（除了给定的观测值 X 和 y 以及设计值 α 和 β 之外）的损失。

将 numpy 导入为 np
从 scipy.stats 导入 gamma
从 module.bayesian_model 导入 BayesianLogisticRegression
def 损失(p, X, y, alpha, beta):
# 计算给定值的损失：
# - 1/sigma²（此处命名为 p 以表示精度）
# - X：特征矩阵
# - y：观察向量
# - alpha：先验 Gamma 分布 alpha 参数超过 1/sigma²
# - beta：先验 Gamma 分布 beta 参数超过 1/sigma²
n_feat = X.shape[1]
m_vec = np.array([0] * n_feat)
p_vec = np.array(p * n_feat)
# 计算 theta*
分辨率 = 最小化（
贝叶斯逻辑回归()._loss,
np.array([0] * n_feat),
args=(X, y, m_vec, p_vec),
方法=“BFGS”，
jac=贝叶斯逻辑回归()._jac,）
theta_star = res.x
# 计算拉普拉斯近似的 Hessian 矩阵
H = 贝叶斯逻辑回归()._hess(theta_star, X, y, m_vec, p_vec)
＃ 损失
损失=0
## 前两项：对数损失和正则化项
损失 += baysian_model._loss(theta_star, X, y, m_vec, p_vec)
## 第三项：sigma 的先验分布，此处写为 p
out -= gamma.logpdf(p, a = alpha, scale = 1 / beta)
## 第四项：拉普拉斯逼近上一项
输出 += 0.5 * np.linalg.slogdet(H)[1] - 0.5 * n_feat * np.log(2 * np.pi)
返回
在我的用例中，我选择通过以下方式对其进行优化

亚当优化器，从中获取了哪些代码回购协议。定义亚当（

乐趣，
x0,
江淮汽车,
参数=(),
学习率=0.001，
β1=0.9，
β2=0.999，
每股收益=1e-8,
启动器=0，
最大=1000，
回调=无，
**夸格
）：
ADAM 的“””``scipy.optimize.minimize`` 兼容实现 -
[http://arxiv.org/pdf/1412.6980.pdf]。
改编自“autograd/misc/optimizers.py”。
”“”
x = x0
m = np.zeros_like(x)
v = np.zeros_like(x)
对于范围内的 i（startiter、startiter + maxiter）：
g = jac(x, *args)
if 回调和回调(x)：
休息
m = (1 - beta1) * g + beta1 * m # 一阶矩估计。
v = (1 - beta2) * (g**2) + beta2 * v # 二阶矩估计。
mhat = m / (1 - beta1**(i + 1)) # 偏差校正。
vhat = v / (1 - beta2**(i + 1))
x = x - 学习率 * mhat / (np.sqrt(vhat) + eps)
我 += 1
返回 OptimizeResult(x=x, fun=fun(x, *args), jac=g, nit=i, nfev=i, success=True)
对于此优化，我们需要先前损失的导数。

我们无法拥有分析形式，因此我决定使用导数的数值近似。

一旦模型在训练数据集上进行训练，就需要对评估数据集进行预测，以评估其性能并比较不同的模型。然而，不可能直接计算新点的实际分布，因为计算很困难。

可以用以下公式来近似结果：

考虑到：

我选择了一个无信息先验而不是精确随机变量。朴素模型表现不佳，对数损失为 0.60，AUC ROC 为 0.50。第二个模型在使用网格搜索和贝叶斯优化进行超优化时表现更好，对数损失为 0.44，AUC ROC 为 0.83。这表明包含因变量的逻辑回归模型优于朴素模型。但是，与网格搜索相比，使用贝叶斯优化没有任何优势，因此我现在将继续使用网格搜索。感谢您的阅读。

……但是等等，我在想。为什么我的参数使用相同的系数进行正则化？我的先验不应该依赖于潜在的因变量吗？也许第一个因变量的参数可以取更高的值，而第二个因变量的参数由于其影响较小，应该更接近于零。让我们探索这些新维度。

到目前为止，我们已经考虑了两种技术：网格搜索和贝叶斯优化。我们可以在更高维度中使用这些相同的技术。

考虑新的维度可能会显着增加网格的节点数量。这就是为什么贝叶斯优化在更高维度上有意义以获得最佳正则化系数。在考虑的用例中，我假设有 3 个正则化参数，每个自变量对应一个。对单个变量进行编码后，我假设生成的新变量都共享相同的正则化参数。因此，即使有超过 3 列作为逻辑回归的输入，总正则化参数也为 3。

我用以下代码更新了之前的损失函数：

将 numpy 导入为 np
从 scipy.stats 导入 gamma
从 module.bayesian_model 导入 BayesianLogisticRegression
def 损失(p, X, y, alpha, beta, X_columns, col_to_p_id):
# 计算给定值的损失：
# - 1/sigma² 向量（此处命名为 p 以表示精度）
# - X：特征矩阵
# - y：观察向量
# - alpha：先验 Gamma 分布 alpha 参数超过 1/sigma²
# - beta：先验 Gamma 分布 beta 参数超过 1/sigma²
# - X_columns：X列的名称列表
# - col_to_p_id：将列名映射到 p 索引的字典
# 因为许多列名可以共享相同的 p 值
n_feat = X.shape[1]
m_vec = np.array([0] * n_feat)
p_列表 = []
对于 X_columns 中的 col：
p_list.append(p[col_to_p_id[col]])
p_vec = np.array(p_list)
# 计算 theta*
分辨率 = 最小化（
贝叶斯逻辑回归()._loss,
np.array([0] * n_feat),
args=(X, y, m_vec, p_vec),
方法=“BFGS”，
jac=贝叶斯逻辑回归()._jac,）
theta_star = res.x
# 计算拉普拉斯近似的 Hessian 矩阵
H = 贝叶斯逻辑回归()._hess(theta_star, X, y, m_vec, p_vec)
＃ 损失
损失=0
## 前两项：对数损失和正则化项
损失 += baysian_model._loss(theta_star, X, y, m_vec, p_vec)
## 第三项：1/sigma² 的先验分布，此处写为 p
## 现在有一个和，因为 p 现在是一个向量
out -= np.sum(gamma.logpdf(p, a = alpha, scale = 1 / beta))
## 第四项：拉普拉斯逼近上一项
输出 += 0.5 * np.linalg.slogdet(H)[1] - 0.5 * n_feat * np.log(2 * np.pi)
返回
使用这种方法，在测试数据集上评估的指标如下：0.39 和 0.88，这比通过网格搜索和贝叶斯方法优化的初始模型要好，所有自变量只有一个先验。

在我的用例中使用不同方法实现的指标。

笔记本。我创建了一个示例来说明该技术的有用性。

然而，我一直无法找到合适的现实世界数据集来充分展示其潜力。当我使用实际数据集时，我无法从应用该技术中获得任何显着的好处。如果您遇到这种情况，请告诉我 - 我很高兴看到这种正则化方法在现实世界中的应用。

总之，使用贝叶斯优化（如果需要，使用拉普拉斯近似）来确定最佳正则化参数可能是传统超参数调整方法的一个很好的替代方法。通过利用概率模型，贝叶斯优化不仅降低了计算成本，而且还提高了找到最佳正则化值的可能性，尤其是在高维情况下。

1. 克里斯托弗·M·毕肖普。（2006）。模式识别和机器学习。施普林格。
2. 贝叶斯岭回归 scikit-learn 用户指南：https://scikit-learn.org/1.5/modules/linear_model.html#bayesian-ridge-regression