logo标签列表

在CppAD中使用导数作为函数

c++gradientautodiff
28

我想修改这里的示例：

# 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)) = something，使得 H 是自动微分的。一个示例可能是 H(X) = jacobian( jacobian(X)[0])，即第一个元素相对于 X 的雅可比矩阵（一种二阶导数）。

问题在于这里写的 jac 是类型为 Vector 的，它是基于原始的 double 参数化的类型，而不是 AD<double>。据我的了解，这意味着输出结果不是自动微分的。

我正在寻找一些建议，看是否有可能在更大的操作中使用雅可比矩阵，并将该大型操作的雅可比矩阵取出（类似于任何算术运算符），或者这不可能。

编辑：这个问题曾经被发起过一次悬赏，但我现在再次发起悬赏，看是否有更好的解决方案，因为我认为这很重要。稍微明确一下，“正确”的答案需要以下元素：

a）计算任意阶导数的方法。

b）不必事先指定导数的阶数的智能方法。如果必须在编译时知道最大阶数导数，则无法通过算法确定导数的阶数。此外，像目前给出的答案中指定一个极大阶数会导致内存分配问题和性能问题。

c）将导数阶数的模板化抽象化为最终用户。这很重要，因为很难跟踪所需的导数的阶数。如果解决了 b），这可能是一种“免费”的东西。

如果有人能解决这个问题，那将是一个很棒的贡献和一个非常有用的操作。

1
我已经添加了所请求的代码。对此感到抱歉！ - user650261
到目前为止，看起来多级可能是正确的答案（例如在 https://coin-or.github.io/CppAD/doc/mul_level.cpp.htm 中），但我不清楚如何将其应用于任意次数，这正是我所需要的。 - user650261
2个回答
6

如果您想嵌套函数，那么您也应该嵌套AD<>。您可以像嵌套其他函数一样嵌套雅各比矩阵，例如请参见下面的代码片段，它通过嵌套雅各比矩阵来计算二阶导数。

#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;
}
谢谢提供示例。这个几乎可以工作 - 除了你必须事先知道你想要的微分级数。有没有一种方法可以在不事先指定n的情况下进行第n阶导数？ - user650261
1
我添加了另一个代码片段，其中计算任意顺序派生。 - Paolo Crosetto
谢谢。我还没有验证这个（希望明天能够），但是我已经给你发放了奖励，因为我在结束之前无法检查它。 - user650261
1
我相信这个结构实际上存在错误，因为对于order=2，它输出的“f'' computed by CppAD = 0”，这是错误的；它应该是2.0。 - user650261
1
哦，对了，有个 bug。我调用 Independent(X) 和 Y[0]=f(X) 的顺序错了。我已经修复了。稍后我会添加一些额外的解释。 - Paolo Crosetto
显示剩余2条评论
1
尽管这个回答在技术上回答了问题，但如果它包含任何相关代码的话，它会更好。 - Ruzihm
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接