说说牛顿迭代 -- 方法篇

说说牛顿迭代 – 方法篇

写这个笔记主要是最近老在考虑最优化问题。今天刚好发现一个不错的手写公式的工具，加上前几天又发现Win10的Windows Ink比我想象得好用，于是来描几笔。主要是想试试这样写起来速度如何。

现在高数也要介绍牛顿法了，一般都是从几何上直接以“切线法”直接引入的。这种引入方式的确很简单，但如果脱离深入一点的分析就不大容易解释后面的事情。所以就在想怎么更直接地从分析的角度来一条线贯穿，把整个过程说一说。

直接开始：

一、牛顿迭代解非线性方程

考虑一个非线性方程：
f(x)=0(*)f(x)=0 \tag{*} f(x)=0(*)
这个方程没有公式可以直接解，于是考虑迭代法。其思想就是构造一个迭代公式：
xk+1=φ(xk)(1)x_{k+1}=\varphi(x_k) \tag{1} xk+1=φ(xk)(1)
如果k→∞k\to \inftyk→∞ 时xk+1→x∗x_{k+1}\to x^*xk+1→x∗，那么就可以不断地利用这个公式以 xk+1x_{k+1}xk+1 去逼近方程的根 x∗x^*x∗。

一般来说更常见的迭代格式是：
xk+1=xk−Δ(xk)(2)x_{k+1}=x_k-\Delta(x_k) \tag{2} xk+1=xk−Δ(xk)(2)
比如梯度下降就是这种格式。

那么问题就很自然了：如何构造这个迭代公式 Δ(xk)\Delta(x_k)Δ(xk)？

先从最简单的方向考虑：如果 xkx_kxk 已经十分接近 x∗x^*x∗ 时，那么这个迭代的长度 Δ(xk)\Delta(x_k)Δ(xk) 应该是多少？

这时候问题就简单得多了，因为当 xkx_kxk 已经十分接近 x∗x^*x∗ 时我们可以认为再移动 Δ(xk)\Delta(x_k)Δ(xk) 就可以使得
f(xk+1)=f(xk−Δ(xk))=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）式得到Δ(xk)\Delta(x_k)Δ(xk)？

