如何在Ceres-solver的EvaluationCallback中使用AutoDiffCostFunction

我是Ceres-solver新手，需要拟合含参数a₁、a₂、a₃、a₄、a₅的非线性曲线f(x; a₁,a₂,a₃,a₄,a₅)到观测数据集y(x)。该函数数学形式复杂，无法手动推导导数闭式表达式，希望使用自动求导，且已能通过模板实现该函数。

我有数千个观测点（如x∈[-2000,…2000]对应的y(x)值），批量计算函数值比单点计算效率更高。

查阅文档后，我认为可参考ceres-solver示例evaluation_callback_example.cc，即：

• 实现继承自EvaluationCallback的类，以高效批量计算所有x对应的残差/雅可比矩阵
• 实现CostAndJacobianCopyingCostFunction，返回单点x对应的残差和雅可比矩阵

但我不确定如何修改该示例，使EvaluationCallback中使用ceres::AutoDiffCostFunction以避免手动计算导数/雅可比矩阵（示例中使用了雅可比矩阵的闭式表达式，而我的函数无法做到）。

示例中的EvaluationCallback实现代码：

// This implementation of the EvaluationCallback interface also stores the
// residuals and Jacobians that the CostFunction copies their values from.
class MyEvaluationCallback : public ceres::EvaluationCallback {
 public:
  // m and c are passed by reference so that we have access to their values as
  // they evolve over time through the course of optimization.
  MyEvaluationCallback(const double& m, const double& c) : m_(m), c_(c) {
    x_ = Eigen::VectorXd::Zero(kNumObservations);
    y_ = Eigen::VectorXd::Zero(kNumObservations);
    residuals_ = Eigen::VectorXd::Zero(kNumObservations);
    jacobians_ = Eigen::MatrixXd::Zero(kNumObservations, 2);
    for (int i = 0; i < kNumObservations; ++i) {
      x_[i] = data[2 * i];
      y_[i] = data[2 * i + 1];
    }
    PrepareForEvaluation(true, true);
  }

  void PrepareForEvaluation(bool evaluate_jacobians,
                            bool new_evaluation_point) final {
    if (new_evaluation_point) {
      ComputeResidualAndJacobian(evaluate_jacobians);
      jacobians_are_stale_ = !evaluate_jacobians;
    } else {
      if (evaluate_jacobians && jacobians_are_stale_) {
        ComputeResidualAndJacobian(evaluate_jacobians);
        jacobians_are_stale_ = false;
      }
    }
  }

  const Eigen::VectorXd& residuals() const { return residuals_; }
  const Eigen::MatrixXd& jacobians() const { return jacobians_; }
  bool jacobians_are_stale() const { return jacobians_are_stale_; }

 private:
  void ComputeResidualAndJacobian(bool evaluate_jacobians) {
    residuals_ = -(m_ * x_.array() + c_).exp();
    if (evaluate_jacobians) {
      jacobians_.col(0) = residuals_.array() * x_.array();
      jacobians_.col(1) = residuals_;
    }
    residuals_ += y_;
  }

  const double& m_;
  const double& c_;
  Eigen::VectorXd x_;
  Eigen::VectorXd y_;
  Eigen::VectorXd residuals_;
  Eigen::MatrixXd jacobians_;

  // jacobians_are_stale_ keeps track of whether the jacobian matrix matches the
  // residuals or not, we only compute it if we know that Solver is going to
  // need access to it.
  bool jacobians_are_stale_ = true;
};

核心思路是在EvaluationCallback的批量计算逻辑中，利用Ceres的自动求导机制，对每个观测点计算残差和雅可比，再统一存储供后续CostFunction拷贝。具体步骤如下：

1. 定义残差模板函数 #

先实现符合Ceres自动求导要求的残差模板函数：

template<typename T>
bool ResidualFunc(const T* params, const T& x, const T& y, T* residual) {
    // 替换为你的非线性函数实现，params[0]-params[4]对应a₁到a₅
    residual[0] = y - YourNonlinearFunction(x, params[0], params[1], params[2], params[3], params[4]);
    return true;
}

2. 重写EvaluationCallback类 #

去掉手动计算雅可比的逻辑，改用AutoDiffCostFunction批量计算：

class BatchEvaluationCallback : public ceres::EvaluationCallback {
public:
    BatchEvaluationCallback(const double* params, const Eigen::VectorXd& x, const Eigen::VectorXd& y)
        : params_(params), x_(x), y_(y) {
        const int num_obs = x_.size();
        residuals_ = Eigen::VectorXd::Zero(num_obs);
        jacobians_ = Eigen::MatrixXd::Zero(num_obs, 5); // 5个参数，每行对应一个观测点的梯度
        PrepareForEvaluation(true, true);
    }

