Describes implementation of non-linear optimization using the Levenberg-Marquardt algorithm by deriving from the base class LevenbergMarquardtBase instead of using the interface from Minpack directly. , Minpack
#include <bias_config.h>
#include <Base/Debug/Debug.hh>
#include <Base/Debug/TimeMeasure.hh>
#include <Base/Math/Vector2.hh>
#include <MathAlgo/LevenbergMarquardtBase.hh>
#include <math.h>
using namespace BIAS;
using namespace std;
namespace BIAS
{
{
public:
ExampleLevenbergMarquardt()
{
NewDebugLevel("ExampleLM_EVAL_FUN");
NewDebugLevel("ExampleLM_EVAL_JAC");
}
~ExampleLevenbergMarquardt() {}
const double a0, const double b0, double &a, double &b,
bool approximateJacobian = false)
{
data_ = data;
numData_ = data.size();
paramsStart[0] = a0;
paramsStart[1] = b0;
cout << "- Levenberg-Marquardt algorithm starts at (a0,b0) = "
<< paramsStart << endl;
int res = -1;
if (!approximateJacobian) {
res = LM_Compute(numData_, 2, paramsStart, paramsResult, fResult);
} else {
res = LM_ComputeWithoutJacobian(numData_, 2, paramsStart,
paramsResult, fResult);
}
if (res == 0)
{
a = paramsResult[0];
b = paramsResult[1];
cout << "- Levenberg-Marquardt algorithm finished successfully." << endl;
cout << "- Solution is (a,b) = " << paramsResult
<< " with final target function value f(a,b) = " << fResult
<< " and sum of squared residuals |f(a,b)|^2 = "
<< LM_GetSumOfSquaredErrors() << " after " << LM_GetNumOfIterations()
<< " iterations (max. " << LM_GetMaxIterations() << ")" << endl << endl;
} else {
cout << "- Levenberg-Marquardt algorithm failed (return code " << res
<< ")!" << endl << endl;
}
return res;
}
protected:
int numData_;
std::vector< BIAS::Vector2<double> > data_;
{
double a = params[0];
double b = params[1];
double x, y;
for (int i = 0; i < numData_; i++) {
x = data_[i][0];
y = data_[i][1];
residuals[i] = a*cos(b*x) + b*sin(a*x) - y;
}
if (DebugLevelIsSet("ExampleLM_EVAL_FUN"))
cout << "--- " << LM_GetNumOfIterations() << " ----------------------------"
<< "----------------------------------" << endl << endl
<< " (a,b) = " << params << ", f(a,b) = " << residuals
<<
" |f(a,b)| = " << residuals.
NormL2() << endl << endl << flush;
return 0;
}
{
double a = params[0];
double b = params[1];
double x;
for (int i = 0; i < numData_; i++) {
x = data_[i][0];
J[i][0] = cos(b*x) + x*b*cos(a*x);
J[i][1] = -x*a*sin(b*x) + sin(a*x);
}
if (DebugLevelIsSet("ExampleLM_EVAL_JAC"))
cout << " (a,b) = " << params << ", J(a,b) = " << J << endl << flush;
return 0;
}
};
}
int main()
{
BIAS::ExampleLevenbergMarquardt LM;
LM.LM_SetTolerance(1e-12);
LM.LM_SetMaxIterations(1000);
double aGT = 100, bGT = 100;
const int numData = 8;
std::vector< BIAS::Vector2<double> > data(numData);
double x, y;
for (int i = 0; i < numData; i++) {
x = i / (double) numData * 1.8 * M_PI / bGT;
y = aGT*cos(bGT*x) + bGT*sin(aGT*x);
data[i].Set(x, y);
}
cout << endl << "--------------------------" << endl
<< " ExampleLevenbergMarquardt" << endl
<< "--------------------------" << endl
<< "- Compute parameters (a,b) minimizing the sum of squared residuals" << endl
<< " of a*cos(b*x) + b*sin(a*x) - y from " << numData << " measurements (x,y)"
<< endl << endl << "- Tolerance is " << LM.LM_GetTolerance() << endl
<< "- Max. number of iterations is " << LM.LM_GetMaxIterations() << endl << endl;
cout << "Computation with analytically given Jacobian" << endl
<< "----------------------------------------------" << endl << endl;
double a0 = 66, b0 = 110;
double aRes, bRes;
int res = LM.Compute(data, a0, b0, aRes, bRes);
if (res == 0) {
cout << "-> final Jacobian J(a,b) is " << LM.LM_GetJacobianResult() << endl;
}
if (res == 0) {
cout << "-> estimation error is for parameter a = " << fabs(aRes-aGT)
<< " (" << 100 * (aRes != aGT ? fabs(aRes-aGT)/aGT : 0) << " %)"
<< endl << " and for parameter b = " << fabs(bRes-bGT)
<< " (" << 100 * (bRes != bGT ? fabs(bRes-bGT)/bGT : 0) << " %)"
<< endl << endl;
}
cout << "Computation with numerically approximated Jacobian" << endl
<< "----------------------------------------------------" << endl << endl;
res = LM.Compute(data, a0, b0, aRes, bRes, true);
if (res == 0) {
cout << "-> numerical approximation of final Jacobian J(a,b) is "
<< LM.LM_GetJacobianResult() << endl;
}
if (res == 0) {
cout << "-> estimation error with approximated Jacobian"
<< " is for parameter a = " << fabs(aRes-aGT)
<< " (" << 100 * (aRes != aGT ? fabs(aRes-aGT)/aGT : 0) << " %)"
<< endl << " and for parameter b = " << fabs(bRes-bGT)
<< " (" << 100 * (bRes != bGT ? fabs(bRes-bGT)/bGT : 0) << " %)"
<< endl << endl;
}
return (res != 0 ? -1 : 0);
}