• 2019-04-18 22:27:08

# 关于LM写的不错的文章：

https://blog.csdn.net/a6333230/article/details/83304098

# 三维重构PMP论文中使用的信赖域算法（LM）matlab版,修改在标定中可用

function [G, val] = rect_rot(G, A, d)
% Rectify rotations
%
% G: Rectifying transform
% A, d: Coefficients for orthonormality constraint
% val: Value of residual error

ik = size(G, 1);
ikk = numel(G);
k = ikk/ik;
Id = eye(ik);
IdG = eye(ikk);

% Parameters for trust-region optimization
delta = 1;
eta1 = 1/4;
eta2 = 3/4;
kappa = 2;
epsilon = 1e-3;

% Trust-region strategy
f = A*reshape(G*G', [], 1)-d;
c = norm(f)^2;
pc = inf;
while pc-c > 1e-10
pG = G;
pf = f;
pc = c;

J = A*kron(G, Id)+reshape(permute(reshape(A*kron(Id, G), [], k, ik), [1 3 2]), [], ikk);
H = J'*J;
g = J'*f;

while true
temp = H+IdG*delta;
r = (temp+trace(temp)*1e-10*IdG)\g;
G(:) = G(:)-r;

f = A*reshape(G*G', [], 1)-d;
c = norm(f)^2;

rho = (pc-c)/(pc-norm(pf-J*r)^2);
rho(isnan(rho)) = 0;
if rho < eta1 || pc < c
delta = delta*kappa;
elseif rho > eta2
delta = delta/kappa;
end

if rho < epsilon || pc < c
G = pG;
if delta > 1e10
c = pc;
break;
end
else
break;
end
end
end

val = mse(f);

end


更多相关内容
• 列文伯格-马夸尔特算法应用于神经网络网络参数优化
• 实现较为粗略的列文伯格-马夸尔特方法代码实现，因为matlab代码简单，很容易修改参数和进行变种。配合https://blog.csdn.net/Pengcode/article/details/107783834进行LM研究将会更加合理。
• #include "LM.h" #include <opencv2/opencv.hpp> int main() { const double a = 1, b = 5, c = 2; double a_ = 0.9, b_ = 5.5, c_ = 2.2;... LM::LevenbergMarquardt lm(a_, b_, c_);... lm.SetParameters(1e-10...
#include "LM.h"
#include <opencv2/opencv.hpp>

int main()
{
const double a = 1, b = 5, c = 2;
double a_ = 0.9, b_ = 5.5, c_ = 2.2;

LM::LevenbergMarquardt lm(a_, b_, c_);
lm.SetParameters(1e-10, 1e-10, 100, true);

const size_t N = 1000;
cv::RNG rng(cv::getTickCount());

for (size_t i = 0; i < N; i++)
{
double x = rng.uniform(0.0, 10.0);
double y = exp(a * x * x + b * x + c) + rng.gaussian(5.0);
}

lm.Slove();
return 0;
}
#pragma once
#include <iostream>
#include <vector>
#include <chrono>
#include <Eigen/Core>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

namespace LM
{
class RunTime
{
public:
void Stat();
void Stop();
double Duration();
private:
};

class LevenbergMarquardt
{
public:
LevenbergMarquardt(double da, double db, double dc);
void SetParameters(double sEpsilon1, double sEpsilin2, int sMaxIter, bool sIsout);
void AddObservation(const double x, const double y);
void CalcJ_fx();
void CalcH_g();
double GetCost();
double F(double a, double b, double c);
double L0_L(Eigen::Vector3d h);
void Slove();
public:
double epsilon1, epsilon2;
int MaxIter;
std::vector<Eigen::Vector2d> obsPt;
double a, b, c;
bool isOut;
private:
MatrixXd fx_;
MatrixXd J_;
Matrix3d H_;
Vector3d g_;
};
}


#include "LM.h"
#include <iomanip>

void LM::RunTime::Stat()
{
}

void LM::RunTime::Stop()
{
}

double LM::RunTime::Duration()
{
return chrono::duration_cast<chrono::duration<double>>(tend - tstat).count();
}

LM::LevenbergMarquardt::LevenbergMarquardt(double da, double db, double dc) :a(da), b(db), c(dc)
{
epsilon1 = 1e-6;
epsilon2 = 1e-6;
MaxIter = 50;
isOut = true;
}

void LM::LevenbergMarquardt::SetParameters(double sEpsilon1, double sEpsilon2, int sMaxIter, bool sIsout)
{
epsilon1 = sEpsilon1;
epsilon2 = sEpsilon2;
MaxIter = sMaxIter;
isOut = sIsout;
}

void LM::LevenbergMarquardt::AddObservation(const double x, const double y)
{
obsPt.push_back(Eigen::Vector2d(x, y));
}

void LM::LevenbergMarquardt::CalcJ_fx()
{
J_.resize(obsPt.size(), 3);
fx_.resize(obsPt.size(), 1);

for (int i = 0; i < obsPt.size(); i++)
{
Eigen::Vector2d obs = obsPt.at(i);
const double x = obs(0);
const double y = obs(1);
double dfda = -x * x * exp(a * x * x + b * x + c);
double dfdb = -x * exp(a * x * x + b * x + c);
double dfdc = -exp(a * x * x + b * x + c);
J_(i, 0) = dfda;
J_(i, 1) = dfdb;
J_(i, 2) = dfdc;
fx_(i, 0) = y - exp(a * x * x + b * x + c);
}
}

void LM::LevenbergMarquardt::CalcH_g()
{
H_ = J_.transpose() * J_;
g_ = - J_.transpose() * fx_;
}

double LM::LevenbergMarquardt::GetCost()
{
MatrixXd cost = fx_.transpose() * fx_;
return cost(0, 0);
}

double LM::LevenbergMarquardt::F(double da, double db, double dc)
{
Eigen::VectorXd fx;
fx.resize(obsPt.size(), 1);

for (int i = 0; i < obsPt.size(); i++)
{
const Eigen::Vector2d obs = obsPt.at(i);
const double x = obs(0);
const double y = obs(1);
fx(i) = y - exp(da * x * x + db * x + dc);
}
Eigen::MatrixXd F = 0.5 * fx.transpose() * fx;
return F(0, 0);
}

double LM::LevenbergMarquardt::L0_L(Eigen::Vector3d h)
{
/*Eigen::MatrixXd L = -fx_.transpose() * J_ * h - 0.5 * h.transpose() * J_ * J_.transpose() * h;*/

Eigen::MatrixXd L = -h.transpose() * J_.transpose() * fx_ - 0.5 * h.transpose() * J_.transpose() * J_ * h;
return L(0, 0);
}

void LM::LevenbergMarquardt::Slove()
{
int k = 0;
double vu = 2.0;
CalcJ_fx();
CalcH_g();
bool found = (g_.lpNorm<Eigen::Infinity>() < epsilon1); //误差足够小
std::vector<double> H_diag;
for (int i = 0; i < 3; i++)
{
H_diag.push_back(H_(i, i));
}
auto max = std::max_element(H_diag.begin(), H_diag.end());
double mu = *max;
double tsum = 0;
while (!found && k < MaxIter) // 误差足够小，增量足够小，达到迭代次数
{
k++;
RunTime runtime;
runtime.Stat();
Eigen::Matrix3d G = H_ + mu * Eigen::Matrix3d::Identity();
Eigen::Vector3d h = G.ldlt().solve(g_);

if (h.norm() < epsilon2 * sqrt(a * a + b * b + c * c) + epsilon2) //增量足够小
{
found = true;
}
else
{
double a_new = a + h(0);
double b_new = b + h(1);
double c_new = c + h(2);
double sig = (F(a, b, c) - F(a_new, b_new, c_new)) / L0_L(h);
if (sig > 0)
{
a = a_new;
b = b_new;
c = c_new;
CalcJ_fx();
CalcH_g();
found = (g_.lpNorm<Eigen::Infinity>() < epsilon1);
mu = mu * std::max<double>(0.333, 1 - std::pow(2 * sig - 1, 3));
vu = 2.0;
}
else
{
mu = mu * vu;
vu *= 2;
}
}

runtime.Stop();
if (isOut)
{
std::cout << "Iter: " << std::left << std::setw(5) << k << std::left << std::setw(10)
<< "Result: " << std::left << std::setw(10) << a << " " << std::left << std::setw(10) << b << " " << std::left << std::setw(10) << c << std::left << std::setw(10)
<< "step: " << std::left << std::setw(15) << h.norm() << std::left << std::setw(10)
<< "cost: " << std::left << std::setw(10) << GetCost() << std::left << std::setw(10)
<< "time: " << std::left << std::setw(10) << (tsum += runtime.Duration()) << std::endl;
}
}

if (found)
{
std::cout << "\nConverged\n\n";
}
else
{
std::cout << "\nDiverged\n\n";
}
}


展开全文
• 列文伯格(1944)和马夸尔特(1963)先后对高斯牛顿法进行了改进，求解过程中引入了阻尼因子，将无约束最小二乘问题转变为有约束最小二乘问题。

# 1、介绍

列文伯格(1944)和马夸尔特(1963)先后对高斯牛顿法进行了改进，求解过程中引入了阻尼因子。
公式 36 的无约束最小二乘问题转变为公式 44 有约束最小二乘问题，其中， 1 2 × ( ∥ D Δ X k ∥ 2 − μ ) ⩽ 0 \frac{1}{2} \times ( \begin{Vmatrix} D\Delta X_k \end{Vmatrix}^2 -\mu )\leqslant 0 是信赖区间对应的条件， D D 是系数矩阵(列文伯格方法和马夸尔特方法的主要区别就是系数矩阵 D D 不同)， μ \mu 是信赖 半径。
F ( X k + Δ X k ) ≈ 1 2 ∥ L ( X k ) + J ( X k ) ⏟ L T Δ X k ∥ 2 , s . t . ( 1 2 × ( ∥ D Δ X k ∥ 2 − μ ) ⩽ 0 ) ( 公 式 44 ) F(X_k+\Delta X_k)\approx \frac{1}{2} \begin{Vmatrix} \\ L(X_k)+\underbrace{J(X_k)}_{L}{^T} \Delta X_k \end{Vmatrix}^2 \qquad ,s.t.( \frac{1}{2} \times ( \begin{Vmatrix} D\Delta X_k \end{Vmatrix}^2 -\mu )\leqslant 0 ) \qquad (公式44)

# 2、数学原理

将公式 44 的约束条件引入，可以定义拉格朗日函数 L a g ( X ) Lag(X)
L a g ( Δ X ) = d e f 1 2 ∥ L ( X k ) + J ( X k ) ⏟ L T Δ X k ∥ 2 + λ k × 1 2 × ( ∥ D Δ X k ∥ 2 − μ ) , s . t . ( λ k ⩾ 0 ) ( 公 式 45 ) Lag(\Delta X)\stackrel{\mathrm{def}}{=} \frac{1}{2} \begin{Vmatrix} \\ L(X_k)+\underbrace{J(X_k)}_{L}{^T} \Delta X_k \end{Vmatrix}^2 +{\lambda}_k \times \frac{1}{2} \times ( \begin{Vmatrix} D\Delta X_k \end{Vmatrix}^2 -\mu ) \qquad ,s.t.( {\lambda}_k \geqslant 0 ) \qquad (公式45)
公式45中 J ( X k ) ⏟ L \underbrace{J(X_k)}_{L} 是残差函数 L ( X k ) L(X_k) 的雅可比矩阵。当 ∥ D Δ X k ∥ 2 ⩾ 0 \begin{Vmatrix} D\Delta X_k \end{Vmatrix}^2 \geqslant 0 的时候，在优化目标函数 F ( X ) F(X) 中引入阻尼因子 λ k {\lambda}_k 作为惩罚项。
该表达式中 L a g ( Δ X k ) Lag(\Delta X_k) L ( X k ) L(X_k) 是一个常数， J ( X k ) ⏟ L \underbrace{J(X_k)}_{L} D D 是一个常数矩阵， Δ X k \Delta X_k 是一个变量矩阵，即函数 L a g ( Δ X k ) Lag(\Delta X_k) 是以 Δ X k \Delta X_k 为自变量的二次函数。综上所述，当 L a g ( Δ X k ) Lag(\Delta X_k) 一阶导数为0的时候，函数 L a g ( Δ X k ) Lag(\Delta X_k) 取得极值，可推得：
L a g Δ X k ′ ( Δ X k ) = 0 ( 公 式 46 ) Lag'_{\Delta X_k}(\Delta X_k)=0 \qquad (公式46)
由公式 45、公式 46 可推得：
0 = λ k D T D Δ X k + L ( X k ) J ( X k ) ⏟ L + J ( X k ) ⏟ L J ( X k ) ⏟ L T Δ X k ( 公 式 47 ) 0= {\lambda}_k D^TD \Delta X_k + L(X_k) \underbrace{J(X_k)}_{L} + \underbrace{J(X_k)}_{L} \underbrace{J(X_k)}_{L}{^T} \Delta X_k \qquad (公式47)
即：
0 = λ k D T D Δ X k + ∑ i = 1 m ( L i ( X k ) J ( X k ) ⏟ L i + J ( X k ) ⏟ L i J ( X k ) ⏟ L i T Δ X k ) ( 公 式 48 ) 0= {\lambda}_k D^TD \Delta X_k+ \sum_{i=1}^m ( L_i(X_k) \underbrace{J(X_k)}_{L_i} + \underbrace{J(X_k)}_{L_i} \underbrace{J(X_k)}_{L_i}{^T} \Delta X_k ) \qquad (公式48)
由公式 47 可推得:
0 = L ( X k ) J ( X k ) ⏟ L + ( J ( X k ) ⏟ L J ( X k ) ⏟ L T + λ k D T D ) Δ X k ( 公 式 49 ) 0= L(X_k) \underbrace{J(X_k)}_{L} + ( \underbrace{J(X_k)}_{L} \underbrace{J(X_k)}_{L}{^T} +{\lambda}_k D^TD ) \Delta X_k \qquad (公式49)
设：由公式49结合公式26结构形式，可近似推得函数 F ( X ) F(X) 的黑塞矩阵 H ( X k ) ⏟ F \underbrace{H(X_k)}_{F} 和雅克比矩阵 J ( X k ) ⏟ F \underbrace{J(X_k)}_{F}
H ( X k ) ⏟ F ≈ d e f J ( X k ) ⏟ L J ( X k ) ⏟ L T ( 公 式 50 ) \underbrace{H(X_k)}_{F} \stackrel{\mathrm{def}}{\approx} \underbrace{J(X_k)}_{L} \underbrace{J(X_k)}_{L}{^T} \qquad (公式50)
J ( X k ) ⏟ F ≈ d e f L ( X k ) J ( X k ) ⏟ L ( 公 式 51 ) \underbrace{J(X_k)}_{F} \stackrel{\mathrm{def}}{\approx} L(X_k) \underbrace{J(X_k)}_{L} \qquad (公式51)
由公式 49、公式 50、公式 51 可推得：
Δ X k = − ( H ( X k ) ⏟ F + λ k D T D ) − 1 J ( X k ) ⏟ F ( 公 式 52 ) \Delta X_k=- ( {\underbrace{H(X_k)}_{F}} + {\lambda}_k D^TD )^{-1} \underbrace{J(X_k)}_{F} \qquad (公式52)
由公式 52 可推得目标函数 F ( X ) F(X) 的最优化迭代公式：
X k + 1 = d e f X k − ( H ( X k ) ⏟ F + λ k D T D ) − 1 J ( X k ) ⏟ F ( 公 式 53 ) X_{k+1}\stackrel{\mathrm{def}}{=} X_{k} -( {\underbrace{H(X_k)}_{F}} + {\lambda}_k D^TD )^{-1} \underbrace{J(X_k)}_{F} \qquad (公式53)
当阻尼因子 λ k > 0 {\lambda}_k>0 的时候，保证 H ( X k ) ⏟ F + λ k D T D {\underbrace{H(X_k)}_{F}} + {\lambda}_k D^TD 正定，迭代朝着下降方向进行；当 λ k {\lambda}_k 趋近于 0 的时候，退化为高斯牛顿算法；当 λ k {\lambda}_k 趋近于正无穷大的时候，退化为最速下降算法；

# 3、阻尼因子更新策略

目标函数 F ( X ) F(X) X = X k X=X_k 处的不含皮亚诺余项的二阶泰勒公式如下：
F ( X k + Δ X k ) ≈ F ( X k ) + J ( X k ) ⏟ F T Δ X k ＋ 1 2 Δ X k T H ( X k ) ⏟ F Δ X k ( 公 式 54 ) F(X_k+\Delta X_k)\approx F(X_k)+ \underbrace{J(X_k)}_{F}{^T} \Delta X_k ＋\frac{1}{2} {\Delta X_k}^T \underbrace{H(X_k)}_{F} \Delta X_k \qquad (公式54)
根据公式 50、公式 51 的近似关系和公式 54 可推得：
U ( Δ X k ) = d e f F ( X k + Δ X k ) ≈ F ( X k ) + J ( X k ) ⏟ L T L ( X k ) Δ X k ＋ 1 2 Δ X k T J ( X k ) ⏟ L T J ( X k ) ⏟ L Δ X k ( 公 式 55 ) U(\Delta X_k)\stackrel{\mathrm{def}}{=} F(X_k+\Delta X_k)\approx F(X_k)+ {\underbrace{J(X_k)}_{L}}^T L(X_k) \Delta X_k ＋\frac{1}{2} {\Delta X_k}^T {\underbrace{J(X_k)}_{L}}^T {\underbrace{J(X_k)}_{L}} \Delta X_k \qquad (公式55)
阻尼因子 λ k {\lambda}_k 的初始值 λ 0 {\lambda}_0 的选取， H 0 ⏟ F ≈ d e f J ( X 0 ) ⏟ L J ( X 0 ) ⏟ L T \underbrace{H_0}_{F} \stackrel{\mathrm{def}}{\approx} \underbrace{J(X_0)}_{L} \underbrace{J(X_0)}_{L}{^T} 对角线上面的数 h i i h_{ii} 有关。
λ 0 = d e f τ × m a x { h i i } , s . t . ( τ ∈ [ 1 0 − 8 , 1 ] ) ( 公 式 56 ) {\lambda}_0\stackrel{\mathrm{def}}{=} \tau \times max\begin{Bmatrix} h_{ii} \end{Bmatrix} \qquad ,s.t.( \tau \in[10^{-8},1] ) \qquad (公式56)
而阻尼因子 λ k {\lambda}_k 的更新可以由增益比 β k {\beta}_k 来定量分析：
β k = d e f F ( X k ) − F ( X k + Δ X k ) U ( 0 ) − U ( Δ X k ) ( 公 式 57 ) {\beta}_k\stackrel{\mathrm{def}}{=} \frac{F(X_k)-F(X_k+\Delta X_k)}{U(0)-U(\Delta X_k)} \qquad (公式57)
其中公式57的分子为目标函数 F ( X k ) F(X_k) 的值在步长 Δ X k \Delta X_k 下的实际变化，分母为目标函数 F ( X k ) F(X_k) 的值二阶近 似的变化，可得出：
如果 β k {\beta}_k 的值越大，那么目标函数的二阶近似变化对实际变化效果越好，可以缩小 λ k {\lambda}_k 的值接近高斯牛顿算法;
如果 β k {\beta}_k 的值越小，那么目标函数的二阶近似变化对实际变化效果越差，可以增大 λ k {\lambda}_k 的值接近最速下降算法;
阻尼因子 λ k {\lambda}_k Nielsen \text {Nielsen} (1991)更新策略如下：
如 果 ( β k > 0 ) : λ k = d e f λ k − 1 × m a x { 1 3 ， 1 − ( 2 β k − 1 ) 3 } v k = d e f 2 否 则 ( β k ⩽ 0 ) : λ k = d e f λ k − 1 × v k − 1 v k = 2 × v k − 1 ( 公 式 58 ) \begin{aligned} 如果(\beta_k>0): & \qquad \\ &\qquad {\lambda}_k\stackrel{\mathrm{def}}{=} {\lambda}_{k-1} \times max\begin{Bmatrix} \frac{1}{3}，1-(2\beta_k-1)^3 \end{Bmatrix}\\ &\qquad v_k\stackrel{\mathrm{def}}{=}2\\ 否则(\beta_k \leqslant 0): &\\ &\qquad {\lambda}_k\stackrel{\mathrm{def}}{=} {\lambda}_{k-1} \times v_{k-1} \\ &\qquad v_k=2 \times v_{k-1} \end{aligned} \qquad (公式58)

# 4、列文伯格方法

在列文伯格方法中，系数矩阵 D D 是单位矩阵 I I ，信赖区域是一个球形状。

# 5、马夸尔特方法

在马夸尔特方法中，系数矩阵 D D H ( X k ) ⏟ F ≈ d e f J ( X k ) ⏟ L J ( X k ) ⏟ L T \underbrace{H(X_k)}_{F} \stackrel{\mathrm{def}}{\approx} \underbrace{J(X_k)}_{L} \underbrace{J(X_k)}_{L}{^T} 对角元素的平方根，信赖区域是一个椭球形状。

# 6、Matlab程序

链接地址

展开全文
• 前言 LM方法是适用于求解方程最小值的一种方法，在非线性优化的框架中，优化方法分为Line Search 和 Trust Region，也就是线搜索和信任域方法，它们是两种不同性质的方法。 不同之处: LIne Search:不管当前迭代点...

# 前言

LM方法是适用于求解方程最小值的一种方法，在非线性优化的框架中，优化方法分为Line Search 和 Trust Region，也就是线搜索和信任域方法，它们是两种不同性质的方法。

不同之处:

LIne Search:不管当前迭代点X(k)到最优解X*之间的路径，每次迭代X(k)得到X(k+1),都是使用该点的反向梯度方向进行值得寻找，这就导致了这样一种可能得问题：‘在靠近X*的时候，X(k)反复震荡，不容易收敛。这种情况特别是会发生在固定迭代步长的时候，Eg:对于求解f(x)=x^2的最优解X*,当X(k)=0.01,而迭代步长为1的时候，这个时候下一个X(k+1)=-0.01，如此这就发生了震荡，无法收敛。这确实是特例，但是只要有这种情况存在，也就是我们能够找出一个情况算法有问题，我们就无法保证所需要优化的函数f(x)不存在这种性质。

代表:梯度下降法，牛顿法，高斯牛顿法。

Trust Region:对于当前迭代点X(k)，我们会找出一个合适的迭代区域||D(X(k+1)-X(k))||<=hk,hk>0，下一个X(k+1)产生于该区域内，然后寻找这个区域的内最佳的迭代向量dX(k)。同时uk会根据每次的X(k)来计算调整，这就i是信任域的普遍思想。

代表:LM法。

# LM法的个人理解

具体的实现方法我不做说明，网络上已经有很多人进行了详细的介绍。

这里需要说明的是，很多博主在描述到LM具体算法实现的伪代码描述时，大多会混淆拉格朗日乘子$\lambda$和信任域hk这两个参量。导致在说到具体实现的时候逻辑不清，在应该让h(k+1)=2*h(k)的时候，说成h(k+1)=h(k)/2。

下面来说明一下算法和LM思想原理之间的对应关系:

高斯牛顿法优化目标是$\bigtriangleup x*=arg min \frac{1}{2}||f(x)+J(x)\bigtriangleup x||$，利用求导数为0，我们得到$\bigtriangleup x*=-(J(x)^{T}J(x))^{-1}J(x)f(x)$,但是这个要求$J(x)^{T}J(x)$必须是可逆矩阵，要求较为严格。

而LM算法优化目标是$\bigtriangleup x*=arg min \frac{1}{2}||f(x)+J(x)\bigtriangleup x||,st.||D\Delta x||\leqslant h(k)$,因为约束条件的h(k)是迭代变化的，为了保证公式一致性，把约束条件改写成$g(\Delta x)=\frac{||D\Delta x||}{h(k)}-1\leq 0$是加了约束条件的优化目标.利用拉格朗日乘子法进行转换成无约束目标$\bigtriangleup x*=arg min \frac{1}{2}||f(x)+J(x)\bigtriangleup x||+\frac{\lambda}{2 }(\frac{||D\Delta x||}{h(k)}-1)$，令$u(k)=\frac{\lambda }{h(k)}$

采用拉格朗日乘数法确实需要解出不同情况下对应的$\lambda$*，但是在LM法中进行了这样的操作，把$\lambda$看成常数，相等于把问题转化成了数学建模中的规划问题，约束条件的权重为$\frac{\lambda }{\lambda +1}$,而高斯牛顿法优化目标为$\frac{1}{\lambda +1}$,当$\lambda =1$时，两者同比重，也就是总优化目标同时考虑到了两者。这样我们就可以减少求解乘子的步骤，得到最终的迭代向量$\bigtriangleup x*=-(J(x)^{T}J(x)+u(k)D^{T}D)^{-1}J(x)f(x)$，特别的当D为单位矩阵$I$时，表示信任域是一个球。

接下来讲述关键问题，也就是如何迭代u(k).

LM定义了这样一个比值参量$\rho(x,\Delta x) =\frac{f(x+\Delta x)-f(x)}{J(x)\Delta x}$,表征了实际下降量与一阶微分量（也就是近似下降量)的比值。

在算法伪代码中有这么一个判断:

如果$\rho(x(k),\Delta x)\leq 0.25$,那么h(k+1)=h(k)/2,也就是减小信任域，此时在具体的算法中u(k+1)=2*u(k)；

如果$\rho(x(k),\Delta x)\geq 0.75$,那么h(k+1)=2*h(k),也就是增大信任域，此时在具体的算法中u(k+1)=u(k)/2；

这其中的逻辑关系从LM迭代公式$\bigtriangleup x*=-(J(x)^{T}J(x)+u(k)D^{T}D)^{-1}J(x)f(x)$可以知道，因为如果实际下降比一阶近似小很多，那么表示二阶以上分量占了比较大的比重，这个时候非线性比较严重，那么就应该缩小信任域，u(k)增大，而$\bigtriangleup x*$的模长变小，也就是变化较小了。同理可说明另外的情况。

最后说明一下如果u(k)的迭代搞混了，也就是之前提到的，该增大的时候变小了，而该变小的时候增大了，这个时候会出现的问题。

最容易遇到的问题就是$\Delta x$迭代过小，通过查看$\Delta x$的模长||$\Delta x$||$\approx$2^(-k)，就是同一个数量级，而且收敛于非最优解，过早收敛了。产生这种情况的逻辑在于:当某一次迭代过程得到$\rho(x(k),\Delta x)\geq 0.75$,本该u(k+1)=u(k)/2,但是搞混的情况下u(k+1)=2*u(k),下次的$\Delta x$迭代将会变得更小，因为u(k+1)=2*u(k)等价于信任域变小了。另外$\lim_{||\Delta x||}\rho(x,\Delta x) =\frac{f(x+\Delta x)-f(x)}{J(x)\Delta x}=1$，这样一次比一次的迭代小了，而且还是指数变小，那么就会使得目标在非最优解上收敛。通过修改程序上的正常逻辑，我们可以验证这样的情况。

# 最后附上LM算法的matlab程序

链接：https://pan.baidu.com/s/1_3AF-ZRkU9DPlcgxoMY-DA
提取码：z7y8
复制这段内容后打开百度网盘手机App，操作更方便哦

展开全文
• 分享一篇写的非常好的博客LM算法，简单明了地介绍了梯度下降法、牛顿法、高斯牛顿法和Levenberg-Marquardt这几种优化算法之间的联系和区别。
• #高斯牛顿法
• 它是使用最广泛的非线性最小二乘算法，中文为列文伯格-马夸尔特法。它是利用梯度求最大（小）值的算法，形象的说，属于“爬山”法的一种。它同时具有 梯度法 和 牛顿法 的优点。当λ很小时，步长等于牛顿法步长，当...
• 列文伯格-马夸尔特 (Levenberg-Marquardt) 算法 ，简称LM 算法 。 下面我们详细探讨一下，高斯牛顿和LM算法的原理以及在VSLAM中的应用。 首先，最小二乘是要解决什么问题？  1 最小二乘算法    1.1 ...
• 版权声明：Davidwang原创文章，严禁用于任何商业途径，授权后方可转载...列文伯格马夸尓法针对这种情况进行了改进，即对$J(x)J^T(x)$矩阵进行改造，添加了一个参数值，保证最后得到的矩阵为正定矩阵，即增量方程变为
• 高斯牛顿(Gauss Newton)、列文伯格-马夸尔特(Levenberg-Marquardt)最优化算法与VSLAM中的具体应用
• 列文伯格-马夸尔特法(Levenberg-Marquardt)进行计算，该算法既有高斯-牛顿法的局部收敛性，又具有梯度下降法的全局性。
• 首先谈一下应用场景——在拟合的时候进行应用 什么是拟合？你有一堆数据点，我有一个函数，但是这个函数的很多参数是未知的，我只知道你的这些数据点都在我的函数上，因此我可以用你的数据点来求我的函数的未知参数...
•   在《SLAM之路－列文伯格马夸尓法Cpp实现(元素视角)》这篇文章中，我们直接通过for循环迭代获取增量方程，看起来很冗余，不直观，在本文中，我们将对上篇文章进行改进，使用矩阵表达方式实现算法。   依然...
• 第5部分介绍列文伯格-马夸特法； 第6部分介绍几种方法的联系。 1、最小二乘 我们想要求解的问题是SLAM过程中的最大似然估计问题。由于最大似然估计问题可以转换为最小二乘问题，所以也就是最小二乘的求解问题。 ...
• 前段时间要用到这个拟合算法，于是在网上搜了下发现了一段matlab代码，写的非常详细，马上改写之。假设我要拟合的函数为a * atan(b * x + c) + d，代码如下： 1 using Accord.Math; 2 public class ...
• 1.状态估计问题 对于SLAM经典模型，我们知道是由一个运动方程和一个观测方程构成，如下方程： {xk=f(xk−1,uk)+wkzk,j=f(yj,xk)+vk,j{xk=f(xk−1,uk)+wkzk,j=f(yj,xk)+vk,j\left\{\begin{matrix} \mathit{x}_{k}=f...
• ## Levenberg-Marquardt(LM算法)

万次阅读 多人点赞 2016-03-30 14:18:21
它是使用最广泛的非线性最小二乘算法，中文为列文伯格-马夸尔特法。它是利用梯度求最大（小）值的算法，形象的说，属于“爬山”法的一种。它同时具有梯度法和牛顿法的优点。当λ很小时，步长等于牛顿法步长，当λ很...
• 在高翔slambook第七讲非线性优化方法中，对于PnP问题选择的优化方法是g2o中的L-M方法，而对于ICP问题选择的是G-N法，... PnP中说g2o中没有G-N方法： ... 而且早在第六讲提出非线性优化时，在用g2o拟合曲线时我尝试选用...
• 用图表示单个点重投影误差是这样的： 求解上述BA问题，最著名的就是高斯-牛顿法（G-M）和列文伯格-马夸尔特法（L-M）。这里我们用高斯牛顿法。 1.2 高斯牛顿法求解BA 高斯牛顿法是最优化算法里面最简单的方法之一。...
• 列文伯格-马夸特算法的例子 采用Levenberg-Marquardt算法实现照相机调焦，包括matlab例子代码 HW5_LM_handout.pdf Figure8.jpg Levenberg-Marquardt算法实现照相机...

...