模拟退火

字数统计: 1.8k字 | 阅读时长: 7分

2020-12-03

leetcode367, 代码模板119, 优化14, 搜索50, 物理8, 随机算法20

摘要: 介绍模拟退火的原理并解决力扣 1515

【对数据分析、人工智能、金融科技、风控服务感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：潮汐朝夕
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

各位好，在文章梯度下降和爬山法中我们学习了解决在目标函数上找最优解的两类算法，并分别解决了力扣 1515。

但是这两种算法只能解决目标函数为凸函数的情况，其中目标函数如果容易求导，可以用梯度下降；如果不容易求导，可以尝试爬山法。

如果目标函数非凸，可以尝试模拟退火算法。本文首先介绍模拟退火解决局部最优问题的思想，然后从热力学的角度理解，给出模拟退火的算法，最后依然以力扣 1515 为例，实现模拟退火算法，可以作为代码模板。

$1 背景

考虑优化问题：

$\min\limits_{x \in X}f(x)$

如果 $X$ 为离散有限值，可以通过搜索算法配合相关的优化手段得的问题的最优解。

如果 $X$ 连续，且 $f(x)$ 是凸的，可以通过梯度下降，爬山法获得最优解。

如果 $X$ 连续，且 $f(x)$ 非凸，可以也可以通过爬山法等近似算法找到解，但由于局部最优的问题，解基本上不会是全局最优。

$2 模拟退火

考虑以下优化问题，$X$ 连续，且 $f(x)$ 非凸：

$\min\limits_{x \in X}f(x)$

使用爬山法等贪心搜索算法可以得到解，会面临的局部最优问题。模拟退火是应对这种局部最优问题的解法之一。

从搜索空间内的某个点开始，每一步先选择一个邻居，在计算从现有位置到达邻居的概率，若概率大于阈值，则跳转到邻居，若概率较小，则停留在原位置不动。

模拟退火在迭代可行解时，以一定概率接受一个比当前解差的解，因此可能会跳出局部最优解。

热力学原理

将搜索空间中的点视为空气分子，这个点也可以像空气分子一样带有能量。由热力学的原理，温度为Ｔ时，出现能量差为 $dE$ 的概率为 $p(dE)$，对这块原理感兴趣的可以去看玻尔兹曼分布相关的内容。

$p(dE) = exp(\frac{dE}{kT})$

$k$ 为玻尔兹曼常数，$k = 1.3806488\times 10^{-23}$，$dE < 0$，$0 < p(dE) < 1$。温度越高，出现一次能量差为 $dE$ 的降温的概率就越大。

算法 (最小值优化)

在搜索问题中，当前状态为 $vec$，下一个状态为 $newvec$，能量差为 $\Delta f = fabs(f(newvec) - f(vec))$。当新状态更好时(例如最小值优化时, f(newvec) < f(vec))，直接跳转到 newvec；当新状态没有更好时，跳转到这个没有比当前解更好的解的概率:

$P(\Delta f) = exp(-\frac{\Delta f}{T}) \\$

初始化：初始温度 T = INF，温度下限 T = EPS，初始状态 vec，每个 T 的迭代次数 C
迭代 vec 直到温度降为 Tmin:
    产生新状态 new_vec，增量为 delta = f(new_vec) - f(vec)
    delta < 0，即邻居比当前点更小，则取 new_vec 为新状态;
    delta >= 0，即概率 exp(-delta/T) 取 new_vec 作为新状态
    check: 检查某些条件，若 check 通过，则改变某些其它参数，例如步长
    温度 T 迭代次数达到 C 时，降温
当温度将为 Tmin 后搜索结束

缺点

参数难调：

初始温度高，搜到全局最优的概率大，但搜索时间变长
退火速度，同一温度下迭代次数越大，搜索越充分，但搜索时间变长

$3 题目

1515. 服务中心的最佳位置

