CppAD中基于Jacobian构建可自动微分高阶函数的实现咨询

你遇到的核心问题是直接调用f.Jacobian(x)得到的是数值结果（double类型），丢失了自动微分所需的AD信息，因此无法对其进一步求导。下面提供两种满足你需求的方案，都能支持任意阶导数、无需预先指定阶数，且可以封装细节对用户隐藏模板复杂度。

方案一：AD嵌套（AD of AD）—— 支持任意复合函数 #

这种方法通过嵌套AD类型（如AD<AD<double>>），将Jacobian的计算过程也纳入自动微分的tape中，从而让Jacobian元素成为可微分的AD变量，进而构建任意复杂的复合函数H并求导。

实现步骤与代码示例 #

1. 先封装基础函数f，使其能接受任意模板类型的输入（支持AD嵌套）：

#include <cppad/cppad.hpp>
#include <vector>
#include <limits>

template <typename T>
CppAD::vector<T> base_f(const CppAD::vector<T>& X) {
    size_t m = 3;
    CppAD::vector<T> Y(m);
    T Square = X[0] * X[0];
    Y[0] = Square * CppAD::exp(X[1]);
    Y[1] = Square * CppAD::sin(X[1]);
    Y[2] = Square * CppAD::cos(X[1]);
    return Y;
}

2. 构建复合函数H（此处以H(X) = df₀/dx₀(X)为例，即Jacobian的第一个元素），并对H求导：

bool nested_ad_example() {
    bool ok = true;
    using CppAD::AD;
    using CppAD::NearEqual;
    double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
    size_t n = 2;

    // 记录H函数的tape：输入是AD<double>类型的X，输出是df0/dx0(X)
    CppAD::vector<AD<double>> X_ad(n);
    X_ad[0] = 2.0;
    X_ad[1] = 1.0;
    CppAD::Independent(X_ad);

    // 构建支持AD<double>输入的f版本，此时f的Jacobian结果是AD<double>类型
    CppAD::vector<AD<double>> X_copy = X_ad;
    CppAD::ADFun<AD<double>> f_ad;
    f_ad.Dependent(X_copy, base_f(X_copy));

    // 计算f在X_ad处的Jacobian（AD<double>类型，保留微分信息）
    CppAD::vector<AD<double>> jac_ad = f_ad.Jacobian(X_ad);

    // 定义H的输出为Jacobian的第一个元素
    CppAD::vector<AD<double>> H_output(1);
    H_output[0] = jac_ad[0];

    // 构建H函数：X → H(X) = df0/dx0(X)
    CppAD::ADFun<double> H(X_ad, H_output);

    // 对H求Jacobian，得到二阶导数d(df0/dx0)/dX
    CppAD::vector<double> x(n);
    x[0] = 2.0;
    x[1] = 1.0;
    CppAD::vector<double> H_jac = H.Jacobian(x);

    // 验证结果：d(df0/dx0)/dx0 = 2*exp(x[1])，d(df0/dx0)/dx1 = 2*x[0]*exp(x[1])
    ok &= NearEqual(2*CppAD::exp(x[1]), H_jac[0], eps99, eps99);
    ok &= NearEqual(2*x[0]*CppAD::exp(x[1]), H_jac[1], eps99, eps99);

    // 若需三阶导数，只需继续嵌套AD类型（如AD<AD<AD<double>>>），重复上述流程即可
    return ok;
}

方案优势 #

• 支持任意复杂的复合函数H（比如Jacobian元素的加减乘除、嵌套其他函数）；
• 无需提前指定最大导数阶数，每次仅处理当前所需阶数，避免内存浪费；
• 模板细节可封装，用户只需调用封装后的接口，无需接触嵌套AD类型。

方案二：Forward/Reverse高阶微分接口—— 高效计算高阶导数 #

CppAD原生提供Forward和Reverse方法，可直接计算任意阶导数，无需显式构建嵌套AD函数，性能更优。

实现步骤与代码示例 #

