2009年4月30日 星期四

[機率論] 期望值 與 Lebesgue 積分

這次要介紹機率論中一個重要的概念:期望值 (Expectation),本質上期望值被視為一個 Lebesgue 積分。更進一步地說就是在較抽象的高等機率論中, 期望值被定義為對某機率測度 (Probability measure, $P$ ) 的 Lebesgue 積分 。亦即 考慮機率空間為 $(\Omega, \mathcal{F}, P)$,則 $E[X]$ 具有如下形式:
\[
E[X] := \int_{\Omega} X d P  = \int_\Omega X(\omega) P(d\omega)
\]

我們由 Simple Function 出發 逐步建構 Lebesgue integral:

Step 1: Simple function 的期望值:

首先,我們定義 $X$ 為一個 simple function。亦即此函數可由 可測集合 (measurable sets $A_i$ ) 的 Indicator function 所組成。我們將其寫作是 Finite sum 如下:
\[
X = \displaystyle \sum_{i=1}^{n} c_i 1_{A_i}
\]其中 $c_i \in \mathbb{R}$, $A_i \in \mathcal{F}$ 且
\[{1_A}\left( x \right): = \left\{ \begin{array}{l}
1,\begin{array}{*{20}{c}}
{}
\end{array}if\begin{array}{*{20}{c}}
{}
\end{array}x \in A\\
0,\begin{array}{*{20}{c}}
{}
\end{array}if\begin{array}{*{20}{c}}
{}
\end{array}x \notin A
\end{array} \right.\]則我們定義 對此 Simple function $X$ 的期望值 (或稱對此 Simple function 的 Integral)為
\[
E[X] := \int X dP := \sum_{i=1}^n c_i P( \{A_i\} )
\]接著,我們定義對 非負隨機變數 的期望值:


Step 2: 非負隨機變數 的期望值:

對任意 非負隨機變數 $Y$ $(Y : \Omega \rightarrow [0,1])$,我們定義 $Y$ 的期望值為上述 simple function 的期望值的 supremum,亦即
\[
E[Y] := \sup \{E[X]: 0 \leq X \leq Y \ \text{with $X$ a simple function}  \}
\] 注意到 $E[X] := \int X dP := \sum_{i=1}^n c_i P(\{A_i\}) $


有了上述定義,我們可以拓展到一般的隨機變數:


Step 3: 一般隨機變數的期望值:

最後,對一般的 隨機變數的情況,該如何定義其期望值呢?

想法如下:
透過將 一般隨機變數 轉換為 上述 非負隨機變數,再用上述非負隨機變數所定義的期望值 即可。

故現在考慮 $Y$ 為任意隨機變數,我們引入兩個 新的非負隨機變數 $Y^+, Y^- \geq 0$ (此法類似線性規劃中,將最佳化問題改寫成標準型式所加入的 slack variable: 有興趣的讀者請參考: [最佳化] 淺談線性規劃(0)- Standard form of Linear Programming ):
\[
Y^+ := Y \cdot 1_{ \{ Y \geq 0 \} }\\
Y^- := -Y \cdot 1_{ \{ Y \leq 0 \}}
\]那麼現在觀察上述 新引入的隨機變數,我們可知其與原本任意隨機變數有如下關係:
\[
Y=Y^+ - Y^-
\]且上述新的隨機變數 (由於 $Y^+, Y^- \geq 0$),故其期望值可由先前所定義的 非負隨機變數期望值定義而得,故 $E[Y^+], E[Y^-]$ 為 well-defined。有了這些結果我們可以定義對一般隨機變數 $Y$ 的期望值為:
\[
E[Y] := E[Y^+] - E[Y^-] \ \ \ \ (*)
\]上式為 well-defined 若 $E[Y^+], E[Y^-]$ 並非全為 $\infty$ (亦即只要不要發生 $\infty - \infty$ 的情況,則上述定法的期望值 $E[Y]$ 都是 well-defined。 )

Comment:
1. 上式 $(*)$ 與下列積分等價:
\[
\int Y dP := \int Y^+ dP - \int Y^- dP
\]
2. 若 $E|X| = E[X^+] + E[X^-] < \infty$ 則我們稱 $X$ 有 finite expectation,(因為 $E[X^+] < \infty$ 且 $E[X^-] < \infty$)。
3. 如果現在給定 $X$ 但想計算 $X$ 在某個子集上的期望值,也就是說我們現在不是取整個樣本空間 $\Omega$ 而是只有某個 $\mathcal{F}$ 中的集合 $A \in \mathcal{F}$ 則我們仍可求其期望值
\[
\int_A X dP = \int_\Omega X 1_A dP = E[X 1_A]
\]上式稱作 integral of $X$ with respect to $P$ over $A$。另外若此積分存在且有限,則我們稱 $X$ 為 integrable with respect to $P$ over $A$


注意到期望值 $E[X] = \int_\Omega X dP $ 仍不容易計算,因為樣本空間 $\Omega$ 可以內含各式各樣奇形怪狀的東西比如說 $\Omega:=\{apple, orange...\}$,此時計算 $E[X]$ 顯然會過於抽象。那麼我們想問什麼時候才可能比較容易計算 $E[X]$?答案是如果 能給出 $X$ 的 累積分佈函數 (Cumulative Distribution Function, CDF) , 記作 $F$, with respect to $P$,亦即 $F(x) := P(X \le x)$  則 $E[X]$ 可計算如下
\[
E[X] = \int_\Omega X dP = \int_{-\infty}^{\infty} x dF(x) \;\;\;\;\; (*)
\]讀者應注意到第二等式 內含的是 $x$ 不是 $X$,且 上式 $(*)$ 為 Riemann-Stigjies Integral。
若 $X$ 有機率密度函數 (Probability Density Function, f) with respect to $P$ 則回憶機率密度函數定義為 $dF(x)/dx = f(x)$ 故我們可以得到更容易計算的期望值如下
\[
E[X] = \int_{-\infty}^{\infty} x dF(x) = \int_{-\infty}^\infty x f(x) dx
\]但上式有一定限制,因為任意隨機變數必定存在 CDF 但不一定有 PDF (因為 CDF 不一定可微,故不一定有 PDF)

Comments: 
1. 若 $g$ 為 measurable 函數 on $\mathbb{R}$ 則我們亦可計算 $E[g(X)]$,亦即
\[
E[g(X)] = \int_\Omega g(X) dP = \int_{-\infty}^\infty g(x) dF(x) = \int_{-\infty}^\infty g(x) f(x)dx
\]
2. 若 $X$ 為離散隨機變數取值為 $x_1,...,x_n$ 且 具有 probability mass function, pmf with respect to $P$, 記作 $P(X = x_n) = p(x_n)$, 則
\[
E[X] = \sum_{n=1}^\infty x_n p(x_n)
\]

\[
E[g(X)] = \sum_{n=1}^\infty g(x_n) p(x_n)
\]

現在看幾個例子:
Example 1:
$\Omega :=\{1,2,3,4\}$ 且 $\mathcal{F} := \sigma ( \{1\}, \{2\}, \{3\}, \{4\})$ 且 令隨機變數 $X = i$ 其中 $i=1,2,3,4$ 且
\[
P(X=1)=1/2;\;\; P(X=2) = 1/4;\;\; P(X=3)=1/6;\;\; P(X=4)=1/12
\]現在令 $X:= 5 \cdot 1_{\{X=1\}} + 2 \cdot 1_{\{X=2\}} - 4 \cdot 1_{\{X=3 \text{ or } X=4\}}$。試求 $E[X]=?$
Solution:
\[\begin{array}{*{20}{l}}
\begin{array}{l}
E[X] = \int_\Omega  {XdP = \sum\limits_n^{} {{x_n}P\left( {X = {x_n}} \right)} } \\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = 5\int_\Omega  {{1_{\left\{ {X = 1} \right\}}}dP}  + 2\int_\Omega  {{1_{\left\{ {X = 2} \right\}}}dP}  - 4\int_\Omega  {{1_{\left\{ {X = 3orX = 4} \right\}}}dP}
\end{array}\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = 5 \cdot P\left( {X = 1} \right) + 2 \cdot P\left( {X = 2} \right) - 4 \cdot \underbrace {P\left( {X = 3orX = 4} \right)}_{ = P\left( {\left\{ {X = 3} \right\} \cup \left\{ {X = 4} \right\}} \right) = P\left( {X = 3} \right) + P\left( {X = 4} \right)}}\\
{\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = 5\left( {\frac{1}{2}} \right) + 2\left( {\frac{1}{4}} \right) - 4\left( {\frac{1}{6} + \frac{1}{{12}}} \right) = 2}
\end{array}\]