一家快递公司希望在新城市建立新的服务中心。公司统计了该城市所有客户在二维地图上的坐标，并希望能够以此为依据为新的服务中心选址：使服务中心到所有客户的欧几里得距离的总和最小。

给你一个数组 positions ，其中 positions[i] = [xi, yi] 表示第 i 个客户在二维地图上的位置，返回到所有客户的欧几里得距离的最小总和。

换句话说，请你为服务中心选址，该位置的坐标 [xcentre, ycentre] 需要使下面的公式取到最小值：

$\sum\limits_{i=0}\limits^{n-1}\sqrt{(x_{centre} - x_{i})^{2} + (y_{centre} - y_{i})^{2}}$

与真实值误差在 1e-5 之内的答案将被视作正确答案。

示例 1：
输入：positions = [[0,1],[1,0],[1,2],[2,1]]
输出：4.00000
解释：如图所示，你可以选 [xcentre, ycentre] = [1, 1] 作为新中心的位置，这样一来到每个客户的距离就都是 1，所有距离之和为 4 ，这也是可以找到的最小值。

示例 2：
输入：positions = [[1,1],[3,3]]
输出：2.82843
解释：欧几里得距离可能的最小总和为 sqrt(2) + sqrt(2) = 2.82843

提示：

1
2
3

1 <= positions.length <= 50
positions[i].length == 2
0 <= xi, yi <= 100

算法：模拟退火

这里使用模拟退火，具体的算法流程就是前面介绍的，下面的代码是对该算法的实现，可以作为代码模板。

代码 (C++)

参数设置：

double step = 1.0; // 步长
double decay = 5e-1; // 步长衰减
double T0 = 5e-1;
double min_T = 1e-8;
double T_decay = 0.9; // 温度衰减
double T = T0;
int C = 10; // 每个温度的迭代次数

struct Function
{
    // 二元函数
    vector<Vector2> *points;
    Function(vector<Vector2>* p):points(p){}
    double operator()(Vector2& p0) const
    {
        double ans = 0;
        for(Vector2 p: *points)
            ans += (p0 - p).norm();
        return ans;
    }
};

class Solution {
public:
    double getMinDistSum(vector<vector<int>>& positions) {
        int n = positions.size();
        Vector2 vec; // 初始值
        vector<Vector2> data; // 训练数据
        for(vector<int>& p: positions)
        {
            vec.x += p[0];
            vec.y += p[1];
            data.emplace_back(p[0], p[1]);
        }
        Function f(&data);
        vec = vec * (1 / (double)n);
        double step = 1.0; // 步长
        double decay = 5e-1; // 步长衰减
        double T0 = 5e-1;
        double min_T = 1e-8;
        double T_decay = 0.9; // 温度衰减
        double T = T0;
        int C = 10; // 每个温度的迭代次数
        int dx[4] = {0, 1, 0, -1};
        int dy[4] = {1, 0, -1, 0};
        Vector2 nxt_vec;
        std::default_random_engine dre;
        std::uniform_real_distribution<double> dr;
        int c = 0;
        while(T > min_T)
        {
            bool flag = false;
            for(int d = 0; d < 4; ++d)
            {
                nxt_vec.x = vec.x + step * dx[d];
                nxt_vec.y = vec.y + step * dy[d];
                double delta = f(nxt_vec) - f(vec);
                if(delta < -EPS)
                {
                    vec = nxt_vec;
                    flag = true;
                    break;
                }
                else
                {
                    double r = dr(dre);
                    if(r < exp(-fabs(delta) / T))
                    {
                        vec = nxt_vec;
                        flag = true;
                        break;
                    }
                }
            }
            if(!flag && step > EPS)
                step *= (1 - decay);
            if(++c == C)
            {
                c = 0;
                T = T * T_decay;
            }
        }
        return f(vec);
    }
};

潮汐朝夕的生活实验室 \Doge 陪伴一个算法工程师的职业生涯

数优化