前言

本篇是牛顿法的最后一篇，Levenberg-Marquardt算法，也就是阻尼牛顿法。

高斯牛顿法

回顾上篇中的高斯牛顿法：
f(x0+Δx)=f(x0)+J(x0)Δxmin⁡Δx12∣∣f(x)∣∣2=min⁡Δx12(f(x0)+J(x0)Δx)T(f(x0)+J(x0)Δx)JTJΔx=−JTff({\bf x_0+\Delta x}) = f({\bf x_0}) + J({\bf x_0}){\bf \Delta x} \\ \quad \\ \min_{\bf \Delta x} \frac{1}{2} ||f({\bf x})||^2 \\ = \min_{\bf \Delta x} \frac{1}{2} (f({\bf x_0})+J({\bf x_0}){\bf \Delta x})^T(f({\bf x_0})+J({\bf x_0}){\bf \Delta x}) \\ \quad \\ J^TJ{\bf \Delta x} = -J^Tf f(x0+Δx)=f(x0)+J(x0)ΔxΔxmin21∣∣f(x)∣∣2=Δxmin21(f(x0)+J(x0)Δx)T(f(x0)+J(x0)Δx)JTJΔx=−JTf
将牛顿法中的海森矩阵HHH，使用一阶梯度JTJJ^TJJTJ进行了替换，因而避免了计算二阶梯度。然而，JTJJ^TJJTJ并非正定，可能是奇异阵或者病态矩阵，求逆不稳定，使得迭代的稳定性和收敛性差。

LM算法

阻尼牛顿法

为了提升迭代的稳定性，在优化函数上添加阻尼（机器学习里有时称为正则项），来降低迭代步长：
min⁡Δx(12∣∣f(x)∣∣2+12λΔxTΔx)=min⁡Δx12(f(x0)+J(x0)Δx)T(f(x0)+J(x0)Δx)+12λΔxTΔx)\min_{\bf \Delta x} (\frac{1}{2} ||f({\bf x})||^2+\frac{1}{2}\lambda{\bf \Delta x^T\Delta x}) \\ = \min_{\bf \Delta x} \frac{1}{2} (f({\bf x_0})+J({\bf x_0}){\bf \Delta x})^T(f({\bf x_0})+J({\bf x_0}){\bf \Delta x}) + \frac{1}{2} \lambda {\bf \Delta x^T\Delta x}) \\ Δxmin(21∣∣f(x)∣∣2+21λΔxTΔx)=Δxmin21(f(x0)+J(x0)Δx)T(f(x0)+J(x0)Δx)+21λΔxTΔx)
在12∣∣f(x)∣∣2\frac{1}{2}||f({\bf x})||^221∣∣f(x)∣∣2后添加12λΔxTΔx\frac {1}{2}\lambda {\bf \Delta x}^T {\bf \Delta x}21λΔxTΔx，当求出的Δx\bf \Delta xΔx比较大时，对函数值进行惩罚。

还是按照高斯牛顿法一样求一阶梯度等于零，获得迭代表达式：
JTJΔx+λΔx=−JTf(JTJ+λI)Δx=−JTfΔx=−(JTJ+λI)−1JTfJ^TJ{\bf \Delta x} +\lambda {\bf \Delta x} = -J^Tf \\ (J^TJ+\lambda I){\bf \Delta x}=-J^Tf \\ {\bf \Delta x} = -(J^TJ+\lambda I)^{-1}J^Tf \\ JTJΔx+λΔx=−JTf(JTJ+λI)Δx=−JTfΔx=−(JTJ+λI)−1JTf
可以看出，阻尼牛顿法通过JTJ+λIJ^TJ+\lambda IJTJ+λI替换HHH，一方面避免二阶梯度计算，另一方面保证了矩阵的可逆性。

阻尼系数λ\lambdaλ

上面的惩罚项中，λ\lambdaλ用于控制Δx\bf \Delta xΔx对迭代步长的影响。当f(x0+Δx)f(\bf x_0 +\Delta x)f(x0+Δx)与一阶泰勒展开f(x0)+J(x0)Δxf({\bf x_0}) +J({\bf x_0}){\bf\Delta x}f(x0)+J(x0)Δx足够近似，迭代步长可以大一些，加速收敛；如果函数值与泰勒展开不够近似，迭代步长要小一些，增加稳定。