Example 2:
同上題,試求 $E[X^2]=?$
Solution:
\[\begin{array}{l}
E[{X^2}] = \int_{ - \infty }^\infty  {{x^2}f\left( x \right)dx}  = \sum\limits_n^{} {{{\left( {{x_n}} \right)}^2}P\left( {X = {x_n}} \right)} \\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = 5P\left( {X = 1} \right) + 2P\left( {X = 2} \right) + 4P\left( {\left\{ {X = 3} \right\} \cup \left\{ {X = 4} \right\}} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = 25\left( {\frac{1}{2}} \right) + 4\left( {\frac{1}{4}} \right) + 16\left( {\frac{1}{6} + \frac{1}{{12}}} \right) = \frac{{35}}{2}
\end{array}\]

Example 3:
假設 $X \sim \mathcal{N}(0,1)$ 亦即我們有 pdf
\[
f_X(x) = \frac{1}{ \sqrt{2 \pi}} \exp(-x^2/2)
\]試證 $E[X]=0$
Proof:
觀察 \[\begin{array}{l}
E[X] = \int_\Omega  {XdP}  = \int_{ - \infty }^\infty  {xdF\left( x \right)}  = \int_{ - \infty }^\infty  {xf\left( x \right)dx} \\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = \int_{ - \infty }^\infty  {x\frac{1}{{\sqrt {2\pi } }}\exp \left( { - {x^2}/2} \right)dx} \\
\begin{array}{*{20}{c}}
{}&{}&{}
\end{array} = \frac{1}{{\sqrt {2\pi } }}\int_{ - \infty }^\infty  {x{e^{\frac{{ - {x^2}}}{2}}}dx}
\end{array}\]現在令\[u = \frac{{{x^2}}}{2} \Rightarrow 2du = 2xdx\]則 $E[X] = 0$
上述積分有更快速的做法如下:因為 被積函數 $x exp(-x^2/2)$ 為奇函數,且積分範圍對原點對稱 $(-\infty, \infty)$ 故積分為零。此法同理可證 $E[X^3]=E[X^5]=E[X^{2k+1}] = 0$ $k\in \mathbb{N} $



Example 4:
令 $a>0$,設 $F$ 為在 $[0, 3a]$ 上的連續函數滿足
\[F\left( x \right): = \left\{ \begin{array}{l}
\pi ,\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array}0 \le x < a\\
4 + a - x,\begin{array}{*{20}{c}}
{}&{}
\end{array}a \le x < 2a\\
{\left( {x - 2a} \right)^2},\begin{array}{*{20}{c}}
{}&{}
\end{array}x \ge 2a
\end{array} \right.
\]且在 $(0, 3a)$ 有一階導數連續 $F' = f$ 且 $f \in L^1(0,3a)$ 試計算 下列 Lebesgue-Stieltjes integral,
\[\int_{\left( {0,3a} \right]}^{} {xdF\left( x \right)}
\]

Solution
Lebesgue-Stieltjes integral 重點在於對於不連續處需給定其測度值。給定之後其餘部分如同一般微積分課程中採用的積分方法計算即可。故我們直接求解
\[\begin{array}{l}
\int_{\left( {0,3a} \right]}^{} {xdF\left( x \right)}  = \int_{\left( {0,a} \right]}^{} {xd\left( \pi  \right)}  + \int_{\left( {a,2a} \right]}^{} {xd\left( {4 + a - x} \right)}  + \int_{\left( {2a,3a} \right]}^{} {xd{{\left( {x - 2a} \right)}^2}} \\
\begin{array}{*{20}{c}}
{}&{}&{}&{}&{}&{}&{}&{}
\end{array} + {\left. {x\left( {F\left( x \right) - F\left( {x - } \right)} \right)} \right|_{x = a}} + {\left. {x\left( {F\left( x \right) - F\left( {x - } \right)} \right)} \right|_{x = 2a}}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}&{}
\end{array} = 0 - \int_{\left( {a,2a} \right]}^{} {xdx}  + 2\int_{\left( {2a,3a} \right]}^{} {x\left( {x - 2a} \right)dx}  + a\left( {4 - \pi } \right) + 2a\left( { - 4 + a} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}&{}
\end{array} = \frac{{{a^2}}}{2} + \frac{{8{a^3}}}{3} - \left( {4 + \pi } \right)a \ \ \ \ \square
\end{array}\]


