docs/powell.tex - platform/external/ceres-solver - Git at Google

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