當(dāng)前位置：首頁(yè) >

LM优化算法

發(fā)布時(shí)間：2023/12/10 39 豆豆

生活随笔收集整理的這篇文章主要介紹了 LM优化算法小編覺(jué)得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.

LM算法

理論知識(shí)
- - 梯度下降
  - 高斯牛頓
  - Levenberg–Marquardt
- 算法框架
- - 算法的整體流程
  - 求解器update流程
  - 說(shuō)明
- 算法實(shí)現(xiàn)
- - - 頭文件
    - cpp
- 算法調(diào)用

LM優(yōu)化算法，是一種非線性優(yōu)化算法，其可以看作是梯度下降和高斯牛頓法的結(jié)合。綜合了梯度下降對(duì)初值不敏感和高斯牛頓在最優(yōu)值附近收斂速度快的特點(diǎn)。

本人非數(shù)學(xué)專業(yè)，且對(duì)算法理解可能不到位，詳細(xì)的算法推導(dǎo)及各個(gè)優(yōu)化算法之間的關(guān)系，非常推薦看METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS ，其介紹更詳細(xì)也更專業(yè)。

以下簡(jiǎn)要介紹LM算法，然后給出opencv中有關(guān)的實(shí)現(xiàn)。

本文代碼也是來(lái)自opencv相關(guān)代碼，位置在opencv-master\modules\calib3d\src\calibration.cpp\cvFindExtrinsicCameraParams2()。

當(dāng)然，有關(guān)LM算法的實(shí)現(xiàn)，在opencv的usac模塊中有更為正式的實(shí)現(xiàn)，在該模塊中，對(duì)許多算法進(jìn)行了集成和優(yōu)化，后續(xù)再進(jìn)行研究。

理論知識(shí)

對(duì)于一個(gè)最小二乘問(wèn)題，有：

$\frac12\sum_{i=1}^m(f_i(x))^2\\ f_i(x) 為殘差$
我們的目的是求解殘差f(x)的一組系數(shù)，使目標(biāo)函數(shù)（或代價(jià)函數(shù)）F(x)達(dá)到最小值。通常的做法是，給定一個(gè)初值x0，優(yōu)化算法不斷的迭代，產(chǎn)生x1,x2,…，最終收斂于X。

因此，在這個(gè)收斂的過(guò)程中，我們需要兩樣?xùn)|西，即方向h和步長(zhǎng)α
$xk+1=xk+αhdhd即為收斂的方向α即為收斂的步長(zhǎng)x_{k+1} = x_k + \alpha h_d \\ h_d即為收斂的方向 \\ \alpha即為收斂的步長(zhǎng)$
所以，我們看到的梯度下降、高斯牛頓以及本文說(shuō)的LM算法，其思想都是一致的。

梯度下降

這里給出梯度下降的公式，有：
$xk+1=xk?αF˙(x)F˙(x)為F(x)一階導(dǎo)數(shù)x_{k+1} = x_k - \alpha \dot{F}(x) \\ \dot{F}(x)為F(x)一階導(dǎo)數(shù)$

高斯牛頓

高斯牛頓是在牛頓迭代的基礎(chǔ)上，用雅可比矩陣的平方來(lái)近似Hessian矩陣，求解起來(lái)將會(huì)更加高效（大多數(shù)情況下，Hessian矩陣太難求）。這樣做的缺點(diǎn)就是，要求迭代初值必須在最優(yōu)解附近，否則以上假設(shè)將不成立。這里給出其公式：
$xk+1=xk+αhgnJTJhgn=?JTfx_{k+1} = x_k + \alpha h_{gn} \\ J^TJh_{gn} = -J^Tf \\$

Levenberg–Marquardt

給出其公式：
$xk+1=xk+αhlm(JTJ+μI)hlm=?JTfx_{k+1} = x_k + \alpha h_{lm} \\ (J^TJ + \mu I)h_{lm} = -J^Tf$
在有些情況下，JtJ不可逆，導(dǎo)致高斯牛頓無(wú)法使用。LM通過(guò)將JtJ加上一個(gè)μI（I為單位陣），保證了正規(guī)方程的左側(cè)一定可逆。

在算法實(shí)際的應(yīng)用中，通過(guò)調(diào)節(jié)μ的大小，可以使算法在梯度下降和高斯牛頓兩者之間切換，如：

使用較大的μ值，算法可以看作梯度下降，適合當(dāng)前優(yōu)化位置距離最優(yōu)解較遠(yuǎn)的情況
使用較小的μ值，算法可以看作高斯牛頓，適合當(dāng)前優(yōu)化位置接近最優(yōu)解。

算法框架

LM算法通過(guò)求解正規(guī)方程，得到每次迭代的步長(zhǎng)和方向。因此，算法的主要目的是求解正規(guī)方程左右側(cè)的項(xiàng)，然后通過(guò)SVD分解即可得到參數(shù)更新的步長(zhǎng)及方向，即：

JtJ - 雅可比矩陣的轉(zhuǎn)置與其自己相乘
JErr - 雅可比矩陣的轉(zhuǎn)置與殘差矩陣（向量）相乘

