use*_*261 28 c++ gradient autodiff
我想在这里修改示例:
# 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;
}
我试图通过以下方式修改它:
设G是jac本例中计算的雅可比行列式:
jac = f.Jacobian(x);
并且,如在示例中,让我们X成为自变量.我想构造一个新函数,H它是一个函数jac,即H(jacobian(X))=某事,使得H是自动不同的.一个例子可以是H(X) = jacobian( jacobian(X)[0]),即jacobian(X)wrt 的第一个元素的jacobian X(各种二阶导数).
问题是,jac这里写的是类型Vector,它是原始的参数化类型double,而不是AD<double>.据我所知,这意味着输出不是自动不同的.
我正在寻找一些建议,如果有可能在更大的操作中使用雅可比行列式,并采用更大的操作的雅可比行列式(与任何算术运算符不同)或者如果这是不可能的话.
编辑:这已经获得了一次赏金,但我再次提出来看看是否有更好的解决方案,因为我觉得这很重要.更清楚一点,"正确"答案所需的要素是:
a)计算任意阶导数的方法.
b)一种不必先验地指定导数顺序的智能方法.如果必须在编译时知道最大阶导数,则无法通过算法确定导数的阶数.此外,指定一个非常大的订单,如当前给出的答案将导致内存分配问题,我想,性能问题.
c)从最终用户中提取衍生订单的模板.这很重要,因为很难跟踪所需衍生物的顺序.如果b)得到解决,这可能是"免费"的.
如果任何人都可以解决这个问题,那将是一项非常有用的贡献和非常有用的操作.
如果要嵌套函数，则还应该嵌套AD <>。您可以将Jacobian嵌套为其他函数，例如，参见下面的代码段，该代码段通过嵌套Jacobian计算双导数
#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;
}
附带说明一下，由于C ++ 14或11实现表达式模板和自动区分变得更加容易，并且可以用更少的精力完成，如此视频中最后所示https://www.youtube.com/watch ？v = cC9MtflQ_nI（对不起，质量不好）。如果必须实现相当简单的符号操作，则可以从现代C ++从零开始：您可以编写更简单的代码，并且会得到容易理解的错误。
编辑： 泛化示例以构建任意阶导数可以是模板元编程练习。以下代码段显示可以使用模板递归
#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;
}