超参数优化之网格优化¶

In [13]:
import numpy as np
import pandas as pd
import sklearn
import matplotlib as mlp
import matplotlib.pyplot as plt
import seaborn as sns
import time
#import re, pip, conda
In [14]:
"""
for package in [sklearn,mlp,np,pd,sns,pip,conda]:
    print(re.findall("([^']*)",str(package))[2],package.__version__)
"""
Out[14]:
'\nfor package in [sklearn,mlp,np,pd,sns,pip,conda]:\n    print(re.findall("([^\']*)",str(package))[2],package.__version__)\n'
In [15]:
#pip install --upgrade scikit-learn
#conda update scikit-learn

目录

一 超参数优化与枚举网格的理论极限
二 随机网格搜索（RandomSearchCV）
  1 随机网格搜索的基本原理
  2 随机网格搜索实现
  3 随机网格的理论极限
三 Halving网格搜索（HalvingSearchCV）
  1 对半网格搜索的基本流程
  2 对半网格搜索的实现
  3 对半随机网格搜索
四 【加餐】AutoML前沿进展
  1 自动化第三方框架
  2 AutoML三大成熟研究领域
  3 AutoML的新兴研究

一 超参数优化与枚举网格的理论极限¶

  • 超参数优化HPO（HyperParameter Optimization）

每一个机器学习算法都会有超参数，而超参数的设置很大程度上影响了算法实际的使用效果，因此调参是机器学习算法工程师最为基础和重要的任务。现代机器学习与深度学习算法的超参数量众多，不仅实现方法异常灵活、算法性能也受到更多的参数的复合影响，因此当人工智能浪潮来临时，可以自动选择超参数的超参数优化HPO领域也迎来了新一轮爆发。

在算法的世界中，我们渴望一切流程最终都走向完美自动化，专门研究机器学习自动化的学科被称为AutoML，而超参数自动优化是AutoML中最成熟、最深入、也是最知名的方向。理论上来说，当算力与数据足够时，HPO的性能一定是超过人类的。HPO能够降低人为工作量，并且HPO得出的结果比认为搜索的复现可能性更高，所以HPO可以极大程度提升科学研究的复现性和公平性。当代超参数优化算法主要可以分为：

基于网格的各类搜索（Grid）——最易入门
基于贝叶斯优化的各类优化算法（Baysian）——最为成熟
基于梯度的各类优化（Gradient-based）——研究中...
基于种群的各类优化（进化算法，遗传算法等）——研究中...

其中，各类网格搜索方法与基于贝叶斯的优化方法是最为盛行的，贝叶斯优化方法甚至可以被称为是当代超参数优化中的SOTA模型。这些模型对于复杂集成算法的调整有极大的作用与意义。

  • 网格搜索的理论极限与缺点

在所有超参数优化的算法当中，枚举网格搜索是最为基础和经典的方法。在搜索开始之前，我们需要人工将每个超参数的备选值一一列出，多个不同超参数的不同取值之间排列组合，最终将组成一个参数空间（parameter space）。枚举网格搜索算法会将这个参数空间当中所有的参数组合带入模型进行训练，最终选出泛化能力最强的组合作为模型的最终超参数。

对网格搜索而言，如果参数空间中的某一个点指向了损失函数真正的最小值，那枚举网格搜索时一定能够捕捉到该最小值以及对应的参数（相对的，假如参数空间中没有任意一点指向损失函数真正的最小值，那网格搜索就一定无法找到最小值对应的参数组合）。

参数空间越大、越密，参数空间中的组合刚好覆盖损失函数最小值点的可能性就会越大。这是说，极端情况下，当参数空间穷尽了所有可能的取值时，网格搜索一定能够找到损失函数的最小值所对应的最优参数组合，且该参数组合的泛化能力一定是强于人工调参的。

但是，参数空间越大，网格搜索所需的算力和时间也会越大，当参数维度上升时，网格搜索所需的计算量更是程指数级上升的。以随机森林为例：

