diff --git a/include/ceres/iteration_callback.h b/include/ceres/iteration_callback.h
index 6bab004..db5d0ef 100644
--- a/include/ceres/iteration_callback.h
+++ b/include/ceres/iteration_callback.h
@@ -69,7 +69,7 @@
   // Step was numerically valid, i.e., all values are finite and the
   // step reduces the value of the linearized model.
   //
-  // Note: step_is_valid is false when iteration = 0.
+  // Note: step_is_valid is always true when iteration = 0.
   bool step_is_valid;
 
   // Step did not reduce the value of the objective function
@@ -77,7 +77,7 @@
   // acceptance criterion used by the non-monotonic trust region
   // algorithm.
   //
-  // Note: step_is_nonmonotonic is false when iteration = 0;
+  // Note: step_is_nonmonotonic is always false when iteration = 0;
   bool step_is_nonmonotonic;
 
   // Whether or not the minimizer accepted this step or not. If the
@@ -89,7 +89,7 @@
   // relative decrease is not sufficient, the algorithm may accept the
   // step and the step is declared successful.
   //
-  // Note: step_is_successful is false when iteration = 0.
+  // Note: step_is_successful is always true when iteration = 0.
   bool step_is_successful;
 
   // Value of the objective function.
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 3a51309..a9fee1a 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -108,6 +108,7 @@
     triplet_sparse_matrix.cc
     trust_region_preprocessor.cc
     trust_region_minimizer.cc
+    trust_region_step_evaluator.cc
     trust_region_strategy.cc
     types.cc
     visibility.cc
diff --git a/internal/ceres/trust_region_minimizer.cc b/internal/ceres/trust_region_minimizer.cc
index 627430c..7a4a775 100644
--- a/internal/ceres/trust_region_minimizer.cc
+++ b/internal/ceres/trust_region_minimizer.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2016 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -43,681 +43,732 @@
 #include "ceres/coordinate_descent_minimizer.h"
 #include "ceres/evaluator.h"
 #include "ceres/file.h"
-#include "ceres/internal/eigen.h"
-#include "ceres/internal/scoped_ptr.h"
 #include "ceres/line_search.h"
-#include "ceres/linear_least_squares_problems.h"
-#include "ceres/sparse_matrix.h"
 #include "ceres/stringprintf.h"
-#include "ceres/trust_region_strategy.h"
 #include "ceres/types.h"
 #include "ceres/wall_time.h"
 #include "glog/logging.h"
 
