diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
index 5e1b3a6..bb8d3b4 100644
--- a/examples/bundle_adjuster.cc
+++ b/examples/bundle_adjuster.cc
@@ -86,6 +86,8 @@
 DEFINE_bool(use_block_amd, true, "Use a block oriented fill reducing ordering.");
 DEFINE_string(trust_region_strategy, "lm", "Options are: lm, dogleg");
 DEFINE_double(max_solver_time, 1e32, "Maximum solve time in seconds.");
+DEFINE_bool(nonmonotonic_steps, false, "Trust region algorithm can use"
+            " nonmonotic steps");
 
 namespace ceres {
 namespace examples {
@@ -221,6 +223,7 @@
   options->num_threads = FLAGS_num_threads;
   options->eta = FLAGS_eta;
   options->max_solver_time_in_seconds = FLAGS_max_solver_time;
+  options->use_nonmonotonic_steps = FLAGS_nonmonotonic_steps;
   if (FLAGS_trust_region_strategy == "lm") {
     options->trust_region_strategy_type = LEVENBERG_MARQUARDT;
   } else if (FLAGS_trust_region_strategy == "dogleg") {
diff --git a/include/ceres/solver.h b/include/ceres/solver.h
index 4d59bde..8a58cb1 100644
--- a/include/ceres/solver.h
+++ b/include/ceres/solver.h
@@ -58,6 +58,8 @@
     // Default constructor that sets up a generic sparse problem.
     Options() {
       trust_region_strategy_type = LEVENBERG_MARQUARDT;
+      use_nonmonotonic_steps = false;
+      max_consecutive_nonmonotonic_steps = 5;
       max_num_iterations = 50;
       max_solver_time_in_seconds = 1e9;
       num_threads = 1;
@@ -119,6 +121,34 @@
 
     TrustRegionStrategyType trust_region_strategy_type;
 
+    // The classical trust region methods are descent methods, in that
+    // they only accept a point if it strictly reduces the value of
+    // the objective function.
+    //
+    // Relaxing this 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 princpled manner allows the algorithm to "jump over
+    // boulders" as the method is not restricted to move into narrow
+    // valleys while preserving its convergence properties.
+    //
+    // Setting use_nonmonotonic_steps to true enables the
+    // non-monotonic trust region algorithm as described by Conn,
+    // Gould & Toint in "Trust Region Methods", Section 10.1.
+    //
+    // The parameter max_consecutive_nonmonotonic_steps controls the
+    // window size used by the step selection algorithm to accept
+    // non-monotonic steps.
+    //
+    // Even though the value of the objective function may be larger
+    // than the minimum value encountered over the course of the
+    // optimization, the final parameters returned to the user are the
+    // ones corresponding to the minimum cost over all iterations.
+    bool use_nonmonotonic_steps;
+    int max_consecutive_nonmonotonic_steps;
+
     // Maximum number of iterations for the minimizer to run for.
     int max_num_iterations;
 
diff --git a/internal/ceres/minimizer.h b/internal/ceres/minimizer.h
index 70b530f..28b76ce 100644
--- a/internal/ceres/minimizer.h
+++ b/internal/ceres/minimizer.h
@@ -68,6 +68,9 @@
       min_relative_decrease = options.min_relative_decrease;
       eta = options.eta;
       jacobi_scaling = options.jacobi_scaling;
+      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;
@@ -96,6 +99,8 @@
     double min_relative_decrease;
     double eta;
     bool jacobi_scaling;
+    bool use_nonmonotonic_steps;
+    bool max_consecutive_nonmonotonic_steps;
     vector<int> lsqp_iterations_to_dump;
     DumpFormatType lsqp_dump_format_type;
     string lsqp_dump_directory;
diff --git a/internal/ceres/trust_region_minimizer.cc b/internal/ceres/trust_region_minimizer.cc
index 357d822..a6cb830 100644
--- a/internal/ceres/trust_region_minimizer.cc
+++ b/internal/ceres/trust_region_minimizer.cc
@@ -132,7 +132,8 @@
   const int num_effective_parameters = evaluator->NumEffectiveParameters();
   const int num_residuals = evaluator->NumResiduals();
 
-  VectorRef x(parameters, num_parameters);
+  VectorRef x_min(parameters, num_parameters);
+  Vector x = x_min;
   double x_norm = x.norm();
 
   Vector residuals(num_residuals);
@@ -145,8 +146,8 @@
 
   IterationSummary iteration_summary;
   iteration_summary.iteration = 0;
-  iteration_summary.step_is_valid=false;
-  iteration_summary.step_is_successful=false;
+  iteration_summary.step_is_valid = false;
+  iteration_summary.step_is_successful = false;
   iteration_summary.cost = summary->initial_cost;
   iteration_summary.cost_change = 0.0;
   iteration_summary.gradient_max_norm = 0.0;
@@ -167,6 +168,13 @@
     return;
   }
 
+  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;
+
   gradient.setZero();
   jacobian->LeftMultiply(residuals.data(), gradient.data());
   iteration_summary.gradient_max_norm = gradient.lpNorm<Eigen::Infinity>();
@@ -366,10 +374,43 @@
         return;
       }
 
-      iteration_summary.relative_decrease =
+      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
+          ? max(relative_decrease, historical_relative_decrease)
+          : relative_decrease;
+
       iteration_summary.step_is_successful =
           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 (relative_decrease <= options_.min_relative_decrease) {
+          VLOG(2) << "Non-monotonic step! "
+                  << " relative_decrease: " << relative_decrease
+                  << " historical_relative_decrease: "
+                  << historical_relative_decrease;
+        }
+      }
     }
 
     if (iteration_summary.step_is_successful) {
@@ -377,6 +418,7 @@
       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(),
@@ -405,6 +447,47 @@
       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(2) << "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(2) << "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) {