【最优化】牛顿法、高斯-牛顿法

一、牛顿法

1、牛顿法在一维搜索中的应用

在一维搜索中我们所要解决的问题是如何找函数f(x)的最小值。
牛顿法的核心思想是用二次函数拟合函数f(x)的某一邻域区间，用二次函数的极小值点作为下一次的迭代点。通过多次迭代使得二次函数的极小值逼近函数f(x)的极小值
g ( x ) = f ( x ( k ) ) + f ′ ( x ( k ) ) ( x − x ( k ) ) + 1 2 f ′ ′ ( x ( k ) ) ( x − x ( k ) ) 2 g ( x ) ≈ f ( x ) , 求 f ( x ) 最小值 ≈ 求 g ( x ) 最小值 g ′ ( x ) = f ′ ( x ( k ) ) + f ′ ′ ( x ( k ) ) ( x − x ( k ) ) 令 g ′ ( x ) = 0 , x = x ( k ) − f ′ ( x ( k ) ) f ′ ′ ( x ( k ) ) 只有在 f ′ ′ ( x ) > 0 时成立， f ′ ( x ) = 0 只能保证该点为极值点， f ′ ′ ( x ) > 0 保证该点为极小值点 g(x)=f(x(k))+f′(x(k))(x−x(k))+12f″(x(k))(x−x(k))2g(x)≈f(x),求f(x)最小值≈求g(x)最小值g′(x)=f′(x(k))+f″(x(k))(x−x(k))令g′(x)=0,x=x(k)−f′(x(k))f″(x(k))只有在f″(x)>0时成立，f′(x)=0只能保证该点为极值点，f″(x)>0保证该点为极小值点

​g(x)=f(x(k))+f′(x(k))(x−x(k))+21​f′′(x(k))(x−x(k))2g(x)≈f(x),求f(x)最小值≈求g(x)最小值g′(x)=f′(x(k))+f′′(x(k))(x−x(k))令g′(x)=0,x=x(k)−f′′(x(k))f′(x(k))​只有在f′′(x)>0时成立，f′(x)=0只能保证该点为极值点，f′′(x)>0保证该点为极小值点​

2、牛顿法在多维函数中的应用

多维的情况与一维类似，如果是二维函数，拟合的是一个二次曲面，用二次曲面的最低点作为下一次的迭代点。
g ( X ) = f ( X ( k ) ) + ( X − X ( k ) ) ∇ f ′ ( X ( k ) ) + 1 2 ( X − X ( k ) ) T ∇ 2 f ( X ( k ) ) ( X − X ( k ) ) g ( X ) ≈ f ( X ) , 求 f ( X ) 最小值 ≈ 求 g ( X ) 最小值 令 ∇ g ( X ) = ∇ f ( X ( k ) ) + ∇ f ( X ( k ) ) ( X − X ( k ) ) = 0 如果 ∇ 2 f ( X ) > 0 ( 正定矩阵 ) , X = X ( k ) − [ ∇ f ( x ( k ) ) ] − 1 ∇ f ( x ( k ) ) g(X)=f(X(k))+(X−X(k))∇f′(X(k))+12(X−X(k))T∇2f(X(k))(X−X(k))g(X)≈f(X),求f(X)最小值≈求g(X)最小值令∇g(X)=∇f(X(k))+∇f(X(k))(X−X(k))=0如果∇2f(X)>0(正定矩阵),X=X(k)−[∇f(x(k))]−1∇f(x(k))

​g(X)=f(X(k))+(X−X(k))∇f′(X(k))+21​(X−X(k))T∇2f(X(k))(X−X(k))g(X)≈f(X),求f(X)最小值≈求g(X)最小值令∇g(X)=∇f(X(k))+∇f(X(k))(X−X(k))=0如果∇2f(X)>0(正定矩阵),X=X(k)−[∇f(x(k))]−1∇f(x(k))​

