高等数学:牛顿迭代发解方程

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

文章目录

  • 一、牛顿迭代解非线性方程
  • 二、牛顿迭代解极值问题

直接开始:

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

考虑一个非线性方程:
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^* kxk+1x,那么就可以不断地利用这个公式以 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))+21f′′(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+21f′′(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.

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/42170.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

【深度学习】PyTorch快速入门

【深度学习】学习PyTorch基础 介绍PyTorch 深度学习框架是一种软件工具,旨在简化和加速构建、训练和部署深度学习模型的过程。深度学习框架提供了一系列的函数、类和工具,用于定义、优化和执行各种深度神经网络模型。这些框架帮助研究人员和开发人员专注…

漏洞+常见漏洞修复建议

一、漏洞级别 高级别漏洞:大部分远程和本地管理员权限漏洞 中级别漏洞:大部分普通用户权限、权限提升、读懂受限文件、远程和本杜漏洞、拒绝服务漏洞 低级别漏洞:大部分远程非授权文件存取、口令恢复、欺骗、信息泄露漏洞 二、漏洞扫描的…

Kotlin Lambda和高阶函数

Lambda和高阶函数 本文链接: 文章目录 Lambda和高阶函数 lambda输出(返回类型)深入探究泛型 inline原理探究 高阶函数集合、泛型自己实现Kotlin内置函数 扩展函数原理companion object 原理 > 静态内部类函数式编程 lambda 1、lambda的由…

Flink流批一体计算(13):PyFlink Tabel API之SQL DDL

1. TableEnvironment 创建 TableEnvironment from pyflink.table import Environmentsettings, TableEnvironment# create a streaming TableEnvironmentenv_settings Environmentsettings.in_streaming_mode()table_env TableEnvironment.create(env_settings)# or create…

嵌入式Linux开发实操(九):CAN接口开发

前言: CAN网络在汽车中的使用可以说相当广泛。而CAN网络需要的收发器最常用的就是NXP 的TJA1042: CAN网络:

605. 种花问题

链接 假设有一个很长的花坛,一部分地块种植了花,另一部分却没有。可是,花不能种植在相邻的地块上,它们会争夺水源,两者都会死去。给你一个整数数组 flowerbed 表示花坛,由若干 0 和 1 组成,其中…

8/16总结

WebSocket是双向通信协议,模拟Socket协议,可以双向发送或者接收信息 而Http是单向的 WebSocket是需要浏览器和服务器握手进行建立连接的 而http是浏览器发起向服务器的连接,服务器预先并不知道这个连接 WebSocket在建立握手时,数…

Python3内置函数大全

吐血整理 Python3内置函数大全 1.abs()函数2.all()函数3.any()函数4.ascii()函数5.bin()函数6.bool()函数7.bytes()函数8.challable()函数9.chr()函数10.classmethod()函数11.complex()函数12.complie()函数13.delattr()函数14.dict()函数15.dir()函数16.divmod()函数17.enumer…

注解@JsonInclude

注解JsonInclude 1. 注解由来 JsonInclude是一个用于Java类中字段或方法的注解,它来自于Jackson库。Jackson库是一个用于处理JSON数据的流行开源库,在Java对象和JSON之间进行序列化和反序列化时经常被使用。 2. 注解示例 下面是JsonInclude注解的一个…

【kubernetes】Pod控制器

目录 Pod控制器及其功用 pod控制器有多种类型 1、ReplicaSet ReplicaSet主要三个组件组成 2、Deployment 3、DaemonSet 4、StatefulSet 5、Job 6、Cronjob Pod与控制器之间的关系 1、Deployment 查看控制器配置 查看历史版本 2、SatefulSet 为什么要有headless&…

2023-08-18力扣每日一题

链接: 1388. 3n 块披萨 题意: 一个长度3n的环,选n次数字,每次选完以后相邻的数字会消失,求选取结果最大值 解: 这波是~~(ctrl)CV工程师了~~ 核心思想是选取n个不相邻的元素一定…

无涯教程-Perl - splice函数

描述 此函数从LENGTH元素的OFFSET元素中删除ARRAY元素,如果指定,则用LIST替换删除的元素。如果省略LENGTH,则从OFFSET开始删除所有内容。 语法 以下是此函数的简单语法- splice ARRAY, OFFSET, LENGTH, LISTsplice ARRAY, OFFSET, LENGTHsplice ARRAY, OFFSET返回值 该函数…

Vue 项目运行 npm install 时,卡在 sill idealTree buildDeps 没有反应

解决方法:切换到淘宝镜像。 以下是之前安装的 xmzs 包,用于控制切换淘宝镜像。 该截图是之前其他项目切换淘宝镜像的截图。 切换镜像后,顺利执行 npm install 。

生成国密密钥对

在线生成国密密钥对 生成的密钥对要妥善保管,丢失是无法找回的。

selinux

一、selinux的说明 二、selinux的工作原理 三、selinux的启动、关闭与查看 Enforcing和permissive都是临时的,重启还是依据配置文件中,禁用selinux,修改配置文件: 之后重启生效 四、selinux对linux服务的影响

SpringBoot 接口调用出现乱码解决 中文乱码

SpringBoot 接口调用出现乱码解决 package com.cxjg.mvc.util;import org.springframework.context.annotation.Bean; import org.springframework.context.annotation.Configuration; import org.springframework.http.converter.HttpMessageConverter; import org.springfra…

相同数字的积木游戏

题目描述 题目描述 小华和小薇一起通过玩积木游戏学习数学。 他们有很多积木,每个积木块上都有一个数字,积木块上的数字可能相同。 小华随机拿一些积木挨着排成一排,请小薇找到这排积木中数字相同目所处位置最远的2块积木块,计算…

【JAVA】我们该如何规避代码中可能出现的错误?(一)

个人主页:【😊个人主页】 系列专栏:【❤️初识JAVA】 文章目录 前言三种类型的异常异常处理JAVA内置异常类Exception 类的层次 前言 异常是程序中的一些错误,但并不是所有的错误都是异常,并且错误有时候是可以避免的&…

【BASH】回顾与知识点梳理(三十三)

【BASH】回顾与知识点梳理 三十三 三十三. 认识系统服务 (daemons)33.1 什么是 daemon 与服务 (service)早期 System V 的 init 管理行为中 daemon 的主要分类 (Optional)systemd 使用的 unit 分类systemd 的配置文件放置目录systemd 的 unit 类型分类说明 33.2 透过 systemctl…

Grounding dino + segment anything + stable diffusion 实现图片编辑

目录 总体介绍总体流程 模块介绍目标检测: grounding dino目标分割:Segment Anything Model (SAM)整体思路模型结构:数据引擎 图片绘制 集成样例 其他问题附录 总体介绍 总体流程 本方案用到了三个步骤,按顺序依次为&#xff1a…