c++ - Usando derivadas como funciones en CppAD

gradients autodiff

Question

Estoy tratando de modificar el ejemplo here :

# include <cppad/cppad.hpp>
namespace { // ---------------------------------------------------------
// define the template function JacobianCases<Vector> in empty namespace
template <typename Vector>
bool JacobianCases()
{ bool ok = true;
using CppAD::AD;
using CppAD::NearEqual;
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
using CppAD::exp;
using CppAD::sin;
using CppAD::cos;
// domain space vector
size_t n = 2;
CPPAD_TESTVECTOR(AD<double>) X(n);
X[0] = 1.;
X[1] = 2.;
// declare independent variables and starting recording
CppAD::Independent(X);
// a calculation between the domain and range values
AD<double> Square = X[0] * X[0];
// range space vector
size_t m = 3;
CPPAD_TESTVECTOR(AD<double>) Y(m);
Y[0] = Square * exp( X[1] );
Y[1] = Square * sin( X[1] );
Y[2] = Square * cos( X[1] );
// create f: X -> Y and stop tape recording
CppAD::ADFun<double> f(X, Y);
// new value for the independent variable vector
Vector x(n);
x[0] = 2.;
x[1] = 1.;
// compute the derivative at this x
Vector jac( m * n );
jac = f.Jacobian(x);
/*
F''(x) = [ 2 * x[0] * exp(x[1]) , x[0] * x[0] * exp(x[1]) ]
[ 2 * x[0] * sin(x[1]) , x[0] * x[0] * cos(x[1]) ]
[ 2 * x[0] * cos(x[1]) , -x[0] * x[0] * sin(x[i]) ]
*/
ok &= NearEqual( 2.*x[0]*exp(x[1]), jac[0*n+0], eps99, eps99);
ok &= NearEqual( 2.*x[0]*sin(x[1]), jac[1*n+0], eps99, eps99);
ok &= NearEqual( 2.*x[0]*cos(x[1]), jac[2*n+0], eps99, eps99);
ok &= NearEqual( x[0] * x[0] *exp(x[1]), jac[0*n+1], eps99, eps99);
ok &= NearEqual( x[0] * x[0] *cos(x[1]), jac[1*n+1], eps99, eps99);
ok &= NearEqual(-x[0] * x[0] *sin(x[1]), jac[2*n+1], eps99, eps99);
return ok;
}
} // End empty namespace
# include <vector>
# include <valarray>
bool Jacobian(void)
{ bool ok = true;
// Run with Vector equal to three different cases
// all of which are Simple Vectors with elements of type double.
ok &= JacobianCases< CppAD::vector <double> >();
ok &= JacobianCases< std::vector <double> >();
ok &= JacobianCases< std::valarray <double> >();
return ok;
}

Estoy tratando de modificarlo de la siguiente manera:

Sea G el jac jacobiano que se calcula en este ejemplo, en la línea:

jac = f.Jacobian(x);

y, como en el ejemplo, sea X las variables independientes. Me gustaría construir una nueva función, H , que es una función de jac , es decir, H(jacobian(X)) = algo, de modo que H sea autodiferenciable. Un ejemplo puede ser H(X) = jacobian( jacobian(X)[0]) , es decir, el jacobian del primer elemento de jacobian(X) wrt X (una segunda derivada de clases).

El problema es que jac tal como está escrito aquí, es de tipo Vector , que es un tipo parametrizado en un double sin double , no en un AD<double> . Que yo sepa, esto significa que la salida no es autodiferenciable.

Estoy buscando algún consejo sobre si es posible usar el jacobiano en una operación más grande, y tomar el jacobiano de esa operación más grande (no diferente a cualquier operador aritmético) o si esto no es posible.

EDITAR: Esto ha sido ofrecido una vez por recompensa, pero lo estoy poniendo de nuevo para ver si hay una mejor solución, porque creo que esto es importante. Para ser un poco más claro, los elementos que necesita la respuesta "correcta" son:

a) Un medio de cálculo de derivados de orden arbitrario.

b) Una forma inteligente de no tener que especificar el orden de los derivados a priori. Si el derivado de orden máximo se debe conocer en el momento de la compilación, el orden de derivado no se puede determinar algorítmicamente. Además, especificar un orden enormemente grande como en la respuesta actual dada dará lugar a problemas de asignación de memoria y, me imagino, problemas de rendimiento.

