力扣1515-服务中心的最佳位置

字数统计: 3.3k字 | 阅读时长: 16分

2020-12-03

摘要: 力扣1515：三分、梯度下降、爬山法、模拟退火

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

此前我们学习了一些优化算法，包括梯度下降、爬山法、模拟退火，都是用的力扣 1515 作为例题，三分查找也是解决力扣 1515 的经典算法。

力扣 1515 是第 197 场周赛 D 题，问题非常经典，解法很多，本文做个总结，首先介绍几何中位数的概念，对本题的目标函数是否为凸函数做个证明，然后罗列一下上面提到的几种解法，此外，用 scipy.optimize.minimaze 也是可以的。

$1 题目

1515. 服务中心的最佳位置

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

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

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

提示：

1
2
3

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

示例 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

示例 3：
输入：positions = [[1,1]]
输出：0.00000

示例 4：
输入：positions = [[1,1],[0,0],[2,0]]
输出：2.73205
解释：乍一看，你可能会将中心定在 [1, 0] 并期待能够得到最小总和，但是如果选址在 [1, 0] 距离总和为 3
如果将位置选在 [1.0, 0.5773502711] ，距离总和将会变为 2.73205
当心精度问题！

示例 5：
输入：positions = [[0,1],[3,2],[4,5],[7,6],[8,9],[11,1],[2,12]]
输出：32.94036
解释：你可以用 [4.3460852395, 4.9813795505] 作为新中心的位置

$2 题解

下面的几种解法中都用到了二维向量的代码模板，参考文章：向量的实现。

几何中位数

求出的 $(x, y)$ 是 ${(x_{i}, y_{i})}$ 的几何中位数。几何中位数没有解析解，无法求全局最小值，但可以证明给定的函数是凸函数，可以用求解局部最小值的方法得到全局最小值。

凸函数证明

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

考察 $f_{i}(x, y) = \sqrt{(x - x_{i})^{2} + (y - y_{i})^{2}}$，$f_{i}$ 连续且可导，一阶导数：

