牛顿法(Newton Method)
牛顿法其实是一个寻找函数的根的方法,而后运用到优化问题中。主要是运用泰勒级数的知识来寻找一个可导函数 f ( x ) = 0 f(x)=0 f ( x ) = 0 的根。
将函数 f ( x ) f(x) f ( x ) 在 x 0 x_0 x 0 处展开成带皮亚诺余项的泰勒展开式:
f ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + o ( x ) f(x) = f(x_0) + f'(x_0)(x-x_0)+o(x) f ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + o ( x )
取前两项,作为 f ( x ) f(x) f ( x ) 的近似,那么就可以用
f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) = 0 f(x_0) + f'(x_0)(x-x_0)=0 f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) = 0
的解来近似 f ( x ) f(x) f ( x ) 的解,求得的解为
x 1 = x 0 − f ( x 0 ) f ′ ( x 0 ) x_1 = x_0-\frac{f(x_0)}{f'(x_0)} x 1 = x 0 − f ′ ( x 0 ) f ( x 0 )
如下图:
由于函数只有一阶展开,求解处的值离真实的根还有一定的差距,但是已经很接近结果了,我们只需要再将 f ( x ) f(x) f ( x ) 在 x 1 x_1 x 1 处泰勒展开,重复以上的步骤,那么就会得到比 x 1 x_1 x 1 更加接近的解 x 2 x_2 x 2 。
只要我们重复上述的步骤,循环迭代 n n n 次,就会得到近似解 x n x_n x n
x n = x n − 1 − f ( x n − 1 ) f ′ ( x n − 1 ) x_n =x_{n-1}-\frac{f(x_{n-1})}{f'(x_{n-1})} x n = x n − 1 − f ′ ( x n − 1 ) f ( x n − 1 )
利用牛顿法求解最优化问题
对于一求个无约束的最优化问题 min x f ( x ) , x ∈ R n \min_{\mathbf x} f(\mathbf x),\mathbf x\in \mathbb R^n min x f ( x ) , x ∈ R n 的最优解,可以将其转换为求解其梯度等于 0 0 0 时的根,即 ∇ f ( x ) = 0 \nabla f(\mathbf x)=0 ∇ f ( x ) = 0 。你看,这就可以利用牛顿法求解最优化问题啦。
采用牛顿法求解的迭代式为:
x k + 1 = x k − H − 1 ( x k ) ⋅ ∇ f ( x k ) \mathbf x_{k+1} = \mathbf x_k -H^{-1}(\mathbf x_k) \cdot \nabla f(\mathbf x_k) x k + 1 = x k − H − 1 ( x k ) ⋅ ∇ f ( x k )
这个 H ( x ) H(\mathbf x) H ( x ) 是海塞矩阵,用来描述一个多元函数二阶导数。
海塞矩阵(Hessian Matrix)
对于一个多元函数 f ( x 1 , x 2 , … , x n ) f(x_1,x_2,\dots,x_n) f ( x 1 , x 2 , … , x n ) ,如果函数的 二阶 导数都存在,则定义 f f f 的海塞矩阵为
H ( f ) = ( ∂ 2 f ∂ x 1 2 ∂ 2 f ∂ x 1 ∂ x 2 ⋯ ∂ 2 f ∂ x 1 ∂ x n ∂ 2 f ∂ x 2 ∂ x 1 ∂ 2 f ∂ x 2 2 ⋯ ∂ 2 f ∂ x 1 ∂ x n ⋮ ⋮ ⋱ ⋮ ∂ 2 f ∂ x n ∂ x 1 ∂ 2 f ∂ x n ∂ x 2 ⋯ ∂ 2 f ∂ x n 2 ) H(f)=
\begin{pmatrix}
\frac{\partial^2f}{\partial x_1^2} & \frac{\partial^2f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2f}{\partial x_1 \partial x_n} \\
\frac{\partial^2f}{\partial x_2 \partial x_1} & \frac{\partial^2f}{\partial x_2^2} & \cdots & \frac{\partial^2f}{\partial x_1 \partial x_n} \\
\vdots & \vdots & \ddots & \vdots
\\
\frac{\partial^2f}{\partial x_n \partial x_1} & \frac{\partial^2f}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2f}{\partial x_n^2} \\
\end{pmatrix} H ( f ) = ⎝ ⎛ ∂ x 1 2 ∂ 2 f ∂ x 2 ∂ x 1 ∂ 2 f ⋮ ∂ x n ∂ x 1 ∂ 2 f ∂ x 1 ∂ x 2 ∂ 2 f ∂ x 2 2 ∂ 2 f ⋮ ∂ x n ∂ x 2 ∂ 2 f ⋯ ⋯ ⋱ ⋯ ∂ x 1 ∂ x n ∂ 2 f ∂ x 1 ∂ x n ∂ 2 f ⋮ ∂ x n 2 ∂ 2 f ⎠ ⎞
如果多元函数 f ( x 1 , x 2 , … , x n ) f(x_1,x_2,\dots,x_n) f ( x 1 , x 2 , … , x n ) 二阶连续可导,并且在某一点 α ( x 1 ′ , x 2 ′ , … , x n ′ ) \alpha(x_1′,x_2′,\dots,x_n’) α ( x 1 ′ , x 2 ′ , … , x n ′ ) 处梯度为 0 0 0 ,即∇ f ( α ) = 0 \nabla f(\alpha)=0 ∇ f ( α ) = 0 ,那么点 α \alpha α 为函数驻点,但不能判断该点是否为极小值。通过之前学习的知识,此时需要对函数求二阶导。
记 f f f 在点 α \alpha α 处的海塞矩阵为 H ( α ) H(\alpha) H ( α ) ,因为
∂ 2 f ∂ x i ∂ x j = ∂ 2 f ∂ x j ∂ x i \frac{\partial^2f}{\partial x_i \partial x_j}=\frac{\partial^2f}{\partial x_j \partial x_i} ∂ x i ∂ x j ∂ 2 f = ∂ x j ∂ x i ∂ 2 f
那么海塞矩阵是一个对称矩阵 ,对于 H ( α ) H(\alpha) H ( α ) ,有以下结论
如果 H ( α ) H(\alpha) H ( α ) 是正定矩阵,则驻点 α \alpha α 为极小值点
如果 H ( α ) H(\alpha) H ( α ) 是负定矩阵,则驻点 α \alpha α 为极大值点
如果 H ( α ) H(\alpha) H ( α ) 是不定矩阵,则驻点 α \alpha α 不是极值点
因为牛顿法在进行迭代时运用了二阶导数,也就是说用二次曲线连模拟点 x k \mathbf x_k x k ,一步便可以找到二次曲线的极值点,所以下降的速度比起梯度下降法更快,能更容易地找到答案;但由于 H − 1 ( x k ) H^{-1}(\mathbf x_k) H − 1 ( x k ) 难于求解,难于操作,于是一种牛顿法的优化方法便出现了,就是下面提到的拟牛顿法。
LM 修正牛顿法
Levenberg-Marquardt 修正牛顿法主要着手于处理 H − 1 ( x k ) H^{-1}(\mathbf x_k) H − 1 ( x k ) 为不可逆的情况,令
G ( x k ) = H ( x k ) + μ I G(\mathbb x_k) = H(\mathbf x_k) + \mu I G ( x k ) = H ( x k ) + μ I
使 μ \mu μ 充分大,则 G ( x k ) G(\mathbf x_k) G ( x k ) 的特征值 λ i + μ \lambda_i + \mu λ i + μ 均大于 0 0 0 。
当 μ → ∞ \mu \to \infty μ → ∞ 时, H − 1 ( x k ) ⋅ ∇ f ( x k ) H^{-1}(\mathbf x_k) \cdot \nabla f(\mathbf x_k) H − 1 ( x k ) ⋅ ∇ f ( x k ) 趋近于 1 μ ∇ f ( x k ) \frac 1 \mu \nabla f(\mathbf x_k) μ 1 ∇ f ( x k ) ,那么该法就趋近于梯度下降法。
高斯-牛顿法(Gauss-Newton Method)
高斯-牛顿法是用来解非线性的最小二乘问题。设有 m m m 个样本点 { ( x ( 1 ) , y ( 1 ) ) , ( x ( 2 ) , y ( 2 ) ) , ⋯ , ( x ( m ) , y ( m ) ) } \{ (x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}),\cdots,(x^{(m)},y^{(m)}) \} {( x ( 1 ) , y ( 1 ) ) , ( x ( 2 ) , y ( 2 ) ) , ⋯ , ( x ( m ) , y ( m ) )} ,x ( i ) ∈ R n x^{(i)}\in \mathbb R^n x ( i ) ∈ R n ,用于训练一个可以用来预测的函数 f θ ( x ) f_\theta(x) f θ ( x ) 。
有别于其他模型,这里我们特别地记 f i ( x ) = f ( x ( i ) ; θ ) f_i(x)=f(x^{(i)};\theta) f i ( x ) = f ( x ( i ) ; θ ) 。有损失函数:
L ( θ ) = ∑ i = 1 m ( f i ( θ ) − y ( i ) ) 2 L(\theta) = \sum_{i=1}^m(f_i(\theta) – y^{(i)})^2 L ( θ ) = i = 1 ∑ m ( f i ( θ ) − y ( i ) ) 2
上式对 θ j \theta_j θ j 进行求导,有
∂ L ∂ θ j = ∑ i = 1 m [ 2 ( f i ( θ ) − y ( i ) ) ∂ f i ∂ θ j ] \frac{\partial L}{\partial \theta_j} = \sum_{i=1}^m\left[2(f_i(\theta) – y^{(i)})\frac{\partial f_i}{\partial\theta_j}\right] ∂ θ j ∂ L = i = 1 ∑ m [ 2 ( f i ( θ ) − y ( i ) ) ∂ θ j ∂ f i ]
写成向量和矩阵形式:
∇ L ( θ ) = ( ∂ L ∂ θ 1 ∂ L ∂ θ 2 ⋮ ∂ L ∂ θ n ) = 2 ( ∂ f 1 ∂ θ 1 ∂ f 2 ∂ θ 1 ⋯ ∂ f m ∂ θ 1 ∂ f 1 ∂ θ 2 ∂ f 2 ∂ θ 2 ⋯ ∂ f m ∂ θ 2 ⋮ ⋮ ⋱ ⋮ ∂ f 1 ∂ θ n ∂ f 2 ∂ θ n ⋯ ∂ f m ∂ θ n ) ( ( f 1 ( θ ) − y ( 1 ) ) ( f 2 ( θ ) − y ( 2 ) ) ⋮ ( f m ( θ ) − y ( m ) ) ) = 2 J T r \begin{align}
\nabla L(\theta) &=
\begin{pmatrix}
\frac{\partial L}{\partial \theta_1}\\
\frac{\partial L}{\partial \theta_2}\\
\vdots\\
\frac{\partial L}{\partial \theta_n}
\end{pmatrix} =2
\begin{pmatrix}
\frac{\partial f_1}{\partial \theta_1} & \frac{\partial f_2}{\partial \theta_1} & \cdots&\frac{\partial f_m}{\partial \theta_1}\\
\frac{\partial f_1}{\partial \theta_2} & \frac{\partial f_2}{\partial \theta_2} & \cdots&\frac{\partial f_m}{\partial \theta_2}\\
\vdots & \vdots & \ddots & \vdots\\
\frac{\partial f_1}{\partial \theta_n} & \frac{\partial f_2}{\partial \theta_n} & \cdots&\frac{\partial f_m}{\partial \theta_n}
\end{pmatrix}
\begin{pmatrix}
(f_1(\theta) – y^{(1)})\\
(f_2(\theta) – y^{(2)})\\
\vdots\\
(f_m(\theta) – y^{(m)})
\end{pmatrix} \\
&=2J^Tr
\end{align} ∇ L ( θ ) = ⎝ ⎛ ∂ θ 1 ∂ L ∂ θ 2 ∂ L ⋮ ∂ θ n ∂ L ⎠ ⎞ = 2 ⎝ ⎛ ∂ θ 1 ∂ f 1 ∂ θ 2 ∂ f 1 ⋮ ∂ θ n ∂ f 1 ∂ θ 1 ∂ f 2 ∂ θ 2 ∂ f 2 ⋮ ∂ θ n ∂ f 2 ⋯ ⋯ ⋱ ⋯ ∂ θ 1 ∂ f m ∂ θ 2 ∂ f m ⋮ ∂ θ n ∂ f m ⎠ ⎞ ⎝ ⎛ ( f 1 ( θ ) − y ( 1 ) ) ( f 2 ( θ ) − y ( 2 ) ) ⋮ ( f m ( θ ) − y ( m ) ) ⎠ ⎞ = 2 J T r
其中 J J J 为 f ( x ; θ ) f(x;\theta) f ( x ; θ ) 关于 θ \theta θ 的雅可比矩阵,r \mathbf r r 为样本与模型的残差。
接下来我们根据已经求得的梯度,更进一步求海森矩阵的元素 h j k h_{jk} h jk
h j k = ∂ 2 L ∂ θ j θ k = ∂ ∂ θ j ( ∑ i = 1 m [ 2 ( f i ( θ ) − y ( i ) ) ∂ f i ∂ θ k ] ) = 2 ∑ i = 1 m ( ∂ f i ∂ θ j ∂ f i ∂ θ k + ( f i ( θ ) − y ( i ) ) ∂ 2 f i ∂ θ j θ k ) \begin{align}
h_{jk} &= \frac{\partial^2 L}{\partial\theta_j\theta_k}\\
&=\frac{\partial}{\partial\theta_j}\left(\sum_{i=1}^m\left[2(f_i(\theta) – y^{(i)})\frac{\partial f_i}{\partial\theta_k}\right]\right)\\
&= 2\sum_{i=1}^m \left( \frac{\partial f_i}{\partial\theta_j}\frac{\partial f_i}{\partial\theta_k} + (f_i(\theta) – y^{(i)})\frac{\partial^2 f_i}{\partial\theta_j\theta_k} \right)
\end{align} h jk = ∂ θ j θ k ∂ 2 L = ∂ θ j ∂ ( i = 1 ∑ m [ 2 ( f i ( θ ) − y ( i ) ) ∂ θ k ∂ f i ] ) = 2 i = 1 ∑ m ( ∂ θ j ∂ f i ∂ θ k ∂ f i + ( f i ( θ ) − y ( i ) ) ∂ θ j θ k ∂ 2 f i )
同样的,写成矩阵形式
H ( θ ) = 2 ( J T J + S ) H(\theta) = 2(J^TJ+S) H ( θ ) = 2 ( J T J + S )
其中
s i j = ∑ i = 1 m ( f i ( θ ) − y ( i ) ) ∂ 2 f i ∂ θ j θ k s_{ij} = \sum_{i=1}^m(f_i(\theta) – y^{(i)})\frac{\partial^2 f_i}{\partial\theta_j\theta_k} s ij = i = 1 ∑ m ( f i ( θ ) − y ( i ) ) ∂ θ j θ k ∂ 2 f i
在迭代过程中,由于残差项很小,往往(?)在迭代过程中可以忽略,则我们找到了一个新的矩阵 J T J J^TJ J T J 来代替海塞矩阵 H H H ,回归到求解 x \mathbf x x 的问题,得到迭代式
x k + 1 = x k − ( J T J ) − 1 J T r \mathbf x_{k+1} = \mathbf x_k -(J^TJ)^{-1} J^T\mathbf r x k + 1 = x k − ( J T J ) − 1 J T r
拟牛顿法(Quasi-Newton Methods)
由于海塞矩阵矩阵不一定可逆,而且求矩阵逆的复杂度较高,所以引入了拟牛顿法。拟牛顿法其实就是准找可以替代矩阵 H − 1 ( x k ) H^{-1}(\mathbf x_k) H − 1 ( x k ) 的方法。
我们对 ∇ f ( x ) \nabla f(\mathbf x) ∇ f ( x ) 在 x k x_k x k 点处做一阶泰勒展开:
∇ f ( x ) = ∇ f ( x k ) + H ( x k ) ⋅ ( x − x k ) \nabla f(\mathbf x) = \nabla f(\mathbf x_k) + H(\mathbf x_{k})\cdot(\mathbf x-\mathbf x_k) ∇ f ( x ) = ∇ f ( x k ) + H ( x k ) ⋅ ( x − x k )
或者说对 f ( x ) f(\mathbf x) f ( x ) 在 x k x_k x k 点处做二阶泰勒展开:
f ( x ) = f ( x k ) + ( x − x k ) T ∇ f ( x k ) + 1 2 ( x − x k ) T H ( x k ) ( x − x k ) f(\mathbf x) = f(\mathbf x_k) + (\mathbf x-\mathbf x_k)^T\nabla f(\mathbf x_k) + \frac 1 2 (\mathbf x-\mathbf x_k)^TH(\mathbf x_{k})(\mathbf x-\mathbf x_k) f ( x ) = f ( x k ) + ( x − x k ) T ∇ f ( x k ) + 2 1 ( x − x k ) T H ( x k ) ( x − x k )
并对其求导。
取 x = x k + 1 \mathbf x = \mathbf x_{k+1} x = x k + 1 有:
∇ f ( x k + 1 ) − ∇ f ( x k ) = H ( x k ) ⋅ ( x k + 1 − x k ) \nabla f(\mathbf x_{k+1})-\nabla f(\mathbf x_k) = H(\mathbf x_k)\cdot(\mathbf x_{k+1}-\mathbf x_k) ∇ f ( x k + 1 ) − ∇ f ( x k ) = H ( x k ) ⋅ ( x k + 1 − x k )
记 y k = ∇ f ( x k + 1 ) − ∇ f ( x k ) \mathbf y_k = \nabla f(\mathbf x_{k+1})-\nabla f(\mathbf x_k) y k = ∇ f ( x k + 1 ) − ∇ f ( x k ) ,δ k = x k + 1 − x k \delta_k = \mathbf x_{k+1}-\mathbf x_k δ k = x k + 1 − x k ,有:
y k = H ( x k ) ⋅ δ k \mathbf y_k = H(\mathbf x_k) \cdot \delta_k y k = H ( x k ) ⋅ δ k
上式被称为拟牛顿条件 。根据对矩阵 H − 1 ( x k ) H^{-1}(\mathbf x_k) H − 1 ( x k ) 或者 H ( x k ) H(\mathbf x_k) H ( x k ) 估计方法的不同,拟牛顿法可以分为以下的几种:
DFP 变尺度法
DFP 方法采用矩阵 G ( x k ) G(\mathbf x_k) G ( x k ) 来对 H − 1 ( x k ) H^{-1}(\mathbf x_k) H − 1 ( x k ) 进行近似,最早是由 Davidon 提出,后经 Fletcher 和 Powell 解释和改进,在命名时以三个人名字的首字母命名。
在 DFP 方法中,假设迭代式:
G ( x k + 1 ) = G ( x k ) + Δ G ( x k ) G(\mathbf x_{k+1}) = G(\mathbf x_{k}) + \Delta G(\mathbf x_{k}) G ( x k + 1 ) = G ( x k ) + Δ G ( x k )
其中 G ( x k ) G(\mathbf x_{k}) G ( x k ) 是正定的,令
Δ G k = α u k u k T + β v k v k T , u k , v k ∈ R n \Delta G_k = \alpha u_ku_k^T + \beta v_kv_k^T, \quad u_k,v_k \in \mathbb R^n Δ G k = α u k u k T + β v k v k T , u k , v k ∈ R n
代入上上式可得
G ( x k + 1 ) = G ( x k ) + α u k u k T + β v k v k T G(\mathbf x_{k+1}) = G(\mathbf x_{k}) + \alpha u_ku_k^T + \beta v_kv_k^T G ( x k + 1 ) = G ( x k ) + α u k u k T + β v k v k T
已知有 H − 1 ( x k ) ⋅ y k = δ k H^{-1}(\mathbf x_k) \cdot \mathbf y_k =\delta_k H − 1 ( x k ) ⋅ y k = δ k ,则
δ k = G ( x k ) y k + α u k u k T y k + β v k v k T y k = G ( x k ) y k + α ( u k T y k ) u k + β ( v k T y k ) v k \begin{align}
\delta_k =& G(\mathbf x_{k})\mathbf y_k + \alpha u_ku_k^T\mathbf y_k + \beta v_kv_k^T\mathbf y_k\\
=& G(\mathbf x_{k})\mathbf y_k + \alpha (u_k^T\mathbf y_k)u_k + \beta (v_k^T\mathbf y_k)v_k
\end{align} δ k = = G ( x k ) y k + α u k u k T y k + β v k v k T y k G ( x k ) y k + α ( u k T y k ) u k + β ( v k T y k ) v k
可知 u k T y k , v k T y k u_k^T\mathbf y_k,v_k^T\mathbf y_k u k T y k , v k T y k 为数,我们设 u k = r δ k u_k = r\delta_k u k = r δ k , v k = θ G ( x k ) y k v_k = \theta G(\mathbf x_k)\mathbf y_k v k = θG ( x k ) y k ,有
δ k = G ( x k ) y k + α ( r δ k T y k ) r δ k + β ( θ y k T G ( x k ) y k ) θ G ( x k ) y k \delta_k = G(\mathbf x_{k})\mathbf y_k+\alpha( r\delta_k^T\mathbf y_k) r\delta_k + \beta (\theta \mathbf y_k^TG(\mathbf x_k)\mathbf y_k)\theta G(\mathbf x_k)\mathbf y_k δ k = G ( x k ) y k + α ( r δ k T y k ) r δ k + β ( θ y k T G ( x k ) y k ) θG ( x k ) y k
整理得:
( α r 2 ( δ k T y k ) − 1 ) δ k + ( β θ 2 ( y k T G ( x k ) y k ) + 1 ) G ( x k ) y k = 0 (\alpha r^2(\delta_k^T\mathbf y_k)-1)\delta_k+(\beta\theta^2(\mathbf y_k^TG(\mathbf x_k)\mathbf y_k)+1) G(\mathbf x_k)\mathbf y_k=0 ( α r 2 ( δ k T y k ) − 1 ) δ k + ( β θ 2 ( y k T G ( x k ) y k ) + 1 ) G ( x k ) y k = 0
令
{ α r 2 ( δ k T y k ) = 1 β θ 2 ( y k T G ( x k ) y k ) = − 1 \begin{cases}
\alpha r^2(\delta_k^T\mathbf y_k) = 1\\
\beta\theta^2(\mathbf y_k^TG(\mathbf x_k)\mathbf y_k)=-1
\end{cases} { α r 2 ( δ k T y k ) = 1 β θ 2 ( y k T G ( x k ) y k ) = − 1
最终代入可得
G ( x k + 1 ) = G ( x k ) + δ k δ k T δ k T y k − G ( x k ) y k y k T G ( x k ) y k T G ( x k ) y k G(\mathbf x_{k+1}) = G(\mathbf x_k) + \frac{\delta_k\delta_k^T}{\delta_k^T\mathbf y_k} – \frac{G(\mathbf x_k)\mathbf y_k\mathbf y_k^TG(\mathbf x_k)}{\mathbf y_k^TG(\mathbf x_k)\mathbf y_k} G ( x k + 1 ) = G ( x k ) + δ k T y k δ k δ k T − y k T G ( x k ) y k G ( x k ) y k y k T G ( x k )
BFGS 方法
BFGS 算法采用矩阵 B ( x k ) B(\mathbf x_k) B ( x k ) 来对 H ( x k ) H(\mathbf x_k) H ( x k ) 进行近似,采取上节 DFP 方法的推导方法,可得
B ( x k + 1 ) = B ( x k ) + y k y k T y k T δ k − B ( x k ) δ k δ k T B ( x k ) δ k T B ( x k ) δ k B(\mathbf x_{k+1}) = B(\mathbf x_k) + \frac{\mathbf y_k\mathbf y_k^T}{\mathbf y_k^T\delta_k} – \frac{B(\mathbf x_k)\delta_k \delta_k^TB(\mathbf x_k)}{\delta_k^TB(\mathbf x_k)\delta_k} B ( x k + 1 ) = B ( x k ) + y k T δ k y k y k T − δ k T B ( x k ) δ k B ( x k ) δ k δ k T B ( x k )
可记 G ( x k ) = B ( x k ) − 1 G(\mathbf x_{k})=B(\mathbf x_{k})^{-1} G ( x k ) = B ( x k ) − 1 ,那么有
G ( x k + 1 ) = ( B ( x k ) + y k y k T y k T δ k − B ( x k ) δ k δ k T B ( x k ) δ k T B ( x k ) δ k ) − 1 G(\mathbf x_{k+1})=\left(B(\mathbf x_k) + \frac{\mathbf y_k\mathbf y_k^T}{\mathbf y_k^T\delta_k} – \frac{B(\mathbf x_k)\delta_k \delta_k^TB(\mathbf x_k)}{\delta_k^TB(\mathbf x_k)\delta_k}\right)^{-1} G ( x k + 1 ) = ( B ( x k ) + y k T δ k y k y k T − δ k T B ( x k ) δ k B ( x k ) δ k δ k T B ( x k ) ) − 1
本式可以运用 Sherman-Morrison-Woodbury 公式 求解
Sherman-Morrison-Woodbury 公式是一种矩阵求逆的方法,公式如下:
( A + U C V ) − 1 = A − 1 − A − 1 U V A − 1 C − 1 + V A − 1 U (A+UCV)^{-1} = A^{-1}-\frac{A^{-1}UVA^{-1}}{C^{-1}+VA^{-1}U} ( A + U C V ) − 1 = A − 1 − C − 1 + V A − 1 U A − 1 U V A − 1
还有一个向量版的,叫 Sherman-Morrison 公式:
( A + u v T ) − 1 = A − 1 − A − 1 u v T A − 1 1 + v T A − 1 u (A+uv^{T})^{-1} = A^{-1}-\frac{A^{-1}uv^{T}A^{-1}}{1+v^{T}A^{-1}u} ( A + u v T ) − 1 = A − 1 − 1 + v T A − 1 u A − 1 u v T A − 1
经过上面两个计算的不停的折腾(完整推倒过程可以 点这里 ),求解得:
G ( x k + 1 ) = ( I − δ k y k T y k T δ k ) G ( x k ) ( I − y k δ k T y k T δ k ) + δ k δ k T y k T δ k G(\mathbf x_{k+1}) = (I-\frac{\delta_k\mathbf y_k^T}{\mathbf y^T_k\delta_k}) G(\mathbf x_{k}) (I-\frac{\mathbf y_k\delta_k^T}{\mathbf y^T_k\delta_k})+\frac{\delta_k\delta_k^T}{\mathbf y^T_k\delta_k} G ( x k + 1 ) = ( I − y k T δ k δ k y k T ) G ( x k ) ( I − y k T δ k y k δ k T ) + y k T δ k δ k δ k T
令 ρ = 1 y k T δ k \rho = \frac 1 {\mathbf y^T_k\delta_k} ρ = y k T δ k 1 ,可以得出最终的 BFGS 方法的迭代方程:
G ( x k + 1 ) = ( I − ρ δ k y k T ) G ( x k ) ( I − ρ y k δ k T ) + ρ δ k δ k T G(\mathbf x_{k+1}) = (I-\rho{\delta_k\mathbf y_k^T}) G(\mathbf x_{k}) (I-\rho{\mathbf y_k\delta_k^T})+\rho{\delta_k\delta_k^T} G ( x k + 1 ) = ( I − ρ δ k y k T ) G ( x k ) ( I − ρ y k δ k T ) + ρ δ k δ k T
L-BFGS 算法
这个算法的英文全称为Limited-memory BFGS(或 Limited-storage BFGS),意思就是上面的 BFGS 方法同样需要消耗计算机的很多计算能力,于是需要一个简洁版的算法。这个 L-BFGS 算法就是这个意思。
想看的看下面这个博客
liuxiaofei.com.cn/blog/lbfgs%…