简介
LAPO,即闪电连接优化(Lightning Attachment Procedure Optimization),听名字就知道是受到了闪电的刺激,而获得灵感的一种算法。
为了走进LAPO的内心,于是上网搜了一下闪电的图片
呃不好意思,是下面这个
发现闪电连接无非是分岔而已,而且这个岔如果非常厚实,那么会继续分叉,否则就会迅速消失。
接下来可以思考,这个岔何以非常强壮?那肯定是这个位置比较适合分岔,在对的地方分岔了,闪电就会继续分岔,否则就面临着消失。
至此,闪电过程被抽象为分支的出现和消失。
原理
初始化
x i = rand ( x L , x R ) x^i = \operatorname{rand}(x_L, x_R) xi=rand(xL,xR)
x L , x R x_L, x_R xL,xR表示随机数的生成范围。
根据目标函数,可以得到 x i x^i xi的适应度
f i = f ( x i ) f^i = f(x^i) fi=f(xi)
下一跳
对于第 i i i个闪电点,周围有很多个可以分岔的位置,其第 j j j个位置记作 x j i x^i_j xji,计算所有可能位置的平均值,及其适应度
计算所有点的平均值,及其平均值的适应度
x ˉ i = 1 N ∑ j N x j i , f ˉ i = f ( x ˉ i ) \bar x^i=\frac{1}{N}\sum_j^Nx_j^i, \bar f^i=f(\bar x^i) xˉi=N1j∑Nxji,fˉi=f(xˉi)
如果 f j i f^i_j fji优于 f ˉ i \bar f^i fˉi,则闪电向这边移动,否则闪电向另一侧移动
x i ∗ = { x j i + ( x ˉ i + x j i ⋅ rand ) ⋅ rand f j i < f ˉ i x j i − ( x ˉ i + x j i ⋅ rand ) ⋅ rand f j i ⩾ f ˉ i x^{i*}=\left\{\begin{aligned} x^i_j+(\bar x^i + x^i_j\cdot\operatorname{rand})\cdot\operatorname{rand}\quad f^i_j<\bar f^i\\ x^i_j-(\bar x^i + x^i_j\cdot\operatorname{rand})\cdot\operatorname{rand}\quad f^i_j\geqslant\bar f^i\\ \end{aligned}\right. xi∗={xji+(xˉi+xji⋅rand)⋅randfji<fˉixji−(xˉi+xji⋅rand)⋅randfji⩾fˉi
分支消失
如果 x i ∗ x^{i*} xi∗的适应度比 x i x^{i} xi要好,那么就这个新点就保留,否则就任其消失。
地面接收
一般来说闪电肯定是要打在避雷针上的,最优解就相当于是避雷针。随着迭代的不断加深,即闪电不断第分岔,也就更加接近避雷针所在的位置。
在LAPO里,通过迭代次数来创建一个概率因子S,其表达式为
S = 1 − t t M exp − t t M S=1-\frac{t}{t_M}\exp-\frac{t}{t_M} S=1−tMtexp−tMt
随着迭代次数增加, S S S的值不断减小,意味着分支点的抖动逐渐降低,其作用在分支上的方式如下
x i ∗ = x i ∗ + rand ⋅ S ⋅ ( x m − x M ) x^{i*}=x^{i*}+\operatorname{rand}\cdot S\cdot(x_m-x_M) xi∗=xi∗+rand⋅S⋅(xm−xM)
其中, x m x_m xm和 x M x_M xM为最优和最差情况下的位置。
Python实现
由于LAPO算法不需要保留上一代闪电点,而且闪电点之间也没有什么关联,所以实现起来比较简单。
首先实现闪电的分支迭代过程
import numpy as np
rand = np.random.rand
# x为当前迭代点;func为迭代函数;N为预选分支点数目
# S为尖端放电系数
def branch(x, func, N, S):
fit = func(x) #x的适应度
# 预选点
stars = [x+rand(*x.shape) for _ in range(N)]
# 所有预选点的适应度
starFits = [func(star) for star in stars]
xBar = np.mean(stars, axis=0)
fBar = np.mean(starFits)
xs = []
for i in range(N):
sign = 1 if starFits[i] < fBar else -1
xs.append(x+sign*(xBar+stars[i]*rand())*rand())
# 上面已经给出了所有可能的分支
xs = np.array(xs)
fits = np.array([func(xi) for xi in xs])
ind = np.where(fits<fit)[0]
if len(ind)==0:
return [x]
xs,fits = xs[ind], fits[ind]
# 如果没有更好的结果,就保留现有结果
S = S*(xs[np.argmin(fits)] - xs[np.argmax(fits)])
return [x+rand()*S for x in xs]
然后实现主函数
from itertools import chain
# nInit为初始化点个数,nDim为数据维度,nBranch为分支点个数
# nIter为迭代次数
def lapo(nInit, nDim, nBranch, nIter, func):
# xs为点集
xs = [rand(nDim) for _ in range(nInit)]
for i in range(nIter):
# 尖端放电系数
S = 1 - (i/nIter)*np.exp(-i/nIter)
xs = [branch(x, func, nBranch, S) for x in xs]
xs = list(chain(*xs))
fits = [func(x) for x in xs]
msg = f"得到{len(xs)}组结果,其中最佳适应度为{np.min(fits)},"
msg += "xs=" + ", ".join([f"{xi}" for xi in xs[np.argmin(fits)]])
print(msg)
最后,找个函数测试一下,测试函数为
f ( x ⃗ ) = ∑ i = 0 cos ( i x i 5 ) ∗ ( i + 1 ) f(\vec x)=\sum_{i=0}\cos(\frac{ix_i}{5})*(i+1) f(x)=i=0∑cos(5ixi)∗(i+1)
def test(xs):
_sum = 0.0
for i in range(len(xs)):
_sum = _sum + np.cos((xs[i]*i)/5)*(i+1)
return _sum
if __name__ == "__main__":
lapo(10, 5, 3, 20, test)
效果为
>python lapo.py
得到1093组结果,其中最佳适应度为-12.570095759883985,xs=-11.260084641709877, -13.431104443084656, -7.471806128776576, -16.196185355184905, -11.894803311699398
当然,这个程序有一个bug,当新种群中没有更优情况的时候,上一代计算结果会被保存,随着迭代次数越来越多,非常容易内存爆炸,看来是要对每一代种群进行以此筛选才行。