+// Helper macro to simplify some of the control flow.
+#define RETURN_IF_ERROR_AND_LOG(expr)                            \
+  do {                                                           \
+    if (!(expr)) {                                               \
+      LOG(ERROR) << "Terminating: " << solver_summary_->message; \
+      return;                                                    \
+    }                                                            \
+  } while (0)
+
 namespace ceres {
 namespace internal {
-namespace {
 
-LineSearch::Summary DoLineSearch(const Minimizer::Options& options,
-                                 const Vector& x,
-                                 const Vector& gradient,
-                                 const double cost,
-                                 const Vector& delta,
-                                 Evaluator* evaluator) {
-  LineSearchFunction line_search_function(evaluator);
+TrustRegionMinimizer::~TrustRegionMinimizer() {}
+
+void TrustRegionMinimizer::Minimize(const Minimizer::Options& options,
+                                    double* parameters,
+                                    Solver::Summary* solver_summary) {
+  start_time_ = WallTimeInSeconds();
+  iteration_start_time_ = start_time_;
+  Init(options, parameters, solver_summary);
+  RETURN_IF_ERROR_AND_LOG(IterationZero());
+
+  // Create the TrustRegionStepEvaluator. The construction needs to be
+  // delayed to this point because we need the cost for the starting
+  // point to initialize the step evaluator.
+  step_evaluator_.reset(new TrustRegionStepEvaluator(
+      x_cost_,
+      options_.use_nonmonotonic_steps
+          ? options_.max_consecutive_nonmonotonic_steps
+          : 0));
+
+  while (FinalizeIterationAndCheckIfMinimizerCanContinue()) {
+    iteration_start_time_ = WallTimeInSeconds();
+    iteration_summary_ = IterationSummary();
+    iteration_summary_.iteration =
+        solver_summary->iterations.back().iteration + 1;
+
+    RETURN_IF_ERROR_AND_LOG(ComputeTrustRegionStep());
+    if (!iteration_summary_.step_is_valid) {
+      RETURN_IF_ERROR_AND_LOG(HandleInvalidStep());
+      continue;
+    }
+
+    if (options_.is_constrained) {
+      // Use a projected line search to enforce the bounds constraints
+      // and improve the quality of the step.
+      DoLineSearch(x_, gradient_, x_cost_, &delta_);
+    }
+
+    ComputeCandidatePointAndEvaluateCost();
+    DoInnerIterationsIfNeeded();
+
+    if (ParameterToleranceReached()) {
+      return;
+    }
+
+    if (FunctionToleranceReached()) {
+      return;
+    }
+
+    if (IsStepSuccessful()) {
+      RETURN_IF_ERROR_AND_LOG(HandleSuccessfulStep());
+      continue;
+    }
+
+    HandleUnsuccessfulStep();
+  }
+}
+
+// Initialize the minimizer, allocate working space and set some of
+// the fields in the solver_summary.
+void TrustRegionMinimizer::Init(const Minimizer::Options& options,
+                                double* parameters,
+                                Solver::Summary* solver_summary) {
+  options_ = options;
+  sort(options_.trust_region_minimizer_iterations_to_dump.begin(),
+       options_.trust_region_minimizer_iterations_to_dump.end());
+
+  parameters_ = parameters;
+
+  solver_summary_ = solver_summary;
+  solver_summary_->termination_type = NO_CONVERGENCE;
+  solver_summary_->num_successful_steps = 0;
+  solver_summary_->num_unsuccessful_steps = 0;
+  solver_summary_->is_constrained = options.is_constrained;
+
+  evaluator_ = CHECK_NOTNULL(options_.evaluator.get());
+  jacobian_ = CHECK_NOTNULL(options_.jacobian.get());
+  strategy_ = CHECK_NOTNULL(options_.trust_region_strategy.get());
+
+  is_not_silent_ = !options.is_silent;
+  inner_iterations_are_enabled_ =
+      options.inner_iteration_minimizer.get() != NULL;
+  inner_iterations_were_useful_ = false;
+
+  num_parameters_ = evaluator_->NumParameters();
+  num_effective_parameters_ = evaluator_->NumEffectiveParameters();
+  num_residuals_ = evaluator_->NumResiduals();
+  num_consecutive_invalid_steps_ = 0;
+
+  x_ = ConstVectorRef(parameters_, num_parameters_);
+  x_norm_ = x_.norm();
+  residuals_.resize(num_residuals_);
+  trust_region_step_.resize(num_effective_parameters_);
+  delta_.resize(num_effective_parameters_);
+  candidate_x_.resize(num_parameters_);
+  gradient_.resize(num_effective_parameters_);
+  model_residuals_.resize(num_residuals_);
+  negative_gradient_.resize(num_effective_parameters_);
+  projected_gradient_step_.resize(num_parameters_);
+
+  // By default scaling is one, if the user requests Jacobi scaling of
+  // the Jacobian, we will compute and overwrite this vector.
+  jacobian_scaling_ = Vector::Ones(num_effective_parameters_);
+
+  x_norm_ = -1;  // Invalid value
+  x_cost_ = std::numeric_limits<double>::max();
+  minimum_cost_ = x_cost_;
+  model_cost_change_ = 0.0;
+}
+
+// 1. Project the initial solution onto the feasible set if needed.
+// 2. Compute the initial cost, jacobian & gradient.
+//
+// Return true if all computations can be performed successfully.
+bool TrustRegionMinimizer::IterationZero() {
+  iteration_summary_ = IterationSummary();
+  iteration_summary_.iteration = 0;
+  iteration_summary_.step_is_valid = false;
+  iteration_summary_.step_is_successful = false;
+  iteration_summary_.cost_change = 0.0;
+  iteration_summary_.gradient_max_norm = 0.0;
+  iteration_summary_.gradient_norm = 0.0;
+  iteration_summary_.step_norm = 0.0;
+  iteration_summary_.relative_decrease = 0.0;
+  iteration_summary_.eta = options_.eta;
+  iteration_summary_.linear_solver_iterations = 0;
+  iteration_summary_.step_solver_time_in_seconds = 0;
+
+  if (options_.is_constrained) {
+    delta_.setZero();
+    if (!evaluator_->Plus(x_.data(), delta_.data(), candidate_x_.data())) {
+      solver_summary_->message =
+          "Unable to project initial point onto the feasible set.";
+      solver_summary_->termination_type = FAILURE;
+      return false;
+    }
+
+    x_ = candidate_x_;
+    x_norm_ = x_.norm();
+  }
+
+  if (!EvaluateGradientAndJacobian()) {
+    return false;
+  }
+
+  solver_summary_->initial_cost = x_cost_ + solver_summary_->fixed_cost;
+  iteration_summary_.step_is_valid = true;
+  iteration_summary_.step_is_successful = true;
+  return true;
+}
+
+// For the current x_, compute
+//
+//  1. Cost
+//  2. Jacobian
+//  3. Gradient
+//  4. Scale the Jacobian if needed (and compute the scaling if we are
+//     in iteration zero).
+//  5. Compute the 2 and max norm of the gradient.
+//
+// Returns true if all computations could be performed
+// successfully. Any failures are considered fatal and the
+// Solver::Summary is updated to indicate this.
+bool TrustRegionMinimizer::EvaluateGradientAndJacobian() {
+  if (!evaluator_->Evaluate(x_.data(),
+                            &x_cost_,
+                            residuals_.data(),
+                            gradient_.data(),
+                            jacobian_)) {
+    solver_summary_->message = "Residual and Jacobian evaluation failed.";
+    solver_summary_->termination_type = FAILURE;
+    return false;
+  }
+
+  iteration_summary_.cost = x_cost_ + solver_summary_->fixed_cost;
+
+  if (options_.jacobi_scaling) {
+    if (iteration_summary_.iteration == 0) {
+      // Compute a scaling vector that is used to improve the
+      // conditioning of the Jacobian.
+      //
+      // jacobian_scaling_ = diag(J'J)^{-1}
+      jacobian_->SquaredColumnNorm(jacobian_scaling_.data());
+      for (int i = 0; i < jacobian_->num_cols(); ++i) {
+        // Add one to the denominator to prevent division by zero.
+        jacobian_scaling_[i] = 1.0 / (1.0 + sqrt(jacobian_scaling_[i]));
+      }
+    }
+
+    // jacobian = jacobian * diag(J'J) ^{-1}
+    jacobian_->ScaleColumns(jacobian_scaling_.data());
+  }
+
+  // The gradient exists in the local tangent space. To account for
+  // the bounds constraints correctly, instead of just computing the
+  // norm of the gradient vector, we compute
+  //
+  // |Plus(x, -gradient) - x|
+  //
+  // Where the Plus operator lifts the negative gradient to the
+  // ambient space, adds it to x and projects it on the hypercube
+  // defined by the bounds.
+  negative_gradient_ = -gradient_;
+  if (!evaluator_->Plus(x_.data(),
+                        negative_gradient_.data(),
+                        projected_gradient_step_.data())) {
+    solver_summary_->message =
+        "projected_gradient_step = Plus(x, -gradient) failed.";
+    solver_summary_->termination_type = FAILURE;
+    return false;
+  }
+
+  iteration_summary_.gradient_max_norm =
+      (x_ - projected_gradient_step_).lpNorm<Eigen::Infinity>();
+  iteration_summary_.gradient_norm = (x_ - projected_gradient_step_).norm();
+  return true;
+}
+
+// 1. Add the final timing information to the iteration summary.
+// 2. Run the callbacks
+// 3. Check for termination based on
+//    a. Run time
+//    b. Iteration count
+//    c. Max norm of the gradient
+//    d. Size of the trust region radius.
+//
+// Returns true if user did not terminate the solver and none of these
+// termination criterion are met.
+bool TrustRegionMinimizer::FinalizeIterationAndCheckIfMinimizerCanContinue() {
+  if (iteration_summary_.step_is_successful) {
+    ++solver_summary_->num_successful_steps;
+    if (x_cost_ < minimum_cost_) {
+      minimum_cost_ = x_cost_;
+      VectorRef(parameters_, num_parameters_) = x_;
+      iteration_summary_.step_is_nonmonotonic = false;
+    } else {
+      iteration_summary_.step_is_nonmonotonic = true;
+    }
+  } else {
+    ++solver_summary_->num_unsuccessful_steps;
+  }
+
+  iteration_summary_.trust_region_radius = strategy_->Radius();
+  iteration_summary_.iteration_time_in_seconds =
+      WallTimeInSeconds() - iteration_start_time_;
+  iteration_summary_.cumulative_time_in_seconds =
+      WallTimeInSeconds() - start_time_ +
+      solver_summary_->preprocessor_time_in_seconds;
+
+  solver_summary_->iterations.push_back(iteration_summary_);
+
+  if (!RunCallbacks(options_, iteration_summary_, solver_summary_)) {
+    return false;
+  }
+
+  if (MaxSolverTimeReached()) {
+    return false;
+  }
+
+  if (MaxSolverIterationsReached()) {
+    return false;
+  }
+
+  if (GradientToleranceReached()) {
+    return false;
+  }
+
+  if (MinTrustRegionRadiusReached()) {
+    return false;
+  }
+
+  return true;
+}
+
+// Compute the trust region step using the TrustRegionStrategy chosen
+// by the user.
+//
+// If the strategy returns with LINEAR_SOLVER_FATAL_ERROR, which
+// indicates an unrecoverable error, return false. This is the only
+// condition that returns false.
+//
+// If the strategy returns with LINEAR_SOLVER_FAILURE, which indicates
+// a numerical failure that could be recovered from by retrying
+// (e.g. by increasing the strength of the regularization), we set
+// iteration_summary_.step_is_valid to false and return true.
+//
+// In all other cases, we compute the decrease in the trust region
+// model problem. In exact arithmetic, this should always be
+// positive, but due to numerical problems in the TrustRegionStrategy
+// or round off error when computing the decrease it may be
+// negative. In which case again, we set
+// iteration_summary_.step_is_valid to false.
+bool TrustRegionMinimizer::ComputeTrustRegionStep() {
+  const double strategy_start_time = WallTimeInSeconds();
+  iteration_summary_.step_is_valid = false;
+  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));
+  }
+
+  TrustRegionStrategy::Summary strategy_summary =
+      strategy_->ComputeStep(per_solve_options,
+                             jacobian_,
+                             residuals_.data(),
+                             trust_region_step_.data());
+
+  if (strategy_summary.termination_type == LINEAR_SOLVER_FATAL_ERROR) {
+    solver_summary_->message =
+        "Linear solver failed due to unrecoverable "
+        "non-numeric causes. Please see the error log for clues. ";
+    solver_summary_->termination_type = FAILURE;
+    return false;
+  }
+
+  iteration_summary_.step_solver_time_in_seconds =
+      WallTimeInSeconds() - strategy_start_time;
+  iteration_summary_.linear_solver_iterations = strategy_summary.num_iterations;
+
+  if (strategy_summary.termination_type == LINEAR_SOLVER_FAILURE) {
+    return true;
+  }
+
+  // new_model_cost
+  //  = 1/2 [f + J * step]^2
+  //  = 1/2 [ f'f + 2f'J * step + step' * J' * J * step ]
+  // model_cost_change
+  //  = cost - new_model_cost
+  //  = f'f/2  - 1/2 [ f'f + 2f'J * step + step' * J' * J * step]
+  //  = -f'J * step - step' * J' * J * step / 2
+  //  = -(J * step)'(f + J * step / 2)
+  model_residuals_.setZero();
+  jacobian_->RightMultiply(trust_region_step_.data(), model_residuals_.data());
+  model_cost_change_ =
+      -model_residuals_.dot(residuals_ + model_residuals_ / 2.0);
+
+  // TODO(sameeragarwal)
+  //
+  //  1. What happens if model_cost_change_ = 0
+  //  2. What happens if -epsilon <= model_cost_change_ < 0 for some
+  //     small epsilon due to round off error.
+  iteration_summary_.step_is_valid = (model_cost_change_ > 0.0);
+  if (iteration_summary_.step_is_valid) {
+    // Undo the Jacobian column scaling.
+    delta_ = (trust_region_step_.array() * jacobian_scaling_.array()).matrix();
+    num_consecutive_invalid_steps_ = 0;
+  }
+
+  VLOG_IF(1, is_not_silent_ && !iteration_summary_.step_is_valid)
+      << "Invalid step: current_cost: " << x_cost_
+      << " absolute model cost change: " << model_cost_change_
+      << " relative model cost change: " << (model_cost_change_ / x_cost_);
+  return true;
+}
+
+// Invalid steps can happen due to a number of reasons, and we allow a
+// limited number of consecutive failures, and return false if this
+// limit is exceeded.
+bool TrustRegionMinimizer::HandleInvalidStep() {
+  // TODO(sameeragarwal): Should we be returning FAILURE or
+  // NO_CONVERGENCE? The solution value is still usable in many cases,
+  // it is not clear if we should declare the solver a failure
+  // entirely. For example the case where model_cost_change ~ 0.0, but
+  // just slightly negative.
+  if (++num_consecutive_invalid_steps_ >=
+      options_.max_num_consecutive_invalid_steps) {
+    solver_summary_->message = StringPrintf(
+        "Number of consecutive invalid steps more "
+        "than Solver::Options::max_num_consecutive_invalid_steps: %d",
+        options_.max_num_consecutive_invalid_steps);
+    solver_summary_->termination_type = FAILURE;
+    return false;
+  }
+
+  strategy_->StepIsInvalid();
+
+  // We are going to try and reduce the trust region radius and
+  // solve again. To do this, we are going to treat this iteration
+  // as an unsuccessful iteration. Since the various callbacks are
+  // still executed, we are going to fill the iteration summary
+  // with data that assumes a step of length zero and no progress.
+  iteration_summary_.cost = x_cost_ + solver_summary_->fixed_cost;
+  iteration_summary_.cost_change = 0.0;
+  iteration_summary_.gradient_max_norm =
+      solver_summary_->iterations.back().gradient_max_norm;
+  iteration_summary_.gradient_norm =
+      solver_summary_->iterations.back().gradient_norm;
+  iteration_summary_.step_norm = 0.0;
+  iteration_summary_.relative_decrease = 0.0;
+  iteration_summary_.eta = options_.eta;
+  return true;
+}
+
+// Use the supplied coordinate descent minimizer to perform inner
+// iterations and compute the improvement due to it. Returns the cost
+// after performing the inner iterations.
+//
+// The optimization is performed with candidate_x_ as the starting
+// point, and if the optimization is successful, candidate_x_ will be
+// updated with the optimized parameters.
+void TrustRegionMinimizer::DoInnerIterationsIfNeeded() {
+  if (!inner_iterations_are_enabled_ ||
+      candidate_cost_ >= std::numeric_limits<double>::max()) {
+    return;
+  }
+
+  double inner_iteration_start_time = WallTimeInSeconds();
+  ++solver_summary_->num_inner_iteration_steps;
+  inner_iteration_x_ = candidate_x_;
+  Solver::Summary inner_iteration_summary;
+  options_.inner_iteration_minimizer->Minimize(
+      options_, inner_iteration_x_.data(), &inner_iteration_summary);
+  double inner_iteration_cost;
+  if (!evaluator_->Evaluate(
+          inner_iteration_x_.data(), &inner_iteration_cost, NULL, NULL, NULL)) {
+    VLOG_IF(2, is_not_silent_) << "Inner iteration failed.";
+    return;
+  }
+
+  VLOG_IF(2, is_not_silent_)
+      << "Inner iteration succeeded; Current cost: " << x_cost_
+      << " Trust region step cost: " << candidate_cost_
+      << " Inner iteration cost: " << inner_iteration_cost;
+
+  candidate_x_ = inner_iteration_x_;
+
+  // Normally, the quality of a trust region step is measured by
+  // the ratio
+  //
+  //              cost_change
+  //    r =    -----------------
+  //           model_cost_change
+  //
+  // All the change in the nonlinear objective is due to the trust
+  // region step so this ratio is a good measure of the quality of
+  // the trust region radius. However, when inner iterations are
+  // being used, cost_change includes the contribution of the
+  // inner iterations and its not fair to credit it all to the
+  // trust region algorithm. So we change the ratio to be
+  //
+  //                              cost_change
+  //    r =    ------------------------------------------------
+  //           (model_cost_change + inner_iteration_cost_change)
+  //
+  // Practically we do this by increasing model_cost_change by
+  // inner_iteration_cost_change.
+
+  const double inner_iteration_cost_change =
+      candidate_cost_ - inner_iteration_cost;
+  model_cost_change_ += inner_iteration_cost_change;
+  inner_iterations_were_useful_ = inner_iteration_cost < x_cost_;
+  const double inner_iteration_relative_progress =
+      1.0 - inner_iteration_cost / candidate_cost_;
+
+  // Disable inner iterations once the relative improvement
+  // drops below tolerance.
+  inner_iterations_are_enabled_ =
+      (inner_iteration_relative_progress > options_.inner_iteration_tolerance);
+  VLOG_IF(2, is_not_silent_ && !inner_iterations_are_enabled_)
+      << "Disabling inner iterations. Progress : "
+      << inner_iteration_relative_progress;
+  candidate_cost_ = inner_iteration_cost;
+
+  solver_summary_->inner_iteration_time_in_seconds +=
+      WallTimeInSeconds() - inner_iteration_start_time;
+}
+
+// Perform a projected line search to improve the objective function
+// value along delta.
+//
+// TODO(sameeragarwal): The current implementation does not do
+// anything illegal but is incorrect and not terribly effective.
+//
+// https://github.com/ceres-solver/ceres-solver/issues/187
+void TrustRegionMinimizer::DoLineSearch(const Vector& x,
+                                        const Vector& gradient,
+                                        const double cost,
+                                        Vector* delta) {
+  LineSearchFunction line_search_function(evaluator_);
 
   LineSearch::Options line_search_options;
   line_search_options.is_silent = true;
   line_search_options.interpolation_type =
-      options.line_search_interpolation_type;
-  line_search_options.min_step_size = options.min_line_search_step_size;
+      options_.line_search_interpolation_type;
+  line_search_options.min_step_size = options_.min_line_search_step_size;
   line_search_options.sufficient_decrease =
-      options.line_search_sufficient_function_decrease;
+      options_.line_search_sufficient_function_decrease;
   line_search_options.max_step_contraction =
-      options.max_line_search_step_contraction;
+      options_.max_line_search_step_contraction;
   line_search_options.min_step_contraction =
-      options.min_line_search_step_contraction;
+      options_.min_line_search_step_contraction;
   line_search_options.max_num_iterations =
-      options.max_num_line_search_step_size_iterations;
+      options_.max_num_line_search_step_size_iterations;
   line_search_options.sufficient_curvature_decrease =
-      options.line_search_sufficient_curvature_decrease;
+      options_.line_search_sufficient_curvature_decrease;
   line_search_options.max_step_expansion =
-      options.max_line_search_step_expansion;
+      options_.max_line_search_step_expansion;
   line_search_options.function = &line_search_function;
 
   std::string message;
-  scoped_ptr<LineSearch> line_search(
-      CHECK_NOTNULL(LineSearch::Create(ceres::ARMIJO,
-                                       line_search_options,
-                                       &message)));
-  LineSearch::Summary summary;
-  line_search_function.Init(x, delta);
-  line_search->Search(1.0, cost, gradient.dot(delta), &summary);
-  return summary;
-}
+  scoped_ptr<LineSearch> line_search(CHECK_NOTNULL(
+      LineSearch::Create(ceres::ARMIJO, line_search_options, &message)));
+  LineSearch::Summary line_search_summary;
+  line_search_function.Init(x, *delta);
+  line_search->Search(1.0, cost, gradient.dot(*delta), &line_search_summary);
 
-}  // namespace
+  solver_summary_->num_line_search_steps += line_search_summary.num_iterations;
+  solver_summary_->line_search_cost_evaluation_time_in_seconds +=
+      line_search_summary.cost_evaluation_time_in_seconds;
+  solver_summary_->line_search_gradient_evaluation_time_in_seconds +=
+      line_search_summary.gradient_evaluation_time_in_seconds;
+  solver_summary_->line_search_polynomial_minimization_time_in_seconds +=
+      line_search_summary.polynomial_minimization_time_in_seconds;
+  solver_summary_->line_search_total_time_in_seconds +=
+      line_search_summary.total_time_in_seconds;
 
