diff --git a/docs/source/bibliography.rst b/docs/source/bibliography.rst
index bab6e61..ab2a29f 100644
--- a/docs/source/bibliography.rst
+++ b/docs/source/bibliography.rst
@@ -94,6 +94,9 @@
    Part II: Implementation and Experiments**, Management Science,
    20(5), 863-874, 1974.
 
+.. [Ridders] C. J. F. Ridders, **Accurate computation of F'(x) and
+   F'(x) F"(x)**, Advances in Engineering Software 4(2), 75-76, 1978.
+
 .. [RuheWedin] A. Ruhe and P.Å. Wedin, **Algorithms for separable
    nonlinear least squares problems**, Siam Review, 22(3):318–337,
    1980.
diff --git a/docs/source/nnls_modeling.rst b/docs/source/nnls_modeling.rst
index 408544d..038684d 100644
--- a/docs/source/nnls_modeling.rst
+++ b/docs/source/nnls_modeling.rst
@@ -367,7 +367,7 @@
     .. code-block:: c++
 
       template <typename CostFunctor,
-                NumericDiffMethod method = CENTRAL,
+                NumericDiffMethodType method = CENTRAL,
                 int kNumResiduals,  // Number of residuals, or ceres::DYNAMIC.
                 int N0,       // Number of parameters in block 0.
                 int N1 = 0,   // Number of parameters in block 1.
@@ -481,10 +481,30 @@
    independent variables, and there is no limit on the dimensionality
    of each of them.
 
-   The ``CENTRAL`` difference method is considerably more accurate at
+   There are three available numeric differentiation schemes in ceres-solver:
+
+   The ``FORWARD`` difference method, which approximates :math:`f'(x)`
+   by computing :math:`\frac{f(x+h)-f(x)}{h}`, computes the cost function
+   one additional time at :math:`x+h`. It is the fastest but least accurate
+   method.
+
+   The ``CENTRAL`` difference method is more accurate at
    the cost of twice as many function evaluations than forward
-   difference. Consider using central differences begin with, and only
-   after that works, trying forward difference to improve performance.
+   difference, estimating :math:`f'(x)` by computing
+   :math:`\frac{f(x+h)-f(x-h)}{2h}`.
+
+   The ``RIDDERS`` difference method[Ridders]_ is an adaptive scheme that
+   estimates derivatives by performing multiple central differences
+   at varying scales. Specifically, the algorithm starts at a certain
+   :math:`h` and as the derivative is estimated, this step size decreases.
+   To conserve function evaluations and estimate the derivative error, the 
+   method performs Richardson extrapolations between the tested step sizes.
+   The algorithm exhibits considerably higher accuracy, but does so by
+   additional evaluations of the cost function.
+
+   Consider using ``CENTRAL`` differences to begin with. Based on the
+   results, either try forward difference to improve performance or
+   Ridders' method to improve accuracy.
 
    **WARNING** A common beginner's error when first using
    NumericDiffCostFunction is to get the sizing wrong. In particular,
@@ -530,7 +550,7 @@
 
    **Alternate Interface**
 
-   For a variety of reason, including compatibility with legacy code,
+   For a variety of reasons, including compatibility with legacy code,
    :class:`NumericDiffCostFunction` can also take
    :class:`CostFunction` objects as input. The following describes
    how.
@@ -571,7 +591,7 @@
 
      .. code-block:: c++
 
-      template <typename CostFunctor, NumericDiffMethod method = CENTRAL>
+      template <typename CostFunctor, NumericDiffMethodType method = CENTRAL>
       class DynamicNumericDiffCostFunction : public CostFunction {
       };
 
diff --git a/examples/more_garbow_hillstrom.cc b/examples/more_garbow_hillstrom.cc
index b0ea78a..835cb7f 100644
--- a/examples/more_garbow_hillstrom.cc
+++ b/examples/more_garbow_hillstrom.cc
@@ -63,34 +63,62 @@
 DEFINE_bool(use_numeric_diff, false,
             "Use numeric differentiation instead of automatic "
             "differentiation.");
