ceres之LM算法

Ceres作为一个优化算法库，在许多领域中有着至关重要的作用，比如slam系统中的优化问题-集束调整BA，就可以通过Ceres去实现，官方文档地址：http://ceres-solver.org/nnls_tutorial.html#bundle-adjustment

本文主要是解析ceres中的LM算法过程，参考代码地址：

https://github.com/ceres-solver/ceres-solver/tree/master/internal/ceres

一、主要流程

先贴个图，LM算法的大概流程如下

可以看到，LM算法的输入为（1）雅可比矩阵J(x)；（2）残差向量f(x)；（3）待优化变量初值x0；（4）控制参数等

LM算法要求解的问题为： $x^{*} = argmin(Fx) = argmin(loss(f(x)^{2})))$

其中 $f(x)$ 为残差函数，它的导函数为 $f'(x) = J$ ,二阶导函数的近似为 $f''(x)=J^T*J$

分为几个步骤：

（1）初始化：首先计算系数矩阵A和残差向量g,初始化参数

（2）while循环：如果达到收敛条件就停止迭代

（3）求解（A+miu*I)dx = -g

(4) xnew = x+dx

（5） $\rho = (Fx - Fx_{new})/(L(0)-L(dx))$ ,其中Fx,Fxnew是所有残差的平方和