一些期望值性質我們將其記錄如下:令 令 $a,b \in \mathbb{R}$ 且 $X,Y$ 為隨機變數,回憶  $E[X \cdot 1_A] = \int_A X dP$ 則
Property (1): Absolutely Integrability
\[
\int_A XdP < \infty \Leftrightarrow \int_A |X| dP < \infty
\]
Property (2): Linearity
\[
\int_A (a X + bY)dP = a \int_A X dP + b \int_A Y dP
\]
Property (3): Countably Additivity over sets
若 $\{A_n\}$ disjoint set sequences 則
\[
\int_{\cup_n A_n} X dP = \sum_n \int_{A_n} X dP
\]
Property (4): Nonnegativity
若 $X \ge 0$ $P$-almost surely (記作 $P$-a.s.),則
\[
\int_A X dP \ge 0
\]其中 $P$-a.s. 表示 $P(X \ge 0) = 1$。

comments: 為何要在意 almost surely? 因為我們只關心測度非零的區域,測度為零的情況比如說 單點的測度為零我們並不關心。

Property (5): Monotonicity
若 $X_1 \le X \le X_2$ almost surely 則
\[
\int_A X_1 dP \le \int_A X dP \le \int_A X_2 dP
\]
Property (6): Modulus Inequality
\[
\left| \int_A X dP \right| \le \int_A |X| dP
\]


