学习黑盒优化算法CMA和RandomSearch，借助阿里达摩院MindOpt的RABBO榜单【系列2/4】

1  黑盒优化的概念

上一篇《新手一步步学习黑盒优化算法，借助达摩院MindOpt的RABBO榜单【系列1/4】》中给大家介绍了黑盒优化的概念、和怎么借助阿里达摩院的RABBO“动手”。这篇将着重学习黑盒优化的算法。
​

先复习一下上篇关于雨黑盒优化的内容：
什么是黑盒优化？
「黑盒优化问题」泛指目标函数难以从数学上解析表达，缺少可直接利用的梯度信息，仅可利用目标函数输入和对应输出函数值进行最优解搜索的优化问题。
这类问题要怎么优化求解？
可以 通过获取不同控制参数(输入变量)对应的系统表现，来推断和搜寻优化解。架构就像下面的图这样，我们假设我们要求解的问题不好描述，就把这个问题做成一个仿真系统，这个系统可以根据输入的变量值（“候选解”）来给出对应的评价（“观测值”)，然后黑盒优化算法，就是接上这个仿真系统，通过不断地去提供候选解和得到观测值，来搜索可能的最优表现的候选解。

图注： 云栖大会上王孟昌博士 https://yunqi.aliyun.com/2021/agenda/session170 直播回放大约20分钟的时候截图

2 RABBO中的算法：CMA方法、randomSearch方法。

根据网页介绍：https://tianchi.aliyun.com/specials/promotion/BlackBox 的指引，下载了RABBO仓库。可以参看上一篇文档《【系列1/4】》指导来下载。RABBO的algorithms文件夹里面有两个内置的算法，在目录 blackbox_starter_kit -> algorithms文件夹下面。

2.1 CMA（CMA-ES, 协方差矩阵自适应进化策略)

2.1.1 cma 背景知识

点开RABBO文件夹下面的blackbox_starter_kit/algorithms/cmaes/cma_algorithm.py文件，可以看到第13行写着：

import cma.purecma as purecma

看来这个代码示例是利用cma工具包来操作的。在安装的requirements.txt文件中需要cma。
可在这儿查看cma的一些介绍：https://www.cnpython.com/pypi/cma （英文的在这儿：https://pypi.org/project/cma/）。

一种求解困难（非凸、病态、多模、粗糙、带噪声）连续优化问题的随机数值优化算法，使用Python实现。（ 参考 ）
其中pycma模块提供了独立的两个CMA-ES算法实现，即cma.CMAEvolutionStrategy类和cma.purecma.CMAES类。

再百度搜索一下CMA-ES（协方差矩阵自适应进化策略）会有很多的理论的教程。

• 对多元高斯分布进行采样得到新解，使用其中较好的解更新高斯分布的参数，最大熵原理。（参考）
• 这是一篇论文里面的CMA-ES的总体方案；第一步取值；第二步挑精英值；第三步，更新均值μ；第四步更新步长σ需要的演化路径；第五步更新步长；第六步更新path2 需要的C的演化路径；第七步我个人觉得不需要？；第八步更新协方差矩阵。（参考）

2.1.2  读代码

来看看RABBO中的示例代码。首先看blackbox_starter_kit/algorithms/submission/submission_algorithm.py ，与cma_algorithm.py相似。submission_algorithm.py中注释掉的代码部分，是4种功能的参考：

1. 提取问题相关接口和参数。
2. 实现一个循环采样的逻辑直到完成指定数量的采样并且消耗完指定的优化时长。
3. 在循环采样过程中使用相应api评估函数值。
4. 在循环采样过程中使用相应api评估约束违背情况。

图注：这个截图的文件里，我添加了cma算法代码进去，所以行号和长度不一样。

然后再看cma_algorithm.py。

• 后面有几个自定义函数def convert_to_real_config(self):、def get_bounds(self, real_config):、def convert_to_origin_param_dict(self, suggest_x):、def update_best(self, param, fval, directs):，做参数数据等的计算和更新工作。
• 算法迭代的内容在def optimize(self, prob_info, query_budget, time_budget):中。

