在計算數學(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
留言
張貼留言