有了上述期望值的概念之後,我們可以看看如果我們取 limit 甚麼時候可以與 積分互換:也就是說 若現在 給定 一組隨機變數的數列 $\{X_n\}$ 我們想問
\[
\int_A \lim_{n \to \infty} X_n dP =?= \lim_{n \to \infty} \int_A X_n dP
\]回憶在高等微積分我們希望積分極限互換的情況的條件是需要 uniform convergence 但在機率論或者測度理論我們有以下三個極為重要的定理可以幫助我們在不需 uniform convergence 條件之下仍可達成積分與極限互換,此三個定理分別為 Dominated Convergence Theorem,Monotone Convergence Theorem、Fatou's Lemma
此三者為機率論 與積分理論的重要基石。分別紀錄如下:

考慮 $\{X_n\}$ 為任意隨機變數的 sequence。
=========================
Dominated Convergence Theorem (DCT)
若 $ P( \lim_{n \to \infty} X_n = X) =1$ 且對 $1 \leq n < \infty$,
\[
|X_n| \leq Y \text{ a.s. } \ \& \ E[Y] < \infty
\],則 $E[|X|]<\infty$ 且
\[
\lim_{n \rightarrow \infty} E[X_n] = E[\lim_{n \rightarrow \infty} X_n] = E[X]
\]========================

=========================
Monotone Convergence Theorem (MCT)
若 $ 0 \leq X_n \leq X_{n+1} (亦即為 Monotone Functions), \ \forall n \geq 1$ 且 $\lim_{n \rightarrow \infty} X_n = X$則
\[
\lim_{n \rightarrow \infty} E[X_n] = E[\lim_{n \rightarrow \infty} X_n] = E[X]
\]========================

