「最优化算法」牛顿法

牛顿法(Newton Method)

牛顿法其实是一个寻找函数的根的方法,而后运用到优化问题中。主要是运用泰勒级数的知识来寻找一个可导函数 f(x)=0f(x)=0 的根。

将函数 f(x)f(x)x0x_0 处展开成带皮亚诺余项的泰勒展开式:

f(x)=f(x0)+f(x0)(xx0)+o(x)f(x) = f(x_0) + f'(x_0)(x-x_0)+o(x)

取前两项,作为 f(x)f(x) 的近似,那么就可以用

f(x0)+f(x0)(xx0)=0f(x_0) + f'(x_0)(x-x_0)=0

的解来近似 f(x)f(x) 的解,求得的解为

x1=x0f(x0)f(x0)x_1 = x_0-\frac{f(x_0)}{f'(x_0)}

如下图:

image.png

由于函数只有一阶展开,求解处的值离真实的根还有一定的差距,但是已经很接近结果了,我们只需要再将 f(x)f(x)x1x_1 处泰勒展开,重复以上的步骤,那么就会得到比 x1x_1 更加接近的解 x2x_2

image.png

只要我们重复上述的步骤,循环迭代 nn 次,就会得到近似解 xnx_n

xn=xn1f(xn1)f(xn1)x_n =x_{n-1}-\frac{f(x_{n-1})}{f'(x_{n-1})}

利用牛顿法求解最优化问题

对于一求个无约束的最优化问题 minxf(x),xRn\min_{\mathbf x} f(\mathbf x),\mathbf x\in \mathbb R^n 的最优解,可以将其转换为求解其梯度等于 00 时的根,即 f(x)=0\nabla f(\mathbf x)=0。你看,这就可以利用牛顿法求解最优化问题啦。

采用牛顿法求解的迭代式为:

xk+1=xkH1(xk)f(xk)\mathbf x_{k+1} = \mathbf x_k -H^{-1}(\mathbf x_k) \cdot \nabla f(\mathbf x_k)

这个 H(x)H(\mathbf x) 是海塞矩阵,用来描述一个多元函数二阶导数。

海塞矩阵(Hessian Matrix)

对于一个多元函数 f(x1,x2,,xn)f(x_1,x_2,\dots,x_n),如果函数的二阶导数都存在,则定义 ff 的海塞矩阵为

H(f)=(2fx122fx1x22fx1xn2fx2x12fx222fx1xn2fxnx12fxnx22fxn2)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}

如果多元函数 f(x1,x2,,xn)f(x_1,x_2,\dots,x_n) 二阶连续可导,并且在某一点 α(x1,x2,,xn)\alpha(x_1′,x_2′,\dots,x_n’) 处梯度为 00,即f(α)=0\nabla f(\alpha)=0,那么点 α\alpha 为函数驻点,但不能判断该点是否为极小值。通过之前学习的知识,此时需要对函数求二阶导。

ff 在点 α\alpha 处的海塞矩阵为 H(α)H(\alpha),因为

2fxixj=2fxjxi\frac{\partial^2f}{\partial x_i \partial x_j}=\frac{\partial^2f}{\partial x_j \partial x_i}

那么海塞矩阵是一个对称矩阵,对于 H(α)H(\alpha),有以下结论

  • 如果 H(α)H(\alpha) 是正定矩阵,则驻点 α\alpha 为极小值点
  • 如果 H(α)H(\alpha) 是负定矩阵,则驻点 α\alpha 为极大值点
  • 如果 H(α)H(\alpha) 是不定矩阵,则驻点 α\alpha 不是极值点

因为牛顿法在进行迭代时运用了二阶导数,也就是说用二次曲线连模拟点 xk\mathbf x_k,一步便可以找到二次曲线的极值点,所以下降的速度比起梯度下降法更快,能更容易地找到答案;但由于 H1(xk)H^{-1}(\mathbf x_k) 难于求解,难于操作,于是一种牛顿法的优化方法便出现了,就是下面提到的拟牛顿法。

LM 修正牛顿法

Levenberg-Marquardt 修正牛顿法主要着手于处理 H1(xk)H^{-1}(\mathbf x_k) 为不可逆的情况,令

G(xk)=H(xk)+μIG(\mathbb x_k) = H(\mathbf x_k) + \mu I

使 μ\mu 充分大,则 G(xk)G(\mathbf x_k) 的特征值 λi+μ\lambda_i + \mu 均大于 00

μ\mu \to \infty 时, H1(xk)f(xk)H^{-1}(\mathbf x_k) \cdot \nabla f(\mathbf x_k) 趋近于 1μf(xk)\frac 1 \mu \nabla f(\mathbf x_k),那么该法就趋近于梯度下降法。

高斯-牛顿法(Gauss-Newton Method)

高斯-牛顿法是用来解非线性的最小二乘问题。设有 mm 个样本点 {(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(i)Rnx^{(i)}\in \mathbb R^n,用于训练一个可以用来预测的函数 fθ(x)f_\theta(x)

有别于其他模型,这里我们特别地记 fi(x)=f(x(i);θ)f_i(x)=f(x^{(i)};\theta) 。有损失函数:

L(θ)=i=1m(fi(θ)y(i))2L(\theta) = \sum_{i=1}^m(f_i(\theta) – y^{(i)})^2

上式对 θj\theta_j 进行求导,有

Lθj=i=1m[2(fi(θ)y(i))fiθ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]

写成向量和矩阵形式:

L(θ)=(Lθ1Lθ2Lθn)=2(f1θ1f2θ1fmθ1f1θ2f2θ2fmθ2f1θnf2θnfmθn)((f1(θ)y(1))(f2(θ)y(2))(fm(θ)y(m)))=2JTr\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}

其中 JJf(x;θ)f(x;\theta) 关于 θ\theta 的雅可比矩阵,r\mathbf r 为样本与模型的残差。

接下来我们根据已经求得的梯度,更进一步求海森矩阵的元素 hjkh_{jk}

hjk=2Lθjθk=θj(i=1m[2(fi(θ)y(i))fiθk])=2i=1m(fiθjfiθk+(fi(θ)y(i))2fiθ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(θ)=2(JTJ+S)H(\theta) = 2(J^TJ+S)

其中

sij=i=1m(fi(θ)y(i))2fiθjθks_{ij} = \sum_{i=1}^m(f_i(\theta) – y^{(i)})\frac{\partial^2 f_i}{\partial\theta_j\theta_k}

在迭代过程中,由于残差项很小,往往(?)在迭代过程中可以忽略,则我们找到了一个新的矩阵 JTJJ^TJ 来代替海塞矩阵 HH,回归到求解 x\mathbf x 的问题,得到迭代式

xk+1=xk(JTJ)1JTr\mathbf x_{k+1} = \mathbf x_k -(J^TJ)^{-1} J^T\mathbf r

拟牛顿法(Quasi-Newton Methods)

由于海塞矩阵矩阵不一定可逆,而且求矩阵逆的复杂度较高,所以引入了拟牛顿法。拟牛顿法其实就是准找可以替代矩阵 H1(xk)H^{-1}(\mathbf x_k) 的方法。

我们对 f(x)\nabla f(\mathbf x)xkx_k 点处做一阶泰勒展开:

f(x)=f(xk)+H(xk)(xxk)\nabla f(\mathbf x) = \nabla f(\mathbf x_k) + H(\mathbf x_{k})\cdot(\mathbf x-\mathbf x_k)

或者说对 f(x)f(\mathbf x)xkx_k 点处做二阶泰勒展开:

f(x)=f(xk)+(xxk)Tf(xk)+12(xxk)TH(xk)(xxk)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)

并对其求导。

x=xk+1\mathbf x = \mathbf x_{k+1} 有:

f(xk+1)f(xk)=H(xk)(xk+1xk)\nabla f(\mathbf x_{k+1})-\nabla f(\mathbf x_k) = H(\mathbf x_k)\cdot(\mathbf x_{k+1}-\mathbf x_k)

yk=f(xk+1)f(xk)\mathbf y_k = \nabla f(\mathbf x_{k+1})-\nabla f(\mathbf x_k)δk=xk+1xk\delta_k = \mathbf x_{k+1}-\mathbf x_k,有:

yk=H(xk)δk\mathbf y_k = H(\mathbf x_k) \cdot \delta_k

上式被称为拟牛顿条件。根据对矩阵 H1(xk)H^{-1}(\mathbf x_k) 或者 H(xk)H(\mathbf x_k) 估计方法的不同,拟牛顿法可以分为以下的几种:

DFP 变尺度法

DFP 方法采用矩阵 G(xk)G(\mathbf x_k) 来对 H1(xk)H^{-1}(\mathbf x_k) 进行近似,最早是由 Davidon 提出,后经 Fletcher 和 Powell 解释和改进,在命名时以三个人名字的首字母命名。

在 DFP 方法中,假设迭代式:

G(xk+1)=G(xk)+ΔG(xk)G(\mathbf x_{k+1}) = G(\mathbf x_{k}) + \Delta G(\mathbf x_{k})

其中 G(xk)G(\mathbf x_{k}) 是正定的,令

ΔGk=αukukT+βvkvkT,uk,vkRn\Delta G_k = \alpha u_ku_k^T + \beta v_kv_k^T, \quad u_k,v_k \in \mathbb R^n

代入上上式可得

G(xk+1)=G(xk)+αukukT+βvkvkTG(\mathbf x_{k+1}) = G(\mathbf x_{k}) + \alpha u_ku_k^T + \beta v_kv_k^T

已知有 H1(xk)yk=δk H^{-1}(\mathbf x_k) \cdot \mathbf y_k =\delta_k,则

δk=G(xk)yk+αukukTyk+βvkvkTyk=G(xk)yk+α(ukTyk)uk+β(vkTyk)vk\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}

可知 ukTyk,vkTyku_k^T\mathbf y_k,v_k^T\mathbf y_k 为数,我们设 uk=rδku_k = r\delta_kvk=θG(xk)ykv_k = \theta G(\mathbf x_k)\mathbf y_k,有

δk=G(xk)yk+α(rδkTyk)rδk+β(θykTG(xk)yk)θG(xk)yk\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

整理得:

(αr2(δkTyk)1)δk+(βθ2(ykTG(xk)yk)+1)G(xk)yk=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

{αr2(δkTyk)=1βθ2(ykTG(xk)yk)=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}

最终代入可得

G(xk+1)=G(xk)+δkδkTδkTykG(xk)ykykTG(xk)ykTG(xk)ykG(\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}

BFGS 方法

BFGS 算法采用矩阵 B(xk)B(\mathbf x_k) 来对 H(xk)H(\mathbf x_k) 进行近似,采取上节 DFP 方法的推导方法,可得

B(xk+1)=B(xk)+ykykTykTδkB(xk)δkδkTB(xk)δkTB(xk)δkB(\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}

可记 G(xk)=B(xk)1G(\mathbf x_{k})=B(\mathbf x_{k})^{-1},那么有

G(xk+1)=(B(xk)+ykykTykTδkB(xk)δkδkTB(xk)δkTB(xk)δk)1G(\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}

本式可以运用 Sherman-Morrison-Woodbury 公式求解

Sherman-Morrison-Woodbury 公式是一种矩阵求逆的方法,公式如下:

(A+UCV)1=A1A1UVA1C1+VA1U(A+UCV)^{-1} = A^{-1}-\frac{A^{-1}UVA^{-1}}{C^{-1}+VA^{-1}U}

还有一个向量版的,叫 Sherman-Morrison 公式:

(A+uvT)1=A1A1uvTA11+vTA1u(A+uv^{T})^{-1} = A^{-1}-\frac{A^{-1}uv^{T}A^{-1}}{1+v^{T}A^{-1}u}

经过上面两个计算的不停的折腾(完整推倒过程可以 点这里 ),求解得:

G(xk+1)=(IδkykTykTδk)G(xk)(IykδkTykTδk)+δkδkTykTδkG(\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}

ρ=1ykTδk\rho = \frac 1 {\mathbf y^T_k\delta_k},可以得出最终的 BFGS 方法的迭代方程:

G(xk+1)=(IρδkykT)G(xk)(IρykδkT)+ρδkδkTG(\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}

L-BFGS 算法

这个算法的英文全称为Limited-memory BFGS(或 Limited-storage BFGS),意思就是上面的 BFGS 方法同样需要消耗计算机的很多计算能力,于是需要一个简洁版的算法。这个 L-BFGS 算法就是这个意思。

想看的看下面这个博客

liuxiaofei.com.cn/blog/lbfgs%…

© 版权声明
THE END
喜欢就支持一下吧
点赞0

Warning: mysqli_query(): (HY000/3): Error writing file '/tmp/MYQAygYr' (Errcode: 28 - No space left on device) in /www/wwwroot/583.cn/wp-includes/class-wpdb.php on line 2345
admin的头像-五八三
评论 抢沙发
头像
欢迎您留下宝贵的见解!
提交
头像

昵称

图形验证码
取消
昵称代码图片