-// Compute a scaling vector that is used to improve the conditioning
-// of the Jacobian.
-void TrustRegionMinimizer::EstimateScale(const SparseMatrix& jacobian,
-                                         double* scale) const {
-  jacobian.SquaredColumnNorm(scale);
-  for (int i = 0; i < jacobian.num_cols(); ++i) {
-    scale[i] = 1.0 / (1.0 + sqrt(scale[i]));
+  if (line_search_summary.success) {
+    *delta *= line_search_summary.optimal_step_size;
   }
 }
 
-void TrustRegionMinimizer::Init(const Minimizer::Options& options) {
-  options_ = options;
-  sort(options_.trust_region_minimizer_iterations_to_dump.begin(),
-       options_.trust_region_minimizer_iterations_to_dump.end());
+bool TrustRegionMinimizer::MaxSolverTimeReached() {
+  const double total_solver_time =
+      WallTimeInSeconds() - start_time_ +
+      solver_summary_->preprocessor_time_in_seconds;
+  if (total_solver_time < options_.max_solver_time_in_seconds) {
+    return false;
+  }
+
+  solver_summary_->message = StringPrintf("Maximum solver time reached. "
+                                          "Total solver time: %e >= %e.",
+                                          total_solver_time,
+                                          options_.max_solver_time_in_seconds);
+  solver_summary_->termination_type = NO_CONVERGENCE;
+  VLOG_IF(1, is_not_silent_) << "Terminating: " << solver_summary_->message;
+  return true;
 }
 
-void TrustRegionMinimizer::Minimize(const Minimizer::Options& options,
-                                    double* parameters,
-                                    Solver::Summary* summary) {
-  double start_time = WallTimeInSeconds();
-  double iteration_start_time =  start_time;
-  Init(options);
+bool TrustRegionMinimizer::MaxSolverIterationsReached() {
+  if (iteration_summary_.iteration < options_.max_num_iterations) {
+    return false;
+  }
 
-  Evaluator* evaluator = CHECK_NOTNULL(options_.evaluator.get());
-  SparseMatrix* jacobian = CHECK_NOTNULL(options_.jacobian.get());
-  TrustRegionStrategy* strategy =
-      CHECK_NOTNULL(options_.trust_region_strategy.get());
+  solver_summary_->message =
+      StringPrintf("Maximum number of iterations reached. "
+                   "Number of iterations: %d.",
+                   iteration_summary_.iteration);
 
-  const bool is_not_silent = !options.is_silent;
+  solver_summary_->termination_type = NO_CONVERGENCE;
+  VLOG_IF(1, is_not_silent_) << "Terminating: " << solver_summary_->message;
+  return true;
+}
 
-  // If the problem is bounds constrained, then enable the use of a
-  // line search after the trust region step has been computed. This
-  // line search will automatically use a projected test point onto
-  // the feasible set, there by guaranteeing the feasibility of the
-  // final output.
+bool TrustRegionMinimizer::GradientToleranceReached() {
+  if (!iteration_summary_.step_is_successful ||
+      iteration_summary_.gradient_max_norm > options_.gradient_tolerance) {
+    return false;
+  }
+
+  solver_summary_->message = StringPrintf(
+      "Gradient tolerance reached. "
+      "Gradient max norm: %e <= %e",
+      iteration_summary_.gradient_max_norm,
+      options_.gradient_tolerance);
+  solver_summary_->termination_type = CONVERGENCE;
+  VLOG_IF(1, is_not_silent_) << "Terminating: " << solver_summary_->message;
+  return true;
+}
+
+bool TrustRegionMinimizer::MinTrustRegionRadiusReached() {
+  if (iteration_summary_.trust_region_radius >
+      options_.min_trust_region_radius) {
+    return false;
+  }
+
+  solver_summary_->message =
+      StringPrintf("Minimum trust region radius reached. "
+                   "Trust region radius: %e <= %e",
+                   iteration_summary_.trust_region_radius,
+                   options_.min_trust_region_radius);
+  solver_summary_->termination_type = CONVERGENCE;
+  VLOG_IF(1, is_not_silent_) << "Terminating: " << solver_summary_->message;
+  return true;
+}
+
+// Solver::Options::parameter_tolerance based convergence check.
+bool TrustRegionMinimizer::ParameterToleranceReached() {
+  iteration_summary_.step_norm = (x_ - candidate_x_).norm();
+  const double step_size_tolerance =
+      options_.parameter_tolerance * (x_norm_ + options_.parameter_tolerance);
+
+  if (iteration_summary_.step_norm > step_size_tolerance) {
+    return false;
+  }
+
+  solver_summary_->message = StringPrintf(
+      "Parameter tolerance reached. "
+      "Relative step_norm: %e <= %e.",
+      (iteration_summary_.step_norm / (x_norm_ + options_.parameter_tolerance)),
+      options_.parameter_tolerance);
+  solver_summary_->termination_type = CONVERGENCE;
+  VLOG_IF(1, is_not_silent_) << "Terminating: " << solver_summary_->message;
+  return true;
+}
+
+// Solver::Options::function_tolerance based convergence check.
+bool TrustRegionMinimizer::FunctionToleranceReached() {
+  iteration_summary_.cost_change = x_cost_ - candidate_cost_;
+  const double absolute_function_tolerance =
+      options_.function_tolerance * x_cost_;
+
+  if (fabs(iteration_summary_.cost_change) > absolute_function_tolerance) {
+    return false;
+  }
+
+  solver_summary_->message = StringPrintf(
+      "Function tolerance reached. "
+      "|cost_change|/cost: %e <= %e",
+      fabs(iteration_summary_.cost_change) / x_cost_,
+      options_.function_tolerance);
+  solver_summary_->termination_type = CONVERGENCE;
+  VLOG_IF(1, is_not_silent_) << "Terminating: " << solver_summary_->message;
+  return true;
+}
+
+// Compute candidate_x_ = Plus(x_, delta_)
+// Evaluate the cost of candidate_x_ as candidate_cost_.
+//
+// Failure to compute the step or the cost mean that candidate_cost_
+// is set to std::numeric_limits<double>::max(). Unlike
+// EvaluateGradientAndJacobian, failure in this function is not fatal
+// as we are only computing and evaluating a candidate point, and if
+// for some reason we are unable to evaluate it, we consider it to be
+// a point with very high cost. This allows the user to deal with edge
+// cases/constraints as part of the LocalParameterization and
+// CostFunction objects.
+void TrustRegionMinimizer::ComputeCandidatePointAndEvaluateCost() {
+  if (!evaluator_->Plus(x_.data(), delta_.data(), candidate_x_.data())) {
+    LOG_IF(WARNING, is_not_silent_)
+        << "x_plus_delta = Plus(x, delta) failed. "
+        << "Treating it as a step with infinite cost";
+    candidate_cost_ = std::numeric_limits<double>::max();
+    return;
+  }
+
+  if (!evaluator_->Evaluate(
+          candidate_x_.data(), &candidate_cost_, NULL, NULL, NULL)) {
+    LOG_IF(WARNING, is_not_silent_)
+        << "Step failed to evaluate. "
+        << "Treating it as a step with infinite cost";
+    candidate_cost_ = std::numeric_limits<double>::max();
+  }
+}
+
+bool TrustRegionMinimizer::IsStepSuccessful() {
+  iteration_summary_.relative_decrease =
+      step_evaluator_->StepQuality(candidate_cost_, model_cost_change_);
+
+  // In most cases, boosting the model_cost_change by the
+  // improvement caused by the inner iterations is fine, but it can
+  // be the case that the original trust region step was so bad that
+  // the resulting improvement in the cost was negative, and the
+  // change caused by the inner iterations was large enough to
+  // improve the step, but also to make relative decrease quite
+  // small.
   //
-  // TODO(sameeragarwal): Make line search available more generally.
-  const bool use_line_search = options.is_constrained;
-
-  summary->termination_type = NO_CONVERGENCE;
-  summary->num_successful_steps = 0;
-  summary->num_unsuccessful_steps = 0;
-  summary->is_constrained = options.is_constrained;
-
-  const int num_parameters = evaluator->NumParameters();
-  const int num_effective_parameters = evaluator->NumEffectiveParameters();
-  const int num_residuals = evaluator->NumResiduals();
-
-  Vector residuals(num_residuals);
-  Vector trust_region_step(num_effective_parameters);
-  Vector delta(num_effective_parameters);
-  Vector x_plus_delta(num_parameters);
-  Vector gradient(num_effective_parameters);
-  Vector model_residuals(num_residuals);
-  Vector scale(num_effective_parameters);
-  Vector negative_gradient(num_effective_parameters);
-  Vector projected_gradient_step(num_parameters);
-
-  IterationSummary iteration_summary;
-  iteration_summary.iteration = 0;
-  iteration_summary.step_is_valid = false;
-  iteration_summary.step_is_successful = false;
-  iteration_summary.cost_change = 0.0;
-  iteration_summary.gradient_max_norm = 0.0;
-  iteration_summary.gradient_norm = 0.0;
-  iteration_summary.step_norm = 0.0;
-  iteration_summary.relative_decrease = 0.0;
-  iteration_summary.trust_region_radius = strategy->Radius();
-  iteration_summary.eta = options_.eta;
-  iteration_summary.linear_solver_iterations = 0;
-  iteration_summary.step_solver_time_in_seconds = 0;
-
-  VectorRef x_min(parameters, num_parameters);
-  Vector x = x_min;
-  // Project onto the feasible set.
-  if (options.is_constrained) {
-    delta.setZero();
-    if (!evaluator->Plus(x.data(), delta.data(), x_plus_delta.data())) {
-      summary->message =
-          "Unable to project initial point onto the feasible set.";
-      summary->termination_type = FAILURE;
-      LOG_IF(WARNING, is_not_silent) << "Terminating: " << summary->message;
-      return;
-    }
-    x_min = x_plus_delta;
-    x = x_plus_delta;
-  }
-
-  double x_norm = x.norm();
-
-  // Do initial cost and Jacobian evaluation.
-  double cost = 0.0;
-  if (!evaluator->Evaluate(x.data(),
-                           &cost,
-                           residuals.data(),
-                           gradient.data(),
-                           jacobian)) {
-    summary->message = "Residual and Jacobian evaluation failed.";
-    summary->termination_type = FAILURE;
-    LOG_IF(WARNING, is_not_silent) << "Terminating: " << summary->message;
-    return;
-  }
-
-  negative_gradient = -gradient;
-  if (!evaluator->Plus(x.data(),
-                       negative_gradient.data(),
-                       projected_gradient_step.data())) {
-    summary->message = "Unable to compute gradient step.";
-    summary->termination_type = FAILURE;
-    LOG(ERROR) << "Terminating: " << summary->message;
-    return;
-  }
-
-  summary->initial_cost = cost + summary->fixed_cost;
-  iteration_summary.cost = cost + summary->fixed_cost;
-  iteration_summary.gradient_max_norm =
-    (x - projected_gradient_step).lpNorm<Eigen::Infinity>();
-  iteration_summary.gradient_norm = (x - projected_gradient_step).norm();
-
-  if (iteration_summary.gradient_max_norm <= options.gradient_tolerance) {
-    summary->message = StringPrintf("Gradient tolerance reached. "
-                                    "Gradient max norm: %e <= %e",
-                                    iteration_summary.gradient_max_norm,
-                                    options_.gradient_tolerance);
-    summary->termination_type = CONVERGENCE;
-    VLOG_IF(1, is_not_silent) << "Terminating: " << summary->message;
-
-    // Ensure that there is an iteration summary object for iteration
-    // 0 in Summary::iterations.
-    iteration_summary.iteration_time_in_seconds =
-        WallTimeInSeconds() - iteration_start_time;
-    iteration_summary.cumulative_time_in_seconds =
-        WallTimeInSeconds() - start_time +
-        summary->preprocessor_time_in_seconds;
-    summary->iterations.push_back(iteration_summary);
-    return;
-  }
-
-  if (options_.jacobi_scaling) {
-    EstimateScale(*jacobian, scale.data());
-    jacobian->ScaleColumns(scale.data());
-  } else {
-    scale.setOnes();
-  }
-
-  iteration_summary.iteration_time_in_seconds =
-      WallTimeInSeconds() - iteration_start_time;
-  iteration_summary.cumulative_time_in_seconds =
-      WallTimeInSeconds() - start_time
-      + summary->preprocessor_time_in_seconds;
-  summary->iterations.push_back(iteration_summary);
-
-  int num_consecutive_nonmonotonic_steps = 0;
-  double minimum_cost = cost;
-  double reference_cost = cost;
-  double accumulated_reference_model_cost_change = 0.0;
-  double candidate_cost = cost;
-  double accumulated_candidate_model_cost_change = 0.0;
-  int num_consecutive_invalid_steps = 0;
-  bool inner_iterations_are_enabled =
-      options.inner_iteration_minimizer.get() != NULL;
-  while (true) {
-    bool inner_iterations_were_useful = false;
-    if (!RunCallbacks(options, iteration_summary, summary)) {
-      return;
-    }
-
-    iteration_start_time = WallTimeInSeconds();
-    if (iteration_summary.iteration >= options_.max_num_iterations) {
-      summary->message = "Maximum number of iterations reached.";
-      summary->termination_type = NO_CONVERGENCE;
-      VLOG_IF(1, is_not_silent) << "Terminating: " << summary->message;
-      return;
-    }
-
-    const double total_solver_time = iteration_start_time - start_time +
-        summary->preprocessor_time_in_seconds;
-    if (total_solver_time >= options_.max_solver_time_in_seconds) {
-      summary->message = "Maximum solver time reached.";
-      summary->termination_type = NO_CONVERGENCE;
-      VLOG_IF(1, is_not_silent) << "Terminating: " << summary->message;
-      return;
-    }
-
-    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());
-
-    if (strategy_summary.termination_type == LINEAR_SOLVER_FATAL_ERROR) {
-      summary->message =
-          "Linear solver failed due to unrecoverable "
-          "non-numeric causes. Please see the error log for clues. ";
-      summary->termination_type = FAILURE;
-      LOG_IF(WARNING, is_not_silent) << "Terminating: " << summary->message;
-      return;
-    }
-
-    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;
-    iteration_summary.step_is_valid = false;
-    iteration_summary.step_is_successful = false;
-
-    double model_cost_change = 0.0;
-    if (strategy_summary.termination_type != LINEAR_SOLVER_FAILURE) {
-      // new_model_cost
-      //  = 1/2 [f + J * step]^2
-      //  = 1/2 [ f'f + 2f'J * step + step' * J' * J * step ]
-      // model_cost_change
-      //  = cost - new_model_cost
-      //  = f'f/2  - 1/2 [ f'f + 2f'J * step + step' * J' * J * step]
-      //  = -f'J * step - step' * J' * J * step / 2
-      model_residuals.setZero();
-      jacobian->RightMultiply(trust_region_step.data(), model_residuals.data());
-      model_cost_change =
-          - model_residuals.dot(residuals + model_residuals / 2.0);
-
-      if (model_cost_change < 0.0) {
-        VLOG_IF(1, is_not_silent)
-            << "Invalid step: current_cost: " << cost
-            << " absolute difference " << model_cost_change
-            << " relative difference " << (model_cost_change / cost);
-      } else {
-        iteration_summary.step_is_valid = true;
-      }
-    }
-
-    if (!iteration_summary.step_is_valid) {
-      // Invalid steps can happen due to a number of reasons, and we
-      // allow a limited number of successive failures, and return with
-      // FAILURE if this limit is exceeded.
-      if (++num_consecutive_invalid_steps >=
-          options_.max_num_consecutive_invalid_steps) {
-        summary->message = StringPrintf(
-            "Number of successive invalid steps more "
-            "than Solver::Options::max_num_consecutive_invalid_steps: %d",
-            options_.max_num_consecutive_invalid_steps);
-        summary->termination_type = FAILURE;
-        LOG_IF(WARNING, is_not_silent) << "Terminating: " << summary->message;
-        return;
-      }
-
-      // We are going to try and reduce the trust region radius and
-      // solve again. To do this, we are going to treat this iteration
-      // as an unsuccessful iteration. Since the various callbacks are
-      // still executed, we are going to fill the iteration summary
-      // with data that assumes a step of length zero and no progress.
-      iteration_summary.cost = cost + summary->fixed_cost;
-      iteration_summary.cost_change = 0.0;
-      iteration_summary.gradient_max_norm =
-          summary->iterations.back().gradient_max_norm;
-      iteration_summary.gradient_norm =
-          summary->iterations.back().gradient_norm;
-      iteration_summary.step_norm = 0.0;
-      iteration_summary.relative_decrease = 0.0;
-      iteration_summary.eta = options_.eta;
-    } else {
-      // The step is numerically valid, so now we can judge its quality.
-      num_consecutive_invalid_steps = 0;
-
-      // Undo the Jacobian column scaling.
-      delta = (trust_region_step.array() * scale.array()).matrix();
-
-      // Try improving the step further by using an ARMIJO line
-      // search.
-      //
-      // TODO(sameeragarwal): What happens to trust region sizing as
-      // it interacts with the line search ?
-      if (use_line_search) {
-        const LineSearch::Summary line_search_summary =
-            DoLineSearch(options, x, gradient, cost, delta, evaluator);
-
-        // Iterations inside the line search algorithm are considered
-        // 'steps' in the broader context, to distinguish these inner
-        // iterations from from the outer iterations of the trust
-        // region minimizer The number of line search steps is the
-        // total number of inner line search iterations (or steps)
-        // across the entire minimization.
-        summary->num_line_search_steps += line_search_summary.num_iterations;
-        summary->line_search_cost_evaluation_time_in_seconds +=
-            line_search_summary.cost_evaluation_time_in_seconds;
-        summary->line_search_gradient_evaluation_time_in_seconds +=
-            line_search_summary.gradient_evaluation_time_in_seconds;
-        summary->line_search_polynomial_minimization_time_in_seconds +=
-            line_search_summary.polynomial_minimization_time_in_seconds;
-        summary->line_search_total_time_in_seconds +=
-            line_search_summary.total_time_in_seconds;
-
-        if (line_search_summary.success) {
-          delta *= line_search_summary.optimal_step_size;
-        }
-      }
-
-      double new_cost = std::numeric_limits<double>::max();
-      if (evaluator->Plus(x.data(), delta.data(), x_plus_delta.data())) {
-        if (!evaluator->Evaluate(x_plus_delta.data(),
-                                 &new_cost,
-                                 NULL,
-                                 NULL,
-                                 NULL)) {
-          LOG_IF(WARNING, is_not_silent)
-              << "Step failed to evaluate. "
-              << "Treating it as a step with infinite cost";
-          new_cost = std::numeric_limits<double>::max();
-        }
-      } else {
-        LOG_IF(WARNING, is_not_silent)
-            << "x_plus_delta = Plus(x, delta) failed. "
-            << "Treating it as a step with infinite cost";
-      }
-
-      if (new_cost < std::numeric_limits<double>::max()) {
-        // Check if performing an inner iteration will make it better.
-        if (inner_iterations_are_enabled) {
-          ++summary->num_inner_iteration_steps;
-          double inner_iteration_start_time = WallTimeInSeconds();
-          const double x_plus_delta_cost = new_cost;
-          Vector inner_iteration_x = x_plus_delta;
-          Solver::Summary inner_iteration_summary;
-          options.inner_iteration_minimizer->Minimize(options,
-                                                      inner_iteration_x.data(),
-                                                      &inner_iteration_summary);
-          if (!evaluator->Evaluate(inner_iteration_x.data(),
-                                   &new_cost,
-                                   NULL, NULL, NULL)) {
-            VLOG_IF(2, is_not_silent) << "Inner iteration failed.";
-            new_cost = x_plus_delta_cost;
-          } else {
-            x_plus_delta = inner_iteration_x;
-            // Boost the model_cost_change, since the inner iteration
-            // improvements are not accounted for by the trust region.
-            model_cost_change +=  x_plus_delta_cost - new_cost;
-            VLOG_IF(2, is_not_silent)
-                << "Inner iteration succeeded; Current cost: " << cost
-                << " Trust region step cost: " << x_plus_delta_cost
-                << " Inner iteration cost: " << new_cost;
-
-            inner_iterations_were_useful = new_cost < cost;
-
-            const double inner_iteration_relative_progress =
-                1.0 - new_cost / x_plus_delta_cost;
-            // Disable inner iterations once the relative improvement
-            // drops below tolerance.
-            inner_iterations_are_enabled =
-                (inner_iteration_relative_progress >
-                 options.inner_iteration_tolerance);
-            VLOG_IF(2, is_not_silent && !inner_iterations_are_enabled)
-                << "Disabling inner iterations. Progress : "
-                << inner_iteration_relative_progress;
-          }
-          summary->inner_iteration_time_in_seconds +=
-              WallTimeInSeconds() - inner_iteration_start_time;
-        }
-      }
-
-      iteration_summary.step_norm = (x - x_plus_delta).norm();
-
-      // Convergence based on parameter_tolerance.
-      const double step_size_tolerance =  options_.parameter_tolerance *
-          (x_norm + options_.parameter_tolerance);
-      if (iteration_summary.step_norm <= step_size_tolerance) {
-        summary->message =
-            StringPrintf("Parameter tolerance reached. "
-                         "Relative step_norm: %e <= %e.",
-                         (iteration_summary.step_norm /
-                          (x_norm + options_.parameter_tolerance)),
-                         options_.parameter_tolerance);
-        summary->termination_type = CONVERGENCE;
-        VLOG_IF(1, is_not_silent) << "Terminating: " << summary->message;
-        return;
-      }
-
-      iteration_summary.cost_change =  cost - new_cost;
-      const double absolute_function_tolerance =
-          options_.function_tolerance * cost;
-      if (fabs(iteration_summary.cost_change) <= absolute_function_tolerance) {
-        summary->message =
-            StringPrintf("Function tolerance reached. "
-                         "|cost_change|/cost: %e <= %e",
-                         fabs(iteration_summary.cost_change) / cost,
-                         options_.function_tolerance);
-        summary->termination_type = CONVERGENCE;
-        VLOG_IF(1, is_not_silent) << "Terminating: " << summary->message;
-        return;
-      }
-
-      const double relative_decrease =
-          iteration_summary.cost_change / model_cost_change;
-
-      const double historical_relative_decrease =
-          (reference_cost - new_cost) /
-          (accumulated_reference_model_cost_change + model_cost_change);
-
-      // If monotonic steps are being used, then the relative_decrease
-      // is the usual ratio of the change in objective function value
-      // divided by the change in model cost.
-      //
-      // If non-monotonic steps are allowed, then we take the maximum
-      // of the relative_decrease and the
-      // historical_relative_decrease, which measures the increase
-      // from a reference iteration. The model cost change is
-      // estimated by accumulating the model cost changes since the
-      // reference iteration. The historical relative_decrease offers
-      // a boost to a step which is not too bad compared to the
-      // reference iteration, allowing for non-monotonic steps.
-      iteration_summary.relative_decrease =
-          options.use_nonmonotonic_steps
-          ? std::max(relative_decrease, historical_relative_decrease)
-          : relative_decrease;
-
-      // Normally, the quality of a trust region step is measured by
-      // the ratio
-      //
-      //              cost_change
-      //    r =    -----------------
-      //           model_cost_change
-      //
-      // All the change in the nonlinear objective is due to the trust
-      // region step so this ratio is a good measure of the quality of
-      // the trust region radius. However, when inner iterations are
-      // being used, cost_change includes the contribution of the
-      // inner iterations and its not fair to credit it all to the
-      // trust region algorithm. So we change the ratio to be
-      //
-      //                              cost_change
-      //    r =    ------------------------------------------------
-      //           (model_cost_change + inner_iteration_cost_change)
-      //
-      // In most cases this is fine, but it can be the case that the
-      // change in solution quality due to inner iterations is so large
-      // and the trust region step is so bad, that this ratio can become
-      // quite small.
-      //
-      // This can cause the trust region loop to reject this step. To
-      // get around this, we expicitly check if the inner iterations
-      // led to a net decrease in the objective function value. If
-      // they did, we accept the step even if the trust region ratio
-      // is small.
-      //
-      // Notice that we do not just check that cost_change is positive
-      // which is a weaker condition and would render the
-      // min_relative_decrease threshold useless. Instead, we keep
-      // track of inner_iterations_were_useful, which is true only
-      // when inner iterations lead to a net decrease in the cost.
-      iteration_summary.step_is_successful =
-          (inner_iterations_were_useful ||
-           iteration_summary.relative_decrease >
-           options_.min_relative_decrease);
-
-      if (iteration_summary.step_is_successful) {
-        accumulated_candidate_model_cost_change += model_cost_change;
-        accumulated_reference_model_cost_change += model_cost_change;
-
-        if (!inner_iterations_were_useful &&
-            relative_decrease <= options_.min_relative_decrease) {
-          iteration_summary.step_is_nonmonotonic = true;
-          VLOG_IF(2, is_not_silent)
-              << "Non-monotonic step! "
-              << " relative_decrease: "
-              << relative_decrease
-              << " historical_relative_decrease: "
-              << historical_relative_decrease;
-        }
-      }
-    }
-
-    if (iteration_summary.step_is_successful) {
-      ++summary->num_successful_steps;
-      strategy->StepAccepted(iteration_summary.relative_decrease);
-
-      x = x_plus_delta;
-      x_norm = x.norm();
-
-      // Step looks good, evaluate the residuals and Jacobian at this
-      // point.
-      if (!evaluator->Evaluate(x.data(),
-                               &cost,
-                               residuals.data(),
-                               gradient.data(),
-                               jacobian)) {
-        summary->message = "Residual and Jacobian evaluation failed.";
-        summary->termination_type = FAILURE;
-        LOG_IF(WARNING, is_not_silent) << "Terminating: " << summary->message;
-        return;
-      }
-
-      negative_gradient = -gradient;
-      if (!evaluator->Plus(x.data(),
-                           negative_gradient.data(),
-                           projected_gradient_step.data())) {
-        summary->message =
-            "projected_gradient_step = Plus(x, -gradient) failed.";
-        summary->termination_type = FAILURE;
-        LOG(ERROR) << "Terminating: " << summary->message;
-        return;
-      }
-
-      iteration_summary.gradient_max_norm =
-          (x - projected_gradient_step).lpNorm<Eigen::Infinity>();
-      iteration_summary.gradient_norm = (x - projected_gradient_step).norm();
-
-      if (options_.jacobi_scaling) {
-        jacobian->ScaleColumns(scale.data());
-      }
-
-      // Update the best, reference and candidate iterates.
-      //
-      // Based on algorithm 10.1.2 (page 357) of "Trust Region
-      // Methods" by Conn Gould & Toint, or equations 33-40 of
-      // "Non-monotone trust-region algorithms for nonlinear
-      // optimization subject to convex constraints" by Phil Toint,
-      // Mathematical Programming, 77, 1997.
-      if (cost < minimum_cost) {
-        // A step that improves solution quality was found.
-        x_min = x;
-        minimum_cost = cost;
-        // Set the candidate iterate to the current point.
-        candidate_cost = cost;
-        num_consecutive_nonmonotonic_steps = 0;
-        accumulated_candidate_model_cost_change = 0.0;
-      } else {
-        ++num_consecutive_nonmonotonic_steps;
-        if (cost > candidate_cost) {
-          // The current iterate is has a higher cost than the
-          // candidate iterate. Set the candidate to this point.
-          VLOG_IF(2, is_not_silent)
-              << "Updating the candidate iterate to the current point.";
-          candidate_cost = cost;
-          accumulated_candidate_model_cost_change = 0.0;
-        }
-
-        // At this point we have made too many non-monotonic steps and
-        // we are going to reset the value of the reference iterate so
-        // as to force the algorithm to descend.
-        //
-        // This is the case because the candidate iterate has a value
-        // greater than minimum_cost but smaller than the reference
-        // iterate.
-        if (num_consecutive_nonmonotonic_steps ==
-            options.max_consecutive_nonmonotonic_steps) {
-          VLOG_IF(2, is_not_silent)
-              << "Resetting the reference point to the candidate point";
-          reference_cost = candidate_cost;
-          accumulated_reference_model_cost_change =
-              accumulated_candidate_model_cost_change;
-        }
-      }
-    } else {
-      ++summary->num_unsuccessful_steps;
-      if (iteration_summary.step_is_valid) {
-        strategy->StepRejected(iteration_summary.relative_decrease);
-      } else {
-        strategy->StepIsInvalid();
-      }
-    }
-
-    iteration_summary.cost = cost + summary->fixed_cost;
-    iteration_summary.trust_region_radius = strategy->Radius();
-    iteration_summary.iteration_time_in_seconds =
-        WallTimeInSeconds() - iteration_start_time;
-    iteration_summary.cumulative_time_in_seconds =
-        WallTimeInSeconds() - start_time
-        + summary->preprocessor_time_in_seconds;
-    summary->iterations.push_back(iteration_summary);
-
-    // If the step was successful, check for the gradient norm
-    // collapsing to zero, and if the step is unsuccessful then check
-    // if the trust region radius has collapsed to zero.
-    //
-    // For correctness (Number of IterationSummary objects, correct
-    // final cost, and state update) these convergence tests need to
-    // be performed at the end of the iteration.
-    if (iteration_summary.step_is_successful) {
-      // Gradient norm can only go down in successful steps.
-      if (iteration_summary.gradient_max_norm <= options.gradient_tolerance) {
-        summary->message = StringPrintf("Gradient tolerance reached. "
-                                        "Gradient max norm: %e <= %e",
-                                        iteration_summary.gradient_max_norm,
-                                        options_.gradient_tolerance);
-        summary->termination_type = CONVERGENCE;
-        VLOG_IF(1, is_not_silent) << "Terminating: " << summary->message;
-        return;
-      }
-    } else {
-      // Trust region radius can only go down if the step if
-      // unsuccessful.
-      if (iteration_summary.trust_region_radius <
-          options_.min_trust_region_radius) {
-        summary->message = "Termination. Minimum trust region radius reached.";
-        summary->termination_type = CONVERGENCE;
-        VLOG_IF(1, is_not_silent) << summary->message;
-        return;
-      }
-    }
-  }
+  // This can cause the trust region loop to reject this step. To
+  // get around this, we expicitly check if the inner iterations
+  // led to a net decrease in the objective function value. If
+  // they did, we accept the step even if the trust region ratio
+  // is small.
+  //
+  // Notice that we do not just check that cost_change is positive
+  // which is a weaker condition and would render the
+  // min_relative_decrease threshold useless. Instead, we keep
+  // track of inner_iterations_were_useful, which is true only
+  // when inner iterations lead to a net decrease in the cost.
+  return (inner_iterations_were_useful_ ||
+          iteration_summary_.relative_decrease >
+              options_.min_relative_decrease);
 }
 
+bool TrustRegionMinimizer::HandleSuccessfulStep() {
+  x_ = candidate_x_;
+  x_norm_ = x_.norm();
+
+  if (!EvaluateGradientAndJacobian()) {
+    return false;
+  }
+
+  iteration_summary_.step_is_successful = true;
+  strategy_->StepAccepted(iteration_summary_.relative_decrease);
+  step_evaluator_->StepAccepted(candidate_cost_, model_cost_change_);
+  return true;
+}
+
+void TrustRegionMinimizer::HandleUnsuccessfulStep() {
+  iteration_summary_.step_is_successful = false;
+  strategy_->StepRejected(iteration_summary_.relative_decrease);
+  iteration_summary_.cost = candidate_cost_ + solver_summary_->fixed_cost;
+}
 
 }  // namespace internal
 }  // namespace ceres