那么很自然地，如果能通过（2）（3）式导出一个关于 Δ(xk)\Delta(x_k)Δ(xk) 的迭代式就行了。从（3）式很容易想到如果把第二个式子展开，就可以得到一个 Δ(xk)\Delta(x_k)Δ(xk) 的解析式，从而就能得到它的显式表达。这样很容易想到Taylor展开（一阶逼近）：
f(xk−Δ(xk))≈f(xk)+f′(xk)⋅(−Δ(xk))(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）式就可以得到牛顿迭代的完整表达式：
Δ(xk)=f(xk)f′(xk)\Delta \left( x_{k}\right) =\dfrac{f\left( x_{k}\right) }{f'\left( x_{k}\right) } Δ(xk)=f′(xk)f(xk)
于是得到牛顿迭代的完整格式：
xk+1=xk−f(xk)f′(xk)(5)x_{k+1}=x_k-\dfrac{f\left( x_{k}\right) }{f'\left( x_{k}\right) \tag{5}} xk+1=xk−f′(xk)f(xk)(5)
当然这里是基于如果 xkx_kxk 已经十分接近 x∗x^*x∗ 时的情况推导，直接利用上面的分析还不太容易解释从何意一个可行的点（根的附近）出发，（5）式一定能让 xkx_kxk 不断地接近 x∗x^*x∗。

那么这时结合几何来简单看看就好说了：

结合上图的标记就很容易知道，无论如何这个迭代公式总会使 xkx_kxk 不断接近 x∗x^*x∗。

这里再记几笔：

1、如果不是上述的单调情况，那么牛顿迭代是不会收敛到图中所示的 x∗x^*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)
那么自然就有它的迭代公式：
xk+1=xk−g(xx)g′(xk)=xk−f′(xk)f′′(xk)(6)x_{k+1}=x_{k}-\dfrac{g\left( x_{x}\right) }{g'\left( x_{k}\right) }=x_{k}-\dfrac{f'\left( x_{k}\right) }{f''\left( x_{k}\right) \tag{6}} xk+1=xk−g′(xk)g(xx)=xk−f′′(xk)f′(xk)(6)
这种理解十分便于记忆，所以还是十分推荐的。

但另外一种思路更适合于推广，再来介绍一下。

基本的思考方式和上面解非线性方程一样，我们仍然考虑 xkx_kxk 在极值点 x∗∗x^{**}x∗∗ 附近，那么问题就变成了：

给多长的增量Δ(xk)\Delta \left( x_{k}\right)Δ(xk)，能使得点 xk+1=xk−Δ(xk)x_{k+1}=x_{k}-\Delta \left( x_{k}\right)xk+1=xk−Δ(xk) 为极值点（或者近似的极值点），换言之希望此时 f(xk+1)≈0f(x_{k+1})\approx 0f(xk+1)≈0。

类似地，我们仍然考虑它的Taylor展开，不过此时展开到二阶：
f(xk+1)≈f(xk)+f′(xk)⋅(−Δ(xk))+12f′′(xk)⋅(−Δ(xk))2f\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))+21f′′(xk)⋅(−Δ(xk))2
为了方便描述我们记 xk+1=xk+zx_{k+1}=x_{k}+zxk+1=xk+z，这样二阶展开就可以写成：
f(xk+1)≈f(xk)+f′(xk)⋅z+12f′′(xk)⋅z2=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+21f′′(xk)⋅z2=g(z)(7)
继续考虑 g(z)g(z)g(z) 的极值：
dg(z)dz=f′(xx)+f′′(xk)⋅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
可以解出：
z=−f′(xx)f′′(xk)z=-\dfrac{f'\left( x_{x}\right) }{f''\left( x_{k}\right) } z=−f′′(xk)f′(xx)
这和（6）式中的迭代长度是一样的。

注1：这里有个很有意思的事情，如果我们只考虑一阶逼近
f(xk+1)≈f(xk)+f′(xk)⋅zf\left( x_{k+1}\right) \approx f\left( x_{k}\right) +f'\left( x_{k}\right)\cdot z f(xk+1)≈f(xk)+f′(xk)⋅z
那么让它对 zzz 求导来求极值就直接得到： f′(xk)=0f'\left( x_{k}\right)=0f′(xk)=0

很显然这种做法就不太可取了。由此也可以看出一阶逼近的粗糙，以至于在讨论极值问题时就不适合直接使用了。

注2： 有同学可能也会问，为什么不直接考虑该函数任意值在 xkx_kxk 附近的展开？

也就是说考虑：
f(x)≈f(xk)+f′(xk)(x−xk)+12f′′(xk)(x−xk)2f\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)≈f(xk)+f′(xk)(x−xk)+21f′′(xk)(x−xk)2
此时直接考虑 f′(x)=0f'(x)=0f′(x)=0，那么容易得到：
f′(xk)+f′′(xk)(x−xk)≈0f'\left( x_{k}\right) +f''\left( x_{k}\right) \left( x-x_{k}\right) \approx 0 f′(xk)+f′′(xk)(x−xk)≈0
于是：
x≈xk−f′(xk)f′′(xk)(8)x \approx x_{k}-\dfrac{f'\left( x_{k}\right) }{f''\left( x_{k}\right) } \tag{8} x≈xk−f′′(xk)f′(xk)(8)
这就有意思了，前面的牛顿迭代式这里都是严格的等号。而上式却是不等号。这里就可以再借题发挥一下：

上面的讨论中（7）式是直接令逼近式严格等于 g(z)g(z)g(z)，因此得到的迭代步长就是严格等于。
而（8）式是直接令 f′(x)=0f'(x)=0f′(x)=0，那么将二阶近似代回计算时自然就要变成约等号。那么得到的值就自然是约等于。不过借助（8）式也方便说明一些问题。比如：此时得到的 xxx 实际上就是极值点，那么如果 xkx_kxk 本身就离极值点很近，这个式子就更加精确；相反，如果 xkx_kxk 离极值本身还很远，那么这个式子就十分粗糙。

注3： 上面的（7）式也是非常有意思的，在极值问题中这个式子的几何意义就十分明确了。先简单移项：
f(xk+1)−f(xk)≈f′(xk)⋅z+12f′′(xk)⋅z2f\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 f(xk+1)−f(xk)≈f′(xk)⋅z+21f′′(xk)⋅z2
这时右端的导数也是 dg(z)dz\dfrac{dg\left( z\right) }{dz}dzdg(z)，那么此时令dg(z)dz=0\dfrac{dg\left( z\right) }{dz}=0dzdg(z)=0 从几何上讲意思就是：取改 zzz 为多少能使得下一步迭代时 f(xk+1)f\left( x_{k+1}\right)f(xk+1) 相对上一步迭代时 f(xk)f\left( x_{k}\right)f(xk) 的改变最大？

这里就能很明显体会到，用上述思路其实就是在整体寻优的过程中构造了一个局部寻优的问题。

接下来讨论多元函数的零点问题

三、多元函数的零点问题

多元的情形自然复杂，要完全写对需要非常仔细。但过于仔细又容易影响总体思路，所以先不管式子是否正确，直接对照一元情形凭感觉把推导写完再说：

另外我们这里考虑的是普通微积分里讨论的多元函数，也就是说f:Rn→Rf:R^n\to Rf:Rn→R。

（1）标量函数

回顾一元的情形，我们仍然是考虑在零点附近的一个 xk\mathbf{x}_kxk ，给它一个增量 z\mathbf{z}z，使得 xk+1=xk+z\mathbf{x}_{k+1}=\mathbf{x}_k+\mathbf{z}xk+1=xk+z 接近零点 x∗∗\mathbf{x}^{**}x∗∗。类似一元的情形，考虑函数的一阶逼近：
f(xk+z)≈f(xk)+∇f(xk)⋅zf(\mathbf{x}_k+\mathbf{z}) \approx f(\mathbf{x}_k)+\nabla f(\mathbf{x}_k)\cdot\mathbf{z} f(xk+z)≈f(xk)+∇f(xk)⋅z
同样地，考虑逼近式等于 000，于是有：
f(xk)+∇f(xk)⋅z=0f(\mathbf{x}_k)+\nabla f(\mathbf{x}_k)\cdot\mathbf{z} = 0 f(xk)+∇f(xk)⋅z=0
到这里就要注意一下了:

f(xk)f(\mathbf{x}_k)f(xk) 是标量
∇f(xk)⋅z\nabla f(\mathbf{x}_k)\cdot\mathbf{z}∇f(xk)⋅z 也必须是标量
但 z\mathbf{z}z 是向量，所以∇f(xk)\nabla f(\mathbf{x}_k)∇f(xk) 也是向量

那么我们先把写法再规范一下：
f(xk)+[∇f(xk)]Tz=0f(\mathbf{x}_k)+ \left[ \nabla f(\mathbf{x}_k) \right ]^T \mathbf{z} = 0 f(xk)+[∇f(xk)]Tz=0

这样写也明确了 ∇f(xk)\nabla f(\mathbf{x}_k)∇f(xk) 是列向量。

于是要解出 z\mathbf{z}z 就可以使用线性代数的工具了：
[∇f(xk)]Tz=−f(xk)\left[ \nabla f(\mathbf{x}_k) \right ]^T \mathbf{z}=-f(\mathbf{x}_k) [∇f(xk)]Tz=−f(xk)
此时上式就是一个线性方程组。而且仔细看会发现非常有意思：它其实只有一个方程。对这种方程我们可以考虑使用广义逆矩阵：
∇f(xk)[∇f(xk)]Tz=−∇f(xk)f(xk)\nabla f(\mathbf{x}_k)\left[ \nabla f(\mathbf{x}_k) \right ]^T \mathbf{z}=-\nabla f(\mathbf{x}_k)f(\mathbf{x}_k) ∇f(xk)[∇f(xk)]Tz=−∇f(xk)f(xk)
此时 ∇f(xk)[∇f(xk)]T\nabla f(\mathbf{x}_k)\left[ \nabla f(\mathbf{x}_k) \right ]^T∇f(xk)[∇f(xk)]T 就是一个方阵了，于是左右乘以它的逆：
z=−{∇f(xk)[∇Df(xk)]T}−1[∇f(xk)]f(xk)\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) z=−{∇f(xk)[∇Df(xk)]T}−1[∇f(xk)]f(xk)
这样一来，多元函数求根的迭代式就可以得到了：
xk+1=xk−{∇f(xk)[∇f(xk)]T}−1[∇f(xk)]f(xk)(M-1)\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} xk+1=xk−{∇f(xk)[∇f(xk)]T}−1[∇f(xk)]f(xk)(M-1)
这时简单回顾一下向量导数：
∇f(xk)=∂f∂x∣x=xk=[∂f∂x1,∂f∂x2,…,∂f∂xn]T∣x=xk\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(xk)=∂x∂f∣∣∣x=xk=[∂x1∂f,∂x2∂f,…,∂xn∂f]T∣∣∣x=xk