=========================
Fatou's Lemma
若 $X_n \ge 0 \text{ a.s. } \ \forall n \geq 1$ 則
\[
E[\liminf_{n \rightarrow \infty} X_n] \leq \liminf_{n \rightarrow \infty} E[X_n]
\]========================

ref:
S. I. Resnick, A Probability Path, Birkhauser
J. M. Steele, Stochastic Calculus and Financial Applications, Springer
J. A. Gubner, Probability and Random Processes for Electrical and Computer Engineers, Cambridge.

2009年4月28日 星期二

[動態規劃] 淺談 離散時間動態規劃 (0) - Principle of Optimality & Bellman equation in Finite Horzion

這次要跟大家介紹離散時間的動態規劃 (Dynamic Programming, D.P.) 問題:
此為最佳化理論中的一種方法,由 Professor Bellman 的 Principle of Optimality 所奠基。

基本想法如下: 從最佳解往回推 (backward),逐步求解最佳值。再從初始位置沿著以求得的最佳路徑前進到最佳解。

我們現在把動態規劃問題寫下:

(有限時間) 動態規劃問題 (Dynamic Programming Problem in Finite Horizon):

考慮 Performance index (cost function)
\[
J(u) := \displaystyle \sum_{k=0}^{N-1} J(x(k), u(k)) + \Phi(x(N))
\]與離散時間的狀態方程( state equation)
\[
x(k+1) = f(x(k), u(k)), \ x(0) \ \text{is given} \\
\]其中 $x(k) \in \mathbb{R}^n$ 為系統在第 $k$ 時刻的 狀態,$u(k) \in \mathbb{R}^m$ 為控制力,函數 $f:\mathbb{R}^n \times \mathbb{R}^m \rightarrow \mathbb{R}^n$,

另外我們仍需考慮控制力拘束條件
\[
u(k) \in \Omega_k(x(k))
\]其中 $\Omega_k$ 為第 $k$ 時刻的拘束條件

Comment:
1. 關於狀態方程:上式中的 $f(x(k),u(k))$ 稱作 autonomous state equation
但事實上DP問題亦可直接考慮 Non-autonomous state equation,(亦即函數中考慮時間 $k$)
\[
f((x(k), u(k), k))
\] Example: Autonomous and Non-autonomous system
Autonomous state equation:  $x(k+1) = x(k) + u(k)$
Non-autonomous state equation: $x(k+1) = k \cdot x(k)$

2. 關於 Performance Index
\[
J(u) := \displaystyle \sum_{k=0}^{N-1} J(x(k), u(k)) + \Psi (x(N))
\] 其中 $J(x(k), u(k))$ 稱為 Branch cost。 $\Psi (x(N))$ 稱為 Terminal cost


現在我們介紹 Bellman's Principle of Optimality

==========================
Principle of Optimality
假設我們有一組 optimal sequence
\[
u^*(0), u^*(1), ..., u^*(N-1)
\]令 $x^*(0), x^*(1), ..., x^*(N)$ 為 induced states,亦即對 $k=0,1,...,N-1$
\[
u(k) \in \Omega_k(x^*(k))
\]且 $x^*(k+1) = f(x^*(k),u^*(k))$  (也就是說上述的 optimal sequence 最佳化整個問題。)
那麼現在考慮一個 子問題 (subproblem)
考慮 $k = l$ (現在考慮由 時刻 $l$ 開始 (不是從 $0$ 開始) ),初始狀態現在為 $x(l)=x^*(l)$ (亦即此時的初始狀態 $x(l)$ 落在最佳解上 $x^*(l)$),此時對應的 Cost function 為
\[
J_l(u) := \displaystyle \sum_{k=l}^{N-1} J(x(k), u(k)) + \Psi (x(N))
\]則 Principle of Optimality 告訴我們
\[
u^*(l), u^*(l+1), ..., u^*(N-1)
\]為對此 subproblem 的 optimal sequence.
==========================
Comment:
上述的 Principle of Optimality 其實表達的想法很簡單,就是如果已有一組最佳解序列在手,那麼從任意時刻開始的解也必須是最佳解。比如說我有一組最佳控制力序列
$u^*(0), u^*(1), ..., u^*(10)$ 那麼考慮從 $k=5, 6,...10$ 的子問題,則對應的控制力序列
\[
u^*(5), u^*(6),...u^*(10)
\]亦為此子問題的最佳解。

