Python實現(xiàn)LM算法的示例代碼

更新時間：2023年06月14日 11:45:59 作者：逃逸的卡路里

L-M方法全稱Levenberg-Marquardt方法，是一種非線性最小二乘優(yōu)化算法，這篇文章整理了該算法的Python和C++實現(xiàn)方法，需要的可以參考一下

L-M方法全稱Levenberg-Marquardt方法，是一種非線性最小二乘優(yōu)化算法，它通過同時利用高斯牛頓方法和梯度下降方法來解決非線性最小二乘問題。其核心思想是在每次迭代中，根據(jù)當前參數(shù)估計計算目標函數(shù)的梯度和海森矩陣，并使用這些信息來更新模型參數(shù)。

LM算法適用于解決各種非線性最小二乘問題，例如數(shù)據(jù)擬合、無約束非線性優(yōu)化等。LM算法相對于其他算法的優(yōu)勢在于其自適應調整步長的能力，使得模型更容易收斂到最優(yōu)解，并且可以避免梯度爆炸或消失的問題。

最近公司一直在研究大氣污染物的檢測技術，對采集數(shù)據(jù)后要進行擬合算法處理，試了很多，在機械機構和電子電路上也做了很多嘗試，昨天對這個LM算法也在軟件上測試了下，被逼無奈的，哈哈哈。其實數(shù)學原理，我是一竅不通。

以下是F = A * exp(t * B) 的LM算法的代碼解釋：

python代碼實現(xiàn)LM算法

python
import numpy as np
from scipy.optimize import leastsq
# 定義目標函數(shù)
def func(x, p):
    A, B = p
    return A * np.exp( * )
# 定義殘差函數(shù)
def residuals(p, y, t):
    return y - func(t, p)
# 初始化參數(shù)
p0 = [1, 1]
t = np.linspace(0, 1, 10)
y = func(t, [2, -1]) + 0.2 * np.random.randn(len(t))
# 使用LM算法擬合模型
plsq = leastsq(residuals, p0, args=(y, t), ftol=1e-15, xtol=1e-15)
# 輸出擬合結果
print(plsq[0])  # [1.99170234, -1.01074562]

上述代碼中，我們首先定義了目標函數(shù) func 和殘差函數(shù) residuals，然后初始化了待擬合的數(shù)據(jù) t 和 y，以及初始參數(shù)值 p0。最后，我們使用 leastsq 函數(shù)對模型進行擬合，并將擬合結果輸出到控制臺中。