不过查了一圈资料似乎很少看到这种迭代格式。

一来可能是这种多元函数的求根问题相对讨论得少，更多肯定还是讨论的多元函数的极值问题。

二来从分析过程也可以看到，关键的一步推导中增量向量仅由一个方程确定，这显然是会导致极大的不确定性。当然这一点也容易理解，比如一个二元函数 f(x)=x12+x22−2=0f(\mathbf{x})=x_1^2+x_2^2-2=0f(x)=x12+x22−2=0，它的解实际上是一个圆。而如果用数值方法来求解，最好的情况也只能是收敛到一个点，自然也就用处不大了。而在实际应用中理应是增加各种限制条件来使它的解唯一或者有限。

（2）向量值函数

很神奇地是对于向量函数这个问题反而变得更自然。考虑一个向量函数：
f:Rn→Rm\mathbf{f}:R^n \to R^m f:Rn→Rm

f(x)=[f1(x),f2(x),⋯,fm(x)]T\mathbf{f}(\mathbf{x})=\left[f_1(\mathbf{x}),f_2(\mathbf{x}),\cdots,f_m(\mathbf{x})\right]^T f(x)=[f1(x),f2(x),⋯,fm(x)]T

也就是说，该函数的值是一个向量，它的每个分量都是一个多元函数。它的零点问题就变成：
f(x)=0\mathbf{f}(\mathbf{x}) = \mathbf{0} f(x)=0
有了上面的经验，再次考虑一阶逼近：
f(xk+z)≈f(xk)+(∂f(x)∂x∣x=xk)z\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} f(xk+z)≈f(xk)+(∂x∂f(x)∣∣∣x=xk)z
这里故意写成 ∂f(x)∂x\dfrac{\partial \mathbf{f}(\mathbf{x})}{\partial \mathbf{x}}∂x∂f(x) 是因为可以很好地对应矩阵导数：
∂f(x)∂x=(∂f(x)∂x1∂f(x)∂x2⋮∂f(x)∂xn)=(∂f1(x)∂x1∂f1(x)∂x2…∂f1(x)∂xn∂f2(x)∂x1∂f2(x)∂x2…∂f2(x)∂xn⋮⋮⋱⋮∂fm(x)∂x1∂fm(x)∂x2…∂fm(x)∂xn)(J)\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} ∂x∂f(x)=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛∂x1∂f(x)∂x2∂f(x)⋮∂xn∂f(x)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛∂x1∂f1(x)∂x1∂f2(x)⋮∂x1∂fm(x)∂x2∂f1(x)∂x2∂f2(x)⋮∂x2∂fm(x)……⋱…∂xn∂f1(x)∂xn∂f2(x)⋮∂xn∂fm(x)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞(J)
这下子就好办了，它就是一个Jacobian矩阵，也可以简记为：
Jk=Jf(xk)=(∂f1(x)∂x1∂f1(x)∂x2…∂f1(x)∂xn∂f2(x)∂x1∂f2(x)∂x2…∂f2(x)∂xn⋮⋮⋱⋮∂fm(x)∂x1∂fm(x)∂x2…∂fm(x)∂xn)\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} Jk=Jf(xk)=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛∂x1∂f1(x)∂x1∂f2(x)⋮∂x1∂fm(x)∂x2∂f1(x)∂x2∂f2(x)⋮∂x2∂fm(x)……⋱…∂xn∂f1(x)∂xn∂f2(x)⋮∂xn∂fm(x)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞
于是令一阶逼近式为 0\mathbf{0}0，得到
Jkz=−f(xk)\mathbf{J}_k \mathbf{z}=-\mathbf{f}(\mathbf{x}_k) Jkz=−f(xk)
注意到这个线性方程组的系数矩阵也不一定是方阵，因此仍然要使用广义逆：
z=−(JkTJk)−1JkTf(xk)\mathbf{z}=-\left ( \mathbf{J}_k^T\mathbf{J}_k \right )^{-1} \mathbf{J}_k^T\mathbf{f}(\mathbf{x}_k) z=−(JkTJk)−1JkTf(xk)
于是完整迭代式可以写出来：
xk+1=xk−(JkTJk)−1JkTf(xk)(M-2)\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} xk+1=xk−(JkTJk)−1JkTf(xk)(M-2)
这个式子看起来就比较眼熟了。

