docs/powell.tex

%!TEX root = ceres-solver.tex
\chapter{Powell's Function}
\label{chapter:tutorial:powell}
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
\begin{align}
   f_1(x) &= x_1 + 10*x_2 \\
   f_2(x) &= \sqrt{5} * (x_3 - x_4)\\
   f_3(x) &= (x_2 - 2*x_3)^2\\
   f_4(x) &= \sqrt{10} * (x_1 - x_4)^2\\
	F(x) & = \left[f_1(x),\ f_2(x),\ f_3(x),\ f_4(x) \right]
\end{align}
$F(x)$ is a function of four parameters, and has four residuals. Now,
one way to solve this problem would be to define four
\texttt{CostFunction} objects that compute the residual and Jacobians. \eg the following code shows the implementation for $f_4(x)$.
\begin{minted}[mathescape]{c++}
class F4 : public ceres::SizedCostFunction<1, 4> {
 public:
  virtual ~F4() {}
  virtual bool Evaluate(double const* const* parameters,
                        double* residuals,
                        double** jacobians) const {
    double x1 = parameters[0][0];
    double x4 = parameters[0][3];
    // $f_4 = \sqrt{10} * (x_1 - x_4)^2$
    residuals[0] = sqrt(10.0) * (x1 - x4) * (x1 - x4)
    if (jacobians != NULL) {
      jacobians[0][0] = 2.0 * sqrt(10.0) * (x1 - x4);   // $\partial_{x_1}f_4(x)$
      jacobians[0][1] = 0.0;                            // $\partial_{x_2}f_4(x)$
      jacobians[0][2] = 0.0;                            // $\partial_{x_3}f_4(x)$
      jacobians[0][3] = -2.0 * sqrt(10.0) * (x1 - x4);  // $\partial_{x_4}f_4(x)$
    }
    return true;
  }
};
\end{minted}

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.

\section{Automatic Differentiation}
\label{sec:tutorial:autodiff}
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
\begin{minted}[frame=lines,mathescape]{c++}
class F4 {
 public:
  template <typename T> bool operator()(const T* const x1,
                                        const T* const x4,
                                        T* residual) const {
    // $f_4 = \sqrt{10} * (x_1 - x_4)^2$
    residual[0] = T(sqrt(10.0)) * (x1[0] - x4[0]) * (x1[0] - x4[0]);
    return true;
  }
};
\end{minted}

The important thing to note here is that \texttt{operator()} is a
templated method, which assumes that all its inputs and outputs are of
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.

Note also that the parameters are not packed
into a single array, they are instead passed as separate arguments to
\texttt{operator()}. Similarly we can define classes \texttt{F1,F2}
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}.}
\begin{minted}[mathescape]{c++}
double x1 =  3.0; double x2 = -1.0; double x3 =  0.0; double x4 =  1.0;
// Add residual terms to the problem using the using the autodiff
// wrapper to get the derivatives automatically. 
problem.AddResidualBlock(
  new ceres::AutoDiffCostFunction<F1, 1, 1, 1>(new F1), NULL, &x1, &x2);
problem.AddResidualBlock(
  new ceres::AutoDiffCostFunction<F2, 1, 1, 1>(new F2), NULL, &x3, &x4);
problem.AddResidualBlock(
  new ceres::AutoDiffCostFunction<F3, 1, 1, 1>(new F3), NULL, &x2, &x3)
problem.AddResidualBlock(
  new ceres::AutoDiffCostFunction<F4, 1, 1, 1>(new F4), NULL, &x1, &x4);
\end{minted}
A few things are worth noting in the code above. First, the object
being added to the \texttt{Problem} is an
\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.


Compiling and running \texttt{powell.cc} gives us:
\begin{minted}{bash}
Initial x1 = 3, x2 = -1, x3 = 0, x4 = 1
   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
   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
   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
   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
   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
   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
   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
   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
   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
   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
  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
  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
Ceres Solver Report: Iterations: 12, Initial cost: 1.075000e+02, \
Final cost: 2.865573e-13, Termination: GRADIENT_TOLERANCE.
Final x1 = 0.000583994, x2 = -5.83994e-05, x3 = 9.55401e-05, x4 = 9.55401e-05
\end{minted}
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}$.

\section{Numeric Differentiation}
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}.

We recommend that if possible,  automatic differentiation should be used. The use of
C++ templates makes automatic differentiation extremely efficient,
whereas numeric differentiation can be quite expensive, prone to
numeric errors and leads to slower convergence.