在使用LM算法時，需要預設收斂精度(ftol和xtol等參數(shù)），以及控制最大迭代次數(shù)的限制來保證計算效率和避免出現(xiàn)意外循環(huán)情況。

C++ 代碼實現(xiàn)LM算法

#include <cmath>
#include <Eigen/Dense>
#include <unsupported/Eigen/NonLinearOptimization>
using namespace Eigen;
// 定義目標函數(shù)
struct Func {
  int operator()(const VectorXd& x, VectorXd& fvec) const {
    double A = x[0], B = x[1];
    for (int i = 0; i < fvec.size(); ++i) {
      fvec[i] = A * exp(x[2] * i) - B;
    }
    return 0;
  }
};
// 初始化參數(shù)和數(shù)據(jù)
int main() {
  VectorXd x(3);
  x[0] = 2.0;
  x[1] = 1.0;
  x[2] = -1.0;
  VectorXd y(10);
  for (int i = 0; i < y.size(); ++i + 0.2 * std::rand()/RAND_MAX;
  }
  // 定義 LM 問題并求解
  NumericalDiff<Func> numerical_diff;
  LevenbergMarquardt<NumericalDiff<Func>> lm(numerical_diff);
  lm.parameters.ftol = 1e-15;
  lm.parameters.xtol = 1e-15;
  lm.parameters.maxfev = 1000;
  lm.minimize(x, y);
  // 輸出擬合結果
  std::cout << "A = " << x[0] << ", B = " << x[1] << ", t = " << x[2] << std::endl;
  return 0;
}

以上代碼中，我們使用 Eigen 庫來計算矩陣的逆和乘積，使用 NumericalDiff 類計算函數(shù)的一階導數(shù)，使用 LevenbergMarquardt 類求解 LM 問題。其中，minimize 函數(shù)接收待擬合數(shù)據(jù) y 和初始參數(shù)向量 x 作為輸入，并可以自動調整 LM 算法的縮放因子 λ 和 μ，同時支持設置迭代次數(shù)限制、收斂精度等參數(shù)來控制算法行為。

以上兩種是實現(xiàn)方法是人工智能生成的，有興趣的小伙伴可以測試跑跑看。

自己測試過得代碼實現(xiàn)LM算法

#include "stdafx.h"
#include <cstdio>
#include <vector>
#include <opencv2/core/core.hpp>
#include <fstream>
#include <string>
using namespace std;
using namespace cv;
const double DERIV_STEP = 1e-4; // 擬合精度
const int MAX_ITER = 100; // 最大循環(huán)次數(shù)
void LM(double(*Func)(const Mat &input, const Mat params), // function pointer
	const Mat &inputs, const Mat &outputs, Mat& params);
double Deriv(double(*Func)(const Mat &input, const Mat params), // function pointer
	const Mat &input, const Mat params, int n);
// The user defines their function here
double Func(const Mat &input, const Mat params);
int main()
{
	// For this demo we're going to try and fit to the function
	// F = A*exp(t*B)
	// There are 2 parameters: A B
	int num_params = 2;
	// Generate random data using these parameters
	int total_data = 1410;
	Mat inputs(total_data, 1, CV_64F);
	Mat outputs(total_data, 1, CV_64F);
// 	//load observation data
// 	for (int i = 0; i < total_data; i++) {
// 		inputs.at<double>(i, 0) = i + 1;  //load year
// 	}
// 	//load America population
// 	outputs.at<double>(0, 0) = 8.3;
// 	outputs.at<double>(1, 0) = 11.0;
// 	outputs.at<double>(2, 0) = 14.7;
// 	outputs.at<double>(3, 0) = 19.7;
// 	outputs.at<double>(4, 0) = 26.7;
// 	outputs.at<double>(5, 0) = 35.2;
// 	outputs.at<double>(6, 0) = 44.4;
// 	outputs.at<double>(7, 0) = 55.9;
	for (int i = 0; i < total_data; i++) 
	{
		inputs.at<double>(i, 0) = i * 0.05; 
	}
	ifstream ifstxt;
	ifstxt.open("shuju.txt");
	string strline;
	int iout = 0;
	while (getline(ifstxt,strline))
	{
		outputs.at<double>(iout, 0) = stoi(strline.c_str());
		iout++;
	}
	ifstxt.close();
	///
	// Guess the parameters, it should be close to the true value, else it can fail for very sensitive functions!
	Mat params(num_params, 1, CV_64F);
	//init guess
	params.at<double>(0, 0) = 14000;
	params.at<double>(1, 0) = -0.05;
	LM(Func, inputs, outputs, params);
	printf("Parameters from LM: %lf %lf\n", params.at<double>(0, 0), params.at<double>(1, 0));
	system("pause");
	return 0;
}
double Func(const Mat &input, const Mat params)
{
	// Assumes input is a single row matrix
	// Assumes params is a column matrix
	double A = params.at<double>(0, 0);
	double B = params.at<double>(1, 0);
	double x = input.at<double>(0, 0);
	return A*exp(x*B);
}
//calc the n-th params' partial derivation ， the params are our  final target
double Deriv(double(*Func)(const Mat &input, const Mat params), const Mat &input, const Mat params, int n)
{
	// Assumes input is a single row matrix
	// Returns the derivative of the nth parameter
	Mat params1 = params.clone();
	Mat params2 = params.clone();
	// Use central difference  to get derivative
	params1.at<double>(n, 0) -= DERIV_STEP;
	params2.at<double>(n, 0) += DERIV_STEP;
	double p1 = Func(input, params1);
	double p2 = Func(input, params2);
	double d = (p2 - p1) / (2 * DERIV_STEP);
	return d;
}
void LM(double(*Func)(const Mat &input, const Mat params),
	const Mat &inputs, const Mat &outputs, Mat& params)
{
	int m = inputs.rows;
	int n = inputs.cols;
	int num_params = params.rows;
	Mat r(m, 1, CV_64F); // residual matrix
	Mat r_tmp(m, 1, CV_64F);
	Mat Jf(m, num_params, CV_64F); // Jacobian of Func()
	Mat input(1, n, CV_64F); // single row input
	Mat params_tmp = params.clone();
	double last_mse = 0;
	float u = 1, v = 2;
	Mat I = Mat::ones(num_params, num_params, CV_64F);//construct identity matrix
	int i = 0;
	for (i = 0; i < MAX_ITER; i++)
	{
		double mse = 0;
		double mse_temp = 0;
		for (int j = 0; j < m; j++)
		{
			for (int k = 0; k < n; k++)
			{//copy Independent variable vector, the year
				input.at<double>(0, k) = inputs.at<double>(j, k);
			}
			r.at<double>(j, 0) = outputs.at<double>(j, 0) - Func(input, params);//diff between previous estimate and observation population
			mse += r.at<double>(j, 0)*r.at<double>(j, 0);
			for (int k = 0; k < num_params; k++) {
				Jf.at<double>(j, k) = Deriv(Func, input, params, k);  //construct jacobian matrix
			}
		}
		mse /= m;
		//printf("%lf\n", mse /= m);
		params_tmp = params.clone();
		Mat hlm = (Jf.t()*Jf + u*I).inv()*Jf.t()*r; //calculate deta
		params_tmp += hlm; //update value
		for (int j = 0; j < m; j++) {
			r_tmp.at<double>(j, 0) = outputs.at<double>(j, 0) - Func(input, params_tmp);//diff between current estimate and observation population
			mse_temp += r_tmp.at<double>(j, 0)*r_tmp.at<double>(j, 0);//diff square sum
		}
		mse_temp /= m;//diff square sum
		Mat q(1, 1, CV_64F);
		q = (mse - mse_temp) / (0.5*hlm.t()*(u*hlm - Jf.t()*r));
		double q_value = q.at<double>(0, 0);
		if (q_value > 0)
		{
			double s = 1.0 / 3.0;
			v = 2;
			mse = mse_temp;
			params = params_tmp;
			double temp = 1 - pow(2 * q_value - 1, 3);
			if (s > temp)
			{
				u = u * s;
			}
			else
			{
				u = u * temp;
			}
		}
		else
		{
			u = u*v;
			v = 2 * v;
			params = params_tmp;
		}
		// The difference in mse is very small, so quit
		if (fabs(mse - last_mse) < 1e-5)  // 這個值得大小，影響循環(huán)跳出，計算精度（考慮計算時間問題）
		{
			printf("%d %lf\n", i, mse);
			break;
		}
		//printf("%d: mse=%f\n", i, mse);
		//printf("%d %lf\n", i, mse);
		last_mse = mse;
	}
}

在配置#include <opencv2/core/core.hpp>，這個的時候會遇到一些問題，小伙伴可以在網(wǎng)上自己搜下解決辦法，

下面是我再vs2015中設置的截圖：電腦上安裝兩個版本OpenCV，設置有點亂