当然如果 Jacobian矩阵是方阵且可逆的话自然就可以直接求逆：
xk+1=xk−Jk−1f(xk)\mathbf{x}_{k+1}=\mathbf{x}_{k}-\mathbf{J}_k^{-1} \mathbf{f}(\mathbf{x}_k) xk+1=xk−Jk−1f(xk)
不过这个情况，显然也就不具有一般性了。

四、多元函数的极值问题

多元函数的极值可能是我们讨论最多的问题了。尤其现在的机器学习训练问题本质上就是一个多元函数的极值问题。这里我们只讨论标量函数的情况。

类似第二小节，最简单的思路是先求出一阶导，再令一阶导为0。这里再明确一下：标量对向量的导数是向量。 于是写出一阶导的完整形式：
g(x)=∇f(x)=∂f(x)∂x=(∂f(x)∂x1∂f(x)∂x2⋮∂f(x)∂xn)\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} g(x)=∇f(x)=∂x∂f(x)=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛∂x1∂f(x)∂x2∂f(x)⋮∂xn∂f(x)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞
仍然考虑一阶导的一阶逼近：
g(xk+z)≈g(xk)+(∂g(x)∂x∣x=xk)z\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} g(xk+z)≈g(xk)+(∂x∂g(x)∣∣∣x=xk)z
那么此时就可以令这个逼近式等于 0\mathbf{0}0，得到：
(∂g(x)∂x∣x=xk)z=−g(xk)\left ( \dfrac{\partial \mathbf{g}(\mathbf{x})}{\partial \mathbf{x}}\Big| _{\mathbf{x}=\mathbf{x}_k} \right ) \mathbf{z}=-\mathbf{g}(\mathbf{x}_k) (∂x∂g(x)∣∣∣x=xk)z=−g(xk)
接下来就很自然了：
z=−(∂g(x)∂x∣x=xk)−1g(xk)\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) z=−(∂x∂g(x)∣∣∣x=xk)−1g(xk)
那么这时候有意思的就来了，一阶导对向量 x\mathbf{x}x 的导数又变成了矩Jacobian阵：
Jg(x)=(∂g1(x)∂x1∂g1(x)∂x2…∂g1(x)∂xn∂g2(x)∂x1∂g2(x)∂x2…∂g2(x)∂xn⋮⋮⋱⋮∂gn(x)∂x1∂gn(x)∂x2…∂gn(x)∂xn)\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} Jg(x)=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛∂x1∂g1(x)∂x1∂g2(x)⋮∂x1∂gn(x)∂x2∂g1(x)∂x2∂g2(x)⋮∂x2∂gn(x)……⋱…∂xn∂g1(x)∂xn∂g2(x)⋮∂xn∂gn(x)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞
而对于函数 f(x)f (\mathbf{x} )f(x) 又最做 Hessian 矩阵
∇2f(x)=Hf(x)=[∂2f∂2x1∂2f∂x1∂x2…∂2f∂x1∂xn∂2f∂x2∂x1∂2f∂x2∂xn⋮⋮∂2f∂xn∂x1……∂2f∂2xn].\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}. ∇2f(x)=Hf(x)=⎣⎢⎢⎢⎢⎡∂2x1∂2f∂x2∂x1∂2f⋮∂xn∂x1∂2f∂x1∂x2∂2f………∂x1∂xn∂2f∂x2∂xn∂2f⋮∂2xn∂2f⎦⎥⎥⎥⎥⎤.

