Add NumericDiffFirstOrderFunction This has been a long requested feature so that users can minimize functions using numeric differentiation. As part of this, I have also redone rosenbrock.cc, which now has three variants. rosenbrock.cc now uses automatic differentiation. rosenbrock_numeric_diff.cc uses numeric differentiation. rosenbrock_analytic_diff.cc uses analytic derivatives. This is analogus to how the helloworld example code is structured. The tutorial for GradientProblemSolver has also been updated to reflect this. https://github.com/ceres-solver/ceres-solver/issues/691 Change-Id: Ib0fb9e35127fe4c8299d4793bea3558722c70dd7
diff --git a/docs/source/gradient_tutorial.rst b/docs/source/gradient_tutorial.rst index 3fef6b6..aea407a 100644 --- a/docs/source/gradient_tutorial.rst +++ b/docs/source/gradient_tutorial.rst
@@ -8,46 +8,46 @@ General Unconstrained Minimization ================================== -While much of Ceres Solver is devoted to solving non-linear least -squares problems, internally it contains a solver that can solve -general unconstrained optimization problems using just their objective -function value and gradients. The ``GradientProblem`` and -``GradientProblemSolver`` objects give the user access to this solver. - -So without much further ado, let us look at how one goes about using -them. +Ceres Solver besides being able to solve non-linear least squares +problem can also solve general unconstrained problems using just their +objective function value and gradients. In this chapter we will see +how to do this. Rosenbrock's Function ===================== -We consider the minimization of the famous `Rosenbrock's function +Consider minimizing the famous `Rosenbrock's function <http://en.wikipedia.org/wiki/Rosenbrock_function>`_ [#f1]_. -We begin by defining an instance of the ``FirstOrderFunction`` -interface. This is the object that is responsible for computing the -objective function value and the gradient (if required). This is the -analog of the :class:`CostFunction` when defining non-linear least -squares problems in Ceres. +The simplest way to minimize is to define a templated functor to +evaluate the objective value of this function and then use Ceres +Solver's automatic differentiation to compute its derivatives. + +We begin by defining a templated functor and then using +``AutoDiffFirstOrderFunction`` to construct an instance of the +``FirstOrderFunction`` interface. This is the object that is +responsible for computing the objective function value and the +gradient (if required). This is the analog of the +:class:`CostFunction` when defining non-linear least squares problems +in Ceres. .. code:: - class Rosenbrock : public ceres::FirstOrderFunction { - public: - virtual bool Evaluate(const double* parameters, - double* cost, - double* gradient) const { - const double x = parameters[0]; - const double y = parameters[1]; - + // f(x,y) = (1-x)^2 + 100(y - x^2)^2; + struct Rosenbrock { + template <typename T> + bool operator()(const T* parameters, T* cost) const { + const T x = parameters[0]; + const T y = parameters[1]; cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x); - if (gradient != nullptr) { - gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x; - gradient[1] = 200.0 * (y - x * x); - } return true; } - virtual int NumParameters() const { return 2; } + static ceres::FirstOrderFunction* Create() { + constexpr int kNumParameters = 2; + return new ceres::AutoDiffFirstOrderFunction<Rosenbrock, kNumParameters>( + new Rosenbrock); + } }; @@ -58,7 +58,7 @@ double parameters[2] = {-1.2, 1.0}; - ceres::GradientProblem problem(new Rosenbrock()); + ceres::GradientProblem problem(Rosenbrock::Create()); ceres::GradientProblemSolver::Options options; options.minimizer_progress_to_stdout = true; @@ -74,43 +74,43 @@ .. code-block:: bash - 0: f: 2.420000e+01 d: 0.00e+00 g: 2.16e+02 h: 0.00e+00 s: 0.00e+00 e: 0 it: 2.00e-05 tt: 2.00e-05 - 1: f: 4.280493e+00 d: 1.99e+01 g: 1.52e+01 h: 2.01e-01 s: 8.62e-04 e: 2 it: 7.32e-05 tt: 2.19e-04 - 2: f: 3.571154e+00 d: 7.09e-01 g: 1.35e+01 h: 3.78e-01 s: 1.34e-01 e: 3 it: 2.50e-05 tt: 2.68e-04 - 3: f: 3.440869e+00 d: 1.30e-01 g: 1.73e+01 h: 1.36e-01 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 2.92e-04 - 4: f: 3.213597e+00 d: 2.27e-01 g: 1.55e+01 h: 1.06e-01 s: 4.59e-01 e: 1 it: 2.86e-06 tt: 3.14e-04 - 5: f: 2.839723e+00 d: 3.74e-01 g: 1.05e+01 h: 1.34e-01 s: 5.24e-01 e: 1 it: 2.86e-06 tt: 3.36e-04 - 6: f: 2.448490e+00 d: 3.91e-01 g: 1.29e+01 h: 3.04e-01 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 3.58e-04 - 7: f: 1.943019e+00 d: 5.05e-01 g: 4.00e+00 h: 8.81e-02 s: 7.43e-01 e: 1 it: 4.05e-06 tt: 3.79e-04 - 8: f: 1.731469e+00 d: 2.12e-01 g: 7.36e+00 h: 1.71e-01 s: 4.60e-01 e: 2 it: 9.06e-06 tt: 4.06e-04 - 9: f: 1.503267e+00 d: 2.28e-01 g: 6.47e+00 h: 8.66e-02 s: 1.00e+00 e: 1 it: 3.81e-06 tt: 4.33e-04 - 10: f: 1.228331e+00 d: 2.75e-01 g: 2.00e+00 h: 7.70e-02 s: 7.90e-01 e: 1 it: 3.81e-06 tt: 4.54e-04 - 11: f: 1.016523e+00 d: 2.12e-01 g: 5.15e+00 h: 1.39e-01 s: 3.76e-01 e: 2 it: 1.00e-05 tt: 4.82e-04 - 12: f: 9.145773e-01 d: 1.02e-01 g: 6.74e+00 h: 7.98e-02 s: 1.00e+00 e: 1 it: 3.10e-06 tt: 5.03e-04 - 13: f: 7.508302e-01 d: 1.64e-01 g: 3.88e+00 h: 5.76e-02 s: 4.93e-01 e: 1 it: 2.86e-06 tt: 5.25e-04 - 14: f: 5.832378e-01 d: 1.68e-01 g: 5.56e+00 h: 1.42e-01 s: 1.00e+00 e: 1 it: 3.81e-06 tt: 5.47e-04 - 15: f: 3.969581e-01 d: 1.86e-01 g: 1.64e+00 h: 1.17e-01 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 5.68e-04 - 16: f: 3.171557e-01 d: 7.98e-02 g: 3.84e+00 h: 1.18e-01 s: 3.97e-01 e: 2 it: 9.06e-06 tt: 5.94e-04 - 17: f: 2.641257e-01 d: 5.30e-02 g: 3.27e+00 h: 6.14e-02 s: 1.00e+00 e: 1 it: 3.10e-06 tt: 6.16e-04 - 18: f: 1.909730e-01 d: 7.32e-02 g: 5.29e-01 h: 8.55e-02 s: 6.82e-01 e: 1 it: 4.05e-06 tt: 6.42e-04 - 19: f: 1.472012e-01 d: 4.38e-02 g: 3.11e+00 h: 1.20e-01 s: 3.47e-01 e: 2 it: 1.00e-05 tt: 6.69e-04 - 20: f: 1.093558e-01 d: 3.78e-02 g: 2.97e+00 h: 8.43e-02 s: 1.00e+00 e: 1 it: 3.81e-06 tt: 6.91e-04 - 21: f: 6.710346e-02 d: 4.23e-02 g: 1.42e+00 h: 9.64e-02 s: 8.85e-01 e: 1 it: 3.81e-06 tt: 7.12e-04 - 22: f: 3.993377e-02 d: 2.72e-02 g: 2.30e+00 h: 1.29e-01 s: 4.63e-01 e: 2 it: 9.06e-06 tt: 7.39e-04 - 23: f: 2.911794e-02 d: 1.08e-02 g: 2.55e+00 h: 6.55e-02 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 7.62e-04 - 24: f: 1.457683e-02 d: 1.45e-02 g: 2.77e-01 h: 6.37e-02 s: 6.14e-01 e: 1 it: 3.81e-06 tt: 7.84e-04 - 25: f: 8.577515e-03 d: 6.00e-03 g: 2.86e+00 h: 1.40e-01 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 8.05e-04 - 26: f: 3.486574e-03 d: 5.09e-03 g: 1.76e-01 h: 1.23e-02 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 8.27e-04 - 27: f: 1.257570e-03 d: 2.23e-03 g: 1.39e-01 h: 5.08e-02 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 8.48e-04 - 28: f: 2.783568e-04 d: 9.79e-04 g: 6.20e-01 h: 6.47e-02 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 8.69e-04 - 29: f: 2.533399e-05 d: 2.53e-04 g: 1.68e-02 h: 1.98e-03 s: 1.00e+00 e: 1 it: 3.81e-06 tt: 8.91e-04 - 30: f: 7.591572e-07 d: 2.46e-05 g: 5.40e-03 h: 9.27e-03 s: 1.00e+00 e: 1 it: 3.81e-06 tt: 9.12e-04 - 31: f: 1.902460e-09 d: 7.57e-07 g: 1.62e-03 h: 1.89e-03 s: 1.00e+00 e: 1 it: 2.86e-06 tt: 9.33e-04 - 32: f: 1.003030e-12 d: 1.90e-09 g: 3.50e-05 h: 3.52e-05 s: 1.00e+00 e: 1 it: 3.10e-06 tt: 9.54e-04 - 33: f: 4.835994e-17 d: 1.00e-12 g: 1.05e-07 h: 1.13e-06 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 9.81e-04 - 34: f: 1.885250e-22 d: 4.84e-17 g: 2.69e-10 h: 1.45e-08 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 1.00e-03 + 0: f: 2.420000e+01 d: 0.00e+00 g: 2.16e+02 h: 0.00e+00 s: 0.00e+00 e: 0 it: 2.29e-05 tt: 2.29e-05 + 1: f: 4.280493e+00 d: 1.99e+01 g: 1.52e+01 h: 2.01e-01 s: 8.62e-04 e: 2 it: 8.39e-05 tt: 1.62e-04 + 2: f: 3.571154e+00 d: 7.09e-01 g: 1.35e+01 h: 3.78e-01 s: 1.34e-01 e: 3 it: 2.22e-05 tt: 1.91e-04 + 3: f: 3.440869e+00 d: 1.30e-01 g: 1.73e+01 h: 1.36e-01 s: 1.00e+00 e: 1 it: 5.01e-06 tt: 2.01e-04 + 4: f: 3.213597e+00 d: 2.27e-01 g: 1.55e+01 h: 1.06e-01 s: 4.59e-01 e: 1 it: 3.81e-06 tt: 2.10e-04 + 5: f: 2.839723e+00 d: 3.74e-01 g: 1.05e+01 h: 1.34e-01 s: 5.24e-01 e: 1 it: 4.05e-06 tt: 2.19e-04 + 6: f: 2.448490e+00 d: 3.91e-01 g: 1.29e+01 h: 3.04e-01 s: 1.00e+00 e: 1 it: 5.01e-06 tt: 2.28e-04 + 7: f: 1.943019e+00 d: 5.05e-01 g: 4.00e+00 h: 8.81e-02 s: 7.43e-01 e: 1 it: 4.05e-06 tt: 2.36e-04 + 8: f: 1.731469e+00 d: 2.12e-01 g: 7.36e+00 h: 1.71e-01 s: 4.60e-01 e: 2 it: 1.22e-05 tt: 2.52e-04 + 9: f: 1.503267e+00 d: 2.28e-01 g: 6.47e+00 h: 8.66e-02 s: 1.00e+00 e: 1 it: 5.96e-06 tt: 2.66e-04 + 10: f: 1.228331e+00 d: 2.75e-01 g: 2.00e+00 h: 7.70e-02 s: 7.90e-01 e: 1 it: 4.05e-06 tt: 2.75e-04 + 11: f: 1.016523e+00 d: 2.12e-01 g: 5.15e+00 h: 1.39e-01 s: 3.76e-01 e: 2 it: 9.06e-06 tt: 2.88e-04 + 12: f: 9.145773e-01 d: 1.02e-01 g: 6.74e+00 h: 7.98e-02 s: 1.00e+00 e: 1 it: 5.01e-06 tt: 2.97e-04 + 13: f: 7.508302e-01 d: 1.64e-01 g: 3.88e+00 h: 5.76e-02 s: 4.93e-01 e: 1 it: 5.01e-06 tt: 3.05e-04 + 14: f: 5.832378e-01 d: 1.68e-01 g: 5.56e+00 h: 1.42e-01 s: 1.00e+00 e: 1 it: 4.77e-06 tt: 3.13e-04 + 15: f: 3.969581e-01 d: 1.86e-01 g: 1.64e+00 h: 1.17e-01 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 3.20e-04 + 16: f: 3.171557e-01 d: 7.98e-02 g: 3.84e+00 h: 1.18e-01 s: 3.97e-01 e: 2 it: 8.82e-06 tt: 3.33e-04 + 17: f: 2.641257e-01 d: 5.30e-02 g: 3.27e+00 h: 6.14e-02 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 3.42e-04 + 18: f: 1.909730e-01 d: 7.32e-02 g: 5.29e-01 h: 8.55e-02 s: 6.82e-01 e: 1 it: 1.00e-05 tt: 4.64e-04 + 19: f: 1.472012e-01 d: 4.38e-02 g: 3.11e+00 h: 1.20e-01 s: 3.47e-01 e: 2 it: 1.29e-05 tt: 4.87e-04 + 20: f: 1.093558e-01 d: 3.78e-02 g: 2.97e+00 h: 8.43e-02 s: 1.00e+00 e: 1 it: 5.01e-06 tt: 4.97e-04 + 21: f: 6.710346e-02 d: 4.23e-02 g: 1.42e+00 h: 9.64e-02 s: 8.85e-01 e: 1 it: 4.05e-06 tt: 5.06e-04 + 22: f: 3.993377e-02 d: 2.72e-02 g: 2.30e+00 h: 1.29e-01 s: 4.63e-01 e: 2 it: 1.00e-05 tt: 5.25e-04 + 23: f: 2.911794e-02 d: 1.08e-02 g: 2.55e+00 h: 6.55e-02 s: 1.00e+00 e: 1 it: 5.01e-06 tt: 5.34e-04 + 24: f: 1.457683e-02 d: 1.45e-02 g: 2.77e-01 h: 6.37e-02 s: 6.14e-01 e: 1 it: 4.05e-06 tt: 5.42e-04 + 25: f: 8.577515e-03 d: 6.00e-03 g: 2.86e+00 h: 1.40e-01 s: 1.00e+00 e: 1 it: 3.81e-06 tt: 5.49e-04 + 26: f: 3.486574e-03 d: 5.09e-03 g: 1.76e-01 h: 1.23e-02 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 5.57e-04 + 27: f: 1.257570e-03 d: 2.23e-03 g: 1.39e-01 h: 5.08e-02 s: 1.00e+00 e: 1 it: 3.81e-06 tt: 5.65e-04 + 28: f: 2.783568e-04 d: 9.79e-04 g: 6.20e-01 h: 6.47e-02 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 5.73e-04 + 29: f: 2.533399e-05 d: 2.53e-04 g: 1.68e-02 h: 1.98e-03 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 5.81e-04 + 30: f: 7.591572e-07 d: 2.46e-05 g: 5.40e-03 h: 9.27e-03 s: 1.00e+00 e: 1 it: 5.96e-06 tt: 6.30e-04 + 31: f: 1.902460e-09 d: 7.57e-07 g: 1.62e-03 h: 1.89e-03 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 6.39e-04 + 32: f: 1.003030e-12 d: 1.90e-09 g: 3.50e-05 h: 3.52e-05 s: 1.00e+00 e: 1 it: 3.81e-06 tt: 6.47e-04 + 33: f: 4.835994e-17 d: 1.00e-12 g: 1.05e-07 h: 1.13e-06 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 6.59e-04 + 34: f: 1.885250e-22 d: 4.84e-17 g: 2.69e-10 h: 1.45e-08 s: 1.00e+00 e: 1 it: 4.05e-06 tt: 6.67e-04 - Solver Summary (v 1.12.0-lapack-suitesparse-cxsparse-no_openmp) + Solver Summary (v 2.0.0-eigen-(3.3.9)-lapack-suitesparse-(5.8.1)-cxsparse-(3.2.0)-acceleratesparse-eigensparse-no_openmp-no_custom_blas) Parameters 2 Line search direction LBFGS (20) @@ -119,20 +119,86 @@ Cost: Initial 2.420000e+01 - Final 1.885250e-22 + Final 1.955192e-27 Change 2.420000e+01 - Minimizer iterations 35 + Minimizer iterations 36 Time (in seconds): - Cost evaluation 0.000 - Gradient evaluation 0.000 - Total 0.003 + Cost evaluation 0.000000 (0) + Gradient & cost evaluation 0.000008 (44) + Polynomial minimization 0.000067 + Total 0.000721 - Termination: CONVERGENCE (Gradient tolerance reached. Gradient max norm: 9.032775e-13 <= 1.000000e-10) + Termination: CONVERGENCE (Parameter tolerance reached. Relative step_norm: 1.890726e-11 <= 1.000000e-08.) + + Initial x: -1.2 y: 1 + Final x: 1 y: 1 + + +If you are unable to use automatic differentiation for some reason +(say beacause you need to call an external library), then you can +use numeric differentiation. In that case the functor is defined as +follows [#f2]_. + +.. code:: + + // f(x,y) = (1-x)^2 + 100(y - x^2)^2; + struct Rosenbrock { + bool operator()(const double* parameters, double* cost) const { + const double x = parameters[0]; + const double y = parameters[1]; + cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x); + return true; + } + + static ceres::FirstOrderFunction* Create() { + constexpr int kNumParameters = 2; + return new ceres::NumericDiffFirstOrderFunction<Rosenbrock, + ceres::CENTRAL, + kNumParameters>( + new Rosenbrock); + } +}; + +And finally, if you would rather compute the derivatives by hand (say +because the size of the parameter vector is too large to be +automatically differentiated). Then you should define an instance of +``FirstOrderFunction``, which is the analog of :class:`CostFunction` +for non-linear least squares problems [#f3]_. + +.. code:: + + // f(x,y) = (1-x)^2 + 100(y - x^2)^2; + class Rosenbrock final : public ceres::FirstOrderFunction { + public: + ~Rosenbrock() override {} + + bool Evaluate(const double* parameters, + double* cost, + double* gradient) const override { + const double x = parameters[0]; + const double y = parameters[1]; + + cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x); + if (gradient) { + gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x; + gradient[1] = 200.0 * (y - x * x); + } + return true; + } + + int NumParameters() const override { return 2; } + }; .. rubric:: Footnotes .. [#f1] `examples/rosenbrock.cc <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/rosenbrock.cc>`_ + +.. [#f2] `examples/rosenbrock_numeric_diff.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/rosenbrock_numeric_diff.cc>`_ + +.. [#f3] `examples/rosenbrock_analytic_diff.cc + <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/rosenbrock_analytic_diff.cc>`_
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 7f9b117..a43c8aa 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt
@@ -50,6 +50,12 @@ add_executable(rosenbrock rosenbrock.cc) target_link_libraries(rosenbrock Ceres::ceres) +add_executable(rosenbrock_analytic_diff rosenbrock_analytic_diff.cc) +target_link_libraries(rosenbrock_analytic_diff Ceres::ceres) + +add_executable(rosenbrock_numeric_diff rosenbrock_numeric_diff.cc) +target_link_libraries(rosenbrock_numeric_diff Ceres::ceres) + add_executable(curve_fitting_c curve_fitting.c) target_link_libraries(curve_fitting_c Ceres::ceres) # Force CMake to link curve_fitting_c using the C linker.
diff --git a/examples/rosenbrock.cc b/examples/rosenbrock.cc index 1b9aef6..8243457 100644 --- a/examples/rosenbrock.cc +++ b/examples/rosenbrock.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2021 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -32,25 +32,20 @@ #include "glog/logging.h" // f(x,y) = (1-x)^2 + 100(y - x^2)^2; -class Rosenbrock : public ceres::FirstOrderFunction { - public: - virtual ~Rosenbrock() {} - - virtual bool Evaluate(const double* parameters, - double* cost, - double* gradient) const { - const double x = parameters[0]; - const double y = parameters[1]; - +struct Rosenbrock { + template <typename T> + bool operator()(const T* parameters, T* cost) const { + const T x = parameters[0]; + const T y = parameters[1]; cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x); - if (gradient != NULL) { - gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x; - gradient[1] = 200.0 * (y - x * x); - } return true; } - virtual int NumParameters() const { return 2; } + static ceres::FirstOrderFunction* Create() { + constexpr int kNumParameters = 2; + return new ceres::AutoDiffFirstOrderFunction<Rosenbrock, kNumParameters>( + new Rosenbrock); + } }; int main(int argc, char** argv) { @@ -62,7 +57,7 @@ options.minimizer_progress_to_stdout = true; ceres::GradientProblemSolver::Summary summary; - ceres::GradientProblem problem(new Rosenbrock()); + ceres::GradientProblem problem(Rosenbrock::Create()); ceres::Solve(options, problem, parameters, &summary); std::cout << summary.FullReport() << "\n";
diff --git a/examples/rosenbrock_analytic_diff.cc b/examples/rosenbrock_analytic_diff.cc new file mode 100644 index 0000000..b133f5a --- /dev/null +++ b/examples/rosenbrock_analytic_diff.cc
@@ -0,0 +1,75 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2021 Google Inc. All rights reserved. +// http://ceres-solver.org/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: sameeragarwal@google.com (Sameer Agarwal) + +#include "ceres/ceres.h" +#include "glog/logging.h" + +// f(x,y) = (1-x)^2 + 100(y - x^2)^2; +class Rosenbrock final : public ceres::FirstOrderFunction { + public: + ~Rosenbrock() override {} + + bool Evaluate(const double* parameters, + double* cost, + double* gradient) const override { + const double x = parameters[0]; + const double y = parameters[1]; + + cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x); + + if (gradient) { + gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x; + gradient[1] = 200.0 * (y - x * x); + } + + return true; + } + + int NumParameters() const override { return 2; } +}; + +int main(int argc, char** argv) { + google::InitGoogleLogging(argv[0]); + + double parameters[2] = {-1.2, 1.0}; + + ceres::GradientProblemSolver::Options options; + options.minimizer_progress_to_stdout = true; + + ceres::GradientProblemSolver::Summary summary; + ceres::GradientProblem problem(new Rosenbrock()); + ceres::Solve(options, problem, parameters, &summary); + + std::cout << summary.FullReport() << "\n"; + std::cout << "Initial x: " << -1.2 << " y: " << 1.0 << "\n"; + std::cout << "Final x: " << parameters[0] << " y: " << parameters[1] + << "\n"; + return 0; +}
diff --git a/examples/rosenbrock_numeric_diff.cc b/examples/rosenbrock_numeric_diff.cc new file mode 100644 index 0000000..f2f6e8a --- /dev/null +++ b/examples/rosenbrock_numeric_diff.cc
@@ -0,0 +1,69 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2021 Google Inc. All rights reserved. +// http://ceres-solver.org/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: sameeragarwal@google.com (Sameer Agarwal) + +#include "ceres/ceres.h" +#include "glog/logging.h" + +// f(x,y) = (1-x)^2 + 100(y - x^2)^2; +struct Rosenbrock { + bool operator()(const double* parameters, double* cost) const { + const double x = parameters[0]; + const double y = parameters[1]; + cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x); + return true; + } + + static ceres::FirstOrderFunction* Create() { + constexpr int kNumParameters = 2; + return new ceres::NumericDiffFirstOrderFunction<Rosenbrock, + ceres::CENTRAL, + kNumParameters>( + new Rosenbrock); + } +}; + +int main(int argc, char** argv) { + google::InitGoogleLogging(argv[0]); + + double parameters[2] = {-1.2, 1.0}; + + ceres::GradientProblemSolver::Options options; + options.minimizer_progress_to_stdout = true; + + ceres::GradientProblemSolver::Summary summary; + ceres::GradientProblem problem(Rosenbrock::Create()); + ceres::Solve(options, problem, parameters, &summary); + + std::cout << summary.FullReport() << "\n"; + std::cout << "Initial x: " << -1.2 << " y: " << 1.0 << "\n"; + std::cout << "Final x: " << parameters[0] << " y: " << parameters[1] + << "\n"; + return 0; +}
diff --git a/include/ceres/ceres.h b/include/ceres/ceres.h index d249351..86a84ab 100644 --- a/include/ceres/ceres.h +++ b/include/ceres/ceres.h
@@ -35,6 +35,7 @@ #define CERES_PUBLIC_CERES_H_ #include "ceres/autodiff_cost_function.h" +#include "ceres/autodiff_first_order_function.h" #include "ceres/autodiff_local_parameterization.h" #include "ceres/conditioned_cost_function.h" #include "ceres/context.h" @@ -47,6 +48,7 @@ #include "ceres/dynamic_cost_function_to_functor.h" #include "ceres/dynamic_numeric_diff_cost_function.h" #include "ceres/evaluation_callback.h" +#include "ceres/first_order_function.h" #include "ceres/gradient_checker.h" #include "ceres/gradient_problem.h" #include "ceres/gradient_problem_solver.h" @@ -55,6 +57,7 @@ #include "ceres/local_parameterization.h" #include "ceres/loss_function.h" #include "ceres/numeric_diff_cost_function.h" +#include "ceres/numeric_diff_first_order_function.h" #include "ceres/numeric_diff_options.h" #include "ceres/ordered_groups.h" #include "ceres/problem.h"
diff --git a/include/ceres/numeric_diff_first_order_function.h b/include/ceres/numeric_diff_first_order_function.h new file mode 100644 index 0000000..8663835 --- /dev/null +++ b/include/ceres/numeric_diff_first_order_function.h
@@ -0,0 +1,162 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2019 Google Inc. All rights reserved. +// http://ceres-solver.org/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: sameeragarwal@google.com (Sameer Agarwal) + +#ifndef CERES_PUBLIC_NUMERIC_DIFF_FIRST_ORDER_FUNCTION_H_ +#define CERES_PUBLIC_NUMERIC_DIFF_FIRST_ORDER_FUNCTION_H_ + +#include <algorithm> +#include <memory> + +#include "ceres/first_order_function.h" +#include "ceres/internal/eigen.h" +#include "ceres/internal/fixed_array.h" +#include "ceres/internal/numeric_diff.h" +#include "ceres/internal/parameter_dims.h" +#include "ceres/internal/variadic_evaluate.h" +#include "ceres/numeric_diff_options.h" +#include "ceres/types.h" + + +namespace ceres { + +// Creates FirstOrderFunctions as needed by the GradientProblem +// framework, with gradients computed via numeric differentiation. For +// more information on numeric differentiation, see the wikipedia +// article at https://en.wikipedia.org/wiki/Numerical_differentiation +// +// To get an numerically differentiated cost function, you must define +// a class with an operator() (a functor) that computes the cost. +// +// The function must write the computed value in the last argument +// (the only non-const one) and return true to indicate success. +// +// For example, consider a scalar error e = x'y - a, where both x and y are +// two-dimensional column vector parameters, the prime sign indicates +// transposition, and a is a constant. +// +// To write an numerically-differentiable cost function for the above model, +// first define the object +// +// class QuadraticCostFunctor { +// public: +// explicit QuadraticCostFunctor(double a) : a_(a) {} +// bool operator()(const double* const xy, double* cost) const { +// constexpr int kInputVectorLength = 2; +// const double* const x = xy; +// const double* const y = xy + kInputVectorLength; +// *cost = x[0] * y[0] + x[1] * y[1] - a_; +// return true; +// } +// +// private: +// double a_; +// }; +// +// +// Note that in the declaration of operator() the input parameters xy +// come first, and are passed as const pointers to array of +// doubles. The output cost is the last parameter. +// +// Then given this class definition, the numerically differentiated +// first order function with central differences used for computing the +// derivative can be constructed as follows. +// +// FirstOrderFunction* function +// = new NumericDiffFirstOrderFunction<MyScalarCostFunctor, CENTRAL, 4>( +// new QuadraticCostFunctor(1.0)); ^ ^ ^ +// | | | +// Finite Differencing Scheme -+ | | +// Dimension of xy ------------------------+ +// +// +// In the instantiation above, the template parameters following +// "QuadraticCostFunctor", "CENTRAL, 4", describe the finite +// differencing scheme as "central differencing" and the functor as +// computing its cost from a 4 dimensional input. +template <typename FirstOrderFunctor, + NumericDiffMethodType method, + int kNumParameters> +class NumericDiffFirstOrderFunction : public FirstOrderFunction { + public: + NumericDiffFirstOrderFunction( + FirstOrderFunctor* functor, + Ownership ownership = TAKE_OWNERSHIP, + const NumericDiffOptions& options = NumericDiffOptions()) + : functor_(functor), ownership_(ownership), options_(options) { + static_assert(kNumParameters > 0, "kNumParameters must be positive"); + } + + ~NumericDiffFirstOrderFunction() override { + if (ownership_ != TAKE_OWNERSHIP) { + functor_.release(); + } + } + + bool Evaluate(const double* const parameters, + double* cost, + double* gradient) const override { + using ParameterDims = internal::StaticParameterDims<kNumParameters>; + constexpr int kNumResiduals = 1; + + // Get the function value (cost) at the the point to evaluate. + if (!internal::VariadicEvaluate<ParameterDims>( + *functor_, ¶meters, cost)) { + return false; + } + + if (gradient == nullptr) { + return true; + } + + // Create a copy of the parameters which will get mutated. + internal::FixedArray<double, 32> parameters_copy(kNumParameters); + std::copy_n(parameters, kNumParameters, parameters_copy.data()); + double* parameters_ptr = parameters_copy.data(); + internal::EvaluateJacobianForParameterBlocks< + ParameterDims>::template Apply<method, kNumResiduals>(functor_.get(), + cost, + options_, + kNumResiduals, + ¶meters_ptr, + &gradient); + return true; + } + + int NumParameters() const override { return kNumParameters; } + + private: + std::unique_ptr<FirstOrderFunctor> functor_; + Ownership ownership_; + NumericDiffOptions options_; +}; + +} // namespace ceres + +#endif // CERES_PUBLIC_NUMERIC_DIFF_FIRST_ORDER_FUNCTION_H_
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index 80f8bdc..18f6090 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt
@@ -486,6 +486,7 @@ ceres_test(minimizer) ceres_test(normal_prior) ceres_test(numeric_diff_cost_function) + ceres_test(numeric_diff_first_order_function) ceres_test(ordered_groups) ceres_test(parallel_for) ceres_test(parallel_utils) @@ -549,4 +550,3 @@ add_subdirectory(autodiff_benchmarks) endif (BUILD_BENCHMARKS) -
diff --git a/internal/ceres/autodiff_first_order_function_test.cc b/internal/ceres/autodiff_first_order_function_test.cc index 7db7835..4870a45 100644 --- a/internal/ceres/autodiff_first_order_function_test.cc +++ b/internal/ceres/autodiff_first_order_function_test.cc
@@ -44,7 +44,7 @@ explicit QuadraticCostFunctor(double a) : a_(a) {} template <typename T> bool operator()(const T* const x, T* cost) const { - cost[0] = x[0] * x[1] + x[2] * x[3] - T(a_); + cost[0] = x[0] * x[1] + x[2] * x[3] - a_; return true; }
diff --git a/internal/ceres/numeric_diff_first_order_function_test.cc b/internal/ceres/numeric_diff_first_order_function_test.cc new file mode 100644 index 0000000..1926caf --- /dev/null +++ b/internal/ceres/numeric_diff_first_order_function_test.cc
@@ -0,0 +1,79 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2021 Google Inc. All rights reserved. +// http://ceres-solver.org/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: sameeragarwal@google.com (Sameer Agarwal) + +#include "ceres/numeric_diff_first_order_function.h" + +#include <memory> + +#include "ceres/array_utils.h" +#include "ceres/first_order_function.h" +#include "gtest/gtest.h" + +namespace ceres { +namespace internal { + +class QuadraticCostFunctor { + public: + explicit QuadraticCostFunctor(double a) : a_(a) {} + bool operator()(const double* const x, double* cost) const { + cost[0] = x[0] * x[1] + x[2] * x[3] - a_; + return true; + } + + private: + double a_; +}; + +TEST(NumericDiffFirstOrderFunction, BilinearDifferentiationTest) { + std::unique_ptr<FirstOrderFunction> function( + new NumericDiffFirstOrderFunction<QuadraticCostFunctor, CENTRAL, 4>( + new QuadraticCostFunctor(1.0))); + + double parameters[4] = {1.0, 2.0, 3.0, 4.0}; + double gradient[4]; + double cost; + + function->Evaluate(parameters, &cost, nullptr); + EXPECT_EQ(cost, 13.0); + + cost = -1.0; + function->Evaluate(parameters, &cost, gradient); + + EXPECT_EQ(cost, 13.0); + + const double kTolerance = 1e-9; + EXPECT_NEAR(gradient[0], parameters[1], kTolerance); + EXPECT_NEAR(gradient[1], parameters[0], kTolerance); + EXPECT_NEAR(gradient[2], parameters[3], kTolerance); + EXPECT_NEAR(gradient[3], parameters[2], kTolerance); +} + +} // namespace internal +} // namespace ceres