How to use cv::LMSolver to solve a nonlinear optimization problem?

C++
1

I want to solve a normal distribution fitting problem, I found cv::LMSolver maybe help me, but I don’t how to use it. I can’t find any example about cv::LMSolver.

I have a series of data(x, y) that satisfies the normal distribution, I want to know The parameter μ and standard deviation σ.
Here’s my code, but it doesn’t work. I don’t know what’s wrong.

#include <opencv2/calib3d.hpp>

struct CvCallbackNDFitting : public cv::LMSolver::Callback
{
	const std::vector<double> xs, ys;
	int npoints;
	CvCallbackNDFitting(const std::vector<double>& x_data, const std::vector<double>& y_data) : xs(x_data), ys(y_data)
	{
		npoints = x_data.size();
	}
	bool compute(cv::InputArray params, cv::OutputArray residuals, cv::OutputArray jacobians) const
	{
		double mu = *(params.getMat().ptr<double>());
		double sigma = *(params.getMat().ptr<double>() + 1);

		if (residuals.empty()) residuals.create(npoints, 1, CV_64F);
		if (jacobians.needed() && jacobians.empty()) jacobians.create(npoints, params.rows(), CV_64F);

		cv::Mat err = residuals.getMat();
		for (size_t i = 0; i < npoints; i++)
		{
			err.at<double>(i, 0) = abs(1.0 / sqrt(2.0 * CV_PI) / sigma * exp(-pow(xs[i] - mu, 2) / 2.0 / pow(sigma, 2)) - ys[i]);
		}
		return true;
	}
};

int main(int argc, char* argv[])
{
	cv::TickMeter tm;
	cv::RNG rng((uint64)time(NULL));
	double mu = 25.0, sigma = 1.0;

	std::vector<double> x_data, y_data;
	for (size_t i = 0; i < 50; i++)
	{
		x_data.push_back(i - 25.0);
		y_data.push_back(1.0 / sqrt(2.0 * CV_PI) / sigma * exp(-pow(x_data[i] - mu, 2) / 2 / pow(sigma, 2)) + rng.gaussian(0.01));
	}

	//LMSolver
	tm.start();
	cv::Mat_<double> miu_sigma = (cv::Mat_<double>(2, 1) << 25.0, 1.0);
	int nIter = cv::LMSolver::create(cv::makePtr<CvCallbackNDFitting>(x_data, y_data), 50)->run(miu_sigma);
	tm.stop();
	
	std::cout << "μ: " << *(miu_sigma[0]) << ", σ: " << *(miu_sigma[1]) << std::endl;
	std::cout << "opencv LMSolver: " << tm.getAvgTimeMilli() << " ms, Iteration: " << nIter << std::endl;
	
	return 0;
}
2

we can only help you, if you tell us the exact problem

mu and sigma are not defined here

and btw, you never calculate any jacobians ?

would a simple meanStdDev() already do it ?

3

jacobians is must be calculate? I don’t need jacobians.

I’m missing a line of code, I have re edited my main post.
double mu = 25.0, sigma = 1.0;

calculate mean and stddev does not seem to be the right way to solve the problem.

1 Like
4

so, your callback error does not change, and it terminates after the 1st iteration.

i’m quite sure you do need valid jacobians, if jacobians.needed()==true
currently they have correct size (2x50), but uninitialized memory, wrecking e.g. this:

5

Oh! thank you, I try to calculate jacobians, the cv::LMSolver is working fine!

struct CvCallbackNDFitting : public cv::LMSolver::Callback
{
	const std::vector<double> xs, ys;
	int npoints;
	CvCallbackNDFitting(const std::vector<double>& x_data, const std::vector<double>& y_data) : xs(x_data), ys(y_data)
	{
		npoints = x_data.size();
	}
	double gaussian(double x, double mu, double sigma) const
	{
		return 1.0 / sqrt(2.0 * CV_PI) / sigma * exp(-pow(x - mu, 2) / 2.0 / pow(sigma, 2));
	}
	bool compute(cv::InputArray params, cv::OutputArray residuals, cv::OutputArray jacobians) const
	{
		//1.ExtractInput
		double mu = *(params.getMat().ptr<double>());
		double sigma = *(params.getMat().ptr<double>() + 1);

		//2.CreatJacAndErr
		if (residuals.empty()) residuals.create(npoints, 1, CV_64F);
		if (jacobians.needed() && jacobians.empty()) jacobians.create(npoints, params.rows(), CV_64F);

		//3.CalcJacAndErr
		cv::Mat err = residuals.getMat();
		cv::Mat jac = jacobians.getMat();
		double eps = 1e-5;
		for (size_t i = 0; i < npoints; i++)
		{
			err.at<double>(i, 0) = abs(gaussian(xs[i], mu, sigma) - ys[i]);

			if (jac.data)
			{
				//partial derivative
				jac.at<double>(i, 0) = abs(gaussian(xs[i], mu + eps, sigma) - ys[i]) - abs(gaussian(xs[i], mu - eps, sigma) - ys[i]);
				jac.at<double>(i, 1) = abs(gaussian(xs[i], mu, sigma + eps) - ys[i]) - abs(gaussian(xs[i], mu, sigma - eps) - ys[i]);
			}
		}

		return true;
	}
};
1 Like