或者更淺白的說,假設今天要坐飛機到美國,最佳航程為
台灣 => 香港 =>美國
如果今天如果不從台灣開始,而是從香港開始(且香港已經在最佳路徑上),則到美國的這條路徑仍維持最佳解。(香港=>美國)


Proof: Bellman's Principle of Optimality

用矛盾證法 (Proof by contradiciton),
考慮一 子問題(subproblem) 為 從狀態 $x^*(l)$ 開始,且此狀態落在最佳解上,但 $u^*(k), \ k=l,l+1,...,N-1$ 並非 Optimal sequence。

既然上述不是 subproblem 的 Optimal sequence,則必定有其他人比此 $u^*$ 這組 sequence 可以得到更低的 cost,也就是說存在一組
\[
\hat{u}(l), \hat{u}(l+1),... \hat{u}(N-1)
\]使得 $J_l(\hat{u}) < J_l (u^*)$

現在我們定義 Optimal sequence: $\tilde{u}(k)$
\[\tilde u(k): = \left\{ \begin{array}{l}
{u^*}\left( k \right),\begin{array}{*{20}{c}}
{}
\end{array}k < l\\
\hat u\left( k \right),\begin{array}{*{20}{c}}
{}
\end{array}l \le k \le N - 1
\end{array} \right.
\] 那麼我們現在觀察 整個最佳化問題 (非子問題 ) 的 cost;亦即對 $k=0$ 開始 (而非 $k=l$),則我們有
\[\begin{array}{l}
J\left( {\tilde u} \right) = \sum\limits_{k = 0}^{N - 1} {J\left( {x\left( k \right),u\left( k \right)} \right)}  + \Psi \left( {x\left( N \right)} \right)\\
 \Rightarrow J\left( {\tilde u} \right) = \underbrace {\sum\limits_{k = 0}^{l - 1} {J\left( {x\left( k \right),u\left( k \right)} \right)} }_{{J_0}\left( {{u^*}} \right)} + \underbrace {\sum\limits_{k = l}^{N - 1} {J\left( {x\left( k \right),u\left( k \right)} \right)} }_{{J_l}\left( {\hat u} \right)} + \Psi \left( {x\left( N \right)} \right)\\
 \Rightarrow J\left( {\tilde u} \right) < {J_0}\left( {{u^*}} \right) + {J_l}\left( {{u^*}} \right) + \Psi \left( {x\left( N \right)} \right)\\
 \Rightarrow J\left( {\tilde u} \right) < J\left( {{u^*}} \right)
\end{array}
\] 上述結果說明 $\tilde{u}$ 的序列為 Optimal sequence,但此結果與 Principle of Optimality 中原本假設的 最佳序列 $u^*(0), u^*(1),..., u^*(N-1)$ 矛盾。故得證。 $\square$


有了上述的 Principle of Optimality,我們便可以建構 Dynamic Programming Equaiton (DPE) (DPE 又稱為 Bellman equation)


Bellman Equation
首先考慮我們的 Subproblem:
由初始狀態 $x(l)$,定義 Optimal cost-to-go 記作 $I(x(l), N-l)$ (亦即考慮 Subproblem 由 $x(l)$ 開始 到 $x(N)$ 的 optimal cost )
\[I(x(l),N - l): = \mathop {\min }\limits_{\scriptstyle u(k) \in {\Omega _k},\atop
\scriptstyle k \ge l} {J_l}\left( u \right)
\]由 $J_l(u)$ 定義可知,${J_l}(u): = \sum\limits_{k = l}^{N - 1} J (x(k),u(k)) + \Psi (x(N))$,故
\[
\begin{array}{l}
I(x(l),N - l): = \mathop {\min }\limits_{\begin{array}{*{20}{c}}
{u(k) \in {\Omega _k},}\\
{k \ge l}
\end{array}} \left[ {\sum\limits_{k = l}^{N - 1} J (x(k),u(k)) + \Psi (x(N))} \right]\\
 \Rightarrow I(x(l),N - l) = \mathop {\min }\limits_{\begin{array}{*{20}{c}}
{u(k) \in {\Omega _k},}\\
{k \ge l}
\end{array}} \left[ \begin{array}{l}
J(x(l),u(l))\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} + \sum\limits_{k = l + 1}^{N - 1} {J(x(k),u(k))}  + \Psi (x(N))
\end{array} \right]
\end{array}
\]由 Principle of Optimality 可知,上式可改寫為
\[
\small{\begin{array}{l}
 \Rightarrow I(x(l),N - l): = \mathop {\min }\limits_{u(k) \in {\Omega _k}} \left[ \begin{array}{l}
J(x(l),u(l))\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} + \mathop {\min }\limits_{\begin{array}{*{20}{c}}
{u(k) \in {\Omega _k},}\\
{k > l}
\end{array}} \left\{ {\sum\limits_{k = l + 1}^{N - 1} {J(x(k),u(k))}  + \Psi (x(N))} \right\}
\end{array} \right]\\
 \Rightarrow I(x(l),N - l): = \mathop {\min }\limits_{u(k) \in {\Omega _k}} \left[ {J(x(l),u(l)) + I\left( {x(l + 1),N - \left( {l + 1} \right)} \right)} \right]
\end{array}}
\]
上式稱為 Bellman Equation 或稱 Dynamic Programming Equation
\[
I(x(l), N-l) = \min_{u(l) \in \Omega_l} \{J(x(l), u(l)) + I(x(l+1), N-(l+1)) \}
\]

