自研贝叶斯优化算法遇到的坑

自研贝叶斯优化算法，如何判断算法能拟合？我目前是在一个tiny的数据集上跑一下，看算法能否收敛到正确的的局部最小值。这里要有两个关键词：

1. 收敛。算法是需要收敛的。黑盒优化的本质就是增加在优势样本附近的采样率。如果算法如同随机搜索不收敛，那么是有问题的。
2. 正确。收敛点是正确的，如果收敛到错误的点，那还不如随机搜索。

自研SMAC

代理模型

SMAC的本质是用随机森林作为代理模型。这个代理模型调包就好了（前提是你熟读开源代码千百遍，知道调什么）

from skopt.learning.forest import RandomForestRegressor, ExtraTreesRegressor

众所周知，RandomForestRegressor不仅会对行做采样，也会对列做采样。ExtraTreesRegressor只会对行做采样。就我使用的经验来看，特征>样本的情况适合用RF，其余情况一般用ET。SMAC的论文提到，他使用的随机森林会用所有的样本，但特征的采样率是83%. （SMAC源码分析->代理模型的构建）

从skopt文档来看，一般来说ET要比RF表现好。

获益函数

就我的经验来看，EI肯定要比PI好，因为EI算的是期望，利用的比重其实比探索要大。PI更注重探索，更发散。

skopt实现了EI, PI, LCB, EIPS等。我目前实现了EI，LogEI。

LogEI的代码可以看SMAC，RoBO

就实验来看，感觉LogEI和EI差别不大。个人感觉上，无loss_transform+EI == log_scaled loss_transform+LogEI

SMAC的LogEI貌似就是搭配 log_scaled loss_transform的。

RoBO的LogEI与SMAC的实现有很大不同。

使用上来看，EI+ log_scaled loss_transform即可，xi($\xi$ξ) 这个参数设0.01，0好像都没什么区别

class EI():
    def __init__(self, xi=0.01):
        # in SMAC, xi=0.0,
        # smac/optimizer/acquisition.py:341
        # par: float=0.0
        # in scikit-optimize, xi=0.01
        # this blog recommend xi=0.01
        # http://krasserm.github.io/2018/03/21/bayesian-optimization/
        self.xi = xi

    def __call__(self, model, X, y_opt):
        mu, std = model.predict(X, return_std=True)
        values = np.zeros_like(mu)
        mask = std > 0
        improve = y_opt - self.xi - mu[mask]
        scaled = improve / std[mask]
        cdf = norm.cdf(scaled)
        pdf = norm.pdf(scaled)
        exploit = improve * cdf
        explore = std[mask] * pdf
        values[mask] = exploit + explore
        # You can find the derivation of the EI formula in this blog
        # http://ash-aldujaili.github.io/blog/2018/02/01/ei/
        return values

class LogEI():
    def __init__(self, xi=0.01):
        self.xi = xi

    def __call__(self, model, X, y_opt):
        mu, std = model.predict(X, return_std=True)
        var = std ** 2
        values = np.zeros_like(mu)
        mask = std > 0
        f_min = y_opt - self.xi
        improve = f_min - mu[mask]
        # in SMAC, v := scaled
        # smac/optimizer/acquisition.py:388
        scaled = improve / std[mask]
        values[mask] = (np.exp(f_min) * norm.cdf(scaled)) - \
                       (np.exp(0.5 * var[mask] + mu[mask]) * norm.cdf(scaled - std[mask]))
        return values

遇到过的坑

不当的loss scale会导致算法陷入错误的局部最优

experiment_id=39

实验代码：

from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split


X, y = load_digits(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)
pipe = SystemClassifier(
    DAG_workflow={
        "num->target": [
            "liblinear_svc",
            "libsvm_svc",
            "logistic_regression",
            "random_forest",
            # "catboost",
        ]
    },
    config_generator="ET",
    config_generator_params={
        "acq_func": "EI",
        "xi": 0,
        "loss_transformer":None,
        "min_points_in_model": 20
    },
    warm_start=False,
    random_state=0,
    min_n_samples_for_SH=50,
    concurrent_type="thread",
    # max_budget=1,
    n_jobs_in_algorithm=3,
    n_workers=1,
    SH_only=True,
    min_budget=1/16,
    max_budget=1/16,
    n_iterations=100,
    # min_budget=1 / 4,
    debug_evaluator=True,
)
pipe.fit(
    X_train, y_train,
    # is_not_realy_run=True,
    fit_ensemble_params=False)
# score = accuracy_score(y_test, y_pred)
score = pipe.score(X_test, y_test)
print(score)

实验结果：
红色部分是warming up（启动过程，即开始时的随机搜索，20次），蓝色是贝叶斯搜索。可以看到启动时最好的样本是 (random_forest, 0.977) ，但是算法却陷入了libsvm_svc的局部最优，真是匪夷所思。

排查后认为，是没有对loss（也就是代理模型拟合的label）做log_scaled 变换。

class LogScaledLossTransformer(LossTransformer):
    def fit_transform(self, y, *args):
        y = super(LogScaledLossTransformer, self).fit_transform(y)
        # Subtract the difference between the percentile and the minimum
        y_min = self.y_min - (self.perc - self.y_min)
        y_min -= 1e-10
        # linear scaling
        if y_min == self.y_max:
            # prevent diving by zero
            y_min *= 1 - (1e-10)
        y = (y - y_min) / (self.y_max - y_min)
        y = np.log(y)
        f_max = y[np.isfinite(y)].max()
        f_min = y[np.isfinite(y)].min()
        y[np.isnan(y)] = f_max
        y[y == -np.inf] = f_min
        y[y == np.inf] = f_max
        return y

TPE

我写的TPE基于一个假设：随机变量相互独立。

遇到的坑

experiment_id=39

不对deactivated的随机变量做impute导致算法前期部分收敛，后期收敛到错误的局部最小值

实验代码：

    config_generator="TPE",
    config_generator_params={
        "fill_deactivated_value":False,
        "min_points_in_model": 40
    },

前期
前期对除了random_forest以外的choices有所探索，这是正常现象
后期收敛到了错误的libsvm_svc（为什么都喜欢翻这个算法的牌子，logistic_regression哭晕在厕所）

在predict函数中加了这段代码后：

 if N_deactivated > 0 and self.fill_deactivated_value:
     good_pdf[~mask, i] = np.random.choice(good_pdf_activated)
     bad_pdf[~mask, i] = np.random.choice(bad_pdf_activated)

experiment_id=35

前期
后期
可以看到随着算法不断迭代，TPE的观测增多，于是算法逐渐从探索转为开发，最后变成了在RF附近做局部搜索。