这个形式很显然看起来就比较讨厌了。

同样地，由这两个记号就可以写出多元标量函数极值问题的牛顿迭代格式了：
xk+1=xk−Hf−1(xk)g(xk)=xk−Hf−1(xk)∇f(xk)(M-3)\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} xk+1=xk−Hf−1(xk)g(xk)=xk−Hf−1(xk)∇f(xk)(M-3)
这里就注意到一个很有意思的事情，Hessian 矩阵一定是方阵，所以可以直接考虑它的逆。当然如果不可逆时也会考虑一些技术处理。而此时广义逆就不一定有用了，更多的情况是在其对角线加一个正数，这种方法称为正则化。
xk+1=xk−(Hf−1(xk)+μI)∇f(xk)(M-4)\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} xk+1=xk−(Hf−1(xk)+μI)∇f(xk)(M-4)
到这还没完，比着一元函数的牛顿迭代我们还有一种思路就是直接考虑二阶逼近。而上面这些铺垫其实已经把二阶展开里要用的麻烦符号定义完了，所以直接写：
f(xk+z)≈f(xk)+(∂f(x)∂x∣x=xk)Tz+zT(∂2f(x)∂x2∣x=xk)zf(\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} f(xk+z)≈f(xk)+(∂x∂f(x)∣∣∣x=xk)Tz+zT(∂x2∂2f(x)∣∣∣x=xk)z
可以写成：
f(xk+z)≈f(xk)+(∇f(xk))Tz+12zT(∇2f(xk))z(M-5)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} f(xk+z)≈f(xk)+(∇f(xk))Tz+21zT(∇2f(xk))z(M-5)
那么如果直接对（M-5）式求导，就可以直接得到其迭代式（M-3）。