Example
考慮系統狀態方程表示如下:
\[x(k+1) = x(k) - u(k)
\]且 cost function 為
\[
J(u) = \displaystyle \sum_{k=0}^{N-1} (x(k+1)-u(k))^2 + x^2(k+1)
\]試求 Optimal $u^*(N-1), u^*(N-2),...$ 與其對應的 optimal cost to go

Solution
為簡便起見,這邊只求 $u(N-1)$ 與其對應的 optimal cost to go。

故首先由 $x(N-1)$ 開始,亦即 $l= N-1$,由 DPE 可知
\[\begin{array}{l}
I(x(l),N - l): = \mathop {\min }\limits_{u(k) \in {\Omega _k}} \left[ {J(x(l),u(l)) + I\left( {x(l + 1),N - \left( {l + 1} \right)} \right)} \right]\\
 \Rightarrow I(x(N - 1),1) = \mathop {\min }\limits_{u(k) \in {\Omega _k}} \left[ {J(x(N - 1),u(N - 1)) + I\left( {x(N),0} \right)} \right]\\
 \Rightarrow I(x(N - 1),1) = \mathop {\min }\limits_{u(k) \in {\Omega _k}} \left[ {{{\left[ {x(N) - u(N - 1)} \right]}^2} + {x^2}\left( N \right) + {\rm{0}}} \right]\\
 \Rightarrow I(x(N - 1),1) = \mathop {\min }\limits_{u(k) \in {\Omega _k}} \left[ {{{\left[ {x(N) - u(N - 1)} \right]}^2} + {x^2}\left( N \right)} \right]
\end{array}
\]現在再由系統狀態方程 $x\left( N \right) = x(N - 1) - u(N - 1) $ 代入上式可得
\[I(x(N - 1),1) = \mathop {\min }\limits_{u(N - 1)} \left[ {{{\left[ {x(N - 1) - 2u(N - 1)} \right]}^2} + {{\left[ {x(N - 1) - u(N - 1)} \right]}^2}} \right]
\] 由 FONC 求最佳解:
\[\begin{array}{l}
\frac{{\partial I(x(N - 1),1)}}{{\partial u(N - 1)}} = 0 \\
 \Rightarrow {u^*}(N - 1) = \frac{3}{5}x(N - 1)
\end{array}
\]其對應的Optimal cost-to-go
\[\begin{array}{l}
 \Rightarrow {\left. {I(x(N - 1),1)} \right|_{{u^*}}} = {\left[ {x(N - 1) - 2\frac{3}{5}x(N - 1)} \right]^2}\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}&{}&{}&{}&{}