$\left\{ \begin{array}{**lr**} \frac{\partial f_{i}}{\partial x} = \frac{x - x_{i}}{\sqrt{(x - x_{i})^{2} + (y - y_{i})^{2}}}\\ \frac{\partial f_{i}}{\partial y} = \frac{y - y_{i}}{\sqrt{(x - x_{i})^{2} + (y - y_{i})^{2}}}\\ \end{array} \right.$

二阶导数：

$\left\{ \begin{array}{**lr**} \frac{\partial^{2} f_{i}}{\partial x^{2}} = \frac{(y - y_{i})^{2}}{((x - x_{i})^{2} + (y - y_{i})^{2})^{\frac{3}{2}}}\\ \frac{\partial^{2} f_{i}}{\partial y^{2}} = \frac{(x - x_{i})^{2}}{((x - x_{i})^{2} + (y - y_{i})^{2})^{\frac{3}{2}}}\\ \frac{\partial^{2} f_{i}}{\partial x\partial y} = \frac{\partial^{2} f_{i}}{\partial y\partial x} = -\frac{(x - x_{i})(y - y_{i})}{((x - x_{i})^{2} + (y - y_{i})^{2})^{\frac{3}{2}}}\\ \end{array} \right.$

$f_{i}$ 对应的 Hessian 矩阵

$H(f_{i}) = \begin{bmatrix} \frac{\partial^{2} f_{i}}{\partial x^{2}} & \frac{\partial^{2} f_{i}}{\partial x\partial y} \\ \frac{\partial^{2} f_{i}}{\partial y\partial x} & \frac{\partial^{2} f_{i}}{\partial y^{2}} \\ \end{bmatrix}$

$H(f_{i})$ 的主子式均大于等于 0

$\left\{ \begin{array}{**lr**} \frac{\partial^{2} f_{i}}{\partial x^{2}} \geq 0 \\ \frac{\partial^{2} f_{i}}{\partial y^{2}} \geq 0 \\ \frac{\partial^{2} f_{i}}{\partial x^{2}}\frac{\partial^{2} f_{i}}{\partial y^{2}} - \frac{\partial^{2} f_{i}}{\partial x\partial y}\frac{\partial^{2} f_{i}}{\partial y\partial x} = 0 \\ \end{array} \right.$

$f_{i}$ 的 Hessian 矩阵半正定，所以 $f_{i}$ 凸函数，凸函数的和还是凸函数，因此 $f$ 是凸函数.

算法1: 实数三分搜索

参考三分查找。

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

要求 $f(x, y)$ 的局部最小值，可以用寻找下降方向的梯度下降或者爬山法。

也三分查找，逐步缩小范围，直到找到最小值点。

给定的函数是二元的，外层的三分确定 x 的范围，分别固定 x 为 xm1, xm2 后，内层的三分进行最小值查找，将所得的最小值进行比较，进而排除对外层三分的某个范围 [xL, xm1], [xm2, xR]。

代码 C++ (模板：二元单峰函数、三分套三分)

struct Function
{
    // 二元函数
    vector<Vector2> points;
    Function(vector<vector<int>>& p)
    {
        int n = p.size();
        points = vector<Vector2>(n);
        for(int i = 0; i < n; ++i)
        {
            points[i].y = p[i][1];
            points[i].x = p[i][0];
        }
    }
    double operator()(double x, double y) const
    {
        double ans = 0;
        Vector2 p0(x, y);
        for(const Vector2 &p: points)
            ans += (p0 - p).norm();
        return ans;
    }
};

struct PartialFunction
{
    Function f; // 二元函数
    double x;
    PartialFunction(Function& f, double x):f(f),x(x){}
    double operator()(double y) const
    {
        return f(x, y);
    }
};

double ternary_search(PartialFunction& f, double L, double R)
{
    double l = L, r = R;
    while(l + EPS < r)
    {
        double m1 = l + (r - l) / 2.0;
        double m2 = m1 + (r - m1) / 2.0;
        if(f(m1) < f(m2))
            r = m2;
        else
            l = m1;
    }
    return l;
}

double ternary_search(Function& f, double Lx, double Rx, double Ly, double Ry)
{
    double lx = Lx, rx = Rx, ly = Ly, ry = Ry;
    double y_max = 0.0;
    while(lx + EPS < rx)
    {
        double xm1 = lx + (rx - lx) / 2.0;
        PartialFunction f1(f, xm1);
        double ym1_max = ternary_search(f1, ly, ry);
        double xm2 = xm1 + (rx - xm1) / 2.0;
        PartialFunction f2(f, xm2);
        double ym2_max = ternary_search(f2, ly, ry);
        if(f(xm1, ym1_max) < f(xm2, ym2_max))
        {
            y_max = ym1_max;
            rx = xm2;
        }
        else
        {
            y_max = ym2_max;
            lx = xm1;
        }
    }
    return f(lx, y_max);
}

class Solution {
public:
    double getMinDistSum(vector<vector<int>>& positions) {
        Function f(positions);
        int Lx = positions[0][0], Rx = positions[0][0];
        int Ly = positions[0][1], Ry = positions[0][1];
        for(vector<int>& p: positions)
        {
            Lx = min(Lx, p[0]);
            Rx = max(Rx, p[0]);
            Ly = min(Ly, p[1]);
            Ry = max(Ry, p[1]);
        }
        if(Lx == Rx && Ly == Ry)
            return 0.0;
        if(Lx == Rx)
        {
            PartialFunction pf(f, Lx);
            return ternary_search(pf, Ly, Ry);
        }
        return ternary_search(f, Lx, Rx, Ly, Ry);
    }
};

算法2: 梯度下降

参考梯度下降。

选定初始点 $\vec{x} = (x, y)$
迭代 $\vec{x}$，直到两次迭代之间的 $f(\vec{x})$ 差值足够小，例如 EPS=1e-9
- 计算梯度 $\nabla$
- 更新 $\vec{x} \leftarrow \vec{x} - \gamma \cdot \nabla$

代码 (C++) 小批量梯度下降

参数设置：

初始点 : 所有带你的算术平均值
学习率 : $\gamma = 1$
学习率衰减 : $\eta = 1e-3$
EPS : 1e-9
batch_size : 16

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;
    }
};

struct Gradient
{
    // 二元函数的梯度
    vector<Vector2> *points;
    Gradient(vector<Vector2>* p):points(p){}
    Vector2 operator()(const Vector2& p0, vector<int>& batch_data_idx) const
    {
        Vector2 nabla;
        for(int idx: batch_data_idx)
        {
            Vector2 p = (*points)[idx];
            nabla.x += (p0.x - p.x) / ((p0 - p).norm() + EPS);
            nabla.y += (p0.y - p.y) / ((p0 - p).norm() + EPS);
        }
        return nabla;
    }
};

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]);
        }
        Gradient grad(&data);
        Function f(&data);
        vec = vec * (1 / (double)n);
        double gamma = 1.0; // 学习率
        double decay = 1e-3; // 学习率衰减
        int batch_size = 16;
        std::default_random_engine dre(6);
        vector<int> data_idx(n);
        for(int i = 0; i < n; ++i)
            data_idx[i] = i;
        Vector2 iter_vec = vec;
        vector<int> batch_data_idx;
        while(true)
        {
            shuffle(data_idx.begin(), data_idx.end(), dre);
            for(int i = 0; i < n; i += batch_size)
            {
                int j = min(i + batch_size, n);
                batch_data_idx.assign(j - i, 0);
                for(int k = i; k < j; ++k)
                    batch_data_idx[k - i] = data_idx[k];
                Vector2 g = grad(iter_vec, batch_data_idx);
                iter_vec = iter_vec - g * gamma; // 更新
                gamma *= (1.0 - decay); // 每一轮学习率衰减 1 次
            }
            if(vec == iter_vec)
                break;
            vec = iter_vec;
        }
        return f(vec);
    }
};

