精华内容
下载资源
问答
  • 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.AddObservation(x, y);
    	}
    
    	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:
    		chrono::steady_clock::time_point tstat;
    		chrono::steady_clock::time_point tend;
    	};
    
    	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()
    {
    	tstat = chrono::steady_clock::now();
    }
    
    void LM::RunTime::Stop()
    {
    	tend = chrono::steady_clock::now();
    }
    
    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 21×(DΔXk2μ)0 是信赖区间对应的条件, D D 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) F(Xk+ΔXk)21L(Xk)+L J(Xk)TΔXk2,s.t.(21×(DΔXk2μ)0)(44)

    2、数学原理

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

    3、阻尼因子更新策略

    目标函数 F ( X ) F(X) F(X) X = X k X=X_k X=Xk处的不含皮亚诺余项的二阶泰勒公式如下:
    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) F(Xk+ΔXk)F(Xk)+F J(Xk)TΔXk21ΔXkTF H(Xk)ΔXk(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) U(ΔXk)=defF(Xk+ΔXk)F(Xk)+L J(Xk)TL(Xk)ΔXk21ΔXkTL J(Xk)TL J(Xk)ΔXk(55)
    阻尼因子 λ k {\lambda}_k λk的初始值 λ 0 {\lambda}_0 λ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} F H0defL J(X0)L J(X0)T对角线上面的数 h i i h_{ii} hii有关。
    λ 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) λ0=defτ×max{hii},s.t.(τ[108,1])(56)
    而阻尼因子 λ k {\lambda}_k λk 的更新可以由增益比 β k {\beta}_k β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) βk=defU(0)U(ΔXk)F(Xk)F(Xk+ΔXk)(57)
    其中公式57的分子为目标函数 F ( X k ) F(X_k) F(Xk) 的值在步长 Δ X k \Delta X_k ΔXk 下的实际变化,分母为目标函数 F ( X k ) F(X_k) F(Xk) 的值二阶近 似的变化,可得出:
    如果 β k {\beta}_k βk 的值越大,那么目标函数的二阶近似变化对实际变化效果越好,可以缩小 λ k {\lambda}_k λk 的值接近高斯牛顿算法;
    如果 β k {\beta}_k βk 的值越小,那么目标函数的二阶近似变化对实际变化效果越差,可以增大 λ k {\lambda}_k λk 的值接近最速下降算法;
    阻尼因子 λ k {\lambda}_k λk Nielsen \text {Nielsen} 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) (βk>0):(βk0):λk=defλk1×max{311(2βk1)3}vk=def2λk=defλk1×vk1vk=2×vk1(58)

    4、列文伯格方法

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

    5、马夸尔特方法

    在马夸尔特方法中,系数矩阵 D 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} F H(Xk)defL J(Xk)L J(Xk)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法的个人理解

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

    这里转载一个我认为没有问题的他人的博客。https://www.codelast.com/%e5%8e%9f%e5%88%9blm%e7%ae%97%e6%b3%95%e7%9a%84%e5%ae%9e%e7%8e%b0/

    这里需要说明的是,很多博主在描述到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||\approx2^(-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 ...
  • SLAM学习——非线性优化

    千次阅读 2017-07-10 23:07:42
    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算法实现照相机...

空空如也

空空如也

1 2 3 4 5 ... 15
收藏数 297
精华内容 118
关键字:

列文伯格马夸尔特