Optimization

而世之奇偉、瑰怪、非常之觀,常在於險遠,         

而人之所罕至焉,故非有志者不能至也。《王安石.游褒禪山記》

Optimization

求根找到零,最佳化找到極值。找到確切的輸入數值與輸出數值,輸出數值是最大值(最小值),稱作「最佳化」。

以函數圖形表達函數的極值:最大值就是最高的地方,最小值就是最低的地方。有時候最大值、最小值會在無限遠的地方。

中學數學談過一點點多項式函數的最佳化,比如一元二次多項式函數的最佳化(求拋物線的頂點)。大學微積分課程也談過多項式函數最佳化,比如一階導數等於零。

以下則是談一般的函數的最佳化。

特殊函數類型

容易找到極值的函數類型

Unimodal Function:單峰函數。先嚴格遞增、再嚴格遞減的函數。只有遞增或者只有遞減,也算是單峰函數。

Convex Function:凸函數。隨便劃一道割線,函數曲線總是凹下去的函數。凸函數是單峰函數。凸函數的斜率是遞增函數。

Concave Function:凹函數。凸函數上下顛倒。注意到凸函數看起來是凹的,凹函數看起來是凸的,不要問我為什麼。

Lipschitz Function:平緩函數。任意割線的斜率的絕對值,小於等於一個定值。

Polynomial Function:多項式函數,無限可微的函數。「梯度等於零」的地方是極值、鞍點,可以手工推導公式解。

Ternary Search(三分搜尋)

求根可用二分搜尋,求極值可用三分搜尋。

適用於單峰函數、凸函數、凹函數。

UVa 1476 10385 11243 11562 11613 ICPC 5009 7598

Ensemble Average(總體平均)

三個臭皮匠勝過一個諸葛亮。實施多種最佳化演算法,得到多種極值。如果這些極值位置們都不是正確位置,那麼這些極值位置們的加權平均值,通常更接近正確位置。權重由人類手動調整。

適用於凸函數、凹函數。

地面勘查類型

概論

Show[ Plot3D[BesselJ[0, Norm[{x, y}]], {x, -4Pi, 4Pi}, {y, -4Pi, 4Pi}, PlotRange -> {-1, 1}, Boxed -> False, Mesh -> None, Axes -> False, PlotPoints -> 50, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-4, 4}]] &)], ParametricPlot3D[{{-5, u, BesselJ[0, Norm[{-5, u}]]}, {u, -8, BesselJ[0, Norm[{u, -8}]]}}, {u, -4Pi, 4Pi}, PlotStyle -> {Red, Red}] ]
ParametricPlot3D[{v Cos[u], v Sin[u], 5 Cos[Pi/10] BesselJ[0, v]}, {u, 0, 4 (Pi/2)}, {v, 0, 15}, Boxed -> False, Axes -> False, Mesh -> None]

選一個起點,一步一腳印,走向極值。

Hill Climbing(登山法)

沿著函數圖形的表面前進。隨便往某個方向跨出一步,確定是往上,就直直走;確定是往下,就不走。最後成功登頂。

步伐太大,會越過山峰,走往低處。只要步伐持續變小,仍可攻頂。但是步伐變小太快,便會在山腰掙扎,無法登頂。取捨兩難。

選擇各種起點,登上各個山峰。由於無法預測最高峰的位置,只好嘗試多種起點、實施多次登山法,依賴運氣尋找最高峰。嘗試越多起點,越能找到最高峰,越能避免落入區域極值。

適用:連續函數。缺陷:可能停在山腰。可能只找到區域極值。

Simulated Annealing(模擬退火法)

登山法改良版。允許往下走,以便翻山越嶺。

隨便往某個方向跨出一步,確定是往上,就走;確定是往下,以機率exp(Δf⋅t)決定走或不走,其中Δf是上升高度(往下走時是負值),t是時刻(大於等於1)。大致上來說,往下越陡就越不想往下,年紀越大就越不想往下。

註:原論文只找山谷、未找山峰。原論文沒有時刻t,而是溫度T;溫度不斷下降,因此機率公式是exp(Δf/T) 。

適用:連續函數。缺陷:可能停在山腰。可能只找到區域極值。

UVa 10228 ICPC 5102

Gradient Descent(梯度下降法)

登山法改良版。直接朝著梯度方向走,不必試誤。

X座標、Y座標、……分開處理,互不干涉。當前位置,求得相鄰高度差,座標大的高度減去座標小的高度,正負值可判別相對高低。當前位置座標,加上相鄰高度差,就是往上走。若相鄰高度差越多,則前後座標差越多。坡越陡、跨越遠、走越快。

