优化算法实战
在遗传算法详解中我们介绍了爬山算法，退火算法和遗传算法的原理。在进行代码实战前让我们先复习下这几种算法。
爬山算法每次在上次解的基础上进行改进，所以最终得到的解极可能是个局部最优解。而退火算法则是在爬山算法的基础上进行一定的改进——爬山算法每次结果都是向着更好地结果迈进但是容易陷入一个局部最优解，而退火算法则是引入了一定的几率可以选择一个较当前更差的解以跳出局部最优解来寻找全局最优。爬山算法和退火算法每次都是优化一个解，而遗传算法则是每次优化一组解——淘汰当前种群中的一组较差的解而保留较好的那些解，然后基于这些最优的解来生成其余的解以保证种群的大小，由于种群的初始值最初都是都是随机化得到的，因此最后得到的这些解可能相似性较小。想象一下，由于各个解有较大差别，因此最终得到的最优种群中各个解可能分布在各个局部最优解中，这样某个解可能就正好是全局最优解。
下面通过一个例子实战更好地理解这些算法。
问题描述
现有10个学生需要安排到5个宿舍中，每个宿舍有两个隔间各可以住一个人。每个学生都有自己最想去的宿舍和次选的。如下：
# 宿舍，每个宿舍两个可用的隔间
dorms = ['Zeus', 'Athena', 'Hercules', 'Bacchus', 'Pluto']
# 代表学生的首选和次选
prefs = [('Toby', ('Bacchus', 'Hercules')),
('Steve', ('Zeus', 'Pluto')),
('Andrea', ('Athena', 'Zeus')),
('Sarah', ('Zeus', 'Pluto')),
('Dave', ('Athena', 'Bacchus')),
('Jeff', ('Hercules', 'Pluto')),
('Fred', ('Pluto', 'Athena')),
('Suzie', ('Bacchus', 'Hercules')),
('Laura', ('Bacchus', 'Hercules')),
('Neil', ('Hercules', 'Athena'))]
# 若学生排队选择宿舍，则各个学生可供选择的隔间
domain = [(0, (len(dorms)*2)-i-1) for i in range(len(dorms)*2)]
domain
# 输出，即第一个学生有10个床位可以选择，第二个学生9个，第三个。。。
[(0, 9),
(0, 8),
(0, 7),
(0, 6),
(0, 5),
(0, 4),
(0, 3),
(0, 2),
(0, 1),
(0, 0)]
现在我们随便选择一个解，即每个学生都选择剩余床位中的第一个，解如下：
def print_solution(vec):
'''给学生分别安排宿舍
'''
slots = []
# 为每个宿舍创建两个槽
# slots=[0,0,1,1,2,2,3,3,4,4] => ['Zeus','Zeus','Athena','Athena','Hercules','Hercules','Bacchus','Bacchus','Pluto','Pluto']
for i in range(len(dorms)):
slots+=[i, i]
# 遍历每一名学生的安置情况
for i in range(len(vec)):
x = int(vec[i])
# 从剩余槽中选择
dorm = dorms[slots[x]]
# 输出学生及其被分配的宿舍
print(prefs[i][0], dorm)
# 删除已被选用的槽
del slots[x]
# 每个学生都选择余下的第一个宿舍
vec = [0]*len(prefs)
print_solution(vec)
# 输出
Toby Zeus
Steve Zeus
Andrea Athena
Sarah Athena
Dave Hercules
Jeff Hercules
Fred Bacchus
Suzie Bacchus
Laura Pluto
Neil Pluto
在上面我们看到了一个随意的分配结果。在下面我们将通过不同的算法来找到更优的解。不过在介绍各种算法之前我们都需要先定义一个用来比较不同解优劣的方法。
成本函数
在成本函数中，我们对不同的分配给予不同的惩罚：若一个学生的分配结果正好是其首选的则罚0分，若是次选的则罚1分，而不再其选择之中则需要罚3分。
def dormcost(vec):
'''每个学生选择不同的宿舍时候的成本
'''
cost = 0
# 建立一个槽序列
slots = [0, 0, 1, 1, 2, 2, 3, 3, 4, 4]
# 遍历每一个学生
for i in range(len(vec)):
x = int(vec[i])
dorm = dorms[slots[x]]
# pref => 可供学生选择的首选和次选
pref = prefs[i][1]
# 首选成本为0，次选成本为1
if pref[0] == dorm:
cost += 0
elif pref[1] == dorm:
cost += 1
else: # 不在选择之列成本值为3
cost += 3
# 删除已被占用的槽
del slots[x]
return cost
既然定义好了用来比较不同解的成本函数，那么我们就可以直接对不同解进行比较了。
随机算法
在实现上面的3个算法之前不妨先看看一个仅生成随机解的算法能达到什么样的水平。
def random_optimize(domain, costf):
'''随机优化算法, 每次随机选择一个解
Args:
domain: 内嵌二元组的列表，指定了每次变量可选的最大值和最小值区间
costf: 成本函数
'''
best = 999999
bestr = None
# 创建1000次随机解，从中选择最好的解
for i in range(1000):
# 生成一个随机解
r = [random.randint(domain[j][0], domain[j][1])
for j in range(len(domain))]
# 得到成本
cost = costf(r)
# 与到目前为止的最优解进行比较
if cost < best:
best = cost
bestr = r
return r
# 运行程序10次，查看每次得到的最优解
for i in range(10):
rs = random_optimize(domain, dormcost)
print(dormcost(rs))
# 输出
23
16
12
26
21
22
22
22
20
23
# 学生全选0的解
vec = [0]*len(prefs)
print(dormcost(vec))
16
从上面可以看到，随机算法得到的解大部分时候还不如学生全选0的解。
爬山算法
def hill_climb(domain, costf):
'''爬山算法
Args:
domain: 内嵌二元组的列表，指定了每次变量可选的最大值和最小值区间
costf: 成本函数
'''
# 先创建一个随机解
sol = [random.randint(domain[j][0], domain[j][1])
for j in range(len(domain))]
# 主循环
while True:
# 创建相邻的列表，也当前位置各个维度小爬一步后的解组成的列表
neighbors = []
for j in range(len(domain)):
# domain中的每个值相当于一个维度（方向），在每个方向上相对于原值偏离一点，生成一个相邻位置的解
if sol[j] > domain[j][0]:
# 大的值变小一点，即 sol[j] = sol[j]-1
neighbors.append(sol[0:j]+[sol[j]-1]+sol[j+1:])
if sol[j] < domain[j][0]:
# 小的值变大一点，即 sol[j] = sol[j]+1
neighbors.append(sol[0:j]+[sol[j]+1]+sol[j+1:])
# 在相邻解中寻找最优解
current = costf(sol)
best = current
for j in range(len(neighbors)):
# 计算在每个相邻解的成本
cost = costf(neighbors[j])
if cost < best:
best = cost
sol = neighbors[j]
# 如果没有更好地解，则退出循环
if best == current:
break
return sol
# 运行程序10次，查看每次得到的最优解
for i in range(10):
hs = hill_climb(domain, dormcost)
print(dormcost(hs))
# 输出
13
16
15
16
17
15
13
12
17
21
可以看到，爬山算法相较于随机算法已经有了较大的改进了。
退火算法
def annealing_optimize(domain, costf, T=1000.0, cool=0.95, step=1):
'''模拟退火算法
Args:
domain: 内嵌二元组的列表，指定了每次变量可选的最大值和最小值区间
costf: 成本函数
T: 初始温度
cool: 降温幅度
step: 每一步走的大小
'''
# 随机初始化值
vec = [float(random.randint(domain[i][0], domain[i][1]))
for i in range(len(domain))]
while T > 0.1:
# 选择一个索引值
i = random.randint(0, len(domain)-1)
# 选择改变索引值的方向
direct = random.randint(-step, step)
# 创建一个代表题解的新列表，改变其中一个值
vecb = vec[:]
vecb[i] += direct
# 若取值超出规定范围，将其拉回来
if vecb[i] < domain[i][0]:
vecb[i] = domain[i][0]
elif vecb[i] > domain[i][1]:
vecb[i] = domain[i][1]
# 计算当前成本和新的成本
ea = costf(vec)
eb = costf(vecb)
# 它是更好地解码？或者是趋向最优解的可能的临界值吗？
if (eb < ea or random.random() < pow(math.e, -(eb-ea)/T)):
vec = vecb
# 降低温度
T = T*cool
return vec
# 运行程序10次，查看每次得到的最优解
for i in range(10):
ass = annealing_optimize(domain, dormcost)
print(dormcost(ass))
# 输出
7
13
9
8
7
13
7
13
8
17
相较于爬山算法，退火算法结果更进一步，更加优秀。
遗传算法
def genetic_optimize(domain, costf, popsize=50, step=1,
mut_prob=0.2, elite=0.2, maxiter=100):
'''遗传算法
Args:
domain: 内嵌二元组的列表，指定了每次变量可选的最大值和最小值区间
costf: 成本函数
popsize: 种群大小
step: 每一步走的大小
mut_prob: 突变概率
elite: 保留的精英数量，种群中被认为是优解且被允许传入下一代的部分
maxiter: 繁殖多少代
'''
def mutate(vec):
'''变异操作(基因突变)
'''
i = random.randint(0, len(domain)-1)
if random.random() < 0.5 and vec[i] > domain[i][0]:
return vec[0:i]+[vec[i]-step]+vec[i+1:]
elif vec[i] < domain[i][1]:
return vec[0:i]+[vec[i]+step]+vec[i+1:]
def crossover(r1, r2):
'''交叉操作(基因重组)
'''
i = random.randint(1, len(domain)-2)
return r1[0:i]+r2[i:]
# 构造初始种群
pop = []
for i in range(popsize):
vec = [random.randint(domain[i][0], domain[i][1])
for i in range(len(domain))]
pop.append(vec)
# 每一代中有多少胜出者
top_elite = int(elite*popsize)
# 主循环
for i in range(maxiter):
scores = sorted([(costf(v), v) for v in pop])
ranked = [v for (s, v) in scores]
# 从纯粹的胜出者开始
pop = ranked[0:top_elite]
# 添加变异和配对基因重组后的胜出者
while len(pop) < popsize:
if random.random() < mut_prob:
# 变异
c = random.randint(0, top_elite)
sol_mut = mutate(ranked[c])
if sol_mut: # mutate的第三种情况else未指定返回的是None，这时候会报错
pop.append(sol_mut)
else: # 基因重组
c1 = random.randint(0, top_elite)
c2 = random.randint(0, top_elite)
pop.append(crossover(ranked[c1], ranked[c2]))
# 打印当前最优值
#print(scores[0][0])
return scores[0][1]
# 运行程序10次，查看每次得到的最优解
for i in range(10):
gs = genetic_optimize(domain, dormcost)
print(dormcost(gs))
# 输出
3
4
2
5
5
4
6
7
4
5
可以看到，遗传算法得到的解已经非常优秀了，最小成本仅为2，也就是说，仅有2个同学选到的是次选项，其余同学都分配到了最想要的结果。
将几种算法结果放在一起进行比较
最后，让我们将这几种方法均运行100次，然后进行可视化，观察一下这几种算法的解如何。
import matplotlib.pyplot as plt
plt.style.use('ggplot')
rand_s = []
hill_s = []
ann_s = []
gen_s = []
for i in range(100):
rand_s.append(dormcost(random_optimize(domain, dormcost)))
hill_s.append(dormcost(hill_climb(domain, dormcost)))
ann_s.append(dormcost(annealing_optimize(domain, dormcost)))
gen_s.append(dormcost(genetic_optimize(domain, dormcost)))
print(sum(rand_s)/100)
print(sum(hill_s)/100)
print(sum(ann_s)/100)
print(sum(gen_s)/100)
# 输出
19.87
14.65
12.09
4.39
# 可视化
df = pd.DataFrame({'rand':rand_s, 'hill':hill_s, 'ann':ann_s, 'gen':gen_s})
df.plot()
从上面可以看出，遗传算法要明显的由于其他的几种算法。
关于“传统”优化算法与梯度下降优化算法的一点思考
首先需要指出的是，将上面的爬山算法，退火算法，遗传算法称为传统优化算法仅是为了方便与梯度优化算法进行比较，在其他地方并未见提及。
让我们继续以上面这个学生宿舍分配的例子来进行说明。在梯度下降中可以将每个参数视为是一个维度，因此梯度下降的目的就是要在一个多维空间中找到找到一个最优解，若图示化的话则是找到一个谷底。我们将这种思想放入传统优化算法中，对于上面的宿舍分配问题，每个学生就可以看作一个维度，因此同样的我们可以认为传统算法也是在一个多维空间中找到最优解，只是其算法不及梯度下降稳定。因为传统优化算法中含有更多的随机性，即便是遗传算法，我们也无法确定每一次突变或基因重组生成的是更优的还是一个更差的解，其更多的是通过加多迭代次数和种群筛选来保证一个结果不会变的更差。而梯度下降则是一个非常稳定的算法，其保证了每次迭代都是向着损失（或者成本）最小的方向走的，撇开计算量不谈，从这个方面也可以看出来为什么现在的机器学习、深度学习算法大都是用的梯度下降算法了。最后，可能有人会提到，梯度下降虽然稳定，但是其无法保证找到全局最优解，从这一点来说，传统算法和梯度下降算法是相似的，两者都是无法保证能找到全局最优解的。而从结果来说，大部分时候一个局部最优解其实已经基本满足我们的需要了。
最后，这只是一个非数学计算机专业算法爱好者的一些思考，不能保证正确性，如有疑问欢迎指正。
参考：
《集体智慧编程》