这个示例代码里，用CMA库去做RABBO库中的问题求解的主要步骤如下：
（1） 首先，初始化参数，参数用处后面会体现。

    def __init__(self, sigma=0.3, penalty_coeff=10e8, penalty_thresh=10e-3, **kwargs):
        """Build wrapper class to use an optimizer in benchmark.

        Parameters
        ----------

        """
        AbstractAlgorithm.__init__(self)
        self._sigma = sigma
        self.name = "CMA"
        self.penalty_coeff = penalty_coeff
        self.penalty_thresh = penalty_thresh # 0

（2） 然后对输入的一番处理。
如_query_budget和_time_budget（次数和时间预算）用于后面迭代次数的控制，在examples/experiment_example.py中会用到。time_budget时间单位在optimize(）中注释的是秒。

（3） 创建一个CMAES的实例

        x0 = np.asarray(lb+(ub-lb)/2)    # initial point for CMA
        sigma = np.asarray(self._sigma)     #  standard deviation (step size)
        es = purecma.CMAES(x0, sigma)    # create CMAES instance

（4）开始迭代优化

• 最外层while循环是用_query_budget和_time_budget（次数和时间预算）控制的。

  • 每次通过es.ask()来获取一组推荐解X；es.tell(X, fitlist)来更新ES。 这里要注意CMA默认的算法direction是minimization的，因此fitlist赋值时候需要转换一下。

    • X循环内，

      • 采用 funcval = self._problem(suggest_x_ori) 来获得待优化问题的评价evaluate，
      • 采用violations, fg = self._constr_eval(suggest_x_ori) 来评估约束的违反情况，并根据惩罚penalty调整评价。
    • 如果没有违反约束（fg)，则把def update_best(self, param, fval, directs)结果更新。

（5）返回结果：return self.best_param, self.best_fval
核心部分的代码copy如下：

        ct = 0
        while self._query_budget - ct > 0 or self._time_budget - tm >0:    # loop until both query budget and time budget are used
            X = es.ask()  # deliver candidate solutions
            fitlist=[]
            for x in X:
                fg = True
                suggest_x_ori = self.convert_to_origin_param_dict(x)
                self.para_list.append(suggest_x_ori)
                funcval = self._problem(suggest_x_ori)
                self.func_list.append(funcval)
                aug_funcval = copy.deepcopy(funcval)

                if self._constraint_types != []:    
                    penalty = 0
                    violations, fg = self._constr_eval(suggest_x_ori)
                    self.constraint_viol.append(violations)
                    
                    for j in range(len(violations)):
                        if self._constraint_types[j] == 'lt':
                            penalty = max(violations[j], 0) * self.penalty_coeff
                        elif self._constraint_types[j] == 'eq':
                            penalty = abs(violations) * self.penalty_coeff
                        for k in range(len(funcval)):
                            if self._directions[k] == "minimize":
                                aug_funcval[k] += penalty
                            elif self._directions[k] == "maximize":
                                aug_funcval[k] -= penalty
                    self.aug_func_list.append(aug_funcval)

                if self._directions[0] == "maximize": # for now only consider single objective
                    fitlist.append(-aug_funcval[0])
                elif self._directions[0] == "minimize":
                    fitlist.append(aug_funcval[0])

                # if penalty > self.penalty_thresh * abs(funcval[k]): # use customized judgement
                if fg: # fg = True means feasible
                    self.update_best(suggest_x_ori, funcval, self._directions)
                else:
                    # print("The current point is infeasible.")
                    print(".") # yoyo add 

                ct += len(funcval)
            es.tell(X, fitlist)  # all the work is done here
            tm = time.time() - st
            
        if self.best_param == []:
            print("No feasible points have been found.")
            self.best_fval = []
        else:
            print("Best value found:\n\tf(x) = {}\nObserved at:\n\tx = {}".format(self.best_fval, self.best_param))
        return self.best_param, self.best_fval

