c++ - 在 CppAD 中使用导数作为函数

Question

我正在尝试修改示例 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;
}

我正在尝试通过以下方式对其进行修改:

设 G 为雅可比 jac在此示例中计算，在以下行中:

jac = f.Jacobian(x);

并且，如示例所示，让 X成为自变量。我想构造一个新函数，H ，它是 jac 的函数，即 H(jacobian(X)) = 某物，使得 H 是可自微的。一个例子可能是 H(X) = jacobian( jacobian(X)[0]) ，即 jacobian(X) 的第一个元素的雅可比w.r.t X (某种二阶导数)。

问题是jac写在这里是类型Vector ，它是原始 double 上的参数化类型，而不是 AD<double> .据我所知，这意味着输出不可自微。

我正在寻找一些关于是否可以在更大的运算中使用雅可比行列式的建议，并采用该较大运算的雅可比行列式(与任何算术运算符不同)或者这是否不可能。

编辑:这已经提出了一次赏金，但我再次提出来看看是否有更好的解决方案，因为我认为这很重要。更清楚一点，“正确”答案需要的元素是:

a) 一种计算任意阶导数的方法。

b) 一种不必先验地指定导数顺序的智能方法。如果在编译时必须知道最大阶导数，则无法通过算法确定导数的阶数。此外，在给出的当前答案中指定一个非常大的订单将导致内存分配问题，我想还会导致性能问题。

c) 将衍生订单的模板从最终用户中抽象出来。这很重要，因为很难跟踪所需的衍生品的顺序。如果 b) 得到解决，这可能是“免费”提供的东西。

如果有人能破解这个，那将是一个了不起的贡献和非常有用的操作。

Answer 1

如果你想嵌套函数，你也应该嵌套 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;
}