跳到主要內容

[最佳化理論] Conjugate Direction Methods (2) - The Conjugate Gradient Algorithm for Quadratic Objective function

OK,現在接續前面兩篇
1. [最佳化理論] Conjugate Direction Methods (0) -Theory
2. [最佳化理論] Conjugate Direction Methods (1) - Basic Algorithm for Quadratic Objective function

這次要來解決如何找到 Conjugate Direction的辦法。幫助我們找到 Conjugate Direction 的方法稱作 共軛梯度演算法 (Conjugate Gradient algorithm) 。

此演算法 藉由 梯度的幫助,使得在任意跌代步驟中,Conjugate Direction 可以由 前一個跌代步的方向 與 現在的 梯度做線性組合來計算出來。且用此法所產生的 Direction 可以保證是 每個方向都是互為 Q-conjugate。故名為 Conjugate Gradient Algorithm;最後我們會給出一個以 二階具有 $n=3$ 個變數的目標函數 作為例子來展示這套方法。

===============

如前所述,考慮標準二階目標函數
\[
J(u) = \frac{1}{2} u^T Q u - u^T b, \ u \in \mathbb{R}^n
\]
其中 $Q = Q^T >0$。

現在我們來看看Conjugate Gradient Algorithm 是如何找到 Conjugate Direction:

對 初始跌代步:$k=0$。
給定任意初始值 $x^{(0)}$ ,且設定 初始方向為最陡坡度(steepest descent)方向,亦即

$d^{(0)} = - \bigtriangledown J(u^{(0)})$

因此 $ u^{(1)} = u^{(0)} + \alpha_0 d^{(0)}$

其中 $\alpha_0$ 一般而言是由line search (一般而言使用MATLAB fminsearch 指令較為簡便)得到,亦即 $\alpha_0 = \arg \min J(x^{(0)} + \alpha d^{(0)})$
不過如果是 考慮 二階目標函數中,則我們可以求得精確的 $\alpha_k$ 如下式
\[
\alpha_0 = \arg \min J(x^{(0)} + \alpha d^{(0)}) =  - \frac{{\bigtriangledown^T J({u^{(k)}}){d^{(k)}}}}{{{d^{(k)T}}Q{d^{(k)}}}}
\] 在算出 $ u^{(1)} $ 之後,接著我們移動到下一個跌代步。
因為我們要找到 $d^{(1)}$ 且 此方向與 $d^{(0)}$ 互為Q-conjugate。所以如前所述我們選擇
\[
d^{(1)} = \text {linear combination of } \bigtriangledown J(u^{(1)}) \ \& \ d^{(0)}
\]亦即,對 $k+1$ 跌代步的時候 我們選
\[ d^{(k+1)} = - \bigtriangledown J(u^{(k+1)}) + \beta_k \cdot d^{(k)} \ \ \ ,k=0,1,2,... \]
其中 $\beta_k = - \frac{{\bigtriangledown^T J({u^{(k+1)}}){d^{(k)}}}}{{{d^{(k)T}}Q{d^{(k)}}}}$

=========================
我們將 Conjugate Gradient Algorithm (for Quadratic Objective Function) 總結如下8個步驟:

  1. 在 $k=0$ 時候給定初始值 $u^{(0)}$
  2. 計算初始梯度 $ \bigtriangledown J(u^{(0)})$ ,且設令其為初始方向 $d^{(0)} = - \bigtriangledown J(u^{(0)})$ :注意,如果 $\bigtriangledown J(u^{(0)}) = 0$ 則跌代停止
  3. 決定 $\alpha_k = - \frac{{\bigtriangledown^T J({u^{(k)}}){d^{(k)}}}}{{{d^{(k)T}}Q{d^{(k)}}}}$
  4. 計算 $ u^{(k+1)} = u^{(k)} + \alpha_k d^{(k)}$
  5. 計算下一個跌代的梯度:  $ \bigtriangledown J(u^{(k+1)})$;如果 $\bigtriangledown J(u^{(0)}) = 0$ 則跌代停止
  6. 計算 $\beta_k = - \frac{{\bigtriangledown^T J({u^{(k+1)}}){d^{(k)}}}}{{{d^{(k)T}}Q{d^{(k)}}}}$
  7. 計算下一個Conjugate 方向 $d^{(k+1)} = - \bigtriangledown J(u^{(k+1)}) + \beta_k \cdot d^{(k)}$
  8. 令 $k=k+1$ 並重覆到step 3

