求解广义Fermat-Torricelli问题的多层邻近梯度算法

马丽丽，谢秋玲，胡清洁

(桂林电子科技大学 数学与计算科学学院，广西 桂林 541004)

Fermat-Torricelli问题最初由法国数学家Pierre de Fermat提出，即给定平面上的3点，找到1个点，使其与3点的欧几里德距离之和最小。这个问题引发了对位置科学的研究兴趣，意大利数学家Evangelista·Torricelli最早研究该问题。随后这个问题被称为Fermat-Torricelli问题。该问题一直是优化领域中的一个研究热点，在优化网络[1]、设施选址[2]、计算几何[3]、定位科学、物流等方面具有广泛应用。l2范数下广义Fermat-Torricelli问题模型为

(1)

其中：ai∈Rn,i=1,2,…,m；Ω为Rn上的非空闭凸集。集合的广义Fermat-Torricelli问题模型为

(2)

其中：Ωi(i=1,2,…,m)为目标集。

近年来，Fermat-Torricelli问题得到广泛关注。Weiszfeld[4]首次提出求解Fermat-Torricelli问题的算法，随后，许多研究者对该算法进行了修正和改进[5-7]。Nam等[8]提出用Nesterov光滑技术求解广义Fermat-Torricelli问题，其主要思想是对原目标函数构造光滑逼近，利用加速梯度法求解，并用数值实验说明算法的有效性。Parpas等[9]基于多层优化方法，提出用多层加速梯度镜面下降算法求解大规模人脸识别问题，并用数值实验验证了该算法的可行性。

为求解广义Fermat-Torricelli问题，基于文献[9-10]，提出多层邻近梯度算法，并给出相关收敛速度分析和数值实验。

1 预备知识

定义1给定集合Ω⊆Rn，则

称为集合Ω的示性函数。

定义2[8]点x到Rn上非空闭凸集Ω的欧式投影为

π(x;Ω)={w∈Ω|d(x;Ω)=‖x-w‖}。

特别地，点x到以δ为中心、r为半径的立方体集合Ω上的投影为

PΩ(x)=max{δi-r,min{xi,δi+r}}。

2 算法描述

对模型(1)的目标函数进行μ光滑[11]和引入示性函数，则原问题可转化为

同样地，模型(2)的问题可转化为

将上述2个问题的一般形式记为

(3)

其中：f:Rn→R为光滑凸函数；g:Rn→R为连续凸函数不一定光滑。

QL(x,y)=f(x)+〈y-x,f(x)〉+

定义PL(x)=argmin{QL(x,y),y∈Rn}。

多层算法需要在层次之间使用适当的信息传递机制。为了在层次之间传递信息，引入线性限制算子R:Rn×nH和延拓算子P:RnH×n，其中nH

σP=RT，

其中σ>0，不失一般性，假设σ=1。算子的选取参见文献[12]。

为求解粗糙子问题，引入一个光滑参数ε>0，对g(x)进行光滑化得到gε(x)，光滑后的目标函数为

Fε(x)=f(x)+gε(x)。

粗糙模型的关键属性是初始点xH,0处2个模型的最优性条件相匹配，可通过在粗略模型的目标函数中引入线性项来实现。粗糙模型的目标函数为

FH,ε(xH)=fH(xH)+gH,ε(xH)+〈vH,xH〉。

其中：fH(xH)、gH,ε(xH)均为光滑函数；vH=RFε(x)-(fH(Rx)+gH,ε(Rx))。

迭代过程中，算法是否使用粗糙模型构建搜索方向取决于当前迭代点xk的最优性条件。进行粗糙迭代应满足以下2个条件：

(4)

为满足相容性，用xH,0=Rxk开始粗糙迭代，通过一阶优化方法求解粗糙模型并构造搜索方向，

dk(xk)=P(xH,NH-Rxk)，

其中xH,NH∈RnH为执行NH次迭代后粗糙模型的近似解。由于通常找不到精确解，利用Nesterov加速梯度算法[8]求解该子问题，得到一个可接受的解xH,NH。

多层邻近梯度算法的步骤为：

1)选取L>0，τ>0,ε>0,0

2)若粗糙条件(4)不满足，转步骤3)，否则执行NH次粗糙迭代得xH,NH，计算dk=P(xH,NH-Rxk)。由Armijo线搜索求得步长sk，

Fε(xk+skdk)≤Fε(xk)+cskFε(xk)Tdk。

3)计算yk=PL(xk)，若‖xk-yk‖<τ，迭代终止，输出yk作为近似极小点，否则转步骤4)。

3 收敛速度分析

为了分析多层邻近梯度算法的收敛速度，首先给出如下4个引理。

引理1若x∈Rn，存在L>0，使得

Fε(PL(x))≤Q(PL(x),x)

成立，则对∀y∈Rn，有

L〈x-y,PL(x)-x〉。

引理2设{xk}和{yk}是由多层邻近梯度算法产生的序列，则对∀k≥1，有

其中