只有1个参数n_estimators，备选范围是[50,100,150,200,250,300]，需要建模6次。
增加参数max_depth，且备选范围是[2,3,4,5,6]，乘以5，需要建模30次。
增加参数min_sample_split，且备选范围为[2,3,4,5]，再乘以4，一共需要建模120次。

同时，参数优化的目标是找出令模型泛化能力最强的组合，因此需要交叉验证来体现模型的泛化能力，假设交叉验证次数为5，则三个参数就需要建模600次。在面对超参数众多、且超参数取值可能无限的人工神经网络、融合模型、集成模型时，伴随着数据和模型的复杂度提升，网格搜索所需要的时间会急剧增加，完成一次枚举网格搜索可能需要耗费几天几夜。考虑到后续实践过程中，算法和数据都将更加复杂，而建模过程中超参数调优是模型训练的必备环节，因此，我们急需寻找到一种更加高效的超参数搜索方法。在本节课中，我们将介绍三种基于网格进行改进的超参数优化方法，并将他们的结果与网格搜索进行时间/空间/效果上的对比。

  • 建立benchmark：随机森林中枚举网格搜索的结果
In [16]:
from sklearn.ensemble import RandomForestRegressor as RFR
from sklearn.model_selection import cross_validate, KFold, GridSearchCV
In [18]:
data = pd.read_csv(r"..\Lesson 09.随机森林模型\datasets\House Price\train_encode.csv",index_col=0)
In [19]:
X = data.iloc[:,:-1] # 特征
y = data.iloc[:,-1]  # 标签
In [20]:
X.shape
Out[20]:
(1460, 80)
In [21]:
X.head()
Out[21]:
Id 住宅类型 住宅区域 街道接触面积(英尺) 住宅面积 街道路面状况 巷子路面状况 住宅形状(大概) 住宅现状 水电气 ... 半开放式门廊面积 泳池面积 泳池质量 篱笆质量 其他配置 其他配置的价值 销售月份 销售年份 销售类型 销售状态
0 0.0 5.0 3.0 36.0 327.0 1.0 0.0 3.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 1.0 2.0 8.0 4.0
1 1.0 0.0 3.0 51.0 498.0 1.0 0.0 3.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.0 1.0 8.0 4.0
2 2.0 5.0 3.0 39.0 702.0 1.0 0.0 0.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 8.0 2.0 8.0 4.0
3 3.0 6.0 3.0 31.0 489.0 1.0 0.0 0.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 8.0 0.0
4 4.0 5.0 3.0 55.0 925.0 1.0 0.0 0.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 11.0 2.0 8.0 4.0

5 rows × 80 columns

In [22]:
y.head() # 标签
Out[22]:
0    208500
1    181500
2    223500
3    140000
4    250000
Name: SalePrice, dtype: int64
In [50]:
y.describe() #RMSE
Out[50]:
count      1460.000000
mean     180921.195890
std       79442.502883
min       34900.000000
25%      129975.000000
50%      163000.000000
75%      214000.000000
max      755000.000000
Name: SalePrice, dtype: float64
In [52]:
#参数空间
param_grid_simple = {"criterion": ["squared_error","poisson"]
                     , 'n_estimators': [*range(20,100,5)]
                     , 'max_depth': [*range(10,25,2)]
                     , "max_features": ["log2","sqrt",16,32,64,"auto"]
                     , "min_impurity_decrease": [*np.arange(0,5,10)]
                    }
In [53]:
#参数空间大小计算

2 * len([*range(20,100,5)]) * len([*range(10,25,2)]) * len(["log2","sqrt",16,32,64,"auto"]) * len([*np.arange(0,5,10)])
Out[53]:
1536
In [54]:
#直接使用循环计算
no_option = 1
for i in param_grid_simple:
    no_option *= len(param_grid_simple[i])
In [55]:
no_option # 最终备选项数目，组合数量
Out[55]:
1536
In [ ]:
#模型，交叉验证，网格搜索
reg = RFR(random_state=1412,verbose=True,n_jobs=-1)
cv = KFold(n_splits=5,shuffle=True,random_state=1412)
search = GridSearchCV(estimator=reg
                     ,param_grid=param_grid_simple
                     ,scoring = "neg_mean_squared_error"
                     ,verbose = True
                     ,cv = cv
                     ,n_jobs=-1)
