當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

LM优化算法

發布時間：2023/12/10 编程问答 26 豆豆

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

LM算法

理論知識
- - 梯度下降
  - 高斯牛頓
  - Levenberg–Marquardt
- 算法框架
- - 算法的整體流程
  - 求解器update流程
  - 說明
- 算法實現
- - - 頭文件
    - cpp
- 算法調用

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

本人非數學專業，且對算法理解可能不到位，詳細的算法推導及各個優化算法之間的關系，非常推薦看METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS ，其介紹更詳細也更專業。

以下簡要介紹LM算法，然后給出opencv中有關的實現。

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

當然，有關LM算法的實現，在opencv的usac模塊中有更為正式的實現，在該模塊中，對許多算法進行了集成和優化，后續再進行研究。

理論知識

對于一個最小二乘問題，有：

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

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

梯度下降

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

高斯牛頓

高斯牛頓是在牛頓迭代的基礎上，用雅可比矩陣的平方來近似Hessian矩陣，求解起來將會更加高效（大多數情況下，Hessian矩陣太難求）。這樣做的缺點就是，要求迭代初值必須在最優解附近，否則以上假設將不成立。這里給出其公式：
$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不可逆，導致高斯牛頓無法使用。LM通過將JtJ加上一個μI（I為單位陣），保證了正規方程的左側一定可逆。

在算法實際的應用中，通過調節μ的大小，可以使算法在梯度下降和高斯牛頓兩者之間切換，如：

使用較大的μ值，算法可以看作梯度下降，適合當前優化位置距離最優解較遠的情況
使用較小的μ值，算法可以看作高斯牛頓，適合當前優化位置接近最優解。

算法框架

LM算法通過求解正規方程，得到每次迭代的步長和方向。因此，算法的主要目的是求解正規方程左右側的項，然后通過SVD分解即可得到參數更新的步長及方向，即：

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

算法框架在外部計算雅可比矩陣和殘差矩陣，然后在算法內部通過求解正規方程，得到每次參數更新的步長及方向。

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

算法的整體流程

求解器update流程

求解器updata()內部根據不同的狀態進行相應的計算，具體流程如下：

說明

J表示雅可比矩陣
err表示殘差矩陣
求解器內外同步表示求解器內部和外部相應的參數指向相同，便于在求解器外部進行雅可比矩陣和殘差矩陣的計算

算法實現

算法的實現主要步驟都體現在上述流程圖中。

頭文件

#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);} }

算法調用

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); // 將優化的變量保存下來if (!proceed || !err)break;if (J){// 計算雅可比矩陣和殘差矩陣}else{//殘差矩陣} }

總結

以上是生活随笔為你收集整理的LM优化算法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。

算法
LM