[最佳化] 淺談 Steepest Descent Method (1) - Optimal step size

延續上篇文章 [最佳化] 淺談 Steepest Descent Method (0) -Why Steepest !?，這次我們要介紹 Steepest Descent Method with Optimal Step size

修訂後的Steepest Descent Algorithm 需要甚麼呢?

1. 初始條件 $u^0$
2. 最大跌代步長上限(Maximum fixed step size): $H$
3. Steepest Descent with Optimal Step size的跌代架構(iterative scheme) \[u^{k+1} = u^k - h_k \frac{ \nabla J(u^k)}{|| \nabla J(u^k)||} \ \ \ \ (*)\]
4. 演算法停止判別機制(stopping criterion) : EX: \[ ||\nabla J(u^k)|| < \varepsilon\]

那麼問題變成 $h_k$ 該怎麼求?

首先我們考慮第 $k$次 跌代，手上有 $u^k$，則我們可以定義

\[
\tilde J(h) := J(u^k - h \frac{ \nabla J(u^k)}{|| \nabla J(u^k)||})
\]接著我們做 Line search 找出一個最佳的 $h \in [0,H] $ ( 亦即 $\min \tilde J(h)$) 把此 $h$ 稱做 $h_k $，也就是說

\[
\tilde J(h_k) = \displaystyle \min_{h \in [0, H]} \tilde J(h)
\]做Line search之後得到的 $h_k$ 再把他放回 $(*)$ 即可!!

\[
u^{k+1} = u^k - h_k \frac{ \nabla J(u^k)}{|| \nabla J(u^k)||} \ \ \ \ (*)
\]上式即稱為 Steepest Descent with Optimal Step size (此Optimal step 由Line search 對 $ \min J(u^k - h \frac{ \nabla J(u^k)}{|| \nabla J(u^k)||})$ 求得)

對於Steepest Descent Algorithm而言，我們有 現在的跌代步的梯度 與下一個跌代步的梯度互為垂直；亦即
\[
\left ( \nabla J(u^k) \right )^T \cdot \nabla J(u^{k+1}) =0
\]現在如果我們考慮更一般的情況，

===============================
Theorem: (Optimal Descent Condition)
考慮 $v \in \mathbb{R}^n$ 為某一個方向 (不必是梯度)，且假設 $h_k$ 把 $\tilde J(u^k + h \cdot v)$最小化，且我們的跌代式為
\[
u^{k+1} = u^k + h_k \cdot v
\]則我們有 $v^T \cdot \nabla J(u^{k+1}) =0$，我們稱此條件為 Optimal Descent Conditon。
===============================

Proof:
我們欲證  $v^T \cdot \nabla J(u^{k+1}) =0$，

故由  $h_k$ 把 $\tilde J(u^k + h \cdot v)$最小化 的假設，我們已知一階必要條件成立，故可利用一階必要條件FONC
\[\frac{{\partial \tilde J({u^k} + h \cdot v)}}{{\partial h}} = 0
\]為了方便起見，現在令 $z: = {u^k} + h \cdot v$，則上式可推得
\[\begin{array}{l}
\frac{{d\tilde J(z)}}{{dh}} = \frac{{\partial \tilde J(z)}}{{\partial {z_1}}}\frac{{\partial {z_1}}}{{\partial h}} + \frac{{\partial \tilde J(z)}}{{\partial {z_2}}}\frac{{\partial {z_2}}}{{\partial h}}... + \frac{{\partial \tilde J(z)}}{{\partial {z_n}}}\frac{{\partial {z_n}}}{{\partial h}} = 0\\
 \Rightarrow \frac{{d\tilde J(z)}}{{dh}} = \frac{{\partial \tilde J(z)}}{{\partial {z_1}}}{v_1} + \frac{{\partial \tilde J(z)}}{{\partial {z_2}}}{v_2}... + \frac{{\partial \tilde J(z)}}{{\partial {z_n}}}{v_n} = 0\\
 \Rightarrow \frac{{d\tilde J(z)}}{{dh}} = \left[ {\begin{array}{*{20}{c}}
{{v_1}}&{{v_2}}& \cdots &{{v_n}}
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{\frac{{\partial \tilde J(z)}}{{\partial {z_1}}}}\\
{\frac{{\partial \tilde J(z)}}{{\partial {z_2}}}}\\
 \vdots \\
{\frac{{\partial \tilde J(z)}}{{\partial {z_n}}}}
\end{array}} \right] = 0\\
 \Rightarrow \frac{{d\tilde J(z)}}{{dh}} = {v^T}\nabla \tilde J(z) = {v^T}\nabla \tilde J({u^k} + h \cdot v) = {v^T}\nabla \tilde J({u^{k + 1}}) = 0
\end{array}
\]最後一個等式成立是由於 $h$ 最小化 $\tilde J(u^k + h \cdot v)$，故稱此$h=h_k$，又由跌代式的假設 $u^{k+1} = u^k + h_k \cdot v$。 $\square$


