在計算數學(Computational Mathematics)，數值分析(Numerical Analysis)領域裡，介紹在$\mathbb{R}^n$ 上解非線性方程組 (System of Nonlinear Equations) 學術上有名的牛頓法(Newton Method) ，它也是很多最佳化演算法的核心概念
現實問題中，我們通常會想要找一組解 $x \in \mathbb{R}^{n}$上，使得 $f(x) = \vec{0}$，
其中 $f : \mathbb{R}^n \longrightarrow \mathbb{R}^{m}$，
詳細一點可以寫成一個有 $n$ 個未知數(unknowns) 跟 $m$ 個限制式(constraints)
$\text{Find } x = (x_1,x_2,...,x_n) \in \mathbb{R}^n \quad \text{such that } \left\{\begin{array}{c}f_1(x_1,x_2,x_3,...,x_n)=0\\f_2(x_1,x_2,x_3,...,x_n)=0\\f_3(x_1,x_2,x_3,...,x_n)=0\\f_4(x_1,x_2,x_3,...,x_n)=0\\...\\f_m(x_1,x_2,x_3,...,x_n)=0\\\end{array}\right.$
精確解我們記做 $x^{*}$
1.如果 $f$ 是線性函數，即為我們熟知的線性聯立方程組 $Ax = b $ 或寫成
$\left(\forall i = 1,2,3,...m \quad f_{i}(x) = \sum^{n}_{j=1}a_{ij}x_j = b_i \right)$
線性系統可以用大家熟知高斯消去法(Gauss Elimination)，經過有限步驟得到精確解。
2.如果$f$ 是非線性但連續且可微，即 Jacobian Matrix $f$ 存在，則我們可以用牛頓迭代法，
[以下是虛擬碼 (pseudo code)]
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
$\text{Given Initial Guess : } x^0 \in \mathbb{R}^n , \quad \text{sufficient small stop criterion : } \epsilon > 0 $
$\text{For }( k=0,1,2,3,... )$
$\qquad$ $\qquad$ $J_f (x^k) d^k = f(x^k) \quad \text{[Solve Linear System to get }d^k \in \mathbb{R}^n] \qquad $
其中 $f : \mathbb{R}^n \longrightarrow \mathbb{R}^{m}$，
詳細一點可以寫成一個有 $n$ 個未知數(unknowns) 跟 $m$ 個限制式(constraints)
$\text{Find } x = (x_1,x_2,...,x_n) \in \mathbb{R}^n \quad \text{such that } \left\{\begin{array}{c}f_1(x_1,x_2,x_3,...,x_n)=0\\f_2(x_1,x_2,x_3,...,x_n)=0\\f_3(x_1,x_2,x_3,...,x_n)=0\\f_4(x_1,x_2,x_3,...,x_n)=0\\...\\f_m(x_1,x_2,x_3,...,x_n)=0\\\end{array}\right.$
精確解我們記做 $x^{*}$
1.如果 $f$ 是線性函數，即為我們熟知的線性聯立方程組 $Ax = b $ 或寫成
$\left(\forall i = 1,2,3,...m \quad f_{i}(x) = \sum^{n}_{j=1}a_{ij}x_j = b_i \right)$
線性系統可以用大家熟知高斯消去法(Gauss Elimination)，經過有限步驟得到精確解。
2.如果$f$ 是非線性但連續且可微，即 Jacobian Matrix $f$ 存在，則我們可以用牛頓迭代法，
[以下是虛擬碼 (pseudo code)]
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
$\text{Given Initial Guess : } x^0 \in \mathbb{R}^n , \quad \text{sufficient small stop criterion : } \epsilon > 0 $
$\text{For }( k=0,1,2,3,... )$
$\qquad$ $\qquad$ $J_f (x^k) d^k = f(x^k) \quad \text{[Solve Linear System to get }d^k \in \mathbb{R}^n] \qquad $
$\qquad$ $\qquad$ $x^{k+1} := x^k -d^k \quad \text{[Newton Iteration]} $
$\qquad$ $\qquad$ $\text{If }(|| f(x^{k+1}) || < \epsilon ) \quad \text{break and output } x^{k+1} \text{ as solution}$
$\text{EndFor }$
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
其中 $x^k$ 代表第 $k$ 次計算後的結果 !!
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
其中 $x^k$ 代表第 $k$ 次計算後的結果 !!
其中 $|| u || = \sqrt{u_1^2+u_2^2+....+u_n^2}$ 為歐式距離 (Euclidean distance)
其中 $J_f$ 為 $m \times n$ Jacobian Matrix (多維空間的微分 $f'$)
其中 $J_f$ 為 $m \times n$ Jacobian Matrix (多維空間的微分 $f'$)
$$J_f =\left[ \begin{array}{cccc}\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & ... & \frac{\partial f_1}{\partial x_n} \ \\ \frac{\partial f_2}{\partial x_1}& \frac{\partial f_2 }{\partial x_2} & ... & \frac{\partial f_2}{\partial x_n} \\ ... & ... & ... & ... \\ \frac{\partial f_m}{\partial x_1}& \frac{\partial f_m}{\partial x_2} & ... & \frac{\partial f_m }{\partial x_n} \\ \end{array}\right] \qquad \frac{\partial f_i}{\partial x_j}(x) = \frac{\partial f_i}{\partial x_j}(x_1,x_2,...x_n) $$
其中 $ d^k := \underbrace{J_f^{-1} (x^k)\cdot f(x^k)}_{\text{Matrix-Vector Multiplication}}$
其中 $ d^k := \underbrace{J_f^{-1} (x^k)\cdot f(x^k)}_{\text{Matrix-Vector Multiplication}}$
但實務上我們是直接解線性聯立方程組。而不會計算 $J_f^{-1}$ ，因為耗時昂貴，只是寫法上好理解。
[一維的例子]
$m=1,n = 1$ 則 Newton Iteration 可以寫成 $ x^{k+1} := x^k - (f')^{-1}(x^k) f(x^k) = x^{k} - \frac{f(x^{k})}{f'(x^{k})}$ ,
[一維的例子]
$m=1,n = 1$ 則 Newton Iteration 可以寫成 $ x^{k+1} := x^k - (f')^{-1}(x^k) f(x^k) = x^{k} - \frac{f(x^{k})}{f'(x^{k})}$ ,
注意: 這邊的 $(f')^{-1}$ 指的是"倒數"在分母的意思，不是"反函數"(inverse)的意思
[概念澄清 I]
牛頓法是一種迭代法(Iterative Method) 即我們是製造數列 $\{x^k\}$ 使得 $$\lim_{k\rightarrow \infty} x^k = x^{*} \quad \text{ s.t } f(x^{*}) = \vec{0}$$ ，我們無法得到精確解，但是我們可以得到無窮靠近的近似解!!有別於高斯消去法的有限步驟得到精確解，直接法(Direct Method)。
[概念澄清II]
由於實際跑演算法，我們無法跑無窮次，而當 $x^{k} \rightarrow x^{*} \overset{f \text{ is continuous}}{\Longrightarrow} f(x^{k}) \longrightarrow \underset{=\vec{0}}{f(x^{*})} \overset{\forall i , f_i(x) \rightarrow 0 }{\Longrightarrow} ||f(x^{k})|| \rightarrow 0 $
所以停止條件設為 $|| f(x^{k+1}) || < \epsilon $ ，當你給定的 $\epsilon $ 越小(解的品質要求越精準)，則所花的計算時間就會拉長，是一種 TradeOff ，$\epsilon$ 需要依照實際需求而定 !!
[概念澄清III]
起始點的選擇 $x_0 \in \mathbb{R}^n$ 是使用牛頓法的大議題，有時候會"收斂"，有時候會"發散"，即只有 "Local Convergence"，必須分析"收斂域"(哪些 Input $x_0$ 會收斂)，但經驗是"收斂域"往往很小!!
[概念澄清 IV]
不同的起始點可能會收斂到不同的解，假如 $x^{*}\in \{x\in \mathbb{R}^n : f(x)=\vec{0} \} $ 有多組解的話
[優點]
理論證明收斂性 Quardratic Convergence ， 即當$ k $ 足夠大，存在Constant $C$ s.t $||x^{k+1}- x^{*}|| \leq C ||x^{k}- x^{*}||^2$ 。以白話來說，假如 $C=1$，這次得到的 $x_k$ 與 $x^{*}$ 距離差距上界為 $10^{-4}$ 則下一次迭代比較$x^{k+1}$與$x^{*}$的距離，可以把差距上界縮小成 $10^{-8}$左右
[進階最佳化應用]
在數值最佳化(Numerical Optimization)領域，即$f:\mathbb{R}^{n}\longrightarrow \mathbb{R}$，
求解無限制式最佳化問題(unconstrained optimization problem) $$\underset{x \in \mathbb{R}^n}{\text{min}} f(x)$$
我們也可以使用牛頓法求解 gradient $f$ 方程組，$\bigtriangledown f(x) = \vec{0}$，即白話文來說"極值發生在微分=0"，詳細一點可以寫成以下方程組 $m = n$ ，p.s $\bigtriangledown f : \mathbb{R}^{n} \longrightarrow \mathbb{R}^{n}$
$$\text{Find } x \in \mathbb{R}^n \text{ such that }\left\{\begin{array}{c}\frac{\partial f}{\partial x_1}(x)=0\\\frac{\partial f}{\partial x_2}(x)=0\\\frac{\partial f}{\partial x_3}(x)=0\\\frac{\partial f}{\partial x_4}(x)=0\\...\\\frac{\partial f}{\partial x_n}(x)=0\\\end{array}\right.$$
同理我們可以計算 Jacobian Matrix of $\bigtriangledown f$ ，剛好會變成 $n \times n$ Hessian Matrix of $f$ (Hessian 可以看成廣義的二次微分)
即 $J_{\bigtriangledown f} = H_{f} = \left[ \begin{array}{cccc}\frac{\partial^2 f}{\partial x^2_1} & \frac{\partial^2 f}{\partial x_1\partial x_2} & ... & \frac{\partial^2 f}{\partial x_1\partial x_n} \ \\ \frac{\partial^2 f}{\partial x_2\partial x_1}& \frac{\partial^2 f }{ \partial x_2^2} & ... & \frac{\partial^2 f}{\partial x_2\partial x_n} \\ ... & ... & ... & ... \\ \frac{\partial^2 f}{\partial x_n \partial x_1}& \frac{\partial^2 f}{\partial x_n \partial x_2} & ... & \frac{\partial^2 f}{\partial x_n^2 } \\ \end{array}\right] $
注意 : $\frac{\partial^2 f}{\partial x_j^2}$ 代表 $f$ 對 $x_j$ 做二次偏微
注意 : 在大多情況下，Hessian 每一項元素(entity)都是連續函數則 Hessian 是對稱矩陣，即 $\frac{\partial^2 f}{\partial x_i\partial x_j} = \frac{\partial^2 f}{\partial x_j\partial x_i} \quad \forall i , j $ ，代表說我們只需要計算上半三角形 $\frac{n(n+1)}{2}$ 項就足夠了 !!
所以 Newton Iteration 可以改寫成
$x^{k+1}:= x^{k} - d^{k}$
$ d^{k} := H^{-1}_f(x^{k})\cdot \bigtriangledown f(x^{k})$
根據微積分知識我們知道以下事情
$\underbrace{\{x \in \mathbb{R}^{n} : \bigtriangledown f(x) = \vec{0} , H_f(x) \succ 0 \}}_{\text{second derivative test}}\subseteq \{x \in \mathbb{R}^{n}: x \text{ is local minima } \} \subseteq \underbrace{\{x \in \mathbb{R}^{n} : \bigtriangledown f(x) = \vec{0} \}}_{\text{first derivative test}}$
根據 Second Derivative Test 在求解最佳化問題使用牛頓法我們還需要多檢查收斂的解是否滿足 $H_f(x^{k}) \succ 0$ ，正定矩陣 (positive definite) ，才能確定是否為 local minima !! 如果 $H_f(x^{k}) \prec 0$，則為 local maxima， 其他狀況代表收斂到任何可能 local minima , local maxima , saddle point
[概念澄清 I]
牛頓法是一種迭代法(Iterative Method) 即我們是製造數列 $\{x^k\}$ 使得 $$\lim_{k\rightarrow \infty} x^k = x^{*} \quad \text{ s.t } f(x^{*}) = \vec{0}$$ ，我們無法得到精確解，但是我們可以得到無窮靠近的近似解!!有別於高斯消去法的有限步驟得到精確解，直接法(Direct Method)。
[概念澄清II]
由於實際跑演算法，我們無法跑無窮次，而當 $x^{k} \rightarrow x^{*} \overset{f \text{ is continuous}}{\Longrightarrow} f(x^{k}) \longrightarrow \underset{=\vec{0}}{f(x^{*})} \overset{\forall i , f_i(x) \rightarrow 0 }{\Longrightarrow} ||f(x^{k})|| \rightarrow 0 $
所以停止條件設為 $|| f(x^{k+1}) || < \epsilon $ ，當你給定的 $\epsilon $ 越小(解的品質要求越精準)，則所花的計算時間就會拉長，是一種 TradeOff ，$\epsilon$ 需要依照實際需求而定 !!
[概念澄清III]
起始點的選擇 $x_0 \in \mathbb{R}^n$ 是使用牛頓法的大議題，有時候會"收斂"，有時候會"發散"，即只有 "Local Convergence"，必須分析"收斂域"(哪些 Input $x_0$ 會收斂)，但經驗是"收斂域"往往很小!!
[概念澄清 IV]
不同的起始點可能會收斂到不同的解，假如 $x^{*}\in \{x\in \mathbb{R}^n : f(x)=\vec{0} \} $ 有多組解的話
[優點]
理論證明收斂性 Quardratic Convergence ， 即當$ k $ 足夠大，存在Constant $C$ s.t $||x^{k+1}- x^{*}|| \leq C ||x^{k}- x^{*}||^2$ 。以白話來說，假如 $C=1$，這次得到的 $x_k$ 與 $x^{*}$ 距離差距上界為 $10^{-4}$ 則下一次迭代比較$x^{k+1}$與$x^{*}$的距離，可以把差距上界縮小成 $10^{-8}$左右
[進階最佳化應用]
在數值最佳化(Numerical Optimization)領域，即$f:\mathbb{R}^{n}\longrightarrow \mathbb{R}$，
求解無限制式最佳化問題(unconstrained optimization problem) $$\underset{x \in \mathbb{R}^n}{\text{min}} f(x)$$
我們也可以使用牛頓法求解 gradient $f$ 方程組，$\bigtriangledown f(x) = \vec{0}$，即白話文來說"極值發生在微分=0"，詳細一點可以寫成以下方程組 $m = n$ ，p.s $\bigtriangledown f : \mathbb{R}^{n} \longrightarrow \mathbb{R}^{n}$
$$\text{Find } x \in \mathbb{R}^n \text{ such that }\left\{\begin{array}{c}\frac{\partial f}{\partial x_1}(x)=0\\\frac{\partial f}{\partial x_2}(x)=0\\\frac{\partial f}{\partial x_3}(x)=0\\\frac{\partial f}{\partial x_4}(x)=0\\...\\\frac{\partial f}{\partial x_n}(x)=0\\\end{array}\right.$$
同理我們可以計算 Jacobian Matrix of $\bigtriangledown f$ ，剛好會變成 $n \times n$ Hessian Matrix of $f$ (Hessian 可以看成廣義的二次微分)
即 $J_{\bigtriangledown f} = H_{f} = \left[ \begin{array}{cccc}\frac{\partial^2 f}{\partial x^2_1} & \frac{\partial^2 f}{\partial x_1\partial x_2} & ... & \frac{\partial^2 f}{\partial x_1\partial x_n} \ \\ \frac{\partial^2 f}{\partial x_2\partial x_1}& \frac{\partial^2 f }{ \partial x_2^2} & ... & \frac{\partial^2 f}{\partial x_2\partial x_n} \\ ... & ... & ... & ... \\ \frac{\partial^2 f}{\partial x_n \partial x_1}& \frac{\partial^2 f}{\partial x_n \partial x_2} & ... & \frac{\partial^2 f}{\partial x_n^2 } \\ \end{array}\right] $
注意 : $\frac{\partial^2 f}{\partial x_j^2}$ 代表 $f$ 對 $x_j$ 做二次偏微
注意 : 在大多情況下，Hessian 每一項元素(entity)都是連續函數則 Hessian 是對稱矩陣，即 $\frac{\partial^2 f}{\partial x_i\partial x_j} = \frac{\partial^2 f}{\partial x_j\partial x_i} \quad \forall i , j $ ，代表說我們只需要計算上半三角形 $\frac{n(n+1)}{2}$ 項就足夠了 !!
所以 Newton Iteration 可以改寫成
$x^{k+1}:= x^{k} - d^{k}$
$ d^{k} := H^{-1}_f(x^{k})\cdot \bigtriangledown f(x^{k})$
根據微積分知識我們知道以下事情
$\underbrace{\{x \in \mathbb{R}^{n} : \bigtriangledown f(x) = \vec{0} , H_f(x) \succ 0 \}}_{\text{second derivative test}}\subseteq \{x \in \mathbb{R}^{n}: x \text{ is local minima } \} \subseteq \underbrace{\{x \in \mathbb{R}^{n} : \bigtriangledown f(x) = \vec{0} \}}_{\text{first derivative test}}$
根據 Second Derivative Test 在求解最佳化問題使用牛頓法我們還需要多檢查收斂的解是否滿足 $H_f(x^{k}) \succ 0$ ，正定矩陣 (positive definite) ，才能確定是否為 local minima !! 如果 $H_f(x^{k}) \prec 0$，則為 local maxima， 其他狀況代表收斂到任何可能 local minima , local maxima , saddle point
[Implementation 加速相關]
牛頓法是求解非線性方程組，甚至最佳化演算法主要工具，因為它的收斂性很好，但是在計算 Jacobian/Hessian 上非常的昂貴耗時，所以計算數學家後來發展 Quasi Newton Method，例如 BFGS 算法，原理為加入矩陣的 "Initial Guess" 在迭代過程中會漸漸近似到 Jacobian/Hessian 矩陣 ，而在牛頓法每次迭代子程序都會解一個線性聯立方程組 $J_f (x^k) d^k = f(x^k)$，即 $Ax=b$，在數值線性代數(Numerical Linear Algebra)領域，有許多其他演算法(有別於高斯消去法)可以更有效率解決大型線性系統，如: LU , CG(Conjugate Gradient) , GMRES 等等。
[特殊最佳化與統計應用]
在最小平方法(least square)與非線性迴歸(nonlinear regression)，也有 Newton 想法，如 Gauss-Newton , Levenberg-Marquardt 的演算法
[缺點與替代方案]
1. 在解非線性方程組裡 $f $ is vector-valued，牛頓法需要 Jacobian 存在，如果不存在無法使用。此時廣義的割線法(Secant Method) Broyden method 可以使用
2. 在解最佳化問題時 $f$ is scalar-valued ，牛頓法需要 gradient 存在，而且 Hessian 也要存在。如果 Hessian 不存在，使用大家熟知的 gradient descent !! ，如果連 gradient 都不存在，使用 subgradient !!
牛頓法是求解非線性方程組，甚至最佳化演算法主要工具，因為它的收斂性很好，但是在計算 Jacobian/Hessian 上非常的昂貴耗時，所以計算數學家後來發展 Quasi Newton Method，例如 BFGS 算法，原理為加入矩陣的 "Initial Guess" 在迭代過程中會漸漸近似到 Jacobian/Hessian 矩陣 ，而在牛頓法每次迭代子程序都會解一個線性聯立方程組 $J_f (x^k) d^k = f(x^k)$，即 $Ax=b$，在數值線性代數(Numerical Linear Algebra)領域，有許多其他演算法(有別於高斯消去法)可以更有效率解決大型線性系統，如: LU , CG(Conjugate Gradient) , GMRES 等等。
[特殊最佳化與統計應用]
在最小平方法(least square)與非線性迴歸(nonlinear regression)，也有 Newton 想法，如 Gauss-Newton , Levenberg-Marquardt 的演算法
[缺點與替代方案]
1. 在解非線性方程組裡 $f $ is vector-valued，牛頓法需要 Jacobian 存在，如果不存在無法使用。此時廣義的割線法(Secant Method) Broyden method 可以使用
2. 在解最佳化問題時 $f$ is scalar-valued ，牛頓法需要 gradient 存在，而且 Hessian 也要存在。如果 Hessian 不存在，使用大家熟知的 gradient descent !! ，如果連 gradient 都不存在，使用 subgradient !!
[以上純為學術經驗交流知識分享，如有錯誤或建議可留言~~]
by Plus & Minus 2017.08
留言
張貼留言