算法框架在外部計(jì)算雅可比矩陣和殘差矩陣，然后在算法內(nèi)部通過(guò)求解正規(guī)方程，得到每次參數(shù)更新的步長(zhǎng)及方向。

然后利用更新的參數(shù)，在外部計(jì)算殘差，然后判斷殘差是否朝著收斂方向進(jìn)行。

算法的整體流程

求解器update流程

求解器updata()內(nèi)部根據(jù)不同的狀態(tài)進(jìn)行相應(yīng)的計(jì)算，具體流程如下：

說(shuō)明

J表示雅可比矩陣
err表示殘差矩陣
求解器內(nèi)外同步表示求解器內(nèi)部和外部相應(yīng)的參數(shù)指向相同，便于在求解器外部進(jìn)行雅可比矩陣和殘差矩陣的計(jì)算

算法實(shí)現(xiàn)

算法的實(shí)現(xiàn)主要步驟都體現(xiàn)在上述流程圖中。

頭文件

#pragma once#include <iostream>#include "opencv2/core/types_c.h" #include "opencv2/core/core_c.h"using namespace cv;struct Iteration {int iters = 0;TermCriteria criteria = TermCriteria(0, 0, 0);int lamda_lg10 = 0;};class LevMarq { public:LevMarq();LevMarq(int nparams, int nerrs,TermCriteria criteria = TermCriteria(TermCriteria::EPS + TermCriteria::COUNT, 30, DBL_EPSILON), bool completeSymmFlag = false);~LevMarq();void clear();void initParam(Mat params);void init(int nparams, int nerrs,TermCriteria criteria = TermCriteria(TermCriteria::EPS + TermCriteria::COUNT, 30, DBL_EPSILON), bool completeSymmFlag = false);bool update(const CvMat*& params, CvMat*& J, CvMat*& err);void step();private:double m_lambda_lg10;enum {DONE,START,CALC_J,CHECK_ERR};int m_state;TermCriteria m_criteria;int m_iters;bool m_complete_symm_flag;double m_err_norm, m_prev_err_norm;int m_solver_method;const double m_log10 = log(10.);Ptr<CvMat> m_mask;Ptr<CvMat> m_params;Ptr<CvMat> m_prev_params;Ptr<CvMat> m_JtJ;Ptr<CvMat> m_JtErr;Ptr<CvMat> m_J;Ptr<CvMat> m_Err;Ptr<CvMat> m_avl_JtJ;Ptr<CvMat> m_avl_JtErr;Ptr<CvMat> m_avl_params; };

cpp

