ceres 实现x=p*3+m的最小二乘估计

假设有n个数据点，其中第i个数据点的观测值为yi，设计矩阵为Xi，参数向量为x，残差向量为ei。那么，可以将x=p*3+m表示为如下形式：

f(x) = x1 * p^3 + x2 * m - y

其中，x = [p, m]，y = [y1, y2, ..., yn]，f(x)和y均为n维向量。

将f(x)进行最小二乘优化，得到残差平方和：

S = e1^2 + e2^2 + ... + en^2
  = (x1 * p^3 + x2 * m - y1)^2 + (x1 * p^3 + x2 * m - y2)^2 + ... + (x1 * p^3 + x2 * m - yn)^2

使用Ceres实现的代码如下：

#include "ceres/ceres.h"

// 自定义残差函数
struct MyResidual {
  MyResidual(double yi, double *Xi) : yi(yi), Xi(Xi) {}

  template <typename T>
  bool operator()(const T *const x, T *residual) const {
    residual[0] = T(Xi[0]) * ceres::pow(x[0], T(3)) + T(Xi[1]) * x[1] - T(yi);
    return true;
  }

  double yi;
  double *Xi;
};

int main() {
  // 数据点
  double y[] = {1.0, 2.0, 3.0, 4.0, 5.0};
  double X[][2] = {{1.0, 1.0}, {2.0, 1.0}, {3.0, 1.0}, {4.0, 1.0}, {5.0, 1.0}};

  // 初始化问题
  ceres::Problem problem;
  for (int i = 0; i < 5; ++i) {
    ceres::CostFunction *cost_function =
        new ceres::AutoDiffCostFunction<MyResidual, 1, 2>(new MyResidual(y[i], X[i]));
    problem.AddResidualBlock(cost_function, nullptr, x);
  }

  // 配置求解器
  ceres::Solver::Options options;
  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.FullReport() << std::endl;
  std::cout << "p: " << x[0] << ", m: " << x[1] << std::endl;

  return 0;
}

其中，MyResidual为自定义的残差函数，用于计算每个数据点的残差。在main函数中，首先构造问题，然后配置求解器，最后调用ceres::Solve函数进行求解。