1 %!TEX root = ceres-solver.tex 2 \chapter{Powell's Function} 3 \label{chapter:tutorial:powell} 4 Consider now a slightly more complicated example -- the minimization of Powell's function. Let $x = \left[x_1, x_2, x_3, x_4 \right]$ and 5 \begin{align} 6 f_1(x) &= x_1 + 10*x_2 \\ 7 f_2(x) &= \sqrt{5} * (x_3 - x_4)\\ 8 f_3(x) &= (x_2 - 2*x_3)^2\\ 9 f_4(x) &= \sqrt{10} * (x_1 - x_4)^2\\ 10 F(x) & = \left[f_1(x),\ f_2(x),\ f_3(x),\ f_4(x) \right] 11 \end{align} 12 $F(x)$ is a function of four parameters, and has four residuals. Now, 13 one way to solve this problem would be to define four 14 \texttt{CostFunction} objects that compute the residual and Jacobians. \eg the following code shows the implementation for $f_4(x)$. 15 \begin{minted}[mathescape]{c++} 16 class F4 : public ceres::SizedCostFunction<1, 4> { 17 public: 18 virtual ~F4() {} 19 virtual bool Evaluate(double const* const* parameters, 20 double* residuals, 21 double** jacobians) const { 22 double x1 = parameters[0][0]; 23 double x4 = parameters[0][3]; 24 // $f_4 = \sqrt{10} * (x_1 - x_4)^2$ 25 residuals[0] = sqrt(10.0) * (x1 - x4) * (x1 - x4) 26 if (jacobians != NULL) { 27 jacobians[0][0] = 2.0 * sqrt(10.0) * (x1 - x4); // $\partial_{x_1}f_4(x)$ 28 jacobians[0][1] = 0.0; // $\partial_{x_2}f_4(x)$ 29 jacobians[0][2] = 0.0; // $\partial_{x_3}f_4(x)$ 30 jacobians[0][3] = -2.0 * sqrt(10.0) * (x1 - x4); // $\partial_{x_4}f_4(x)$ 31 } 32 return true; 33 } 34 }; 35 \end{minted} 36 37 But this can get painful very quickly, especially for residuals involving complicated multi-variate terms. Ceres provides two ways around this problem. Numeric and automatic symbolic differentiation. 38 39 \section{Automatic Differentiation} 40 \label{sec:tutorial:autodiff} 41 With its automatic differentiation support, Ceres allows you to define templated objects/functors that will compute the residual and it takes care of computing the Jacobians as needed and filling the \texttt{jacobians} arrays with them. For example, for $f_4(x)$ we define 42 \begin{minted}[frame=lines,mathescape]{c++} 43 class F4 { 44 public: 45 template <typename T> bool operator()(const T* const x1, 46 const T* const x4, 47 T* residual) const { 48 // $f_4 = \sqrt{10} * (x_1 - x_4)^2$ 49 residual[0] = T(sqrt(10.0)) * (x1[0] - x4[0]) * (x1[0] - x4[0]); 50 return true; 51 } 52 }; 53 \end{minted} 54 55 The important thing to note here is that \texttt{operator()} is a 56 templated method, which assumes that all its inputs and outputs are of 57 some type \texttt{T}. The reason for using templates here is because Ceres will call \texttt{F4::operator<T>()}, with $\texttt{T=double}$ when just the residual is needed, and with a special type $T=\texttt{Jet}$ when the Jacobians are needed. 58 59 Note also that the parameters are not packed 60 into a single array, they are instead passed as separate arguments to 61 \texttt{operator()}. Similarly we can define classes \texttt{F1,F2} 62 and \texttt{F4}. Then let us consider the construction and solution of the problem. For brevity we only describe the relevant bits of code~\footnote{The full source code for this example can be found in \texttt{examples/powell.cc}.} 63 \begin{minted}[mathescape]{c++} 64 double x1 = 3.0; double x2 = -1.0; double x3 = 0.0; double x4 = 1.0; 65 // Add residual terms to the problem using the using the autodiff 66 // wrapper to get the derivatives automatically. 67 problem.AddResidualBlock( 68 new ceres::AutoDiffCostFunction<F1, 1, 1, 1>(new F1), NULL, &x1, &x2); 69 problem.AddResidualBlock( 70 new ceres::AutoDiffCostFunction<F2, 1, 1, 1>(new F2), NULL, &x3, &x4); 71 problem.AddResidualBlock( 72 new ceres::AutoDiffCostFunction<F3, 1, 1, 1>(new F3), NULL, &x2, &x3) 73 problem.AddResidualBlock( 74 new ceres::AutoDiffCostFunction<F4, 1, 1, 1>(new F4), NULL, &x1, &x4); 75 \end{minted} 76 A few things are worth noting in the code above. First, the object 77 being added to the \texttt{Problem} is an 78 \texttt{AutoDiffCostFunction} with \texttt{F1}, \texttt{F2}, \texttt{F3} and \texttt{F4} as template parameters. Second, each \texttt{ResidualBlock} only depends on the two parameters that the corresponding residual object depends on and not on all four parameters. 79 80 81 Compiling and running \texttt{powell.cc} gives us: 82 \begin{minted}{bash} 83 Initial x1 = 3, x2 = -1, x3 = 0, x4 = 1 84 0: f: 1.075000e+02 d: 0.00e+00 g: 1.55e+02 h: 0.00e+00 rho: 0.00e+00 mu: 1.00e-04 li: 0 85 1: f: 5.036190e+00 d: 1.02e+02 g: 2.00e+01 h: 2.16e+00 rho: 9.53e-01 mu: 3.33e-05 li: 1 86 2: f: 3.148168e-01 d: 4.72e+00 g: 2.50e+00 h: 6.23e-01 rho: 9.37e-01 mu: 1.11e-05 li: 1 87 3: f: 1.967760e-02 d: 2.95e-01 g: 3.13e-01 h: 3.08e-01 rho: 9.37e-01 mu: 3.70e-06 li: 1 88 4: f: 1.229900e-03 d: 1.84e-02 g: 3.91e-02 h: 1.54e-01 rho: 9.37e-01 mu: 1.23e-06 li: 1 89 5: f: 7.687123e-05 d: 1.15e-03 g: 4.89e-03 h: 7.69e-02 rho: 9.37e-01 mu: 4.12e-07 li: 1 90 6: f: 4.804625e-06 d: 7.21e-05 g: 6.11e-04 h: 3.85e-02 rho: 9.37e-01 mu: 1.37e-07 li: 1 91 7: f: 3.003028e-07 d: 4.50e-06 g: 7.64e-05 h: 1.92e-02 rho: 9.37e-01 mu: 4.57e-08 li: 1 92 8: f: 1.877006e-08 d: 2.82e-07 g: 9.54e-06 h: 9.62e-03 rho: 9.37e-01 mu: 1.52e-08 li: 1 93 9: f: 1.173223e-09 d: 1.76e-08 g: 1.19e-06 h: 4.81e-03 rho: 9.37e-01 mu: 5.08e-09 li: 1 94 10: f: 7.333425e-11 d: 1.10e-09 g: 1.49e-07 h: 2.40e-03 rho: 9.37e-01 mu: 1.69e-09 li: 1 95 11: f: 4.584044e-12 d: 6.88e-11 g: 1.86e-08 h: 1.20e-03 rho: 9.37e-01 mu: 5.65e-10 li: 1 96 Ceres Solver Report: Iterations: 12, Initial cost: 1.075000e+02, \ 97 Final cost: 2.865573e-13, Termination: GRADIENT_TOLERANCE. 98 Final x1 = 0.000583994, x2 = -5.83994e-05, x3 = 9.55401e-05, x4 = 9.55401e-05 99 \end{minted} 100 It is easy to see that the optimal solution to this problem is at $x_1=0, x_2=0, x_3=0, x_4=0$ with an objective function value of $0$. In 10 iterations, Ceres finds a solution with an objective function value of $4\times 10^{-12}$. 101 102 \section{Numeric Differentiation} 103 If a templated implementation is not possible then a \texttt{NumericDiffCostFunction} object can be used. The user defines a \texttt{CostFunction} object whose \texttt{Evaluate} method is only computes the residuals. A wrapper object \texttt{NumericDiffCostFunction} then uses it to compute the residuals and the Jacobian using finite differencing. \texttt{examples/quadratic\_numeric\_diff.cc} shows a numerically differentiated implementation of \texttt{examples/quadratic.cc}. 104 105 We recommend that if possible, automatic differentiation should be used. The use of 106 C++ templates makes automatic differentiation extremely efficient, 107 whereas numeric differentiation can be quite expensive, prone to 108 numeric errors and leads to slower convergence.