vk=Fε(yk)-Fε(y*),uk=tkyk-(tk-1)yk-1-1。

2〈xk+1-yk,yk+1-xk+1〉。

(5)

2〈xk+1-y*,yk+1-xk+1〉。

(6)

‖tk+1(yk+1-xk+1)‖2+2tk+1×

〈tk+1xk+1-(tk+1-1)yk-y*,yk+1-xk+1〉=

‖tk+1yk+1-(tk+1-1)yk-y*‖2-

‖tk+1xk+1-(tk+1-1)yk-y*‖2=

‖tk+1yk+1-(tk+1-1)yk-y*‖2-

‖tkyk-(tk-1)yk-1-y*‖2。

由uk的定义得

引理得证。

引理3[10]设{ak}和{bk}都是正实数序列，假设对于∀k≥1，a1+b1≤c，c>0，ak-ak+1≥bk+1-bk成立，则对∀k≥1，ak≤c也成立。

引理4[10]设{tk}是由多层邻近梯度算法产生的序列，且t1=1，则对∀k≥1，tk≥(k+1)/2成立。

基于上述4个引理，估计该算法的收敛速度。

定理1设{yk}是由多层邻近梯度算法产生的序列，则对∀k≥1，a1+b1≤c，c>0，有

Fε(y*)-Fε(y1)=Fε(y*)-Fε(PL(y1))≥

等价于

‖y1-y*‖2-‖x1-y*‖2。

即a1+b1≤c成立，由引理2得ak-ak+1≥bk+1-bk成立，由引理3得

从而可得

定理得证。

4 数值实验

为检验算法的有效性，考虑文献[8]的数值算例，采用Matlab编制程序，其测试环境为Window 10操作系统，Intel (R) Core (TM) i5-6500 CPU @ 3.20 GHz，8.00 GB RAM。

算例1给定1217个美国城市的经纬度，其数据可从文献[8]提供的网址找到。为更符合实际分布情况，将其经度乘以-1由正转负。现在需要找到一个点，使得该点到给定美国城市点的距离之和最小。选取初始点y0=(20,-100)，限制算子R=[1/4 1/2]，粗糙参数η=0.2，θ=0.8，M=30，利普希茨常数L=70，精度τ=10-5，得到近似最优解y*≈(38.63,-97.35)。美国各城市及其最优点的分布结果如图1所示。将多层邻近梯度算法(MPGA)与迭代收缩阈值算法(FISTA)算法进行比较，结果如表1所示。从表1可看出，在相同的初始值和最优解达到相同精度时，多层邻近梯度算法迭代次数和运行时间较少。

图1 美国城市及其最优点的分布结果图

算法迭代数运行时间/s近似最优解y*MPGA161.43(38.63,-97.35)FISTA222.21(38.63,-97.35)

算例2与算例1的设置相同，以上述1217个美国城市为中心点，半径r=1.5的1217个正方形，现在需要在直线x-y=180上找到一个点，使得该点到1217个正方形的距离之和最小。选取初始点y0=(0,180)，其他参数的选取同算例1，最后得到近似最优解y*≈(56.84,-123.16)，各广场及其最优点的分布结果如图2所示。将多层临近梯度算法(MPGA)与迭代收缩阈值算法(FISTA)进行比较，结果如表2所示。从表2可看出，在相同的初始值和最优解达到相同精度时，多层邻近梯度算法迭代次数较少，但运行时间略多于FISTA算法。

图2 美国广场及其最优点的分布结果图

算法迭代数运行时间/s近似最优解y*MPGA381.51(56.84,-123.16)FISTA401.36(56.84,-123.16)

算例3考虑以(-6,6,-4)，(-5,-3,-6)，(2,3,4)，(4,-4,-5)，(5,6,-6)，(-5,-2,4)六个点为中心，半径r=1.5的6个正方体，现在需要找到一个点，使得该点到6个正方体的距离之和最小。选取初始点y0=(-6,6,2)，限制算子R=[1/4 1/2 1/4]，粗糙参数和精度的选取与算例1相同，利普希茨常数L=0.8，得到近似最优解y*≈(-1.040 5,0.840 2,-1.432 2)，正方体及其最优点的分布结果如图3所示。将多层邻近梯度算法(MPGA)与迭代收缩阈值算法(FISTA)进行比较，结果如表3所示。从表3可看出，在相同的初始值和最优解达到相同精度时，多层邻近梯度算法迭代次数和运行时间较少。

图3 正方体及其最优点的分布结果图

算法迭代数运行时间/s近似最优解y*MPGA110.027 6(-1.040 5,0.840 2,-1.432 2)FISTA160.028 1(-1.040 5,0.802 4,-1.432 2)

5 结束语

针对点和集合的广义Fermat-Torricelli问题，提出一种多层邻近梯度算法，采用μ光滑方法将原始目标函数转化为光滑函数，再利用多层邻近梯度算法求解，并从理论上证明该算法具有O(1/k2)的收敛速度。数值实验表明，多层邻近梯度算法求解广义Fermat-Torricelli问题是有效的。