本节主要依据Ceres官方指南翻译而成。
最小二乘法(least square)和非线性最小二乘分析的本来目的就是对一组数据进行曲线拟合。本节将介绍曲线拟合的问题。本节所用的采样点根据
生成,并且加入标准差为
高斯噪声。这
个数据,存入data[]
当中。下面我们用下列带未知参数的方程来拟合这些采样点:
同样的,我们定义一个用来计算残差的结构体。注意,对应每个采样点(观测点)都要计算一个残差。
struct ExponentialResidual {
ExponentialResidual(double x, double y)
: x_(x), y_(y) {}
//结构体的构造函数。把x赋值给x_,把y赋值给y_。也就是说在建立一个对象的时候x_和y_是赋完值的。
template <typename T>
bool operator()(const T* const m, const T* const c, T* residual) const {
//计算残差,观测值-理论值。
residual[0] = T(y_) - exp(m[0] * T(x_) + c[0]);
return true;
}
private:
// 一个样本数据的观测值。
const double x_;
const double y_;
};
下面构造Problem
。
double m = 0.0;
double c = 0.0; //初始值
Problem problem;
for (int i = 0; i < kNumObservations; ++i) {
CostFunction* cost_function =
new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>
new ExponentialResidual(data[2 * i], data[2 * i + 1]));
problem.AddResidualBlock(cost_function, NULL, &m, &c);
}
这里,不妨再与Hello World ( )对比一下:
struct CostFunctor { template <typename T> bool operator()(const T* const x, T* residual) const { residual[0] = T(10.0) - x[0]; return true; } };
CostFunction* cost_function = new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor); problem.AddResidualBlock(cost_function, NULL, &x);
通过对比,可以发现。在Hello World (HW)中,CostFunctor中是没有(显式)构造函数的,也就同样没有了初始值。所以在构造对象时,可以直接New CostFunctor
。而在本节的曲线拟合例子中,构造对象时还要加上初始值。
new ExponentialResidual(data[2 * i], data[2 * i + 1]));
在方括号中的参数分别对应残差函数名<ExponentialResidual,
,以及输出值(residual)的维度1,
,还有残差函数各个输入值(m,c)维度1,1>
。所以在本例中一共有三个1,而在HW中,只有两个1,即residual和x的维度。注意先是残差,后是输入参数,而且一一对应。
最后一点就是在把残差块加入problem的过程中,要把输入变量一一带入。比如&m,&c,&x等。以上就是在构建Problem的时候需要顾及到的三个方面。再就是在使用Numeric算法时,还要在方括号中指定计算机如何计算导数,如ceres::CENTRAL
。
最后运行程序(examples/curve_fitting.cc
),得到下列结果:
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.34e-04 2.56e-03
1 1.211734e+02 -2.21e+03 0.00e+00 7.52e-01 -1.87e+01 5.00e+03 1 4.29e-05 3.25e-03
2 1.211734e+02 -2.21e+03 0.00e+00 7.51e-01 -1.86e+01 1.25e+03 1 1.10e-05 3.28e-03
3 1.211734e+02 -2.19e+03 0.00e+00 7.48e-01 -1.85e+01 1.56e+02 1 1.41e-05 3.31e-03
4 1.211734e+02 -2.02e+03 0.00e+00 7.22e-01 -1.70e+01 9.77e+00 1 1.00e-05 3.34e-03
5 1.211734e+02 -7.34e+02 0.00e+00 5.78e-01 -6.32e+00 3.05e-01 1 1.00e-05 3.36e-03
6 3.306595e+01 8.81e+01 4.10e+02 3.18e-01 1.37e+00 9.16e-01 1 2.79e-05 3.41e-03
7 6.426770e+00 2.66e+01 1.81e+02 1.29e-01 1.10e+00 2.75e+00 1 2.10e-05 3.45e-03
8 3.344546e+00 3.08e+00 5.51e+01 3.05e-02 1.03e+00 8.24e+00 1 2.10e-05 3.48e-03
9 1.987485e+00 1.36e+00 2.33e+01 8.87e-02 9.94e-01 2.47e+01 1 2.10e-05 3.52e-03
10 1.211585e+00 7.76e-01 8.22e+00 1.05e-01 9.89e-01 7.42e+01 1 2.10e-05 3.56e-03
11 1.063265e+00 1.48e-01 1.44e+00 6.06e-02 9.97e-01 2.22e+02 1 2.60e-05 3.61e-03
12 1.056795e+00 6.47e-03 1.18e-01 1.47e-02 1.00e+00 6.67e+02 1 2.10e-05 3.64e-03
13 1.056751e+00 4.39e-05 3.79e-03 1.28e-03 1.00e+00 2.00e+03 1 2.10e-05 3.68e-03
Ceres Solver Report: Iterations: 13, Initial cost: 1.211734e+02, Final cost: 1.056751e+00, Termination: CONVERGENCE
Initial m: 0 c: 0
Final m: 0.291861 c: 0.131439
最终经过计算,结果是
。这个值和一开始的设定值(
)略有偏差。因为额外加入了高斯噪声,所以这个偏差的存在是合理的。拟合结果如下图所示,红圈即为拟合出的曲线。
鲁棒的曲线拟合Robust Curve Fitting
如果我们在数据集合中加入几个非常夸张的外围点Outliers,那么拟合结果会受到这几个点的明显影响。在这个时候需要应用损失函数(Loss function)来对异常数据进行过滤。比如在上文的例子中,我们对代码进行以下修改:
problem.AddResidualBlock(cost_function, NULL , &m, &c);
改为:
problem.AddResidualBlock(cost_function, new CauchyLoss(0.5) , &m, &c);
CauchyLoss是Ceres Solver附带的损失函数之一。 参数0.5指定了损失函数的规模。通过对下面两个图的对比,我们可以明显的发现损失函数的作用。