In [7]:
#=====【TIME WARNING: 7mins】=====#
start = time.time()
search.fit(X,y)
print(time.time() - start)

Fitting 5 folds for each of 1536 candidates, totalling 7680 fits
381.6039867401123

[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=-1)]: Done  18 tasks      | elapsed:    0.0s
[Parallel(n_jobs=-1)]: Done  85 out of  85 | elapsed:    0.0s finished
In [11]:
381.6039/60 # 运行分钟
Out[11]:
6.3600650000000005
In [10]:
search.best_estimator_
Out[10]:
RandomForestRegressor(max_depth=23, max_features=16, min_impurity_decrease=0,
                      n_estimators=85, n_jobs=-1, random_state=1412,
                      verbose=True)
In [9]:
abs(search.best_score_)**0.5 # 最优RMSE
Out[9]:
29179.698261599166
In [11]:
# 按最优参数重建模型，查看效果
ad_reg = RFR(n_estimators=85, max_depth=23, max_features=16, random_state=1412)
In [15]:
cv = KFold(n_splits=5,shuffle=True,random_state=1412)
result_post_adjusted = cross_validate(ad_reg,X,y,cv=cv,scoring="neg_mean_squared_error"
                          ,return_train_score=True
                          ,verbose=True
                          ,n_jobs=-1)

[Parallel(n_jobs=-1)]: Using backend LokyBackend with 16 concurrent workers.
[Parallel(n_jobs=-1)]: Done   5 out of   5 | elapsed:    0.2s finished
In [56]:
def RMSE(cvresult,key):
    return (abs(cvresult[key])**0.5).mean()
In [16]:
RMSE(result_post_adjusted,"train_score")
Out[16]:
11000.81099038192
In [17]:
RMSE(result_post_adjusted,"test_score")
Out[17]:
28572.070208366855
HPO方法 默认参数 网格搜索
搜索空间/全域空间 - 1536/1536
运行时间（分钟） - 6.36
搜索最优（RMSE） 30571.266 29179.698
重建最优（RMSE） - 28572.070
In [12]:
# 打包成函数供后续使用
# 评估指标RMSE
def RMSE(cvresult,key):
    return (abs(cvresult[key])**0.5).mean()

# 计算参数空间大小
def count_space(param):
    no_option = 1
    for i in param_grid_simple:
        no_option *= len(param_grid_simple[i])
    print(no_option)
    
# 在最优参数上进行重新建模验证结果
def rebuild_on_best_param(ad_reg):
    cv = KFold(n_splits=5,shuffle=True,random_state=1412)
    result_post_adjusted = cross_validate(ad_reg,X,y,cv=cv,scoring="neg_mean_squared_error"
                                          ,return_train_score=True
                                          ,verbose=True
                                          ,n_jobs=-1)
    print("训练RMSE:{:.3f}".format(RMSE(result_post_adjusted,"train_score")))
    print("测试RMSE:{:.3f}".format(RMSE(result_post_adjusted,"test_score")))

二 随机网格搜索（RandomizedSearchCV）¶

  • 基本原理

在讲解网格搜索时我们提到，伴随着数据和模型的复杂度提升，网格搜索所需要的时间急剧增加。以随机森林算法为例，如果使用过万的数据，搜索时间则会立刻上升好几个小时。因此，我们急需寻找到一种更加高效的超参数搜索方法。

首先，当所使用的算法确定时，决定枚举网格搜索运算速度的因子一共有两个：

1 参数空间的大小：参数空间越大，需要建模的次数越多
2 数据量的大小：数据量越大，每次建模时需要的算力和时间越多

因此，sklearn中的网格搜索优化方法主要包括两类，其一是调整搜索空间——即随机网格搜索，其二是调整每次训练的数据。其中，调整参数空间的具体方法，是放弃原本的搜索中必须使用的全域超参数空间，改为挑选出部分参数组合，构造超参数子空间，并只在子空间中进行参数搜索。

