Ceres曲线拟合

最新推荐文章于 2023-04-02 19:48:44 发布

唐僧洗头用飘柔dp

最新推荐文章于 2023-04-02 19:48:44 发布

阅读量130

点赞数 1

分类专栏： Ceres Solver 教程文章标签：算法 c++

本文链接：https://blog.csdn.net/weixin_45860565/article/details/126541333

版权

Ceres Solver 教程专栏收录该内容

4 篇文章 0 订阅

订阅专栏

Ceres Solver 曲线拟合

$\quad$ 到目前为止，我们看到的示例都是没有数据的简单优化问题。最小二乘和非线性最小二乘分析的最初目的是拟合数据曲线。

$\quad$ 这次拟合的曲线为 $y=e^{0.3x+0.1}$ 并且加入了 $0.2$ 的高斯噪声。

$\quad$ 我们首先定义一个模板化对象来评估残差。每次观察都会有一个残差。

struct ExponentialResidual {
  ExponentialResidual(double x, double y)
      : x_(x), y_(y) {}

  template <typename T>
  bool operator()(const T* const m, const T* const c, T* residual) const {
    residual[0] = y_ - exp(m[0] * x_ + c[0]);
    return true;
  }

 private:
  // Observations for a sample.
  const double x_;
  const double y_;
};

$\quad$ 构建残差块:

double m = 0.0;
  double c = 0.0;

  Problem problem;
  for (int i = 0; i < kNumObservations; ++i) {
    problem.AddResidualBlock(
        new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
            new ExponentialResidual(data[2 * i], data[2 * i + 1])),
        nullptr,
        &m,
        &c);
  }

$\quad$ 这里说一下个人理解：

$\quad$ 首先本次出现的问题和前面的两个问题略有不同，前面两个问题是没有观测这一说的。无论是在计算 $10 - x$ 还是在计算 Powell's Function 时，其表达式都是明确知道的。但是在这里却不是的，我们要估计的就是表达式的 $m$ 和 $c$ 的值。

$\quad$ 在第一个例子中需要计算的是 $10-x)^2$ 的最小值，在这个例子中是只有一个残差块的，因此在添加残差块的时候也只添加了一次，并且没有所谓的观测值，调用 Ceres 的自动求导，找到其梯度下降的方向，然后进行优化即可。

$\quad$ 在第二个例子中优化的是 Powell's Function 的最小值，其中 $x$ 是一个多维变量，在只有一个约束的情况下是没办法找到最小值的，因此该问题中也是有四个约束，则分别构成了四个残差块，但是也没有所谓的观测值。

$\quad$ 在这个问题中我们有大量的 $x_i,y_i]$ 每个观测可以构成一个残差约束块：
$residual_i=y_i-e^{mx_i+c}$
$\quad$ 笔者在最开始的时候没能理解为什么有的需要观测值，有的则不需要。后来仔细想了想也算是理解了。从这三个问题出发，首先明确一点，有多少个观测就会有多少个残差块，并且残差块的个数是不能少于变量的维度的。

$\quad$ 例 $1$ 不需要观测是因为本身就求的是该函数的最小值，就是求个导数即可，可以将其想象成在拥有观测问题中的某一次观测，变量为 $x$ ，其他部分常数部分则为观测值。

$\quad$ 例 $2$ 中不需要观测则可以这样理解，四个方程构成了四个观测，每个观测对其中部分变量进行了约束。

$\quad$ 在曲线拟合的例子中， $x_i,y_i]$ 均为观测值，为常数，每一对构成了一个方程，其中的变量，或者说待优化变量则为 $m ， c$ ，求导数的对象就是这两个参数才对。

$\quad$ 最终优化结果如下，与真实值的接近程度和噪声的大小和数据量的大小都有一定的关系。

iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time
   0  1.211734e+02    0.00e+00    3.61e+02   0.00e+00   0.00e+00  1.00e+04        0    5.57e-04    2.30e-03
   1  2.334822e+03   -2.21e+03    0.00e+00   7.52e-01  -1.87e+01  5.00e+03        1    7.49e-05    3.14e-03
   2  2.331438e+03   -2.21e+03    0.00e+00   7.51e-01  -1.86e+01  1.25e+03        1    2.20e-05    3.17e-03
   3  2.311313e+03   -2.19e+03    0.00e+00   7.48e-01  -1.85e+01  1.56e+02        1    2.48e-05    3.21e-03
   4  2.137268e+03   -2.02e+03    0.00e+00   7.22e-01  -1.70e+01  9.77e+00        1    1.83e-05    3.23e-03
   5  8.553131e+02   -7.34e+02    0.00e+00   5.78e-01  -6.32e+00  3.05e-01        1    2.12e-05    3.27e-03
   6  3.306595e+01    8.81e+01    4.10e+02   3.18e-01   1.37e+00  9.16e-01        1    6.54e-04    3.93e-03
   7  6.426770e+00    2.66e+01    1.81e+02   1.29e-01   1.10e+00  2.75e+00        1    5.70e-04    4.52e-03
   8  3.344546e+00    3.08e+00    5.51e+01   3.05e-02   1.03e+00  8.24e+00        1    5.57e-04    5.09e-03
   9  1.987485e+00    1.36e+00    2.33e+01   8.87e-02   9.94e-01  2.47e+01        1    5.52e-04    5.66e-03
  10  1.211585e+00    7.76e-01    8.22e+00   1.05e-01   9.89e-01  7.42e+01        1    6.53e-04    6.33e-03
  11  1.063265e+00    1.48e-01    1.44e+00   6.06e-02   9.97e-01  2.22e+02        1    5.98e-04    6.95e-03
  12  1.056795e+00    6.47e-03    1.18e-01   1.47e-02   1.00e+00  6.67e+02        1    5.50e-04    7.51e-03
  13  1.056751e+00    4.39e-05    3.79e-03   1.28e-03   1.00e+00  2.00e+03        1    5.52e-04    8.08e-03
