現在高數也要介紹牛頓法了，一般都是從幾何上直接以“切線法”直接引入的。這種引入方式的確很簡單，但如果脫離深入一點的分析就不大容易解釋后面的事情。所以就在想怎么更直接地從分析的角度來一條線貫穿，把整個過程說一說。
文章目錄
- 一、牛頓迭代解非線性方程
- 二、牛頓迭代解極值問題
直接開始：
一、牛頓迭代解非線性方程
考慮一個非線性方程：
f ( x ) = 0 (*) f(x)=0 \tag{*} f(x)=0(*)
這個方程沒有公式可以直接解，于是考慮迭代法。其思想就是構造一個迭代公式：
x k + 1 = φ ( x k ) (1) x_{k+1}=\varphi(x_k) \tag{1} xk+1?=φ(xk?)(1)
如果 k → ∞ 時 x k + 1 → x ? k\to \infty 時 x_{k+1}\to x^* k→∞時xk+1?→x?，那么就可以不斷地利用這個公式以 x k + 1 x_{k+1} xk+1?去逼近方程的根 x ? x^* x?。
一般來說更常見的迭代格式是：
x k + 1 = x k ? Δ ( x k ) (2) x_{k+1}=x_k-\Delta(x_k) \tag{2} xk+1?=xk??Δ(xk?)(2)
比如梯度下降就是這種格式。
那么問題就很自然了：如何構造這個迭代公式 Δ ( x k ) \Delta(x_k) Δ(xk?)？
先從最簡單的方向考慮：如果 x k x_k xk? 已經十分接近 x ? x^* x? 時，那么這個迭代的長度 Δ ( x k ) \Delta(x_k) Δ(xk?)應該是多少？
這時候問題就簡單得多了，因為當 x k x_k xk? 已經十分接近 x ? x^* x?時我們可以認為再移動 Δ ( x k ) \Delta(x_k) Δ(xk?) 就可以使得
f ( x k + 1 ) = f ( x k ? Δ ( x k ) ) = 0 (3) f(x_{k+1})=f(x_k-\Delta(x_k))=0 \tag{3} f(xk+1?)=f(xk??Δ(xk?))=0(3)
這下問題就可以轉化為：如何通過（3）式和（2）式得到 Δ ( x k ) \Delta(x_k) Δ(xk?)？
那么很自然地，如果能通過（2）（3）式導出一個關于 Δ ( x k ) \Delta(x_k) Δ(xk?)的迭代式就行了。從（3）式很容易想到如果把第二個式子展開，就可以得到一個 \Delta(x_k) 的解析式，從而就能得到它的顯式表達。這樣很容易想到 Taylor 展開（一階逼近）：
f ( x k ? Δ ( x k ) ) ≈ f ( x k ) + f ′ ( x k ) ? ( ? Δ ( x k ) ) (4) f\left( x_{k}-\Delta \left( x_{k}\right) \right) \approx f\left( x_{k}\right) +f'\left( x_{k}\right) \cdot \left( -\Delta \left( x_{k}\right) \right) \tag{4} f(xk??Δ(xk?))≈f(xk?)+f′(xk?)?(?Δ(xk?))(4)
把（4）式代回（3）式就可以得到牛頓迭代的完整表達式：
Δ ( x k ) = f ( x k ) f ′ ( x k ) \Delta \left( x_{k}\right) =\dfrac{f\left( x_{k}\right) }{f'\left( x_{k}\right) } \\ Δ(xk?)=f′(xk?)f(xk?)?
于是得到牛頓迭代的完整格式：
當然這里是基于如果 x k x_k xk?已經十分接近 x ? x^* x? 時的情況推導，直接利用上面的分析還不太容易解釋從何意一個可行的點（根的附近）出發，（5）式一定能讓 x k x_k xk? 不斷地接近 x ? x^* x?。
那么這時結合幾何來簡單看看就好說了：
結合上圖的標記就很容易知道，無論如何這個迭代公式總會使 x k x_k xk? 不斷接近 x ? x^* x?。
這里再記幾筆：
1、如果不是上述的單調情況，那么牛頓迭代是不會收斂到圖中所示的 x^* 的。這從局部收斂性的描述很容易理解。
2、嚴格的證明可以借助另一個定理。在理論篇再來討論。
二、牛頓迭代解極值問題
解方程是牛頓法的第一種情況，而這種情況其實可以很容易地推廣到極值問題。因為極值的必要條件是：
f ′ ( x ) = 0 (**) f'\left( x\right) =0 \tag{**} f′(x)=0(**)
這里最簡單的一種思路就是直接令
g ( x ) = f ′ ( x ) g(x)=f'\left( x\right) \\ g(x)=f′(x)
那么自然就有它的迭代公式：
這種理解十分便于記憶，所以還是十分推薦的。
但另外一種思路更適合于推廣，再來介紹一下。
基本的思考方式和上面解非線性方程一樣，我們仍然考慮 x k x_k xk? 在極值點 x ? ? x^{**} x?? 附近，那么問題就變成了：
給多長的增量 Δ ( x k ) \Delta \left( x_{k}\right) Δ(xk?)，能使得點 x k + 1 = x k ? Δ ( x k ) x_{k+1}=x_{k}-\Delta \left( x_{k}\right) xk+1?=xk??Δ(xk?) 為極值點（或者近似的極值點），換言之希望此時 f ( x k + 1 ) ≈ 0 f(x_{k+1})\approx 0 f(xk+1?)≈0。
類似地，我們仍然考慮它的 T a y l o r Taylor Taylor 展開，不過此時展開到二階：
f ( x k + 1 ) ≈ f ( x k ) + f ′ ( x k ) ? ( ? Δ ( x k ) ) + 1 2 f ′ ′ ( x k ) ? ( ? Δ ( x k ) ) 2 f\left( x_{k+1}\right) \approx f\left( x_{k}\right) +f'\left( x_{k}\right) \cdot \left( -\Delta \left( x_{k}\right)\right) +\frac{1}{2}f''\left( x_{k}\right) \cdot \left( -\Delta \left( x_{k}\right) \right) ^{2} \\ f(xk+1?)≈f(xk?)+f′(xk?)?(?Δ(xk?))+21?f′′(xk?)?(?Δ(xk?))2
為了方便描述我們記 x k + 1 = x k + z x_{k+1}=x_{k}+z xk+1?=xk?+z，這樣二階展開就可以寫成：
f ( x k + 1 ) ≈ f ( x k ) + f ′ ( x k ) ? z + 1 2 f ′ ′ ( x k ) ? z 2 = g ( z ) (7) f\left( x_{k+1}\right) \approx f\left( x_{k}\right) +f'\left( x_{k}\right) \cdot z +\frac{1}{2}f''\left( x_{k}\right) \cdot z^2=g(z) \tag{7} f(xk+1?)≈f(xk?)+f′(xk?)?z+21?f′′(xk?)?z2=g(z)(7)
繼續考慮 g ( z ) g(z) g(z) 的極值：
d g ( z ) d z = f ′ ( x x ) + f ′ ′ ( x k ) ? z = 0 \dfrac{dg\left( z\right) }{dz}=f'\left( x_{x}\right) + f''\left( x_{k}\right) \cdot z=0 \\ dzdg(z)?=f′(xx?)+f′′(xk?)?z=0
可以解出：
\bbox[pink]{z=-\dfrac{f’\left( x_{x}\right) }{f’'\left( x_{k}\right) }} \
這和（6）式中的迭代長度是一樣的。
注 1：這里有個很有意思的事情，如果我們只考慮一階逼近
f\left( x_{k+1}\right) \approx f\left( x_{k}\right) +f’\left( x_{k}\right)\cdot z \
那么讓它對 z 求導來求極值就直接得到： f’\left( x_{k}\right)=0
很顯然這種做法就不太可取了。由此也可以看出一階逼近的粗糙，以至于在討論極值問題時就不適合直接使用了。
注 2： 有同學可能也會問，為什么不直接考慮該函數任意值在 x_k 附近的展開？
也就是說考慮：
f\left( x\right) \approx f\left( x_{k}\right) +f’\left( x_{k}\right) \left( x-x_{k}\right) +\dfrac{1}{2}f’'\left( x_{k}\right) \left( x-x_{k}\right) ^{2} \
此時直接考慮 f’(x)=0，那么容易得到：
f’\left( x_{k}\right) +f’'\left( x_{k}\right) \left( x-x_{k}\right) \approx 0 \
于是：
x \bbox[red]{\approx} x_{k}-\dfrac{f’\left( x_{k}\right) }{f’'\left( x_{k}\right) } \tag{8}
這就有意思了，前面的牛頓迭代式這里都是嚴格的等號。而上式卻是不等號。這里就可以再借題發揮一下：
上面的討論中（7）式是直接令逼近式嚴格等于g(z)，因此得到的迭代步長就是嚴格等于。
而（8）式是直接令 f’(x)=0，那么將二階近似代回計算時自然就要變成約等號。那么得到的值就自然是約等于。不過借助（8）式也方便說明一些問題。比如：此時得到的 x 實際上就是極值點，那么如果 x_k 本身就離極值點很近，這個式子就更加精確；相反，如果 x_k 離極值本身還很遠，那么這個式子就十分粗糙。
注 3： 上面的（7）式也是非常有意思的，在極值問題中這個式子的幾何意義就十分明確了。先簡單移項：
f\left( x_{k+1}\right)- f\left( x_{k}\right) \approx f’\left( x_{k}\right) \cdot z +\frac{1}{2}f’'\left( x_{k}\right) \cdot z^2 \
這時右端的導數也是 \dfrac{dg\left( z\right) }{dz}，那么此時令 \dfrac{dg\left( z\right) }{dz}=0 從幾何上講意思就是：取改 z 為多少能使得下一步迭代時 f\left( x_{k+1}\right) 相對上一步迭代時 f\left( x_{k}\right) 的改變最大？
這里就能很明顯體會到，用上述思路其實就是在整體尋優的過程中構造了一個局部尋優的問題。
接下來討論多元函數的零點問題
三、多元函數的零點問題
多元的情形自然復雜，要完全寫對需要非常仔細。但過于仔細又容易影響總體思路，所以先不管式子是否正確，直接對照一元情形憑感覺把推導寫完再說：
另外我們這里考慮的是普通微積分里討論的多元函數，也就是說 f:R^n\to R。
（1）標量函數
回顧一元的情形，我們仍然是考慮在零點附近的一個 \mathbf{x}k ，給它一個增量 \mathbf{z}，使得 \mathbf{x}{k+1}=\mathbf{x}_k+\mathbf{z} 接近零點 \mathbf{x}^{**}。類似一元的情形，考慮函數的一階逼近：
\bbox[yellow]{f(\mathbf{x}_k+\mathbf{z}) \approx f(\mathbf{x}_k)+\nabla f(\mathbf{x}_k)\cdot\mathbf{z}} \
同樣地，考慮逼近式等于 0，于是有：
f(\mathbf{x}_k)+\nabla f(\mathbf{x}_k)\cdot\mathbf{z} = 0 \
到這里就要注意一下了:
f(\mathbf{x}_k) 是標量
\nabla f(\mathbf{x}_k)\cdot\mathbf{z} 也必須是標量
但 \mathbf{z} 是向量，所以 \nabla f(\mathbf{x}_k) 也是向量
那么我們先把寫法再規范一下：
f(\mathbf{x}_k)+ \left[ \nabla f(\mathbf{x}_k) \right ]^T \mathbf{z} = 0 \
這樣寫也明確了 \nabla f(\mathbf{x}_k) 是列向量。
于是要解出 \mathbf{z} 就可以使用線性代數的工具了：
\left[ \nabla f(\mathbf{x}_k) \right ]^T \mathbf{z}=-f(\mathbf{x}_k) \
此時上式就是一個線性方程組。而且仔細看會發現非常有意思：它其實只有一個方程。對這種方程我們可以考慮使用廣義逆矩陣：
\nabla f(\mathbf{x}_k)\left[ \nabla f(\mathbf{x}_k) \right ]^T \mathbf{z}=-\nabla f(\mathbf{x}_k)f(\mathbf{x}_k) \
此時 \nabla f(\mathbf{x}_k)\left[ \nabla f(\mathbf{x}_k) \right ]^T 就是一個方陣了，于是左右乘以它的逆：
\mathbf{z}= - \left { \nabla f(\mathbf{x}_k)\left[ \nabla D f(\mathbf{x}_k) \right ]^T \right }^{-1} \left [ \nabla f(\mathbf{x}_k) \right ] f(\mathbf{x}_k) \
這樣一來，多元函數求根的迭代式就可以得到了：
\bbox[pink]{\mathbf{x}_{k+1}=\mathbf{x}_k-\left { \nabla f(\mathbf{x}_k)\left[ \nabla f(\mathbf{x}_k) \right ]^T \right }^{-1} \left [\nabla f(\mathbf{x}_k) \right ] f(\mathbf{x}_k) \tag{M-1}}
這時簡單回顧一下向量導數：
\nabla f\left( \mathbf{x}k\right) = \dfrac{\partial f}{\partial \mathbf{x}}\Big| {\mathbf{x}=\mathbf{x}{k}}= \left[ \dfrac{\partial f}{\partial x{1}},\dfrac{\partial f}{\partial x_{2}},\ldots ,\dfrac{\partial f}{\partial x_{n}}\right] ^{T}\Big|{\mathbf{x}=\mathbf{x}{k}} \
不過查了一圈資料似乎很少看到這種迭代格式。
一來可能是這種多元函數的求根問題相對討論得少，更多肯定還是討論的多元函數的極值問題。
二來從分析過程也可以看到，關鍵的一步推導中增量向量僅由一個方程確定，這顯然是會導致極大的不確定性。當然這一點也容易理解，比如一個二元函數 f(\mathbf{x})=x_12+x_22-2=0，它的解實際上是一個圓。而如果用數值方法來求解，最好的情況也只能是收斂到一個點，自然也就用處不大了。而在實際應用中理應是增加各種限制條件來使它的解唯一或者有限。
（2）向量值函數
很神奇地是對于向量函數這個問題反而變得更自然。考慮一個向量函數：
\mathbf{f}:R^n \to R^m \\mathbf{f}(\mathbf{x})=\left[f_1(\mathbf{x}),f_2(\mathbf{x}),\cdots,f_m(\mathbf{x})\right]^T \
也就是說，該函數的值是一個向量，它的每個分量都是一個多元函數。它的零點問題就變成：
\mathbf{f}(\mathbf{x}) = \mathbf{0} \
有了上面的經驗，再次考慮一階逼近：
\mathbf{f}(\mathbf{x}_k+\mathbf{z}) \approx \mathbf{f}(\mathbf{x}_k)+ \left ( \dfrac{\partial \mathbf{f}(\mathbf{x})}{\partial \mathbf{x}}\Big| _{\mathbf{x}=\mathbf{x}_k} \right ) \mathbf{z} \
這里故意寫成 \dfrac{\partial \mathbf{f}(\mathbf{x})}{\partial \mathbf{x}} 是因為可以很好地對應矩陣導數：
\dfrac{\partial \mathbf{f}(\mathbf{x})}{\partial \mathbf{x}}=\begin{pmatrix} \dfrac{\partial \mathbf{f}\left( \mathbf{x}\right) }{\partial x_{1}} \ \dfrac{\partial \mathbf{f}\left( \mathbf{x}\right) }{\partial x_{2}} \ \vdots \ \dfrac{\partial \mathbf{f}\left( \mathbf{x}\right) }{\partial x_{n}} \end{pmatrix} = \begin{pmatrix} \dfrac{\partial f_{1}\left( \mathbf{x}\right) }{\partial x_{1}} & \dfrac{\partial f_{1}\left( \mathbf{x}\right) }{\partial x_{2}}&\ldots & \dfrac{\partial f_{1}\left( \mathbf{x}\right) }{\partial x_{n}} \ \dfrac{\partial f_{2}\left( \mathbf{x}\right) }{\partial x_{1}} & \dfrac{\partial f_{2}\left( \mathbf{x}\right) }{\partial x_{2}}&\ldots & \dfrac{\partial f_{2}\left( \mathbf{x}\right) }{\partial x_{n}} \ \vdots & \vdots& \ddots & \vdots \ \dfrac{\partial f_{m}\left( \mathbf{x}\right) }{\partial x_{1}} & \dfrac{\partial f_{m}\left( \mathbf{x}\right) }{\partial x_{2}}&\ldots & \dfrac{\partial f_{m}\left( \mathbf{x}\right) }{\partial x_{n}} \end{pmatrix} \tag{J}
這下子就好辦了，它就是一個 Jacobian 矩陣，也可以簡記為：
\mathbf{J}k=\mathbf{J}{\mathbf{f}}(\mathbf{x}k)= \begin{pmatrix} \dfrac{\partial f{1}\left( \mathbf{x}\right) }{\partial x_{1}} & \dfrac{\partial f_{1}\left( \mathbf{x}\right) }{\partial x_{2}}&\ldots & \dfrac{\partial f_{1}\left( \mathbf{x}\right) }{\partial x_{n}} \ \dfrac{\partial f_{2}\left( \mathbf{x}\right) }{\partial x_{1}} & \dfrac{\partial f_{2}\left( \mathbf{x}\right) }{\partial x_{2}}&\ldots & \dfrac{\partial f_{2}\left( \mathbf{x}\right) }{\partial x_{n}} \ \vdots & \vdots& \ddots & \vdots \ \dfrac{\partial f_{m}\left( \mathbf{x}\right) }{\partial x_{1}} & \dfrac{\partial f_{m}\left( \mathbf{x}\right) }{\partial x_{2}}&\ldots & \dfrac{\partial f_{m}\left( \mathbf{x}\right) }{\partial x_{n}} \end{pmatrix} \
于是令一階逼近式為 \mathbf{0}，得到
\mathbf{J}_k \mathbf{z}=-\mathbf{f}(\mathbf{x}_k) \
注意到這個線性方程組的系數矩陣也不一定是方陣，因此仍然要使用廣義逆：
\mathbf{z}=-\left ( \mathbf{J}_k^T\mathbf{J}_k \right )^{-1} \mathbf{J}_k^T\mathbf{f}(\mathbf{x}_k) \
于是完整迭代式可以寫出來：
\bbox[pink]{\mathbf{x}{k+1}=\mathbf{x}{k}-\left ( \mathbf{J}_k^T\mathbf{J}_k \right )^{-1} \mathbf{J}_k^T\mathbf{f}(\mathbf{x}_k) \tag{M-2}}
這個式子看起來就比較眼熟了。
當然如果 Jacobian 矩陣是方陣且可逆的話自然就可以直接求逆：
\mathbf{x}{k+1}=\mathbf{x}{k}-\mathbf{J}_k^{-1} \mathbf{f}(\mathbf{x}_k) \
不過這個情況，顯然也就不具有一般性了。
四、多元函數的極值問題
多元函數的極值可能是我們討論最多的問題了。尤其現在的機器學習訓練問題本質上就是一個多元函數的極值問題。這里我們只討論標量函數的情況。
類似第二小節，最簡單的思路是先求出一階導，再令一階導為 0。這里再明確一下：標量對向量的導數是向量。 于是寫出一階導的完整形式：
\mathbf{g}(\mathbf{x})=\nabla f(\mathbf{x})=\dfrac{\partial f(\mathbf{x})}{\partial \mathbf{x}}=\begin{pmatrix} \dfrac{\partial f\left( \mathbf{x}\right) }{\partial x_{1}} \ \dfrac{\partial f\left( \mathbf{x}\right) }{\partial x_{2}} \ \vdots \ \dfrac{\partial f\left( \mathbf{x}\right) }{\partial x_{n}} \end{pmatrix} \
仍然考慮一階導的一階逼近：
\mathbf{g}(\mathbf{x}_k+\mathbf{z}) \approx \mathbf{g}(\mathbf{x}_k)+ \left ( \dfrac{\partial \mathbf{g}(\mathbf{x})}{\partial \mathbf{x}}\Big| _{\mathbf{x}=\mathbf{x}_k} \right ) \mathbf{z} \
那么此時就可以令這個逼近式等于 \mathbf{0}，得到：
\left ( \dfrac{\partial \mathbf{g}(\mathbf{x})}{\partial \mathbf{x}}\Big| _{\mathbf{x}=\mathbf{x}_k} \right ) \mathbf{z}=-\mathbf{g}(\mathbf{x}_k) \
接下來就很自然了：
\mathbf{z}=-\left ( \dfrac{\partial \mathbf{g}(\mathbf{x})}{\partial \mathbf{x}}\Big| _{\mathbf{x}=\mathbf{x}_k} \right )^{-1}\mathbf{g}(\mathbf{x}_k) \
那么這時候有意思的就來了，一階導對向量 \mathbf{x} 的導數又變成了矩 Jacobian 陣：
\mathbf{J}{\mathbf{g}}(\mathbf{x})=\begin{pmatrix} \dfrac{\partial g{1}\left( \mathbf{x}\right) }{\partial x_{1}} & \dfrac{\partial g_{1}\left( \mathbf{x}\right) }{\partial x_{2}}&\ldots & \dfrac{\partial g_{1}\left( \mathbf{x}\right) }{\partial x_{n}} \ \dfrac{\partial g_{2}\left( \mathbf{x}\right) }{\partial x_{1}} & \dfrac{\partial g_{2}\left( \mathbf{x}\right) }{\partial x_{2}}&\ldots & \dfrac{\partial g_{2}\left( \mathbf{x}\right) }{\partial x_{n}} \ \vdots & \vdots& \ddots & \vdots \ \dfrac{\partial g_{n}\left( \mathbf{x}\right) }{\partial x_{1}} & \dfrac{\partial g_{n}\left( \mathbf{x}\right) }{\partial x_{2}}&\ldots & \dfrac{\partial g_{n}\left( \mathbf{x}\right) }{\partial x_{n}} \end{pmatrix} \
而對于函數 f (\mathbf{x} ) 又最做 Hessian 矩陣
\nabla^2 f(\mathbf{x})=\mathbf{H}_f (\mathbf{x} ) = \begin{bmatrix} \frac{\partial^2 f}{\partial^2 x_1} & \frac{\partial^2 f}{\partial x_1\partial x_2} & \ldots & \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 \partial x_n} \ \vdots &&&\vdots \ \frac{\partial^2 f}{\partial x_n \partial x_1} & \ldots & \ldots & \frac{\partial^2 f}{\partial^2 x_n} \end{bmatrix}. \
這個形式很顯然看起來就比較討厭了。
同樣地，由這兩個記號就可以寫出多元標量函數極值問題的牛頓迭代格式了：
\bbox[pink]{\mathbf{x}{k+1}=\mathbf{x}{k}-\mathbf{H}f ^{-1} (\mathbf{x}{k} ) \mathbf{g}(\mathbf{x}k)=\mathbf{x}{k}-\mathbf{H}f ^{-1} (\mathbf{x}{k} ) \nabla f(\mathbf{x}_k) \tag{M-3}}
這里就注意到一個很有意思的事情，Hessian 矩陣一定是方陣，所以可以直接考慮它的逆。當然如果不可逆時也會考慮一些技術處理。而此時廣義逆就不一定有用了，更多的情況是在其對角線加一個正數，這種方法稱為正則化。
\bbox[pink]{\mathbf{x}{k+1}=\mathbf{x}{k}-\left ( \mathbf{H}f ^{-1} (\mathbf{x}{k}) +\mu \mathbf{I} \right ) \nabla f(\mathbf{x}_k) \tag{M-4}}
到這還沒完，比著一元函數的牛頓迭代我們還有一種思路就是直接考慮二階逼近。而上面這些鋪墊其實已經把二階展開里要用的麻煩符號定義完了，所以直接寫：
f(\mathbf{x}_k+\mathbf{z})\approx f(\mathbf{x}k)+\left ( \dfrac{\partial f(\mathbf{x})}{\partial \mathbf{x}}\Big|{\mathbf{x}=\mathbf{x}k} \right )^T \mathbf{z} +\mathbf{z}^T\left ( \dfrac{\partial^2 f(\mathbf{x})}{\partial \mathbf{x}^2 }\Big|{\mathbf{x}=\mathbf{x}_k} \right ) \mathbf{z} \
可以寫成：
\bbox[pink]{f(\mathbf{x}_k+\mathbf{z})\approx f(\mathbf{x}_k)+\left ( \nabla f(\mathbf{x}_k) \right )^T \mathbf{z} +\frac{1}{2}\mathbf{z}^T\left ( \nabla^2 f(\mathbf{x}_k) \right ) \mathbf{z} \tag{M-5}}
那么如果直接對（M-5）式求導，就可以直接得到其迭代式（M-3）。
剛剛為了盡快進入牛頓迭代格式就直接跳過了另外一條線。事實上，上面在令其二階逼近的導數為零時，求其迭代格式的本質已經轉化成了一個線性方程組，我們簡單記成這樣：
\bbox[pink]{\mathbf{H}_k \mathbf{z}=-\mathbf{g}_k \tag{M-6}}
其實質就是如何快速求出 \mathbf{z} 的值。那么這樣一來，各種快速的求解方法就可以用了。對于較大的矩陣最常用的方法就是 CG 梯度，在《Matrix Computation》一書中就有介紹。不過在一些常見的數值計算工具中（比如 Matlab），似乎對于不太大的矩陣直接用 LU 分解還快一些。
五、最后說兩種簡化版
要說簡化，自然是要把最討厭的東西給簡化掉。上面已經說了，Hessian 矩陣 \mathbf{H} 最討厭。在一些特殊的問題中它是可以被簡化的。這種問題就是
(1) 非線性最小二乘問題
\ell(\mathbf{w}) = \frac{1}{2}| f(\mathbf{X};\mathbf{w})-\mathbf{Y} |^2= \frac{1}{2}\mathbf{e}^2 \
這里簡單作一下說明，用符號 \ell 對應現在流行的說法“損失函數”，同時根據其物理意義也可以叫 “最小平方誤差函數”，由于要求導所以乘以 \frac{1}{2} 后面討論起來比較方便。\mathbf{X} 是輸入變量組成的矩陣，\mathbf{Y} 是輸出變量組成的矩陣，\mathbf{w} 是全體參數組成的向量，\mathbf{e} 是誤差值組成的向量。 那么這時就直接套用一下二階展開（M-5）式：
\ell(\mathbf{w}_k+\mathbf{z})\approx \ell(\mathbf{w}_k)+\left ( \nabla \ell(\mathbf{w}_k) \right )^T \mathbf{z} +\frac{1}{2}\mathbf{z}^T\left ( \nabla^2 \ell(\mathbf{w}_k) \right ) \mathbf{z} \
由于這個最小平方誤差函數的定義比較特殊，因此它的二階展開里的兩個偏導可以寫得比較簡單：
\nabla \ell(\mathbf{w})=\dfrac{\partial \left ( \frac{1}{2}\mathbf{e}^2 \right ) }{\partial \mathbf{w}} = \left ( \dfrac{\partial \mathbf{e} }{\partial \mathbf{w}} \right )^T\mathbf{e} =\mathbf{J}^T_{\mathbf{e}}(\mathbf{w})\mathbf{e} \
這個 Jacobian 矩陣就是（J）式中的定義。這里注意下向量的形狀，始終保證最后的結果是列向量。
二階導自然就是對一階再求一次導：
\nabla^2 \ell(\mathbf{w})=\nabla \left ( \nabla\ell(\mathbf{w}) \right )= \frac{\partial\left ( \left ( \dfrac{\partial \mathbf{e} }{\partial \mathbf{w}} \right )^T\mathbf{e} \right ) }{\partial \mathbf{w}} \
這時很有意思的就來了，向量對向量求導的乘法法則（分配律）仍然成立：
\frac{\partial\left ( \left ( \dfrac{\partial \mathbf{e} }{\partial \mathbf{w}} \right )^T\mathbf{e} \right ) }{\partial \mathbf{w}} = \left ( \dfrac{\partial \mathbf{e} }{\partial \mathbf{w}} \right )^T \left ( \dfrac{\partial \mathbf{e} }{\partial \mathbf{w}} \right ) +\dfrac{\partial^2 \mathbf{e} }{\partial \mathbf{w}^2}\mathbf{e} \
于是它的二階導又可以寫成：
\nabla^2 \ell(\mathbf{w})= \mathbf{J}^T_{\mathbf{e}}(\mathbf{w})\mathbf{J}_{\mathbf{e}}(\mathbf{w}) +\dfrac{\partial^2 \mathbf{e} }{\partial \mathbf{w}^2}\mathbf{e} \
那么直接對應代回（M-3）式就有：
\bbox[pink]{\mathbf{w}{k+1}=\mathbf{w}{k}- \left ( (\mathbf{J}^T_{\mathbf{e}}(\mathbf{w})\mathbf{J}{\mathbf{e}}(\mathbf{w}) +\dfrac{\partial^2 \mathbf{e} }{\partial \mathbf{w}^2}\mathbf{e}) ^{-1} \mathbf{J}^T{\mathbf{e}}(\mathbf{w})\mathbf{e} \right ) \Big|_{\mathbf{w}=\mathbf{w}_k} \tag{LS-1}}
這玩意誰看誰煩！
（2）高斯-牛頓法
一位法國兒童可能也是懷著同樣的心情，于是發了一篇文章：Theoria motus corporum coelestium in sectionibus conicis solem ambientum (1809), pp. 179-180.
在這篇文章里，他大概提到：\nabla^2 \ell(\mathbf{w}) 中的二階項應該很小，所以直接省略，于是再對符號簡化一下，這個 Hessian 近似就可以寫成：
\nabla^2 \ell(\mathbf{w})= \mathbf{J}^T\mathbf{J} \
這下子就非常簡潔了：
KaTeX parse error: Undefined control sequence: \bbox at position 2: \?b?b?o?x?[pink]{\mathbf{…
這個式子就是著名的 Gauss-Newton 法的迭代公式。
（3）Levenberg-Marquardt
在討論一般的極值問題時就提到過一個問題，那個 Hessian 矩陣是不一定可逆的。用 Gauss-Newton 方法之后，其實這種不可逆可以在一定程度上得到緩解。因為不難證明對于形如 \mathbf{J}^T\mathbf{J} 的矩陣通常是半正定的。但在實際計算過程中也難免會有數值上的病態出現，于是就出現了各種的改進。
最簡單直接的思路就是在 J T J \mathbf{J}^T\mathbf{J} JTJ的對角線上加上一個正數：
而這個格式就是后來更常用的 Levenberg-Marquardt 算法。
剩下還有一些其它的加權方式基本都是在這個基礎上進行改進，不再多說了。
寫完了發現其實大部分時間還是在直接敲。主要還是手寫工具反應太慢了。初次寫復雜的式子的時候還有點用，如果改動的話還不如直接整代碼。加起來還是整了好幾個小時。
至于理論篇，再說吧。想想就覺得累。
關于牛頓迭代有個比較好的綜述，懶得翻譯了將就看看：
The name “Newton’s method” is derived from Isaac Newton’s description of a special case of the method in De analysi per aequationes numero terminorum infinitas (written in 1669, published in 1711 by William Jones) and in De metodis fluxionum et serierum infinitarum (written in 1671, translated and published as Method of Fluxions in 1736 by John Colson). However, his method differs substantially from the modern method given above. Newton applied the method only to polynomials, starting with an initial root estimate and extracting a sequence of error corrections. He used each correction to rewrite the polynomial in terms of the remaining error, and then solved for a new correction by neglecting higher-degree terms. He did not explicitly connect the method with derivatives or present a general formula. Newton applied this method to both numerical and algebraic problems, producing Taylor series in the latter case.
Newton may have derived his method from a similar but less precise method by Vieta. The essence of Vieta’s method can be found in the work of the Persian mathematician Sharaf al-Din al-Tusi, while his successor Jamshīd al-Kāshī used a form of Newton’s method to solve xP ? N = 0 to find roots of N (Ypma 1995). A special case of Newton’s method for calculating square roots was known since ancient times and is often called the Babylonian method.
Newton’s method was used by 17th-century Japanese mathematician Seki Kōwa to solve single-variable equations, though the connection with calculus was missing.[1]
Newton’s method was first published in 1685 in A Treatise of Algebra both Historical and Practical by John Wallis.[2] In 1690, Joseph Raphson published a simplified description in Analysis aequationum universalis.[3] Raphson also applied the method only to polynomials, but he avoided Newton’s tedious rewriting process by extracting each successive correction from the original polynomial. This allowed him to derive a reusable iterative expression for each problem. Finally, in 1740, Thomas Simpson described Newton’s method as an iterative method for solving general nonlinear equations using calculus, essentially giving the description above. In the same publication, Simpson also gives the generalization to systems of two equations and notes that Newton’s method can be used for solving optimization problems by setting the gradient to zero.
Arthur Cayley in 1879 in The Newton–Fourier imaginary problem was the first to notice the difficulties in generalizing Newton’s method to complex roots of polynomials with degree greater than 2 and complex initial values. This opened the way to the study of the theory of iterations of rational functions.
