diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
index 0331f32..718fde9 100644
--- a/examples/bundle_adjuster.cc
+++ b/examples/bundle_adjuster.cc
@@ -64,15 +64,12 @@
 #include "ceres/ceres.h"
 
 DEFINE_string(input, "", "Input File name");
-
 DEFINE_string(solver_type, "sparse_schur", "Options are: "
               "sparse_schur, dense_schur, iterative_schur, cholesky, "
               "dense_qr, and conjugate_gradients");
-
 DEFINE_string(preconditioner_type, "jacobi", "Options are: "
               "identity, jacobi, schur_jacobi, cluster_jacobi, "
               "cluster_tridiagonal");
-
 DEFINE_int32(num_iterations, 5, "Number of iterations");
 DEFINE_int32(num_threads, 1, "Number of threads");
 DEFINE_double(eta, 1e-2, "Default value for eta. Eta determines the "
@@ -201,9 +198,9 @@
 
 void SetSolverOptionsFromFlags(BALProblem* bal_problem,
                                Solver::Options* options) {
+  SetMinimizerOptions(options);
   SetLinearSolver(options);
   SetOrdering(bal_problem, options);
-  SetMinimizerOptions(options);
 }
 
 void BuildProblem(BALProblem* bal_problem, Problem* problem) {
diff --git a/include/ceres/iteration_callback.h b/include/ceres/iteration_callback.h
index 88da992..6d7720c 100644
--- a/include/ceres/iteration_callback.h
+++ b/include/ceres/iteration_callback.h
@@ -28,8 +28,9 @@
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
 //
-// When an iteration callback is specified, Ceres calls the callback after each
-// optimizer step and pass it an IterationSummary object, defined below.
+// When an iteration callback is specified, Ceres calls the callback
+// after each minimizer step (if the minimizer has not converged) and
+// passes it an IterationSummary object, defined below.
 
 #ifndef CERES_PUBLIC_ITERATION_CALLBACK_H_
 #define CERES_PUBLIC_ITERATION_CALLBACK_H_
@@ -44,7 +45,15 @@
   // Current iteration number.
   int32 iteration;
 
+  // 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.
+  bool step_is_valid;
+
   // Whether or not the algorithm made progress in this iteration.
+  //
+  // Note: step_is_successful is false when iteration = 0.
   bool step_is_successful;
 
   // Value of the objective function.
@@ -66,9 +75,10 @@
   // cost and the change in the cost of the linearized approximation.
   double relative_decrease;
 
-  // Value of the regularization parameter for Levenberg-Marquardt
-  // algorithm at the end of the current iteration.
-  double mu;
+  // Size of the trust region at the end of the current iteration. For
+  // the Levenberg-Marquardt algorithm, the regularization parameter
+  // mu = 1.0 / trust_region_radius.
+  double trust_region_radius;
 
   // For the inexact step Levenberg-Marquardt algorithm, this is the
   // relative accuracy with which the Newton(LM) step is solved. This
@@ -81,13 +91,15 @@
   // Newton step.
   int linear_solver_iterations;
 
-  // TODO(sameeragarwal): Change to use a higher precision timer using
-  // clock_gettime.
-  // Time (in seconds) spent inside the linear least squares solver.
+  // TODO(sameeragarwal): Change the following two to use a higher
+  // precision timer using clock_gettime.
+  //
+  // Time (in seconds) spent inside the minimizer loop in the current
+  // iteration.
   int iteration_time_sec;
 
-  // Time (in seconds) spent inside the linear least squares solver.
-  int linear_solver_time_sec;
+  // Time (in seconds) spent inside the trust region step solver.
+  int step_solver_time_sec;
 };
 
 // Interface for specifying callbacks that are executed at the end of
@@ -133,7 +145,7 @@
 //                                    summary.gradient_max_norm,
 //                                    summary.step_norm,
 //                                    summary.relative_decrease,
-//                                    summary.mu,
+//                                    summary.trust_region_radius,
 //                                    summary.eta,
 //                                    summary.linear_solver_iterations);
 //       if (log_to_stdout_) {
diff --git a/include/ceres/solver.h b/include/ceres/solver.h
index 701f0b2..ed4c9b8 100644
--- a/include/ceres/solver.h
+++ b/include/ceres/solver.h
@@ -57,13 +57,17 @@
   struct Options {
     // Default constructor that sets up a generic sparse problem.
     Options() {
-      minimizer_type = LEVENBERG_MARQUARDT;
+      trust_region_strategy_type = LEVENBERG_MARQUARDT;
       max_num_iterations = 50;
-      max_solver_time_sec = 1.0e9;
+      max_solver_time_sec = 1e9;
       num_threads = 1;
-      tau = 1e-4;
-      min_mu = 1e-20;
+      initial_trust_region_radius = 1e4;
+      max_trust_region_radius = 1e16;
+      min_trust_region_radius = 1e-32;
       min_relative_decrease = 1e-3;
+      lm_min_diagonal = 1e-6;
+      lm_max_diagonal = 1e32;
+      max_num_consecutive_invalid_steps = 5;
       function_tolerance = 1e-6;
       gradient_tolerance = 1e-10;
       parameter_tolerance = 1e-8;
@@ -93,7 +97,6 @@
       return_final_residuals = false;
       lsqp_dump_directory = "/tmp";
       lsqp_dump_format_type = TEXTFILE;
-      crash_and_dump_lsqp_on_failure = false;
       check_gradients = false;
       gradient_check_relative_precision = 1e-8;
       numeric_derivative_relative_step_size = 1e-6;
@@ -102,7 +105,7 @@
 
     // Minimizer options ----------------------------------------
 
-    MinimizerType minimizer_type;
+    TrustRegionStrategyType trust_region_strategy_type;
 
     // Maximum number of iterations for the minimizer to run for.
     int max_num_iterations;
@@ -114,27 +117,35 @@
     // jacobians.
     int num_threads;
 
-    // For Levenberg-Marquardt, the initial value for the
-    // regularizer. This is the inversely related to the size of the
-    // initial trust region.
-    double tau;
+    // Trust region minimizer settings.
+    double initial_trust_region_radius;
+    double max_trust_region_radius;
 
-    // For Levenberg-Marquardt, the minimum value of the
-    // regularizer. For well constrained problems there shold be no
-    // need to modify the default value, but in some cases, going
-    // below a certain minimum reliably triggers rank deficiency in
-    // the normal equations. In such cases, the LM solver can
-    // oscillate between lowering the value of mu, seeing a numerical
-    // failure, and then increasing it making some progress and then
-    // reducing it again.
-    //
-    // In such cases, it is useful to set a higher value for min_mu.
-    double min_mu;
+    // Minimizer terminates when the trust region radius becomes
+    // smaller than this value.
+    double min_trust_region_radius;
 
-    // For trust region methods, this is lower threshold for the
-    // relative decrease before a step is accepted.
+    // Lower bound for the relative decrease before a step is
+    // accepted.
     double min_relative_decrease;
 
+    // For the Levenberg-Marquadt algorithm, the scaled diagonal of
+    // the normal equations J'J is used to control the size of the
+    // trust region. Extremely small and large values along the
+    // diagonal can make this regularization scheme
+    // fail. lm_max_diagonal and lm_min_diagonal, clamp the values of
+    // diag(J'J) from above and below. In the normal course of
+    // operation, the user should not have to modify these parameters.
+    double lm_min_diagonal;
+    double lm_max_diagonal;
+
+    // Sometimes due to numerical conditioning problems or linear
+    // solver flakiness, the trust region strategy may return a
+    // numerically invalid step that can be fixed by reducing the
+    // trust region size. So the TrustRegionMinimizer allows for a few
+    // successive invalid steps before it declares NUMERICAL_FAILURE.
+    int max_num_consecutive_invalid_steps;
+
     // Minimizer terminates when
     //
     //   (new_cost - old_cost) < function_tolerance * old_cost;
@@ -243,15 +254,6 @@
     string lsqp_dump_directory;
     DumpFormatType lsqp_dump_format_type;
 
-    // Dump the linear least squares problem to disk if the minimizer
-    // fails due to NUMERICAL_FAILURE and crash the process. This flag
-    // is useful for generating debugging information. The problem is
-    // dumped in a file whose name is determined by
-    // Solver::Options::lsqp_dump_format.
-    //
-    // Note: This requires a version of Ceres built with protocol buffers.
-    bool crash_and_dump_lsqp_on_failure;
-
     // Finite differences options ----------------------------------------------
 
     // Check all jacobians computed by each residual block with finite
@@ -299,10 +301,15 @@
     bool update_state_every_iteration;
 
     // Callbacks that are executed at the end of each iteration of the
-    // Minimizer. They are executed in the order that they are
-    // specified in this vector. By default, parameter blocks are
-    // updated only at the end of the optimization, i.e when the
-    // Minimizer terminates. This behaviour is controlled by
+    // Minimizer. An iteration may terminate midway, either due to
+    // numerical failures or because one of the convergence tests has
+    // been satisfied. In this case none of the callbacks are
+    // executed.
+
+    // Callbacks are executed in the order that they are specified in
+    // this vector. By default, parameter blocks are updated only at
+    // the end of the optimization, i.e when the Minimizer
+    // terminates. This behaviour is controlled by
     // update_state_every_variable. If the user wishes to have access
     // to the update parameter blocks when his/her callbacks are
     // executed, then set update_state_every_iteration to true.
diff --git a/include/ceres/types.h b/include/ceres/types.h
index 61a9a0d..705097f 100644
--- a/include/ceres/types.h
+++ b/include/ceres/types.h
@@ -158,8 +158,8 @@
   PER_MINIMIZER_ITERATION
 };
 
-enum MinimizerType {
-  LEVENBERG_MARQUARDT
+enum TrustRegionStrategyType {
+  LEVENBERG_MARQUARDT,
 };
 
 enum SolverTerminationType {
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 4cdf055..2a0a6af 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -53,7 +53,7 @@
     gradient_checking_cost_function.cc
     implicit_schur_complement.cc
     iterative_schur_complement_solver.cc
-    levenberg_marquardt.cc
+    levenberg_marquardt_strategy.cc
     linear_least_squares_problems.cc
     linear_operator.cc
     linear_solver.cc
@@ -78,6 +78,8 @@
     split.cc
     stringprintf.cc
     suitesparse.cc
+    trust_region_minimizer.cc
+    trust_region_strategy.cc
     triplet_sparse_matrix.cc
     types.cc
     visibility_based_preconditioner.cc
@@ -176,7 +178,7 @@
   CERES_TEST(implicit_schur_complement)
   CERES_TEST(iterative_schur_complement_solver)
   CERES_TEST(jet)
-  CERES_TEST(levenberg_marquardt)
+  CERES_TEST(levenberg_marquardt_strategy)
   CERES_TEST(local_parameterization)
   CERES_TEST(loss_function)
   CERES_TEST(normal_prior)
@@ -196,6 +198,7 @@
   CERES_TEST(solver_impl)
   CERES_TEST(symmetric_linear_solver)
   CERES_TEST(triplet_sparse_matrix)
+  CERES_TEST(trust_region_minimizer)
   CERES_TEST(unsymmetric_linear_solver)
   CERES_TEST(visibility)
   IF (GFLAGS)
diff --git a/internal/ceres/levenberg_marquardt.cc b/internal/ceres/levenberg_marquardt.cc
deleted file mode 100644
index 1eae624..0000000
--- a/internal/ceres/levenberg_marquardt.cc
+++ /dev/null
@@ -1,588 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
-// http://code.google.com/p/ceres-solver/
-//
-// 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)
-//
-// Implementation of a simple LM algorithm which uses the step sizing
-// rule of "Methods for Nonlinear Least Squares" by K. Madsen,
-// H.B. Nielsen and O. Tingleff. Available to download from
-//
-// http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
-//
-// The basic algorithm described in this note is an exact step
-// algorithm that depends on the Newton(LM) step being solved exactly
-// in each iteration. When a suitable iterative solver is available to
-// solve the Newton(LM) step, the algorithm will automatically switch
-// to an inexact step solution method. This trades some slowdown in
-// convergence for significant savings in solve time and memory
-// usage. Our implementation of the Truncated Newton algorithm follows
-// the discussion and recommendataions in "Stephen G. Nash, A Survey
-// of Truncated Newton Methods, Journal of Computational and Applied
-// Mathematics, 124(1-2), 45-59, 2000.
-
-#include "ceres/levenberg_marquardt.h"
-
-#include <algorithm>
-#include <cstdlib>
-#include <cmath>
-#include <cstring>
-#include <string>
-#include <vector>
-
-#include <glog/logging.h>
-#include "Eigen/Core"
-#include "ceres/array_utils.h"
-#include "ceres/evaluator.h"
-#include "ceres/file.h"
-#include "ceres/linear_least_squares_problems.h"
-#include "ceres/linear_solver.h"
-#include "ceres/matrix_proto.h"
-#include "ceres/sparse_matrix.h"
-#include "ceres/stringprintf.h"
-#include "ceres/internal/eigen.h"
-#include "ceres/internal/scoped_ptr.h"
-#include "ceres/types.h"
-
-namespace ceres {
-namespace internal {
-namespace {
-
-// Numbers for clamping the size of the LM diagonal. The size of these
-// numbers is heuristic. We will probably be adjusting them in the
-// future based on more numerical experience. With jacobi scaling
-// enabled, these numbers should be all but redundant.
-const double kMinLevenbergMarquardtDiagonal = 1e-6;
-const double kMaxLevenbergMarquardtDiagonal = 1e32;
-
-// Small constant for various floating point issues.
-const double kEpsilon = 1e-12;
-
-// Number of times the linear solver should be retried in case of
-// numerical failure. The retries are done by exponentially scaling up
-// mu at each retry. This leads to stronger and stronger
-// regularization making the linear least squares problem better
-// conditioned at each retry.
-const int kMaxLinearSolverRetries = 5;
-
-// D = 1/sqrt(diag(J^T * J))
-void EstimateScale(const SparseMatrix& jacobian, double* D) {
-  CHECK_NOTNULL(D);
-  jacobian.SquaredColumnNorm(D);
-  for (int i = 0; i < jacobian.num_cols(); ++i) {
-    D[i] = 1.0 / (kEpsilon + sqrt(D[i]));
-  }
-}
-
-// D = diag(J^T * J)
-void LevenbergMarquardtDiagonal(const SparseMatrix& jacobian,
-                                double* D) {
-  CHECK_NOTNULL(D);
-  jacobian.SquaredColumnNorm(D);
-  for (int i = 0; i < jacobian.num_cols(); ++i) {
-    D[i] = min(max(D[i], kMinLevenbergMarquardtDiagonal),
-               kMaxLevenbergMarquardtDiagonal);
-  }
-}
-
-bool RunCallback(IterationCallback* callback,
-                 const IterationSummary& iteration_summary,
-                 Solver::Summary* summary) {
-  const CallbackReturnType status = (*callback)(iteration_summary);
-  switch (status) {
-    case SOLVER_TERMINATE_SUCCESSFULLY:
-      summary->termination_type = USER_SUCCESS;
-      VLOG(1) << "Terminating on USER_SUCCESS.";
-      return false;
-    case SOLVER_ABORT:
-      summary->termination_type = USER_ABORT;
-      VLOG(1) << "Terminating on USER_ABORT.";
-      return false;
-    case SOLVER_CONTINUE:
-      return true;
-    default:
-      LOG(FATAL) << "Unknown status returned by callback: "
-                 << status;
-  }
-}
-
-}  // namespace
-
-LevenbergMarquardt::~LevenbergMarquardt() {}
-
-void LevenbergMarquardt::Minimize(const Minimizer::Options& options,
-                                  Evaluator* evaluator,
-                                  LinearSolver* linear_solver,
-                                  const double* initial_parameters,
-                                  double* final_parameters,
-                                  Solver::Summary* summary) {
-  time_t start_time = time(NULL);
-  const int num_parameters = evaluator->NumParameters();
-  const int num_effective_parameters = evaluator->NumEffectiveParameters();
-  const int num_residuals = evaluator->NumResiduals();
-
-  summary->termination_type = NO_CONVERGENCE;
-  summary->num_successful_steps = 0;
-  summary->num_unsuccessful_steps = 0;
-
-  // Allocate the various vectors needed by the algorithm.
-  memcpy(final_parameters, initial_parameters,
-         num_parameters * sizeof(*initial_parameters));
-
-  VectorRef x(final_parameters, num_parameters);
-  Vector x_new(num_parameters);
-
-  Vector lm_step(num_effective_parameters);
-  Vector gradient(num_effective_parameters);
-  Vector scaled_gradient(num_effective_parameters);
-  // Jacobi scaling vector
-  Vector scale(num_effective_parameters);
-
-  Vector f_model(num_residuals);
-  Vector f(num_residuals);
-  Vector f_new(num_residuals);
-  Vector D(num_parameters);
-  Vector muD(num_parameters);
-
-  // Ask the Evaluator to create the jacobian matrix. The sparsity
-  // pattern of this matrix is going to remain constant, so we only do
-  // this once and then re-use this matrix for all subsequent Jacobian
-  // computations.
-  scoped_ptr<SparseMatrix> jacobian(evaluator->CreateJacobian());
-
-  double x_norm  = x.norm();
-
-  double cost = 0.0;
-  D.setOnes();
-  f.setZero();
-
-  // Do initial cost and Jacobian evaluation.
-  if (!evaluator->Evaluate(x.data(), &cost, f.data(), jacobian.get())) {
-    LOG(WARNING) << "Failed to compute residuals and Jacobian. "
-                 << "Terminating.";
-    summary->termination_type = NUMERICAL_FAILURE;
-    return;
-  }
-
-  if (options.jacobi_scaling) {
-    EstimateScale(*jacobian, scale.data());
-    jacobian->ScaleColumns(scale.data());
-  } else {
-    scale.setOnes();
-  }
-
-  // This is a poor way to do this computation. Even if fixed_cost is
-  // zero, because we are subtracting two possibly large numbers, we
-  // are depending on exact cancellation to give us a zero here. But
-  // initial_cost and cost have been computed by two different
-  // evaluators. One which runs on the whole problem (in
-  // solver_impl.cc) in single threaded mode and another which runs
-  // here on the reduced problem, so fixed_cost can (and does) contain
-  // some numerical garbage with a relative magnitude of 1e-14.
-  //
-  // The right way to do this, would be to compute the fixed cost on
-  // just the set of residual blocks which are held constant and were
-  // removed from the original problem when the reduced problem was
-  // constructed.
-  summary->fixed_cost = summary->initial_cost - cost;
-
-  double model_cost = f.squaredNorm() / 2.0;
-  double total_cost = summary->fixed_cost + cost;
-
-  scaled_gradient.setZero();
-  jacobian->LeftMultiply(f.data(), scaled_gradient.data());
-  gradient = scaled_gradient.array() / scale.array();
-
-  double gradient_max_norm = gradient.lpNorm<Eigen::Infinity>();
-  // We need the max here to guard againt the gradient being zero.
-  const double gradient_max_norm_0 = max(gradient_max_norm, kEpsilon);
-  double gradient_tolerance = options.gradient_tolerance * gradient_max_norm_0;
-
-  double mu = options.tau;
-  double nu = 2.0;
-  int iteration = 0;
-  double actual_cost_change = 0.0;
-  double step_norm = 0.0;
-  double relative_decrease = 0.0;
-
-  // Insane steps are steps which are not sane, i.e. there is some
-  // numerical kookiness going on with them. There are various reasons
-  // for this kookiness, some easier to diagnose then others. From the
-  // point of view of the non-linear solver, they are steps which
-  // cannot be used. We return with NUMERICAL_FAILURE after
-  // kMaxLinearSolverRetries consecutive insane steps.
-  bool step_is_sane = false;
-  int num_consecutive_insane_steps = 0;
-
-  // Whether the step resulted in a sufficient decrease in the
-  // objective function when compared to the decrease in the value of
-  // the lineariztion.
-  bool step_is_successful = false;
-
-  // Parse the iterations for which to dump the linear problem.
-  vector<int> iterations_to_dump = options.lsqp_iterations_to_dump;
-  sort(iterations_to_dump.begin(), iterations_to_dump.end());
-
-  IterationSummary iteration_summary;
-  iteration_summary.iteration = iteration;
-  iteration_summary.step_is_successful = false;
-  iteration_summary.cost = total_cost;
-  iteration_summary.cost_change = actual_cost_change;
-  iteration_summary.gradient_max_norm = gradient_max_norm;
-  iteration_summary.step_norm = step_norm;
-  iteration_summary.relative_decrease = relative_decrease;
-  iteration_summary.mu = mu;
-  iteration_summary.eta = options.eta;
-  iteration_summary.linear_solver_iterations = 0;
-  iteration_summary.linear_solver_time_sec = 0.0;
-  iteration_summary.iteration_time_sec = (time(NULL) - start_time);
-  if (options.logging_type >= PER_MINIMIZER_ITERATION) {
-    summary->iterations.push_back(iteration_summary);
-  }
-
-  // Check if the starting point is an optimum.
-  VLOG(2) << "Gradient max norm: " << gradient_max_norm
-          << " tolerance: " << gradient_tolerance
-          << " ratio: " << gradient_max_norm / gradient_max_norm_0
-          << " tolerance: " << options.gradient_tolerance;
-  if (gradient_max_norm <= gradient_tolerance) {
-    summary->termination_type = GRADIENT_TOLERANCE;
-    VLOG(1) << "Terminating on GRADIENT_TOLERANCE. "
-            << "Relative gradient max norm: "
-            << gradient_max_norm / gradient_max_norm_0
-            << " <= " << options.gradient_tolerance;
-    return;
-  }
-
-  // Call the various callbacks.
-  for (int i = 0; i < options.callbacks.size(); ++i) {
-    if (!RunCallback(options.callbacks[i], iteration_summary, summary)) {
-      return;
-    }
-  }
-
-  // We only need the LM diagonal if we are actually going to do at
-  // least one iteration of the optimization. So we wait to do it
-  // until now.
-  LevenbergMarquardtDiagonal(*jacobian, D.data());
-
-  while ((iteration < options.max_num_iterations) &&
-         (time(NULL) - start_time) <= options.max_solver_time_sec) {
-    time_t iteration_start_time = time(NULL);
-    step_is_sane = false;
-    step_is_successful = false;
-
-    IterationSummary iteration_summary;
-    // The while loop here is just to provide an easily breakable
-    // control structure. We are guaranteed to always exit this loop
-    // at the end of one iteration or before.
-    while (1) {
-      muD = (mu * D).array().sqrt();
-      LinearSolver::PerSolveOptions solve_options;
-      solve_options.D = muD.data();
-      solve_options.q_tolerance = options.eta;
-      // Disable r_tolerance checking. Since we only care about
-      // termination via the q_tolerance. As Nash and Sofer show,
-      // r_tolerance based termination is essentially useless in
-      // Truncated Newton methods.
-      solve_options.r_tolerance = -1.0;
-
-      // Invalidate the output array lm_step, so that we can detect if
-      // the linear solver generated numerical garbage.  This is known
-      // to happen for the DENSE_QR and then DENSE_SCHUR solver when
-      // the Jacobin is severly rank deficient and mu is too small.
-      InvalidateArray(num_effective_parameters, lm_step.data());
-      const time_t linear_solver_start_time = time(NULL);
-      LinearSolver::Summary linear_solver_summary =
-          linear_solver->Solve(jacobian.get(),
-                               f.data(),
-                               solve_options,
-                               lm_step.data());
-      iteration_summary.linear_solver_time_sec =
-          (time(NULL) - linear_solver_start_time);
-      iteration_summary.linear_solver_iterations =
-          linear_solver_summary.num_iterations;
-
-      if (binary_search(iterations_to_dump.begin(),
-                        iterations_to_dump.end(),
-                        iteration)) {
-        CHECK(DumpLinearLeastSquaresProblem(options.lsqp_dump_directory,
-                                            iteration,
-                                            options.lsqp_dump_format_type,
-                                            jacobian.get(),
-                                            muD.data(),
-                                            f.data(),
-                                            lm_step.data(),
-                                            options.num_eliminate_blocks))
-            << "Tried writing linear least squares problem: " 
-            << options.lsqp_dump_directory
-            << " but failed.";
-      }
-
-      // We ignore the case where the linear solver did not converge,
-      // since the partial solution computed by it still maybe of use,
-      // and there is no reason to ignore it, especially since we
-      // spent so much time computing it.
-      if ((linear_solver_summary.termination_type != TOLERANCE) &&
-          (linear_solver_summary.termination_type != MAX_ITERATIONS)) {
-        VLOG(1) << "Linear solver failure: retrying with a higher mu";
-        break;
-      }
-
-      if (!IsArrayValid(num_effective_parameters, lm_step.data())) {
-        LOG(WARNING) << "Linear solver failure. Failed to compute a finite "
-                     << "step. Terminating. Please report this to the Ceres "
-                     << "Solver developers.";
-        summary->termination_type = NUMERICAL_FAILURE;
-        return;
-      }
-
-      step_norm = (lm_step.array() * scale.array()).matrix().norm();
-
-      // Check step length based convergence. If the step length is
-      // too small, then we are done.
-      const double step_size_tolerance =  options.parameter_tolerance *
-          (x_norm + options.parameter_tolerance);
-
-      VLOG(2) << "Step size: " << step_norm
-              << " tolerance: " <<  step_size_tolerance
-              << " ratio: " << step_norm / step_size_tolerance
-              << " tolerance: " << options.parameter_tolerance;
-      if (step_norm <= options.parameter_tolerance *
-          (x_norm + options.parameter_tolerance)) {
-        summary->termination_type = PARAMETER_TOLERANCE;
-        VLOG(1) << "Terminating on PARAMETER_TOLERANCE."
-                << "Relative step size: " << step_norm / step_size_tolerance
-            << " <= " << options.parameter_tolerance;
-        return;
-      }
-
-      Vector delta =  -(lm_step.array() * scale.array()).matrix();
-      if (!evaluator->Plus(x.data(), delta.data(), x_new.data())) {
-        LOG(WARNING) << "Failed to compute Plus(x, delta, x_plus_delta). "
-                     << "Terminating.";
-        summary->termination_type = NUMERICAL_FAILURE;
-        return;
-      }
-
-      double cost_new = 0.0;
-      if (!evaluator->Evaluate(x_new.data(), &cost_new, NULL, NULL)) {
-        LOG(WARNING) << "Failed to compute the value of the objective "
-                     << "function. Terminating.";
-        summary->termination_type = NUMERICAL_FAILURE;
-        return;
-      }
-
-      f_model.setZero();
-      jacobian->RightMultiply(lm_step.data(), f_model.data());
-      const double model_cost_new =
-          (f.segment(0, num_residuals) - f_model).squaredNorm() / 2;
-
-      actual_cost_change = cost - cost_new;
-      double model_cost_change = model_cost - model_cost_new;
-
-      VLOG(2) << "[Model cost] current: " << model_cost
-              << " new : " << model_cost_new
-              << " change: " << model_cost_change;
-
-      VLOG(2) << "[Nonlinear cost] current: " << cost
-              << " new : " << cost_new
-              << " change: " << actual_cost_change
-              << " relative change: " << fabs(actual_cost_change) / cost
-              << " tolerance: " << options.function_tolerance;
-
-      // In exact arithmetic model_cost_change should never be
-      // negative. But due to numerical precision issues, we may end up
-      // with a small negative number. model_cost_change which are
-      // negative and large in absolute value are indicative of a
-      // numerical failure in the solver.
-      if (model_cost_change < -kEpsilon) {
-        VLOG(1) << "Model cost change is negative.\n"
-                << "Current : " << model_cost
-                << " new : " << model_cost_new
-                << " change: " << model_cost_change << "\n";
-        break;
-      }
-
-      // If we have reached this far, then we are willing to trust the
-      // numerical quality of the step.
-      step_is_sane = true;
-      num_consecutive_insane_steps = 0;
-
-      // Check function value based convergence.
-      if (fabs(actual_cost_change) < options.function_tolerance * cost) {
-        VLOG(1) << "Termination on FUNCTION_TOLERANCE."
-                << " Relative cost change: " << fabs(actual_cost_change) / cost
-                << " tolerance: " << options.function_tolerance;
-        summary->termination_type = FUNCTION_TOLERANCE;
-        return;
-      }
-
-      // Clamp model_cost_change at kEpsilon from below.
-      if (model_cost_change < kEpsilon) {
-        VLOG(1) << "Clamping model cost change " << model_cost_change
-                << " to " << kEpsilon;
-        model_cost_change = kEpsilon;
-      }
-
-      relative_decrease = actual_cost_change / model_cost_change;
-      VLOG(2) << "actual_cost_change / model_cost_change = "
-              << relative_decrease;
-
-      if (relative_decrease < options.min_relative_decrease) {
-        VLOG(2) << "Unsuccessful step.";
-        break;
-      }
-
-      VLOG(2) << "Successful step.";
-
-      ++summary->num_successful_steps;
-      x = x_new;
-      x_norm = x.norm();
-
-      if (!evaluator->Evaluate(x.data(), &cost, f.data(), jacobian.get())) {
-        LOG(WARNING) << "Failed to compute residuals and jacobian. "
-                     << "Terminating.";
-        summary->termination_type = NUMERICAL_FAILURE;
-        return;
-      }
-
-      if (options.jacobi_scaling) {
-        jacobian->ScaleColumns(scale.data());
-      }
-
-      model_cost = f.squaredNorm() / 2.0;
-      LevenbergMarquardtDiagonal(*jacobian, D.data());
-      scaled_gradient.setZero();
-      jacobian->LeftMultiply(f.data(), scaled_gradient.data());
-      gradient = scaled_gradient.array() / scale.array();
-      gradient_max_norm = gradient.lpNorm<Eigen::Infinity>();
-
-      // Check gradient based convergence.
-      VLOG(2) << "Gradient max norm: " << gradient_max_norm
-              << " tolerance: " << gradient_tolerance
-              << " ratio: " << gradient_max_norm / gradient_max_norm_0
-              << " tolerance: " << options.gradient_tolerance;
-      if (gradient_max_norm <= gradient_tolerance) {
-        summary->termination_type = GRADIENT_TOLERANCE;
-        VLOG(1) << "Terminating on GRADIENT_TOLERANCE. "
-                << "Relative gradient max norm: "
-                << gradient_max_norm / gradient_max_norm_0
-                << " <= " << options.gradient_tolerance
-                << " (tolerance).";
-        return;
-      }
-
-      mu = mu * max(1.0 / 3.0, 1 - pow(2 * relative_decrease - 1, 3));
-      mu = std::max(options.min_mu, mu);
-      nu = 2.0;
-      step_is_successful = true;
-      break;
-    }
-
-    if (!step_is_sane) {
-      ++num_consecutive_insane_steps;
-    }
-
-    if (num_consecutive_insane_steps == kMaxLinearSolverRetries) {
-      summary->termination_type = NUMERICAL_FAILURE;
-      VLOG(1) << "Too many consecutive retries; ending with numerical fail.";
-
-      if (!options.crash_and_dump_lsqp_on_failure) {
-        return;
-      }
-
-      // Dump debugging information to disk.
-      CHECK(options.lsqp_dump_format_type == TEXTFILE ||
-            options.lsqp_dump_format_type == PROTOBUF)
-          << "Dumping the linear least squares problem on crash "
-          << "requires Solver::Options::lsqp_dump_format_type to be "
-          << "PROTOBUF or TEXTFILE.";
-
-      if (DumpLinearLeastSquaresProblem(options.lsqp_dump_directory,
-                                        iteration,
-                                        options.lsqp_dump_format_type,
-                                        jacobian.get(),
-                                        muD.data(),
-                                        f.data(),
-                                        lm_step.data(),
-                                        options.num_eliminate_blocks)) {
-        LOG(FATAL) << "Linear least squares problem saved to: " 
-                   << options.lsqp_dump_directory
-                   << ". Please provide this to the Ceres developers for "
-                   << " debugging along with the v=2 log.";
-      } else {
-        LOG(FATAL) << "Tried writing linear least squares problem: " 
-                   << options.lsqp_dump_directory
-                   << " but failed.";
-      }
-    }
-
-    if (!step_is_successful) {
-      // Either the step did not lead to a decrease in cost or there
-      // was numerical failure. In either case we will scale mu up and
-      // retry. If it was a numerical failure, we hope that the
-      // stronger regularization will make the linear system better
-      // conditioned. If it was numerically sane, but there was no
-      // decrease in cost, then increasing mu reduces the size of the
-      // trust region and we look for a decrease closer to the
-      // linearization point.
-      ++summary->num_unsuccessful_steps;
-      mu = mu * nu;
-      nu = 2 * nu;
-    }
-
-    ++iteration;
-
-    total_cost = summary->fixed_cost + cost;
-
-    iteration_summary.iteration = iteration;
-    iteration_summary.step_is_successful = step_is_successful;
-    iteration_summary.cost = total_cost;
-    iteration_summary.cost_change = actual_cost_change;
-    iteration_summary.gradient_max_norm = gradient_max_norm;
-    iteration_summary.step_norm = step_norm;
-    iteration_summary.relative_decrease = relative_decrease;
-    iteration_summary.mu = mu;
-    iteration_summary.eta = options.eta;
-    iteration_summary.iteration_time_sec = (time(NULL) - iteration_start_time);
-
-    if (options.logging_type >= PER_MINIMIZER_ITERATION) {
-      summary->iterations.push_back(iteration_summary);
-    }
-
-    // Call the various callbacks.
-    for (int i = 0; i < options.callbacks.size(); ++i) {
-      if (!RunCallback(options.callbacks[i], iteration_summary, summary)) {
-        return;
-      }
-    }
-  }
-}
-
-}  // namespace internal
-}  // namespace ceres
diff --git a/internal/ceres/levenberg_marquardt_strategy.cc b/internal/ceres/levenberg_marquardt_strategy.cc
new file mode 100644
index 0000000..54ff783
--- /dev/null
+++ b/internal/ceres/levenberg_marquardt_strategy.cc
@@ -0,0 +1,103 @@
+#include "ceres/levenberg_marquardt_strategy.h"
+
+#include <cmath>
+#include "glog/logging.h"
+#include "ceres/array_utils.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/linear_solver.h"
+#include "ceres/sparse_matrix.h"
+#include "ceres/trust_region_strategy.h"
+#include "ceres/types.h"
+#include "Eigen/Core"
+
+namespace ceres {
+namespace internal {
+
+LevenbergMarquardtStrategy::LevenbergMarquardtStrategy(
+    const TrustRegionStrategy::Options& options)
+    : linear_solver_(options.linear_solver),
+      radius_(options.initial_radius),
+      max_radius_(options.max_radius),
+      min_diagonal_(options.lm_min_diagonal),
+      max_diagonal_(options.lm_max_diagonal),
+      decrease_factor_(2.0),
+      reuse_diagonal_(false) {
+  CHECK_NOTNULL(linear_solver_);
+  CHECK_GT(min_diagonal_, 0.0);
+  CHECK_LT(min_diagonal_, max_diagonal_);
+  CHECK_GT(max_radius_, 0.0);
+}
+
+LevenbergMarquardtStrategy::~LevenbergMarquardtStrategy() {
+}
+
+LinearSolver::Summary LevenbergMarquardtStrategy::ComputeStep(
+    const TrustRegionStrategy::PerSolveOptions& per_solve_options,
+    SparseMatrix* jacobian,
+    const double* residuals,
+    double* step) {
+  CHECK_NOTNULL(jacobian);
+  CHECK_NOTNULL(residuals);
+  CHECK_NOTNULL(step);
+
+  const int num_parameters = jacobian->num_cols();
+  if (!reuse_diagonal_) {
+    if (diagonal_.rows() != num_parameters) {
+      diagonal_.resize(num_parameters, 1);
+    }
+
+    jacobian->SquaredColumnNorm(diagonal_.data());
+    for (int i = 0; i < num_parameters; ++i) {
+      diagonal_[i] = min(max(diagonal_[i], min_diagonal_), max_diagonal_);
+    }
+  }
+
+  lm_diagonal_ = (diagonal_ / radius_).array().sqrt();
+
+  LinearSolver::PerSolveOptions solve_options;
+  solve_options.D = lm_diagonal_.data();
+  solve_options.q_tolerance = per_solve_options.eta;
+  // Disable r_tolerance checking. Since we only care about
+  // termination via the q_tolerance. As Nash and Sofer show,
+  // r_tolerance based termination is essentially useless in
+  // Truncated Newton methods.
+  solve_options.r_tolerance = -1.0;
+
+  // Invalidate the output array lm_step, so that we can detect if
+  // the linear solver generated numerical garbage.  This is known
+  // to happen for the DENSE_QR and then DENSE_SCHUR solver when
+  // the Jacobin is severly rank deficient and mu is too small.
+  InvalidateArray(num_parameters, step);
+  LinearSolver::Summary linear_solver_summary =
+      linear_solver_->Solve(jacobian, residuals, solve_options, step);
+  if (linear_solver_summary.termination_type == FAILURE ||
+      !IsArrayValid(num_parameters, step)) {
+    LOG(WARNING) << "Linear solver failure. Failed to compute a finite step.";
+    linear_solver_summary.termination_type = FAILURE;
+  }
+
+  reuse_diagonal_ = true;
+  return linear_solver_summary;
+}
+
+void LevenbergMarquardtStrategy::StepAccepted(double step_quality) {
+  CHECK_GT(step_quality, 0.0);
+  radius_ = radius_ / std::max(1.0 / 3.0,
+                               1.0 - pow(2.0 * step_quality - 1.0, 3));
+  radius_ = std::min(max_radius_, radius_);
+  decrease_factor_ = 2.0;
+  reuse_diagonal_ = false;
+}
+
+void LevenbergMarquardtStrategy::StepRejected(double step_quality) {
+  radius_ = radius_ / decrease_factor_;
+  decrease_factor_ *= 2.0;
+  reuse_diagonal_ = true;
+}
+
+double LevenbergMarquardtStrategy::Radius() const {
+  return radius_;
+}
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/internal/ceres/levenberg_marquardt_strategy.h b/internal/ceres/levenberg_marquardt_strategy.h
new file mode 100644
index 0000000..274e972
--- /dev/null
+++ b/internal/ceres/levenberg_marquardt_strategy.h
@@ -0,0 +1,78 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// 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_LEVENBERG_MARQUARDT_STRATEGY_H_
+#define CERES_INTERNAL_LEVENBERG_MARQUARDT_STRATEGY_H_
+
+#include "ceres/linear_solver.h"
+#include "ceres/trust_region_strategy.h"
+
+namespace ceres {
+namespace internal {
+
+// Levenberg-Marquardt step computation and trust region sizing
+// strategy based on on "Methods for Nonlinear Least Squares" by
+// K. Madsen, H.B. Nielsen and O. Tingleff. Available to download from
+//
+// http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
+class LevenbergMarquardtStrategy : public TrustRegionStrategy {
+public:
+  LevenbergMarquardtStrategy(const TrustRegionStrategy::Options& options);
+  virtual ~LevenbergMarquardtStrategy();
+
+  // TrustRegionStrategy interface
+  virtual LinearSolver::Summary ComputeStep(
+      const TrustRegionStrategy::PerSolveOptions& per_solve_options,
+      SparseMatrix* jacobian,
+      const double* residuals,
+      double* step);
+  virtual void StepAccepted(double step_quality);
+  virtual void StepRejected(double step_quality);
+  virtual double Radius() const;
+
+ private:
+  LinearSolver* linear_solver_;
+  double radius_;
+  double max_radius_;
+  const double min_diagonal_;
+  const double max_diagonal_;
+  double decrease_factor_;
+  bool reuse_diagonal_;
+  Vector diagonal_;   // diagonal_ =  diag(J'J)
+  // Scaled copy of diagonal_. Stored here as optimization to prevent
+  // allocations in every iteration and reuse when a step fails and
+  // ComputeStep is called again.
+  Vector lm_diagonal_;  // lm_diagonal_ = diagonal_ / radius_;
+};
+
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // CERES_INTERNAL_LEVENBERG_MARQUARDT_STRATEGY_H_
diff --git a/internal/ceres/levenberg_marquardt_strategy_test.cc b/internal/ceres/levenberg_marquardt_strategy_test.cc
new file mode 100644
index 0000000..370d0a5
--- /dev/null
+++ b/internal/ceres/levenberg_marquardt_strategy_test.cc
@@ -0,0 +1,152 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// 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 <glog/logging.h>
+#include "gmock/gmock.h"
+#include "gtest/gtest.h"
+#include "ceres/mock_log.h"
+#include "ceres/trust_region_strategy.h"
+#include "ceres/levenberg_marquardt_strategy.h"
+#include "ceres/internal/eigen.h"
+
+using testing::AllOf;
+using testing::AnyNumber;
+using testing::HasSubstr;
+using testing::ScopedMockLog;
+using testing::_;
+
+namespace ceres {
+namespace internal {
+
+const double kTolerance = 1e-16;
+
+// Linear solver that takes as input a vector and checks that the
+// caller passes the same vector as LinearSolver::PerSolveOptions.D.
+class RegularizationCheckingLinearSolver : public DenseSparseMatrixSolver {
+ public:
+  RegularizationCheckingLinearSolver(const int num_cols, const double* diagonal)
+      : num_cols_(num_cols),
+        diagonal_(diagonal) {
+  }
+
+  virtual ~RegularizationCheckingLinearSolver(){}
+
+ private:
+  virtual LinearSolver::Summary SolveImpl(
+      DenseSparseMatrix* A,
+      const double* b,
+      const LinearSolver::PerSolveOptions& per_solve_options,
+      double* x) {
+    CHECK_NOTNULL(per_solve_options.D);
+    for (int i = 0; i < num_cols_; ++i) {
+      EXPECT_NEAR(per_solve_options.D[i], diagonal_[i], kTolerance)
+          << i << " " << per_solve_options.D[i] << " " << diagonal_[i];
+    }
+  }
+  const int num_cols_;
+  const double* diagonal_;
+};
+
+TEST(LevenbergMarquardtStrategy, AcceptRejectStepRadiusScaling) {
+  TrustRegionStrategy::Options options;
+  options.initial_radius = 2.0;
+  options.max_radius = 20.0;
+  options.lm_min_diagonal = 1e-8;
+  options.lm_max_diagonal = 1e8;
+
+  // We need a non-null pointer here, so anything should do.
+  options.linear_solver = new RegularizationCheckingLinearSolver(0, NULL);
+
+  LevenbergMarquardtStrategy lms(options);
+  EXPECT_EQ(lms.Radius(), options.initial_radius);
+  lms.StepRejected(0.0);
+  EXPECT_EQ(lms.Radius(), 1.0);
+  lms.StepRejected(-1.0);
+  EXPECT_EQ(lms.Radius(), 0.25);
+  lms.StepAccepted(1.0);
+  EXPECT_EQ(lms.Radius(), 0.25 * 3.0);
+  lms.StepAccepted(1.0);
+  EXPECT_EQ(lms.Radius(), 0.25 * 3.0 * 3.0);
+  lms.StepAccepted(0.25);
+  EXPECT_EQ(lms.Radius(), 0.25 * 3.0 * 3.0 / 1.125);
+  lms.StepAccepted(1.0);
+  EXPECT_EQ(lms.Radius(), 0.25 * 3.0 * 3.0 / 1.125 * 3.0);
+  lms.StepAccepted(1.0);
+  EXPECT_EQ(lms.Radius(), 0.25 * 3.0 * 3.0 / 1.125 * 3.0 * 3.0);
+  lms.StepAccepted(1.0);
+  EXPECT_EQ(lms.Radius(), options.max_radius);
+}
+
+TEST(LevenbergMarquardtStrategy, CorrectDiagonalToLinearSolver) {
+  Matrix jacobian(2,3);
+  jacobian.setZero();
+  jacobian(0,0) = 0.0;
+  jacobian(0,1) = 1.0;
+  jacobian(1,1) = 1.0;
+  jacobian(0,2) = 100.0;
+
+  double residual = 1.0;
+  double x[2];
+  DenseSparseMatrix dsm(jacobian);
+
+  TrustRegionStrategy::Options options;
+  options.initial_radius = 2.0;
+  options.max_radius = 20.0;
+  options.lm_min_diagonal = 1e-2;
+  options.lm_max_diagonal = 1e2;
+
+  double diagonal[3];
+  diagonal[0] = options.lm_min_diagonal;
+  diagonal[1] = 2.0;
+  diagonal[2] = options.lm_max_diagonal;
+  for (int i = 0; i < 3; ++i) {
+    diagonal[i] = sqrt(diagonal[i] / options.initial_radius);
+  }
+
+  RegularizationCheckingLinearSolver linear_solver(3, diagonal);
+  options.linear_solver = &linear_solver;
+
+  LevenbergMarquardtStrategy lms(options);
+  TrustRegionStrategy::PerSolveOptions pso;
+
+  {
+    ScopedMockLog log;
+    EXPECT_CALL(log, Log(_, _, _)).Times(AnyNumber());
+    EXPECT_CALL(log, Log(WARNING, _,
+                         HasSubstr("Failed to compute a finite step.")));
+
+    LinearSolver::Summary summary = lms.ComputeStep(pso, &dsm, &residual, x);
+    EXPECT_EQ(summary.termination_type, FAILURE);
+  }
+}
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/internal/ceres/minimizer.h b/internal/ceres/minimizer.h
index 19b0344..eeda298 100644
--- a/internal/ceres/minimizer.h
+++ b/internal/ceres/minimizer.h
@@ -40,6 +40,8 @@
 
 class Evaluator;
 class LinearSolver;
+class SparseMatrix;
+class TrustRegionStrategy;
 
 // Interface for non-linear least squares solvers.
 class Minimizer {
@@ -48,55 +50,87 @@
   // see solver.h for detailed information about the meaning and
   // default values of each of these parameters.
   struct Options {
+    Options() {
+      Init(Solver::Options());
+    }
+
     explicit Options(const Solver::Options& options) {
+      Init(options);
+    }
+
+    void Init(const Solver::Options& options) {
       max_num_iterations = options.max_num_iterations;
       max_solver_time_sec = options.max_solver_time_sec;
+      max_step_solver_retries = 5;
       gradient_tolerance = options.gradient_tolerance;
       parameter_tolerance = options.parameter_tolerance;
       function_tolerance = options.function_tolerance;
       min_relative_decrease = options.min_relative_decrease;
       eta = options.eta;
-      tau = options.tau;
-      min_mu = options.min_mu;
       jacobi_scaling = options.jacobi_scaling;
-      crash_and_dump_lsqp_on_failure = options.crash_and_dump_lsqp_on_failure;
       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;
       num_eliminate_blocks = options.num_eliminate_blocks;
-      logging_type = options.logging_type;
+      max_num_consecutive_invalid_steps =
+          options.max_num_consecutive_invalid_steps;
+      min_trust_region_radius = options.min_trust_region_radius;
+      evaluator = NULL;
+      trust_region_strategy = NULL;
+      jacobian = NULL;
     }
 
     int max_num_iterations;
     int max_solver_time_sec;
+
+    // Number of times the linear solver should be retried in case of
+    // numerical failure. The retries are done by exponentially scaling up
+    // mu at each retry. This leads to stronger and stronger
+    // regularization making the linear least squares problem better
+    // conditioned at each retry.
+    int max_step_solver_retries;
     double gradient_tolerance;
     double parameter_tolerance;
     double function_tolerance;
     double min_relative_decrease;
     double eta;
-    double tau;
-    double min_mu;
     bool jacobi_scaling;
-    bool crash_and_dump_lsqp_on_failure;
     vector<int> lsqp_iterations_to_dump;
     DumpFormatType lsqp_dump_format_type;
     string lsqp_dump_directory;
     int num_eliminate_blocks;
-    LoggingType logging_type;
+    int max_num_consecutive_invalid_steps;
+    int min_trust_region_radius;
 
     // List of callbacks that are executed by the Minimizer at the end
     // of each iteration.
     //
-    // Client owns these pointers.
+    // The Options struct does not own these pointers.
     vector<IterationCallback*> callbacks;
+
+    // Object responsible for evaluating the cost, residuals and
+    // Jacobian matrix. The Options struct does not own this pointer.
+    Evaluator* evaluator;
+
+    // Object responsible for actually computing the trust region
+    // step, and sizing the trust region radius. The Options struct
+    // does not own this pointer.
+    TrustRegionStrategy* trust_region_strategy;
+
+    // Object holding the Jacobian matrix. It is assumed that the
+    // sparsity structure of the matrix has already been initialized
+    // and will remain constant for the life time of the
+    // optimization. The Options struct does not own this pointer.
+    SparseMatrix* jacobian;
   };
 
   virtual ~Minimizer() {}
+
+  // Note: The minimizer is expected to update the state of the
+  // parameters array every iteration. This is required for the
+  // StateUpdatingCallback to work.
   virtual void Minimize(const Options& options,
-                        Evaluator* evaluator,
-                        LinearSolver* linear_solver,
-                        const double* initial_parameters,
-                        double* final_parameters,
+                        double* parameters,
                         Solver::Summary* summary) = 0;
 };
 
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc
index d642a54..f5f5d91 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -32,7 +32,6 @@
 #include "ceres/solver.h"
 
 #include <vector>
-#include "ceres/levenberg_marquardt.h"
 #include "ceres/program.h"
 #include "ceres/solver_impl.h"
 #include "ceres/stringprintf.h"
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc
index 6f8fa17..4eb3d1f 100644
--- a/internal/ceres/solver_impl.cc
+++ b/internal/ceres/solver_impl.cc
@@ -34,18 +34,19 @@
 #include <numeric>
 #include "ceres/evaluator.h"
 #include "ceres/gradient_checking_cost_function.h"
-#include "ceres/levenberg_marquardt.h"
+#include "ceres/iteration_callback.h"
+#include "ceres/levenberg_marquardt_strategy.h"
 #include "ceres/linear_solver.h"
 #include "ceres/map_util.h"
 #include "ceres/minimizer.h"
 #include "ceres/parameter_block.h"
+#include "ceres/problem.h"
 #include "ceres/problem_impl.h"
 #include "ceres/program.h"
 #include "ceres/residual_block.h"
 #include "ceres/schur_ordering.h"
 #include "ceres/stringprintf.h"
-#include "ceres/iteration_callback.h"
-#include "ceres/problem.h"
+#include "ceres/trust_region_minimizer.h"
 
 namespace ceres {
 namespace internal {
@@ -104,7 +105,7 @@
                                  summary.gradient_max_norm,
                                  summary.step_norm,
                                  summary.relative_decrease,
-                                 summary.mu,
+                                 summary.trust_region_radius,
                                  summary.linear_solver_iterations);
     if (log_to_stdout_) {
       cout << output << endl;
@@ -124,32 +125,40 @@
                           Program* program,
                           Evaluator* evaluator,
                           LinearSolver* linear_solver,
-                          double* initial_parameters,
-                          double* final_parameters,
+                          double* parameters,
                           Solver::Summary* summary) {
   Minimizer::Options minimizer_options(options);
-
   LoggingCallback logging_callback(options.minimizer_progress_to_stdout);
   if (options.logging_type != SILENT) {
     minimizer_options.callbacks.push_back(&logging_callback);
   }
 
-  StateUpdatingCallback updating_callback(program, final_parameters);
+  StateUpdatingCallback updating_callback(program, parameters);
   if (options.update_state_every_iteration) {
     minimizer_options.callbacks.push_back(&updating_callback);
   }
 
-  LevenbergMarquardt levenberg_marquardt;
+  minimizer_options.evaluator = evaluator;
+  scoped_ptr<SparseMatrix> jacobian(evaluator->CreateJacobian());
+  minimizer_options.jacobian = jacobian.get();
 
-  time_t start_minimizer_time_seconds = time(NULL);
-  levenberg_marquardt.Minimize(minimizer_options,
-                               evaluator,
-                               linear_solver,
-                               initial_parameters,
-                               final_parameters,
-                               summary);
-  summary->minimizer_time_in_seconds =
-      time(NULL) - start_minimizer_time_seconds;
+  TrustRegionStrategy::Options trust_region_strategy_options;
+  trust_region_strategy_options.linear_solver = linear_solver;
+  trust_region_strategy_options.initial_radius =
+      options.initial_trust_region_radius;
+  trust_region_strategy_options.max_radius = options.max_trust_region_radius;
+  trust_region_strategy_options.lm_min_diagonal = options.lm_min_diagonal;
+  trust_region_strategy_options.lm_max_diagonal = options.lm_max_diagonal;
+  trust_region_strategy_options.trust_region_strategy_type =
+      options.trust_region_strategy_type;
+  scoped_ptr<TrustRegionStrategy> strategy(
+      TrustRegionStrategy::Create(trust_region_strategy_options));
+  minimizer_options.trust_region_strategy = strategy.get();
+
+  TrustRegionMinimizer minimizer;
+  time_t minimizer_start_time = time(NULL);
+  minimizer.Minimize(minimizer_options, parameters, summary);
+  summary->minimizer_time_in_seconds = time(NULL) - minimizer_start_time;
 }
 
 void SolverImpl::Solve(const Solver::Options& original_options,
@@ -259,19 +268,17 @@
   }
 
   // The optimizer works on contiguous parameter vectors; allocate some.
-  Vector initial_parameters(reduced_program->NumParameters());
-  Vector optimized_parameters(reduced_program->NumParameters());
+  Vector parameters(reduced_program->NumParameters());
 
   // Collect the discontiguous parameters into a contiguous state vector.
-  reduced_program->ParameterBlocksToStateVector(&initial_parameters[0]);
+  reduced_program->ParameterBlocksToStateVector(parameters.data());
 
   // Run the optimization.
   Minimize(options,
            reduced_program.get(),
            evaluator.get(),
            linear_solver.get(),
-           initial_parameters.data(),
-           optimized_parameters.data(),
+           parameters.data(),
            summary);
 
   // If the user aborted mid-optimization or the optimization
@@ -283,7 +290,7 @@
   }
 
   // Push the contiguous optimized parameters back to the user's parameters.
-  reduced_program->StateVectorToParameterBlocks(&optimized_parameters[0]);
+  reduced_program->StateVectorToParameterBlocks(parameters.data());
   reduced_program->CopyParameterBlockStateToUserState();
 
   // Return the final cost and residuals for the original problem.
diff --git a/internal/ceres/solver_impl.h b/internal/ceres/solver_impl.h
index 957ebcc..7dee03c 100644
--- a/internal/ceres/solver_impl.h
+++ b/internal/ceres/solver_impl.h
@@ -92,8 +92,7 @@
                        Program* program,
                        Evaluator* evaluator,
                        LinearSolver* linear_solver,
-                       double* initial_parameters,
-                       double* final_parameters,
+                       double* parameters,
                        Solver::Summary* summary);
 
   // Remove the fixed or unused parameter blocks and residuals
diff --git a/internal/ceres/trust_region_minimizer.cc b/internal/ceres/trust_region_minimizer.cc
new file mode 100644
index 0000000..bb6c4a1
--- /dev/null
+++ b/internal/ceres/trust_region_minimizer.cc
@@ -0,0 +1,445 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// 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 "ceres/trust_region_minimizer.h"
+
+#include <algorithm>
+#include <cstdlib>
+#include <cmath>
+#include <cstring>
+#include <string>
+#include <vector>
+#include <glog/logging.h>
+
+#include "Eigen/Core"
+#include "ceres/array_utils.h"
+#include "ceres/evaluator.h"
+#include "ceres/linear_least_squares_problems.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/sparse_matrix.h"
+#include "ceres/trust_region_strategy.h"
+#include "ceres/types.h"
+
+namespace ceres {
+namespace internal {
+namespace {
+// Small constant for various floating point issues.
+const double kEpsilon = 1e-12;
+}  // namespace
+
+// Execute the list of IterationCallbacks sequentially. If any one of
+// the callbacks does not return SOLVER_CONTINUE, then stop and return
+// its status.
+CallbackReturnType TrustRegionMinimizer::RunCallbacks(
+    const Minimizer::Options& options_,
+    const IterationSummary& iteration_summary) {
+  for (int i = 0; i < options_.callbacks.size(); ++i) {
+    const CallbackReturnType status =
+        (*options_.callbacks[i])(iteration_summary);
+    if (status != SOLVER_CONTINUE) {
+      return status;
+    }
+  }
+  return SOLVER_CONTINUE;
+}
+
+// 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 / (kEpsilon + sqrt(scale[i]));
+  }
+}
+
+void TrustRegionMinimizer::Init(const Minimizer::Options& options) {
+  options_ = options;
+  sort(options_.lsqp_iterations_to_dump.begin(),
+       options_.lsqp_iterations_to_dump.end());
+}
+
+bool TrustRegionMinimizer::MaybeDumpLinearLeastSquaresProblem(
+    const int iteration,
+    const SparseMatrix* jacobian,
+    const double* residuals,
+    const double* step) const  {
+  // TODO(sameeragarwal): Since the use of trust_region_radius has
+  // moved inside TrustRegionStrategy, its not clear how we dump the
+  // regularization vector/matrix anymore.
+  //
+  // Doing this right requires either an API change to the
+  // TrustRegionStrategy and/or how LinearLeastSquares problems are
+  // stored on disk.
+  //
+  // For now, we will just not dump the regularizer.
+  return (!binary_search(options_.lsqp_iterations_to_dump.begin(),
+                         options_.lsqp_iterations_to_dump.end(),
+                         iteration) ||
+          DumpLinearLeastSquaresProblem(options_.lsqp_dump_directory,
+                                        iteration,
+                                        options_.lsqp_dump_format_type,
+                                        jacobian,
+                                        NULL,
+                                        residuals,
+                                        step,
+                                        options_.num_eliminate_blocks));
+}
+
+void TrustRegionMinimizer::Minimize(const Minimizer::Options& options,
+                                    double* parameters,
+                                    Solver::Summary* summary) {
+  time_t start_time = time(NULL);
+  time_t iteration_start_time =  start_time;
+  Init(options);
+
+  summary->termination_type = NO_CONVERGENCE;
+  summary->num_successful_steps = 0;
+  summary->num_unsuccessful_steps = 0;
+
+  Evaluator* evaluator = CHECK_NOTNULL(options_.evaluator);
+  SparseMatrix* jacobian = CHECK_NOTNULL(options_.jacobian);
+  TrustRegionStrategy* strategy = CHECK_NOTNULL(options_.trust_region_strategy);
+
+  const int num_parameters = evaluator->NumParameters();
+  const int num_effective_parameters = evaluator->NumEffectiveParameters();
+  const int num_residuals = evaluator->NumResiduals();
+
+  VectorRef x(parameters, num_parameters);
+  double x_norm = x.norm();
+
+  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);
+
+  IterationSummary iteration_summary;
+  iteration_summary.iteration = 0;
+  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;
+  iteration_summary.step_norm = 0.0;
+  iteration_summary.relative_decrease = 0.0;
+  iteration_summary.trust_region_radius = strategy->Radius();
+  // TODO(sameeragarwal): Rename eta to linear_solver_accuracy or
+  // something similar across the board.
+  iteration_summary.eta = options_.eta;
+  iteration_summary.linear_solver_iterations = 0;
+  iteration_summary.step_solver_time_sec = 0;
+
+  // Do initial cost and Jacobian evaluation.
+  double cost = 0.0;
+  if (!evaluator->Evaluate(x.data(), &cost, residuals.data(), jacobian)) {
+    LOG(WARNING) << "Terminating: Residual and Jacobian evaluation failed.";
+    summary->termination_type = NUMERICAL_FAILURE;
+    return;
+  }
+
+  // Compute the fixed part of the cost.
+  //
+  // This is a poor way to do this computation. Even if fixed_cost is
+  // zero, because we are subtracting two possibly large numbers, we
+  // are depending on exact cancellation to give us a zero here. But
+  // initial_cost and cost have been computed by two different
+  // evaluators. One which runs on the whole problem (in
+  // solver_impl.cc) in single threaded mode and another which runs
+  // here on the reduced problem, so fixed_cost can (and does) contain
+  // some numerical garbage with a relative magnitude of 1e-14.
+  //
+  // The right way to do this, would be to compute the fixed cost on
+  // just the set of residual blocks which are held constant and were
+  // removed from the original problem when the reduced problem was
+  // constructed.
+  summary->fixed_cost = summary->initial_cost - cost;
+
+  gradient.setZero();
+  jacobian->LeftMultiply(residuals.data(), gradient.data());
+  iteration_summary.gradient_max_norm = gradient.lpNorm<Eigen::Infinity>();
+
+  if (options_.jacobi_scaling) {
+    EstimateScale(*jacobian, scale.data());
+    jacobian->ScaleColumns(scale.data());
+  } else {
+    scale.setOnes();
+  }
+
+  const double absolute_gradient_tolerance =
+      options_.gradient_tolerance *
+      max(iteration_summary.gradient_max_norm, kEpsilon);
+  if (iteration_summary.gradient_max_norm <= absolute_gradient_tolerance) {
+    summary->termination_type = GRADIENT_TOLERANCE;
+    VLOG(1) << "Terminating: Gradient tolerance reached."
+            << "Relative gradient max norm: "
+            << iteration_summary.gradient_max_norm /
+        summary->iterations[0].gradient_max_norm
+            << " <= " << options_.gradient_tolerance;
+    return;
+  }
+
+  // Call the various callbacks.
+  iteration_summary.iteration_time_sec = time(NULL) - iteration_start_time;
+  summary->iterations.push_back(iteration_summary);
+  switch (RunCallbacks(options_, iteration_summary)) {
+    case SOLVER_TERMINATE_SUCCESSFULLY:
+      summary->termination_type = USER_SUCCESS;
+      VLOG(1) << "Terminating: User callback returned USER_SUCCESS.";
+      return;
+    case SOLVER_ABORT:
+      summary->termination_type = USER_ABORT;
+      VLOG(1) << "Terminating: User callback returned  USER_ABORT.";
+      return;
+    case SOLVER_CONTINUE:
+      break;
+    default:
+      LOG(FATAL) << "Unknown type of user callback status";
+  }
+
+  int num_consecutive_invalid_steps = 0;
+  while (true) {
+    iteration_start_time = time(NULL);
+    if (iteration_summary.iteration >= options_.max_num_iterations) {
+      summary->termination_type = NO_CONVERGENCE;
+      VLOG(1) << "Terminating: Maximum number of iterations reached.";
+      break;
+    }
+
+    if ((start_time - iteration_start_time) >= options_.max_solver_time_sec) {
+      summary->termination_type = NO_CONVERGENCE;
+      VLOG(1) << "Terminating: Maximum solver time reached.";
+      break;
+    }
+
+    iteration_summary = IterationSummary();
+    iteration_summary = summary->iterations.back();
+    iteration_summary.iteration = summary->iterations.back().iteration + 1;
+    iteration_summary.step_is_valid = false;
+    iteration_summary.step_is_successful = false;
+
+    const time_t strategy_start_time = time(NULL);
+    TrustRegionStrategy::PerSolveOptions per_solve_options;
+    per_solve_options.eta = options_.eta;
+    LinearSolver::Summary strategy_summary =
+        strategy->ComputeStep(per_solve_options,
+                              jacobian,
+                              residuals.data(),
+                              trust_region_step.data());
+
+    iteration_summary.step_solver_time_sec = time(NULL) - strategy_start_time;
+    iteration_summary.linear_solver_iterations =
+        strategy_summary.num_iterations;
+
+    if (!MaybeDumpLinearLeastSquaresProblem(iteration_summary.iteration,
+                                            jacobian,
+                                            residuals.data(),
+                                            trust_region_step.data())) {
+      LOG(FATAL) << "Tried writing linear least squares problem: "
+                 << options.lsqp_dump_directory << "but failed.";
+    }
+
+    double new_model_cost = 0.0;
+    if (strategy_summary.termination_type != FAILURE) {
+      // new_model_cost = 1/2 |J * step - f|^2
+      model_residuals = -residuals;
+      jacobian->RightMultiply(trust_region_step.data(), model_residuals.data());
+      new_model_cost = model_residuals.squaredNorm() / 2.0;
+
+      // In exact arithmetic, this would never be the case. But poorly
+      // conditioned matrices can give rise to situations where the
+      // new_model_cost can actually be larger than half the squared
+      // norm of the residual vector. We allow for small tolerance
+      // around cost and beyond that declare the step to be invalid.
+      if (cost < (new_model_cost - kEpsilon)) {
+        VLOG(1) << "Invalid step: current_cost: " << cost
+                << " new_model_cost " << new_model_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
+      // NUMERICAL_FAILURE if this limit is exceeded.
+      if (++num_consecutive_invalid_steps >=
+          options_.max_num_consecutive_invalid_steps) {
+        summary->termination_type = NUMERICAL_FAILURE;
+        LOG(WARNING) << "Terminating. Number of successive invalid steps more "
+                     << "than "
+                     << "Solver::Options::max_num_consecutive_invalid_steps: "
+                     << options_.max_num_consecutive_invalid_steps;
+        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;
+      iteration_summary.cost_change = 0.0;
+      iteration_summary.gradient_max_norm =
+          summary->iterations.back().gradient_max_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;
+
+      // We allow some slop around 0, and clamp the model_cost_change
+      // at kEpsilon from below.
+      //
+      // There is probably a better way to do this, as it is going to
+      // create problems for problems where the objective function is
+      // kEpsilon close to zero.
+      const double model_cost_change = max(kEpsilon, cost - new_model_cost);
+
+      // Undo the Jacobian column scaling.
+      delta =  -(trust_region_step.array() * scale.array()).matrix();
+      iteration_summary.step_norm = 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) {
+        VLOG(1) << "Terminating. Parameter tolerance reached. "
+                << "relative step_norm: "
+                << iteration_summary.step_norm /
+            (x_norm + options_.parameter_tolerance)
+                << " <= " << options_.parameter_tolerance;
+        summary->termination_type = PARAMETER_TOLERANCE;
+        return;
+      }
+
+      if (!evaluator->Plus(x.data(), delta.data(), x_plus_delta.data())) {
+        summary->termination_type = NUMERICAL_FAILURE;
+        LOG(WARNING) << "Terminating. Failed to compute "
+                     << "Plus(x, delta, x_plus_delta).";
+        return;
+      }
+
+      // Try this step.
+      double new_cost;
+      if (!evaluator->Evaluate(x_plus_delta.data(), &new_cost, NULL, NULL)) {
+        summary->termination_type = NUMERICAL_FAILURE;
+        LOG(WARNING) << "Terminating: Cost evaluation failed.";
+        return;
+      }
+
+      VLOG(2) << "old cost: " << cost << "new cost: " << new_cost;
+      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) {
+        VLOG(1) << "Terminating. Function tolerance reached. "
+                << "|cost_change|/cost: " << fabs(iteration_summary.cost_change) / cost
+                << " <= " << options_.function_tolerance;
+        summary->termination_type = FUNCTION_TOLERANCE;
+        return;
+      }
+
+      iteration_summary.relative_decrease =
+          iteration_summary.cost_change / model_cost_change;
+      iteration_summary.step_is_successful =
+          iteration_summary.relative_decrease > options_.min_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(), jacobian)) {
+        summary->termination_type = NUMERICAL_FAILURE;
+        LOG(WARNING) << "Terminating: Residual and Jacobian evaluation failed.";
+        return;
+      }
+
+      gradient.setZero();
+      jacobian->LeftMultiply(residuals.data(), gradient.data());
+      iteration_summary.gradient_max_norm = gradient.lpNorm<Eigen::Infinity>();
+
+      if (iteration_summary.gradient_max_norm <= absolute_gradient_tolerance) {
+        summary->termination_type = GRADIENT_TOLERANCE;
+        VLOG(1) << "Terminating: Gradient tolerance reached."
+                << "Relative gradient max norm: "
+                << iteration_summary.gradient_max_norm /
+            summary->iterations[0].gradient_max_norm
+                << " <= " << options_.gradient_tolerance;
+        return;
+      }
+
+      if (options_.jacobi_scaling) {
+        jacobian->ScaleColumns(scale.data());
+      }
+    } else {
+      ++summary->num_unsuccessful_steps;
+      strategy->StepRejected(iteration_summary.relative_decrease);
+    }
+
+    iteration_summary.cost = cost + summary->fixed_cost;
+    iteration_summary.trust_region_radius = strategy->Radius();
+    if (iteration_summary.trust_region_radius <
+        options_.min_trust_region_radius) {
+      summary->termination_type = PARAMETER_TOLERANCE;
+      VLOG(1) << "Termination. Minimum trust region radius reached.";
+      return;
+    }
+
+    iteration_summary.iteration_time_sec = time(NULL) - iteration_start_time;
+    summary->iterations.push_back(iteration_summary);
+    switch (RunCallbacks(options_, iteration_summary)) {
+      case SOLVER_TERMINATE_SUCCESSFULLY:
+        summary->termination_type = USER_SUCCESS;
+        VLOG(1) << "Terminating: User callback returned USER_SUCCESS.";
+        return;
+      case SOLVER_ABORT:
+        summary->termination_type = USER_ABORT;
+        VLOG(1) << "Terminating: User callback returned  USER_ABORT.";
+        return;
+      case SOLVER_CONTINUE:
+        break;
+      default:
+        LOG(FATAL) << "Unknown type of user callback status";
+    }
+  }
+}
+
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/internal/ceres/levenberg_marquardt.h b/internal/ceres/trust_region_minimizer.h
similarity index 62%
rename from internal/ceres/levenberg_marquardt.h
rename to internal/ceres/trust_region_minimizer.h
index d00bb90..4337b18 100644
--- a/internal/ceres/levenberg_marquardt.h
+++ b/internal/ceres/trust_region_minimizer.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
+// Copyright 2012 Google Inc. All rights reserved.
 // http://code.google.com/p/ceres-solver/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -27,39 +27,42 @@
 // POSSIBILITY OF SUCH DAMAGE.
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
-//
-// Implmentation of Levenberg Marquardt algorithm based on "Methods for
-// Nonlinear Least Squares" by K. Madsen, H.B. Nielsen and
-// O. Tingleff. Available to download from
-//
-// http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
-//
 
-#ifndef CERES_INTERNAL_LEVENBERG_MARQUARDT_H_
-#define CERES_INTERNAL_LEVENBERG_MARQUARDT_H_
+#ifndef CERES_INTERNAL_TRUST_REGION_MINIMIZER_H_
+#define CERES_INTERNAL_TRUST_REGION_MINIMIZER_H_
 
 #include "ceres/minimizer.h"
 #include "ceres/solver.h"
+#include "ceres/types.h"
 
 namespace ceres {
 namespace internal {
 
-class Evaluator;
-class LinearSolver;
-
-class LevenbergMarquardt : public Minimizer {
+// Generic trust region minimization algorithm. The heavy lifting is
+// done by a TrustRegionStrategy object passed in as one of the
+// arguments to the Minimize method.
+//
+// For example usage, see SolverImpl::Minimize.
+class TrustRegionMinimizer : public Minimizer {
  public:
-  virtual ~LevenbergMarquardt();
-
+  ~TrustRegionMinimizer() {}
   virtual void Minimize(const Minimizer::Options& options,
-                        Evaluator* evaluator,
-                        LinearSolver* linear_solver,
-                        const double* initial_parameters,
-                        double* final_parameters,
+                        double* parameters,
                         Solver::Summary* summary);
+
+ private:
+  void Init(const Minimizer::Options& options);
+  void EstimateScale(const SparseMatrix& jacobian, double* scale) const;
+  CallbackReturnType RunCallbacks(const Minimizer::Options& options,
+                                  const IterationSummary& iteration_summary);
+  bool MaybeDumpLinearLeastSquaresProblem( const int iteration,
+                                           const SparseMatrix* jacobian,
+                                           const double* residuals,
+                                           const double* step) const;
+
+  Minimizer::Options options_;
 };
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif  // CERES_INTERNAL_LEVENBERG_MARQUARDT_H_
+#endif  // CERES_INTERNAL_TRUST_REGION_MINIMIZER_H_
diff --git a/internal/ceres/levenberg_marquardt_test.cc b/internal/ceres/trust_region_minimizer_test.cc
similarity index 72%
rename from internal/ceres/levenberg_marquardt_test.cc
rename to internal/ceres/trust_region_minimizer_test.cc
index 020abfa..b4ef601 100644
--- a/internal/ceres/levenberg_marquardt_test.cc
+++ b/internal/ceres/trust_region_minimizer_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
+// Copyright 2012 Google Inc. All rights reserved.
 // http://code.google.com/p/ceres-solver/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -27,19 +27,21 @@
 // POSSIBILITY OF SUCH DAMAGE.
 //
 // Author: keir@google.com (Keir Mierle)
+//         sameeragarwal@google.com (Sameer Agarwal)
 //
-// This tests the Levenberg-Marquardt loop using a direct Evaluator
-// implementation, rather than having a test that goes through all the Program
-// and Problem machinery.
+// This tests the TrustRegionMinimizer loop using a direct Evaluator
+// implementation, rather than having a test that goes through all the
+// Program and Problem machinery.
 
 #include <cmath>
 #include "ceres/dense_qr_solver.h"
 #include "ceres/dense_sparse_matrix.h"
 #include "ceres/evaluator.h"
-#include "ceres/levenberg_marquardt.h"
+#include "ceres/internal/port.h"
+#include "ceres/levenberg_marquardt_strategy.h"
 #include "ceres/linear_solver.h"
 #include "ceres/minimizer.h"
-#include "ceres/internal/port.h"
+#include "ceres/trust_region_minimizer.h"
 #include "gtest/gtest.h"
 
 namespace ceres {
@@ -181,67 +183,74 @@
 // Templated function to hold a subset of the columns fixed and check
 // if the solver converges to the optimal values or not.
 template<bool col1, bool col2, bool col3, bool col4>
-void IsSolveSuccessful() {
-  LevenbergMarquardt lm;
+void IsLevenbergMarquardtSolveSuccessful() {
   Solver::Options solver_options;
+  LinearSolver::Options linear_solver_options;
+  DenseQRSolver linear_solver(linear_solver_options);
+
+  double parameters[4] = { 3, -1, 0, 1.0 };
+
+  // If the column is inactive, then set its value to the optimal
+  // value.
+  parameters[0] = (col1 ? parameters[0] : 0.0);
+  parameters[1] = (col2 ? parameters[1] : 0.0);
+  parameters[2] = (col3 ? parameters[2] : 0.0);
+  parameters[3] = (col4 ? parameters[3] : 0.0);
+
+  PowellEvaluator2<col1, col2, col3, col4> powell_evaluator;
+  scoped_ptr<SparseMatrix> jacobian(powell_evaluator.CreateJacobian());
+
   Minimizer::Options minimizer_options(solver_options);
   minimizer_options.gradient_tolerance = 1e-26;
   minimizer_options.function_tolerance = 1e-26;
   minimizer_options.parameter_tolerance = 1e-26;
-  LinearSolver::Options linear_solver_options;
-  DenseQRSolver linear_solver(linear_solver_options);
+  minimizer_options.evaluator = &powell_evaluator;
+  minimizer_options.jacobian = jacobian.get();
 
-  double initial_parameters[4] = { 3, -1, 0, 1.0 };
-  double final_parameters[4] = { -1.0, -1.0, -1.0, -1.0 };
+  TrustRegionStrategy::Options trust_region_strategy_options;
+  trust_region_strategy_options.linear_solver = &linear_solver;
+  trust_region_strategy_options.initial_radius = 1e4;
+  trust_region_strategy_options.max_radius = 1e20;
+  trust_region_strategy_options.lm_min_diagonal = 1e-6;
+  trust_region_strategy_options.lm_max_diagonal = 1e32;
+  scoped_ptr<TrustRegionStrategy> strategy(
+      TrustRegionStrategy::Create(trust_region_strategy_options));
+  minimizer_options.trust_region_strategy = strategy.get();
 
-  // If the column is inactive, then set its value to the optimal
-  // value.
-  initial_parameters[0] = (col1 ? initial_parameters[0] : 0.0);
-  initial_parameters[1] = (col2 ? initial_parameters[1] : 0.0);
-  initial_parameters[2] = (col3 ? initial_parameters[2] : 0.0);
-  initial_parameters[3] = (col4 ? initial_parameters[3] : 0.0);
-
-  PowellEvaluator2<col1, col2, col3, col4> powell_evaluator;
-
+  TrustRegionMinimizer minimizer;
   Solver::Summary summary;
-  lm.Minimize(minimizer_options,
-              &powell_evaluator,
-              &linear_solver,
-              initial_parameters,
-              final_parameters,
-              &summary);
+  minimizer.Minimize(minimizer_options, parameters, &summary);
 
   // The minimum is at x1 = x2 = x3 = x4 = 0.
-  EXPECT_NEAR(0.0, final_parameters[0], 0.001);
-  EXPECT_NEAR(0.0, final_parameters[1], 0.001);
-  EXPECT_NEAR(0.0, final_parameters[2], 0.001);
-  EXPECT_NEAR(0.0, final_parameters[3], 0.001);
+  EXPECT_NEAR(0.0, parameters[0], 0.001);
+  EXPECT_NEAR(0.0, parameters[1], 0.001);
+  EXPECT_NEAR(0.0, parameters[2], 0.001);
+  EXPECT_NEAR(0.0, parameters[3], 0.001);
 };
 
-TEST(LevenbergMarquardt, PowellsSingularFunction) {
+TEST(TrustRegionMinimizer, PowellsSingularFunction) {
   // This case is excluded because this has a local minimum and does
   // not find the optimum. This should not affect the correctness of
   // this test since we are testing all the other 14 combinations of
   // column activations.
+  //
+  //   IsSolveSuccessful<true, true, false, true>();
 
-  // IsSolveSuccessful<true, true, false, true>();
-
-  IsSolveSuccessful<true,  true,  true,  true>();
-  IsSolveSuccessful<true,  true,  true,  false>();
-  IsSolveSuccessful<true,  false, true,  true>();
-  IsSolveSuccessful<false, true,  true,  true>();
-  IsSolveSuccessful<true,  true,  false, false>();
-  IsSolveSuccessful<true,  false, true,  false>();
-  IsSolveSuccessful<false, true,  true,  false>();
-  IsSolveSuccessful<true,  false, false, true>();
-  IsSolveSuccessful<false, true,  false, true>();
-  IsSolveSuccessful<false, false, true,  true>();
-  IsSolveSuccessful<true,  false, false, false>();
-  IsSolveSuccessful<false, true,  false, false>();
-  IsSolveSuccessful<false, false, true,  false>();
-  IsSolveSuccessful<false, false, false, true>();
+  IsLevenbergMarquardtSolveSuccessful<true,  true,  true,  true>();
+  IsLevenbergMarquardtSolveSuccessful<true,  true,  true,  false>();
+  IsLevenbergMarquardtSolveSuccessful<true,  false, true,  true>();
+  IsLevenbergMarquardtSolveSuccessful<false, true,  true,  true>();
+  IsLevenbergMarquardtSolveSuccessful<true,  true,  false, false>();
+  IsLevenbergMarquardtSolveSuccessful<true,  false, true,  false>();
+  IsLevenbergMarquardtSolveSuccessful<false, true,  true,  false>();
+  IsLevenbergMarquardtSolveSuccessful<true,  false, false, true>();
+  IsLevenbergMarquardtSolveSuccessful<false, true,  false, true>();
+  IsLevenbergMarquardtSolveSuccessful<false, false, true,  true>();
+  IsLevenbergMarquardtSolveSuccessful<true,  false, false, false>();
+  IsLevenbergMarquardtSolveSuccessful<false, true,  false, false>();
+  IsLevenbergMarquardtSolveSuccessful<false, false, true,  false>();
+  IsLevenbergMarquardtSolveSuccessful<false, false, false, true>();
 }
 
-
 }  // namespace internal
 }  // namespace ceres
diff --git a/internal/ceres/trust_region_strategy.cc b/internal/ceres/trust_region_strategy.cc
new file mode 100644
index 0000000..5c69bd4
--- /dev/null
+++ b/internal/ceres/trust_region_strategy.cc
@@ -0,0 +1,23 @@
+#include "ceres/trust_region_strategy.h"
+#include "ceres/levenberg_marquardt_strategy.h"
+
+namespace ceres {
+namespace internal {
+
+TrustRegionStrategy::~TrustRegionStrategy() {}
+
+TrustRegionStrategy* TrustRegionStrategy::Create(const Options& options) {
+  switch (options.trust_region_strategy_type) {
+    case LEVENBERG_MARQUARDT:
+      return new LevenbergMarquardtStrategy(options);
+    default:
+      LOG(FATAL) << "Unknown trust region strategy: "
+                 << options.trust_region_strategy_type;
+  }
+
+  LOG(FATAL) << "Unknown trust region strategy: "
+             << options.trust_region_strategy_type;
+}
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/internal/ceres/trust_region_strategy.h b/internal/ceres/trust_region_strategy.h
new file mode 100644
index 0000000..afa4400
--- /dev/null
+++ b/internal/ceres/trust_region_strategy.h
@@ -0,0 +1,109 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// 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_STRATEGY_H_
+#define CERES_INTERNAL_TRUST_REGION_STRATEGY_H_
+
+#include "ceres/linear_solver.h"
+
+namespace ceres {
+namespace internal {
+
+// Interface for classes implementing various trust region strategies
+// for nonlinear least squares problems.
+//
+// The object is expected to maintain and update a trust region
+// radius, which it then uses to solve for the trust region step using
+// the jacobian matrix and residual vector.
+//
+// Here the term trust region radius is used loosely, as the strategy
+// is free to treat it as guidance and violate it as need be. e.g.,
+// the LevenbergMarquardtStrategy uses the inverse of the trust region
+// radius to scale the damping term, which controls the step size, but
+// does not set a hard limit on its size.
+class TrustRegionStrategy {
+public:
+  struct Options {
+    Options()
+        : trust_region_strategy_type(LEVENBERG_MARQUARDT),
+          initial_radius(1e4),
+          max_radius(1e32),
+          lm_min_diagonal(1e-6),
+          lm_max_diagonal(1e32) {
+    }
+
+    TrustRegionStrategyType trust_region_strategy_type;
+    // Linear solver used for actually solving the trust region step.
+    LinearSolver* linear_solver;
+    double initial_radius;
+    double max_radius;
+
+    // Minimum and maximum values of the diagonal damping matrix used
+    // by LevenbergMarquardtStrategy.
+    double lm_min_diagonal;
+    double lm_max_diagonal;
+  };
+
+  // Per solve options.
+  struct PerSolveOptions {
+    // Forcing sequence for inexact solves.
+    double eta;
+  };
+
+  ~TrustRegionStrategy();
+
+  // Use the current radius to solve for the trust region step.
+  virtual LinearSolver::Summary ComputeStep(
+      const PerSolveOptions& per_solve_options,
+      SparseMatrix* jacobian,
+      const double* residuals,
+      double* step) = 0;
+
+  // Inform the strategy that the current step has been accepted, and
+  // that the ratio of the decrease in the non-linear objective to the
+  // decrease in the trust region model is step_quality.
+  virtual void StepAccepted(double step_quality) = 0;
+
+  // Inform the strategy that the current step has been rejected, and
+  // that the ratio of the decrease in the non-linear objective to the
+  // decrease in the trust region model is step_quality.
+  virtual void StepRejected(double step_quality) = 0;
+
+  // Current trust region radius.
+  virtual double Radius() const = 0;
+
+  // Factory.
+  static TrustRegionStrategy* Create(const Options& options);
+};
+
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // CERES_INTERNAL_TRUST_REGION_STRATEGY_H_