$L(0) - L(dx) = Fx - (Fx+F'(x)*dx+0.5*F''(x)*dx^{2}) = -dx^{T}*J^{T}*f-0.5*(J*dx)^T*(J*dx)=-(J*dx)^T*((J*dx)/2+f)$

(6)如果第五步结果大于零，表示这个迭代是有效的，可以接受

然后更新x = xnew;Fx =Fnew ,A = J'*J, g = J'*f

(7)否则这个dx得到的结果是无效的，收缩搜索半径，相当于增大 $\mu$

ceres对应代码：https://github.com/ceres-solver/ceres-solver/blob/master/internal/ceres/trust_region_minimizer.cc

以及：https://github.com/ceres-solver/ceres-solver/blob/master/internal/ceres/levenberg_marquardt_strategy.cc

void TrustRegionMinimizer::Minimize(const Minimizer::Options& options,double* parameters,Solver::Summary* solver_summary) {start_time_in_secs_ = WallTimeInSeconds();iteration_start_time_in_secs_ = start_time_in_secs_;
//初始化Init(options, parameters, solver_summary);
//得到初始的雅可比矩阵以及残差矩阵RETURN_IF_ERROR_AND_LOG(IterationZero());// Create the TrustRegionStepEvaluator. The construction needs to be// delayed to this point because we need the cost for the starting// point to initialize the step evaluator.step_evaluator_.reset(new TrustRegionStepEvaluator(x_cost_,options_.use_nonmonotonic_steps? options_.max_consecutive_nonmonotonic_steps: 0));
//while循环：是否达到迭代终止条件while (FinalizeIterationAndCheckIfMinimizerCanContinue()) {iteration_start_time_in_secs_ = WallTimeInSeconds();iteration_summary_ = IterationSummary();iteration_summary_.iteration =solver_summary->iterations.back().iteration + 1;
//计算(J'*J+D'*D/radius)*dx = -J'*fRETURN_IF_ERROR_AND_LOG(ComputeTrustRegionStep());if (!iteration_summary_.step_is_valid) {
//如果当前迭代不是有效的，则收缩搜索半径RETURN_IF_ERROR_AND_LOG(HandleInvalidStep());continue;}if (options_.is_constrained &&options_.max_num_line_search_step_size_iterations > 0) {// Use a projected line search to enforce the bounds constraints// and improve the quality of the step.DoLineSearch(x_, gradient_, x_cost_, &delta_);}
//计算Xnew = x+dxComputeCandidatePointAndEvaluateCost();DoInnerIterationsIfNeeded();
//是否||dx||太小了，太小了直接终止if (ParameterToleranceReached()) {return;}
//是否||Fx-Fxnew||太小了，太小了意味着残差平方和更新不动直接终止if (FunctionToleranceReached()) {return;}
//如果是有效的if (IsStepSuccessful()) {
//更新x,雅可比矩阵，残差向量等RETURN_IF_ERROR_AND_LOG(HandleSuccessfulStep());continue;}
//否则收缩搜索半径HandleUnsuccessfulStep();}
}

上面为什么是信赖域算法过程，Ceres里面把LM算法和dogleg算法（也叫狗腿算法）集成到统一的框架下--信赖域算法框架，不同的是LM算法求解dx的过程和狗腿算法不同，下面是LM算法求解dx的过程以及搜索半径的更新

TrustRegionStrategy::Summary LevenbergMarquardtStrategy::ComputeStep(const TrustRegionStrategy::PerSolveOptions& per_solve_options,SparseMatrix* jacobian,const double* residuals,double* step) {CHECK(jacobian != nullptr);CHECK(residuals != nullptr);CHECK(step != nullptr);//计算对角矩阵Dconst int num_parameters = jacobian->num_cols();if (!reuse_diagonal_) {if (diagonal_.rows() != num_parameters) {diagonal_.resize(num_parameters, 1);}jacobian->SquaredColumnNorm(diagonal_.data());for (int i = 0; i < num_parameters; ++i) {diagonal_[i] = std::min(std::max(diagonal_[i], min_diagonal_),max_diagonal_);}}
// D = diag{J'*J}/radiuslm_diagonal_ = (diagonal_ / radius_).array().sqrt();LinearSolver::PerSolveOptions solve_options;solve_options.D = lm_diagonal_.data();solve_options.q_tolerance = per_solve_options.eta;solve_options.r_tolerance = -1.0;InvalidateArray(num_parameters, step);
//求解 （A+D'*D)dx = -gLinearSolver::Summary linear_solver_summary =linear_solver_->Solve(jacobian, residuals, solve_options, step);if (linear_solver_summary.termination_type == LINEAR_SOLVER_FATAL_ERROR) {LOG(WARNING) << "Linear solver fatal error: "<< linear_solver_summary.message;} else if (linear_solver_summary.termination_type == LINEAR_SOLVER_FAILURE)  {LOG(WARNING) << "Linear solver failure. Failed to compute a step: "<< linear_solver_summary.message;} else if (!IsArrayValid(num_parameters, step)) {LOG(WARNING) << "Linear solver failure. Failed to compute a finite step.";linear_solver_summary.termination_type = LINEAR_SOLVER_FAILURE;} else {VectorRef(step, num_parameters) *= -1.0;}reuse_diagonal_ = true;if (per_solve_options.dump_format_type == CONSOLE ||(per_solve_options.dump_format_type != CONSOLE &&!per_solve_options.dump_filename_base.empty())) {if (!DumpLinearLeastSquaresProblem(per_solve_options.dump_filename_base,per_solve_options.dump_format_type,jacobian,solve_options.D,residuals,step,0)) {LOG(ERROR) << "Unable to dump trust region problem."<< " Filename base: " << per_solve_options.dump_filename_base;}}TrustRegionStrategy::Summary summary;summary.residual_norm = linear_solver_summary.residual_norm;summary.num_iterations = linear_solver_summary.num_iterations;summary.termination_type = linear_solver_summary.termination_type;return summary;
}//如果当前迭代是有效的，更新搜索半径等参数
void LevenbergMarquardtStrategy::StepAccepted(double step_quality) {CHECK_GT(step_quality, 0.0);radius_ = radius_ / std::max(1.0 / 3.0,1.0 - pow(2.0 * step_quality - 1.0, 3));radius_ = std::min(max_radius_, radius_);decrease_factor_ = 2.0;reuse_diagonal_ = false;
}//如果当前迭代是无效的，收缩搜索半径等参数
void LevenbergMarquardtStrategy::StepRejected(double step_quality) {radius_ = radius_ / decrease_factor_;decrease_factor_ *= 2.0;reuse_diagonal_ = true;
}

ceres之LM算法相关推荐

非线性优化：徒手实现LM算法
点击上方"3D视觉工坊",选择"星标" 干货第一时间送达本文由知乎作者David LEE授权转载,不得擅自二次转载. 原文链接:https://zhuanla ...
lm opencv 算法_相机模型与标定（七）--LM算法在相机标定中的使用
LM算法在相机标定的应用共有三处. (1)单目标定或双目标定中,在内参固定的情况下,计算最佳外参.OpenCV中对应的函数为findExtrinsicCameraParams2. (2)单目标定中,在 ...
lm opencv 算法_Levenberg–Marquardt算法学习（和matlab的LM算法对比）
回顾高斯牛顿算法,引入LM算法惩罚因子的计算(迭代步子的计算) 完整的算法流程及代码样例 1. 回顾高斯牛顿,引入LM算法根据之前的博文:Gauss-Newton算法学习假设我们研究如 ...
matlab lm计算方式,lm算法(lm算法原理介绍)
请问 MATLAB中 LM算法(Levenberg-Marquard-algorithm)的函数是什么?. http://www.mathworks.com/matlabcentral/fileexc ...
信赖域(Trust Region)算法和L-M算法
近期在将样条数据进行公式化. 这时候发现多项高斯函数函数很好用,里面发现了两种优化高斯函数参数的算法一种是信赖域算法(Trust Region),另一种是L-M算法(Levenberg-Marquad ...
相机校正、张氏标定法、极大似然估计/极大似然参数估计、牛顿法、高斯牛顿法、LM算法、sin/cos/tan/cot
日萌社人工智能AI:Keras PyTorch MXNet TensorFlow PaddlePaddle 深度学习实战(不定时更新) CNN:RCNN.SPPNet.Fast RCNN.Faste ...
梯度下降、牛顿法、高斯牛顿L-M算法比较
本文梳理一下常见的几种优化算法:梯度下降法,牛顿法.高斯-牛顿法和L-M算法,优化算法根据阶次可以分为一阶方法(梯度下降法),和二阶方法(牛顿法等等),不同优化算法的收敛速度,鲁棒性都有所不同.一般来 ...
Levenberg-Marquardt(LM算法)的理解
Levenberg-Marquardt LM算法的理解 1. convex optimization 1.1 convex set 1.2 convex function 1.3 optimizat ...
视觉orb_slam中LM算法的分析
orb_slam中LM算法的分析刘康 2017年 1月在orb_slam的PoseOptimization()函数中的LM算法,可以分解成几个部分,首先是线性方程的建立和求解,然后是位姿矩阵的迭代 ...
sba(sparse bundle adjustment):一个基于Levenberg-Marquardt(LM)算法的通用稀疏光束法平差C/C++软件包
[转]sba(sparse bundle adjustment):一个基于Levenberg-Marquardt(LM)算法的通用稀疏光束法平差C/C++软件包 2011-01-04 10:44 转载 ...

ceres之LM算法

ceres之LM算法相关推荐

最新文章

热门文章