因此，可以设计以下策略调整阻尼系数：
ρ=f(x0+Δx)−f(x0)J(x0)Δxifρ>34:λ=12λifρ<12:λ=2λ\rho = \frac {f({\bf x_0 + \Delta x})-f({\bf x_0})} {J({\bf x_0}){\bf \Delta x}} \\ \quad \\ if \quad \rho > \frac{3}{4}: \quad \lambda=\frac{1}{2}\lambda \\ \quad \\ if \quad \rho < \frac{1}{2}: \quad \lambda= 2\lambda ρ=J(x0)Δxf(x0+Δx)−f(x0)ifρ>43:λ=21λifρ<21:λ=2λ

置信域法

LM算法还可以通过置信域思想实现：
min⁡Δx12∣∣f(x)∣∣2s.t.ΔxTΔx≤d\min_{\bf \Delta x} \frac{1}{2} ||f({\bf x})||^2 \\ \quad \\ s.t. \quad {\bf \Delta x}^T{\bf \Delta x} \le d Δxmin21∣∣f(x)∣∣2s.t.ΔxTΔx≤d
通过限制Δx\bf \Delta xΔx的长度来约束迭代步长。这时可以通过拉格朗日乘子法把有约束的优化问题转换成无约束的优化问题：

12∣∣f(x0)+J(x0)Δx∣∣2+12λ(ΔxTΔx−d)\frac{1}{2}||f({\bf x_0})+J({\bf x_0}){\bf \Delta x}||^2+\frac{1}{2}\lambda({\bf \Delta x}^T{\bf \Delta x} - d) 21∣∣f(x0)+J(x0)Δx∣∣2+21λ(ΔxTΔx−d)
这个优化问题涉及到KKT收敛性等等判定，最后的迭代表示式形式与阻尼牛顿法相同

置信域ddd

置信域法中的λ\lambdaλ是求解拉格朗日乘子法计算出来的，而我们需要调节的是Δx\bf \Delta xΔx的置信域ddd：
ρ=f(x0+Δx)−f(x0)J(x0)Δx0ifρ>34:d=2difρ<12:d=12d\rho = \frac {f({\bf x_0 + \Delta x})-f({\bf x_0})} {J({\bf x_0}){\bf \Delta x_0}} \\ \quad \\ if \quad \rho > \frac{3}{4}: \quad d=2d \\ \quad \\ if \quad \rho < \frac{1}{2}: \quad d=\frac{1}{2}d ρ=J(x0)Δx0f(x0+Δx)−f(x0)ifρ>43:d=2difρ<21:d=21d

LM算法的流程

最后给出LM算法的整体流程：

根据迭代表达式计算Δx\bf \Delta xΔx
计算ρ\rhoρ，判断是否满足终止条件，不满足则计算x=Δx+x0\bf x=\Delta x+ x_0x=Δx+x0，否则迭代结束
调节阻尼系数λ\lambdaλ，或者调节置信域ddd，返回步骤1

代码示例

给出阻尼牛顿思想的LM算法代码：

import numpy as np
import scipy.optimize
import time
import mathdef partial_derivate_xy(x, y, F):dx = (F(x + 0.001, y) - F(x, y))/0.001dy = (F(x, y + 0.001) - F(x, y))/0.001return dx, dydef non_linear_func(x, y):fxy = 0.5 * (x ** 2 + y ** 2)return fxydef non_linear_func_2(x, y):fxy = x*x + 2*y*y + 2*x*y + 3*x - y - 2return fxydef non_linear_func_3(x, y):fxy = 0.5 * (x ** 2 - y ** 2)return fxydef non_linear_func_4(x, y):fxy = x**4 + 2*y**4 + 3*x**2*y**2 + 4*x*y**2 + x*y + x + 2*y + 0.5return fxydef non_linear_func_5(x, y):fxy = math.exp(x) + math.exp(0.5 * y) + xreturn fxydef non_linear_func_5_least_square(x, y):fxy = math.pow(math.exp(x) + math.exp(0.5 * y) + x, 2)return fxydef damping_newton(x, y, F, l):dx, dy = partial_derivate_xy(x, y, F)fx = F(x, y)grad = np.array([[dx], [dy]])H = np.matmul(grad, grad.T) + l * np.eye(2)g = - grad * fxvec_delta = np.matmul(np.linalg.inv(H), g)vec_opt = np.array([[x], [y]]) + vec_deltax_opt = vec_opt[0][0]y_opt = vec_opt[1][0]rho = (F(x_opt, y_opt) - F(x, y)) / (np.matmul(grad.T, vec_delta))if rho < 0.25:l *= 2if rho > 0.75:l *= 0.5return x_opt, y_opt, vec_delta, ldef optimizer(x0, y0, F, th=0.0001):x = x0y = y0counter = 0l = 1while True:x_opt, y_opt, vec_delta, l = damping_newton(x, y, F, l)if np.linalg.norm(vec_delta) < th:breakx = x_opty = y_optcounter = counter + 1print('iter: {}'.format(counter), 'optimized (x, y) = ({}, {})'.format(x, y), 'lambda: {}'.format(l))return x, ydef verify_min(x, y, F):fx = F(x, y)deltax = np.linspace(-0.1, 0.1, 100)deltay = np.linspace(-0.1, 0.1, 100)x_range = x + deltaxy_range = y + deltaycounter = 0for i in range(100):for j in range(100):f_range = F(x_range[i], y_range[j])f_delta = fx - f_rangeif f_delta < 0:counter += 1print('counter: {}'.format(counter))if __name__ == '__main__':x0 = 2.y0 = 2.start = time.time()for i in range(1000):result_x, result_y = optimizer(x0, y0, non_linear_func_5)if i == 0:breakend = time.time()print(result_x, result_y, 'cost time: {}'.format(end - start))print(partial_derivate_xy(result_x, result_y, non_linear_func_5_least_square))verify_min(result_x, result_y, non_linear_func_5_least_square)# scipyRes = scipy.optimize.fmin_cg(scipyF, np.array([0, 0]))# print(scipyRes)

