從幾個簡單例子談隨機優化技術
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每次實際選擇的前進迭代方式又具備了一定的隨機性。
&n