以下图的二维空间为例，在这个n_estimators与max_depth共同组成的参数空间中，n_estimators的取值假设为[50,100,150,200,250,300]，max_depth的取值假设为[2,3,4,5,6]，则枚举网格搜索必须对30种参数组合都进行搜索。当我们调整搜索空间，我们可以只抽样出橙色的参数组合作为“子空间”，并只对橙色参数组合进行搜索。如此一来，整体搜索所需的计算量就大大下降了，原本需要30次建模，现在只需要8次建模。

In [26]:
fig, [ax1, ax2] = plt.subplots(1,2,dpi=300)
n_e_list = [*range(50,350,50)]
m_d_list = [*range(2,7)]
comb = pd.DataFrame([(n_estimators, max_depth) for n_estimators in n_e_list for max_depth in m_d_list])

ax1.scatter(comb.iloc[:,0],comb.iloc[:,1],cmap="Blues")
ax1.set_xticks([*range(50,350,50)])
ax1.set_yticks([*range(2,7)])
ax1.set_xlabel("n_estimators")
ax1.set_ylabel("max_depth")
ax1.set_title("GridSearch")

ax2.scatter(comb.iloc[:,0],comb.iloc[:,1],cmap="Blues")
ax2.scatter([50,250,200,200,300,100,150,150],[4,2,6,3,2,3,2,5],cmap="red",s=20,linewidths=5)
ax2.set_xticks([*range(50,350,50)])
ax2.set_yticks([*range(2,7)])
ax2.set_xlabel("n_estimators")
ax2.set_ylabel("max_depth")
ax2.set_title("RandomSearch");

在sklearn中，随机抽取参数子空间并在子空间中进行搜索的方法叫做随机网格搜索RandomizedSearchCV。由于搜索空间的缩小，需要枚举和对比的参数组的数量也对应减少，整体搜索耗时也将随之减少，因此我们得到三个结论：

当设置相同的全域空间时，随机搜索的运算速度比枚举网格搜索快很多。

当设置相同的训练次数时，随机搜索可以覆盖的空间比枚举网格搜索大很多。

同时，绝妙的是，随机网格搜索得出的最小损失与枚举网格搜索得出的最小损失很接近。

可以说，是提升了运算速度，又没有过多地伤害搜索的精度。

不过，需要注意的是，随机网格搜索在实际运行时并不是先抽样出子空间，再对子空间进行搜索，而是仿佛“循环迭代”一般，在这一次迭代中随机抽取1组参数进行建模，下一次迭代再随机抽取1组参数进行建模，由于这种随机抽样是不放回的，因此不会出现两次抽中同一组参数的问题。我们可以控制随机网格搜索的迭代次数，来控制整体被抽出的参数子空间的大小，这种做法往往被称为“赋予随机网格搜索固定的计算量，当全部计算量被消耗完毕之后，随机网格搜索就停止”。

  • 随机网格搜索的实现
In [9]:
from sklearn.model_selection import RandomizedSearchCV

class sklearn.model_selection.RandomizedSearchCV(estimator, param_distributions, *, n_iter=10, scoring=None, n_jobs=None, refit=True, cv=None, verbose=0, pre_dispatch='2*n_jobs', random_state=None, error_score=nan, return_train_score=False)

全部参数解读如下，其中加粗的是随机网格搜索独有的参数：

Name Description
estimator 调参对象，某评估器
param_distributions 全域参数空间，可以是字典或者字典构成的列表
n_iter 迭代次数，迭代次数越多，抽取的子参数空间越大
scoring 评估指标，支持同时输出多个参数
n_jobs 设置工作时参与计算的线程数
refit 挑选评估指标和最佳参数，在完整数据集上进行训练
cv 交叉验证的折数
verbose 输出工作日志形式
pre_dispatch 多任务并行时任务划分数量
random_state 随机数种子
error_score 当网格搜索报错时返回结果，选择'raise'时将直接报错并中断训练过程，其他情况会显示警告信息后继续完成训练
return_train_score 在交叉验证中是否显示训练集中参数得分

我们依然借用之前在网格搜索上见过的X和y，以及随机森林回归器，来实现随机网格搜索：