c) abstraer la plantilla de orden derivada lejos del usuario final. Esto es importante, ya que puede ser difícil hacer un seguimiento del orden de los derivados necesarios. Esto es probablemente algo que viene "gratis" si b) se resuelve.

Si alguien puede resolver esto, sería una contribución increíble y una operación extremadamente útil.

Answer 1

Hay una nueva función en CppAD que elimina la necesidad de AD <AD>, consulte https://coin-or.github.io/CppAD/doc/base2ad.cpp.htm

Answer 2

Si desea anidar funciones, también debe anidar AD <>. Puede anidar a los jacobianos como a otras funciones, por ejemplo, vea el fragmento de código a continuación, que calcula la doble derivada anidando a los jacobianos.

#include <cstring>
#include <iostream> // standard input/output
#include <vector> // standard vector
#include <cppad/cppad.hpp> // the CppAD package http://www.coin-or.org/CppAD/
// main program
int main(void)
{ using CppAD::AD; // use AD as abbreviation for CppAD::AD
using std::vector; // use vector as abbreviation for std::vector
size_t i; // a temporary index
// domain space vector
auto Square = [](auto t){return t*t;};
vector< AD<AD<double>> > X(1); // vector of domain space variables
// declare independent variables and start recording operation sequence
CppAD::Independent(X);
// range space vector
vector< AD<AD<double>> > Y(1); // vector of ranges space variables
Y[0] = Square(X[0]); // value during recording of operations
// store operation sequence in f: X -> Y and stop recording
CppAD::ADFun<AD<double>> f(X, Y);
// compute derivative using operation sequence stored in f
vector<AD<double>> jac(1); // Jacobian of f (m by n matrix)
vector<AD<double>> x(1); // domain space vector
CppAD::Independent(x);
jac = f.Jacobian(x); // Jacobian for operation sequence
CppAD::ADFun<double> f2(x, jac);
vector<double> result(1);
vector<double> x_res(1);
x_res[0]=15.;
result=f2.Jacobian(x_res);
// print the results
std::cout << "f'''' computed by CppAD = " << result[0] << std::endl;
}

Como nota al margen, desde C ++ 14 u 11 la implementación de plantillas de expresión y la diferenciación automática se hizo más fácil y se puede hacer con mucho menos esfuerzo, como se muestra, por ejemplo, en este video hacia el final https://www.youtube.com/watch?v=cC9MtflQ_nI (perdón por la mala calidad). Si tuviera que implementar operaciones simbólicas razonablemente simples, comenzaría desde cero con C ++ moderno: puede escribir código más simple y obtener errores que puede comprender fácilmente.

Edición: Generalizar el ejemplo para construir derivados de orden arbitrarios puede ser un ejercicio de metaprogramación de plantilla. El siguiente fragmento de código muestra que es posible usar la recursión de la plantilla.

#include <cstring>
#include <iostream>
#include <vector>
#include <cppad/cppad.hpp>
using CppAD::AD;
using std::vector;
template<typename T>
struct remove_ad{
using type=T;
};
template<typename T>
struct remove_ad<AD<T>>{
using type=T;
};
template<int N>
struct derivative{
using type = AD<typename derivative<N-1>::type >;
static constexpr int order = N;
};
template<>
struct derivative<0>{
using type = double;
static constexpr int order = 0;
};
template<typename T>
struct Jac{
using value_type = typename remove_ad<typename T::type>::type;
template<typename P, typename Q>
auto operator()(P & X, Q & Y){
CppAD::ADFun<value_type> f(X, Y);
vector<value_type> jac(1);
vector<value_type> x(1);
CppAD::Independent(x);
jac = f.Jacobian(x);
return Jac<derivative<T::order-1>>{}(x, jac);
}
};
template<>
struct Jac<derivative<1>>{
using value_type = derivative<0>::type;
template<typename P, typename Q>
auto operator()(P & x, Q & jac){
CppAD::ADFun<value_type> f2(x, jac);
vector<value_type> res(1);
vector<value_type> x_res(1);
x_res[0]=15.;
return f2.Jacobian(x_res);
}
};
int main(void)
{
constexpr int order=4;
auto Square = [](auto t){return t*t;};
vector< typename derivative<order>::type > X(1);
vector< typename derivative<order>::type > Y(1);
CppAD::Independent(X);
Y[0] = Square(X[0]);
auto result = Jac<derivative<order>>{}(X, Y);
std::cout << "f'''' computed by CppAD = " << result[0] << std::endl;
}