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

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) << ", σ: " << *(miu_sigma) << std::endl;
std::cout << "opencv LMSolver: " << tm.getAvgTimeMilli() << " ms, Iteration: " << nIter << std::endl;

return 0;
}
``````

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 ?

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

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:

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