从几个简单例子谈随机优化技术
阅读目录(Content)
1. 关于随机优化（stochastic optimization）
2. 日常和工程中常见的需要用到优化技术的典型场景
0x1：无约束搜索空间问题的随机优化 - 组团旅游问题
1. 描述题解
2. 成本函数（cost function）
0x2：带约束搜索空间问题的随机优化 - 学生宿舍安排问题
1. 描述题解
2. 成本函数
0x3：网络可视化布局问题
1. 布局问题（题解描述）
2. 计算交叉线（损失函数）
3. 随机搜索
4. 爬山法
5. 模拟退火算法（simulated annealing algorithm）
6. 遗传算法（genetic algorithm）
7. 极大似然估计
8. SGD（随机梯度下降法）
回到顶部(go to top)
1. 关于随机优化（stochastic optimization）
随机优化技术常被用来处理协作类问题，它特别擅长处理：受多种变量的影响，存在许多可能解的问题，以及结果因这些变量的组合而产生很大变化的问题。例如：
在物理学中，研究分子的运动
在生物学中，预测蛋白质的结构
在计算机科学中，预测算法的最坏可能运行时间
NASA甚至使用优化技术来设计具有正确操作特性的天线，而这些天线看起似乎不像是人类设计者创造出来的
优化算法是通过尝试许多不同解并给这些解打分以确定其质量的方式，来找到一个问题的最优解。优化算法的典型应用场景是，存在大量可能的解（解搜索空间）以至于我们无法对其进行一一尝试的情况。
最粗暴简单但也是最低效的求解方法，是去尝试随机猜测上千个题解，并从中找出最佳解。而相对的，更有效率的方法，则是一种对题解可能有改进的方式来对其进行智能地修正。
优化算法是否管用的最关键因素：
文章的前面提到了优化算法的几个核心因素，但这里笔者希望强调的是，一种优化算法是否管用很大程度取决于问题本身。大多数优化算法都依赖于这样一个事实：最优解应该接近于其他的优解（即损失函数连续可微），来看一个可能不起作用的例子，
在图的右边，成本的最低点实际上处在一个非常陡峭的峡谷区域内，接近它的任何解都有可能被排除在外，因为这些解的成本都很高，所以我们可能永远都无法找到通往全局最小值的途径。大多数优化算法都会陷入到图中左边某个局部最小化的区域里。
Relevant Link:
《集体智慧编程》Toby segaran - 第5章
回到顶部(go to top)
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章
回到顶部(go to top)
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,