下面我們將給個例子來驗證上述的方法。

=================

Example:
考慮下列二階 目標函數
\[
J({u_1},{u_2},{u_3}) = \frac{3}{2}u_1^2 + 2u_2^2 + \frac{3}{2}u_3^2 + u_1u_3 + 2u_2u_3 - 3{u_1} - {u_3}
\] 令初始值為 $u^{(0) } = [0, \ 0, \ 0]^T$ 求使用 Conjugate Gradient Method 找出最佳解。

Solution

首先觀察二階 目標函數,我們可以將其改寫成矩陣型式如下

$J(u) = \frac{1}{2} u^T Q u - u^T b$

$\Rightarrow J(u)= \frac{1}{2}\left[ {\begin{array}{*{20}{c}}{{u_1}}&{{u_2}}&{{u_3}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}3&0&1\\0&4&2\\1&2&3\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{u_1}}\\{{u_2}}\\{{u_3}}\end{array}} \right] - \left[ {\begin{array}{*{20}{c}}{{u_1}}&{{u_2}}&{{u_3}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}3\\0\\1\end{array}} \right]$

接著我們計算梯度:

$ \bigtriangledown J(u) = Q u - b = \left[ {\begin{array}{*{20}{c}}{3{u_1} + {u_3} - 3}\\{4{u_2} + 2{u_3}}\\{{u_1} + 2{u_2} + 3{u_3} - 1}\end{array}} \right]$

在使用 Conjugate Gradient Method之前我們可以先看看其最佳解落在哪裡:
由一階必要條件  $\bigtriangledown J(u^*)=0 $ 與二階充分條件 $\bigtriangledown^2 J(u^*) >0 $ 可知最佳解為

 $u^* = Q^{-1}b = [1, 0, 0]^T$

現在我們開始用 Conjugate Gradient Method 來看看會
現在由步驟2可知初始梯度為

$ \bigtriangledown J(u^{(0)}) = Q u^{(0)} - b = \left[ {\begin{array}{*{20}{c}}{ - 3}\\{0}\\{ - 1}\end{array}} \right]$

故 初始Conjugate direction

$d^{(0)} = -\bigtriangledown J(u^{(0)}) =\left[ {\begin{array}{*{20}{c}}{  3}\\{0}\\{  1}\end{array}} \right] $

接著由步驟3 我們可以計算 初始 $\alpha_k$

 $\alpha_0 = - \frac{{\bigtriangledown^T J({u^{(0)}}){d^{(0)}}}}{{{d^{(0)T}}Q{d^{(0)}}}} = \frac{10}{36} = 0.2778$

接著由步驟4 : $ u^{(k+1)} = u^{(k)} + \alpha_k d^{(k)}$ 可以算出

$ u^{(1)} = u^{(0)} + \alpha_0 d^{(0)} = [0.8333,\ 0, \ 0.2778]^T$

接著由步驟5:計算下一個跌代的梯度:  $ \bigtriangledown J(u^{(k+1)})$

$ \bigtriangledown J(u^{(1)}) = Q u^{(1)} - b =  [-0.2222,\ 0.5566, \ 0.6667]^T$

接著由步驟6:計算  $\beta_k $

$\beta_0 = - \frac{{\bigtriangledown^T J({u^{(1)}}){d^{(0)}}}}{{{d^{(0)T}}Q{d^{(0)}}}} = 0.08025$

再者由步驟7 :計算下一個Conjugate 方向 $d^{(k+1)} $

 $d^{(1)} = - \bigtriangledown J(u^{(1)}) + \beta_0 \cdot d^{(0)} = [0.4630, -0.5556, -0.5864]^T$

再來回到步驟3 重複前述步驟

最後會得到 $u^{(3)} = u^{(2)} + \alpha_2 d^{(2)} = [1, \ 0, \ 0]^T$ 且 $\bigtriangledown J(u^{(3)}) = 0$