\end{array} + {\left[ {x(N - 1) - \frac{3}{5}x(N - 1)} \right]^2}\\
 \Rightarrow {\left. {I(x(N - 1),1)} \right|_{{u^*}}} = \underbrace {\frac{1}{5}}_{{K_{N - 1}}}{x^2}(N - 1)
\end{array}
\] 上式子中的 $K_{N-1}$ 被稱作 feedback-gain at $N-1$ stage。

接著我們計算 $x(N-2)$ 亦即 $l=N-2$ 的 Optimal cost-to-go:
\[\begin{array}{*{20}{l}}
{I(x(l),N - l): = \mathop {\min }\limits_{u(k) \in {\Omega _k}} \left[ {J(x(l),u(l)) + I\left( {x(l + 1),N - \left( {l + 1} \right)} \right)} \right]}\\
{ \Rightarrow I(x(N - 2),2) = \mathop {\min }\limits_{u(N - 2) \in {\Omega _{N - 2}}} \left[ {J(x(N - 2),u(N - 2)) + I\left( {x(N - 1),1} \right)} \right]}\\
{ \Rightarrow I(x(N - 2),2) = \mathop {\min }\limits_{u(N - 2)} \left[ \begin{array}{l}
{\left( {x\left( {N - 1} \right) - u\left( {N - 2} \right)} \right)^2}\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} + {x^2}\left( {N - 1} \right) + \frac{1}{5}{x^2}\left( {N - 1} \right)
\end{array} \right]}\\
{ \Rightarrow I(x(N - 2),2) = \mathop {\min }\limits_{u(N - 2)} \left[ \begin{array}{l}
{\left( {x\left( {N - 2} \right) - 2u\left( {N - 2} \right)} \right)^2}\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} + {\left[ {x\left( {N - 2} \right) - u\left( {N - 2} \right)} \right]^2}\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} + \frac{1}{5}{\left[ {x\left( {N - 2} \right) - u\left( {N - 2} \right)} \right]^2}
\end{array} \right]}
\end{array}
\] 由FONC: $\frac{\partial }{{\partial u(N - 2)}} = 0$,可知
\[\begin{array}{l}
\left\{ \begin{array}{l}
 - 4\left( {x\left( {N - 2} \right) - 2u\left( {N - 2} \right)} \right)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} - 2\left[ {x\left( {N - 2} \right) - u\left( {N - 2} \right)} \right] - \frac{2}{5}\left[ {x\left( {N - 2} \right) - u\left( {N - 2} \right)} \right]
\end{array} \right\} = 0\\
 \Rightarrow {u^*}\left( {N - 2} \right) = \frac{8}{{13}}x\left( {N - 2} \right)
\end{array}
\] 將上式代回 $I(x(N-2),2)$ 可得對應的 Optimal Cost to go 為
\[\begin{array}{l}
I(x(N - 2),2) = \left[ {\frac{9}{{{{13}^2}}}{x^2}\left( {N - 2} \right) + \frac{{25}}{{{{13}^2}}}{x^2}\left( {N - 2} \right) + \frac{5}{{{{13}^2}}}{x^2}\left( {N - 2} \right)} \right]\\
 \Rightarrow I(x(N - 2),2) = \frac{{39}}{{{{13}^2}}}{x^2}\left( {N - 2} \right)
\end{array}
\] $\square$

延伸閱讀:
[動態規劃] 淺談 離散時間動態規劃 (1) - Bellman equation in Infinite Horzion

2009年4月16日 星期四

[最佳化] 淺談 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
\]

2009年4月15日 星期三

[最佳化] 淺談 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