那麼現在我們來看看如果目標函數是標準二次的情況

考慮如下標準二次函數
\[
J(u) = u^T A u + b^T u + c, \ A=A^T, \ A \succ 0
\]注意到上述二次函數可以直接用 一階必要條件FONC ( $\nabla J(u^k) =0$) 與 二階充分條件SOSC ($\nabla^2 J(u^k) \succ 0$) 直接求解，可得最佳解為 $u^* = \frac{1}{2} A^{-1}b$

那麼如果我們現在採用 Steepest Descent Algorithm with Optimal Step size $h_k$ ，我們想要知道選怎樣的 $h_k$ 可以得到同樣的最佳解呢?

故首先給定初始條件 $u^0$, 且給定足夠大的步長上限 $H$，

接著我們寫下 Line Search需要的式子
\[
\tilde J (h) := J(u^k  - h \cdot \nabla J(u^k)) \ \ \ \ (*)
\] ，目標是要找出 $h =?$

首先觀察 $(*)$，我們可以先計算上式的梯度部分，由 FONC 可知梯度為
\[
\nabla J(u^k) = 2 A u^k + b
\]故將其代入 $(*)$ 可得
\[
\tilde J(h): = J\left( {{u^k} - h\left( {2A{u^k} + b} \right)} \right)
\]又因為
\[
J(u) = u^T A u + b^T u + c
\]故可推得
\[\begin{array}{l}
 \Rightarrow \tilde J(h) = {\left( {{u^k} - h\left( {2A{u^k} + b} \right)} \right)^T}A\left( {{u^k} - h\left( {2A{u^k} + b} \right)} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} + {b^T}\left( {{u^k} - h\left( {2A{u^k} + b} \right)} \right) + c
\end{array}
\]由於 $h$ 為最小化 $\tilde J(h)$ 故由FONC對 $h$ 可知 $\frac{d \tilde J(h)}{dh} =0$ 亦即
\[\begin{array}{l}
 \Rightarrow \frac{{d\tilde J(h)}}{{dh}} = 0\\
 \Rightarrow  - {u^k}^TA\left( {2A{u^k} + b} \right) - {\left( {2A{u^k} + b} \right)^T}A{u^k}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} + 2h{\left( {2A{u^k} + b} \right)^T}A\left( {2A{u^k} + b} \right) - {b^T}\left( {2A{u^k} + b} \right) = 0\\
 \Rightarrow {h_k}: = h = \frac{1}{2}\frac{{{{\left( {2A{u^k} + b} \right)}^T}\left( {2A{u^k} + b} \right)}}{{{{\left( {2A{u^k} + b} \right)}^T}A\left( {2A{u^k} + b} \right)}}
\end{array}
\]
Comments:
1. 注意到上式中分母為 ${{{\left( {2A{u^k} + b} \right)}^T}A\left( {2A{u^k} + b} \right)}$ 為 $1 \times 1$ 此時分母不再是矩陣或者向量，故可以直接執行除法。且由於我們的假設 $A$ 矩陣為正定矩陣，亦即 $x^T A x >0, \forall x \neq 0$，仔細觀察上式分母，若令 $x:={\left( {2A{u^k} + b} \right)}$，我們確實得到
\[
{{{\left( {2A{u^k} + b} \right)}^T}A\left( {2A{u^k} + b} \right)} = x^T A x >0
\]

[最佳化] 淺談 Steepest Descent Method (0) -Why Steepest !?

這次要介紹的是最陡坡度法(Steepest Descent Method)，又稱 Gradient descent method：

想法：透過負梯度 (negative gradient) 作為最陡坡度，逐步找到 (局部)最小值 (最佳解 $u^*$)

這個演算法需要甚麼呢?

1. 初始條件 $u^0$
2. 固定的跌代步長(fixed step size): $h$
3. Steepest Descent 的跌代架構(iterative scheme) \[u^{k+1} = u^k - h \frac{ \nabla J(u^k)}{|| \nabla J(u^k)||}\]
4. 演算法停止判別機制(stopping criterion) : EX: 給定誤差 $\varepsilon>0$，檢驗 \[ ||\nabla J(u^k)|| < \varepsilon\]

那麼現在我們來解決一個問題：

為什麼此法被稱作 "最陡" 坡度? 
也就是說 為什麼Iterative scheme 中的方向 $ \nabla J(u^k)$ 被稱做是最陡(Steepest)方向??

考慮目標函數 $J: \mathbb{R}^n \rightarrow \mathbb{R}$，其在某點 $u^0 \in \mathbb{R}$ 與方向$v$ 的方向導數(Directional derivative at point $u^0$ in direction $v$)定義如下：