3、Levenberg-Marquardt修正

上述方法只有在$Hess$矩阵正定是成立，如果$Hess$矩阵不是正定的要怎么办？
$Hsee$矩阵是实对称矩阵($\frac{\partial f^2(X)}{\partial x_j\partial x_i}=\frac{\partial f^2(X)}{\partial x_i\partial x_j}$​​)，而实对称矩阵一定可以三角化
∇ 2 f ( X ( k ) ) = U T Λ U = [ λ 1 0 ⋯ 0 0 λ 2 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ λ n ] , U T U = I 如果 ∇ 2 f ( X ( k ) ) 非正定，说明 λ 1 ∼ λ n 中有若干个特征值小于 0

∇2f(X(k))=UTΛU=[λ10⋯00λ2⋯0⋮⋮⋱⋮00⋯λn],UTU=I如果∇2f(X(k))非正定，说明λ1∼λn中有若干个特征值小于0 ​​∇2f(X(k))=UTΛU= ​λ1​0⋮0​0λ2​⋮0​⋯⋯⋱⋯​00⋮λn​​ ​,UTU=I如果∇2f(X(k))非正定，说明λ1​∼λn​中有若干个特征值小于0​
用最小的特征值$\lambda min(\lambda min<0)$对$Hess$​矩阵进行修正
$$\begin{array}{r}\end{array}\begin{array}{rl}& {\nabla }^{2}f(X^{(k)})=U^T\Lambda U+(\epsilon -{\lambda }_{min})I\\ & =U^T\left[\begin{array}{cccc}{\lambda }_1& 0& \cdots & 0\\ 0& {\lambda }_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\lambda }_n\end{array}\right]U+\left[\begin{array}{cccc}\epsilon -{\lambda }_{min}& 0& \cdots & 0\\ 0& \epsilon -{\lambda }_{min}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \epsilon -{\lambda }_{min}\end{array}\right]U^{T}U\\ & =U^T\left[\begin{array}{cccc}{\lambda }_1+\epsilon -{\lambda }_{min}& 0& \cdots & 0\\ 0& {\lambda }_2+\epsilon -{\lambda }_{min}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\lambda }_n+\epsilon -{\lambda }_{min}\end{array}\right]U\\ & =U^T[\Lambda +(\epsilon -{\lambda }_{min})I]U\end{array}$$
为了不和原始的$Hess$矩阵偏差太大，$\epsilon$越小越好

​ Levenberg-Marquardt修正后，即保证了特征值都是正数，也保留原$Hess$矩阵尽可能多的信息

​ 在工程中，为了减少算法计算的复杂度，不会计算特征值特征向量，而是根据经验值手动设置一个${\mu }_k$,同时还会引入一个步长因子$\alpha$
$$X=X^{(k)}-{\alpha }^{(k)}[{\nabla }^{2}f(X^{(k)})+{\mu }_{k}I{]}^{-1}\nabla f(X^{(k)})$$
通过手动调节${\mu }_k$的值，使得$[{\nabla }^{2}f(X^{(k)})+{\mu }_{k}I]>0$

​(${\mu }_k\to 0$:趋近原牛顿法；$\mu_k \rightarrow \infty $:趋近步长很小的梯度下降法)

二、高斯-牛顿法

1、应用范围

高斯-牛顿法用于解决什么问题？

有一个函数 $y=Asin(\alpha t+\beta )$ ,其中$A、\alpha 、\beta$未知，已知一些输入输出数据$[t_1,y_1],[t_2,y_2],\cdots ,[t_n,y_n]$​

高斯-牛顿法想要解决的问题是如何根据已知数据，估计未知参数

这是一个非线性最小二乘问题，${\min }_{\hat{A},\hat{\alpha },\hat{\beta }}{\sum }_{i=1}^n(\hat{A}sin(\hat{\alpha }{t}_i+\hat{\beta })-y_i{)}^2$

2、高斯-牛顿法原理

