How to use AutoDiffCostFunction within EvaluationCallback with Ceres-solver

Question

I am a newbie trying to understand how to use ceres-solver for my problem, and was wondering if someone could point to an example etc.. that may be readily adapted to my use-case. My problem is as follows;

I wish to fit a non-linear, parameterised curve f(x; a_1, a_2, a_3, a_4, a_5) to an observed dataset y(x), where a_1...a_5 are parameters to be estimated. The curve is mathematically complicated, and difficult to determine closed-form expressions for the derivatives, so am hoping to use auto-derivatives. I am able to implement f(x; a_1, a_2, a_3, a_4, a_5) with templates.

There are several thousand observations (eg. values of y(x) for x in [-2000,...2000] ) with which to fit f(x; a_1, a_2, a_3, a_4, a_5). It is much more efficient to evaluate the function for all values of x at the same time, rather than evaluating individual points in isolation.

As I currently understand (after first reading of the documentation), the best way forward would be similar to ceres-solver-directory/examples/evaluation_callback_example.cc. ie;

• Implement a class derived from EvaluationCallback (to efficiently evaluate the residuals/Jacobian of the function at all values of 'x' at the same time)
• Implement a 'CostAndJacobianCopyingCostFunction' to return residuals and Jacobians evaluated at individual values of 'x'.

The part I am not quite sure of is;

Q: How to modify the evaluation_callback_example so that ceres::AutoDiffCostFunction is used in the EvaluationCallback function to avoid calculating derivatives/Jacobian by hand? (The evaluation_callback_example uses a closed form expression for the Jacobians, but my function is too complicated for that...)

For convenience, here is a copy of the EvaluationCallback function in evaluation_callback_example

// 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;
};

Answer 1

Update: Does not exactly answer my original question, but I think I will be able to solve my problem. I have reworked the curve_fitting.cc example code so it was structured similar to my needs (ie. to use autodiff rather than calculating Jacobians in closed form, and calculating all values at the same time).

Can anyone see issues with this restructure that I am oblivious to?

const int kNUMOBSERVATIONS = 67;

// Separated interleaved x,y data into individual arrays for convenience.
const double x_data[] = {
  0.000000e+00, 
  7.500000e-02, 
  1.500000e-01, 
  ...etc...
};
const double y_data[] = {
  1.133898e+00,
  1.334902e+00,
  1.213546e+00,
  ...etc...
};


struct ExponentialResidual 
{
  ExponentialResidual(double* x, double* y)  
  {
      x_ =  Eigen::ArrayXd::Zero(N_OBSERVATIONS);
      y_ =  Eigen::ArrayXd::Zero(N_OBSERVATIONS);
      
      for (int ii = 0; ii < N_OBSERVATIONS; ii++) 
      {
          x_[ii] = x[ii];
          y_[ii] = y[ii];
      }
  }

  // Calculate all residuals at once as a vector.
  // This allows Jacobians to be calculated using Autodiff.
  template <typename T>
  bool operator()(const T* const m, const T* const c, T* residual) const {
      Eigen::Array<T,kNUMOBSERVATIONS,1> temp  = y_ - exp(m[0] * x_ + c[0]);
      for (int ii = 0; ii < kNUMOBSERVATIONS; ii++) 
        residual[ii] = temp(ii);
    return true;
  }

 private:  
  Eigen::ArrayXd x_;
  Eigen::ArrayXd y_;
};

int main(int argc, char** argv) {
  google::InitGoogleLogging(argv[0]);

  const double initial_m = 0.0;
  const double initial_c = 0.0;
  double m = initial_m;
  double c = initial_c;

  ceres::Problem problem;
  // Rather than lots of one-dimensional residual blocks,
  // we have one vector residual block.
  // This allows the residuals to be calculated in one place. 
  problem.AddResidualBlock(
    new ceres::AutoDiffCostFunction<ExponentialResidual, kNUMOBSERVATIONS, 1, 1>(
        new ExponentialResidual((double*) x_data, (double*) y_data)),
    nullptr,
    &m,
    &c);
  

  ceres::Solver::Options options;
  options.max_num_iterations = 25;
  options.linear_solver_type = ceres::DENSE_QR;
  options.minimizer_progress_to_stdout = true;

  ceres::Solver::Summary summary;
  ceres::Solve(options, &problem, &summary);
  std::cout << summary.BriefReport() << "\n";
  std::cout << "Initial m: " << initial_m << " c: " << initial_c << "\n";
  std::cout << "Final   m: " << m << " c: " << c << "\n";
  return 0;
}