From c4a329155c66816998a5313adbac052dcf0f0ef7 Mon Sep 17 00:00:00 2001 From: Sameer Agarwal Date: Thu, 13 Jun 2013 22:00:48 -0700 Subject: [PATCH] Enable support for dumping trust region minimizer problems. This support was broken due to the TrustRegionMinimizer refactoring. It is now enabled again, with the responsibilty for dumping the problem shifted to the individual TrustRegionStrategy. There is however one wrinkle, which is perhaps an indication of poor design to start with. The LinearLeastSquaresProblemProto carries in it num_eliminate_blocks, something which does not exist anymore. More importantly, the TrustRegionStrategy does not have access to this quantity anymore. Dealing with this will be the subject of a future change. Change-Id: I358adf6a2e386f4940b617bf950d6c7e87d2635d --- docs/source/solving.rst | 40 ++++++++------ include/ceres/solver.h | 23 ++++---- internal/ceres/dogleg_strategy.cc | 19 ++++++- internal/ceres/dogleg_strategy.h | 6 +- .../ceres/levenberg_marquardt_strategy.cc | 18 +++++- .../ceres/linear_least_squares_problems.cc | 55 +++++++------------ .../ceres/linear_least_squares_problems.h | 3 +- internal/ceres/minimizer.h | 15 +++-- internal/ceres/solver_impl.cc | 10 +++- internal/ceres/trust_region_minimizer.cc | 55 ++++++++----------- internal/ceres/trust_region_strategy.h | 16 ++++++ 11 files changed, 149 insertions(+), 111 deletions(-) diff --git a/docs/source/solving.rst b/docs/source/solving.rst index 80aee8636..d3f2d0faf 100644 --- a/docs/source/solving.rst +++ b/docs/source/solving.rst @@ -1172,29 +1172,32 @@ elimination group [LiSaad]_. #. ``it`` is the time take by the current iteration. #. ``tt`` is the the total time taken by the minimizer. -.. member:: vector Solver::Options::lsqp_iterations_to_dump +.. member:: vector Solver::Options::trust_region_minimizer_iterations_to_dump Default: ``empty`` - List of iterations at which the optimizer should dump the linear - least squares problem to disk. Useful for testing and - benchmarking. If ``empty``, no problems are dumped. + List of iterations at which the trust region minimizer should dump + the trust region problem. Useful for testing and benchmarking. If + ``empty``, no problems are dumped. -.. member:: string Solver::Options::lsqp_dump_directory +.. member:: string Solver::Options::trust_region_problem_dump_directory Default: ``/tmp`` - If :member:`Solver::Options::lsqp_iterations_to_dump` is non-empty, then - this setting determines the directory to which the files containing - the linear least squares problems are written to. + Directory to which the problems should be written to. Should be + non-empty if + :member:`Solver::Options::trust_region_minimizer_iterations_to_dump` is + non-empty and + :member:`Solver::Options::trust_region_problem_dump_format_type` is not + ``CONSOLE``. -.. member:: DumpFormatType Solver::Options::lsqp_dump_format +.. member:: DumpFormatType Solver::Options::trust_region_problem_dump_format Default: ``TEXTFILE`` - The format in which linear least squares problems should be logged - when :member:`Solver::Options::lsqp_iterations_to_dump` is non-empty. - There are three options: + The format in which trust region problems should be logged when + :member:`Solver::Options::trust_region_minimizer_iterations_to_dump` + is non-empty. There are three options: * ``CONSOLE`` prints the linear least squares problem in a human readable format to ``stderr``. The Jacobian is printed as a @@ -1203,22 +1206,25 @@ elimination group [LiSaad]_. problems. * ``PROTOBUF`` Write out the linear least squares problem to the - directory pointed to by :member:`Solver::Options::lsqp_dump_directory` as + directory pointed to by + :member:`Solver::Options::trust_region_problem_dump_directory` as a protocol buffer. ``linear_least_squares_problems.h/cc`` contains routines for loading these problems. For details on the on disk format used, see ``matrix.proto``. The files are named - ``lm_iteration_???.lsqp``. This requires that ``protobuf`` be + ``ceres_solver_iteration_???.bin``. This requires that ``protobuf`` be linked into Ceres Solver. * ``TEXTFILE`` Write out the linear least squares problem to the - directory pointed to by member::`Solver::Options::lsqp_dump_directory` as + directory pointed to by + :member:`Solver::Options::trust_region_problem_dump_directory` as text files which can be read into ``MATLAB/Octave``. The Jacobian is dumped as a text file containing :math:`(i,j,s)` triplets, the vectors :math:`D`, `x` and `f` are dumped as text files containing a list of their values. - A ``MATLAB/Octave`` script called ``lm_iteration_???.m`` is also - output, which can be used to parse and load the problem into memory. + A ``MATLAB/Octave`` script called + ``ceres_solver_iteration_???.m`` is also output, which can be + used to parse and load the problem into memory. .. member:: bool Solver::Options::check_gradients diff --git a/include/ceres/solver.h b/include/ceres/solver.h index 8bf87f9d2..a86a58629 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h @@ -106,8 +106,8 @@ class Solver { inner_iteration_ordering = NULL; logging_type = PER_MINIMIZER_ITERATION; minimizer_progress_to_stdout = false; - lsqp_dump_directory = "/tmp"; - lsqp_dump_format_type = TEXTFILE; + trust_region_problem_dump_directory = "/tmp"; + trust_region_problem_dump_format_type = TEXTFILE; check_gradients = false; gradient_check_relative_precision = 1e-8; numeric_derivative_relative_step_size = 1e-6; @@ -490,14 +490,17 @@ class Solver { // is sent to STDOUT. bool minimizer_progress_to_stdout; - // List of iterations at which the optimizer should dump the - // linear least squares problem to disk. Useful for testing and - // benchmarking. If empty (default), no problems are dumped. - // - // This is ignored if protocol buffers are disabled. - vector lsqp_iterations_to_dump; - string lsqp_dump_directory; - DumpFormatType lsqp_dump_format_type; + // List of iterations at which the minimizer should dump the trust + // region problem. Useful for testing and benchmarking. If empty + // (default), no problems are dumped. + vector trust_region_minimizer_iterations_to_dump; + + // Directory to which the problems should be written to. Should be + // non-empty if trust_region_minimizer_iterations_to_dump is + // non-empty and trust_region_problem_dump_format_type is not + // CONSOLE. + string trust_region_problem_dump_directory; + DumpFormatType trust_region_problem_dump_format_type; // Finite differences options ---------------------------------------------- diff --git a/internal/ceres/dogleg_strategy.cc b/internal/ceres/dogleg_strategy.cc index a330ad2c7..03c0d78c9 100644 --- a/internal/ceres/dogleg_strategy.cc +++ b/internal/ceres/dogleg_strategy.cc @@ -34,6 +34,7 @@ #include "Eigen/Dense" #include "ceres/array_utils.h" #include "ceres/internal/eigen.h" +#include "ceres/linear_least_squares_problems.h" #include "ceres/linear_solver.h" #include "ceres/polynomial.h" #include "ceres/sparse_matrix.h" @@ -127,7 +128,7 @@ TrustRegionStrategy::Summary DoglegStrategy::ComputeStep( ComputeCauchyPoint(jacobian); LinearSolver::Summary linear_solver_summary = - ComputeGaussNewtonStep(jacobian, residuals); + ComputeGaussNewtonStep(per_solve_options, jacobian, residuals); TrustRegionStrategy::Summary summary; summary.residual_norm = linear_solver_summary.residual_norm; @@ -507,6 +508,7 @@ bool DoglegStrategy::FindMinimumOnTrustRegionBoundary(Vector2d* minimum) const { } LinearSolver::Summary DoglegStrategy::ComputeGaussNewtonStep( + const PerSolveOptions& per_solve_options, SparseMatrix* jacobian, const double* residuals) { const int n = jacobian->num_cols(); @@ -561,6 +563,21 @@ LinearSolver::Summary DoglegStrategy::ComputeGaussNewtonStep( solve_options, gauss_newton_step_.data()); + if (per_solve_options.dump_format_type == CONSOLE || + (per_solve_options.dump_format_type != CONSOLE && + !per_solve_options.dump_filename_base.empty())) { + if (!DumpLinearLeastSquaresProblem(per_solve_options.dump_filename_base, + per_solve_options.dump_format_type, + jacobian, + solve_options.D, + residuals, + gauss_newton_step_.data(), + 0)) { + LOG(ERROR) << "Unable to dump trust region problem." + << " Filename base: " << per_solve_options.dump_filename_base; + } + } + if (linear_solver_summary.termination_type == FAILURE || !IsArrayValid(n, gauss_newton_step_.data())) { mu_ *= mu_increase_factor_; diff --git a/internal/ceres/dogleg_strategy.h b/internal/ceres/dogleg_strategy.h index 7131467d6..71c785cc3 100644 --- a/internal/ceres/dogleg_strategy.h +++ b/internal/ceres/dogleg_strategy.h @@ -79,8 +79,10 @@ class DoglegStrategy : public TrustRegionStrategy { typedef Eigen::Matrix Vector2d; typedef Eigen::Matrix Matrix2d; - LinearSolver::Summary ComputeGaussNewtonStep(SparseMatrix* jacobian, - const double* residuals); + LinearSolver::Summary ComputeGaussNewtonStep( + const PerSolveOptions& per_solve_options, + SparseMatrix* jacobian, + const double* residuals); void ComputeCauchyPoint(SparseMatrix* jacobian); void ComputeGradient(SparseMatrix* jacobian, const double* residuals); void ComputeTraditionalDoglegStep(double* step); diff --git a/internal/ceres/levenberg_marquardt_strategy.cc b/internal/ceres/levenberg_marquardt_strategy.cc index 9e6a59e38..6b29244e9 100644 --- a/internal/ceres/levenberg_marquardt_strategy.cc +++ b/internal/ceres/levenberg_marquardt_strategy.cc @@ -34,6 +34,7 @@ #include "Eigen/Core" #include "ceres/array_utils.h" #include "ceres/internal/eigen.h" +#include "ceres/linear_least_squares_problems.h" #include "ceres/linear_solver.h" #include "ceres/sparse_matrix.h" #include "ceres/trust_region_strategy.h" @@ -111,9 +112,24 @@ TrustRegionStrategy::Summary LevenbergMarquardtStrategy::ComputeStep( } else { VectorRef(step, num_parameters) *= -1.0; } - reuse_diagonal_ = true; + if (per_solve_options.dump_format_type == CONSOLE || + (per_solve_options.dump_format_type != CONSOLE && + !per_solve_options.dump_filename_base.empty())) { + if (!DumpLinearLeastSquaresProblem(per_solve_options.dump_filename_base, + per_solve_options.dump_format_type, + jacobian, + solve_options.D, + residuals, + step, + 0)) { + LOG(ERROR) << "Unable to dump trust region problem." + << " Filename base: " << per_solve_options.dump_filename_base; + } + } + + TrustRegionStrategy::Summary summary; summary.residual_norm = linear_solver_summary.residual_norm; summary.num_iterations = linear_solver_summary.num_iterations; diff --git a/internal/ceres/linear_least_squares_problems.cc b/internal/ceres/linear_least_squares_problems.cc index 6c886a1be..d1ee7f0c3 100644 --- a/internal/ceres/linear_least_squares_problems.cc +++ b/internal/ceres/linear_least_squares_problems.cc @@ -574,9 +574,7 @@ LinearLeastSquaresProblem* LinearLeastSquaresProblem3() { } namespace { -bool DumpLinearLeastSquaresProblemToConsole(const string& directory, - int iteration, - const SparseMatrix* A, +bool DumpLinearLeastSquaresProblemToConsole(const SparseMatrix* A, const double* D, const double* b, const double* x, @@ -602,8 +600,7 @@ bool DumpLinearLeastSquaresProblemToConsole(const string& directory, }; #ifndef CERES_NO_PROTOCOL_BUFFERS -bool DumpLinearLeastSquaresProblemToProtocolBuffer(const string& directory, - int iteration, +bool DumpLinearLeastSquaresProblemToProtocolBuffer(const string& filename_base, const SparseMatrix* A, const double* D, const double* b, @@ -632,18 +629,14 @@ bool DumpLinearLeastSquaresProblemToProtocolBuffer(const string& directory, } lsqp.set_num_eliminate_blocks(num_eliminate_blocks); - string format_string = JoinPath(directory, - "lm_iteration_%03d.lsqp"); - string filename = - StringPrintf(format_string.c_str(), iteration); - LOG(INFO) << "Dumping least squares problem for iteration " << iteration - << " to disk. File: " << filename; + + const string filename = filename_base + ".bin"; + LOG(INFO) << "Dumping least squares problem to disk. File: " << filename; WriteStringToFileOrDie(lsqp.SerializeAsString(), filename); return true; } #else -bool DumpLinearLeastSquaresProblemToProtocolBuffer(const string& directory, - int iteration, +bool DumpLinearLeastSquaresProblemToProtocolBuffer(const string& filename_base, const SparseMatrix* A, const double* D, const double* b, @@ -669,31 +662,25 @@ void WriteArrayToFileOrDie(const string& filename, fclose(fptr); } -bool DumpLinearLeastSquaresProblemToTextFile(const string& directory, - int iteration, +bool DumpLinearLeastSquaresProblemToTextFile(const string& filename_base, const SparseMatrix* A, const double* D, const double* b, const double* x, int num_eliminate_blocks) { CHECK_NOTNULL(A); - string format_string = JoinPath(directory, - "lm_iteration_%03d"); - string filename_prefix = - StringPrintf(format_string.c_str(), iteration); - - LOG(INFO) << "writing to: " << filename_prefix << "*"; + LOG(INFO) << "writing to: " << filename_base << "*"; string matlab_script; StringAppendF(&matlab_script, - "function lsqp = lm_iteration_%03d()\n", iteration); + "function lsqp = load_trust_region_problem()\n"); StringAppendF(&matlab_script, "lsqp.num_rows = %d;\n", A->num_rows()); StringAppendF(&matlab_script, "lsqp.num_cols = %d;\n", A->num_cols()); { - string filename = filename_prefix + "_A.txt"; + string filename = filename_base + "_A.txt"; FILE* fptr = fopen(filename.c_str(), "w"); CHECK_NOTNULL(fptr); A->ToTextFile(fptr); @@ -709,34 +696,33 @@ bool DumpLinearLeastSquaresProblemToTextFile(const string& directory, if (D != NULL) { - string filename = filename_prefix + "_D.txt"; + string filename = filename_base + "_D.txt"; WriteArrayToFileOrDie(filename, D, A->num_cols()); StringAppendF(&matlab_script, "lsqp.D = load('%s', '-ascii');\n", filename.c_str()); } if (b != NULL) { - string filename = filename_prefix + "_b.txt"; + string filename = filename_base + "_b.txt"; WriteArrayToFileOrDie(filename, b, A->num_rows()); StringAppendF(&matlab_script, "lsqp.b = load('%s', '-ascii');\n", filename.c_str()); } if (x != NULL) { - string filename = filename_prefix + "_x.txt"; + string filename = filename_base + "_x.txt"; WriteArrayToFileOrDie(filename, x, A->num_cols()); StringAppendF(&matlab_script, "lsqp.x = load('%s', '-ascii');\n", filename.c_str()); } - string matlab_filename = filename_prefix + ".m"; + string matlab_filename = filename_base + ".m"; WriteStringToFileOrDie(matlab_script, matlab_filename); return true; } } // namespace -bool DumpLinearLeastSquaresProblem(const string& directory, - int iteration, +bool DumpLinearLeastSquaresProblem(const string& filename_base, DumpFormatType dump_format_type, const SparseMatrix* A, const double* D, @@ -745,19 +731,16 @@ bool DumpLinearLeastSquaresProblem(const string& directory, int num_eliminate_blocks) { switch (dump_format_type) { case CONSOLE: - return DumpLinearLeastSquaresProblemToConsole(directory, - iteration, - A, D, b, x, + return DumpLinearLeastSquaresProblemToConsole(A, D, b, x, num_eliminate_blocks); case PROTOBUF: return DumpLinearLeastSquaresProblemToProtocolBuffer( - directory, - iteration, + filename_base, A, D, b, x, num_eliminate_blocks); + case TEXTFILE: - return DumpLinearLeastSquaresProblemToTextFile(directory, - iteration, + return DumpLinearLeastSquaresProblemToTextFile(filename_base, A, D, b, x, num_eliminate_blocks); default: diff --git a/internal/ceres/linear_least_squares_problems.h b/internal/ceres/linear_least_squares_problems.h index c76ae91c7..0b550ff9d 100644 --- a/internal/ceres/linear_least_squares_problems.h +++ b/internal/ceres/linear_least_squares_problems.h @@ -73,8 +73,7 @@ LinearLeastSquaresProblem* LinearLeastSquaresProblem3(); // Write the linear least squares problem to disk. The exact format // depends on dump_format_type. -bool DumpLinearLeastSquaresProblem(const string& directory, - int iteration, +bool DumpLinearLeastSquaresProblem(const string& filename_base, DumpFormatType dump_format_type, const SparseMatrix* A, const double* D, diff --git a/internal/ceres/minimizer.h b/internal/ceres/minimizer.h index 18ce64d2d..a2ebe2f1b 100644 --- a/internal/ceres/minimizer.h +++ b/internal/ceres/minimizer.h @@ -73,9 +73,12 @@ class Minimizer { use_nonmonotonic_steps = options.use_nonmonotonic_steps; max_consecutive_nonmonotonic_steps = options.max_consecutive_nonmonotonic_steps; - lsqp_dump_directory = options.lsqp_dump_directory; - lsqp_iterations_to_dump = options.lsqp_iterations_to_dump; - lsqp_dump_format_type = options.lsqp_dump_format_type; + trust_region_problem_dump_directory = + options.trust_region_problem_dump_directory; + trust_region_minimizer_iterations_to_dump = + options.trust_region_minimizer_iterations_to_dump; + trust_region_problem_dump_format_type = + options.trust_region_problem_dump_format_type; max_num_consecutive_invalid_steps = options.max_num_consecutive_invalid_steps; min_trust_region_radius = options.min_trust_region_radius; @@ -110,9 +113,9 @@ class Minimizer { bool jacobi_scaling; bool use_nonmonotonic_steps; int max_consecutive_nonmonotonic_steps; - vector lsqp_iterations_to_dump; - DumpFormatType lsqp_dump_format_type; - string lsqp_dump_directory; + vector trust_region_minimizer_iterations_to_dump; + DumpFormatType trust_region_problem_dump_format_type; + string trust_region_problem_dump_directory; int max_num_consecutive_invalid_steps; double min_trust_region_radius; LineSearchDirectionType line_search_direction_type; diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc index f0ac2f6c1..3d48700c6 100644 --- a/internal/ceres/solver_impl.cc +++ b/internal/ceres/solver_impl.cc @@ -391,9 +391,13 @@ void SolverImpl::TrustRegionSolve(const Solver::Options& original_options, summary->num_threads_given = original_options.num_threads; summary->num_threads_used = options.num_threads; - if (options.lsqp_iterations_to_dump.size() > 0) { - LOG(WARNING) << "Dumping linear least squares problems to disk is" - " currently broken. Ignoring Solver::Options::lsqp_iterations_to_dump"; + if (options.trust_region_minimizer_iterations_to_dump.size() > 0 && + options.trust_region_problem_dump_format_type != CONSOLE && + options.trust_region_problem_dump_directory.empty()) { + summary->error = + "Solver::Options::trust_region_problem_dump_directory is empty."; + LOG(ERROR) << summary->error; + return; } event_logger.AddEvent("Init"); diff --git a/internal/ceres/trust_region_minimizer.cc b/internal/ceres/trust_region_minimizer.cc index 3206decb2..d6ae0ab1a 100644 --- a/internal/ceres/trust_region_minimizer.cc +++ b/internal/ceres/trust_region_minimizer.cc @@ -41,6 +41,7 @@ #include "Eigen/Core" #include "ceres/array_utils.h" #include "ceres/evaluator.h" +#include "ceres/file.h" #include "ceres/internal/eigen.h" #include "ceres/internal/scoped_ptr.h" #include "ceres/linear_least_squares_problems.h" @@ -70,25 +71,8 @@ void TrustRegionMinimizer::EstimateScale(const SparseMatrix& jacobian, void TrustRegionMinimizer::Init(const Minimizer::Options& options) { options_ = options; - sort(options_.lsqp_iterations_to_dump.begin(), - options_.lsqp_iterations_to_dump.end()); -} - -bool TrustRegionMinimizer::MaybeDumpLinearLeastSquaresProblem( - const int iteration, - const SparseMatrix* jacobian, - const double* residuals, - const double* step) const { - // TODO(sameeragarwal): Since the use of trust_region_radius has - // moved inside TrustRegionStrategy, its not clear how we dump the - // regularization vector/matrix anymore. - // - // Also num_eliminate_blocks is not visible to the trust region - // minimizer either. - // - // Both of these indicate that this is the wrong place for this - // code, and going forward this should needs fixing/refactoring. - return true; + sort(options_.trust_region_minimizer_iterations_to_dump.begin(), + options_.trust_region_minimizer_iterations_to_dump.end()); } void TrustRegionMinimizer::Minimize(const Minimizer::Options& options, @@ -212,33 +196,38 @@ void TrustRegionMinimizer::Minimize(const Minimizer::Options& options, break; } - iteration_summary = IterationSummary(); - iteration_summary = summary->iterations.back(); - iteration_summary.iteration = summary->iterations.back().iteration + 1; - iteration_summary.step_is_valid = false; - iteration_summary.step_is_successful = false; - const double strategy_start_time = WallTimeInSeconds(); TrustRegionStrategy::PerSolveOptions per_solve_options; per_solve_options.eta = options_.eta; + if (find(options_.trust_region_minimizer_iterations_to_dump.begin(), + options_.trust_region_minimizer_iterations_to_dump.end(), + iteration_summary.iteration) != + options_.trust_region_minimizer_iterations_to_dump.end()) { + per_solve_options.dump_format_type = + options_.trust_region_problem_dump_format_type; + per_solve_options.dump_filename_base = + JoinPath(options_.trust_region_problem_dump_directory, + StringPrintf("ceres_solver_iteration_%03d", + iteration_summary.iteration)); + } else { + per_solve_options.dump_format_type = TEXTFILE; + per_solve_options.dump_filename_base.clear(); + } + TrustRegionStrategy::Summary strategy_summary = strategy->ComputeStep(per_solve_options, jacobian, residuals.data(), trust_region_step.data()); + iteration_summary = IterationSummary(); + iteration_summary.iteration = summary->iterations.back().iteration + 1; iteration_summary.step_solver_time_in_seconds = WallTimeInSeconds() - strategy_start_time; iteration_summary.linear_solver_iterations = strategy_summary.num_iterations; - - if (!MaybeDumpLinearLeastSquaresProblem(iteration_summary.iteration, - jacobian, - residuals.data(), - trust_region_step.data())) { - LOG(FATAL) << "Tried writing linear least squares problem: " - << options.lsqp_dump_directory << "but failed."; - } + iteration_summary.step_is_valid = false; + iteration_summary.step_is_successful = false; double model_cost_change = 0.0; if (strategy_summary.termination_type != FAILURE) { diff --git a/internal/ceres/trust_region_strategy.h b/internal/ceres/trust_region_strategy.h index f150594bb..f7fd25a81 100644 --- a/internal/ceres/trust_region_strategy.h +++ b/internal/ceres/trust_region_strategy.h @@ -31,6 +31,8 @@ #ifndef CERES_INTERNAL_TRUST_REGION_STRATEGY_H_ #define CERES_INTERNAL_TRUST_REGION_STRATEGY_H_ +#include +#include "ceres/internal/port.h" #include "ceres/types.h" namespace ceres { @@ -82,8 +84,22 @@ class TrustRegionStrategy { // Per solve options. struct PerSolveOptions { + PerSolveOptions() + : eta(0), + dump_filename_base(""), + dump_format_type(TEXTFILE) { + } + // Forcing sequence for inexact solves. double eta; + + // If non-empty and dump_format_type is not CONSOLE, the trust + // regions strategy will write the linear system to file(s) with + // name starting with dump_filename_base. If dump_format_type is + // CONSOLE then dump_filename_base will be ignored and the linear + // system will be written to the standard error. + string dump_filename_base; + DumpFormatType dump_format_type; }; struct Summary {