数值计算之 Levenberg-Marquardt算法相关推荐

LM(Levenberg–Marquardt)算法原理及其python自定义实现
LM算法原理及其python自定义实现 LM(Levenberg–Marquardt)算法原理 LM算法python实现实现步骤: 代码: 运行结果: LM(Levenberg–Marquardt) ...
Levenberg–Marquardt算法学习
本次是对Levenberg–Marquardt的学习总结,是为之后看懂sparse bundle ajdustment打基础.这篇笔记包含如下内容: 回顾高斯牛顿算法,引入LM算法惩罚因子的计算(迭 ...
非线性最小二乘问题与Levenberg–Marquardt算法详解
1 最小二乘问题给定一组连续函数 f:Rn→Rm,m⩾n{\mathbf{f}}:{\mathbb{R}^n} \to {\mathbb{R}^m},{\text{ }}m \geqslant nf ...
手写系列之手写LM(Levenberg–Marquardt)算法(基于eigen)
紧接上次的手写高斯牛顿,今天顺便将LM算法也进行了手写实现,并且自己基于eigen的高斯牛顿进行了对比,可以很明显的看到,LM的算法收敛更快,精度也略胜一筹,这次高博的书不够用了,参考网上伪代码进行实 ...
Levenberg–Marquardt算法
Levenberg-Marquardt 算法最优估计算法中通常的做法都是写一个代价函数,并迭代计算最小化代价函数. 解决非线性最小二乘问题的方法有高斯牛顿法等,下面介绍Levenberg-Marqu ...
lm opencv 算法_Levenberg–Marquardt算法学习（和matlab的LM算法对比）
回顾高斯牛顿算法,引入LM算法惩罚因子的计算(迭代步子的计算) 完整的算法流程及代码样例 1. 回顾高斯牛顿,引入LM算法根据之前的博文:Gauss-Newton算法学习假设我们研究如 ...
Levenberg–Marquardt algorithm
Levenberg-Marquardt又称莱文伯格-马夸特方法(Levenberg–Marquardt algorithm)能提供数非线性最小化(局部最小)的数值解. 此算法能借由执行时修改参数达到结 ...
高斯牛顿算法matlab代码,matlab实现高斯牛顿法、Levenberg–Marquardt方法
高斯牛顿法: function [ x_ans ] = GaussNewton( xi, yi, ri) % input : x = the x vector of 3 points % y = th ...
Levenberg–Marquardt（LM）
Levenberg–Marquardt(LM)详解 1.基础概念 1.1.信赖域法 1.2.泰勒展开 1.2.正定矩阵(positive definite matrix) 1.3.雅克比矩阵(Jaco ...
如何高效的通过BP算法来训练CNN
< Neural Networks Tricks of the Trade.2nd>这本书是收录了1998-2012年在NN上面的一些技巧.原理.算法性文章,对于初学者或者是正在学习NN的 ...

数值计算之 Levenberg-Marquardt算法