In [11]:
X.shape
Out[11]:
(1460, 80)
In [12]:
X.head()
Out[12]:
Id 住宅类型 住宅区域 街道接触面积(英尺) 住宅面积 街道路面状况 巷子路面状况 住宅形状(大概) 住宅现状 水电气 ... 半开放式门廊面积 泳池面积 泳池质量 篱笆质量 其他配置 其他配置的价值 销售月份 销售年份 销售类型 销售状态
0 0.0 5.0 3.0 36.0 327.0 1.0 0.0 3.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 1.0 2.0 8.0 4.0
1 1.0 0.0 3.0 51.0 498.0 1.0 0.0 3.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 4.0 1.0 8.0 4.0
2 2.0 5.0 3.0 39.0 702.0 1.0 0.0 0.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 8.0 2.0 8.0 4.0
3 3.0 6.0 3.0 31.0 489.0 1.0 0.0 0.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 8.0 0.0
4 4.0 5.0 3.0 55.0 925.0 1.0 0.0 0.0 3.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 11.0 2.0 8.0 4.0

5 rows × 80 columns

In [13]:
y.describe()
Out[13]:
count      1460.000000
mean     180921.195890
std       79442.502883
min       34900.000000
25%      129975.000000
50%      163000.000000
75%      214000.000000
max      755000.000000
Name: SalePrice, dtype: float64
  • 相同的全域参数空间
In [57]:
# 定义全域参数空间 - 使用与网格搜索时完全一致的空间，以便于对比
param_grid_simple = {"criterion": ["squared_error","poisson"]
                     , 'n_estimators': [*range(20,100,5)]
                     , 'max_depth': [*range(10,25,2)]
                     , "max_features": ["log2","sqrt",16,32,64,"auto"]
                     , "min_impurity_decrease": [*np.arange(0,5,10)]
                    }

# 建立用于训练的回归器、交叉验证
reg = RFR(random_state=1412, verbose=True, n_jobs=-1)
cv = KFold(n_splits=5, shuffle=True, random_state=1412)
In [58]:
# 计算全域参数空间大小，这是我们能够抽样的最大值
count_space(param_grid_simple)

1536
In [59]:
# 定义随机搜索
search = RandomizedSearchCV(estimator=reg
                            ,param_distributions=param_grid_simple
                            ,n_iter = 800 #设置的子空间的大小是全域空间的一半左右
                            ,scoring = "neg_mean_squared_error"
                            ,verbose = True
                            ,cv = cv
                            ,random_state=1412
                            ,n_jobs=-1
                           )
In [17]:
# 训练随机搜索评估器
#=====【TIME WARNING: 5~10min】=====#
start = time.time()
search.fit(X,y)
print(time.time() - start)

Fitting 5 folds for each of 800 candidates, totalling 4000 fits
170.16785073280334

