internal/ceres/trust_region_preprocessor.cc

   1 // Ceres Solver - A fast non-linear least squares minimizer
   2 // Copyright 2015 Google Inc. All rights reserved.
   3 // http://ceres-solver.org/
   4 //
   5 // Redistribution and use in source and binary forms, with or without
   6 // modification, are permitted provided that the following conditions are met:
   7 //
   8 // * Redistributions of source code must retain the above copyright notice,
   9 //   this list of conditions and the following disclaimer.
  10 // * Redistributions in binary form must reproduce the above copyright notice,
  11 //   this list of conditions and the following disclaimer in the documentation
  12 //   and/or other materials provided with the distribution.
  13 // * Neither the name of Google Inc. nor the names of its contributors may be
  14 //   used to endorse or promote products derived from this software without
  15 //   specific prior written permission.
  16 //
  17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  27 // POSSIBILITY OF SUCH DAMAGE.
  28 //
  29 // Author: sameeragarwal@google.com (Sameer Agarwal)
  30
  31 #include "ceres/trust_region_preprocessor.h"
  32
  33 #include <numeric>
  34 #include <string>
  35 #include "ceres/callbacks.h"
  36 #include "ceres/evaluator.h"
  37 #include "ceres/linear_solver.h"
  38 #include "ceres/minimizer.h"
  39 #include "ceres/parameter_block.h"
  40 #include "ceres/preconditioner.h"
  41 #include "ceres/preprocessor.h"
  42 #include "ceres/problem_impl.h"
  43 #include "ceres/program.h"
  44 #include "ceres/reorder_program.h"
  45 #include "ceres/suitesparse.h"
  46 #include "ceres/trust_region_strategy.h"
  47 #include "ceres/wall_time.h"
  48
  49 namespace ceres {
  50 namespace internal {
  51
  52 using std::vector;
  53
  54 namespace {
  55
  56 ParameterBlockOrdering* CreateDefaultLinearSolverOrdering(
  57     const Program& program) {
  58   ParameterBlockOrdering* ordering = new ParameterBlockOrdering;
  59   const vector<ParameterBlock*>& parameter_blocks =
  60       program.parameter_blocks();
  61   for (int i = 0; i < parameter_blocks.size(); ++i) {
  62     ordering->AddElementToGroup(
  63         const_cast<double*>(parameter_blocks[i]->user_state()), 0);
  64   }
  65   return ordering;
  66 }
  67
  68 // Check if all the user supplied values in the parameter blocks are
  69 // sane or not, and if the program is feasible or not.
  70 bool IsProgramValid(const Program& program, std::string* error) {
  71   return (program.ParameterBlocksAreFinite(error) &&
  72           program.IsFeasible(error));
  73 }
  74
  75 void AlternateLinearSolverAndPreconditionerForSchurTypeLinearSolver(
  76     Solver::Options* options) {
  77   if (!IsSchurType(options->linear_solver_type)) {
  78     return;
  79   }
  80
  81   const LinearSolverType linear_solver_type_given = options->linear_solver_type;
  82   const PreconditionerType preconditioner_type_given =
  83       options->preconditioner_type;
  84   options->linear_solver_type = LinearSolver::LinearSolverForZeroEBlocks(
  85       linear_solver_type_given);
  86
  87   std::string message;
  88   if (linear_solver_type_given == ITERATIVE_SCHUR) {
  89     options->preconditioner_type = Preconditioner::PreconditionerForZeroEBlocks(
  90         preconditioner_type_given);
  91
  92     message =
  93         StringPrintf(
  94             "No E blocks. Switching from %s(%s) to %s(%s).",
  95             LinearSolverTypeToString(linear_solver_type_given),
  96             PreconditionerTypeToString(preconditioner_type_given),
  97             LinearSolverTypeToString(options->linear_solver_type),
  98             PreconditionerTypeToString(options->preconditioner_type));
  99   } else {
 100     message =
 101         StringPrintf(
 102             "No E blocks. Switching from %s to %s.",
 103             LinearSolverTypeToString(linear_solver_type_given),
 104             LinearSolverTypeToString(options->linear_solver_type));
 105   }
 106
 107   VLOG_IF(1, options->logging_type != SILENT) << message;
 108 }
 109
 110 // For Schur type and SPARSE_NORMAL_CHOLESKY linear solvers, reorder
 111 // the program to reduce fill-in and increase cache coherency.
 112 bool ReorderProgram(PreprocessedProblem* pp) {
 113   Solver::Options& options = pp->options;
 114   if (IsSchurType(options.linear_solver_type)) {
 115     return ReorderProgramForSchurTypeLinearSolver(
 116         options.linear_solver_type,
 117         options.sparse_linear_algebra_library_type,
 118         pp->problem->parameter_map(),
 119         options.linear_solver_ordering.get(),
 120         pp->reduced_program.get(),
 121         &pp->error);
 122   }
 123
 124
 125   if (options.linear_solver_type == SPARSE_NORMAL_CHOLESKY &&
 126       !options.dynamic_sparsity) {
 127     return ReorderProgramForSparseNormalCholesky(
 128         options.sparse_linear_algebra_library_type,
 129         *options.linear_solver_ordering,
 130         pp->reduced_program.get(),
 131         &pp->error);
 132   }
 133
 134   return true;
 135 }
 136
 137 // Configure and create a linear solver object. In doing so, if a
 138 // sparse direct factorization based linear solver is being used, then
 139 // find a fill reducing ordering and reorder the program as needed
 140 // too.
 141 bool SetupLinearSolver(PreprocessedProblem* pp) {
 142   Solver::Options& options = pp->options;
 143   if (options.linear_solver_ordering.get() == NULL) {
 144     // If the user has not supplied a linear solver ordering, then we
 145     // assume that they are giving all the freedom to us in choosing
 146     // the best possible ordering. This intent can be indicated by
 147     // putting all the parameter blocks in the same elimination group.
 148     options.linear_solver_ordering.reset(
 149         CreateDefaultLinearSolverOrdering(*pp->reduced_program));
 150   } else {
 151     // If the user supplied an ordering, then check if the first
 152     // elimination group is still non-empty after the reduced problem
 153     // has been constructed.
 154     //
 155     // This is important for Schur type linear solvers, where the
 156     // first elimination group is special -- it needs to be an
 157     // independent set.
 158     //
 159     // If the first elimination group is empty, then we cannot use the
 160     // user's requested linear solver (and a preconditioner as the
 161     // case may be) so we must use a different one.
 162     ParameterBlockOrdering* ordering = options.linear_solver_ordering.get();
 163     const int min_group_id = ordering->MinNonZeroGroup();
 164     ordering->Remove(pp->removed_parameter_blocks);
 165     if (IsSchurType(options.linear_solver_type) &&
 166         min_group_id != ordering->MinNonZeroGroup()) {
 167       AlternateLinearSolverAndPreconditionerForSchurTypeLinearSolver(
 168           &options);
 169     }
 170   }
 171
 172   // Reorder the program to reduce fill in and improve cache coherency
 173   // of the Jacobian.
 174   if (!ReorderProgram(pp)) {
 175     return false;
 176   }
 177
 178   // Configure the linear solver.
 179   pp->linear_solver_options = LinearSolver::Options();
 180   pp->linear_solver_options.min_num_iterations =
 181       options.min_linear_solver_iterations;
 182   pp->linear_solver_options.max_num_iterations =
 183       options.max_linear_solver_iterations;
 184   pp->linear_solver_options.type = options.linear_solver_type;
 185   pp->linear_solver_options.preconditioner_type = options.preconditioner_type;
 186   pp->linear_solver_options.visibility_clustering_type =
 187       options.visibility_clustering_type;
 188   pp->linear_solver_options.sparse_linear_algebra_library_type =
 189       options.sparse_linear_algebra_library_type;
 190   pp->linear_solver_options.dense_linear_algebra_library_type =
 191       options.dense_linear_algebra_library_type;
 192   pp->linear_solver_options.use_explicit_schur_complement =
 193       options.use_explicit_schur_complement;
 194   pp->linear_solver_options.dynamic_sparsity = options.dynamic_sparsity;
 195   pp->linear_solver_options.num_threads = options.num_linear_solver_threads;
 196   pp->linear_solver_options.use_postordering = options.use_postordering;
 197
 198   if (IsSchurType(pp->linear_solver_options.type)) {
 199     OrderingToGroupSizes(options.linear_solver_ordering.get(),
 200                          &pp->linear_solver_options.elimination_groups);
 201
 202     // Schur type solvers expect at least two elimination groups. If
 203     // there is only one elimination group, then it is guaranteed that
 204     // this group only contains e_blocks. Thus we add a dummy
 205     // elimination group with zero blocks in it.
 206     if (pp->linear_solver_options.elimination_groups.size() == 1) {
 207       pp->linear_solver_options.elimination_groups.push_back(0);
 208     }
 209
 210     if (options.linear_solver_type == SPARSE_SCHUR) {
 211       // When using SPARSE_SCHUR, we ignore the user's postordering
 212       // preferences in certain cases.
 213       //
 214       // 1. SUITE_SPARSE is the sparse linear algebra library requested
 215       //    but cholmod_camd is not available.
 216       // 2. CX_SPARSE is the sparse linear algebra library requested.
 217       //
 218       // This ensures that the linear solver does not assume that a
 219       // fill-reducing pre-ordering has been done.
 220       //
 221       // TODO(sameeragarwal): Implement the reordering of parameter
 222       // blocks for CX_SPARSE.
 223       if ((options.sparse_linear_algebra_library_type == SUITE_SPARSE &&
 224            !SuiteSparse::
 225            IsConstrainedApproximateMinimumDegreeOrderingAvailable()) ||
 226           (options.sparse_linear_algebra_library_type == CX_SPARSE)) {
 227         pp->linear_solver_options.use_postordering = true;
 228       }
 229     }
 230   }
 231
 232   pp->linear_solver.reset(LinearSolver::Create(pp->linear_solver_options));
 233   return (pp->linear_solver.get() != NULL);
 234 }
 235
 236 // Configure and create the evaluator.
 237 bool SetupEvaluator(PreprocessedProblem* pp) {
 238   const Solver::Options& options = pp->options;
 239   pp->evaluator_options = Evaluator::Options();
 240   pp->evaluator_options.linear_solver_type = options.linear_solver_type;
 241   pp->evaluator_options.num_eliminate_blocks = 0;
 242   if (IsSchurType(options.linear_solver_type)) {
 243     pp->evaluator_options.num_eliminate_blocks =
 244         options
 245         .linear_solver_ordering
 246         ->group_to_elements().begin()
 247         ->second.size();
 248   }
 249
 250   pp->evaluator_options.num_threads = options.num_threads;
 251   pp->evaluator_options.dynamic_sparsity = options.dynamic_sparsity;
 252   pp->evaluator.reset(Evaluator::Create(pp->evaluator_options,
 253                                         pp->reduced_program.get(),
 254                                         &pp->error));
 255
 256   return (pp->evaluator.get() != NULL);
 257 }
 258
 259 // If the user requested inner iterations, then find an inner
 260 // iteration ordering as needed and configure and create a
 261 // CoordinateDescentMinimizer object to perform the inner iterations.
 262 bool SetupInnerIterationMinimizer(PreprocessedProblem* pp) {
 263   Solver::Options& options = pp->options;
 264   if (!options.use_inner_iterations) {
 265     return true;
 266   }
 267
 268   // With just one parameter block, the outer iteration of the trust
 269   // region method and inner iterations are doing exactly the same
 270   // thing, and thus inner iterations are not needed.
 271   if (pp->reduced_program->NumParameterBlocks() == 1) {
 272     LOG(WARNING) << "Reduced problem only contains one parameter block."
 273                  << "Disabling inner iterations.";
 274     return true;
 275   }
 276
 277   if (options.inner_iteration_ordering.get() != NULL) {
 278     // If the user supplied an ordering, then remove the set of
 279     // inactive parameter blocks from it
 280     options.inner_iteration_ordering->Remove(pp->removed_parameter_blocks);
 281     if (options.inner_iteration_ordering->NumElements() == 0) {
 282       LOG(WARNING) << "No remaining elements in the inner iteration ordering.";
 283       return true;
 284     }
 285
 286     // Validate the reduced ordering.
 287     if (!CoordinateDescentMinimizer::IsOrderingValid(
 288             *pp->reduced_program,
 289             *options.inner_iteration_ordering,
 290             &pp->error)) {
 291       return false;
 292     }
 293   } else {
 294     // The user did not supply an ordering, so create one.
 295     options.inner_iteration_ordering.reset(
 296         CoordinateDescentMinimizer::CreateOrdering(*pp->reduced_program));
 297   }
 298
 299   pp->inner_iteration_minimizer.reset(new CoordinateDescentMinimizer);
 300   return pp->inner_iteration_minimizer->Init(*pp->reduced_program,
 301                                              pp->problem->parameter_map(),
 302                                              *options.inner_iteration_ordering,
 303                                              &pp->error);
 304 }
 305
 306 // Configure and create a TrustRegionMinimizer object.
 307 void SetupMinimizerOptions(PreprocessedProblem* pp) {
 308   const Solver::Options& options = pp->options;
 309
 310   SetupCommonMinimizerOptions(pp);
 311   pp->minimizer_options.is_constrained =
 312       pp->reduced_program->IsBoundsConstrained();
 313   pp->minimizer_options.jacobian.reset(pp->evaluator->CreateJacobian());
 314   pp->minimizer_options.inner_iteration_minimizer =
 315       pp->inner_iteration_minimizer;
 316
 317   TrustRegionStrategy::Options strategy_options;
 318   strategy_options.linear_solver = pp->linear_solver.get();
 319   strategy_options.initial_radius =
 320       options.initial_trust_region_radius;
 321   strategy_options.max_radius = options.max_trust_region_radius;
 322   strategy_options.min_lm_diagonal = options.min_lm_diagonal;
 323   strategy_options.max_lm_diagonal = options.max_lm_diagonal;
 324   strategy_options.trust_region_strategy_type =
 325       options.trust_region_strategy_type;
 326   strategy_options.dogleg_type = options.dogleg_type;
 327   pp->minimizer_options.trust_region_strategy.reset(
 328       CHECK_NOTNULL(TrustRegionStrategy::Create(strategy_options)));
 329 }
 330
 331 }  // namespace
 332
 333 TrustRegionPreprocessor::~TrustRegionPreprocessor() {
 334 }
 335
 336 bool TrustRegionPreprocessor::Preprocess(const Solver::Options& options,
 337                                          ProblemImpl* problem,
 338                                          PreprocessedProblem* pp) {
 339   CHECK_NOTNULL(pp);
 340   pp->options = options;
 341   ChangeNumThreadsIfNeeded(&pp->options);
 342
 343   pp->problem = problem;
 344   Program* program = problem->mutable_program();
 345   if (!IsProgramValid(*program, &pp->error)) {
 346     return false;
 347   }
 348
 349   pp->reduced_program.reset(
 350       program->CreateReducedProgram(&pp->removed_parameter_blocks,
 351                                     &pp->fixed_cost,
 352                                     &pp->error));
 353
 354   if (pp->reduced_program.get() == NULL) {
 355     return false;
 356   }
 357
 358   if (pp->reduced_program->NumParameterBlocks() == 0) {
 359     // The reduced problem has no parameter or residual blocks. There
 360     // is nothing more to do.
 361     return true;
 362   }
 363
 364   if (!SetupLinearSolver(pp) ||
 365       !SetupEvaluator(pp) ||
 366       !SetupInnerIterationMinimizer(pp)) {
 367     return false;
 368   }
 369
 370   SetupMinimizerOptions(pp);
 371   return true;
 372 }
 373
 374 }  // namespace internal
 375 }  // namespace ceres