Plot[BesselJ[0, Norm[{-5, y}]], {y, -4Pi, 4Pi}, PlotRange -> {-1,1}] Plot[BesselJ[0, Norm[{x, -8}]], {x, -4Pi, 4Pi}, PlotRange -> {-1,1}]

梯度是相鄰高度差再除以dx,功效差不多。因為數學家喜歡梯度,所以採用梯度。為了得到梯度,函數必須一次可微。

梯度大小是傾斜程度。梯度方向是最陡的方向,是等高線的垂直方向。沿梯度方向走會上升、得極大值,沿梯度反方向走會下降、得極小值。步伐太大,會走之字路線,無傷大雅。

註:原論文只找山谷、未找山峰,因而取名為「下降」。

適用:一次可微函數。優點:方向明確,不必隨機亂走試誤,攻頂速度快。缺陷:可能停在山腰。可能只找到區域極值、鞍點。

Nonlinear Conjugate Gradient Method(非線性共軛梯度法)

梯度下降法改良版。前進方向是「本次的梯度」與「上次的前進方向」的加權平均。權重有許多種選擇,每個人自有一套理論:

http://people.cs.vt.edu/~asandu/Public/Qual2011/Optim/Hager_2006_CG-survey.pdf

http://sebastianruder.com/optimizing-gradient-descent/

最近走紅的演算法有AdaDelta、ADAM、Momentum,層出不窮。

註:演算法名稱雖源自「Conjugate Gradient Method」,內容卻毫無關聯。古代數學家考慮欠周的結果。

Coordinate Descent(座標下降法)

直接朝著座標軸方向走,走到最高處。每個座標軸輪流一遍。

偏微分等於零的地方就是該方向的極值與鞍點。

如果偏微分等於零、可以推導公式解,那麼就能一步到位。

如果偏微分等於零、無法推導公式解,那麼只能步步登高。

適用:一次可微函數。優點:有了公式解,不必煩惱步伐大小。缺陷:故意繞遠路,攻頂時間更久。

Coordinate Gradient Descent(座標梯度下降法)

直接朝著座標軸方向走一步。每個座標軸輪流一遍。

劣於梯度下降法,毫無可取之處。

過程宛如「Gibbs Sampling」,最佳化與取樣的橋樑!

Newton's Method(牛頓法)

梯度等於零的地方是極值與鞍點(不考慮無限遠處)。

牛頓法原先是求根的方法,找到函數等於零的地方;但只要事先求得梯度,就可以套用牛頓法求得極值(連同鞍點)。

求得梯度需要一次微分,牛頓法的過程需要再一次微分,前後合計兩次,等同二次微分。

輸入變數只有一個的情況下,牛頓法非常實用。牛頓法可以視作梯度下降法的改良版,步伐大小是二次微分的倒數(曲率半徑越大、路線越直、步伐越大),攻頂速度更快。

輸入變數有許多個的情況下,牛頓法不太實用。當函數有N個輸入變數,其梯度有N個方向,得到N個函數。N個函數同時實施牛頓法,每一步需時O(N²),計算時間更久。

適用:二次可微函數。缺陷:繼承全部缺陷,計算時間更久。

Quasi-Newton Method(擬牛頓法)(類牛頓法)

牛頓法改良版。二次微分的部分,改成其他類似的東西。有許多種選擇,每個人自有一套理論:

https://en.wikipedia.org/wiki/Quasi-Newton_method

最近走紅的演算法有BHHH、BFGS,層出不窮。

註:古時候尚未發明程式語言。古人遇到迴圈,習慣寫成向量運算、矩陣運算、級數。簡單明瞭的數學概念,經過古人的轉換,往往讓現代人摸不著頭緒。請讀者好自為之。

Fixed Point Iteration(不動點遞推法)

求根、求不動點是同一件事。

不動點遞推法原先是求不動點的方法,找到輸出等於輸入的地方;但只要事先求得梯度、加x,就可以套用不動點遞推法求得極值(連同鞍點)。

  stationary points of f(x)
= roots of ∇f(x)
= roots of λ ∇f(x)
= fixed points of λ ∇f(x) + x

不動點遞推法就是梯度下降法!不動點遞推法是正統根源,擁有清楚的收斂條件:λ ∇f(x) + x是平緩函數。平緩函數保證收斂,不是平緩函數則可能收斂、可能不收斂。

  fixed point iteration of λ ∇f(x) + x
