/* This is an example on how automatic differentiation can be used in
a multidimensional Newton iteration, which eliminates need for
(often painful) symbolic derivatives.
Copyright (c) 2009 by Wenzel Jakob
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include
#include "autodiff.h"
#include "cholesky.h"
using namespace std;
int main(int argc, char **argv) {
/* Instantiate a differentiation environment with 2 independent variables.
(this is placed on the heap - alternatively, the typedef
could be changed to
typedef DiffEnvFixed DEnv;
in which case the variable count is fixed at compile time and the
DEnv instance is placed on the stack instead)
*/
typedef DiffEnv DEnv;
DEnv env(2);
/* Do a multidimensional newton iteration for the function
f(x1, x2) = -log(1-x1-x2) - log(x1) - log(x2)
starting at x1 = .85 and x2 = .05
*/
env[0] = 0.85; env[1] = 0.05;
for (int i=0; i<10; ++i) {
/* Print the current iterate -- it should end
up converging to [1/3, 1/3]
*/
cout << "Iteration " << i << " (x1=" << env[0] << ", x2="
<< env[1] << ")" << endl;
/* Get a reference to the independent variables. These are basically
large data structures, which contain the value, gradient and
hessian of the two functions
f(x1,x2) = x1 (=> gradient=[1; 0], hessian = [0 0; 0 0])
f(x1,x2) = x2 (=> gradient=[0; 1], hessian = [0 0; 0 0])
evaluated at the current x1,x2 values.
*/
DEnv::Scalar x1 = env.getIndepScalar(0);
DEnv::Scalar x2 = env.getIndepScalar(1);
/* Compute f(x1, x2). This additionally propagates first and
second derivative information through builtin operators and
math library calls.
*/
DEnv::Scalar Fx = -log(1-x1-x2) - log(x1) - log(x2);
/* Print the local derivative information */
cout << " - Local information: " << Fx << endl << endl;
/* Compute the Cholesky decomposition of f's Hessian
and use it to take a full Newton step */
DEnv::Hessian hess = Fx.getHessian();
DEnv::Gradient grad = Fx.getGradient();
ublas::triangular_matrix tri(2, 2);
if (!cholesky_decompose(hess, tri))
throw std::range_error("Matrix is not symmetric positive definite!");
cholesky_solve(tri, grad);
/* Update the differentiation environment with new coefficients */
env[0] -= grad(0);
env[1] -= grad(1);
}
return 0;
}