從幾個簡單例子談隨機優化技術

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每次實際選擇的前進迭代方式又具有了必定的隨機性。