bool forward_reverse_example() {
    bool ok = true;
    using CppAD::AD;
    using CppAD::NearEqual;
    double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
    size_t n = 2;
    size_t m = 3;

    // 构建原始函数f的tape
    CppAD::vector<AD<double>> X(n);
    X[0] = 1.0;
    X[1] = 2.0;
    CppAD::Independent(X);
    AD<double> Square = X[0] * X[0];
    CppAD::vector<AD<double>> Y(m);
    Y[0] = Square * CppAD::exp(X[1]);
    Y[1] = Square * CppAD::sin(X[1]);
    Y[2] = Square * CppAD::cos(X[1]);
    CppAD::ADFun<double> f(X, Y);

    // 目标：计算d(df0/dx0)/dx0 和 d(df0/dx0)/dx1（二阶导数）
    CppAD::vector<double> x(n);
    x[0] = 2.0;
    x[1] = 1.0;

    // 先将f的输入设置为x
    f.Forward(0, x);

    // 用Reverse模式计算二阶导数：权重向量w指定我们关心f0的导数
    CppAD::vector<double> w(m, 0.0);
    w[0] = 1.0;
    size_t p = 2; // 求到二阶导数
    CppAD::vector<double> r = f.Reverse(p, w);

    // 解析结果：r[j*(p+1)+k]表示x_j的k阶导数系数
    double d2f0_dx0dx0 = r[0*(p+1) + 2]; // 二阶导数d²f0/dx0²
    double d2f0_dx0dx1 = r[0*(p+1) + 2 + 1]; // 二阶导数d²f0/(dx0 dx1)

    // 验证结果
    ok &= NearEqual(d2f0_dx0dx0, 2*CppAD::exp(x[1]), eps99, eps99);
    ok &= NearEqual(d2f0_dx0dx1, 2*x[0]*CppAD::exp(x[1]), eps99, eps99);

    // 若需三阶导数，只需将p设为3，解析对应的r元素即可
    return ok;
}

方案优势 #

• 无需嵌套AD类型，内存占用更低，计算效率更高；
• 直接通过参数p指定所需导数阶数，无需提前配置最大阶数；
• 适合直接计算原始函数的高阶导数，逻辑简洁。

封装模板细节的示例 #

你可以将上述逻辑封装为类，对终端用户隐藏所有AD模板和高阶接口细节：

class HighOrderDerivCalculator {
private:
    CppAD::ADFun<double> f_;
    size_t n_;
    size_t m_;

public:
    template <typename Func>
    HighOrderDerivCalculator(size_t n, size_t m, Func func) : n_(n), m_(m) {
        CppAD::vector<AD<double>> X(n);
        for (size_t i = 0; i < n; ++i) X[i] = 0.0;
        CppAD::Independent(X);
        CppAD::vector<AD<double>> Y = func(X);
        f_.Dependent(X, Y);
    }

    // 计算output_idx对应的输出函数，关于input_idxs序列的k阶导数
    double calc(size_t k, size_t output_idx, const std::vector<size_t>& input_idxs, const std::vector<double>& x) {
        f_.Forward(0, x);
        CppAD::vector<double> w(m_, 0.0);
        w[output_idx] = 1.0;
        CppAD::vector<double> r = f_.Reverse(k, w);

        // 此处仅实现二阶导数的解析，更高阶可扩展
        if (k == 2 && input_idxs.size() == 2) {
            size_t j = input_idxs[0];
            size_t l = input_idxs[1];
            if (j == l) return r[j*(k+1) + k];
            else return r[j*(k+1) + k + (l - j)];
        }
        return 0.0;
    }
};

// 使用示例
bool encapsulated_example() {
    bool ok = true;
    size_t n = 2, m = 3;
    auto func = [](const CppAD::vector<AD<double>>& X) {
        CppAD::vector<AD<double>> Y(3);
        AD<double> Square = X[0] * X[0];
        Y[0] = Square * CppAD::exp(X[1]);
        Y[1] = Square * CppAD::sin(X[1]);
        Y[2] = Square * CppAD::cos(X[1]);
        return Y;
    };

    HighOrderDerivCalculator calc(n, m, func);
    std::vector<double> x = {2.0, 1.0};
    double d2f0_dx0dx0 = calc.calc(2, 0, {0,0}, x);
    ok &= CppAD::NearEqual(d2f0_dx0dx0, 2*CppAD::exp(1.0), 99*std::numeric_limits<double>::epsilon(), 99*std::numeric_limits<double>::epsilon());
    return ok;
}

内容的提问来源于stack exchange，提问作者user650261