diff --git a/examples/nist.cc b/examples/nist.cc index 696bd675b..f78d8eb8f 100644 --- a/examples/nist.cc +++ b/examples/nist.cc @@ -264,7 +264,7 @@ NIST_END // y = b1 / ((1+exp[b2-b3*x])**(1/b4)) + e NIST_BEGIN(Rat43) - b[0] / pow(T(1.0) + exp(b[1] - b[2] * x), T(1.0) / b[4]) + b[0] / pow(T(1.0) + exp(b[1] - b[2] * x), T(1.0) / b[3]) NIST_END // y = (b1 + b2*x + b3*x**2 + b4*x**3) / @@ -320,13 +320,14 @@ int RegressionDriver(const std::string& filename, Matrix response = nist_problem.response(); Matrix final_parameters = nist_problem.final_parameters(); std::vector summaries(nist_problem.num_starts() + 1); + std::cerr << filename << std::endl; // Each NIST problem comes with multiple starting points, so we // construct the problem from scratch for each case and solve it. for (int start = 0; start < nist_problem.num_starts(); ++start) { Matrix initial_parameters = nist_problem.initial_parameters(start); - ceres::Problem problem; + ceres::Problem problem; for (int i = 0; i < nist_problem.num_observations(); ++i) { problem.AddResidualBlock( new ceres::AutoDiffCostFunction( @@ -340,44 +341,51 @@ int RegressionDriver(const std::string& filename, } // Ugly hack to get the objective function value at the certified - // optimal parameter values. - Matrix initial_parameters = nist_problem.final_parameters(); - ceres::Problem problem; - for (int i = 0; i < nist_problem.num_observations(); ++i) { - problem.AddResidualBlock( - new ceres::AutoDiffCostFunction( - new Model(predictor.data() + nist_problem.predictor_size() * i, - response.data() + nist_problem.response_size() * i)), - NULL, - initial_parameters.data()); + // optimal parameter values. So we build the problem and call Ceres + // with zero iterations to get the initial_cost. + { + Matrix initial_parameters = nist_problem.final_parameters(); + ceres::Problem problem; + for (int i = 0; i < nist_problem.num_observations(); ++i) { + problem.AddResidualBlock( + new ceres::AutoDiffCostFunction( + new Model(predictor.data() + nist_problem.predictor_size() * i, + response.data() + nist_problem.response_size() * i)), + NULL, + initial_parameters.data()); + } + + ceres::Solver::Options options; + options.max_num_iterations = 0; + Solve(options, &problem, &summaries[nist_problem.num_starts()]); } - Solve(options, &problem, &summaries.back()); - double certified_cost = summaries.back().initial_cost; - std::cout << filename << std::endl; + + double certified_cost = summaries[nist_problem.num_starts()].initial_cost; int num_success = 0; - for (int i = 0; i < nist_problem.num_starts(); ++i) { + for (int start = 0; start < nist_problem.num_starts(); ++start) { + const ceres::Solver::Summary& summary = summaries[start]; const int num_matching_digits = -std::log10(1e-18 + - fabs(summaries[i].final_cost - certified_cost) + fabs(summary.final_cost - certified_cost) / certified_cost); - std::cout << "start " << i + 1 << " " ; + std::cerr << "start " << start + 1 << " " ; if (num_matching_digits > 4) { ++num_success; - std::cout << "SUCCESS"; + std::cerr << "SUCCESS"; } else { - std::cout << "FAILURE"; + std::cerr << "FAILURE"; } - std::cout << " digits: " << num_matching_digits; - std::cout << " termination: " - << ceres::SolverTerminationTypeToString(summaries[i].termination_type) + std::cerr << " digits: " << num_matching_digits; + std::cerr << " summary: " + << summary.BriefReport() << std::endl; } return num_success; } void SolveNISTProblems(const ceres::Solver::Options& options) { - std::cout << "Lower Difficulty\n"; + std::cerr << "Lower Difficulty\n"; int easy_success = 0; easy_success += RegressionDriver("Misra1a.dat", options); easy_success += RegressionDriver("Chwirut1.dat", options); @@ -388,7 +396,7 @@ void SolveNISTProblems(const ceres::Solver::Options& options) { easy_success += RegressionDriver("DanWood.dat", options); easy_success += RegressionDriver("Misra1b.dat", options); - std::cout << "\nMedium Difficulty\n"; + std::cerr << "\nMedium Difficulty\n"; int medium_success = 0; medium_success += RegressionDriver("Kirby2.dat", options); medium_success += RegressionDriver("Hahn1.dat", options); @@ -402,22 +410,23 @@ void SolveNISTProblems(const ceres::Solver::Options& options) { medium_success += RegressionDriver("Roszman1.dat", options); medium_success += RegressionDriver("ENSO.dat", options); - std::cout << "\nHigher Difficulty\n"; + std::cerr << "\nHigher Difficulty\n"; int hard_success = 0; hard_success += RegressionDriver("MGH09.dat", options); hard_success += RegressionDriver("Thurber.dat", options); hard_success += RegressionDriver("BoxBOD.dat", options); hard_success += RegressionDriver("Rat42.dat", options); hard_success += RegressionDriver("MGH10.dat", options); + hard_success += RegressionDriver("Eckerle4.dat", options); hard_success += RegressionDriver("Rat43.dat", options); hard_success += RegressionDriver("Bennett5.dat", options); - 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"; + 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"; } int main(int argc, char** argv) { @@ -428,7 +437,7 @@ int main(int argc, char** argv) { // linear solvers. ceres::Solver::Options options; options.linear_solver_type = ceres::DENSE_QR; - options.max_num_iterations = 1000; + options.max_num_iterations = 2000; options.function_tolerance *= 1e-10; options.gradient_tolerance *= 1e-10; options.parameter_tolerance *= 1e-10; diff --git a/internal/ceres/residual_block.cc b/internal/ceres/residual_block.cc index 03867891d..bdb88b1dd 100644 --- a/internal/ceres/residual_block.cc +++ b/internal/ceres/residual_block.cc @@ -102,8 +102,11 @@ bool ResidualBlock::Evaluate(double* cost, InvalidateEvaluation(*this, cost, residuals, eval_jacobians); - if (!cost_function_->Evaluate(parameters.get(), residuals, eval_jacobians) || - !IsEvaluationValid(*this, + if (!cost_function_->Evaluate(parameters.get(), residuals, eval_jacobians)) { + return false; + } + + if (!IsEvaluationValid(*this, parameters.get(), cost, residuals, diff --git a/internal/ceres/trust_region_minimizer.cc b/internal/ceres/trust_region_minimizer.cc index 19e3f3583..76c4f8a75 100644 --- a/internal/ceres/trust_region_minimizer.cc +++ b/internal/ceres/trust_region_minimizer.cc @@ -34,6 +34,7 @@ #include #include #include +#include #include #include @@ -281,9 +282,11 @@ void TrustRegionMinimizer::Minimize(const Minimizer::Options& options, // new_model_cost can actually be larger than half the squared // norm of the residual vector. We allow for small tolerance // around cost and beyond that declare the step to be invalid. - if (cost < (new_model_cost - kEpsilon)) { + if ((1.0 - new_model_cost / cost) < -kEpsilon) { VLOG(1) << "Invalid step: current_cost: " << cost - << " new_model_cost " << new_model_cost; + << " new_model_cost " << new_model_cost + << " absolute difference " << (cost - new_model_cost) + << " relative difference " << (1.0 - new_model_cost/cost); } else { iteration_summary.step_is_valid = true; } @@ -320,12 +323,23 @@ void TrustRegionMinimizer::Minimize(const Minimizer::Options& options, num_consecutive_invalid_steps = 0; // We allow some slop around 0, and clamp the model_cost_change - // at kEpsilon from below. + // at kEpsilon * min(1.0, cost) from below. // - // There is probably a better way to do this, as it is going to - // create problems for problems where the objective function is - // kEpsilon close to zero. - const double model_cost_change = max(kEpsilon, cost - new_model_cost); + // In exact arithmetic this should never be needed, as we are + // guaranteed to new_model_cost <= cost. However, due to various + // numerical issues, it is possible that new_model_cost is + // nearly equal to cost, and the difference is a small negative + // number. To make sure that the relative_decrease computation + // remains sane, as clamp the difference (cost - new_model_cost) + // from below at a small positive number. + // + // This number is the minimum of kEpsilon * (cost, 1.0), which + // ensures that it will never get too large in absolute value, + // while scaling down proportionally with the magnitude of the + // cost. This is important for problems where the minimum of the + // objective function is near zero. + const double model_cost_change = + max(kEpsilon * min(1.0, cost), cost - new_model_cost); // Undo the Jacobian column scaling. delta = (trust_region_step.array() * scale.array()).matrix(); @@ -356,9 +370,11 @@ void TrustRegionMinimizer::Minimize(const Minimizer::Options& options, if (!evaluator->Evaluate(x_plus_delta.data(), &new_cost, NULL, NULL, NULL)) { - summary->termination_type = NUMERICAL_FAILURE; - LOG(WARNING) << "Terminating: Cost evaluation failed."; - return; + // If the evaluation of the new cost fails, treat it as a step + // with high cost. + LOG(WARNING) << "Step failed to evaluate. " + << "Treating it as step with infinite cost"; + new_cost = numeric_limits::max(); } VLOG(2) << "old cost: " << cost << " new cost: " << new_cost;