-
+DEFINE_string(numeric_diff_method, "ridders", "When using numeric "
+              "differentiation, selects algorithm. Options are: central, "
+              "forward, ridders.");
+DEFINE_int32(ridders_extrapolations, 3, "Maximal number of extrapolations in "
+             "Ridders' method.");
 
 namespace ceres {
 namespace examples {
 
 const double kDoubleMax = std::numeric_limits<double>::max();
 
-#define BEGIN_MGH_PROBLEM(name, num_parameters, num_residuals)          \
-  struct name {                                                         \
-    static const int kNumParameters = num_parameters;                   \
-    static const double initial_x[kNumParameters];                      \
-    static const double lower_bounds[kNumParameters];                   \
-    static const double upper_bounds[kNumParameters];                   \
-    static const double constrained_optimal_cost;                       \
-    static const double unconstrained_optimal_cost;                     \
-    static CostFunction* Create() {                                     \
-      if (FLAGS_use_numeric_diff) {                                     \
-        return new NumericDiffCostFunction<name,                        \
-                                           CENTRAL,                     \
-                                           num_residuals,               \
-                                           num_parameters>(new name);   \
-      } else {                                                          \
-        return new AutoDiffCostFunction<name,                           \
-                                        num_residuals,                  \
-                                        num_parameters>(new name);      \
-      }                                                                 \
-    }                                                                   \
-    template <typename T>                                               \
+static void SetNumericDiffOptions(ceres::NumericDiffOptions* options) {
+  options->max_num_ridders_extrapolations = FLAGS_ridders_extrapolations;
+}
+
+#define BEGIN_MGH_PROBLEM(name, num_parameters, num_residuals)            \
+  struct name {                                                           \
+    static const int kNumParameters = num_parameters;                     \
+    static const double initial_x[kNumParameters];                        \
+    static const double lower_bounds[kNumParameters];                     \
+    static const double upper_bounds[kNumParameters];                     \
+    static const double constrained_optimal_cost;                         \
+    static const double unconstrained_optimal_cost;                       \
+    static CostFunction* Create() {                                       \
+      if (FLAGS_use_numeric_diff) {                                       \
+        ceres::NumericDiffOptions options;                                \
+        SetNumericDiffOptions(&options);                                  \
+        if (FLAGS_numeric_diff_method == "central") {                     \
+          return new NumericDiffCostFunction<name,                        \
+                                             ceres::CENTRAL,              \
+                                             num_residuals,               \
+                                             num_parameters>(             \
+              new name, ceres::TAKE_OWNERSHIP, num_residuals, options);   \
+        } else if (FLAGS_numeric_diff_method == "forward") {              \
+          return new NumericDiffCostFunction<name,                        \
+                                             ceres::FORWARD,              \
+                                             num_residuals,               \
+                                             num_parameters>(             \
+              new name, ceres::TAKE_OWNERSHIP, num_residuals, options);   \
+        } else if (FLAGS_numeric_diff_method == "ridders") {              \
+          return new NumericDiffCostFunction<name,                        \
+                                             ceres::RIDDERS,              \
+                                             num_residuals,               \
+                                             num_parameters>(             \
+              new name, ceres::TAKE_OWNERSHIP, num_residuals, options);   \
+        } else {                                                          \
+          LOG(ERROR) << "Invalid numeric diff method specified";          \
+          return NULL;                                                    \
+        }                                                                 \
+      } else {                                                            \
+        return new AutoDiffCostFunction<name,                             \
+                                        num_residuals,                    \
+                                        num_parameters>(new name);        \
+      }                                                                   \
+    }                                                                     \
+    template <typename T>                                                 \
     bool operator()(const T* const x, T* residual) const {
 
 #define END_MGH_PROBLEM return true; } };  // NOLINT
diff --git a/examples/nist.cc b/examples/nist.cc
index ceeb457..f7bd38e 100644
--- a/examples/nist.cc
+++ b/examples/nist.cc
@@ -65,7 +65,7 @@
 // solver had the highest LRE.
 
 // In this file, we implement the same evaluation methodology using
-// Ceres. Currently using Levenberg-Marquard with DENSE_QR, we get
+// Ceres. Currently using Levenberg-Marquardt with DENSE_QR, we get
 //
 //               Excel  Gnuplot  GaussFit  HBN  MinPack  Ceres
 // Average LRE     2.3      4.3       4.0  6.8      4.4    9.4
@@ -118,6 +118,13 @@
 DEFINE_bool(use_numeric_diff, false,
             "Use numeric differentiation instead of automatic "
             "differentiation.");
+DEFINE_string(numeric_diff_method, "ridders", "When using numeric "
+              "differentiation, selects algorithm. Options are: central, "
+              "forward, ridders.");
+DEFINE_double(ridders_step_size, 1e-9, "Initial step size for Ridders "
+              "numeric differentiation.");
+DEFINE_int32(ridders_extrapolations, 3, "Maximal number of Ridders "
+             "extrapolations.");
 
 namespace ceres {
 namespace examples {
@@ -413,6 +420,11 @@
   double y_;
 };
 
+static void SetNumericDiffOptions(ceres::NumericDiffOptions* options) {
+  options->max_num_ridders_extrapolations = FLAGS_ridders_extrapolations;
+  options->ridders_relative_initial_step_size = FLAGS_ridders_step_size;
+}
+
 template <typename Model, int num_residuals, int num_parameters>
 int RegressionDriver(const string& filename,
                      const ceres::Solver::Options& options) {
@@ -439,11 +451,30 @@
           response.data() + nist_problem.response_size() * i);
       ceres::CostFunction* cost_function = NULL;
       if (FLAGS_use_numeric_diff) {
-        cost_function =
-            new ceres::NumericDiffCostFunction<Model,
-                                               ceres::CENTRAL,
-                                               num_residuals,
-                                               num_parameters>(model);
+        ceres::NumericDiffOptions options;
+        SetNumericDiffOptions(&options);
+        if (FLAGS_numeric_diff_method == "central") {
+          cost_function = new NumericDiffCostFunction<Model,
+                                                      ceres::CENTRAL,
+                                                      num_residuals,
+                                                      num_parameters>(
+              model, ceres::TAKE_OWNERSHIP, num_residuals, options);
+        } else if (FLAGS_numeric_diff_method == "forward") {
+          cost_function = new NumericDiffCostFunction<Model,
+                                                      ceres::FORWARD,
+                                                      num_residuals,
+                                                      num_parameters>(
+              model, ceres::TAKE_OWNERSHIP, num_residuals, options);
+        } else if (FLAGS_numeric_diff_method == "ridders") {
+          cost_function = new NumericDiffCostFunction<Model,
+                                                      ceres::RIDDERS,
+                                                      num_residuals,
+                                                      num_parameters>(
+              model, ceres::TAKE_OWNERSHIP, num_residuals, options);
+        } else {
+          LOG(ERROR) << "Invalid numeric diff method specified";
+          return 0;
+        }
       } else {
          cost_function =
              new ceres::AutoDiffCostFunction<Model,
diff --git a/include/ceres/dynamic_numeric_diff_cost_function.h b/include/ceres/dynamic_numeric_diff_cost_function.h
index 72fa13b..c852d57 100644
--- a/include/ceres/dynamic_numeric_diff_cost_function.h
+++ b/include/ceres/dynamic_numeric_diff_cost_function.h
@@ -68,19 +68,37 @@
 #include "ceres/internal/scoped_ptr.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/internal/numeric_diff.h"
+#include "ceres/numeric_diff_options.h"
 #include "glog/logging.h"
 
 namespace ceres {
 
-template <typename CostFunctor, NumericDiffMethod method = CENTRAL>
+template <typename CostFunctor, NumericDiffMethodType method = CENTRAL>
 class DynamicNumericDiffCostFunction : public CostFunction {
  public:
-  explicit DynamicNumericDiffCostFunction(const CostFunctor* functor,
-                                          Ownership ownership = TAKE_OWNERSHIP,
-                                          double relative_step_size = 1e-6)
+  explicit DynamicNumericDiffCostFunction(
+      const CostFunctor* functor,
+      Ownership ownership = TAKE_OWNERSHIP,
+      const NumericDiffOptions& options = NumericDiffOptions())
       : functor_(functor),
         ownership_(ownership),
-        relative_step_size_(relative_step_size) {
+        options_(options) {
+  }
+
+  // Deprecated. New users should avoid using this constructor. Instead, use the
+  // constructor with NumericDiffOptions.
+  DynamicNumericDiffCostFunction(
+      const CostFunctor* functor,
+      Ownership ownership,
+      double relative_step_size)
+      : functor_(functor),
+        ownership_(ownership),
+        options_() {
+    LOG(WARNING) << "This constructor is deprecated and will be removed in "
+                    "a future version. Please use the NumericDiffOptions "
+                    "constructor instead.";
+
+    options_.relative_step_size = relative_step_size;
   }
 
   virtual ~DynamicNumericDiffCostFunction() {
@@ -140,7 +158,7 @@
                        DYNAMIC, DYNAMIC>::EvaluateJacobianForParameterBlock(
                                              functor_.get(),
                                              residuals,
-                                             relative_step_size_,
+                                             options_,
                                              this->num_residuals(),
                                              block,
                                              block_sizes[block],
@@ -179,7 +197,7 @@
 
   internal::scoped_ptr<const CostFunctor> functor_;
   Ownership ownership_;
-  const double relative_step_size_;
+  NumericDiffOptions options_;
 };
 
 }  // namespace ceres
diff --git a/include/ceres/internal/numeric_diff.h b/include/ceres/internal/numeric_diff.h
index b9aab88..2fffdff 100644
--- a/include/ceres/internal/numeric_diff.h
+++ b/include/ceres/internal/numeric_diff.h
@@ -28,8 +28,9 @@
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
 //         mierle@gmail.com (Keir Mierle)
+//         tbennun@gmail.com (Tal Ben-Nun)
 //
-// Finite differencing routine used by NumericDiffCostFunction.
+// Finite differencing routines used by NumericDiffCostFunction.
 
 #ifndef CERES_PUBLIC_INTERNAL_NUMERIC_DIFF_H_
 #define CERES_PUBLIC_INTERNAL_NUMERIC_DIFF_H_
@@ -37,9 +38,12 @@
 #include <cstring>
 
 #include "Eigen/Dense"
+#include "Eigen/StdVector"
 #include "ceres/cost_function.h"
+#include "ceres/internal/fixed_array.h"
 #include "ceres/internal/scoped_ptr.h"
 #include "ceres/internal/variadic_evaluate.h"
+#include "ceres/numeric_diff_options.h"
 #include "ceres/types.h"
 #include "glog/logging.h"
 
@@ -78,7 +82,7 @@
 // specializations for member functions. The alternative is to repeat the main
 // class for differing numbers of parameters, which is also unfortunate.
 template <typename CostFunctor,
-          NumericDiffMethod kMethod,
+          NumericDiffMethodType kMethod,
           int kNumResiduals,
           int N0, int N1, int N2, int N3, int N4,
           int N5, int N6, int N7, int N8, int N9,
@@ -88,8 +92,8 @@
   // Mutates parameters but must restore them before return.
   static bool EvaluateJacobianForParameterBlock(
       const CostFunctor* functor,
-      double const* residuals_at_eval_point,
-      const double relative_step_size,
+      const double* residuals_at_eval_point,
+      const NumericDiffOptions& options,
       int num_residuals,
       int parameter_block_index,
       int parameter_block_size,
@@ -126,67 +130,285 @@
                                            num_residuals_internal,
                                            parameter_block_size_internal);
 
-    // Mutate 1 element at a time and then restore.
     Map<ParameterVector> x_plus_delta(
         parameters[parameter_block_index_internal],
         parameter_block_size_internal);
     ParameterVector x(x_plus_delta);
-    ParameterVector step_size = x.array().abs() * relative_step_size;
+    ParameterVector step_size = x.array().abs() *
+        ((kMethod == RIDDERS) ? options.ridders_relative_initial_step_size :
+        options.relative_step_size);
 
     // It is not a good idea to make the step size arbitrarily
     // small. This will lead to problems with round off and numerical
     // instability when dividing by the step size. The general
     // recommendation is to not go down below sqrt(epsilon).
-    const double min_step_size =
-        std::sqrt(std::numeric_limits<double>::epsilon());
+    double min_step_size = std::sqrt(std::numeric_limits<double>::epsilon());
+
+    // For Ridders' method, the initial step size is required to be large,
+    // thus ridders_relative_initial_step_size is used.
+    if (kMethod == RIDDERS) {
+      min_step_size = std::max(min_step_size,
+                               options.ridders_relative_initial_step_size);
+    }
 
     // For each parameter in the parameter block, use finite differences to
     // compute the derivative for that parameter.
-    ResidualVector residuals(num_residuals_internal);
+    FixedArray<double> temp_residual_array(num_residuals_internal);
+    FixedArray<double> residual_array(num_residuals_internal);
+    Map<ResidualVector> residuals(residual_array.get(),
+                                  num_residuals_internal);
+
     for (int j = 0; j < parameter_block_size_internal; ++j) {
       const double delta = std::max(min_step_size, step_size(j));
-      x_plus_delta(j) = x(j) + delta;
+
+      if (kMethod == RIDDERS) {
+        if (!EvaluateRiddersJacobianColumn(functor, j, delta,
+                                           options,
+                                           num_residuals_internal,
+                                           parameter_block_size_internal,
+                                           x.data(),
+                                           residuals_at_eval_point,
+                                           parameters,
+                                           x_plus_delta.data(),
+                                           temp_residual_array.get(),
+                                           residual_array.get())) {
+          return false;
+        }
+      } else {
+        if (!EvaluateJacobianColumn(functor, j, delta, 
+                                    num_residuals_internal,
+                                    parameter_block_size_internal,
+                                    x.data(),
+                                    residuals_at_eval_point,
+                                    parameters, 
+                                    x_plus_delta.data(),
+                                    temp_residual_array.get(),
+                                    residual_array.get())) {
+          return false;
+        }
+      }
+
+      parameter_jacobian.col(j).matrix() = residuals;
+    }
+    return true;
+  }
+
+  static bool EvaluateJacobianColumn(const CostFunctor* functor,
+                                     int parameter_index, 
+                                     double delta,
+                                     int num_residuals,
+                                     int parameter_block_size,
+                                     const double* x_ptr,
+                                     const double* residuals_at_eval_point,
+                                     double** parameters,
+                                     double* x_plus_delta_ptr,
+                                     double* temp_residuals_ptr,
+                                     double* residuals_ptr) {
+    using Eigen::Map;
+    using Eigen::Matrix;
+
+    typedef Matrix<double, kNumResiduals, 1> ResidualVector;
+    typedef Matrix<double, kParameterBlockSize, 1> ParameterVector;
+
+    Map<const ParameterVector> x(x_ptr, parameter_block_size);
+    Map<ParameterVector> x_plus_delta(x_plus_delta_ptr,
+                                      parameter_block_size);
+
+    Map<ResidualVector> residuals(residuals_ptr, num_residuals);
+    Map<ResidualVector> temp_residuals(temp_residuals_ptr, num_residuals);
+
+    // Mutate 1 element at a time and then restore.
+    x_plus_delta(parameter_index) = x(parameter_index) + delta;
+
+    if (!EvaluateImpl<CostFunctor, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>(
+            functor, parameters, residuals.data(), functor)) {
+      return false;
+    }
+
+    // Compute this column of the jacobian in 3 steps:
+    // 1. Store residuals for the forward part.
+    // 2. Subtract residuals for the backward (or 0) part.
+    // 3. Divide out the run.
+    double one_over_delta = 1.0 / delta;
+    if (kMethod == CENTRAL || kMethod == RIDDERS) {
+      // Compute the function on the other side of x(parameter_index).
+      x_plus_delta(parameter_index) = x(parameter_index) - delta;
 
       if (!EvaluateImpl<CostFunctor, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>(
-              functor, parameters, residuals.data(), functor)) {
+              functor, parameters, temp_residuals.data(), functor)) {
         return false;
       }
 
-      // Compute this column of the jacobian in 3 steps:
-      // 1. Store residuals for the forward part.
-      // 2. Subtract residuals for the backward (or 0) part.
-      // 3. Divide out the run.
-      parameter_jacobian.col(j) = residuals;
+      residuals -= temp_residuals;
+      one_over_delta /= 2;
+    } else {
+      // Forward difference only; reuse existing residuals evaluation.
+      residuals -=
+          Map<const ResidualVector>(residuals_at_eval_point,
+                                    num_residuals);
+    }
 
-      double one_over_delta = 1.0 / delta;
-      if (kMethod == CENTRAL) {
-        // Compute the function on the other side of x(j).
-        x_plus_delta(j) = x(j) - delta;
+    // Restore x_plus_delta.
+    x_plus_delta(parameter_index) = x(parameter_index);
 
-        if (!EvaluateImpl<CostFunctor, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>(
-                functor, parameters, residuals.data(), functor)) {
-          return false;
-        }
+    // Divide out the run to get slope.
+    residuals *= one_over_delta;
 
-        parameter_jacobian.col(j) -= residuals;
-        one_over_delta /= 2;
-      } else {
-        // Forward difference only; reuse existing residuals evaluation.
-        parameter_jacobian.col(j) -=
-            Map<const ResidualVector>(residuals_at_eval_point,
-                                      num_residuals_internal);
+    return true;
+  }
+
+  // This numeric difference implementation uses adaptive differentiation
+  // on the parameters to obtain the Jacobian matrix. The adaptive algorithm
+  // is based on Ridders' method for adaptive differentiation, which creates
+  // a Romberg tableau from varying step sizes and extrapolates the
+  // intermediate results to obtain the current computational error.
+  //
+  // References:
+  // C.J.F. Ridders, Accurate computation of F'(x) and F'(x) F"(x), Advances
+  // in Engineering Software (1978), Volume 4, Issue 2, April 1982,
+  // Pages 75-76, ISSN 0141-1195,
+  // http://dx.doi.org/10.1016/S0141-1195(82)80057-0.
+  static bool EvaluateRiddersJacobianColumn(
+      const CostFunctor* functor,
+      int parameter_index, 
+      double delta,
+      const NumericDiffOptions& options,
+      int num_residuals,
+      int parameter_block_size,
+      const double* x_ptr,
+      const double* residuals_at_eval_point,
+      double** parameters,
+      double* x_plus_delta_ptr,
+      double* temp_residuals_ptr,
+      double* residuals_ptr) {
+    using Eigen::Map;
+    using Eigen::Matrix;
+    using Eigen::aligned_allocator;
+
+    typedef Matrix<double, kNumResiduals, 1> ResidualVector;
+    typedef Matrix<double, kNumResiduals, DYNAMIC> ResidualCandidateMatrix;
+    typedef Matrix<double, kParameterBlockSize, 1> ParameterVector;
+
+    Map<const ParameterVector> x(x_ptr, parameter_block_size);
+    Map<ParameterVector> x_plus_delta(x_plus_delta_ptr,
+                                      parameter_block_size);
+
+    Map<ResidualVector> residuals(residuals_ptr, num_residuals);
+    Map<ResidualVector> temp_residuals(temp_residuals_ptr, num_residuals);
+
+    // In order for the algorithm to converge, the step size should be
+    // initialized to a value that is large enough to produce a significant
+    // change in the function.
+    // As the derivative is estimated, the step size decreases.
+    // By default, the step sizes are chosen so that the middle column
+    // of the Romberg tableau uses the input delta.
+    double current_step_size = delta *
+        pow(options.ridders_step_shrink_factor,
+            options.max_num_ridders_extrapolations / 2);
+
+    // Double-buffering temporary differential candidate vectors
+    // from previous step size.
+    ResidualCandidateMatrix stepsize_candidates_a(
+        num_residuals,
+        options.max_num_ridders_extrapolations);
+    ResidualCandidateMatrix stepsize_candidates_b(
+        num_residuals,
+        options.max_num_ridders_extrapolations);
+    ResidualCandidateMatrix* current_candidates = &stepsize_candidates_a;
+    ResidualCandidateMatrix* previous_candidates = &stepsize_candidates_b;
+
+    // Represents the computational error of the derivative. This variable is
+    // initially set to a large value, and is set to the difference between
+    // current and previous finite difference extrapolations.
+    // norm_error is supposed to decrease as the finite difference tableau
+    // generation progresses, serving both as an estimate for differentiation
+    // error and as a measure of differentiation numerical stability.
+    double norm_error = std::numeric_limits<double>::max();
+
+    // Loop over decreasing step sizes until:
+    //  1. Error is smaller than a given value (ridders_epsilon),
+    //  2. Maximal order of extrapolation reached, or
+    //  3. Extrapolation becomes numerically unstable.
+    for (int i = 0; i < options.max_num_ridders_extrapolations; ++i) {
+      // Compute the numerical derivative at this step size.
+      if (!EvaluateJacobianColumn(functor, parameter_index, current_step_size,
+                                  num_residuals,
+                                  parameter_block_size,
+                                  x.data(),
+                                  residuals_at_eval_point,
+                                  parameters,
+                                  x_plus_delta.data(),
+                                  temp_residuals.data(),
+                                  current_candidates->col(0).data())) {
+        // Something went wrong; bail.
+        return false;
       }
-      x_plus_delta(j) = x(j);  // Restore x_plus_delta.
 
-      // Divide out the run to get slope.
-      parameter_jacobian.col(j) *= one_over_delta;
+      // Store initial results.
+      if (i == 0) {
+        residuals = current_candidates[0];
+      }
+
+      // Shrink differentiation step size.
+      current_step_size /= options.ridders_step_shrink_factor;
+
+      // Extrapolation factor for Richardson acceleration method (see below).
+      double richardson_factor = options.ridders_step_shrink_factor *
+          options.ridders_step_shrink_factor;
+      for (int k = 1; k <= i; ++k) {
+        // Extrapolate the various orders of finite differences using
+        // the Richardson acceleration method.
+        current_candidates->col(k) =
+            (richardson_factor * current_candidates->col(k - 1) -
+             previous_candidates->col(k - 1)) / (richardson_factor - 1.0);
+
+        richardson_factor *= options.ridders_step_shrink_factor *
+            options.ridders_step_shrink_factor;
+
+        // Compute the difference between the previous value and the current.
+        double candidate_error = std::max(
+            (current_candidates->col(k) -
+             current_candidates->col(k - 1)).norm(),
+            (current_candidates->col(k) -
+             previous_candidates->col(k - 1)).norm());
+
+        // If the error has decreased, update results.
+        if (candidate_error <= norm_error) {
+          norm_error = candidate_error;
+          residuals = current_candidates->col(k);
+
+          // If the error is small enough, stop.
+          if (norm_error < options.ridders_epsilon) {
+            break;
+          }
+        }
+      }
+
+      // After breaking out of the inner loop, declare convergence.
+      if (norm_error < options.ridders_epsilon) {
+        break;
+      }
+
+      // Check to see if the current gradient estimate is numerically unstable.
+      // If so, bail out and return the last stable result.
+      if (i > 0) {
+        double tableau_error = (current_candidates->col(i) -
+            previous_candidates->col(i - 1)).norm();
+
+        // Compare current error to the chosen candidate's error.
+        if (tableau_error >= 2 * norm_error) {
+          break;
+        }
+      }
+
+      std::swap(current_candidates, previous_candidates);
     }
     return true;
   }
 };
 
 template <typename CostFunctor,
-          NumericDiffMethod kMethod,
+          NumericDiffMethodType kMethod,
           int kNumResiduals,
           int N0, int N1, int N2, int N3, int N4,
           int N5, int N6, int N7, int N8, int N9,
@@ -197,8 +419,8 @@
   // Mutates parameters but must restore them before return.
   static bool EvaluateJacobianForParameterBlock(
       const CostFunctor* functor,
-      double const* residuals_at_eval_point,
-      const double relative_step_size,
+      const double* residuals_at_eval_point,
+      const NumericDiffOptions& options,
       const int num_residuals,
       const int parameter_block_index,
       const int parameter_block_size,
@@ -207,7 +429,7 @@
     // Silence unused parameter compiler warnings.
     (void)functor;
     (void)residuals_at_eval_point;
-    (void)relative_step_size;
+    (void)options;
     (void)num_residuals;
     (void)parameter_block_index;
     (void)parameter_block_size;
diff --git a/include/ceres/numeric_diff_cost_function.h b/include/ceres/numeric_diff_cost_function.h
index c15e06d..fa96078 100644
--- a/include/ceres/numeric_diff_cost_function.h
+++ b/include/ceres/numeric_diff_cost_function.h
@@ -132,7 +132,7 @@
 //
 // ALTERNATE INTERFACE
 //
-// For a variety of reason, including compatibility with legacy code,
+// For a variety of reasons, including compatibility with legacy code,
 // NumericDiffCostFunction can also take CostFunction objects as
 // input. The following describes how.
 //
@@ -165,6 +165,7 @@
 #include "ceres/cost_function.h"
 #include "ceres/internal/numeric_diff.h"
 #include "ceres/internal/scoped_ptr.h"
+#include "ceres/numeric_diff_options.h"
 #include "ceres/sized_cost_function.h"
 #include "ceres/types.h"
 #include "glog/logging.h"
@@ -172,7 +173,7 @@
 namespace ceres {
 
 template <typename CostFunctor,
-          NumericDiffMethod method = CENTRAL,
+          NumericDiffMethodType method = CENTRAL,
           int kNumResiduals = 0,  // Number of residuals, or ceres::DYNAMIC
           int N0 = 0,   // Number of parameters in block 0.
           int N1 = 0,   // Number of parameters in block 1.
@@ -189,13 +190,14 @@
                                N0, N1, N2, N3, N4,
                                N5, N6, N7, N8, N9> {
  public:
-  NumericDiffCostFunction(CostFunctor* functor,
-                          Ownership ownership = TAKE_OWNERSHIP,
-                          int num_residuals = kNumResiduals,
-                          const double relative_step_size = 1e-6)
-      :functor_(functor),
-       ownership_(ownership),
-       relative_step_size_(relative_step_size) {
+  NumericDiffCostFunction(
+      CostFunctor* functor,
+      Ownership ownership = TAKE_OWNERSHIP,
+      int num_residuals = kNumResiduals,
+      const NumericDiffOptions& options = NumericDiffOptions())
+      : functor_(functor),
+        ownership_(ownership),
+        options_(options) {
     if (kNumResiduals == DYNAMIC) {
       SizedCostFunction<kNumResiduals,
                         N0, N1, N2, N3, N4,
@@ -204,6 +206,29 @@
     }
   }
 
+  // Deprecated. New users should avoid using this constructor. Instead, use the
+  // constructor with NumericDiffOptions.
+  NumericDiffCostFunction(CostFunctor* functor,
+                          Ownership ownership,
+                          int num_residuals,
+                          const double relative_step_size)
+      :functor_(functor),
+       ownership_(ownership),
+       options_() {
+    LOG(WARNING) << "This constructor is deprecated and will be removed in "
+                    "a future version. Please use the NumericDiffOptions "
+                    "constructor instead.";
+
+    if (kNumResiduals == DYNAMIC) {
+      SizedCostFunction<kNumResiduals,
+                        N0, N1, N2, N3, N4,
+                        N5, N6, N7, N8, N9>
+          ::set_num_residuals(num_residuals);
+    }
+
+    options_.relative_step_size = relative_step_size;
+  }
+
   ~NumericDiffCostFunction() {
     if (ownership_ != TAKE_OWNERSHIP) {
       functor_.release();
@@ -278,7 +303,7 @@
                        N ## block >::EvaluateJacobianForParameterBlock( \
                            functor_.get(),                              \
                            residuals,                                   \
-                           relative_step_size_,                         \
+                           options_,                                    \
                           SizedCostFunction<kNumResiduals,              \
                            N0, N1, N2, N3, N4,                          \
                            N5, N6, N7, N8, N9>::num_residuals(),        \
@@ -309,7 +334,7 @@
  private:
   internal::scoped_ptr<CostFunctor> functor_;
   Ownership ownership_;
-  const double relative_step_size_;
+  NumericDiffOptions options_;
 };
 
 }  // namespace ceres
diff --git a/include/ceres/numeric_diff_options.h b/include/ceres/numeric_diff_options.h
new file mode 100644
index 0000000..119c8a8
--- /dev/null
+++ b/include/ceres/numeric_diff_options.h
@@ -0,0 +1,79 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2015 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: tbennun@gmail.com (Tal Ben-Nun)
+//
+
+#ifndef CERES_PUBLIC_NUMERIC_DIFF_OPTIONS_H_
+#define CERES_PUBLIC_NUMERIC_DIFF_OPTIONS_H_
+
+namespace ceres {
+
+// Options pertaining to numeric differentiation (e.g., convergence criteria,
+// step sizes).
+struct CERES_EXPORT NumericDiffOptions {
+  NumericDiffOptions() {
+    relative_step_size = 1e-6;
+    ridders_relative_initial_step_size = 1e-2;
+    max_num_ridders_extrapolations = 10;
+    ridders_epsilon = 1e-12;
+    ridders_step_shrink_factor = 2.0;
+  }
+
+  // Numeric differentiation step size (multiplied by parameter block's
+  // order of magnitude). If parameters are close to zero, the step size
+  // is set to sqrt(machine_epsilon).
+  double relative_step_size;
+
+  // Initial step size for Ridders adaptive numeric differentiation (multiplied
+  // by parameter block's order of magnitude).
+  // If parameters are close to zero, Ridders' method sets the step size
+  // directly to this value. This parameter is separate from
+  // "relative_step_size" in order to set a different default value.
+  //
+  // Note: For Ridders' method to converge, the step size should be initialized
+  // to a value that is large enough to produce a significant change in the
+  // function. As the derivative is estimated, the step size decreases.
+  double ridders_relative_initial_step_size;
+
+  // Maximal number of adaptive extrapolations (sampling) in Ridders' method.
+  int max_num_ridders_extrapolations;
+
+  // Convergence criterion on extrapolation error for Ridders adaptive
+  // differentiation. The available error estimation methods are defined in
+  // NumericDiffErrorType and set in the "ridders_error_method" field.
+  double ridders_epsilon;
+
+  // The factor in which to shrink the step size with each extrapolation in
+  // Ridders' method.
+  double ridders_step_shrink_factor;
+};
+
+}  // namespace ceres
+
+#endif  // CERES_PUBLIC_NUMERIC_DIFF_OPTIONS_H_
diff --git a/include/ceres/types.h b/include/ceres/types.h
index 9b01af2..2ea4180 100644
--- a/include/ceres/types.h
+++ b/include/ceres/types.h
@@ -397,9 +397,19 @@
   DYNAMIC = -1
 };
 
-enum NumericDiffMethod {
+// The differentiation method used to compute numerical derivatives in
+// NumericDiffCostFunction and DynamicNumericDiffCostFunction.
+enum NumericDiffMethodType {
+  // Compute central finite difference: f'(x) ~ (f(x+h) - f(x-h)) / 2h.
   CENTRAL,
-  FORWARD
+
+  // Compute forward finite difference: f'(x) ~ (f(x+h) - f(x)) / h.
+  FORWARD,
+
+  // Adaptive numerical differentiation using Ridders' method. Provides more
+  // accurate and robust derivatives at the expense of additional cost
+  // function evaluations.
+  RIDDERS
 };
 
 enum LineSearchInterpolationType {
@@ -477,6 +487,12 @@
     std::string value,
     CovarianceAlgorithmType* type);
 
+CERES_EXPORT const char* NumericDiffMethodTypeToString(
+    NumericDiffMethodType type);
+CERES_EXPORT bool StringToNumericDiffMethodType(
+    std::string value,
+    NumericDiffMethodType* type);
+
 CERES_EXPORT const char* TerminationTypeToString(TerminationType type);
 
 CERES_EXPORT bool IsSchurType(LinearSolverType type);
diff --git a/internal/ceres/gradient_checking_cost_function.cc b/internal/ceres/gradient_checking_cost_function.cc
index e3ad3d4..580fd26 100644
--- a/internal/ceres/gradient_checking_cost_function.cc
+++ b/internal/ceres/gradient_checking_cost_function.cc
@@ -86,7 +86,7 @@
 class GradientCheckingCostFunction : public CostFunction {
  public:
   GradientCheckingCostFunction(const CostFunction* function,
-                               double relative_step_size,
+                               const NumericDiffOptions& options,
                                double relative_precision,
                                const string& extra_info)
       : function_(function),
@@ -97,7 +97,7 @@
         new DynamicNumericDiffCostFunction<CostFunction, CENTRAL>(
             function,
             DO_NOT_TAKE_OWNERSHIP,
-            relative_step_size);
+            options);
 
     const vector<int32>& parameter_block_sizes =
         function->parameter_block_sizes();
@@ -235,8 +235,11 @@
     double relative_step_size,
     double relative_precision,
     const string& extra_info) {
+  NumericDiffOptions numeric_diff_options;
+  numeric_diff_options.relative_step_size = relative_step_size;
+
   return new GradientCheckingCostFunction(cost_function,
-                                          relative_step_size,
+                                          numeric_diff_options,
                                           relative_precision,
                                           extra_info);
 }
diff --git a/internal/ceres/numeric_diff_cost_function_test.cc b/internal/ceres/numeric_diff_cost_function_test.cc
index cb74a86..13ab106 100644
--- a/internal/ceres/numeric_diff_cost_function_test.cc
+++ b/internal/ceres/numeric_diff_cost_function_test.cc
@@ -27,6 +27,7 @@
 // POSSIBILITY OF SUCH DAMAGE.
 //
 // Author: keir@google.com (Keir Mierle)
+//         tbennun@gmail.com (Tal Ben-Nun)
 
 #include "ceres/numeric_diff_cost_function.h"
 
@@ -71,6 +72,19 @@
   functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, FORWARD);
 }
 
+TEST(NumericDiffCostFunction, EasyCaseFunctorRidders) {
+  internal::scoped_ptr<CostFunction> cost_function;
+  cost_function.reset(
+      new NumericDiffCostFunction<EasyFunctor,
+                                  RIDDERS,
+                                  3,  /* number of residuals */
+                                  5,  /* size of x1 */
+                                  5   /* size of x2 */>(
+          new EasyFunctor));
+  EasyFunctor functor;
+  functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS);
+}
+
 TEST(NumericDiffCostFunction, EasyCaseCostFunctionCentralDifferences) {
   internal::scoped_ptr<CostFunction> cost_function;
   cost_function.reset(
@@ -97,7 +111,21 @@
   functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, FORWARD);
 }
 
-TEST(NumericDiffCostFunction, TranscendentalCaseFunctorCentralDifferences) {
+TEST(NumericDiffCostFunction, EasyCaseCostFunctionRidders) {
+  internal::scoped_ptr<CostFunction> cost_function;
+  cost_function.reset(
+      new NumericDiffCostFunction<EasyCostFunction,
+                                  RIDDERS,
+                                  3,  /* number of residuals */
+                                  5,  /* size of x1 */
+                                  5   /* size of x2 */>(
+          new EasyCostFunction, TAKE_OWNERSHIP));
+  EasyFunctor functor;
+  functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS);
+}
+
+TEST(NumericDiffCostFunction,
+     TranscendentalCaseFunctorCentralDifferences) {
   internal::scoped_ptr<CostFunction> cost_function;
   cost_function.reset(
       new NumericDiffCostFunction<TranscendentalFunctor,
@@ -110,7 +138,8 @@
   functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, CENTRAL);
 }
 
-TEST(NumericDiffCostFunction, TranscendentalCaseFunctorForwardDifferences) {
+TEST(NumericDiffCostFunction,
+     TranscendentalCaseFunctorForwardDifferences) {
   internal::scoped_ptr<CostFunction> cost_function;
   cost_function.reset(
       new NumericDiffCostFunction<TranscendentalFunctor,
@@ -123,7 +152,28 @@
   functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, FORWARD);
 }
 
-TEST(NumericDiffCostFunction, TranscendentalCaseCostFunctionCentralDifferences) {
+TEST(NumericDiffCostFunction,
+     TranscendentalCaseFunctorRidders) {
+  NumericDiffOptions options;
+
+  // Using a smaller initial step size to overcome oscillatory function
+  // behavior.
+  options.ridders_relative_initial_step_size = 1e-3;
+
+  internal::scoped_ptr<CostFunction> cost_function;
+  cost_function.reset(
+      new NumericDiffCostFunction<TranscendentalFunctor,
+                                  RIDDERS,
+                                  2,  /* number of residuals */
+                                  5,  /* size of x1 */
+                                  5   /* size of x2 */>(
+          new TranscendentalFunctor, TAKE_OWNERSHIP, 2, options));
+  TranscendentalFunctor functor;
+  functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS);
+}
+
+TEST(NumericDiffCostFunction,
+     TranscendentalCaseCostFunctionCentralDifferences) {
   internal::scoped_ptr<CostFunction> cost_function;
   cost_function.reset(
       new NumericDiffCostFunction<TranscendentalCostFunction,
@@ -136,7 +186,8 @@
   functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, CENTRAL);
 }
 
-TEST(NumericDiffCostFunction, TranscendentalCaseCostFunctionForwardDifferences) {
+TEST(NumericDiffCostFunction,
+     TranscendentalCaseCostFunctionForwardDifferences) {
   internal::scoped_ptr<CostFunction> cost_function;
   cost_function.reset(
       new NumericDiffCostFunction<TranscendentalCostFunction,
@@ -149,6 +200,26 @@
   functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, FORWARD);
 }
 
+TEST(NumericDiffCostFunction,
+     TranscendentalCaseCostFunctionRidders) {
+  NumericDiffOptions options;
+
+  // Using a smaller initial step size to overcome oscillatory function
+  // behavior.
+  options.ridders_relative_initial_step_size = 1e-3;
+
+  internal::scoped_ptr<CostFunction> cost_function;
+  cost_function.reset(
+      new NumericDiffCostFunction<TranscendentalCostFunction,
+                                  RIDDERS,
+                                  2,  /* number of residuals */
+                                  5,  /* size of x1 */
+                                  5   /* size of x2 */>(
+          new TranscendentalCostFunction, TAKE_OWNERSHIP, 2, options));
+  TranscendentalFunctor functor;
+  functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS);
+}
+
 template<int num_rows, int num_cols>
 class SizeTestingCostFunction : public SizedCostFunction<num_rows, num_cols> {
  public:
@@ -208,7 +279,8 @@
           new EasyFunctor, TAKE_OWNERSHIP, 2));
 }
 
-TEST(NumericDiffCostFunction, EasyCaseFunctorCentralDifferencesAndDynamicNumResiduals) {
+TEST(NumericDiffCostFunction,
+     EasyCaseFunctorCentralDifferencesAndDynamicNumResiduals) {
   internal::scoped_ptr<CostFunction> cost_function;
   cost_function.reset(
       new NumericDiffCostFunction<EasyFunctor,
@@ -221,5 +293,66 @@
   functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, CENTRAL);
 }
 
+TEST(NumericDiffCostFunction, ExponentialFunctorRidders) {
+  internal::scoped_ptr<CostFunction> cost_function;
+  cost_function.reset(
+      new NumericDiffCostFunction<ExponentialFunctor,
+                                  RIDDERS,
+                                  1,  /* number of residuals */
+                                  1   /* size of x1 */>(
+             new ExponentialFunctor));
+  ExponentialFunctor functor;
+  functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function);
+}
+
+TEST(NumericDiffCostFunction, ExponentialCostFunctionRidders) {
+  internal::scoped_ptr<CostFunction> cost_function;
+  cost_function.reset(
+      new NumericDiffCostFunction<ExponentialCostFunction,
+                                  RIDDERS,
+                                  1,  /* number of residuals */
+                                  1   /* size of x1 */>(
+             new ExponentialCostFunction));
+  ExponentialFunctor functor;
+  functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function);
+}
+
+TEST(NumericDiffCostFunction, RandomizedFunctorRidders) {
+  internal::scoped_ptr<CostFunction> cost_function;
+  NumericDiffOptions options;
+  // Larger initial step size is chosen to produce robust results in the
+  // presence of random noise.
+  options.ridders_relative_initial_step_size = 10.0;
+
+  cost_function.reset(
+      new NumericDiffCostFunction<RandomizedFunctor,
+                                  RIDDERS,
+                                  1,  /* number of residuals */
+                                  1   /* size of x1 */>(
+             new RandomizedFunctor(kNoiseFactor, kRandomSeed), TAKE_OWNERSHIP,
+             1, options));
+  RandomizedFunctor functor (kNoiseFactor, kRandomSeed);
+  functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function);
+}
+
+TEST(NumericDiffCostFunction, RandomizedCostFunctionRidders) {
+  internal::scoped_ptr<CostFunction> cost_function;
+  NumericDiffOptions options;
+  // Larger initial step size is chosen to produce robust results in the
+  // presence of random noise.
+  options.ridders_relative_initial_step_size = 10.0;
+
+  cost_function.reset(
+      new NumericDiffCostFunction<RandomizedCostFunction,
+                                  RIDDERS,
+                                  1,  /* number of residuals */
+                                  1   /* size of x1 */>(
+             new RandomizedCostFunction(kNoiseFactor, kRandomSeed),
+             TAKE_OWNERSHIP, 1, options));
+  RandomizedFunctor functor (kNoiseFactor, kRandomSeed);
+  functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function);
+}
+
+
 }  // namespace internal
 }  // namespace ceres
diff --git a/internal/ceres/numeric_diff_test_utils.cc b/internal/ceres/numeric_diff_test_utils.cc
index 5595449..bab186c 100644
--- a/internal/ceres/numeric_diff_test_utils.cc
+++ b/internal/ceres/numeric_diff_test_utils.cc
@@ -27,6 +27,7 @@
 // POSSIBILITY OF SUCH DAMAGE.
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
+//         tbennun@gmail.com (Tal Ben-Nun)
 
 #include "ceres/numeric_diff_test_utils.h"
 
@@ -56,7 +57,7 @@
 
 void EasyFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
     const CostFunction& cost_function,
-    NumericDiffMethod method) const {
+    NumericDiffMethodType method) const {
   // The x1[0] is made deliberately small to test the performance near
   // zero.
   double x1[] = { 1e-64, 2.0, 3.0, 4.0, 5.0 };
@@ -80,7 +81,21 @@
   EXPECT_EQ(expected_residuals[1], residuals[1]);
   EXPECT_EQ(expected_residuals[2], residuals[2]);
 
-  const double tolerance = (method == CENTRAL)? 3e-9 : 2e-5;
+  double tolerance = 0.0;
+  switch (method) {
+    default:
+    case CENTRAL:
+      tolerance = 3e-9;
+      break;
+
+    case FORWARD:
+      tolerance = 2e-5;
+      break;
+
+    case RIDDERS:
+      tolerance = 1e-13;
+      break;
+  }
 
   for (int i = 0; i < 5; ++i) {
     ExpectClose(x2[i],                    dydx1[5 * 0 + i], tolerance);  // y1
@@ -106,7 +121,7 @@
 
 void TranscendentalFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
     const CostFunction& cost_function,
-    NumericDiffMethod method) const {
+    NumericDiffMethodType method) const {
   struct {
     double x1[5];
     double x2[5];
@@ -150,7 +165,21 @@
       x1x2 += x1[i] * x2[i];
     }
 
-    const double tolerance = (method == CENTRAL)? 2e-7 : 2e-5;
+    double tolerance = 0.0;
+    switch (method) {
+      default:
+      case CENTRAL:
+        tolerance = 2e-7;
+        break;
+
+      case FORWARD:
+        tolerance = 2e-5;
+        break;
+
+      case RIDDERS:
+        tolerance = 3e-12;
+        break;
+    }
 
     for (int i = 0; i < 5; ++i) {
       ExpectClose( x2[i] * cos(x1x2),              dydx1[5 * 0 + i], tolerance);
@@ -161,5 +190,81 @@
   }
 }
 
+bool ExponentialFunctor::operator()(const double* x1,
+                                    double* residuals) const {
+  residuals[0] = exp(x1[0]);
+  return true;
+}
+
+void ExponentialFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
+    const CostFunction& cost_function) const {
+  // Evaluating the functor at specific points for testing.
+  double kTests[] = { 1.0, 2.0, 3.0, 4.0, 5.0 };
+
+  // Minimal tolerance w.r.t. the cost function and the tests.
+  const double kTolerance = 2e-14;
+
+  for (int k = 0; k < CERES_ARRAYSIZE(kTests); ++k) {
+    double *parameters[] = { &kTests[k] };
+    double dydx;
+    double *jacobians[1] = { &dydx };
+    double residual;
+
+    ASSERT_TRUE(cost_function.Evaluate(&parameters[0],
+                                       &residual,
+                                       &jacobians[0]));
+
+
+    double expected_result = exp(kTests[k]);
+
+    // Expect residual to be close to exp(x).
+    ExpectClose(residual, expected_result, kTolerance);
+
+    // Check evaluated differences. dydx should also be close to exp(x).
+    ExpectClose(dydx, expected_result, kTolerance);
+  }
+}
+
+bool RandomizedFunctor::operator()(const double* x1,
+                                   double* residuals) const {
+  double random_value = static_cast<double>(rand()) /
+      static_cast<double>(RAND_MAX);
+
+  // Normalize noise to [-factor, factor].
+  random_value *= 2.0;
+  random_value -= 1.0;
+  random_value *= noise_factor_;
+
+  residuals[0] = x1[0] * x1[0] + random_value;
+  return true;
+}
+
+void RandomizedFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
+    const CostFunction& cost_function) const {
+  double kTests[] = { 0.0, 1.0, 3.0, 4.0, 50.0 };
+
+  const double kTolerance = 2e-4;
+
+  // Initialize random number generator with given seed.
+  srand(random_seed_);
+
+  for (int k = 0; k < CERES_ARRAYSIZE(kTests); ++k) {
+    double *parameters[] = { &kTests[k] };
+    double dydx;
+    double *jacobians[1] = { &dydx };
+    double residual;
+
+    ASSERT_TRUE(cost_function.Evaluate(&parameters[0],
+                                       &residual,
+                                       &jacobians[0]));
+
+    // Expect residual to be close to x^2 w.r.t. noise factor.
+    ExpectClose(residual, kTests[k] * kTests[k], noise_factor_);
+
+    // Check evaluated differences. (dy/dx = ~2x)
+    ExpectClose(dydx, 2 * kTests[k], kTolerance);
+  }
+}
+
 }  // namespace internal
 }  // namespace ceres