刚刚为了尽快进入牛顿迭代格式就直接跳过了另外一条线。事实上，上面在令其二阶逼近的导数为零时，求其迭代格式的本质已经转化成了一个线性方程组，我们简单记成这样：
Hkz=−gk(M-6)\mathbf{H}_k \mathbf{z}=-\mathbf{g}_k \tag{M-6} Hkz=−gk(M-6)
其实质就是如何快速求出 z\mathbf{z}z 的值。那么这样一来，各种快速的求解方法就可以用了。对于较大的矩阵最常用的方法就是CG梯度，在《Matrix Computation》一书中就有介绍。不过在一些常见的数值计算工具中（比如Matlab），似乎对于不太大的矩阵直接用LU分解还快一些。

五、最后说两种简化版

要说简化，自然是要把最讨厌的东西给简化掉。上面已经说了，Hessian 矩阵 H\mathbf{H}H 最讨厌。在一些特殊的问题中它是可以被简化的。这种问题就是

(1) 非线性最小二乘问题

ℓ(w)=12∥f(X;w)−Y∥2=12e2\ell(\mathbf{w}) = \frac{1}{2}\| f(\mathbf{X};\mathbf{w})-\mathbf{Y} \|^2= \frac{1}{2}\mathbf{e}^2 ℓ(w)=21∥f(X;w)−Y∥2=21e2