= gradient descent of f(x)

Proximal Point Algorithm(鄰近點演算法)

不動點遞推法改良版。故事本應劃下句點,結果有人故意搞事。首先故意把λ ∇f(x) + x拿去積分,發現了拋物線x²。

∇g(x) = λ ∇f(x) + x
g(x) = λ f(x) + 0.5 x²   (ignore constant)

接著故意創造prox函數:額外創造維度z,令拋物線x²自由平移,找到最小值位置;又除以λ,不影響最小值位置。最後把不動點遞推法的λ ∇f(x) + x,等價換成prox函數,重新稱作鄰近點演算法。應該是為了裝逼。

                             1
prox  (x) = arg min { f(z) + —— ‖z - x‖² } = λ ∇f(x) + x = ∇g(x)
    λf           z           2λ

空中勘查類型

概論

地面勘查經常癱在山腰、卡在山丘,只好改為飛行模式。像偵察機一樣飛來飛去,不被山峰峽谷阻擋。

這類的演算法,完全憑感覺,一個比一個唬爛,連一個數學性質都不必證明。但是我們也沒有更好的方法了。面對亂七八糟的函數,只好用亂七八糟的算法。

Particle Swarm Optimization(粒子演算法)

https://www.youtube.com/watch?v=_bzRHqmpwvo

一、散佈N個粒子。一個粒子對應一個可行解。
  position[1...N];	// x
  solution[1...N];	// f(x)
二、以個人過去最佳解、團體過去最佳解,決定粒子的速度。
  然後移動粒子,讓粒子飛往最佳解。
  velocity[i] =
   2 * rand(0 ~ 1) * (best_position[i] - position[i]) +
   2 * rand(0 ~ 1) * (best_position[best_particle_index] - position[i]);
三、更新個人最佳解、團體最佳解。
  best_solution[i] = max(best_solution[i], solution[i])
    and record best_position[i];
  best_particle_solution = max_arg(best_solution[1...N])
    and record best_particle_index;
四、重複二三,直到解夠好。

附帶一提,這與粒子的真實行為完全無關。

Bee Colony Optimization(蜜蜂演算法)

https://www.youtube.com/watch?v=zxcb6ZBj5PE

一、散佈N個食物。一個食物對應一個可行解。
  position[1...N];	// x
  solution[1...N];	// f(x)
二、N隻蜜蜂分頭前往N個食物並返回,但是腦中記得的位置會有偏差。
  position[i] +=
   rand(-1 ~ +1) * position[rand(1 ~ N except i)];
三、每隻蜜蜂皆從N個食物中挑一個去,機率由解的好壞決定。
  probability[i] = solution[i] / ( solution[1] + ... + solution[N] )
  返回時腦中記得的位置一樣會有偏差。公式同二。
四、如果某個食物連續K個回合沒去,食物自動消滅。
  隨機散佈食物,補足食物直到N個。
五、重複二三四,直到解夠好。

附帶一提,這與蜜蜂的真實行為完全無關。

Fruit Fly Optimization(果蠅演算法)

《果蠅最佳化演算法:最新演化式計算技術》

排列組合類型

概論

隨意拼湊函數輸入,試著讓函數輸出是極值。

這類的演算法,不僅憑感覺,而且還是隨機亂猜的,唬爛無上限。一個演算法,就能開一學期的課,令人不得不嘖嘖稱奇。

Genetic Algorithm(基因演算法)

https://www.youtube.com/watch?v=ejxfTy4lI6I

靈感來自於染色體減數分裂的過程,優良的基因不斷遺傳下去,逐代演化出更適應環境的基因。基因演算法把答案比擬成染色體,把好的答案不斷分裂再結合,成為更好的答案。

附帶一提,這與基因的真實行為完全無關。

1.
[初始化]
一開始先隨便弄出幾個x。本例是四個。

	1010101010
	1011001011
	1110101011
	0010101000

2.
[fitness function]
根據問題特性,定義好壞程度。

	f(1010101010) = 678

3.
[selection]
隨便找個位置切一刀,每個x都被分成兩段。

	1010101  010
	1011001  011
	1110101  011
	0010101  000

4.
[crossover]
隨便找兩組你覺得夠優良的x,交叉相接變成新答案。
重複一直做,直到x數目跟原先一樣多。本例是四個。

	1010101 \/ 010  ->  1010101 -- 011
	1011001 /\ 011      1011001 -- 010 


	1010101011
	1011001010
	1110101010
	1010101000

