-
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
更多相关内容 -
列文伯格-马夸尔特算法LM算法(基于Python编程语言实现)
2022-04-06 21:09:21列文伯格-马夸尔特算法应用于神经网络网络参数优化 -
LM(列文伯格-马夸尔特)方法的matlab实现代码
2020-08-04 22:15:12实现较为粗略的列文伯格-马夸尔特方法代码实现,因为matlab代码简单,很容易修改参数和进行变种。配合https://blog.csdn.net/Pengcode/article/details/107783834进行LM研究将会更加合理。 -
基于C++ Eigen的列文伯格马夸尔特(LM)优化算法
2021-12-16 21:31:20#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"; } }
-
机器学习训练算法十(列文伯格-马夸尔特法(LM 法))
2022-01-09 15:28:56列文伯格(1944)和马夸尔特(1963)先后对高斯牛顿法进行了改进,求解过程中引入了阻尼因子,将无约束最小二乘问题转变为有约束最小二乘问题。连续函数的最优化方法-LM法
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ΔXk∥∥2−μ)⩽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)≈21∥∥∥∥∥∥L(Xk)+L J(Xk)TΔXk∥∥∥∥∥∥2,s.t.(21×(∥∥DΔXk∥∥2−μ)⩽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)=def21∥∥∥∥∥∥L(Xk)+L J(Xk)TΔXk∥∥∥∥∥∥2+λk×21×(∥∥DΔXk∥∥2−μ),s.t.(λk⩾0)(公式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ΔXk∥∥2⩾0的时候,在优化目标函数 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=1∑m(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ΔXk+21Δ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)ΔXk+21Δ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 H0≈defL 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.(τ∈[10−8,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):否则(βk⩽0):λk=defλk−1×max{31,1−(2βk−1)3}vk=def2λk=defλk−1×vk−1vk=2×vk−1(公式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(列文伯格-马夸尔特)方法的个人理解,以及实现问题
2020-08-04 21:56:51前言 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具体算法实现的伪代码描述时,大多会混淆拉格朗日乘子
和信任域hk这两个参量。导致在说到具体实现的时候逻辑不清,在应该让h(k+1)=2*h(k)的时候,说成h(k+1)=h(k)/2。
下面来说明一下算法和LM思想原理之间的对应关系:
高斯牛顿法优化目标是
,利用求导数为0,我们得到
,但是这个要求
必须是可逆矩阵,要求较为严格。
而LM算法优化目标是
,因为约束条件的h(k)是迭代变化的,为了保证公式一致性,把约束条件改写成
是加了约束条件的优化目标.利用拉格朗日乘子法进行转换成无约束目标
,令
。
采用拉格朗日乘数法确实需要解出不同情况下对应的
*,但是在LM法中进行了这样的操作,把
看成常数,相等于把问题转化成了数学建模中的规划问题,约束条件的权重为
,而高斯牛顿法优化目标为
,当
时,两者同比重,也就是总优化目标同时考虑到了两者。这样我们就可以减少求解乘子的步骤,得到最终的迭代向量
,特别的当D为单位矩阵
时,表示信任域是一个球。
接下来讲述关键问题,也就是如何迭代u(k).
LM定义了这样一个比值参量
,表征了实际下降量与一阶微分量(也就是近似下降量)的比值。
在算法伪代码中有这么一个判断:
如果
,那么h(k+1)=h(k)/2,也就是减小信任域,此时在具体的算法中u(k+1)=2*u(k);
如果
,那么h(k+1)=2*h(k),也就是增大信任域,此时在具体的算法中u(k+1)=u(k)/2;
这其中的逻辑关系从LM迭代公式
可以知道,因为如果实际下降比一阶近似小很多,那么表示二阶以上分量占了比较大的比重,这个时候非线性比较严重,那么就应该缩小信任域,u(k)增大,而
的模长变小,也就是变化较小了。同理可说明另外的情况。
最后说明一下如果u(k)的迭代搞混了,也就是之前提到的,该增大的时候变小了,而该变小的时候增大了,这个时候会出现的问题。
最容易遇到的问题就是
迭代过小,通过查看
的模长||
||
2^(-k),就是同一个数量级,而且收敛于非最优解,过早收敛了。产生这种情况的逻辑在于:当某一次迭代过程得到
,本该u(k+1)=u(k)/2,但是搞混的情况下u(k+1)=2*u(k),下次的
迭代将会变得更小,因为u(k+1)=2*u(k)等价于信任域变小了。另外
,这样一次比一次的迭代小了,而且还是指数变小,那么就会使得目标在非最优解上收敛。通过修改程序上的正常逻辑,我们可以验证这样的情况。
最后附上LM算法的matlab程序
链接:https://pan.baidu.com/s/1_3AF-ZRkU9DPlcgxoMY-DA
提取码:z7y8
复制这段内容后打开百度网盘手机App,操作更方便哦 -
Levenberg-Marquardt列文伯格-马夸尔特算法
2019-04-25 19:47:53分享一篇写的非常好的博客LM算法,简单明了地介绍了梯度下降法、牛顿法、高斯牛顿法和Levenberg-Marquardt这几种优化算法之间的联系和区别。 -
高斯牛顿法与列文伯格--马夸尔特法优缺点
2020-06-01 12:01:50#高斯牛顿法 -
Levenberg-Marquardt(列文伯格-马夸尔特)算法
2017-12-25 20:13:57它是使用最广泛的非线性最小二乘算法,中文为列文伯格-马夸尔特法。它是利用梯度求最大(小)值的算法,形象的说,属于“爬山”法的一种。它同时具有 梯度法 和 牛顿法 的优点。当λ很小时,步长等于牛顿法步长,当... -
高斯牛顿(Gauss Newton)与列文伯格-马夸尔特(Levenberg-Marquardt)迭代算法
2019-08-28 11:33:23列文伯格-马夸尔特 (Levenberg-Marquardt) 算法 ,简称LM 算法 。 下面我们详细探讨一下,高斯牛顿和LM算法的原理以及在VSLAM中的应用。 首先,最小二乘是要解决什么问题? 1 最小二乘算法 1.1 ... -
SLAM之路-列文伯格马夸尓特法Cpp实现(元素视角)
2020-09-01 21:39:51版权声明:Davidwang原创文章,严禁用于任何商业途径,授权后方可转载...列文伯格马夸尓特法针对这种情况进行了改进,即对$J(x)J^T(x)$矩阵进行改造,添加了一个参数值,保证最后得到的矩阵为正定矩阵,即增量方程变为 -
高斯牛顿(Gauss Newton)、列文伯格-马夸尔特(Levenberg-Marquardt)最优化算法与VSLAM
2017-07-11 15:01:37高斯牛顿(Gauss Newton)、列文伯格-马夸尔特(Levenberg-Marquardt)最优化算法与VSLAM中的具体应用 -
Model_weibo_matlab_Levenberg-Marquardt_
2021-10-01 13:44:06用列文伯格-马夸尔特法(Levenberg-Marquardt)进行计算,该算法既有高斯-牛顿法的局部收敛性,又具有梯度下降法的全局性。 -
LM算法——列文伯格-马夸尔特算法(最速下降法,牛顿法,高斯牛顿法)(完美解释负梯度方向)
2018-10-23 11:49:00首先谈一下应用场景——在拟合的时候进行应用 什么是拟合?你有一堆数据点,我有一个函数,但是这个函数的很多参数是未知的,我只知道你的这些数据点都在我的函数上,因此我可以用你的数据点来求我的函数的未知参数... -
SLAM之路-列文伯格马夸尓特法Cpp实现(矩阵视角)
2020-10-13 08:56:22在《SLAM之路-列文伯格马夸尓特法Cpp实现(元素视角)》这篇文章中,我们直接通过for循环迭代获取增量方程,看起来很冗余,不直观,在本文中,我们将对上篇文章进行改进,使用矩阵表达方式实现算法。 依然... -
最小二乘问题的四种解法——牛顿法,梯度下降法,高斯牛顿法和列文伯格-马夸特法的区别和联系
2020-07-18 18:31:57第5部分介绍列文伯格-马夸特法; 第6部分介绍几种方法的联系。 1、最小二乘 我们想要求解的问题是SLAM过程中的最大似然估计问题。由于最大似然估计问题可以转换为最小二乘问题,所以也就是最小二乘的求解问题。 ... -
列文伯格-马夸尔特拟合算法(Levenberg Marquardt Fitting)的C#实现
2019-09-20 22:41:02前段时间要用到这个拟合算法,于是在网上搜了下发现了一段matlab代码,写的非常详细,马上改写之。假设我要拟合的函数为a * atan(b * x + c) + d,代码如下: 1 using Accord.Math; 2 public class ... -
SLAM学习——非线性优化
2017-07-10 23:07:421.状态估计问题 对于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中g2o优化方法L-M(列文伯格-马夸尔特方法),G-N(高斯牛顿法)选择的困惑
2019-09-09 22:44:29在高翔slambook第七讲非线性优化方法中,对于PnP问题选择的优化方法是g2o中的L-M方法,而对于ICP问题选择的是G-N法,... PnP中说g2o中没有G-N方法: ... 而且早在第六讲提出非线性优化时,在用g2o拟合曲线时我尝试选用... -
继续视觉里程计求解问题之PnP算法
2021-08-21 12:58:20用图表示单个点重投影误差是这样的: 求解上述BA问题,最著名的就是高斯-牛顿法(G-M)和列文伯格-马夸尔特法(L-M)。这里我们用高斯牛顿法。 1.2 高斯牛顿法求解BA 高斯牛顿法是最优化算法里面最简单的方法之一。... -
Matlab有关LevenbergMarquardt算法实现照相机调焦-HW5_LM_handout.pdf
2019-08-13 02:24:30列文伯格-马夸特算法的例子 采用Levenberg-Marquardt算法实现照相机调焦,包括matlab例子代码 HW5_LM_handout.pdf Figure8.jpg Levenberg-Marquardt算法实现照相机...