考虑更加一般的情况：
min ⁡ ∑ i = 1 m ( r i ( X ) ) 2 令 r = [ r 1 , r 2 , ⋯   , r m ] T , 则目标函数为 f ( X ) = r ( X ) T r ( X ) , 为了使用牛顿法求解，需要计算梯度和 H e s s 矩阵 梯度 ∇ f ( X ) 的第 j 个元素为 : ( ∇ f ( X ) ) j = ∂ f ∂ x j ( X ) = 2 ∑ i = 1 m r i ( X ) ∂ r i ∂ x i ( X ) r 的 J a c o b i 矩阵为： J ( X ) = [ ∂ r 1 ∂ x 1 ( X ) ∂ r 1 ∂ x 2 ( X ) ⋯ ∂ r 1 ∂ x n ( X ) ∂ r 2 ∂ x 1 ( X ) ∂ r 2 ∂ x 2 ( X ) ⋯ ∂ r 2 ∂ x n ( X ) ⋮ ⋮ ⋱ ⋮ ∂ r m ∂ x 1 ( X ) ∂ r m ∂ x 2 ( X ) ⋯ ∂ r m ∂ x n ( X ) ] 因此，函数 f 的梯度可表示为： ∇ f ( X ) = 2 J ( X ) T r ( X ) 函数 f 的 H e s s 矩阵的第 ( k , j ) 个元素为： ∂ 2 f ∂ x k ∂ x j ( X ) = ∂ ∂ x k ( ∂ f ∂ x j ( X ) ) = ∂ ∂ x k ( 2 ∑ i = 1 m r i ( X ) ∂ r i ∂ x i ( X ) ) = 2 ∑ i = 1 m ( ∂ r i ∂ x k ( X ) ∂ r i ∂ x j ( X ) + r i ( X ) ∂ 2 r i ∂ x k ∂ x j ( X ) ) 令 S ( X ) 表示一个矩阵其中 ( k , j ) 的元素为 ∑ i = 1 m r i ( X ) ∂ 2 r i ∂ x k ∂ x j ( X ) f ( x ) 的 H e s s 矩阵可以表示为： ∇ 2 f = 2 ( J ( X ) T J ( X ) + S ( X ) ) 迭代公式为： X = X ( k ) − ( J ( X ) T J ( X ) + S ( X ) ) − 1 J ( X ) T r ( X ) 由于 S ( X ) 包含函数 r 的二阶导，数值较小可以忽略，所以迭代公式可变为： X = X ( k ) − ( J ( X ) T J ( X ) ) − 1 J ( X ) T r ( X ) \begin{aligned} &\min \sum_{i=1}^{m}(r_i(X))^2 \\ &令r=[r_1,r_2,\cdots,r_m]^T,则目标函数为f(X)=r(X)^Tr(X),为了使用牛顿法求解，需要计算梯度和Hess矩阵\\ &梯度\nabla f(X)的第j个元素为:(\nabla f(X))_j = \frac{\partial f}{\partial x_j}(X) = 2\sum_{i=1}^{m}r_i(X)\frac{\partial r_i}{\partial x_i}(X)\\ &r的Jacobi矩阵为：J(X)= \begin{bmatrix} \frac{\partial r_1}{\partial x_1}(X) & \frac{\partial r_1}{\partial x_2}(X) & \cdots &\frac{\partial r_1}{\partial x_n}(X) \\ \frac{\partial r_2}{\partial x_1}(X) & \frac{\partial r_2}{\partial x_2}(X) & \cdots &\frac{\partial r_2}{\partial x_n}(X) \\ \vdots &\vdots &\ddots &\vdots\\ \frac{\partial r_m}{\partial x_1}(X) & \frac{\partial r_m}{\partial x_2}(X) & \cdots &\frac{\partial r_m}{\partial x_n}(X) \\ \end{bmatrix}\\ &因此，函数f的梯度可表示为：\nabla f(X) = 2J(X)^Tr(X) \\ \\ &函数f的Hess矩阵的第(k,j)个元素为：\\ & \frac{\partial^2f}{\partial x_k \partial x_j}(X) =\frac{\partial}{\partial x_k}\left ( \frac{\partial f}{\partial x_j}(X) \right ) = \frac{\partial}{\partial x_k}\left ( 2\sum_{i=1}^{m}r_i(X)\frac{\partial r_i}{\partial x_i}(X)\right ) = 2\sum_{i=1}^{m}\left( \frac{\partial r_i}{\partial x_k}(X)\frac{\partial r_i}{\partial x_j}(X) +\color{blue} r_i(X)\frac{\partial^2r_i}{\partial x_k \partial x_j}(X) \color{black} \right) \\ &令\color{blue}S(X)\color{black}表示一个矩阵其中(k,j)的元素为\color{blue}\sum_{i=1}^{m} r_i(X)\frac{\partial^2r_i}{\partial x_k \partial x_j}(X) \color{black} \\ &f(x)的Hess矩阵可以表示为：\nabla^2f=2\left(J(X)^TJ(X)+\color{blue} S(X)\color{black} \right)\\ &迭代公式为：X = X^{(k)}-\left(J(X)^TJ(X)+\color{blue} S(X)\color{black} \right)^{-1}J(X)^Tr(X)\\ &由于S(X)包含函数r的二阶导，数值较小可以忽略，所以迭代公式可变为：\\ &X = X^{(k)}-\left(J(X)^TJ(X)\right)^{-1}J(X)^Tr(X)\\ \end{aligned}