上述表示在第三步的時候收斂到最佳解。此結果確實說明了 Conjugate Direction Method 保證對二階 具有 n變數的 目標函數 能夠在 n次跌代之內收斂。此例 n=3 故在第3次跌代的時候即收斂到最佳解。

ref: E. K. P. Chong, S. H. Zak, An Introduction to Optimization 2nd.

---
延伸閱讀

留言

這個網誌中的熱門文章

[數學分析] 什麼是若且唯若 "if and only if"

數學上的 if and only if  ( 此文不討論邏輯學中的 if and only if,只討論數學上的 if and only if。) 中文翻譯叫做  若且唯若 (or 當且僅當) , 記得當初剛接觸這個詞彙的時候,我是完全不明白到底是甚麼意思,查了翻譯也是愛莫能助,畢竟有翻跟沒翻一樣,都是有看沒有懂。 在數學上如果看到 if and only if  這類的句子,其實是表示一種 雙條件句 ,通常可以直接將其視為" 定義(Definition)" 待之,今天要分享的是這樣的一個句子如何用比較直觀的方法去看他 假設我們現在有 兩個邏輯陳述句 A 與  B. 注意到,在此我們不必考慮這兩個陳述句到底是什麼,想表達什麼,或者到底是否為真(true),這些都不重要。只要知道是兩個陳述即可。 現在,考慮新的陳述:  "A if and only if B" 好了,現在主角登場,我們可以怎麼看待這個句子呢? 事實上我們可以很直覺的把這句子拆成兩部分看待,也就是 "( A if B ) and ( A only if B )" 那麼先針對第一個部分  A if B  來看, 其實這句就是說  if B then A, 更直白一點就是 "if B is true, then A is also true".  在數學上等價可以寫為 "B implies A" .  或者更常用一個箭頭符號來表示 "B $\Rightarrow$  A"  現在針對第二個部分  A only if B 此句意指  "If B is not true, then A is also not true". 所以如果已知 A is true,  那麼按照上句不難推得 B is also true 也就是說  A only if B  等價為 "If A is true then B is also true". 同樣,也可以寫作   "A implies B"   或者用箭頭表示  "A   $\Rightarrow$     B".

[數學分析] 淺談各種基本範數 (Norm)

這次要介紹的是數學上一個重要的概念: Norm: 一般翻譯成 範數 (在英語中 norm 有規範的意思,比如我們說normalization就是把某種東西/物品/事件 做 正規化,也就是加上規範使其正常化),不過個人認為其實翻譯成 範數 也是看不懂的...這邊建議把 Norm 想成長度就好 (事實上norm是長度的抽象推廣), 也許讀者會認為好端端的長度不用,為何又要發明一個 norm 來自討苦吃?? 既抽象又艱澀。 事實上想法是這樣的: 比如說現在想要比較兩個數字 $3$ , $5$ 之間的大小,則我們可以馬上知道 $ 3 < 5 $;同樣的,如果再考慮小數與無理數如 $1.8753$ 與 $\pi$,我們仍然可以比較大小 $1.8753 < \pi = 3.1415...$ 故可以發現我們有辦法對 "純量" 做明確的比大小,WHY? 因為前述例子中 $3$, $5$, $1.8753$ or $\pi$ 其各自的大小有辦法被 "measure "! 但是如果是現在考慮的是一組數字 我們如何去measure 其大小呢?? 比如說 \[x:=[1, -2, 0.1, 0 ]^T \]上式的大小該是多少? 是 $1$? $-2$? $0.1$??? 再者如果更過分一點,我們考慮一個矩陣 \[A = \left[ {\begin{array}{*{20}{c}} 1&2\\ 3&4 \end{array}} \right] \],想要知道這個矩陣的大小又該怎麼辦?? 是 $1$ ? $2$ 還是 $4$ ?..其實現階段我們說不清楚。 也正是如此,可以發現我們確實需要新的 "長度" 的定義來幫助我們如何去 measure 矩陣/向量/甚至是函數的大小。 故此,我們首先定義甚麼是Norm,(也就是把 "長度" or "大小" 的本質抽離出來) ================== Definition: Norm 考慮 $V$ 為一個向量空間(Vector space),則我們說  Norm 為一個函數 $||\cdot|| : V \rightarrow \mathbb{R}$ 且滿足下列性質