%!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.