算法3: 爬山法

参考: 爬山法。

当凸函数很难求导的时候，可以考虑用爬山法找最值。

对于二元函数，梯度可以分为 x 和 y 两个方向：$\nabla = (\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y})$

当前点在 (x, y)，下一步考虑四个方向，每个方向向前移动一点，移动的长度记为 step。

1 2	dx = {0, 1, 0, -1} dy = {1, 0, -1, 0}

这样 (x + step * dx + y + step * dy) 就是一个候选点，如果该点的函数值更小，就走到这个点，这样就可以像梯度下降一样得到本轮迭代的下一个点。

如果候选点的函数值没有更小，就考察下一个方向，如果四个方向都没有找到函数值更小的点，就衰减步长。

当步长小于 EPS 时，搜索结束。

代码 (C++)

参数设置：

初始点 : 所有带你的算术平均值
步长：step = 1
步长衰减 : $\eta = 5e-1$
EPS : 1e-9

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; // 步长衰减
        int dx[4] = {0, 1, 0, -1};
        int dy[4] = {1, 0, -1, 0};
        Vector2 nxt_vec;
        while(step > EPS)
        {
            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];
                if(f(nxt_vec) < f(vec))
                {
                    vec = nxt_vec;
                    flag = true;
                    break;
                }
            }
            if(!flag)
                step *= (1 - decay);
        }
        return f(vec);
    }
};

算法4: 模拟退火

在爬山法的基础上，迭代可行解时，以一定概率接受一个比当前解差的解，因此可能会跳出局部最优解。

注：本题实际没有局部最优的问题。

参考: 模拟退火。

代码 (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);
    }
};

算法5: `scipy.optimize.minimize`

scipy.optimize.minimize 是 Python 里常用的非线性规划求极值的方法。

非线性规划是指目标函数或约束条件中有一个或几个非线性函数的最优化问题的方法

minimaze 的使用方法如下：

1	minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)

其中各个参数的简要含义如下：

fun：目标函数，返回单值，
x0：初始迭代值，
args：要输入到目标函数中的参数。
method：求解的算法，例如 BFGS、SLSQP 等
jac：目标函数的雅可比矩阵。
hess，hessp：目标函数的 Hessian（二阶导数矩阵）或目标函数的 Hessian 乘以任意向量 p。
bounds：变量的边界。
constraints：约束条件（只对 COBYLA 和 SLSQP）。dict 类型。
    type : str， ‘eq’ 表示等于0，‘ineq’ 表示不小于0
    fun : 定义约束的目标函数
    jac : 函数的雅可比矩阵 (只用于 SLSQP)，可选项。
    args : fun 和 雅可比矩阵的入参，可选项。
tol：迭代停止的精度。
callback(xk)：每次迭代要回调的函数，需要有参数 xk
options：其他选项
    maxiter :  最大迭代次数
    disp :  是否显示过程信息

代码 (Python)

from scipy.optimize import minimize
class Solution:
    def getMinDistSum(self, positions: List[List[int]]) -> float:
        def getSum(x):
            return sum([sqrt((x[0] - x1) ** 2 + (x[1] - y1) ** 2) for x1, y1 in positions])
        return minimize(getSum, [50, 50]).fun

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

数优化

力扣1515-服务中心的最佳位置

$1 题目

$2 题解

几何中位数

凸函数证明

算法1: 实数三分搜索

代码 C++ (模板：二元单峰函数、三分套三分)

算法2: 梯度下降

代码 (C++) 小批量梯度下降

算法3: 爬山法

代码 (C++)

算法4: 模拟退火

代码 (C++)

算法5: `scipy.optimize.minimize`

代码 (Python)

数优化

$1 题目

$2 题解

几何中位数

凸函数证明

算法1: 实数三分搜索

代码 C++ (模板：二元单峰函数、三分套三分)

算法2: 梯度下降

代码 (C++) 小批量梯度下降

算法3: 爬山法

代码 (C++)

算法4: 模拟退火

代码 (C++)

算法5: scipy.optimize.minimize

代码 (Python)

算法5: `scipy.optimize.minimize`