diff --git a/internal/ceres/trust_region_minimizer.h b/internal/ceres/trust_region_minimizer.h
index ed52c26..ac4a6ed 100644
--- a/internal/ceres/trust_region_minimizer.h
+++ b/internal/ceres/trust_region_minimizer.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2016 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -31,35 +31,106 @@
 #ifndef CERES_INTERNAL_TRUST_REGION_MINIMIZER_H_
 #define CERES_INTERNAL_TRUST_REGION_MINIMIZER_H_
 
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/scoped_ptr.h"
 #include "ceres/minimizer.h"
 #include "ceres/solver.h"
+#include "ceres/sparse_matrix.h"
+#include "ceres/trust_region_step_evaluator.h"
+#include "ceres/trust_region_strategy.h"
 #include "ceres/types.h"
 
 namespace ceres {
 namespace internal {
 
-// Generic trust region minimization algorithm. The heavy lifting is
-// done by a TrustRegionStrategy object passed in as part of options.
+// Generic trust region minimization algorithm.
 //
 // For example usage, see SolverImpl::Minimize.
 class TrustRegionMinimizer : public Minimizer {
  public:
-  ~TrustRegionMinimizer() {}
+  ~TrustRegionMinimizer();
+
+  // This method is not thread safe.
   virtual void Minimize(const Minimizer::Options& options,
                         double* parameters,
-                        Solver::Summary* summary);
+                        Solver::Summary* solver_summary);
 
  private:
-  void Init(const Minimizer::Options& options);
-  void EstimateScale(const SparseMatrix& jacobian, double* scale) const;
-  bool MaybeDumpLinearLeastSquaresProblem(const int iteration,
-                                          const SparseMatrix* jacobian,
-                                          const double* residuals,
-                                          const double* step) const;
+  void Init(const Minimizer::Options& options,
+            double* parameters,
+            Solver::Summary* solver_summary);
+  bool IterationZero();
+  bool FinalizeIterationAndCheckIfMinimizerCanContinue();
+  bool ComputeTrustRegionStep();
+
+  bool EvaluateGradientAndJacobian();
+  void ComputeCandidatePointAndEvaluateCost();
+
+  void DoLineSearch(const Vector& x,
+                    const Vector& gradient,
+                    const double cost,
+                    Vector* delta);
+  void DoInnerIterationsIfNeeded();
+
+  bool ParameterToleranceReached();
+  bool FunctionToleranceReached();
+  bool GradientToleranceReached();
+  bool MaxSolverTimeReached();
+  bool MaxSolverIterationsReached();
+  bool MinTrustRegionRadiusReached();
+
+  bool IsStepSuccessful();
+  void HandleUnsuccessfulStep();
+  bool HandleSuccessfulStep();
+  bool HandleInvalidStep();
 
   Minimizer::Options options_;
+
+  // These pointers are shortcuts to objects passed to the
+  // TrustRegionMinimizer. The TrustRegionMinimizer does not own them.
+  double* parameters_;
+  Solver::Summary* solver_summary_;
+  Evaluator* evaluator_;
+  SparseMatrix* jacobian_;
+  TrustRegionStrategy* strategy_;
+
+  scoped_ptr<TrustRegionStepEvaluator> step_evaluator_;
+
+  bool is_not_silent_;
+  bool inner_iterations_are_enabled_;
+  bool inner_iterations_were_useful_;
+
+  // Summary of the current iteration.
+  IterationSummary iteration_summary_;
+
+  int num_parameters_;
+  int num_effective_parameters_;
+  int num_residuals_;
+
+  Vector delta_;
+  Vector gradient_;
+  Vector inner_iteration_x_;
+  Vector model_residuals_;
+  Vector negative_gradient_;
+  Vector projected_gradient_step_;
+  Vector residuals_;
+  Vector trust_region_step_;
+  Vector x_;
+  Vector candidate_x_;
+  Vector jacobian_scaling_;
+
+  double x_norm_;
+  double x_cost_;
+  double minimum_cost_;
+  double model_cost_change_;
+  double candidate_cost_;
+
+  double start_time_;
+  double iteration_start_time_;
+  int num_consecutive_invalid_steps_;
 };
 
 }  // namespace internal
 }  // namespace ceres
