从几个简单例子谈随机优化技术
1. 关于随机优化（stochastic optimization）
随机优化技术常被用来处理协作类问题，它特别擅长处理：受多种变量的影响，存在许多可能解的问题，以及结果因这些变量的组合而产生很大变化的问题。例如：
- 在物理学中，研究分子的运动
- 在生物学中，预测蛋白质的结构
- 在计算机科学中，预测算法的最坏可能运行时间
- NASA甚至使用优化技术来设计具有正确操作特性的天线，而这些天线看起似乎不像是人类设计者创造出来的
优化算法是通过尝试许多不同解并给这些解打分以确定其质量的方式，来找到一个问题的最优解。优化算法的典型应用场景是，存在大量可能的解（解搜索空间）以至于我们无法对其进行一一尝试的情况。
最粗暴简单但也是最低效的求解方法，是去尝试随机猜测上千个题解，并从中找出最佳解。而相对的，更有效率的方法，则是一种对题解可能有改进的方式来对其进行智能地修正。
优化算法是否管用的最关键因素：
文章的前面提到了优化算法的几个核心因素，但这里笔者希望强调的是，一种优化算法是否管用很大程度取决于问题本身。大多数优化算法都依赖于这样一个事实：最优解应该接近于其他的优解（即损失函数连续可微），来看一个可能不起作用的例子，
在图的右边，成本的最低点实际上处在一个非常陡峭的峡谷区域内，接近它的任何解都有可能被排除在外，因为这些解的成本都很高，所以我们可能永远都无法找到通往全局最小值的途径。大多数优化算法都会陷入到图中左边某个局部最小化的区域里。
Relevant Link:
《集体智慧编程》Toby segaran - 第5章
2. 日常和工程中常见的需要用到优化技术的典型场景
这一章节，笔者这里会先描述定义出2种典型场景，分别是无约束搜索问题和带约束搜索问题，它们的区别如下，
- 无约束搜索问题：参数的搜索空间充满了所有随机变量的整个定义域，每一个随机变量的取值都对其他随机变量的取值没有任何影响
- 带约束搜索问题：参数的搜索空间是所有随机变量定义域的一个子集，某个随机变量的取值确定后会对其他余下随机变量的取值产生约束影响
我们日常生活中和工程中遇到的很多问题，都可以抽象为上述两类问题。本文我们会通过3个具体的例子，来从不同方面考察优化算法的灵活性。
0x1：无约束搜索空间问题的随机优化 - 组团旅游问题
第一个例子是关于为一个团队安排去某地开会的往返航班，这种场景在现实生活中很常见，例如笔者所在的公司在全球多地有不同的分部，而每年的BlackHat会议会在一个估计的地点（例如拉斯维加斯）和时间召开，这个时候，就需要为各个不同地方的员工安排航班，使他们在同一天到达，并且同样安排其在同一天返航会各自的城市。
计划的制定要求有许多不同的输入，比如：
- 每个人的航班时间表应该是什么？
- 需要租用多少辆汽车？
- 哪个飞机场是最通畅的？
许多输出结果也必须考虑，比如：
- 总的成本
- 候机时间
- 起飞的时间
很显然，我们对输出的期望是总成本和总候机时间越短越好，但是我们无法将这些输入用一个简单的公式映射到输出。要解决这个问题，就需要转换思维，将公式思维转换到计算机的优化思维上。
1. 描述题解
为来自不同地方去往同一个地点的人们（本例是Glass一家）安排一次旅游是一件极富挑战的事情。下面虚构出一个家庭成员及其所在地，
people = [('Seymour', 'BOS'),
('Franny', 'DAL'),
('Zooey', 'CAK'),
('Walt', 'MIA'),
('Buddy', 'ORD'),
('Les', 'OMA')]
# Laguardia airport in NewYork
destination = 'LGA'家庭成员们来自全美各地，并且它们希望在纽约会面。他们将在同一天到达，并在同一天离开，而且他们想搭乘相同的交通工具往返（到达接机，离开送机）纽约机场，从各自所在地前往所在地的机场的交通不用安排。每天有许多航班从任何一位家庭成员的所在地飞往纽约，飞机的起飞时间都是不同的，这些航班在价格和飞行时间上也都不尽相同。数据示例如下，
LGA,OMA,6:19,8:13,239
OMA,LGA,6:11,8:31,249
LGA,OMA,8:04,10:59,136
OMA,LGA,7:39,10:24,219
LGA,OMA,9:31,11:43,210
OMA,LGA,9:15,12:03,99
LGA,OMA,11:07,13:24,171
OMA,LGA,11:08,13:07,175
LGA,OMA,12:31,14:02,234
OMA,LGA,12:18,14:56,172
LGA,OMA,14:05,15:47,226
数据中每一列代表分别为：起点、终点、起飞时间、到达时间、价格。
接下来的问题是，家庭成员的每个成员应该乘坐哪个航班呢？当然，使总的票价降下来是一个目标，但是还有其他因素要考虑，例如：总的候机时间、等地其他成员到达的时间、总的飞行时间等。
描述题解的第一步要做的是，抽象化表达。
当处理类似这样的问题时，我们有必要明确潜在的题解将如何抽象的表达，使之不局限于某个具体的业务问题场景。有一种非常通用的表达方式，就是数字序列，在本例中，一个数字可以代表某人选择乘坐的航班，0代表第一次航班、1代表第二次，以此类推。因为每个人都要往返两个航班，所以列表长度是人数的两倍。
例如，
[1, 4, 3, 2, 7, 3, 6, 3, 2, 4, 5, 3]
Seymour BOS 8:04-10:11 $ 95 12:08-14:05 $142
Franny DAL 12:19-15:25 $342 10:51-14:16 $256
Zooey CAK 10:53-13:36 $189 9:58-12:56 $249
Walt MIA 9:15-12:29 $225 16:50-19:26 $304
Buddy ORD 16:43-19:00 $246 10:33-13:11 $132
Les OMA 11:08-13:07 $175 15:07-17:21 $129
上述列表代表了一种题解：Seymour搭乘当天的第2次航班从Boston飞往New York，并搭乘当天的第5次航班返回到Boston。Franny搭乘第4次航班从Dallas飞往New York，并搭乘第3次航班返回。
即使不考虑价格，上述安排仍然是有问题的。例如Les的航班是下午3点的，但是因为Zooey的航班是早上9点的，所以所有人都必须在早晨6点到达机场。
基本上看，所有人的航班”到达纽约时间“和”从纽约出发时间“应该尽可能接近，这样总体的候机时间会减少，但是这没有将总票价和总飞行时间考虑进去。所以，为了确定最佳组合，程序需要一种方法来为不同日程安排的各种属性进行量化地评估，从而决定哪一个方案是最好的。
2. 成本函数（cost function）
成本函数是用优化算法解决问题的关键，任何优化算法的目标，就是要寻找一组能够使成本函数的返回结果达到最小化的输入（本例中即为所有人的航班安排序列），因此成本函数需要返回一个值，用以表示方案的好坏。对于好坏程度的度量没有固定不变的衡量尺度，但是有以下几点基本准则：
- 损失函数返回的值越大，表示该方案越差
- 损失函数需要连续可微，这样保证较低损失的方案是接近于其他低损失方案，让优化的过程是连续渐进的
- 坏方案（无效解）的损失值不宜过大，因为这会导致优化算法很难找到一个次优的解（better solution），因为算法无法确定返回结果是否接近于其他优解，甚或是有效的解
- 尽可能让最优解的损失为零，这样当算法找到一个最优解时，就可以让优化算法停止搜寻更优的解
通常来说，对多个随机变量进行综合评估方案的好坏是比较困难的，主要问题在于不同随机变量之间的量纲维度不一致，我们来考察一些在组团旅游的例子中能够被度量的变量，
- 价格：所有航班的总票价，或者有可能是考虑财务因素之后的加权平均
- 旅行时间：每个人在飞机上花费的总时间
- 等待时间：在机场等待其他成员到达的时间（在达纽约机场等其他人一起拼车前往市区）
- 出发时间：早晨太早起飞的航班也许会产生额外的成本，因为这要求旅行者减少睡眠时间
- 汽车租用时间：如果集体租用一辆汽车，那么他们必须在一天内早于起租时刻之前将车归还，否则将多付一天的租金
确定了对成本产生影响的所有变量因素之后，我们就必须要找到办法将他们”组合“在一起行程一个值。
例如在本例中，我们就需要做如下转换：
- 飞机上时间的时间价值：假定旅途中每一分钟值1美元，这相当于，再加90美元选择直达航班，就可以节省1个半小时的时间
- 机场等待时间的时间价值：机场等待中每节省1分钟则价值0.5美元
- 如果汽车是在租用时间之后归还的，还会追加50美元的罚款
def schedulecost(sol):
totalprice = 0
latestarrival = 0
earliestdep = 24 * 60
for d in range(len(sol) / 2):
# Get the inbound and outbound flights
origin = people[d][1]
outbound = flights[(origin, destination)][int(sol[d])]
returnf = flights[(destination, origin)][int(sol[d + 1])]
# Total price is the price of all outbound and return flights
totalprice += outbound[2]
totalprice += returnf[2]
# Track the latest arrival and earliest departure
if latestarrival < getminutes(outbound[1]): latestarrival = getminutes(outbound[1])
if earliestdep > getminutes(returnf[0]): earliestdep = getminutes(returnf[0])
# Every person must wait at the airport until the latest person arrives.
# They also must arrive at the same time and wait for their flights.
totalwait = 0
for d in range(len(sol) / 2):
origin = people[d][1]
outbound = flights[(origin, destination)][int(sol[d])]
returnf = flights[(destination, origin)][int(sol[d + 1])]
totalwait += latestarrival - getminutes(outbound[1])
totalwait += getminutes(returnf[0]) - earliestdep
# Does this solution require an extra day of car rental? That'll be $50!
if latestarrival > earliestdep: totalprice += 50
return totalprice + totalwait
对前面假设的一个航班序列打印其总损失结果为：
Seymour BOS 8:04-10:11 $ 95 12:08-14:05 $142
Franny DAL 12:19-15:25 $342 10:51-14:16 $256
Zooey CAK 10:53-13:36 $189 9:58-12:56 $249
Walt MIA 9:15-12:29 $225 16:50-19:26 $304
Buddy ORD 16:43-19:00 $246 10:33-13:11 $132
Les OMA 11:08-13:07 $175 15:07-17:21 $129
5285
成本函数建立后，我们的目标就非常明确了，就是要通过选择正确地数字序列来最小化该成本。
在文章的后半部分，我们会集中介绍代表当前主流核心思想的集中优化算法，并通过比较不同优化算法之间的区别，来更加深刻理解随机优化的内涵。现在，我们暂时先继续我们对典型问题场景的讨论上。
0x2：带约束搜索空间问题的随机优化 - 学生宿舍安排问题
本节我们来考查另一个不同的问题，这个问题明显要借助优化算法来加以解决。其一般的表述是：如何将有限的资源分配给多个表达了偏好的人，并尽可能使他们都满意，或尽可能使所有人都满意。
1. 描述题解
本节的示例问题是，依据学生的首选和次选，为其分配宿舍。有5间宿舍，每间宿舍只有2个隔间（只能住2人），由10名学生来竞争住所。每一名学生都有一个首选和一个次选，如下：
# The dorms, each of which has two available spaces
dorms=['Zeus','Athena','Hercules','Bacchus','Pluto']
# People, along with their first and second choices
prefs=[('Toby', ('Bacchus', 'Hercules')),
('Steve', ('Zeus', 'Pluto')),
('Karen', ('Athena', 'Zeus')),
('Sarah', ('Zeus', 'Pluto')),
('Dave', ('Athena', 'Bacchus')),
('Jeff', ('Hercules', 'Pluto')),
('Fred', ('Pluto', 'Athena')),
('Suzie', ('Bacchus', 'Hercules')),
('Laura', ('Bacchus', 'Hercules')),
('James', ('Hercules', 'Athena'))]
通过观察数据我们发现，不可能满足所有人的首选。实际上这也是工程中最普遍的情况，全局最优解很少能直接得到，我们大部分时候是在求全局次优解。
现在来思考如何对这个问题进行抽象建模的问题，理论上，我们可以构造一个数字序列，让每个学生对应于一个学生，表示我们将某个学生安排在了某个宿舍。但是问题在于，这种抽象表达方式无法在题解中体现每间宿舍仅限两名学生居住的约束条件。一个全零序列代表了所有人都安置在了Zeus宿舍，但这不是一个有效的解。
解决这个问题的正确办法是，通过在抽象建模时增加一个约束条件，即找到一种让每个解都有效的题解表示法。这等价于缩小了题解搜索空间，强制优化算法在一个受约束的子空间进行最优化搜索。
将每间宿舍抽象出拥有两个”槽“，如此，在本例中共有10个槽，我们将每名学生依序安置在各个空槽内，第一位学生可置于10个槽中的任何一个内，第二位学生则可置于剩余9个槽中的任何一个位置，依次类推。
# [(0, 9), (0, 8), (0, 7), (0, 6), (0, 5), (0, 4), (0, 3), (0, 2), (0, 1), (0, 0)]
domain=[(0,(len(dorms)*2)-i-1) for i in range(0,len(dorms)*2)]
下面代码示范了槽的分配过程，5间宿舍总共有10个槽，每个槽被分配后会从队列中删除，如此，其他学生就不会再被安置于该槽中了。
def printsolution(vec):
slots=[]
# Create two slots for each dorm
for i in range(len(dorms)): slots+=[i,i]
print "slots: ", slots
# Loop over each students assignment
for i in range(len(vec)):
x=int(vec[i])
# Choose the slot from the remaining ones
dorm=dorms[slots[x]]
# Show the student and assigned dorm
print prefs[i][0],dorm
# Remove this slot
del slots[x]
这里以一种安排序列为例，即给每一个学生都安排当前所剩槽位的第一个，
printsolution([0,0,0, 0,0,0, 0,0,0, 0])
Toby Zeus
Steve Zeus
Karen Athena
Sarah Athena
Dave Hercules
Jeff Hercules
Fred Bacchus
Suzie Bacchus
Laura Pluto
James Pluto
显然这是一个有效解，但不是一个最优解。可以看到，槽的分配程序是一个中立的程序，如何分配槽位，完全取决于输入的学生安排序列，优化算法的任务就是找到一个最佳的安排序列，使槽分配后的结果尽可能满足学生的需求。
笔者提醒：
尽管这是一个非常具体的例子，但是将这种情况推广到其他问题是非常容易的，例如纸牌游戏中玩家的牌桌分配、家庭成员中的家务分配。要理解领会的是，这类问题的目的是为了从个体中提取信息，并将其组合起来产生出优化的结果。
2. 成本函数
成本函数的计算过程，是通过将学生的当前宿舍安置情况，与他的两项期望选择进行对比而得到。
- 如果学生当前被安置的宿舍即是首先宿舍，则总成本为0
- 如果是次选宿舍，则成本加1
- 如果不在其选择之内，则加3
def dormcost(vec):
cost=0
# Create list a of slots
slots=[0,0,1,1,2,2,3,3,4,4]
# Loop over each student
for i in range(len(vec)):
x=int(vec[i])
dorm=dorms[slots[x]]
pref=prefs[i][1]
# First choice costs 0, second choice costs 1
if pref[0]==dorm: cost+=0
elif pref[1]==dorm: cost+=1
else: cost+=3
# Not on the list costs 3
# Remove selected slot
del slots[x]
return cost
以前面举例的[0,0,0, 0,0,0, 0,0,0, 0]为例，
s = [0,0,0, 0,0,0, 0,0,0, 0]
printsolution(s)
print dormcost(s)
18
0x3：网络可视化布局问题
这个例子讨论的是网络的可视化问题，这里的网络指的是任何彼此相连的一组事物，像Myspace、Facebook、linkedIn这样的社交应用便是社交网络的一个典型例子，在这些应用中，人们因为朋友或具备特定关系而彼此相连。网站的每一位成员可以选择与他们相连的其他成员，共同构筑一个人际关系网络。将这样的网络可视化输出，以明确人们彼此间的社会关系结构，是一件很有意义和研究价值的事情。
1. 布局问题（题解描述）
为了展示一大群人及其彼此间的关联，我们将网络在一个二维画布上进行绘制，绘制时会遇到一个问题，我们应该如何安置图中的每个人名的空间布局呢？以下图为例，
混乱的随机布局
在图中，我们可以看到Augustus是Willy、Violet和Miranda的朋友。但是，网络的布局有点杂乱，连接线彼此交错的情况比较多，一个更为清晰的布局如下图所示，
本例的目标就是讨论如何运用优化算法来构建更清晰结构的网络图。首先，我们虚构一些数据，这些数据代表着社会网络的某一小部分，
people=['Charlie','Augustus','Veruca','Violet','Mike','Joe','Willy','Miranda']
links=[('Augustus', 'Willy'),
('Mike', 'Joe'),
('Miranda', 'Mike'),
('Violet', 'Augustus'),
('Miranda', 'Willy'),
('Charlie', 'Mike'),
('Veruca', 'Joe'),
('Miranda', 'Augustus'),
('Willy', 'Augustus'),
('Joe', 'Charlie'),
('Veruca', 'Augustus'),
('Miranda', 'Joe')]
此处，我们的目标是建立一个程序，令其能够读取一组有关谁是谁的朋友的事实数据，并生成一个易于理解的网络图。
要完成这项工作，最基本的需要借助质点弹簧算法（mass-and-spring algorithm），这一算法是从物理学中建模而来的，各节点彼此向对方施以推力并试图分离，而节点间的连接则试图将关联节点彼此拉近。如此一来，网络便会逐渐呈现出这样一个布局，未关联的节点被推离，而关联的节点则被彼此拉近，却又不会靠的很近。
但遗憾的是，原始的质点弹簧算法无法避免交叉线，这使得我们追踪一个大型的社会网络变得很困难。解决问题的办法是使用优化算法来构建布局，将比较交叉作为成本构建一个成本函数，并尝试令其返回值尽可能小。
2. 计算交叉线（损失函数）
将每个节点都抽象为一个（x，y）坐标，因此可以将所有节点的坐标放入一个列表中。随后，成本函数只需对彼此交叉的连线进行计数即可，
def crosscount(v):
# Convert the number list into a dictionary of person:(x,y)
loc=dict([(people[i],(v[i*2],v[i*2+1])) for i in range(0,len(people))])
total=0
# Loop through every pair of links
for i in range(len(links)):
for j in range(i+1,len(links)):
# Get the locations
(x1,y1),(x2,y2)=loc[links[i][0]],loc[links[i][1]]
(x3,y3),(x4,y4)=loc[links[j][0]],loc[links[j][1]]
den=(y4-y3)*(x2-x1)-(x4-x3)*(y2-y1)
# den==0 if the lines are parallel
if den==0: continue
# Otherwise ua and ub are the fraction of the
# line where they cross
ua=((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3))/den
ub=((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3))/den
# If the fraction is between 0 and 1 for both lines
# then they cross each other
if ua>0 and ua<1 and ub>0 and ub<1:
total+=1
for i in range(len(people)):
for j in range(i+1,len(people)):
# Get the locations of the two nodes
(x1,y1),(x2,y2)=loc[people[i]],loc[people[j]]
# Find the distance between them
dist=math.sqrt(math.pow(x1-x2,2)+math.pow(y1-y2,2))
# Penalize any nodes closer than 50 pixels
if dist<50:
total+=(1.0-(dist/50.0))
return total
笔者思考：
优化算法就像一个忠实满足我们愿望的”精灵“，你通过损失函数设定什么样的目标，优化的结果就会朝什么方向发展。从这点来说，在利用机器学习解决生活和工程问题的时候，清楚我们想要什么是非常重要的。时常会有题解满足原本的”最优“条件，但却并非是我们想要的结果，这个时候就要反思一下是不是损失函数的目标设置有问题。
Relevant Link:
《集体智慧编程》Toby segaran - 第5章
3. 随机搜索
前面的章节，我们已经列举了几种常见的典型问题场景，笔者这么排版的目的是想明确传达出一个意思：
优化算法是通用独立的，它不依赖于具体问题，而是一种通用的抽象化接口，只要目标问题能够抽象出一个明确的损失函数，优化算法就可以在损失函数的搜索空间中寻找到一个相对最优解。
这节我们从随机搜索开始讨论，随机搜索不是一种非常好的优化算法，但是它却使我们很容易领会所有算法的真正意图，并且它也是我们评估其他算法优劣的基线（baseline）。
def randomoptimize(domain, costf):
best = 999999999
bestr = None
for i in range(0, 1000):
# Create a random solution
r = [float(random.randint(domain[i][0], domain[i][1]))
for i in range(len(domain))]
# Get the cost
cost = costf(r)
# Compare it to the best one so far
if cost < best:
best = cost
bestr = r
return r
这个函数接受两个参数，
- domain：是一个由二元组（2-tuples）构成的列表，它指定了每个变量的最小和最大值。题解的长度与此列表的长度相同。例如旅游团安排的例子中，每个人都有10个往返航班和，因此传入的domain即是由(0，9)构成的list。
- costf：是成本函数，用于计算成本值
此函数会在domain的定义域内随机产生1000次猜测，并对每一次猜测调用costf。并统计最终的最优结果。
import time
import random
import math
people = [('Seymour', 'BOS'),
('Franny', 'DAL'),
('Zooey', 'CAK'),
('Walt', 'MIA'),
('Buddy', 'ORD'),
('Les', 'OMA')]
# Laguardia
destination = 'LGA'
flights = {}
#
for line in file('schedule.txt'):
origin, dest, depart, arrive, price = line.strip().split(',')
flights.setdefault((origin, dest), [])
# Add details to the list of possible flights
flights[(origin, dest)].append((depart, arrive, int(price)))
def getminutes(t):
x = time.strptime(t, '%H:%M')
return x[3] * 60 + x[4]
def printschedule(r):
for d in range(len(r) / 2):
name = people[d][0]
origin = people[d][1]
out = flights[(origin, destination)][int(r[d])]
ret = flights[(destination, origin)][int(r[d + 1])]
print '%10s%10s %5s-%5s $%3s %5s-%5s $%3s' % (name, origin,
out[0], out[1], out[2],
ret[0], ret[1], ret[2])
def schedulecost(sol):
totalprice = 0
latestarrival = 0
earliestdep = 24 * 60
for d in range(len(sol) / 2):
# Get the inbound and outbound flights
origin = people[d][1]
outbound = flights[(origin, destination)][int(sol[d])]
returnf = flights[(destination, origin)][int(sol[d + 1])]
# Total price is the price of all outbound and return flights
totalprice += outbound[2]
totalprice += returnf[2]
# Track the latest arrival and earliest departure
if latestarrival < getminutes(outbound[1]): latestarrival = getminutes(outbound[1])
if earliestdep > getminutes(returnf[0]): earliestdep = getminutes(returnf[0])
# Every person must wait at the airport until the latest person arrives.
# They also must arrive at the same time and wait for their flights.
totalwait = 0
for d in range(len(sol) / 2):
origin = people[d][1]
outbound = flights[(origin, destination)][int(sol[d])]
returnf = flights[(destination, origin)][int(sol[d + 1])]
totalwait += latestarrival - getminutes(outbound[1])
totalwait += getminutes(returnf[0]) - earliestdep
# Does this solution require an extra day of car rental? That'll be $50!
if latestarrival > earliestdep: totalprice += 50
return totalprice + totalwait
def randomoptimize(domain, costf):
best = 999999999
bestr = None
for i in range(0, 100000):
# Create a random solution
r = [float(random.randint(domain[i][0], domain[i][1]))
for i in range(len(domain))]
# Get the cost
cost = costf(r)
# Compare it to the best one so far
if cost < best:
best = cost
bestr = r
return r
def hillclimb(domain, costf):
# Create a random solution
sol = [random.randint(domain[i][0], domain[i][1])
for i in range(len(domain))]
# Main loop
while 1:
# Create list of neighboring solutions
neighbors = []
for j in range(len(domain)):
# One away in each direction
if sol[j] > domain[j][0]:
neighbors.append(sol[0:j] + [sol[j] + 1] + sol[j + 1:])
if sol[j] < domain[j][1]:
neighbors.append(sol[0:j] + [sol[j] - 1] + sol[j + 1:])
# See what the best solution amongst the neighbors is
current = costf(sol)
best = current
for j in range(len(neighbors)):
cost = costf(neighbors[j])
if cost < best:
best = cost
sol = neighbors[j]
# If there's no improvement, then we've reached the top
if best == current:
break
return sol
def annealingoptimize(domain, costf, T=10000.0, cool=0.95, step=1):
# Initialize the values randomly
vec = [float(random.randint(domain[i][0], domain[i][1]))
for i in range(len(domain))]
while T > 0.1:
# Choose one of the indices
i = random.randint(0, len(domain) - 1)
# Choose a direction to change it
dir = random.randint(-step, step)
# Create a new list with one of the values changed
vecb = vec[:]
vecb[i] += dir
if vecb[i] < domain[i][0]:
vecb[i] = domain[i][0]
elif vecb[i] > domain[i][1]:
vecb[i] = domain[i][1]
# Calculate the current cost and the new cost
ea = costf(vec)
eb = costf(vecb)
p = pow(math.e, (-eb - ea) / T)
# Is it better, or does it make the probability
# cutoff?
if (eb < ea or random.random() < p):
vec = vecb
# Decrease the temperature
T = T * cool
return vec
def geneticoptimize(domain, costf, popsize=50, step=1,
mutprob=0.2, elite=0.2, maxiter=100):
# Mutation Operation
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:]
# Crossover Operation
def crossover(r1, r2):
i = random.randint(1, len(domain) - 2)
return r1[0:i] + r2[i:]
# Build the initial population
pop = []
for i in range(popsize):
vec = [random.randint(domain[i][0], domain[i][1])
for i in range(len(domain))]
pop.append(vec)
# How many winners from each generation?
topelite = int(elite * popsize)
# Main loop
for i in range(maxiter):
scores = [(costf(v), v) for v in pop]
scores.sort()
ranked = [v for (s, v) in scores]
# Start with the pure winners
pop = ranked[0:topelite]
# Add mutated and bred forms of the winners
while len(pop) < popsize:
if random.random() < mutprob:
# Mutation
c = random.randint(0, topelite)
pop.append(mutate(ranked[c]))
else:
# Crossover
c1 = random.randint(0, topelite)
c2 = random.randint(0, topelite)
pop.append(crossover(ranked[c1], ranked[c2]))
# Print current best score
print scores[0][0]
return scores[0][1]
if __name__ == '__main__':
#s = [1, 4, 3, 2, 7, 3, 6, 3, 2, 4, 5, 3]
#printschedule(s)
#print schedulecost(s)
domain = [(0, 9)] * (len(people) * 2)
s = randomoptimize(domain, schedulecost)
print schedulecost(s)
printschedule(s)
# s = hillclimb(domain, schedulecost)
# print schedulecost(s)
# printschedule(s)
# s = annealingoptimize(domain, schedulecost)
# print schedulecost(s)
# printschedule(s)
#s = geneticoptimize(domain, schedulecost)
#printschedule(s)
随机的结果并不算太好，当然也不算太差，而且最大的问题在于，优化结果不可预期，每次运行都可能得到不同的结果。我们希望优化算法可以稳定地获得成功。
Relevant Link:
《集体智慧编程》Toby segaran - 第5章
4. 爬山法
随机尝试的问题在于，没有充分利用已经发现的优解。具体来说，拥有较低总成本的时间安排很可能接近于其他低成本的安排。但是因为随机优化是到处跳跃的（jump around），所以它不会自动去寻找与已经被发现的优解相接近的题解。
本章要介绍的方法叫爬山法（hill climbing），爬山法以一个随机解开始，然后在其临近的解集中寻找更好的题解，这类似于从斜坡向上下走，如下图所示，
这种爬山法的合理性在于，充分利用了损失函数本身的连续可微特性，每次的尝试都建立在充分利用历史尝试的信息之上，步步为营，逐渐找到一个局部最优解。
我们可以应用这种爬山法为Glass一家找到最好的旅行安排方案，
s = hillclimb(domain, schedulecost)
print schedulecost(s)
printschedule(s)
爬山法找到的解通常比随机搜索要好，但是，爬山法有一个较大的缺陷，如下图，
简单从斜坡滑下不一定会产生全局最优解。爬山法得到的只能是一个局部范围内的最小值，它比邻近解的表现逗号，但却不是全局最优的。
一个缓解的手段被称为随机重复爬山法（random-restart hill climbing），即让爬山法以多个随机生成的初始解为起点运行若干次，并取其中最优解，借此希望在多次的随机尝试中，有一个解能逼近全局的最小值。
当然，爬山法现在已经没有使用了，但是我们却能从中领悟到优化算法进化的几个核心要素，
- 渐进优化思想：在定义域内进行局部的小范围修改，充分利用历史优化结果，可以至少保证得到局部最优
- 随机扰动思想：随机扰动可以扩大优化算法的搜索范围，避免落入局部最优的鞍点
- 稳定、可预期、可重入的获得最优解（或次优解）是优化算法的一个良好特性
Relevant Link:
《集体智慧编程》Toby segaran - 第5章
5. 模拟退火算法（simulated annealing algorithm）
模拟退火算法是受物理学领域启发而来的一种优化算法。退火是指将合金加热后再慢慢冷却的过程。
大量的原子因为受到激发而向周围跳跃，然后又逐渐稳定到一个低能阶的状态（例如熔炼武器的过程），所有这些原子能够找到一个低能阶的配置（configuration）。
退火跃迁公式如下，退火算法以一个问题的随机解开始，它用一个变量来表示温度，这一温度初始时很高，之后逐渐变低：
- 每一次迭代中，算法会随机选中题解中的某个数字（例如Seymour的返程航班从第二趟转移到第三趟），这个随机的策略保证了算法具备速记扰动性。
- 新选出的题解会和上一次题解进行损失比较，
- 如果新的成本值更低，则新题解将成为当前题解，这和爬山法的渐进优化特性一致。
- 但是新的成本更高，则按照上面退火跃迁公式计算概率，根据概率结果决定是否要采集该题解为当前题解。可以看到，即使成本值更高，新题解仍然有可能成为当前题解，这是为了避免限于局部最优的一种概率性尝试。可以看到，在开始时，温度非常高，概率几乎为1。
- 不断迭代，直到退火完毕，温度将为0，模拟退火算法完全退化为爬山法，进入最后的局部渐进收敛
在某些情况下，在我们能够得到一个更优的解之前转向一个更差的解是有必要的，这就相当于黎明前的黑暗。模拟退火算法之所以管用，不仅因为它总是接受一个更优的解而且还因为它在退火过程的开始阶段会接受表现较差的解。而随机退火过程的不断进行，算法越来越不可能接受较差的解，收敛逐渐稳定，
def annealingoptimize(domain, costf, T=10000.0, cool=0.95, step=1):
# Initialize the values randomly
vec = [float(random.randint(domain[i][0], domain[i][1]))
for i in range(len(domain))]
while T > 0.1:
# Choose one of the indices
i = random.randint(0, len(domain) - 1)
# Choose a direction to change it
dir = random.randint(-step, step)
# Create a new list with one of the values changed
vecb = vec[:]
vecb[i] += dir
if vecb[i] < domain[i][0]:
vecb[i] = domain[i][0]
elif vecb[i] > domain[i][1]:
vecb[i] = domain[i][1]
# Calculate the current cost and the new cost
ea = costf(vec)
eb = costf(vecb)
p = pow(math.e, (-eb - ea) / T)
# Is it better, or does it make the probability
# cutoff?
if (eb < ea or random.random() < p):
vec = vecb
# Decrease the temperature
T = T * cool
return vec
使用模拟退火算法来对旅行规划问题进行优化，
s = annealingoptimize(domain, schedulecost)
print schedulecost(s)
printschedule(s)
可以看到，模拟退火的优化结果由于随机搜索。
笔者思考：
模拟退火有效的兼顾了”随机扰动“和”渐进收敛“这两个目标的平衡，对于一个优化算法来说，如果只是一味地渐进优化（例如爬山法），则会陷入局部最优的陷阱中。另一方面，如果像随机重复爬山法那样一味地不断地扰动，则又无法保证”可重复收敛“到某一个全局优解上。模拟退火通过一个递减的函数（退火）来动态的调整每一轮优化的扰动幅度，逐步减小扰动幅度，这样做的好处是使得优化算法能稳定的保证经过一定次数的优化后，稳定地收敛到某个全局优解上。
Relevant Link:
《集体智慧编程》Toby segaran - 第5章
6. 遗传算法（genetic algorithm）
遗传算法也是受自然科学的启发，这类算法的运行过程是先随机生成一组解，我们称之为种群（population），在优化过程中的每一步，算法会计算整个种群的成本函数，从而得到种群每一个题解的有序列表，如下列表所示，
在对题解进行排序之后，选出排名前列的一部分（例如10%）题解作为下一代种群的种子，这个过程被称为精英选拔法（elitism）。
接下来，新的种群的种子题解，要经过一次”繁衍“得到一个新的种群。具体的繁衍方式有两种：
- 变异（mutation）：对一个既有的解进行微小的、简单的、随机的改变。例如从题解中选择一个数字，对其进行递增或递减即可。
- 交叉（crossover）或配对（breeding）： 选择最优解中的两个解，将它们按某种方式结。例如从一个解中随机取出一个数字作为新题解中的某个元素，而剩余元素则来自于另一个题解。
这样，一个新的种群（我们称之为下一代）被创建出来。它的大小与旧的种群相同，而后，这个过程会一直重复进行。达到迭代次数，或者连续数代后题解都没有得到改善，则结束优化。
def geneticoptimize(domain, costf, popsize=50, step=1,
mutprob=0.2, elite=0.2, maxiter=100):
# Mutation Operation
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:]
# Crossover Operation
def crossover(r1, r2):
i = random.randint(1, len(domain) - 2)
return r1[0:i] + r2[i:]
# Build the initial population
pop = []
for i in range(popsize):
vec = [random.randint(domain[i][0], domain[i][1])
for i in range(len(domain))]
pop.append(vec)
# How many winners from each generation?
topelite = int(elite * popsize)
# Main loop
for i in range(maxiter):
scores = [(costf(v), v) for v in pop]
scores.sort()
ranked = [v for (s, v) in scores]
# Start with the pure winners
pop = ranked[0:topelite]
# Add mutated and bred forms of the winners
while len(pop) < popsize:
if random.random() < mutprob:
# Mutation
c = random.randint(0, topelite)
pop.append(mutate(ranked[c]))
else:
# Crossover
c1 = random.randint(0, topelite)
c2 = random.randint(0, topelite)
pop.append(crossover(ranked[c1], ranked[c2]))
# Print current best score
print scores[0][0]
return scores[0][1]
运用遗传算法对旅行计划进行优化，
2956
Seymour BOS 9:45-11:50 $172 9:58-11:18 $130
Franny DAL 9:08-12:12 $364 9:49-13:51 $229
Zooey CAK 9:15-12:14 $247 8:19-11:16 $122
Walt MIA 7:34- 9:40 $324 8:23-11:07 $143
Buddy ORD 8:25-10:34 $157 9:11-10:42 $172
Les OMA 9:15-12:03 $ 99 8:04-10:59 $136
可以看到，遗传算法比模拟退火算法的优化效果要好。
笔者提醒：
需要注意的是，遗传算法中，同样体现了优化算法的几个核心要素，
- 随机扰动性：通过变异或配对的方式生成新的种群过程，本质就是引入了一个随机扰动，而随机扰动的加入，一定程度上缓解了陷入局部最优的问题。
- 渐进收敛性：每轮种群的精英拔过程都是选取当前种群中top优的题解，充分利用了上一轮优化的结果，是一种渐进优化思想的体现。
Relevant Link:
《集体智慧编程》Toby segaran - 第5章
7. 极大似然估计
多元非线性损失函数，直接求解其逆矩阵很困难，因此不适合用极大似然估计直接得到全局最优解
8. SGD（随机梯度下降法）
随机梯度下降思想上类似于爬上法，最大的区别在于，
- 爬山法每次都所有的随机变量都进行修改，在各个方向上都相对于原始值偏离相同的step size。
- 而随机梯度会通过计算多元随机变量的梯度，找到当前下降最快的一个梯度方向（例如可能是某几个随机变量偏离相同的step size）。
多元随机函数的梯度计算需要用到Hessian矩阵，这里不再详述。
梯度下降这块遵循了优化算法的渐进收敛特性。而在随机扰动方面，SGD每轮迭代只是随机地选取一部分样本计算损失函数的梯度，随机选取的样本子集可能不一定代表了真实全集的梯度方向，因此，SGD每次实际选择的前进迭代方式又具备了一定的随机性。