5.
[mutation]
每個x都隨便找一個地方把數字改掉,也可以不改。

	1010111011
	1011001000
	1110101010
	1010101001

6.
重複3. 4. 5.,直到裡面有一個x是你滿意的,令f(x)最大的那個x。
1. 隨機產生N個x。
2. 計算fitness function。
3. 重複以下步驟,直到有一個x讓人滿意。
 甲、selection。
 乙、crossover。
 丙、mutation。
 丁、計算fitness function。

演算法的過程,大致就是如此。細部的實作方式、參數的調校方式,隨人話虎爛。

一開始的x的足夠豐富,多演化幾次,就可以得到不錯的結果。一開始的x足夠豐富,可以避免進入區域極值。mutation用於增加x的豐富性,以跳脫區域極值。

範例:0/1 Knapsack Problem

N個物品。選出其中幾個。

[初始化]
答案設計成N個二進位數字,
0代表對應的物品不在背包內,
1代表對應的物品在背包內。

[fitness function]
是背包內物品總價值,數值越大越好。

當「特定組合」具有深邃的影響力,造成最佳解、次佳解們都包含著同一(幾)套「特定組合」,此時就適合使用基因演算法。以0/1背包問題為例,一套好的物品組合方式可以有效節省背包空間,只要依賴幾套好的物品組合方式,就留有足夠的背包空間,能夠隨心所欲的放入其他物品;所有接近最佳解的答案,答案的其中一小段,都會是那幾套固定的物品組合方式──這種情況就非常適合使用基因演算法。

UVa 10715

範例:Minimum Spanning Tree

E條邊。選出V-1條。

[初始化]
答案設計成E個二進位數字,
0代表對應的邊不是MST。
1代表對應的邊是MST。

[fitness function]
是MST的權重,數值越大越好。

用人眼觀察,很直覺的就能在小區域找出最佳的子樹,那便是一套「特定組合」。

範例:Travelling Salesman Problem

N個城市,C(N,2)條路。選出N條。

[初始化]
答案被迫設計成0到N-1的數字,代表一條路徑。

[fitness function]
是路徑的權重,數值越小越好。

[crossover]
哈哈哈,很難搞。

[mutation]
可以有好幾種方式:
1. 隨便交換兩個數字。
2. 隨便找一段數字,顛倒前後順序。
3. 隨便拿出一個數字,隨便插入到其他地方。

旅行推銷員問題也具有「特定組合」的性質,只不過它的答案並不適合分裂再結合,最好不要堅持使用基因演算法,自尋煩惱。

範例:Video Game

Ant Colony Optimization(螞蟻演算法)

https://www.youtube.com/watch?v=hXUCCRiNBOc

靈感來自於螞蟻覓食的過程,螞蟻發現食物後會沿途留下強烈的費洛蒙,驅使其他螞蟻沿著路線前進,不斷留下更多費洛蒙,吸引更多螞蟻;也有螞蟻會另闢新路,而找到更簡潔的路線。螞蟻演算法把答案比擬成螞蟻覓食的路徑,把好的答案不斷做局部調整,成為更好的答案。

螞蟻演算法有各式各樣的版本,這裡介紹其中一個經典版本Ant Colony System,主要是計算一條最短的覓食路徑。

附帶一提,這與螞蟻的真實行為完全無關。

1.
[初始化]
一開始先建立一個地圖,地圖上有P個地點。
有一些地點是食物,有一個地點是蟻窩。

地點與地點間皆有道路,
所有道路的費洛蒙都預設為1.0。

2.
[state transition rule]
一隻螞蟻從蟻窩出發。
每次踏上一個地點,螞蟻有兩種選擇,
  ◎探索:隨便走一條路。機率為q。
  ◎追蹤:走費洛蒙最強的道路。機率為1-q。

q是螞蟻選擇探索的機率,自行設定在0到1之間。

為了不讓螞蟻打轉繞圈,所以會讓螞蟻避開去過的地點。
在探索和追蹤前,要先過濾掉去過的地點。

所有食物都找到後就直接回蟻窩,
沒找到所有食物就不准回蟻窩。
總之就是要刻意弄出一條「嘗遍天下美食」的路線。

3.
[local updating rule]
螞蟻回巢之後,
剛剛走過的每一條道路,費洛蒙都會加強,
道路的費洛蒙會變成下列兩者相加後的數值,
  ◎自然蒸發,餘下:原本費洛蒙值,乘上1-ρ。
  ◎螞蟻路過,添加:道路起點所連接的道路當中,最大的那個費洛蒙值,乘上p。