Ceres Solver Report: Iterations: 14, Initial cost: 1.211734e+02, Final cost: 1.056751e+00, Termination: CONVERGENCE
Initial m: 0 c: 0
Final   m: 0.291861 c: 0.131439

$\quad$ 完整代码如下：

#include "ceres/ceres.h"
#include "glog/logging.h"

using ceres::AutoDiffCostFunction;
using ceres::CostFunction;
using ceres::Problem;
using ceres::Solve;
using ceres::Solver;

// Data generated using the following octave code.
//   randn('seed', 23497);
//   m = 0.3;
//   c = 0.1;
//   x=[0:0.075:5];
//   y = exp(m * x + c);
//   noise = randn(size(x)) * 0.2;
//   y_observed = y + noise;
//   data = [x', y_observed'];

const int kNumObservations = 67;
// clang-format off
const double data[] = {
  0.000000e+00, 1.133898e+00,
  7.500000e-02, 1.334902e+00,
  1.500000e-01, 1.213546e+00,
  2.250000e-01, 1.252016e+00,
  3.000000e-01, 1.392265e+00,
  3.750000e-01, 1.314458e+00,
  4.500000e-01, 1.472541e+00,
  5.250000e-01, 1.536218e+00,
  6.000000e-01, 1.355679e+00,
  6.750000e-01, 1.463566e+00,
  7.500000e-01, 1.490201e+00,
  8.250000e-01, 1.658699e+00,
  9.000000e-01, 1.067574e+00,
  9.750000e-01, 1.464629e+00,
  1.050000e+00, 1.402653e+00,
  1.125000e+00, 1.713141e+00,
  1.200000e+00, 1.527021e+00,
  1.275000e+00, 1.702632e+00,
  1.350000e+00, 1.423899e+00,
  1.425000e+00, 1.543078e+00,
  1.500000e+00, 1.664015e+00,
  1.575000e+00, 1.732484e+00,
  1.650000e+00, 1.543296e+00,
  1.725000e+00, 1.959523e+00,
  1.800000e+00, 1.685132e+00,
  1.875000e+00, 1.951791e+00,
  1.950000e+00, 2.095346e+00,
  2.025000e+00, 2.361460e+00,
  2.100000e+00, 2.169119e+00,
  2.175000e+00, 2.061745e+00,
  2.250000e+00, 2.178641e+00,
  2.325000e+00, 2.104346e+00,
  2.400000e+00, 2.584470e+00,
  2.475000e+00, 1.914158e+00,
  2.550000e+00, 2.368375e+00,
  2.625000e+00, 2.686125e+00,
  2.700000e+00, 2.712395e+00,
  2.775000e+00, 2.499511e+00,
  2.850000e+00, 2.558897e+00,
  2.925000e+00, 2.309154e+00,
  3.000000e+00, 2.869503e+00,
  3.075000e+00, 3.116645e+00,
  3.150000e+00, 3.094907e+00,
  3.225000e+00, 2.471759e+00,
  3.300000e+00, 3.017131e+00,
  3.375000e+00, 3.232381e+00,
  3.450000e+00, 2.944596e+00,
  3.525000e+00, 3.385343e+00,
  3.600000e+00, 3.199826e+00,
  3.675000e+00, 3.423039e+00,
  3.750000e+00, 3.621552e+00,
  3.825000e+00, 3.559255e+00,
  3.900000e+00, 3.530713e+00,
  3.975000e+00, 3.561766e+00,
  4.050000e+00, 3.544574e+00,
  4.125000e+00, 3.867945e+00,
  4.200000e+00, 4.049776e+00,
  4.275000e+00, 3.885601e+00,
  4.350000e+00, 4.110505e+00,
  4.425000e+00, 4.345320e+00,
  4.500000e+00, 4.161241e+00,
  4.575000e+00, 4.363407e+00,
  4.650000e+00, 4.161576e+00,
  4.725000e+00, 4.619728e+00,
  4.800000e+00, 4.737410e+00,
  4.875000e+00, 4.727863e+00,
  4.950000e+00, 4.669206e+00,
};
// clang-format on

struct ExponentialResidual {
  ExponentialResidual(double x, double y) : x_(x), y_(y) {}

  template <typename T>
  bool operator()(const T* const m, const T* const c, T* residual) const {
    residual[0] = y_ - exp(m[0] * x_ + c[0]);
    return true;
  }

 private:
  const double x_;
  const double y_;
};

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

  double m = 0.0;
  double c = 0.0;

  Problem problem;
  for (int i = 0; i < kNumObservations; ++i) {
    problem.AddResidualBlock(
        new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
            new ExponentialResidual(data[2 * i], data[2 * i + 1])),
        nullptr,
        &m,
        &c);
  }

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

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