#include "LevMarq.h"static void subMatrix(const cv::Mat& src, cv::Mat& dst, const std::vector<uchar>& cols,const std::vector<uchar>& rows) {int nonzeros_cols = cv::countNonZero(cols);cv::Mat tmp(src.rows, nonzeros_cols, CV_64FC1);for (int i = 0, j = 0; i < (int)cols.size(); i++){if (cols[i]){src.col(i).copyTo(tmp.col(j++));}}int nonzeros_rows = cv::countNonZero(rows);dst.create(nonzeros_rows, nonzeros_cols, CV_64FC1);for (int i = 0, j = 0; i < (int)rows.size(); i++){if (rows[i]){tmp.row(i).copyTo(dst.row(j++));}} }LevMarq::LevMarq() {// only do some initm_lambda_lg10 = 0;m_state = DONE;m_criteria = TermCriteria(0, 0, 0);m_iters = 0;m_complete_symm_flag = false;m_err_norm = m_prev_err_norm = DBL_MAX;m_solver_method = DECOMP_SVD; }LevMarq::LevMarq(int nparams, int nerrs, TermCriteria criteria, bool completeSymmFlag) {init(nparams, nerrs, criteria, completeSymmFlag); }LevMarq::~LevMarq() {clear(); }void LevMarq::clear() {m_mask.release();m_params.release();m_prev_params.release();m_JtJ.release();m_JtErr.release();m_J.release();m_Err.release();m_avl_JtJ.release();m_avl_JtErr.release();m_avl_params.release(); }void LevMarq::initParam(Mat params) {if (params.empty())return;double* data = m_params->data.db;assert(params.cols == 1, "params dim must be N*1");for (int i = 0; i < params.rows; i++){data[i] = params.at<double>(i, 0);} }/// <summary> /// create some mat used for store related /// </summary> /// <param name="nparams">the nums of parameters to be optimized</param> /// <param name="nerrs">the nums of residual function</param> /// <param name="criteria">condition of convergence </param> /// <param name="completeSymmFlag"></param> void LevMarq::init(int nparams, int nerrs, TermCriteria criteria, bool completeSymmFlag) {if (!m_params || m_params->rows != nparams || nerrs != (m_Err ? m_Err->rows : 0))clear();m_mask.reset(cvCreateMat(nparams, 1, CV_8U));cvSet(m_mask, cvScalarAll(1));m_prev_params.reset(cvCreateMat(nparams, 1, CV_64F));m_params.reset(cvCreateMat(nparams, 1, CV_64F));m_JtJ.reset(cvCreateMat(nparams, nparams, CV_64F));m_JtErr.reset(cvCreateMat(nparams, 1, CV_64F));if (nerrs > 0){m_J.reset(cvCreateMat(nerrs, nparams, CV_64F));m_Err.reset(cvCreateMat(nerrs, 1, CV_64F));}m_err_norm = m_prev_err_norm = DBL_MAX;m_lambda_lg10 = -3;m_criteria = criteria;if (m_criteria.type & TermCriteria::COUNT){m_criteria.maxCount = MIN(MAX(m_criteria.maxCount, 1), 100000);}else{m_criteria.maxCount = 30;}if (m_criteria.type & TermCriteria::EPS){m_criteria.epsilon = MAX(m_criteria.epsilon, 0);}else{m_criteria.epsilon = DBL_EPSILON;}m_state = START;m_iters = 0;m_complete_symm_flag = completeSymmFlag;m_solver_method = DECOMP_SVD; }bool LevMarq::update(const CvMat*& params, CvMat*& J, CvMat*& err) {J = err = 0;assert(!m_Err.empty());if (m_state == DONE){params = m_params;return false;}if (m_state == START){params = m_params;cvZero(m_J);cvZero(m_Err);J = m_J;err = m_Err;m_state = CALC_J;return true;}if (m_state == CALC_J){cvMulTransposed(m_J, m_JtJ, 1);cvGEMM(m_J, m_Err, 1, 0, 0, m_JtErr, CV_GEMM_A_T);cvCopy(m_params, m_prev_params);step();if (m_iters == 0){m_prev_err_norm = cvNorm(m_Err, 0, CV_L2);}params = m_params;cvZero(m_Err);err = m_Err;m_state = CHECK_ERR;return true;}assert(m_state == CHECK_ERR);m_err_norm = cvNorm(m_Err, 0, CV_L2);if (m_err_norm > m_prev_err_norm){if (++m_lambda_lg10 <= 16){step();params = m_params;cvZero(m_Err);err = m_Err;m_state = CHECK_ERR;return true;}}m_lambda_lg10 = MAX(m_lambda_lg10 - 1, -16);if (++m_iters >= m_criteria.maxCount ||cvNorm(m_params, m_prev_params, CV_RELATIVE_L2) <= m_criteria.epsilon){std::cout << "criteria.epsilon: " << cvNorm(m_params, m_prev_params, CV_RELATIVE_L2) << std::endl;params = m_params;m_state = DONE;return true;}params = m_params;m_prev_err_norm = m_err_norm;cvZero(m_J);cvZero(m_Err);J = m_J;err = m_Err;m_state = CALC_J;return true; }void LevMarq::step() {double miu = exp(m_lambda_lg10 * m_log10);int _nparams = m_params->rows;Mat _JtJ = cvarrToMat(m_JtJ);Mat _mask = cvarrToMat(m_mask);int _avl_nparams = countNonZero(_mask);if (!m_avl_JtJ || m_avl_JtJ->rows != _avl_nparams){m_avl_JtJ.reset(cvCreateMat(_avl_nparams, _avl_nparams, CV_64F));m_avl_JtErr.reset(cvCreateMat(_avl_nparams, 1, CV_64F));m_avl_params.reset(cvCreateMat(_avl_nparams, 1, CV_64F));}Mat _avl_JtJ = cvarrToMat(m_avl_JtJ);Mat _avl_JtErr = cvarrToMat(m_avl_JtErr);Mat_<double> _avl_params = cvarrToMat(m_avl_params);subMatrix(cvarrToMat(m_JtErr), _avl_JtErr, std::vector<uchar>(1, 1), _mask);subMatrix(_JtJ, _avl_JtJ, _mask, _mask);if (!m_Err)completeSymm(_avl_JtErr, m_complete_symm_flag);_avl_JtJ.diag() *= 1.0 + miu;std::cout << miu << std::endl;solve(_avl_JtJ, _avl_JtErr, _avl_params, m_solver_method);int j = 0;for (int i = 0; i < _nparams; i++){m_params->data.db[i] = m_prev_params->data.db[i] - (m_mask->data.ptr[i] ? _avl_params(j++) : 0);} }

算法調(diào)用

double param[opt_num] = { 0,0 }; // k1 p1 CvMat _param = cvMat(opt_num, 1, CV_64F, param); CvTermCriteria criteria = cvTermCriteria(CV_TERMCRIT_EPS + CV_TERMCRIT_ITER, 10000, DBL_EPSILON); myCvLevMarq solver(opt_num, total, criteria); solver.initParam(distCoffes.rowRange(Range(0,2)));int cnt = 0; for (;;) {const CvMat* __param = 0;CvMat* J = 0, * err = 0;bool proceed = solver.update(__param, J, err);cvCopy(__param, &_param); // 將優(yōu)化的變量保存下來(lái)if (!proceed || !err)break;if (J){// 計(jì)算雅可比矩陣和殘差矩陣}else{//殘差矩陣} }

總結(jié)

以上是生活随笔為你收集整理的LM优化算法的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。

算法
LM