ρ是費洛蒙蒸發的程度,自行設定在0到1之間。
通常會把參數設的很好,讓相加之後,費洛蒙會比原本增加一些,而不是變更少。

4.
有N隻螞蟻,依序做2. 3.這件事。

5.
[global updating rule]
N隻螞蟻回巢之後,
地圖上所有道路的費洛蒙都會蒸發一定比例。
  ◎自然蒸發,餘下:原本費洛蒙值,乘上1-α。

而目前的最佳路線,在蒸發之後,竟會神奇地額外增加一些費洛蒙。
  ◎最佳路線,添加:目前最佳解的數值,倒數後,再乘上α。

α是費洛蒙蒸發的程度,自行設定在0到1之間。

6.
2. 3. 4. 5.,重複執行R次。
途中不斷記錄最好的路線。
1. 初始化地圖與費洛蒙。
2. 以下動作執行R次:
  1. N隻螞蟻依序找路線,不是同時找路線:
    1. state transition rule,一隻螞蟻找一條路線。 
       此路線是由蟻窩出發,經過所有食物各一次,然後回到蟻窩。
    2. 記錄目前最佳路線。
    3. local updating rule,更新路線費洛蒙。
  2. global updating rule,更新所有道路費洛蒙。

演算法的過程,大致就是如此。細部的實作方式、參數的調校方式,隨人話虎爛。

螞蟻數量足夠豐富,可以避免進入區域極值。隨機探索用於增加答案豐富性,跳脫區域極值。

範例:Travelling Salesman Problem

N個城市,C(N,2)條路。選出N條。

[初始化]
答案設計成0到N-1的數字,代表一條路徑。
地圖上每個地點都有食物。
地圖上可以任意挑一地點作為蟻窩。

當答案具有「特定排列」的性質,就適合採用螞蟻演算法。

Dynamic Optimization

故兵無常勢,水無常形。能因敵變化而取勝者,謂之神。《孫子》

Dynamic Optimization

「動態最佳化」。函數時時變化,有如波浪,隨波逐流。

Plot3D[(Sin[x/2-y] + Sin[y/4])/2, {x, -4Pi, 4Pi}, {y, -4Pi, 4Pi}, PlotRange -> {-1, 1}, Boxed -> False, Mesh -> None, Axes -> False, PlotPoints -> 50, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-4, 4}]] &)]

Mean Shift(平均值移動)

一、均勻散布粒子。
二、以上次的極值座標為中心,
  取得鄰近範圍內所有粒子。(範圍自訂,通常是規定半徑長度。)
三、範圍內所有粒子各自求函數值。
四、範圍內所有粒子的座標,求加權平均(權重是函數值),得到新的極值座標。

另一種觀點是取樣。樣本(粒子)的密度,當作是函數值的大小。可以想成是找到粒子密度最高的地方。【待補文字】

Particle Filter(粒子濾波器)

一、以上次的極值座標為中心,
  散布n個粒子,呈高斯分布。(範圍自訂,範圍即變異數)
  粒子不用飛來飛去了,直接散布。
二、n個粒子各自求函數值。
三、n個粒子的座標,求加權平均(權重是函數值),得到新的極值座標。
四、或者不採用加權平均。
  想像n個粒子各有磁場(高斯分布),強弱由函數值大小決定,
  形成joint distribution(高斯混和分布)。
  下一回合依此分布散布粒子。

Follow The Leader(跟隨領導)

Ensemble Average改良版,適用於凸函數、凹函數。

前幾次的函數相加後,找到極值位置,推定為本次的極值位置、或者本次的登山起點。找極值位置時,只找幾個重點位置。

重點位置隨便找,指引我們向前進──儘管非常唬爛,但是效果出眾,於是有人認真進行數學分析:

http://courses.cs.washington.edu/courses/cse599s/14sp/scribes.html

還有其他進階版本FTRL、FTRL-Proximal。

Multiobjective Optimization

魚,我所欲也;熊掌,亦我所欲也。    

二者不可得兼,舍魚而取熊掌者也。《孟子》

Multiobjective Optimization

「多目標最佳化」。找到一個輸入,同時最佳化多個函數。