+
 #endif  // CERES_INTERNAL_TRUST_REGION_MINIMIZER_H_
diff --git a/internal/ceres/trust_region_step_evaluator.cc b/internal/ceres/trust_region_step_evaluator.cc
new file mode 100644
index 0000000..c9167e6
--- /dev/null
+++ b/internal/ceres/trust_region_step_evaluator.cc
@@ -0,0 +1,107 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2016 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include <algorithm>
+#include "ceres/trust_region_step_evaluator.h"
+#include "glog/logging.h"
+
+namespace ceres {
+namespace internal {
+
+TrustRegionStepEvaluator::TrustRegionStepEvaluator(
+    const double initial_cost,
+    const int max_consecutive_nonmonotonic_steps)
+    : max_consecutive_nonmonotonic_steps_(max_consecutive_nonmonotonic_steps),
+      minimum_cost_(initial_cost),
+      current_cost_(initial_cost),
+      reference_cost_(initial_cost),
+      candidate_cost_(initial_cost),
+      accumulated_reference_model_cost_change_(0.0),
+      accumulated_candidate_model_cost_change_(0.0),
+      num_consecutive_nonmonotonic_steps_(0){
+}
+
+double TrustRegionStepEvaluator::StepQuality(
+    const double cost,
+    const double model_cost_change) const {
+  const double relative_decrease = (current_cost_ - cost) / model_cost_change;
+  const double historical_relative_decrease =
+      (reference_cost_ - cost) /
+      (accumulated_reference_model_cost_change_ + model_cost_change);
+  return std::max(relative_decrease, historical_relative_decrease);
+}
+
+void TrustRegionStepEvaluator::StepAccepted(
+    const double cost,
+    const double model_cost_change) {
+  // Algorithm 10.1.2 from Trust Region Methods by Conn, Gould &
+  // Toint.
+  //
+  // Step 3a
+  current_cost_ = cost;
+  accumulated_candidate_model_cost_change_ += model_cost_change;
+  accumulated_reference_model_cost_change_ += model_cost_change;
+
+  // Step 3b.
+  if (current_cost_ < minimum_cost_) {
+    minimum_cost_ = current_cost_;
+    num_consecutive_nonmonotonic_steps_ = 0;
+    candidate_cost_ = current_cost_;
+    accumulated_candidate_model_cost_change_ = 0.0;
+  } else {
+    // Step 3c.
+    ++num_consecutive_nonmonotonic_steps_;
+    if (current_cost_ > candidate_cost_) {
+      candidate_cost_ = current_cost_;
+      accumulated_candidate_model_cost_change_ = 0.0;
+    }
+  }
+
+  // Step 3d.
+  //
+  // At this point we have made too many non-monotonic steps and
+  // we are going to reset the value of the reference iterate so
+  // as to force the algorithm to descend.
+  //
+  // Note: In the original algorithm by Toint, this step was only
+  // executed if the step was non-monotonic, but that would not handle
+  // the case of max_consecutive_nonmonotonic_steps = 0. The small
+  // modification of doing this always handles that corner case
+  // correctly.
+  if (num_consecutive_nonmonotonic_steps_ ==
+      max_consecutive_nonmonotonic_steps_) {
+    reference_cost_ = candidate_cost_;
+    accumulated_reference_model_cost_change_ =
+        accumulated_candidate_model_cost_change_;
+  }
+}
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/internal/ceres/trust_region_step_evaluator.h b/internal/ceres/trust_region_step_evaluator.h
new file mode 100644
index 0000000..1f14954
--- /dev/null
+++ b/internal/ceres/trust_region_step_evaluator.h
@@ -0,0 +1,98 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2016 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#ifndef CERES_INTERNAL_TRUST_REGION_STEP_EVALUATOR_H_
+#define CERES_INTERNAL_TRUST_REGION_STEP_EVALUATOR_H_
+
+namespace ceres {
+namespace internal {
+
+// The job of the TrustRegionStepEvaluator is to evaluate the quality
+// of a step, i.e., how the cost of a step compares with the reduction
+// in the objective of the trust region problem.
+//
+// Classic trust region methods are descent methods, in that they only
+// accept a point if it strictly reduces the value of the objective
+// function. They do this by measuring the quality of a step as
+//
+//   cost_change / model_cost_change.
+//
+// Relaxing the monotonic descent requirement allows the algorithm to
+// be more efficient in the long term at the cost of some local
+// increase in the value of the objective function.
+//
+// This is because allowing for non-decreasing objective function
+// values in a principaled manner allows the algorithm to "jump over
+// boulders" as the method is not restricted to move into narrow
+// valleys while preserving its convergence properties.
+//
+// The parameter max_consecutive_nonmonotonic_steps controls the
+// window size used by the step selection algorithm to accept
+// non-monotonic steps. Setting this parameter to zero, recovers the
+// classic montonic descent algorithm.
+//
+// Based on algorithm 10.1.2 (page 357) of "Trust Region
+// Methods" by Conn Gould & Toint, or equations 33-40 of
+// "Non-monotone trust-region algorithms for nonlinear
+// optimization subject to convex constraints" by Phil Toint,
+// Mathematical Programming, 77, 1997.
+//
+// Example usage:
+//
+// TrustRegionStepEvaluator* step_evaluator = ...
+//
+// cost = ... // Compute the non-linear objective function value.
+// model_cost_change = ... // Change in the value of the trust region objective.
+// if (step_evaluator->StepQuality(cost, model_cost_change) > threshold) {
+//   x = x + delta;
+//   step_evaluator->StepAccepted(cost, model_cost_change);
+// }
+class TrustRegionStepEvaluator {
+ public:
+  TrustRegionStepEvaluator(double initial_cost,
+                           int max_consecutive_nonmonotonic_steps);
+  double StepQuality(double cost, double model_cost_change) const;
+  void StepAccepted(double cost, double model_cost_change);
+
+ private:
+  const int max_consecutive_nonmonotonic_steps_;
+  double minimum_cost_;
+  double current_cost_;
+  double reference_cost_;
+  double candidate_cost_;
+  double accumulated_reference_model_cost_change_;
+  double accumulated_candidate_model_cost_change_;
+  int num_consecutive_nonmonotonic_steps_;
+};
+
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // CERES_INTERNAL_TRUST_REGION_STEP_EVALUATOR_H_
diff --git a/internal/ceres/trust_region_strategy.h b/internal/ceres/trust_region_strategy.h
index 9560e67..36e8e98 100644
--- a/internal/ceres/trust_region_strategy.h
+++ b/internal/ceres/trust_region_strategy.h
@@ -86,20 +86,20 @@
   struct PerSolveOptions {
     PerSolveOptions()
         : eta(0),
-          dump_filename_base(""),
           dump_format_type(TEXTFILE) {
     }
 
     // Forcing sequence for inexact solves.
     double eta;
 
+    DumpFormatType dump_format_type;
+
     // 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.
     std::string dump_filename_base;
-    DumpFormatType dump_format_type;
   };
 
   struct Summary {