2.2 Random Search（随机搜索）

2.2.1 Random Search 背景知识

RandomSearch的代码在RABBO的blackbox_starter_kit/algorithms/random_search/random_search_algorithm.py文件。
​

继续百度一下Random Search，又会有一堆材料。随机搜索，如字面意思，就是随机撒点来搜索。

• 可参考深度学习里面超参调优的random Search方法例子（参考）。
• 查询资料的时候就会看到GridSearch和RandomSearch等算法的科普。目前业界用得比较多的分别是Grid search、Random search和Bayesian Optimization，业界公认的Random search效果会比Grid search好。（参考）

2.2.2 读代码

继续读RABBO的random_search_algorithm.py文件。挺多代码内容是一样的。类似前面的内容分解代码：
（1） 首先，初始化参数，用处后面会体现。random search与前面的CMA不同的地方是sigma换成了n_suggestion，意思是采样点的个数，在后面的代码里可以看出用途。

    def __init__(self, n_suggestion=1, penalty_coeff=10e8, penalty_thresh=10e-3, **kwargs):
        """Build wrapper class to use an optimizer in benchmark.

        Parameters
        ----------

        """
        AbstractAlgorithm.__init__(self)
        self._n_suggestions = n_suggestion
        self.name = "Random search"
        self.penalty_coeff = penalty_coeff
        self.penalty_thresh = penalty_thresh # 0

（2）然后对输入的一番处理。类似之前CMA的代码。
（3）开始迭代优化，大部分内容与前面CMA代码解析相似，除了红色标记部分。

• 最外层while循环是用_query_budget和_time_budget（次数和时间预算）控制的。

  • _n_suggestions（多个个采样点）的循环内，默认是1个

    • 采用numpy生成随机推荐值cand_x = np.random.uniform(lb, ub, size=(1, len(lb)))
    • 采用 funcval = self._problem(suggest_x_ori) 来获得待优化问题的评价evaluate，
    • 采用violations, fg = self._constr_eval(suggest_x_ori) 来评估约束的违反情况，并根据惩罚penalty调整评价。
    • 如果没有违反约束（fg)，则把def update_best(self, param, fval, directs)结果更新。

（4）返回结果：return self.best_param, self.best_fval
核心部分的代码copy如下：

        while self._query_budget - ct > 0 or self._time_budget - tm >0:
            ct += self._n_suggestions
            for i in range(self._n_suggestions):
                fg = True
                cand_x = np.random.uniform(lb, ub, size=(1, len(lb)))
                suggest_x_ori = self.convert_to_origin_param_dict(cand_x.tolist()[0])
                self.para_list.append(suggest_x_ori)
                funcval = self._problem(suggest_x_ori)
                self.func_list.append(funcval)
                
                if self._constraint_types != []:
                    aug_funcval = copy.deepcopy(funcval)
                    penalty = 0
                    violations, fg = self._constr_eval(suggest_x_ori)
                    print('fg','feasible' if fg==True else 'infeasible')
                    self.cosntraint_viol.append(violations)
                    
                    for j in range(len(violations)):
                        if self._constraint_types[j] == 'lt':
                            penalty = max(violations[j], 0) * self.penalty_coeff
                        elif self._constraint_types[j] == 'eq':
                            penalty = abs(violations) * self.penalty_coeff
                        for k in range(len(funcval)):
                            if self._directions[k] == "minimize":
                                aug_funcval[k] += penalty
                            elif self._directions[k] == "maximize":
                                aug_funcval[k] -= penalty
                    self.aug_func_list.append(aug_funcval)


                # print(self.cosntraint_viol)
                if fg:
                    self.update_best(suggest_x_ori, funcval, self._directions)
                else:
                    print("The point is infeasible.")
            tm = time.time() - st
        if self.best_param == []:
            print("No feasible points have been found.")
            self.best_fval = []
        else:
            print("Best value found:\n\tf(x) = {}\nObserved at:\n\tx = {}".format(self.best_fval, self.best_param))      
        return self.best_param, self.best_fval