​mini=1∑m​(ri​(X))2令r=[r1​,r2​,⋯,rm​]T,则目标函数为f(X)=r(X)Tr(X),为了使用牛顿法求解，需要计算梯度和Hess矩阵梯度∇f(X)的第j个元素为:(∇f(X))j​=∂xj​∂f​(X)=2i=1∑m​ri​(X)∂xi​∂ri​​(X)r的Jacobi矩阵为：J(X)= ​∂x1​∂r1​​(X)∂x1​∂r2​​(X)⋮∂x1​∂rm​​(X)​∂x2​∂r1​​(X)∂x2​∂r2​​(X)⋮∂x2​∂rm​​(X)​⋯⋯⋱⋯​∂xn​∂r1​​(X)∂xn​∂r2​​(X)⋮∂xn​∂rm​​(X)​ ​因此，函数f的梯度可表示为：∇f(X)=2J(X)Tr(X)函数f的Hess矩阵的第(k,j)个元素为：∂xk​∂xj​∂2f​(X)=∂xk​∂​(∂xj​∂f​(X))=∂xk​∂​(2i=1∑m​ri​(X)∂xi​∂ri​​(X))=2i=1∑m​(∂xk​∂ri​​(X)∂xj​∂ri​​(X)+ri​(X)∂xk​∂xj​∂2ri​​(X))令S(X)表示一个矩阵其中(k,j)的元素为i=1∑m​ri​(X)∂xk​∂xj​∂2ri​​(X)f(x)的Hess矩阵可以表示为：∇2f=2(J(X)TJ(X)+S(X))迭代公式为：X=X(k)−(J(X)TJ(X)+S(X))−1J(X)Tr(X)由于S(X)包含函数r的二阶导，数值较小可以忽略，所以迭代公式可变为：X=X(k)−(J(X)TJ(X))−1J(X)Tr(X)​

3、Levenberg-Marquardt修正

高斯-牛顿法在使用的过程中也会出现$J(X{)}^{T}J(X)$不是正定矩阵的问题，所以也可以使用Levenberg-Marquardt修正就加以解决:
$X=X^{(k)}-{\left(J(X{)}^{T}J(X)+{\mu }_{k}I\right)}^{-1}J(X{)}^{T}r(X)$