\[
{\left. {\frac{{\partial J\left( u \right)}}{{\partial v}}} \right|_{u = {u^0}}}: = {\left[ {\nabla J({u^0})} \right]^T} \cdot \frac{v}{{\left\| v \right\|}}
\]由Cauchy-Schwarz inequality $\left| {{x^T}y} \right| \le \left\| x \right\|\left\| y \right\|$，可推得上式如下：
\[
\left| {{{\left[ {\nabla J({u^0})} \right]}^T} \cdot \frac{v}{{\left\| v \right\|}}} \right| \le \left\| {\nabla J({u^0})} \right\|\frac{{\left\| v \right\|}}{{\left\| v \right\|}}
\]現在如果我們把方向 $v$ 選成梯度方向，亦即

\[
v = \nabla J(u^0)
\]則可發現上述不等式變成
\[\begin{array}{l}
 \Rightarrow \left| {{{\left[ {\nabla J({u^0})} \right]}^T} \cdot \frac{{\nabla J({u^0})}}{{\left\| {\nabla J({u^0})} \right\|}}} \right| \le \left\| {\nabla J({u^0})} \right\|\\
 \Rightarrow \frac{{{{\left\| {\nabla J({u^0})} \right\|}^2}}}{{\left\| {\nabla J({u^0})} \right\|}} \le \left\| {\nabla J({u^0})} \right\|\\
 \Rightarrow \left\| {\nabla J({u^0})} \right\| = \left\| {\nabla J({u^0})} \right\|
\end{array}
\]故可知當我們選 $v = \nabla J(u^0) $ Cauchy-Schwarz inequality 的 "等" 式成立，故 $\nabla J(u^0) $ 為使方向導數最大的值! 亦即 最陡方向(Steepest)!

至於為什麼我們說Steepest Descent (最陡坡度下降)，是因為注意到Steepest Descent 演算法中跌代式子

\[u^{k+1} = u^k - h \frac{ \nabla J(u^k)}{|| \nabla J(u^k)||}
\] 的方向為負，亦即 $- \nabla J(u^k)$ !!，故我們是朝著最陡的方向往下逐步得到最佳解(最小值) $u^*$

現在我們給個例子：

Example
考慮如下目標函數
\[
J(u) = {\left( {11 - {u_1} - {u_2}} \right)^2} + {\left( {1 + 10{u_2} + {u_1} - {u_1}{u_2}} \right)^2}
\] 1. 試證上述目標函數最佳解為 $[13, 4]^T$
2. 繪製其 $0 \le u_1 \le 20$ 與 $ 0 \le u_2 \le 15$ 的範圍
3. 給定初始值 $u^0 = [8, 12]^T $使用上述 Steepest Descent algorithm 與不同的固定步長 $h=0.01, 1.0$看看發生甚麼事情

Solution:
1. 透過一階必要條件(FONC) 與 二階充分條件(SOSC)即可求得最佳解 $u^* = [13,
4]^T$。在此不贅述。

2. 透過MATLAB contour 指令可以繪製目標函數的等高線圖如下

3. 考慮$u^0 = [8, 12]^T$ ，並考慮 $h=0.01$的情況，透過MATLAB實現上述Steepest Descent Algorithm並限制停止判別為跌代步驟不超過兩千步 $k_{max} := 2,000$。
在約 $k=1000$ 步之後，可收斂到 $u = [13.01, 3.993]$。如下圖所示： (點圖放大)；

若現在考慮 $h=1.0$ 則Steepest Descent 展示了Zig-Zag的現象，最終落在 $u=[12.6, 3.601] to [13.31, 4.08]$之間，且跌代步如下圖所示 (點圖放大)

注意到上述例子中，對於較大的固定步e.g., $h=1.0$ ，Steepest Descent 表現出來回震盪的情況，對於較小的固定步 e.g., $h=0.01$，Steepest Descent 收斂緩慢 (超過一千步才收斂)。

Summary: 
Steepest Descent Algorithm 雖然想法很直覺，但事實上本質有兩個重大的缺點：
1. 注意到Steepest Descent的跌代式子中除了要計算梯度之後，仍需要給定固定步長 $h$，如果固定步長 $h$ 太大!，則會演算法產生衝過頭的情形。也就是說假設 $h=100$ (100步長單位) 當我可能今天只需要1步就到達最佳解，但Steepest Descent Method卻被迫每次都要走步長單位為 $100$ ，則永遠只能在最佳解附近震盪永遠到不了最佳解，

2. 如果固定步長 $h$ 太小，則雖然在某種程度上解決了震盪問題，但此時收斂速度會變得異常緩慢。

如何解決上述的問題!??
我們需要徹底地拔除固定步長的限制，此法稱做 Steepest Descent with Optimal Step Algorithm。
亦即我們將原本的Steepest Descent Algorithm的跌代式中的固定步長 $h$改成動態步長 $h_k$
\[
u^{k+1} = u^k - h_k \frac{\nabla J(u^k)}{|| \nabla J(u^k)||}
\]至於此 $h_k$該怎麼選? 有興趣的讀者可以參考下篇
[最佳化] 淺談 Steepest Descent Method (1) - Optimal step size