{ minimize f(x)
{ maximize g(x)
{ maximize h(x)

通常不存在如此完美的輸入,所以只好折衷。

各取所長,以求得更均衡的答案。

Scalarization(純量化)

多個函數,極值位置通常不相等。多個函數的加減乘除,極值位置毫無規律。我們難以制定任何策略。就連最佳化專家孟子也自認只能選一個函數來最佳化。這種情況下,大家只好尋求心靈勵志書籍、感恩師父感恩上人了。

{ minimize f(x)   Mencius bear paw
{ maximize g(x)  ------------------>  maximize g(x)
{ maximize h(x)

然而值得一提的是,凸函數相加之後,仍是凸函數;而且新極值位置,夾在原極值位置之間,不會歪太多。得到折衷方案:最佳化凸函數的和,找到極值位置。

{ minimize f(x)   f,g,h are convex
{ minimize g(x)  ------------------>  minimize f(x) + g(x) + h(x)
{ minimize h(x)

為了調整極值位置,可以添上權重。為了保持是凸函數,權重不能是負數。權重可以自由設定,沒有所謂最適當的權重。

注意到,除了權重,函數本身的斜率變化也會影響極值位置。即使權重均等,極值位置也通常不在正中央。因此這個手法是任憑感覺、碰碰運氣。

minimize α f(x) + β g(x) + h(x)  (α ≥ 0, β ≥ 0, γ ≥ 0)

衡量輕重、化做權重,故稱之為「純量化」。

Regularization(正則化)

純量化的進階應用。追加其他的最佳化目標。

minimize f(x) + α g(x) + β h(x) + ...   (α ≥ 0, β ≥ 0, ...)

f(x)是原本的目標,g(x) h(x) ...是追加的目標。

大家視問題需要,訂立適當的函數。例如x的長度平方越小越好g(x) = ‖x‖²;物理學的做功越少越好g(x) = F x。

追加凸函數,導致函數整體變窄、鞍點周遭變陡,梯度下降法更快找到極值、不易停在鞍點。

大家視問題需要,訂立適當的權重。權重越小,函數的影響力就越小,類似於加權平均的概念。

求極值,權重的絕對大小不重要,考慮相對大小即可。原始函數的權重定為1,節省一個權重數值。

明定準則、施予規範,故稱之為「正則化」。

Constrained Optimization

天將降大任於斯人也,必先苦其心志,勞其筋骨,餓其體膚,   

空乏其身,行拂亂其所為,所以動心忍性,曾益其所不能。《孟子》

Constrained Optimization

「約束最佳化」。最佳化,限制輸入範圍。

輸入只能是受規範的、被指定的數值。

maximize f(x)
 x ∈ C

「約束最佳化」。最佳化、等式與不等式組,兩者合體。

輸入必須同時滿足許多道等式與不等式。

maximize f(x)   subject to g(x) = 0, h(x) ≥ 0, ...

約束最佳化當中,函數稱做「目標函數objective function」。等式與不等式組稱做「約束條件constraint」。等式與不等式可以重新整理成函數的模樣,稱做「約束函數constraint function」。

通常會存在太過完美的輸入,所以只好抑制。

施加限制,以求得更對味的答案。

Show[ Plot3D[(Cos[(x-y)/2 - 3] + Cos[y/4])/2 - 1/2, {x, -4Pi, 4Pi}, {y, -4Pi, 4Pi}, PlotRange -> {-1, 1}, Boxed -> False, Mesh -> None, Axes -> False, PlotPoints -> 50, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-4, 4}]] &)], ParametricPlot3D[{{u,-4Pi,0}, {-4Pi,u,0}}, {u, -4Pi, 4Pi}, PlotStyle -> {Black, Black}] ]
Show[ Plot3D[BesselJ[0, Norm[{x, y}]], {x, -4Pi, 4Pi}, {y, -4Pi, 4Pi}, PlotRange -> {-1, 1}, Boxed -> False, Mesh -> None, Axes -> False, PlotPoints -> 50, ColorFunction -> (ColorData["GrayTones"][Rescale[#3, {-6, 2}]] &)], Plot3D[BesselJ[0, Norm[{x, y}]] + 0.01, {x, -4Pi, 4Pi}, {y, -4Pi, 4Pi}, PlotRange -> {-1, 1}, Boxed -> False, Mesh -> None, Axes -> False, PlotPoints -> 50, RegionFunction -> Function[{x, y, z}, Cos[(x-y)/2 - 3] + Cos[y/4] > 1], ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-4, 4}]] &)] ]
Show[ Plot3D[0, {x, -4Pi, 4Pi}, {y, -4Pi, 4Pi}, PlotRange -> {-1, 1}, Boxed -> False, Mesh -> None, Axes -> False, PlotPoints -> 50, ColorFunction -> (ColorData["GrayTones"][Rescale[#3, {-6, 2}]] &)], Plot3D[0.01, {x, -4Pi, 4Pi}, {y, -4Pi, 4Pi}, PlotRange -> {-1, 1}, Boxed -> False, Mesh -> None, Axes -> False, PlotPoints -> 50, RegionFunction -> Function[{x, y, z}, Cos[(x-y)/2 - 3] + Cos[y/4] > 1], ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-4, 4}]] &)] ]
ContourPlot[(Cos[(x-y)/2 - 3] + Cos[y/4])/2 - 1/2 == 0, {x, -4Pi, 4Pi}, {y, -4Pi, 4Pi}]

Conditional Gradient Method(條件梯度法)

梯度下降法改良版。總是朝向約束條件。

p = argmin [ ∇f(x) dot (p - x) ] = argmin [ ∇f(x) dot p ]
     p∈C                            p∈C

xnext = x + step ⋅ (p - x)

從約束條件當中,找出最符合當前梯度方向、最靠近當前位置的地點(以點積大小來判斷),然後朝該地點走一步。

找到該地點,是另一個最佳化問題。兩層最佳化。

如果約束函數很單純,則推導公式解,得到該地點。如果約束函數很複雜,則以梯度下降法,找到該地點。

Projected Gradient Method(投影梯度法)

梯度下降法改良版。每當偏離約束條件,立即歸位。

p = x + step ⋅ ∇f(x)

xnext = argmin ‖ p - xnext ‖
       xnext∈C

朝著梯度方向走一步,然後立即垂直投影到約束條件上面。換句話說,若離開約束條件,則走最短直線抵達約束條件;若未離開約束條件,則維持現狀。

垂直投影,是另一個最佳化問題。兩層最佳化。

如果約束函數很單純,則推導垂直投影公式,一步到位。如果約束函數很複雜,則以梯度下降法慢慢走入約束條件:約束函數的梯度方向,權且做為垂直投影方向。

Lagrange Multiplier(拉格朗日乘數)

當約束條件是等式,限制最佳化可以化作微分方程組求解。

f(x,y) 求極值。必須滿足 g(x,y) = 0。

湊得 f(x,y) = f(x,y) + λ g(x,y)
定義 s(x,y,λ) = f(x,y) + λ g(x,y)

f(x,y) 的極值,等同 s(x,y,λ) = f(x,y) + λ g(x,y) 的極值。
欲求極值:
對 x 偏微分,讓斜率是 0。
對 y 偏微分,讓斜率是 0。

不管 λ 如何變化,λ g(x,y) 都是零,s(x,y,λ) 永遠不變。
欲求永遠不變的地方:
對 λ 偏微分,讓斜率是 0。

三道偏微分方程式聯立之後,其解涵蓋了(不全是)所有符合約束條件的極值。
{ ∂/∂x s(x,y,λ) = 0
{ ∂/∂y s(x,y,λ) = 0
{ ∂/∂λ s(x,y,λ) = 0

微分方程組的解,必定涵蓋所有正確答案,卻也可能包含錯誤答案。事後需要驗證答案,比對一下鄰近函數值即可。

可惜的是,微分方程組求解,沒有簡單的演算法。大家都是手工推導公式解。

另外還可以進一步得到非常漂亮的數學性質:極值所在之處,目標函數與約束函數,兩者梯度方向共線,而梯度大小不明。(梯度大小設定成未知倍率λ,稱作拉格朗日乘數。)

三式分別偏微分。
{ ∂/∂x s(x,y,λ) = 0             { ∂/∂x f(x,y) + λ ⋅ ∂/∂x g(x,y) = 0
{ ∂/∂y s(x,y,λ) = 0   ------>   { ∂/∂y f(x,y) + λ ⋅ ∂/∂y g(x,y) = 0
{ ∂/∂λ s(x,y,λ) = 0             { g(x,y) = 0

一式、二式合起來,即是梯度。
∇f(x,y) + λ ∇g(x,y) = 0
∇f(x,y) = -λ ∇g(x,y)
∇f(x,y) = λ ∇g(x,y)      變數替換 λ := -λ

如果不喜歡變數替換,那麼修改定義。
s(x,y,λ) = f(x,y) - λ g(x,y)

另外還可以進一步整理成「梯度等於零」、「求駐點」的格式。

∇f(x,y) + λ ∇g(x,y) = 0
∇[f(x) + λ g(x,y)] = 0       求駐點。先微後加、先加後微,結果相同。
optimize [f(x) + λ g(x,y)]   ✘ 注意到,求駐點,不等價於求極值。

當約束條件是大量等式,則通通乘上倍率再加總。

f(x,y) 求極值。必須滿足 g₁(x,y) = 0 與 g₂(x,y) = 0 與……。
定義 s(x,y,λ₁,λ₂,...) = f(x,y) + λ₁ g₁(x,y) + λ₂ g₂(x,y) + ...

Karush-Kuhn-Tucker Conditions(KKT條件)

約束條件是等式:稱作拉格朗日乘數。極值位於約束條件上面。

極值所在之處,目標函數與約束函數,兩者梯度方向共線,而梯度大小不明。(梯度大小設定成未知倍率λ。)

maximize f(x)   subject to g(x) = 0

---> solve { ∇f(x) = λ ∇g(x)
           { g(x) = 0

約束條件是不等式:稱做KKT條件。極值可能位於約束條件邊界或者內部。

邊界:目標函數與約束函數,兩者梯度方向共線,而梯度大小不明。(梯度大小設定成未知倍率μ。)

內部:目標函數的極值位置,已經滿足約束條件。約束條件無關緊要。(μ設定成零,讓約束條件失效。)

計算學家將兩種情況分開處理,邊界是微分方程組求解,內部是目標函數最佳化。數學家則是將兩種情況合併成一個方程組。

maximize f(x)   subject to h(x) ≥ 0

---> solve { ∇f(x) = μ ∇h(x)
           { μ ≤ 0
           { h(x) ≥ 0

極值位於約束條件邊界時,梯度方向的判斷方式:

一、最大值位於邊界:目標函數的梯度方向,離開邊界。

二、約束函數大於等於零:約束函數的梯度方向,進入邊界。

梯度方向恰好相反,梯度大小必須異號。(μ必須小於零。)

                      | ∇f(x) = μ ∇h(x) | ∇f(x) + μ ∇h(x) = 0
----------------------|------------------|---------------------
max f(x) st h(x) ≥ 0  | μ ≤ 0            | μ ≥ 0
max f(x) st h(x) ≤ 0  | μ ≥ 0            | μ ≤ 0
min f(x) st h(x) ≥ 0  | μ ≥ 0            | μ ≤ 0
min f(x) st h(x) ≤ 0  | μ ≤ 0            | μ ≥ 0

約束條件是大量等式與不等式:各有倍率,通通加總。

maximize f(x)   subject to g(x) = 0, h(x) ≥ 0

---> solve { ∇f(x) = λ ∇g(x) + μ ∇h(x)
           { g(x) = 0
           { μ ≤ 0
           { h(x) ≥ 0

Regularization與KKT Conditions沒有任何關聯,只不過是數學式子有點像。別忘了KKT Conditions還得聯立其他式子。

Regularization:
minimize f(x) + α g(x) + β h(x) + ...   (α ≥ 0, β ≥ 0, ...)

KKT Conditions:
∇[f(x) + λ g(x) + μ h(x) + ...] = 0     (μ ≥ 0, ...)

最佳化與物理現象

物理  丨數學
一一一一十一一一一一一一一一一
力   丨向量         
合力  丨線性組合(加權總和) 如果有很多力,可以直接乘上倍率再相加。
力>位能丨積分         力對位置積分=位能。
位能>力丨微分         位能對位置微分=力。
位能  丨半正定函數(拋物線) 屬於凸函數,有唯一最小值。
物理  丨數學    丨位能         力        
一一一一十一一一一一一十一一一一一一一一一一一一一一一一一一一一
彈力  丨微分(梯度)丨0.5 k (x - l)²  丨k (x - l)
約束力 丨微分(梯度)丨0.5 k c(x)²     丨k c(x) c′(x)
慣性運動丨梯度等於零 丨min 位能       丨力 = 0
靜力平衡丨梯度等於零 丨min 位能總和     丨合力 = 0
施加外力丨拉格朗日乘數丨min 位能 st 約束能丨力 + λ ⋅ 約束力 = 0

例如彈簧系統。統計所有彈簧的彈力位能,求最小值,求最小值所在位置(每個彈簧的端點座標)。梯度等於零,就是靜力平衡,就是穩定狀態。實際應用諸如建築結構分析、化學結構分析。

例如軌道運動。物體位置限制在軌道上面,動能與位能的總和是常數。根據拉格朗日乘數,限制能量大小,等同施加外力;此例當中就是施加向心力。實際應用諸如太空衛星。