Update nist.cc to better evaluate results.
Ceres beats Minpack and HBN handily.
Change-Id: I7df8a47b753202ed0b53ab128ce48466bf9f8083
diff --git a/examples/nist.cc b/examples/nist.cc
index 440ab5c..7504e51 100644
--- a/examples/nist.cc
+++ b/examples/nist.cc
@@ -28,18 +28,48 @@
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
//
-// NIST non-linear regression problems solved using Ceres.
+// The National Institute of Standards and Technology has released a
+// set of problems to test non-linear least squares solvers.
//
-// The data was obtained from
-// http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml, where more
-// background on these problems can also be found.
+// More information about the background on these problems and
+// suggested evaluation methodology can be found at:
//
-// Currently not all problems are solved successfully. Some of the
-// failures are due to convergence to a local minimum, and some fail
-// because of numerical issues.
+// http://www.itl.nist.gov/div898/strd/nls/nls_info.shtml
//
-// TODO(sameeragarwal): Fix numerical issues so that all the problems
-// converge and then look at convergence to the wrong solution issues.
+// The problem data themselves can be found at
+//
+// http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml
+//
+// The problems are divided into three levels of difficulty, Easy,
+// Medium and Hard. For each problem there are two starting guesses,
+// the first one far away from the global minimum and the second
+// closer to it.
+//
+// A problem is considered successfully solved, if every components of
+// the solution matches the globally optimal solution in at least 4
+// digits or more.
+//
+// This dataset was used for an evaluation of Non-linear least squares
+// solvers:
+//
+// P. F. Mondragon & B. Borchers, A Comparison of Nonlinear Regression
+// Codes, Journal of Modern Applied Statistical Methods, 4(1):343-351,
+// 2005.
+//
+// The results from Mondragon & Borchers can be summarized as
+// Excel Gnuplot GaussFit HBN MinPack
+// Average LRE 2.3 4.3 4.0 6.8 4.4
+// Winner 1 5 12 29 12
+//
+// Where the row Winner counts, the number of problems for which the
+// solver had the highest LRE.
+
+// In this file, we implement the same evaluation methodology using
+// Ceres. Currently using Levenberg-Marquard with DENSE_QR, we get
+//
+// Excel Gnuplot GaussFit HBN MinPack Ceres
+// Average LRE 2.3 4.3 4.0 6.8 4.4 9.4
+// Winner 0 0 5 11 2 41
#include <iostream>
#include <fstream>
@@ -347,11 +377,12 @@
Matrix predictor = nist_problem.predictor();
Matrix response = nist_problem.response();
Matrix final_parameters = nist_problem.final_parameters();
- std::vector<ceres::Solver::Summary> summaries(nist_problem.num_starts() + 1);
- std::cerr << filename << std::endl;
+
+ printf("%s\n", filename.c_str());
// Each NIST problem comes with multiple starting points, so we
// construct the problem from scratch for each case and solve it.
+ int num_success = 0;
for (int start = 0; start < nist_problem.num_starts(); ++start) {
Matrix initial_parameters = nist_problem.initial_parameters(start);
@@ -365,39 +396,41 @@
initial_parameters.data());
}
- Solve(options, &problem, &summaries[start]);
- }
+ ceres::Solver::Summary summary;
+ Solve(options, &problem, &summary);
- const double certified_cost = nist_problem.certified_cost();
-
- int num_success = 0;
- const int kMinNumMatchingDigits = 4;
- for (int start = 0; start < nist_problem.num_starts(); ++start) {
- const ceres::Solver::Summary& summary = summaries[start];
-
- int num_matching_digits = 0;
- if (IsSuccessfulTermination(summary.termination_type)
- && summary.final_cost < certified_cost) {
- num_matching_digits = kMinNumMatchingDigits + 1;
- } else {
- num_matching_digits =
- -std::log10(fabs(summary.final_cost - certified_cost) / certified_cost);
+ // Compute the LRE by comparing each component of the solution
+ // with the ground truth, and taking the minimum.
+ Matrix final_parameters = nist_problem.final_parameters();
+ const double kMaxNumSignificantDigits = 11;
+ double log_relative_error = kMaxNumSignificantDigits + 1;
+ for (int i = 0; i < num_parameters; ++i) {
+ const double tmp_lre =
+ -std::log10(std::fabs(final_parameters(i) - initial_parameters(i)) /
+ std::fabs(final_parameters(i)));
+ // The maximum LRE is capped at 11 - the precision at which the
+ // ground truth is known.
+ //
+ // The minimum LRE is capped at 0 - no digits match between the
+ // computed solution and the ground truth.
+ log_relative_error =
+ std::min(log_relative_error,
+ std::max(0.0, std::min(kMaxNumSignificantDigits, tmp_lre)));
}
- std::cerr << "start " << start + 1 << " " ;
- if (num_matching_digits <= kMinNumMatchingDigits) {
- std::cerr << "FAILURE";
- } else {
- std::cerr << "SUCCESS";
+ const int kMinNumMatchingDigits = 4;
+ if (log_relative_error >= kMinNumMatchingDigits) {
++num_success;
}
- std::cerr << " summary: "
- << summary.BriefReport()
- << " Certified cost: " << certified_cost
- << std::endl;
+ printf("start: %d status: %s lre: %4.1f initial cost: %e final cost:%e certified cost: %e\n",
+ start + 1,
+ log_relative_error < kMinNumMatchingDigits ? "FAILURE" : "SUCCESS",
+ log_relative_error,
+ summary.initial_cost,
+ summary.final_cost,
+ nist_problem.certified_cost());
}
-
return num_success;
}
@@ -427,7 +460,7 @@
ceres::Solver::Options options;
SetMinimizerOptions(&options);
- std::cerr << "Lower Difficulty\n";
+ std::cout << "Lower Difficulty\n";
int easy_success = 0;
easy_success += RegressionDriver<Misra1a, 1, 2>("Misra1a.dat", options);
easy_success += RegressionDriver<Chwirut, 1, 3>("Chwirut1.dat", options);
@@ -438,7 +471,7 @@
easy_success += RegressionDriver<DanWood, 1, 2>("DanWood.dat", options);
easy_success += RegressionDriver<Misra1b, 1, 2>("Misra1b.dat", options);
- std::cerr << "\nMedium Difficulty\n";
+ std::cout << "\nMedium Difficulty\n";
int medium_success = 0;
medium_success += RegressionDriver<Kirby2, 1, 5>("Kirby2.dat", options);
medium_success += RegressionDriver<Hahn1, 1, 7>("Hahn1.dat", options);
@@ -452,7 +485,7 @@
medium_success += RegressionDriver<Roszman1, 1, 4>("Roszman1.dat", options);
medium_success += RegressionDriver<ENSO, 1, 9>("ENSO.dat", options);
- std::cerr << "\nHigher Difficulty\n";
+ std::cout << "\nHigher Difficulty\n";
int hard_success = 0;
hard_success += RegressionDriver<MGH09, 1, 4>("MGH09.dat", options);
hard_success += RegressionDriver<Thurber, 1, 7>("Thurber.dat", options);
@@ -464,11 +497,11 @@
hard_success += RegressionDriver<Rat43, 1, 4>("Rat43.dat", options);
hard_success += RegressionDriver<Bennet5, 1, 3>("Bennett5.dat", options);
- std::cerr << "\n";
- std::cerr << "Easy : " << easy_success << "/16\n";
- std::cerr << "Medium : " << medium_success << "/22\n";
- std::cerr << "Hard : " << hard_success << "/16\n";
- std::cerr << "Total : " << easy_success + medium_success + hard_success << "/54\n";
+ std::cout << "\n";
+ std::cout << "Easy : " << easy_success << "/16\n";
+ std::cout << "Medium : " << medium_success << "/22\n";
+ std::cout << "Hard : " << hard_success << "/16\n";
+ std::cout << "Total : " << easy_success + medium_success + hard_success << "/54\n";
}
int main(int argc, char** argv) {