diff --git a/internal/ceres/numeric_diff_test_utils.h b/internal/ceres/numeric_diff_test_utils.h
index 3ae83fa..2a551d3 100644
--- a/internal/ceres/numeric_diff_test_utils.h
+++ b/internal/ceres/numeric_diff_test_utils.h
@@ -38,6 +38,12 @@
 namespace ceres {
 namespace internal {
 
+// Noise factor for randomized cost function.
+static const double kNoiseFactor = 0.01;
+
+// Default random seed for randomized cost function.
+static const unsigned int kRandomSeed = 1234;
+
 // y1 = x1'x2      -> dy1/dx1 = x2,               dy1/dx2 = x1
 // y2 = (x1'x2)^2  -> dy2/dx1 = 2 * x2 * (x1'x2), dy2/dx2 = 2 * x1 * (x1'x2)
 // y3 = x2'x2      -> dy3/dx1 = 0,                dy3/dx2 = 2 * x2
@@ -46,7 +52,7 @@
   bool operator()(const double* x1, const double* x2, double* residuals) const;
   void ExpectCostFunctionEvaluationIsNearlyCorrect(
       const CostFunction& cost_function,
-      NumericDiffMethod method) const;
+      NumericDiffMethodType method) const;
 };
 
 class EasyCostFunction : public SizedCostFunction<3, 5, 5> {
@@ -71,7 +77,7 @@
   bool operator()(const double* x1, const double* x2, double* residuals) const;
   void ExpectCostFunctionEvaluationIsNearlyCorrect(
       const CostFunction& cost_function,
-      NumericDiffMethod method) const;
+      NumericDiffMethodType method) const;
 };
 
 class TranscendentalCostFunction : public SizedCostFunction<2, 5, 5> {
@@ -85,6 +91,61 @@
   TranscendentalFunctor functor_;
 };
 