这里简单作一下说明，用符号 ℓ\ellℓ 对应现在流行的说法“损失函数”，同时根据其物理意义也可以叫 “最小平方误差函数”，由于要求导所以乘以 12\frac{1}{2}21 后面讨论起来比较方便。X\mathbf{X}X 是输入变量组成的矩阵，Y\mathbf{Y}Y 是输出变量组成的矩阵，w\mathbf{w}w 是全体参数组成的向量，e\mathbf{e}e 是误差值组成的向量。那么这时就直接套用一下二阶展开（M-5）式：
ℓ(wk+z)≈ℓ(wk)+(∇ℓ(wk))Tz+12zT(∇2ℓ(wk))z\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} ℓ(wk+z)≈ℓ(wk)+(∇ℓ(wk))Tz+21zT(∇2ℓ(wk))z
由于这个最小平方误差函数的定义比较特殊，因此它的二阶展开里的两个偏导可以写得比较简单：
∇ℓ(w)=∂(12e2)∂w=(∂e∂w)Te=JeT(w)e\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} ∇ℓ(w)=∂w∂(21e2)=(∂w∂e)Te=JeT(w)e
这个Jacobian矩阵就是（J）式中的定义。这里注意下向量的形状，始终保证最后的结果是列向量。

二阶导自然就是对一阶再求一次导：
∇2ℓ(w)=∇(∇ℓ(w))=∂((∂e∂w)Te)∂w\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}} ∇2ℓ(w)=∇(∇ℓ(w))=∂w∂((∂w∂e)Te)
这时很有意思的就来了，向量对向量求导的乘法法则（分配律）仍然成立：
∂((∂e∂w)Te)∂w=(∂e∂w)T(∂e∂w)+∂2e∂w2e\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} ∂w∂((∂w∂e)Te)=(∂w∂e)T(∂w∂e)+∂w2∂2ee
于是它的二阶导又可以写成：
∇2ℓ(w)=JeT(w)Je(w)+∂2e∂w2e\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} ∇2ℓ(w)=JeT(w)Je(w)+∂w2∂2ee
那么直接对应代回（M-3）式就有：
wk+1=wk−((JeT(w)Je(w)+∂2e∂w2e)−1JeT(w)e)∣w=wk(LS-1)\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} wk+1=wk−((JeT(w)Je(w)+∂w2∂2ee)−1JeT(w)e)∣∣∣w=wk(LS-1)
这玩意谁看谁烦！

（2）高斯-牛顿法

一位法国儿童可能也是怀着同样的心情，于是发了一篇文章：Theoria motus corporum coelestium in sectionibus conicis solem ambientum (1809), pp. 179-180.

在这篇文章里，他大概提到：∇2ℓ(w)\nabla^2 \ell(\mathbf{w})∇2ℓ(w) 中的二阶项应该很小，所以直接省略，于是再对符号简化一下，这个Hessian近似就可以写成：
∇2ℓ(w)=JTJ\nabla^2 \ell(\mathbf{w})= \mathbf{J}^T\mathbf{J} ∇2ℓ(w)=JTJ
这下子就非常简洁了：
wk+1=wk−(JTJ)−1JTe(LS-2)\mathbf{w}_{k+1}=\mathbf{w}_{k}-\left ( \mathbf{J}^T\mathbf{J} \right )^{-1} \mathbf{J}^T\mathbf{e} \tag{LS-2} wk+1=wk−(JTJ)−1JTe(LS-2)
这个式子就是著名的Gauss-Newton法的迭代公式。

（3）Levenberg-Marquardt

在讨论一般的极值问题时就提到过一个问题，那个Hessian 矩阵是不一定可逆的。用Gauss-Newton方法之后，其实这种不可逆可以在一定程度上得到缓解。因为不难证明对于形如 JTJ\mathbf{J}^T\mathbf{J}JTJ 的矩阵通常是半正定的。但在实际计算过程中也难免会有数值上的病态出现，于是就出现了各种的改进。

最简单直接的思路就是在 JTJ\mathbf{J}^T\mathbf{J}JTJ 的对角线上加上一个正数：
wk+1=wk−(JTJ+μI)−1JTe(LS-3)\mathbf{w}_{k+1}=\mathbf{w}_{k}-\left ( \mathbf{J}^T\mathbf{J}+\mu\mathbf{I} \right )^{-1} \mathbf{J}^T\mathbf{e} \tag{LS-3} wk+1=wk−(JTJ+μI)−1JTe(LS-3)
而这个格式就是后来更常用的 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.

一些比较有用的资料：

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.

一些比较有用的资料：

牛顿法求向量函数的根