Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 1 | .. highlight:: c++ |
| 2 | |
| 3 | .. default-domain:: cpp |
| 4 | |
| 5 | .. _chapter-tutorial: |
| 6 | |
| 7 | ======== |
| 8 | Tutorial |
| 9 | ======== |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 10 | |
Carlos Hernandez | 79397c2 | 2014-08-07 17:51:38 -0700 | [diff] [blame] | 11 | Ceres solves robustified non-linear bounds constrained least squares |
| 12 | problems of the form |
| 13 | |
| 14 | .. math:: :label: ceresproblem |
| 15 | |
| 16 | \min_{\mathbf{x}} &\quad \frac{1}{2}\sum_{i} \rho_i\left(\left\|f_i\left(x_{i_1}, ... ,x_{i_k}\right)\right\|^2\right) \\ |
| 17 | \text{s.t.} &\quad l_j \le x_j \le u_j |
| 18 | |
| 19 | Problems of this form comes up in a broad range of areas across |
| 20 | science and engineering - from `fitting curves`_ in statistics, to |
| 21 | constructing `3D models from photographs`_ in computer vision. |
| 22 | |
| 23 | .. _fitting curves: http://en.wikipedia.org/wiki/Nonlinear_regression |
| 24 | .. _3D models from photographs: http://en.wikipedia.org/wiki/Bundle_adjustment |
| 25 | |
| 26 | In this chapter we will learn how to solve :eq:`ceresproblem` using |
| 27 | Ceres Solver. Full working code for all the examples described in this |
| 28 | chapter and more can be found in the `examples |
| 29 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/>`_ |
| 30 | directory. |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 31 | |
| 32 | The expression |
| 33 | :math:`\rho_i\left(\left\|f_i\left(x_{i_1},...,x_{i_k}\right)\right\|^2\right)` |
| 34 | is known as a ``ResidualBlock``, where :math:`f_i(\cdot)` is a |
| 35 | :class:`CostFunction` that depends on the parameter blocks |
| 36 | :math:`\left[x_{i_1},... , x_{i_k}\right]`. In most optimization |
| 37 | problems small groups of scalars occur together. For example the three |
| 38 | components of a translation vector and the four components of the |
| 39 | quaternion that define the pose of a camera. We refer to such a group |
| 40 | of small scalars as a ``ParameterBlock``. Of course a |
Carlos Hernandez | 79397c2 | 2014-08-07 17:51:38 -0700 | [diff] [blame] | 41 | ``ParameterBlock`` can just be a single parameter. :math:`l_j` and |
| 42 | :math:`u_j` are bounds on the parameter block :math:`x_j`. |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 43 | |
| 44 | :math:`\rho_i` is a :class:`LossFunction`. A :class:`LossFunction` is |
| 45 | a scalar function that is used to reduce the influence of outliers on |
Carlos Hernandez | 79397c2 | 2014-08-07 17:51:38 -0700 | [diff] [blame] | 46 | the solution of non-linear least squares problems. |
| 47 | |
| 48 | As a special case, when :math:`\rho_i(x) = x`, i.e., the identity |
| 49 | function, and :math:`l_j = -\infty` and :math:`u_j = \infty` we get |
| 50 | the more familiar `non-linear least squares problem |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 51 | <http://en.wikipedia.org/wiki/Non-linear_least_squares>`_. |
| 52 | |
Carlos Hernandez | 79397c2 | 2014-08-07 17:51:38 -0700 | [diff] [blame] | 53 | .. math:: \frac{1}{2}\sum_{i} \left\|f_i\left(x_{i_1}, ... ,x_{i_k}\right)\right\|^2. |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 54 | :label: ceresproblem2 |
| 55 | |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 56 | .. _section-hello-world: |
| 57 | |
| 58 | Hello World! |
| 59 | ============ |
| 60 | |
| 61 | To get started, consider the problem of finding the minimum of the |
| 62 | function |
| 63 | |
| 64 | .. math:: \frac{1}{2}(10 -x)^2. |
| 65 | |
| 66 | This is a trivial problem, whose minimum is located at :math:`x = 10`, |
| 67 | but it is a good place to start to illustrate the basics of solving a |
| 68 | problem with Ceres [#f1]_. |
| 69 | |
| 70 | The first step is to write a functor that will evaluate this the |
| 71 | function :math:`f(x) = 10 - x`: |
| 72 | |
| 73 | .. code-block:: c++ |
| 74 | |
| 75 | struct CostFunctor { |
| 76 | template <typename T> |
| 77 | bool operator()(const T* const x, T* residual) const { |
| 78 | residual[0] = T(10.0) - x[0]; |
| 79 | return true; |
| 80 | } |
| 81 | }; |
| 82 | |
| 83 | The important thing to note here is that ``operator()`` is a templated |
| 84 | method, which assumes that all its inputs and outputs are of some type |
Carlos Hernandez | 79397c2 | 2014-08-07 17:51:38 -0700 | [diff] [blame] | 85 | ``T``. The use of templating here allows Ceres to call |
| 86 | ``CostFunctor::operator<T>()``, with ``T=double`` when just the value |
| 87 | of the residual is needed, and with a special type ``T=Jet`` when the |
| 88 | Jacobians are needed. In :ref:`section-derivatives` we will discuss the |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 89 | various ways of supplying derivatives to Ceres in more detail. |
| 90 | |
| 91 | Once we have a way of computing the residual function, it is now time |
| 92 | to construct a non-linear least squares problem using it and have |
| 93 | Ceres solve it. |
| 94 | |
| 95 | .. code-block:: c++ |
| 96 | |
| 97 | int main(int argc, char** argv) { |
| 98 | google::InitGoogleLogging(argv[0]); |
| 99 | |
| 100 | // The variable to solve for with its initial value. |
| 101 | double initial_x = 5.0; |
| 102 | double x = initial_x; |
| 103 | |
| 104 | // Build the problem. |
| 105 | Problem problem; |
| 106 | |
| 107 | // Set up the only cost function (also known as residual). This uses |
| 108 | // auto-differentiation to obtain the derivative (jacobian). |
| 109 | CostFunction* cost_function = |
| 110 | new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor); |
| 111 | problem.AddResidualBlock(cost_function, NULL, &x); |
| 112 | |
| 113 | // Run the solver! |
| 114 | Solver::Options options; |
| 115 | options.linear_solver_type = ceres::DENSE_QR; |
| 116 | options.minimizer_progress_to_stdout = true; |
| 117 | Solver::Summary summary; |
| 118 | Solve(options, &problem, &summary); |
| 119 | |
| 120 | std::cout << summary.BriefReport() << "\n"; |
| 121 | std::cout << "x : " << initial_x |
| 122 | << " -> " << x << "\n"; |
| 123 | return 0; |
| 124 | } |
| 125 | |
| 126 | :class:`AutoDiffCostFunction` takes a ``CostFunctor`` as input, |
| 127 | automatically differentiates it and gives it a :class:`CostFunction` |
| 128 | interface. |
| 129 | |
| 130 | Compiling and running `examples/helloworld.cc |
| 131 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/helloworld.cc>`_ |
| 132 | gives us |
| 133 | |
| 134 | .. code-block:: bash |
| 135 | |
Carlos Hernandez | 79397c2 | 2014-08-07 17:51:38 -0700 | [diff] [blame] | 136 | iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time |
| 137 | 0 4.512500e+01 0.00e+00 9.50e+00 0.00e+00 0.00e+00 1.00e+04 0 5.33e-04 3.46e-03 |
| 138 | 1 4.511598e-07 4.51e+01 9.50e-04 9.50e+00 1.00e+00 3.00e+04 1 5.00e-04 4.05e-03 |
| 139 | 2 5.012552e-16 4.51e-07 3.17e-08 9.50e-04 1.00e+00 9.00e+04 1 1.60e-05 4.09e-03 |
| 140 | Ceres Solver Report: Iterations: 2, Initial cost: 4.512500e+01, Final cost: 5.012552e-16, Termination: CONVERGENCE |
| 141 | x : 0.5 -> 10 |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 142 | |
| 143 | Starting from a :math:`x=5`, the solver in two iterations goes to 10 |
| 144 | [#f2]_. The careful reader will note that this is a linear problem and |
| 145 | one linear solve should be enough to get the optimal value. The |
| 146 | default configuration of the solver is aimed at non-linear problems, |
| 147 | and for reasons of simplicity we did not change it in this example. It |
| 148 | is indeed possible to obtain the solution to this problem using Ceres |
| 149 | in one iteration. Also note that the solver did get very close to the |
| 150 | optimal function value of 0 in the very first iteration. We will |
| 151 | discuss these issues in greater detail when we talk about convergence |
| 152 | and parameter settings for Ceres. |
| 153 | |
| 154 | .. rubric:: Footnotes |
| 155 | |
| 156 | .. [#f1] `examples/helloworld.cc |
| 157 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/helloworld.cc>`_ |
| 158 | |
| 159 | .. [#f2] Actually the solver ran for three iterations, and it was |
| 160 | by looking at the value returned by the linear solver in the third |
| 161 | iteration, it observed that the update to the parameter block was too |
| 162 | small and declared convergence. Ceres only prints out the display at |
| 163 | the end of an iteration, and terminates as soon as it detects |
| 164 | convergence, which is why you only see two iterations here and not |
| 165 | three. |
| 166 | |
| 167 | .. _section-derivatives: |
| 168 | |
| 169 | |
| 170 | Derivatives |
| 171 | =========== |
| 172 | |
| 173 | Ceres Solver like most optimization packages, depends on being able to |
| 174 | evaluate the value and the derivatives of each term in the objective |
| 175 | function at arbitrary parameter values. Doing so correctly and |
| 176 | efficiently is essential to getting good results. Ceres Solver |
| 177 | provides a number of ways of doing so. You have already seen one of |
| 178 | them in action -- |
| 179 | Automatic Differentiation in `examples/helloworld.cc |
| 180 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/helloworld.cc>`_ |
| 181 | |
| 182 | We now consider the other two possibilities. Analytic and numeric |
| 183 | derivatives. |
| 184 | |
| 185 | |
| 186 | Numeric Derivatives |
| 187 | ------------------- |
| 188 | |
| 189 | In some cases, its not possible to define a templated cost functor, |
| 190 | for example when the evaluation of the residual involves a call to a |
| 191 | library function that you do not have control over. In such a |
| 192 | situation, numerical differentiation can be used. The user defines a |
| 193 | functor which computes the residual value and construct a |
| 194 | :class:`NumericDiffCostFunction` using it. e.g., for :math:`f(x) = 10 - x` |
| 195 | the corresponding functor would be |
| 196 | |
| 197 | .. code-block:: c++ |
| 198 | |
| 199 | struct NumericDiffCostFunctor { |
| 200 | bool operator()(const double* const x, double* residual) const { |
| 201 | residual[0] = 10.0 - x[0]; |
| 202 | return true; |
| 203 | } |
| 204 | }; |
| 205 | |
| 206 | Which is added to the :class:`Problem` as: |
| 207 | |
| 208 | .. code-block:: c++ |
| 209 | |
| 210 | CostFunction* cost_function = |
| 211 | new NumericDiffCostFunction<NumericDiffCostFunctor, ceres::CENTRAL, 1, 1, 1>( |
| 212 | new NumericDiffCostFunctor) |
| 213 | problem.AddResidualBlock(cost_function, NULL, &x); |
| 214 | |
| 215 | Notice the parallel from when we were using automatic differentiation |
| 216 | |
| 217 | .. code-block:: c++ |
| 218 | |
| 219 | CostFunction* cost_function = |
| 220 | new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor); |
| 221 | problem.AddResidualBlock(cost_function, NULL, &x); |
| 222 | |
| 223 | The construction looks almost identical to the one used for automatic |
| 224 | differentiation, except for an extra template parameter that indicates |
| 225 | the kind of finite differencing scheme to be used for computing the |
| 226 | numerical derivatives [#f3]_. For more details see the documentation |
| 227 | for :class:`NumericDiffCostFunction`. |
| 228 | |
| 229 | **Generally speaking we recommend automatic differentiation instead of |
| 230 | numeric differentiation. The use of C++ templates makes automatic |
| 231 | differentiation efficient, whereas numeric differentiation is |
| 232 | expensive, prone to numeric errors, and leads to slower convergence.** |
| 233 | |
| 234 | |
| 235 | Analytic Derivatives |
| 236 | -------------------- |
| 237 | |
| 238 | In some cases, using automatic differentiation is not possible. For |
| 239 | example, it may be the case that it is more efficient to compute the |
| 240 | derivatives in closed form instead of relying on the chain rule used |
| 241 | by the automatic differentiation code. |
| 242 | |
| 243 | In such cases, it is possible to supply your own residual and jacobian |
| 244 | computation code. To do this, define a subclass of |
| 245 | :class:`CostFunction` or :class:`SizedCostFunction` if you know the |
| 246 | sizes of the parameters and residuals at compile time. Here for |
| 247 | example is ``SimpleCostFunction`` that implements :math:`f(x) = 10 - |
| 248 | x`. |
| 249 | |
| 250 | .. code-block:: c++ |
| 251 | |
| 252 | class QuadraticCostFunction : public ceres::SizedCostFunction<1, 1> { |
| 253 | public: |
| 254 | virtual ~QuadraticCostFunction() {} |
| 255 | virtual bool Evaluate(double const* const* parameters, |
| 256 | double* residuals, |
| 257 | double** jacobians) const { |
| 258 | const double x = parameters[0][0]; |
| 259 | residuals[0] = 10 - x; |
| 260 | |
| 261 | // Compute the Jacobian if asked for. |
| 262 | if (jacobians != NULL && jacobians[0] != NULL) { |
| 263 | jacobians[0][0] = -1; |
| 264 | } |
| 265 | return true; |
| 266 | } |
| 267 | }; |
| 268 | |
| 269 | |
| 270 | ``SimpleCostFunction::Evaluate`` is provided with an input array of |
| 271 | ``parameters``, an output array ``residuals`` for residuals and an |
| 272 | output array ``jacobians`` for Jacobians. The ``jacobians`` array is |
| 273 | optional, ``Evaluate`` is expected to check when it is non-null, and |
| 274 | if it is the case then fill it with the values of the derivative of |
| 275 | the residual function. In this case since the residual function is |
| 276 | linear, the Jacobian is constant [#f4]_ . |
| 277 | |
| 278 | As can be seen from the above code fragments, implementing |
| 279 | :class:`CostFunction` objects is a bit tedious. We recommend that |
| 280 | unless you have a good reason to manage the jacobian computation |
| 281 | yourself, you use :class:`AutoDiffCostFunction` or |
| 282 | :class:`NumericDiffCostFunction` to construct your residual blocks. |
| 283 | |
| 284 | More About Derivatives |
| 285 | ---------------------- |
| 286 | |
| 287 | Computing derivatives is by far the most complicated part of using |
| 288 | Ceres, and depending on the circumstance the user may need more |
| 289 | sophisticated ways of computing derivatives. This section just |
| 290 | scratches the surface of how derivatives can be supplied to |
| 291 | Ceres. Once you are comfortable with using |
| 292 | :class:`NumericDiffCostFunction` and :class:`AutoDiffCostFunction` we |
| 293 | recommend taking a look at :class:`DynamicAutoDiffCostFunction`, |
| 294 | :class:`CostFunctionToFunctor`, :class:`NumericDiffFunctor` and |
| 295 | :class:`ConditionedCostFunction` for more advanced ways of |
| 296 | constructing and computing cost functions. |
| 297 | |
| 298 | .. rubric:: Footnotes |
| 299 | |
| 300 | .. [#f3] `examples/helloworld_numeric_diff.cc |
| 301 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/helloworld_numeric_diff.cc>`_. |
| 302 | |
| 303 | .. [#f4] `examples/helloworld_analytic_diff.cc |
| 304 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/helloworld_analytic_diff.cc>`_. |
| 305 | |
| 306 | |
| 307 | .. _section-powell: |
| 308 | |
| 309 | Powell's Function |
| 310 | ================= |
| 311 | |
| 312 | Consider now a slightly more complicated example -- the minimization |
| 313 | of Powell's function. Let :math:`x = \left[x_1, x_2, x_3, x_4 \right]` |
| 314 | and |
| 315 | |
| 316 | .. math:: |
| 317 | |
| 318 | \begin{align} |
| 319 | f_1(x) &= x_1 + 10x_2 \\ |
| 320 | f_2(x) &= \sqrt{5} (x_3 - x_4)\\ |
| 321 | f_3(x) &= (x_2 - 2x_3)^2\\ |
| 322 | f_4(x) &= \sqrt{10} (x_1 - x_4)^2\\ |
| 323 | F(x) &= \left[f_1(x),\ f_2(x),\ f_3(x),\ f_4(x) \right] |
| 324 | \end{align} |
| 325 | |
| 326 | |
| 327 | :math:`F(x)` is a function of four parameters, has four residuals |
| 328 | and we wish to find :math:`x` such that :math:`\frac{1}{2}\|F(x)\|^2` |
| 329 | is minimized. |
| 330 | |
| 331 | Again, the first step is to define functors that evaluate of the terms |
| 332 | in the objective functor. Here is the code for evaluating |
| 333 | :math:`f_4(x_1, x_4)`: |
| 334 | |
| 335 | .. code-block:: c++ |
| 336 | |
| 337 | struct F4 { |
| 338 | template <typename T> |
| 339 | bool operator()(const T* const x1, const T* const x4, T* residual) const { |
| 340 | residual[0] = T(sqrt(10.0)) * (x1[0] - x4[0]) * (x1[0] - x4[0]); |
| 341 | return true; |
| 342 | } |
| 343 | }; |
| 344 | |
| 345 | |
| 346 | Similarly, we can define classes ``F1``, ``F2`` and ``F4`` to evaluate |
| 347 | :math:`f_1(x_1, x_2)`, :math:`f_2(x_3, x_4)` and :math:`f_3(x_2, x_3)` |
| 348 | respectively. Using these, the problem can be constructed as follows: |
| 349 | |
| 350 | |
| 351 | .. code-block:: c++ |
| 352 | |
| 353 | double x1 = 3.0; double x2 = -1.0; double x3 = 0.0; double x4 = 1.0; |
| 354 | |
| 355 | Problem problem; |
| 356 | |
| 357 | // Add residual terms to the problem using the using the autodiff |
| 358 | // wrapper to get the derivatives automatically. |
| 359 | problem.AddResidualBlock( |
| 360 | new AutoDiffCostFunction<F1, 1, 1, 1>(new F1), NULL, &x1, &x2); |
| 361 | problem.AddResidualBlock( |
| 362 | new AutoDiffCostFunction<F2, 1, 1, 1>(new F2), NULL, &x3, &x4); |
| 363 | problem.AddResidualBlock( |
| 364 | new AutoDiffCostFunction<F3, 1, 1, 1>(new F3), NULL, &x2, &x3) |
| 365 | problem.AddResidualBlock( |
| 366 | new AutoDiffCostFunction<F4, 1, 1, 1>(new F4), NULL, &x1, &x4); |
| 367 | |
| 368 | |
| 369 | Note that each ``ResidualBlock`` only depends on the two parameters |
| 370 | that the corresponding residual object depends on and not on all four |
| 371 | parameters. Compiling and running `examples/powell.cc |
| 372 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/powell.cc>`_ |
| 373 | gives us: |
| 374 | |
| 375 | .. code-block:: bash |
| 376 | |
Carlos Hernandez | 79397c2 | 2014-08-07 17:51:38 -0700 | [diff] [blame] | 377 | Initial x1 = 3, x2 = -1, x3 = 0, x4 = 1 |
| 378 | iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time |
| 379 | 0 1.075000e+02 0.00e+00 1.55e+02 0.00e+00 0.00e+00 1.00e+04 0 4.95e-04 2.30e-03 |
| 380 | 1 5.036190e+00 1.02e+02 2.00e+01 2.16e+00 9.53e-01 3.00e+04 1 4.39e-05 2.40e-03 |
| 381 | 2 3.148168e-01 4.72e+00 2.50e+00 6.23e-01 9.37e-01 9.00e+04 1 9.06e-06 2.43e-03 |
| 382 | 3 1.967760e-02 2.95e-01 3.13e-01 3.08e-01 9.37e-01 2.70e+05 1 8.11e-06 2.45e-03 |
| 383 | 4 1.229900e-03 1.84e-02 3.91e-02 1.54e-01 9.37e-01 8.10e+05 1 6.91e-06 2.48e-03 |
| 384 | 5 7.687123e-05 1.15e-03 4.89e-03 7.69e-02 9.37e-01 2.43e+06 1 7.87e-06 2.50e-03 |
| 385 | 6 4.804625e-06 7.21e-05 6.11e-04 3.85e-02 9.37e-01 7.29e+06 1 5.96e-06 2.52e-03 |
| 386 | 7 3.003028e-07 4.50e-06 7.64e-05 1.92e-02 9.37e-01 2.19e+07 1 5.96e-06 2.55e-03 |
| 387 | 8 1.877006e-08 2.82e-07 9.54e-06 9.62e-03 9.37e-01 6.56e+07 1 5.96e-06 2.57e-03 |
| 388 | 9 1.173223e-09 1.76e-08 1.19e-06 4.81e-03 9.37e-01 1.97e+08 1 7.87e-06 2.60e-03 |
| 389 | 10 7.333425e-11 1.10e-09 1.49e-07 2.40e-03 9.37e-01 5.90e+08 1 6.20e-06 2.63e-03 |
| 390 | 11 4.584044e-12 6.88e-11 1.86e-08 1.20e-03 9.37e-01 1.77e+09 1 6.91e-06 2.65e-03 |
| 391 | 12 2.865573e-13 4.30e-12 2.33e-09 6.02e-04 9.37e-01 5.31e+09 1 5.96e-06 2.67e-03 |
| 392 | 13 1.791438e-14 2.69e-13 2.91e-10 3.01e-04 9.37e-01 1.59e+10 1 7.15e-06 2.69e-03 |
| 393 | |
| 394 | Ceres Solver v1.10.0 Solve Report |
| 395 | ---------------------------------- |
| 396 | Original Reduced |
| 397 | Parameter blocks 4 4 |
| 398 | Parameters 4 4 |
| 399 | Residual blocks 4 4 |
| 400 | Residual 4 4 |
| 401 | |
| 402 | Minimizer TRUST_REGION |
| 403 | |
| 404 | Dense linear algebra library EIGEN |
| 405 | Trust region strategy LEVENBERG_MARQUARDT |
| 406 | |
| 407 | Given Used |
| 408 | Linear solver DENSE_QR DENSE_QR |
| 409 | Threads 1 1 |
| 410 | Linear solver threads 1 1 |
| 411 | |
| 412 | Cost: |
| 413 | Initial 1.075000e+02 |
| 414 | Final 1.791438e-14 |
| 415 | Change 1.075000e+02 |
| 416 | |
| 417 | Minimizer iterations 14 |
| 418 | Successful steps 14 |
| 419 | Unsuccessful steps 0 |
| 420 | |
| 421 | Time (in seconds): |
| 422 | Preprocessor 0.002 |
| 423 | |
| 424 | Residual evaluation 0.000 |
| 425 | Jacobian evaluation 0.000 |
| 426 | Linear solver 0.000 |
| 427 | Minimizer 0.001 |
| 428 | |
| 429 | Postprocessor 0.000 |
| 430 | Total 0.005 |
| 431 | |
| 432 | Termination: CONVERGENCE (Gradient tolerance reached. Gradient max norm: 3.642190e-11 <= 1.000000e-10) |
| 433 | |
| 434 | Final x1 = 0.000292189, x2 = -2.92189e-05, x3 = 4.79511e-05, x4 = 4.79511e-05 |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 435 | |
| 436 | It is easy to see that the optimal solution to this problem is at |
| 437 | :math:`x_1=0, x_2=0, x_3=0, x_4=0` with an objective function value of |
| 438 | :math:`0`. In 10 iterations, Ceres finds a solution with an objective |
| 439 | function value of :math:`4\times 10^{-12}`. |
| 440 | |
| 441 | |
| 442 | .. rubric:: Footnotes |
| 443 | |
| 444 | .. [#f5] `examples/powell.cc |
| 445 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/powell.cc>`_. |
| 446 | |
| 447 | |
| 448 | .. _section-fitting: |
| 449 | |
| 450 | Curve Fitting |
| 451 | ============= |
| 452 | |
| 453 | The examples we have seen until now are simple optimization problems |
| 454 | with no data. The original purpose of least squares and non-linear |
| 455 | least squares analysis was fitting curves to data. It is only |
| 456 | appropriate that we now consider an example of such a problem |
| 457 | [#f6]_. It contains data generated by sampling the curve :math:`y = |
| 458 | e^{0.3x + 0.1}` and adding Gaussian noise with standard deviation |
| 459 | :math:`\sigma = 0.2`. Let us fit some data to the curve |
| 460 | |
| 461 | .. math:: y = e^{mx + c}. |
| 462 | |
| 463 | We begin by defining a templated object to evaluate the |
| 464 | residual. There will be a residual for each observation. |
| 465 | |
| 466 | .. code-block:: c++ |
| 467 | |
| 468 | struct ExponentialResidual { |
| 469 | ExponentialResidual(double x, double y) |
| 470 | : x_(x), y_(y) {} |
| 471 | |
| 472 | template <typename T> |
| 473 | bool operator()(const T* const m, const T* const c, T* residual) const { |
| 474 | residual[0] = T(y_) - exp(m[0] * T(x_) + c[0]); |
| 475 | return true; |
| 476 | } |
| 477 | |
| 478 | private: |
| 479 | // Observations for a sample. |
| 480 | const double x_; |
| 481 | const double y_; |
| 482 | }; |
| 483 | |
| 484 | Assuming the observations are in a :math:`2n` sized array called |
| 485 | ``data`` the problem construction is a simple matter of creating a |
| 486 | :class:`CostFunction` for every observation. |
| 487 | |
| 488 | |
| 489 | .. code-block:: c++ |
| 490 | |
| 491 | double m = 0.0; |
| 492 | double c = 0.0; |
| 493 | |
| 494 | Problem problem; |
| 495 | for (int i = 0; i < kNumObservations; ++i) { |
| 496 | CostFunction* cost_function = |
| 497 | new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>( |
| 498 | new ExponentialResidual(data[2 * i], data[2 * i + 1])); |
| 499 | problem.AddResidualBlock(cost_function, NULL, &m, &c); |
| 500 | } |
| 501 | |
| 502 | Compiling and running `examples/curve_fitting.cc |
| 503 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/curve_fitting.cc>`_ |
| 504 | gives us: |
| 505 | |
| 506 | .. code-block:: bash |
| 507 | |
Carlos Hernandez | 79397c2 | 2014-08-07 17:51:38 -0700 | [diff] [blame] | 508 | iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time |
| 509 | 0 1.211734e+02 0.00e+00 3.61e+02 0.00e+00 0.00e+00 1.00e+04 0 5.34e-04 2.56e-03 |
| 510 | 1 1.211734e+02 -2.21e+03 0.00e+00 7.52e-01 -1.87e+01 5.00e+03 1 4.29e-05 3.25e-03 |
| 511 | 2 1.211734e+02 -2.21e+03 0.00e+00 7.51e-01 -1.86e+01 1.25e+03 1 1.10e-05 3.28e-03 |
| 512 | 3 1.211734e+02 -2.19e+03 0.00e+00 7.48e-01 -1.85e+01 1.56e+02 1 1.41e-05 3.31e-03 |
| 513 | 4 1.211734e+02 -2.02e+03 0.00e+00 7.22e-01 -1.70e+01 9.77e+00 1 1.00e-05 3.34e-03 |
| 514 | 5 1.211734e+02 -7.34e+02 0.00e+00 5.78e-01 -6.32e+00 3.05e-01 1 1.00e-05 3.36e-03 |
| 515 | 6 3.306595e+01 8.81e+01 4.10e+02 3.18e-01 1.37e+00 9.16e-01 1 2.79e-05 3.41e-03 |
| 516 | 7 6.426770e+00 2.66e+01 1.81e+02 1.29e-01 1.10e+00 2.75e+00 1 2.10e-05 3.45e-03 |
| 517 | 8 3.344546e+00 3.08e+00 5.51e+01 3.05e-02 1.03e+00 8.24e+00 1 2.10e-05 3.48e-03 |
| 518 | 9 1.987485e+00 1.36e+00 2.33e+01 8.87e-02 9.94e-01 2.47e+01 1 2.10e-05 3.52e-03 |
| 519 | 10 1.211585e+00 7.76e-01 8.22e+00 1.05e-01 9.89e-01 7.42e+01 1 2.10e-05 3.56e-03 |
| 520 | 11 1.063265e+00 1.48e-01 1.44e+00 6.06e-02 9.97e-01 2.22e+02 1 2.60e-05 3.61e-03 |
| 521 | 12 1.056795e+00 6.47e-03 1.18e-01 1.47e-02 1.00e+00 6.67e+02 1 2.10e-05 3.64e-03 |
| 522 | 13 1.056751e+00 4.39e-05 3.79e-03 1.28e-03 1.00e+00 2.00e+03 1 2.10e-05 3.68e-03 |
| 523 | Ceres Solver Report: Iterations: 13, Initial cost: 1.211734e+02, Final cost: 1.056751e+00, Termination: CONVERGENCE |
| 524 | Initial m: 0 c: 0 |
| 525 | Final m: 0.291861 c: 0.131439 |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 526 | |
| 527 | Starting from parameter values :math:`m = 0, c=0` with an initial |
| 528 | objective function value of :math:`121.173` Ceres finds a solution |
| 529 | :math:`m= 0.291861, c = 0.131439` with an objective function value of |
| 530 | :math:`1.05675`. These values are a a bit different than the |
| 531 | parameters of the original model :math:`m=0.3, c= 0.1`, but this is |
| 532 | expected. When reconstructing a curve from noisy data, we expect to |
| 533 | see such deviations. Indeed, if you were to evaluate the objective |
| 534 | function for :math:`m=0.3, c=0.1`, the fit is worse with an objective |
| 535 | function value of :math:`1.082425`. The figure below illustrates the fit. |
| 536 | |
| 537 | .. figure:: least_squares_fit.png |
| 538 | :figwidth: 500px |
| 539 | :height: 400px |
| 540 | :align: center |
| 541 | |
| 542 | Least squares curve fitting. |
| 543 | |
| 544 | |
| 545 | .. rubric:: Footnotes |
| 546 | |
| 547 | .. [#f6] `examples/curve_fitting.cc |
| 548 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/curve_fitting.cc>`_ |
| 549 | |
| 550 | |
| 551 | Robust Curve Fitting |
| 552 | ===================== |
| 553 | |
| 554 | Now suppose the data we are given has some outliers, i.e., we have |
| 555 | some points that do not obey the noise model. If we were to use the |
| 556 | code above to fit such data, we would get a fit that looks as |
| 557 | below. Notice how the fitted curve deviates from the ground truth. |
| 558 | |
| 559 | .. figure:: non_robust_least_squares_fit.png |
| 560 | :figwidth: 500px |
| 561 | :height: 400px |
| 562 | :align: center |
| 563 | |
| 564 | To deal with outliers, a standard technique is to use a |
| 565 | :class:`LossFunction`. Loss functions, reduce the influence of |
| 566 | residual blocks with high residuals, usually the ones corresponding to |
| 567 | outliers. To associate a loss function in a residual block, we change |
| 568 | |
| 569 | .. code-block:: c++ |
| 570 | |
| 571 | problem.AddResidualBlock(cost_function, NULL , &m, &c); |
| 572 | |
| 573 | to |
| 574 | |
| 575 | .. code-block:: c++ |
| 576 | |
| 577 | problem.AddResidualBlock(cost_function, new CauchyLoss(0.5) , &m, &c); |
| 578 | |
| 579 | :class:`CauchyLoss` is one of the loss functions that ships with Ceres |
| 580 | Solver. The argument :math:`0.5` specifies the scale of the loss |
| 581 | function. As a result, we get the fit below [#f7]_. Notice how the |
| 582 | fitted curve moves back closer to the ground truth curve. |
| 583 | |
| 584 | .. figure:: robust_least_squares_fit.png |
| 585 | :figwidth: 500px |
| 586 | :height: 400px |
| 587 | :align: center |
| 588 | |
| 589 | Using :class:`LossFunction` to reduce the effect of outliers on a |
| 590 | least squares fit. |
| 591 | |
| 592 | |
| 593 | .. rubric:: Footnotes |
| 594 | |
| 595 | .. [#f7] `examples/robust_curve_fitting.cc |
| 596 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/robust_curve_fitting.cc>`_ |
| 597 | |
| 598 | |
| 599 | Bundle Adjustment |
| 600 | ================= |
| 601 | |
| 602 | One of the main reasons for writing Ceres was our need to solve large |
| 603 | scale bundle adjustment problems [HartleyZisserman]_, [Triggs]_. |
| 604 | |
| 605 | Given a set of measured image feature locations and correspondences, |
| 606 | the goal of bundle adjustment is to find 3D point positions and camera |
| 607 | parameters that minimize the reprojection error. This optimization |
| 608 | problem is usually formulated as a non-linear least squares problem, |
| 609 | where the error is the squared :math:`L_2` norm of the difference between |
| 610 | the observed feature location and the projection of the corresponding |
| 611 | 3D point on the image plane of the camera. Ceres has extensive support |
| 612 | for solving bundle adjustment problems. |
| 613 | |
| 614 | Let us solve a problem from the `BAL |
| 615 | <http://grail.cs.washington.edu/projects/bal/>`_ dataset [#f8]_. |
| 616 | |
| 617 | The first step as usual is to define a templated functor that computes |
| 618 | the reprojection error/residual. The structure of the functor is |
| 619 | similar to the ``ExponentialResidual``, in that there is an |
| 620 | instance of this object responsible for each image observation. |
| 621 | |
| 622 | Each residual in a BAL problem depends on a three dimensional point |
| 623 | and a nine parameter camera. The nine parameters defining the camera |
| 624 | are: three for rotation as a Rodriques' axis-angle vector, three |
| 625 | for translation, one for focal length and two for radial distortion. |
| 626 | The details of this camera model can be found the `Bundler homepage |
| 627 | <http://phototour.cs.washington.edu/bundler/>`_ and the `BAL homepage |
| 628 | <http://grail.cs.washington.edu/projects/bal/>`_. |
| 629 | |
| 630 | .. code-block:: c++ |
| 631 | |
| 632 | struct SnavelyReprojectionError { |
| 633 | SnavelyReprojectionError(double observed_x, double observed_y) |
| 634 | : observed_x(observed_x), observed_y(observed_y) {} |
| 635 | |
| 636 | template <typename T> |
| 637 | bool operator()(const T* const camera, |
| 638 | const T* const point, |
| 639 | T* residuals) const { |
| 640 | // camera[0,1,2] are the angle-axis rotation. |
| 641 | T p[3]; |
| 642 | ceres::AngleAxisRotatePoint(camera, point, p); |
| 643 | // camera[3,4,5] are the translation. |
| 644 | p[0] += camera[3]; p[1] += camera[4]; p[2] += camera[5]; |
| 645 | |
| 646 | // Compute the center of distortion. The sign change comes from |
| 647 | // the camera model that Noah Snavely's Bundler assumes, whereby |
| 648 | // the camera coordinate system has a negative z axis. |
| 649 | T xp = - p[0] / p[2]; |
| 650 | T yp = - p[1] / p[2]; |
| 651 | |
| 652 | // Apply second and fourth order radial distortion. |
| 653 | const T& l1 = camera[7]; |
| 654 | const T& l2 = camera[8]; |
| 655 | T r2 = xp*xp + yp*yp; |
| 656 | T distortion = T(1.0) + r2 * (l1 + l2 * r2); |
| 657 | |
| 658 | // Compute final projected point position. |
| 659 | const T& focal = camera[6]; |
| 660 | T predicted_x = focal * distortion * xp; |
| 661 | T predicted_y = focal * distortion * yp; |
| 662 | |
| 663 | // The error is the difference between the predicted and observed position. |
| 664 | residuals[0] = predicted_x - T(observed_x); |
| 665 | residuals[1] = predicted_y - T(observed_y); |
| 666 | return true; |
| 667 | } |
| 668 | |
| 669 | // Factory to hide the construction of the CostFunction object from |
| 670 | // the client code. |
| 671 | static ceres::CostFunction* Create(const double observed_x, |
| 672 | const double observed_y) { |
| 673 | return (new ceres::AutoDiffCostFunction<SnavelyReprojectionError, 2, 9, 3>( |
| 674 | new SnavelyReprojectionError(observed_x, observed_y))); |
| 675 | } |
| 676 | |
| 677 | double observed_x; |
| 678 | double observed_y; |
| 679 | }; |
| 680 | |
| 681 | |
| 682 | Note that unlike the examples before, this is a non-trivial function |
| 683 | and computing its analytic Jacobian is a bit of a pain. Automatic |
| 684 | differentiation makes life much simpler. The function |
| 685 | :func:`AngleAxisRotatePoint` and other functions for manipulating |
| 686 | rotations can be found in ``include/ceres/rotation.h``. |
| 687 | |
| 688 | Given this functor, the bundle adjustment problem can be constructed |
| 689 | as follows: |
| 690 | |
| 691 | .. code-block:: c++ |
| 692 | |
| 693 | ceres::Problem problem; |
| 694 | for (int i = 0; i < bal_problem.num_observations(); ++i) { |
| 695 | ceres::CostFunction* cost_function = |
Carlos Hernandez | 79397c2 | 2014-08-07 17:51:38 -0700 | [diff] [blame] | 696 | SnavelyReprojectionError::Create( |
| 697 | bal_problem.observations()[2 * i + 0], |
| 698 | bal_problem.observations()[2 * i + 1]); |
Sascha Haeberling | 1d2624a | 2013-07-23 19:00:21 -0700 | [diff] [blame] | 699 | problem.AddResidualBlock(cost_function, |
| 700 | NULL /* squared loss */, |
| 701 | bal_problem.mutable_camera_for_observation(i), |
| 702 | bal_problem.mutable_point_for_observation(i)); |
| 703 | } |
| 704 | |
| 705 | |
| 706 | Notice that the problem construction for bundle adjustment is very |
| 707 | similar to the curve fitting example -- one term is added to the |
| 708 | objective function per observation. |
| 709 | |
| 710 | Since this large sparse problem (well large for ``DENSE_QR`` anyways), |
| 711 | one way to solve this problem is to set |
| 712 | :member:`Solver::Options::linear_solver_type` to |
| 713 | ``SPARSE_NORMAL_CHOLESKY`` and call :member:`Solve`. And while this is |
| 714 | a reasonable thing to do, bundle adjustment problems have a special |
| 715 | sparsity structure that can be exploited to solve them much more |
| 716 | efficiently. Ceres provides three specialized solvers (collectively |
| 717 | known as Schur-based solvers) for this task. The example code uses the |
| 718 | simplest of them ``DENSE_SCHUR``. |
| 719 | |
| 720 | .. code-block:: c++ |
| 721 | |
| 722 | ceres::Solver::Options options; |
| 723 | options.linear_solver_type = ceres::DENSE_SCHUR; |
| 724 | options.minimizer_progress_to_stdout = true; |
| 725 | ceres::Solver::Summary summary; |
| 726 | ceres::Solve(options, &problem, &summary); |
| 727 | std::cout << summary.FullReport() << "\n"; |
| 728 | |
| 729 | For a more sophisticated bundle adjustment example which demonstrates |
| 730 | the use of Ceres' more advanced features including its various linear |
| 731 | solvers, robust loss functions and local parameterizations see |
| 732 | `examples/bundle_adjuster.cc |
| 733 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/bundle_adjuster.cc>`_ |
| 734 | |
| 735 | |
| 736 | .. rubric:: Footnotes |
| 737 | |
| 738 | .. [#f8] `examples/simple_bundle_adjuster.cc |
| 739 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/simple_bundle_adjuster.cc>`_ |
| 740 | |
| 741 | |
| 742 | Other Examples |
| 743 | ============== |
| 744 | |
| 745 | Besides the examples in this chapter, the `example |
| 746 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/>`_ |
| 747 | directory contains a number of other examples: |
| 748 | |
| 749 | #. `bundle_adjuster.cc |
| 750 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/bundle_adjuster.cc>`_ |
| 751 | shows how to use the various features of Ceres to solve bundle |
| 752 | adjustment problems. |
| 753 | |
| 754 | #. `circle_fit.cc |
| 755 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/circle_fit.cc>`_ |
| 756 | shows how to fit data to a circle. |
| 757 | |
| 758 | #. `denoising.cc |
| 759 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/denoising.cc>`_ |
| 760 | implements image denoising using the `Fields of Experts |
| 761 | <http://www.gris.informatik.tu-darmstadt.de/~sroth/research/foe/index.html>`_ |
| 762 | model. |
| 763 | |
| 764 | #. `nist.cc |
| 765 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/nist.cc>`_ |
| 766 | implements and attempts to solves the `NIST |
| 767 | <http://www.itl.nist.gov/div898/strd/nls/nls_main.shtm>`_ |
| 768 | non-linear regression problems. |
| 769 | |
| 770 | #. `libmv_bundle_adjuster.cc |
| 771 | <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/libmv_bundle_adjuster.cc>`_ |
| 772 | is the bundle adjustment algorithm used by `Blender <www.blender.org>`_/libmv. |