[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 16 concurrent workers.
[Parallel(n_jobs=-1)]: Done  18 tasks      | elapsed:    0.0s
[Parallel(n_jobs=-1)]: Done  85 out of  85 | elapsed:    0.0s finished
In [21]:
170.1678/60 # 耗时分钟
Out[21]:
2.83613
In [19]:
# 查看模型结果，即最优参数
search.best_estimator_
Out[19]:
RandomForestRegressor(max_depth=24, max_features=16, min_impurity_decrease=0,
                      n_estimators=85, n_jobs=-1, random_state=1412,
                      verbose=True)
In [20]:
abs(search.best_score_)**0.5 # 最好的RMSE
Out[20]:
29251.284326350575
In [23]:
# 根据最优参数重建模型
ad_reg = RFR(max_depth=24, max_features=16, min_impurity_decrease=0,
                      n_estimators=85, n_jobs=-1, random_state=1412,
                      verbose=True)
# 调用自定义的函数
rebuild_on_best_param(ad_reg)

[Parallel(n_jobs=-1)]: Using backend LokyBackend with 16 concurrent workers.

训练RMSE:11031.299
测试RMSE:28639.969

[Parallel(n_jobs=-1)]: Done   5 out of   5 | elapsed:    0.1s finished
HPO方法 默认参数 网格搜索 随机搜索
搜索空间/全域空间 - 1536/1536 800/1536
运行时间（分钟） - 6.36 **2.83(↓)**
搜索最优（RMSE） 30571.266 29179.698 29251.284
重建最优（RMSE） - 28572.070 **28639.969(↑)**

很明显，在相同参数空间、相同模型的情况下，随机网格搜索的运算速度是普通网格搜索的一半，当然，这与子空间是全域空间的一半有很大的联系。由于随机搜索只是降低搜索的次数，并非影响搜索过程本身，因此其运行时间基本就等于n_iter/全域空间组合数 * 网格搜索的运行时间。

  • 随机网格搜索的理论极限

虽然通过缩小子空间可以提升搜索的速度，但是随机网格搜索的精度看起来并没有削减太多，甚至是可以原谅的，那随机网格搜索可以得到和网格搜索一样好的结果吗？它也像网格搜索一样，可以得到最优的参数组合吗？为什么缩小参数空间之后，随机网格搜索的结果还与网格搜索一致？

理论上来说，枚举网格搜索的上限和随机网格搜索的上限哪个高？

从直觉上来说，我们很难回答这些问题，但我们可以从数学的随机过程的角度来理解这个问题。在机器学习算法当中，有非常多通过随机来提升运算速度（比如Kmeans，随机挑选样本构建簇心，小批量随机梯度下降，通过随机来减少每次迭代需要的样本）、或通过随机来提升模型效果的操作（比如随机森林，比如极度随机树）。两种随机背后的原理完全不同，而随机网格搜索属于前者，这一类机器学习方法总是伴随着“从某个全数据集/全域中进行抽样”的操作，而这种操作能够有效的根本原因在于：

  • 抽样出的子空间可以一定程度上反馈出全域空间的分布，且子空间相对越大（含有的参数组合数越多），子空间的分布越接近全域空间的分布

  • 当全域空间本身足够密集时，很小的子空间也能获得与全域空间相似的分布

  • 如果全域空间包括了理论上的损失函数最小值，那一个与全域空间分布高度相似的子空间很可能也包括损失函数的最小值，或包括非常接近最小值的一系列次小值

上述数学事实听起来比较抽象，但其实我们可以通过绘制图像来直观地呈现这些事实。许多在数学上比较抽象的概念都可以被可视化。在这里，我们借助matplotlib工具库mplot3d中的一组默认数据。

In [23]:
from mpl_toolkits.mplot3d import axes3d
In [24]:
# 自动获取数据的功能get_test_data，可以自动生成复合某一分布的数据
p1, p2, MSE = axes3d.get_test_data(0.05)

# 我们现在假设这一组数据中有两个参数，p1与p2，两个参数组成的参数组合对应着损失函数值MSE
# 参数0.05是指参数空间中，点与点之间的距离
# 因此该数字越小，取出来的样本越多
In [25]:
len(p1) #参数1的取值有120个
Out[25]:
120
In [26]:
len(p2) #参数2的取值也有120个
Out[26]:
120

请问现在参数空间当中一共有多少个参数组合？120*120=14400种组合，所以参数空间中一共有14400个点。

In [27]:
MSE.shape #损失函数值，总共14400个点
Out[27]:
(120, 120)
In [28]:
#绘制P1与P2的参数空间 - 这是一个呈现出14400个点的密集空间
plt.figure(dpi=300)
plt.scatter(p1,p2,s=0.2)
plt.xticks(fontsize=9)
plt.yticks(fontsize=9);
In [30]:
# 参数与损失共同构建的函数图像
p1, p2, MSE = axes3d.get_test_data(0.05)
plt.figure(dpi=300)
ax = plt.axes(projection="3d")
ax.plot_wireframe(p1,p2,MSE,rstride=2,cstride=2,linewidth=0.5)
#ax.view_init(2, -15) # 旋转下3D图，看的更清楚
ax.zaxis.set_tick_params(labelsize=7)
ax.xaxis.set_tick_params(labelsize=7)
ax.yaxis.set_tick_params(labelsize=7);