+// y = exp(x), dy/dx = exp(x)
+class ExponentialFunctor {
+ public:
+  bool operator()(const double* x1, double* residuals) const;
+  void ExpectCostFunctionEvaluationIsNearlyCorrect(
+      const CostFunction& cost_function) const;
+};
+
+class ExponentialCostFunction : public SizedCostFunction<1, 1> {
+ public:
+  virtual bool Evaluate(double const* const* parameters,
+                        double* residuals,
+                        double** /* not used */) const {
+    return functor_(parameters[0], residuals);
+  }
+
+ private:
+  ExponentialFunctor functor_;
+};
+
+// Test adaptive numeric differentiation by synthetically adding random noise
+// to a functor.
+// y = x^2 + [random noise], dy/dx ~ 2x
+class RandomizedFunctor {
+ public:
+  RandomizedFunctor(double noise_factor, unsigned int random_seed)
+      : noise_factor_(noise_factor), random_seed_(random_seed) {
+  }
+
+  bool operator()(const double* x1, double* residuals) const;
+  void ExpectCostFunctionEvaluationIsNearlyCorrect(
+      const CostFunction& cost_function) const;
+
+ private:
+  double noise_factor_;
+  unsigned int random_seed_;
+};
+
+class RandomizedCostFunction : public SizedCostFunction<1, 1> {
+ public:
+  RandomizedCostFunction(double noise_factor, unsigned int random_seed)
+      : functor_(noise_factor, random_seed) {
+  }
+
+  virtual bool Evaluate(double const* const* parameters,
+                        double* residuals,
+                        double** /* not used */) const {
+    return functor_(parameters[0], residuals);
+  }
+
+ private:
+  RandomizedFunctor functor_;
+};
+
+
 }  // namespace internal
 }  // namespace ceres
 
diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc
index 275d7f8..f86fb78 100644
--- a/internal/ceres/types.cc
+++ b/internal/ceres/types.cc
@@ -284,6 +284,27 @@
   return false;
 }
 
+const char* NumericDiffMethodTypeToString(
+    NumericDiffMethodType type) {
+  switch (type) {
+    CASESTR(CENTRAL);
+    CASESTR(FORWARD);
+    CASESTR(RIDDERS);
+    default:
+      return "UNKNOWN";
+  }
+}
+
+bool StringToNumericDiffMethodType(
+    string value,
+    NumericDiffMethodType* type) {
+  UpperCase(&value);
+  STRENUM(CENTRAL);
+  STRENUM(FORWARD);
+  STRENUM(RIDDERS);
+  return false;
+}
+
 const char* VisibilityClusteringTypeToString(
     VisibilityClusteringType type) {
   switch (type) {