2.3 两种方法结果对比

2.3.1 肉眼对比

怎么评判算法的好坏呢？RABBO里面提供了个评测思路:对于单一问题，从采样次数和优化时间两个维度考察算法表现，分别记录：

1. 算法在给定采样次数内得到的最优解，
2. 算法在给定优化时间内得到的最优解。

也就是根据最后的目标值f(x) 和运算耗时来评估。

这里我们把RABBO自带的两种方法的输出结果整理一下，得到如下的结果对比表格。问题的f(x)值不是都是越大越好的，需要根据优化方向来定的，可以去对应问题的Python文件看到directions 。
跑实验给定的参数可在examples/experiment_example.py中修改，我用的是默认的（根据后面RABBO线上测评的结果看，是与线上测评的次数设置不一样的，大家根据需要修改）：

exp.run(n_exps=2, query_budget=200, time_budget=1.5)

​

表注：CMA 和 RandomSearch算法效果对比

从上表可以看出，目标值f(x)结果是CMA的结果更优秀呢，虽然在有些问题上耗时久一些。不过风场的这个Windfarm问题两个算法都解不出来，可以感觉到这两种方法都有些弱呢。

2.3.2 提交到RABBO线上测评

上一篇《【系列1/4】》中给大家介绍了怎么打包镜像评测，
最近RABBO的线上评测题有修改，测评的时间非常的长，要好几个小时。
根据评分规则，totalScore是越高越好，据说是不知道上限有多高的，因为有些问题没有人知道最优解是多少。windfarm那个问题没有解出来，拉低了平均分很多呀。smelting的问题只有在线上测评时候有，看起来也很难解呢…………

RABBO线上测评结果截图如下（结果里面的日志是有测评平台输出和我们自己编代码里的print输出）：

RABBO线上评测给出的成绩score整理成表格如下：

3 其他的算法

从前面的测评结果看，CMA（CMA-ES, 协方差矩阵自适应进化策略)要比random Search（随机搜索）要优秀，不过对于RABBO中给出的smelting和windfarm问题都不能好的求解，根据云栖大会MindOpt团队的视频分享（ https://yunqi.aliyun.com/2021/agenda/session170 直播回放），他们的求解器性能应该很不错的。还有更好的算法值得去发掘、学习和研究。

从前面的搜索中，大家应该能看到很多其他的优化算法。有一篇写的挺不错的：《贝叶斯优化: 一种更好的超参数调优方式》https://zhuanlan.zhihu.com/p/29779000

自动调参算法：说到自动调参算法，大家可能已经知道了Grid search（网格搜索）、Random search（随机搜索），还有Genetic algorithm（遗传算法）、Paticle Swarm Optimization（粒子群优化）、Bayesian Optimization（贝叶斯优化）、TPE、SMAC等。

像遗传算法和PSO这些经典黑盒优化算法，我归类为群体优化算法，也不是特别适合模型超参数调优场景，因为需要有足够多的初始样本点，并且优化效率不是特别高，本文也不再详细叙述。

目前业界用得比较多的分别是Grid search、Random search和Bayesian Optimization。

大家也可以多利用互联网搜索，多学习一下。
后面我也会找一个经典的方法改造一下参与测评，欢迎交流。

4 其他文

后面，立个flag，将补充更多文章，学会黑盒优化。
借助阿里达摩院的MindOpt的RABBO榜单系列：
1.《新手一步步学习黑盒优化算法，借助达摩院MindOpt的RABBO榜单【系列1/4】》
2.（本篇）《学习黑盒优化榜单RABBO提供的两个参考算法CMA和RandomSearch》

1. todo：《改编一个自己的黑盒优化算法或开源算法》
2. todo：《将我要解决的业务问题建模成RABBO的问题格式，然后去求解》