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}[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}
In some cases, its not possible to define a templated cost functor. In such a situation, numerical differentiation can be used. The user defines a functor which computes the residual value and construct a \texttt{NumericDiffCostFunction} using it. e.g., for \texttt{F4}, the corresponding functor would be
\begin{minted}[mathescape]{c++}
class F4 {
 public:
  bool operator()(const double* const x1,
                  const double* const x4,
                  double* residual) const {
    // $f_4 = \sqrt{10} * (x_1 - x_4)^2$
    residual[0] = sqrt(10.0) * (x1[0] - x4[0]) * (x1[0] - x4[0]);
    return true;
  }
};
\end{minted}

Which can then be wrapped \texttt{NumericDiffCostFunction} and added to the \texttt{Problem} object as follows

\begin{minted}[mathescape]{c++}
problem.AddResidualBlock(
  new ceres::NumericDiffCostFunction<F4, ceres::CENTRAL, 1, 1, 1>(new F4), NULL, &x1, &x4);
\end{minted}

The construction looks almost identical to the used for automatic differentiation, except for an extra template parameter that indicates the kind of finite differencing scheme to be used for computing the numerical derivatives. \texttt{examples/quadratic\_numeric\_diff.cc} shows a numerically differentiated implementation of \texttt{examples/quadratic.cc}.

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