    void PrepareForEvaluation(bool evaluate_jacobians, bool new_evaluation_point) final {
        if (new_evaluation_point || (evaluate_jacobians && jacobians_are_stale_)) {
            ComputeBatchResidualAndJacobian(evaluate_jacobians);
            jacobians_are_stale_ = !evaluate_jacobians;
        }
    }

    const Eigen::VectorXd& residuals() const { return residuals_; }
    const Eigen::MatrixXd& jacobians() const { return jacobians_; }

private:
    void ComputeBatchResidualAndJacobian(bool evaluate_jacobians) {
        const int num_obs = x_.size();
        using CostFunc = ceres::AutoDiffCostFunction<decltype(&ResidualFunc), 1, 5>;
        
        for (int i = 0; i < num_obs; ++i) {
            // 为单个观测点构造自动求导CostFunction
            CostFunc cost_func(&ResidualFunc, x_[i], y_[i]);
            
            // 计算残差
            double residual[1];
            cost_func.Evaluate(&params_[0], residual, nullptr);
            residuals_[i] = residual[0];
            
            // 自动计算雅可比（如果需要）
            if (evaluate_jacobians) {
                double jacobian[1][5];
                cost_func.Evaluate(&params_[0], nullptr, reinterpret_cast<double**>(jacobian));
                jacobians_.row(i) = Eigen::Map<Eigen::RowVectorXd>(jacobian[0], 5);
            }
        }
    }

    const double* params_;
    const Eigen::VectorXd x_;
    const Eigen::VectorXd y_;
    Eigen::VectorXd residuals_;
    Eigen::MatrixXd jacobians_;
    bool jacobians_are_stale_ = true;
};

关键修改说明 #

• 构造函数直接传入所有观测数据和优化参数引用
• 批量计算时对每个观测点使用AutoDiffCostFunction自动求导，无需手动推导雅可比
• 计算结果统一存储到residuals_和jacobians_中

3. 实现拷贝式CostFunction #

每个单点CostFunction从批量结果中拷贝数据：

class CopyingCostFunction : public ceres::CostFunction {
public:
    CopyingCostFunction(BatchEvaluationCallback* callback, int observation_index)
        : callback_(callback), obs_index_(observation_index) {
        set_num_residuals(1);
        mutable_parameter_block_sizes()->push_back(5);
    }

    bool Evaluate(double const* const* parameters,
                  double* residuals,
                  double** jacobians) const override {
        residuals[0] = callback_->residuals()[obs_index_];
        
        if (jacobians != nullptr && jacobians[0] != nullptr) {
            Eigen::Map<Eigen::RowVectorXd>(jacobians[0], 5) = callback_->jacobians().row(obs_index_);
        }
        return true;
    }

private:
    BatchEvaluationCallback* callback_;
    const int obs_index_;
};

4. 组装优化流程 #

将所有组件整合并运行优化：

int main() {
    // 替换为你的实际观测数据
    const int num_obs = 1000;
    Eigen::VectorXd x(num_obs), y(num_obs);
    // ... 填充x和y的数据 ...

    // 初始化优化参数初始值
    double params[5] = {1.0, 2.0, 3.0, 4.0, 5.0};

    // 创建批量回调实例
    BatchEvaluationCallback callback(params, x, y);

    // 构建优化问题
    ceres::Problem problem;
    for (int i = 0; i < num_obs; ++i) {
        ceres::CostFunction* cost_func = new CopyingCostFunction(&callback, i);
        problem.AddResidualBlock(cost_func, nullptr, params);
    }

    // 配置求解器并运行
    ceres::Solver::Options options;
    options.evaluation_callback = &callback;
    options.minimizer_progress_to_stdout = true;
    ceres::Solver::Summary summary;
    ceres::Solve(options, &problem, &summary);

    std::cout << summary.BriefReport() << std::endl;
    std::cout << "优化后参数: " << params[0] << ", " << params[1] << ", " << params[2] << ", " << params[3] << ", " << params[4] << std::endl;

    return 0;
}

效率优化建议 #

若观测点数量极大，可改用批量自动求导模板函数减少循环开销：

template<typename T>
bool BatchResidualFunc(const T* params, const Eigen::Matrix<T, Eigen::Dynamic, 1>& x, const Eigen::Matrix<T, Eigen::Dynamic, 1>& y, Eigen::Matrix<T, Eigen::Dynamic, 1>* residuals) {
    for (int i = 0; i < x.size(); ++i) {
        (*residuals)[i] = y[i] - YourNonlinearFunction(x[i], params[0], params[1], params[2], params[3], params[4]);
    }
    return true;
}

结合Eigen向量化操作，可一次性计算所有残差和雅可比，进一步提升批量计算效率。

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