用Ceres-Solver求解非线性最小二乘问题

简介

Ceres-Solver是谷歌推出的开源C++库，用于求解优化问题。本文主要介绍用Ceres求解优化问题中最常见的最小二乘问题。例程来源于官网中Tutorial: Non-linear Least Squares章节, Ceres版本为1.14.

最小二乘问题

从官方文档看，Ceres可以求解如下问题：

$$\begin{gathered} x^{\ast }=\underset{x}{\mathrm{argmin}}\sum _i{\rho }_i(||f_i(x_{i_1},x_{i_2},...,x_{i_k})|{|}^2) \\ s.t.l_j\le x_j\le u_j \end{gathered}$$

传统的优化问题可以描述为：选择一组参数（变量），在满足一系列有关的限制条件（约束）下，使设计指标（目标）达到最优值。据此来分析上述公式：

变量

ceres可以优化多个变量，式中的$x$就是待优化变量(下文简称变量或自变量)的集合。再次强调，$x$是变量的集合，每个变量都是向量，它们的维度可以不同。

约束

公式中仅仅对单个变量的上下限进行了约束，这看上去并不完备，因为传统优化问题可能会对多个变量做整体约束。但是在官方的Tutorial: Non-linear Least Squares章节中，约束不作为讨论的重点。

目标

从式中可以看出，待优化的目标是是多个项目的累加。沿用Ceres的命名方式，我们称每个项为"残差项"。

在Ceres中，残差项是处理的最小单元。每个残差项可以写成:

$$R_i={\rho }_i(||f_i(x_i)|{|}^2)$$

其中
$$x_i=\{x_{i_1},x_{i_2},...,x_{i_k}\},x_i\in x$$

可以看出，残差项具有如下形式：代价函数$f$作用于自变量上，对结果求二范数平方，然后用损失函数$\rho$来处理，得到最终结果。

值得注意的是，式中每个残差项都用$i$作为下标，这暗示这每个项都有自己独立的参数，包括：

• 变量：每个残差项都有自己的自变量$x_i$，它是多个变量的集合，是全体待优化变量$x$的子集。并且，各个残差项的$x_i$的元素个数可以不同。换言之，所有残差项的自变量构成了整体自变量集合$x$。

• 代价函数：代价函数作用于变量（以及可能的常量）上。从式中的二范数可以看出，代价函数返回的是一个向量。

• 损失函数：通常用于减少外点对模型拟合的影响，此处不去深究数学特性，仅仅认为是一个函数即可。

流程示例

接下来以官方例程说明Ceres整体流程和结构.

求解问题:

$$\underset{x}{\mathrm{argmin}}(10-x{)}^2,x\in R$$

1. 要素分析

首先明确待优化变量，分解各个残差项，并剥离出每个残差项的待优化变量、代价函数和损失函数。这里待优化变量是$x$，只有一个残差项。

残差项具有"损失函数作用在代价函数的二范数平方上"的运算形式，而二范数通常是一项或多项的平方和，以此为特征，易得代价函数$f(x)=10-x$，损失函数为空，同时存在常数项10。

2. 定义代价函数类

代价函数类通过重载括号运算符，以实现代价函数$f(\cdot )$的运算逻辑。代价函数以自变量作为输入，输出一个向量作为结果。代码中函数的形参列表声明了自变量和输出向量，函数的主体实现了核心计算逻辑。

struct CostFunctor {
   template <typename T>
   bool operator()(const T* const x, T* residual) const {
     residual[0] = T(10.0) - x[0];
     return true;
   }
};

注意，本例中常量10被硬编码到函数主体中，这样缺乏灵活性。后续会讲到如何从外部传入常量。

3. 求解过程调用

调用时包括步骤：

1. 初始化GLog。实际测试本步骤可省。
2. 创建Problem对象，并添加各个残差项。
3. 配置求解器选项并求解。

其中，添加残差项的过程包括：

• 指定(我们自己定义的)代价函数类
• 分别指定代价函数中输入向量与输出向量的维度，其中输入向量可以有多个
• 构造代价函数类的实例，并代价函数的常量(如果有)
• 绑定自变量实体，用以传入初值并传出结果

int main(int argc, char** argv) 
{
  google::InitGoogleLogging(argv[0]);
  // The variable to solve for with its initial value.
  double x = 5.0;

  // Build the problem.
  Problem problem;             
                                  
  CostFunction* cost_function =   
      new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
  // 此处指定了代价函数类名，各个变量的维度，函数输出的维度，并构造代价函数实例
  problem.AddResidualBlock(cost_function, NULL, &x);
  // 指定损失函数，绑定自变量

  // 进行求解
  Solver::Options options;
  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 << "x: "<< x << "\n";
  return 0;
}

总结与要点分析

问题分析与要素提取

对具体问题明确自变量有哪些，划分各个残差项，对每个残差项提取它的自变量、代价函数和损失函数以及常量(如果有)。由于每个残差项都具备"损失函数作用于代价函数的二范数平方上"的形式，而二范数平方通常是一项或多项平方和的形式，以此提取代价函数以及其他各个要素。

初始化操作

初始化操作包括：

1. 初始化glog(实测可省)
2. 创建自变量，并赋予初值
3. 创建Problem对象

  google::InitGoogleLogging(argv[0]);
  double x = 5.0;
  Problem problem;    

然后就是为各个残差项定义代价函数类，并添加到Problem对象中，再进行求解了。

为残差项定义代价函数类

代价函数类实现了代价函数的计算逻辑。它以自变量为输入，输出一个向量。自变量可以是多个，维度可以不同。注意此处没有指定自变量和输出向量的维度。

格式

类名可以自定义，核心函数重载了括号运算符，内部定义了代价函数运算规则。类可以有其它方法或成员，并没有严格要求。

struct MyFunc{
    template <typename T>
    bool operator()(const T* const x1, const T* const x2, ..., T* residual) const {
     residual[0] = T(10.0) - x1[0] + ceres::sin(x2[0]);
     ...
     return true;
   }
};

注意重载函数的格式，不要漏掉尾部const。参数类型也要严格按照范例。参数列表中前面的形参表示输入的自变量(可以是多个)，最后形参表示输出向量。内部实现了输入到输出的运算逻辑。

注意运算的对象都是向量。标量也视为长度为1的向量，别忘了下标。

常数项

实际问题中可能代价函数需要常数项参与运算，此时常数项可以作为类的成员被函数访问，并通过构造函数来从外界传入。常数项的数目和类型等没有强制规范。

struct MyFunc{
	MyFunc(double x): x(x) {}
	template<typename T> bool operator()(....){...} [核心函数]
	private:
		const double x;
};

模板类型

函数主体拥有模板类型T. 代价函数的输入变量和输出量都是T类型的向量，此处的T可以按需实例化为double或Jet(暂不明)类型。注意T不能是用户自己定义的其他数据类型。

根据C++模板函数的规则，T类型会根据实际调用时的情况进行实例化。如果应用场合是针对double的，那么T只会实例化为double。然而笔者建议尽量把操作的数据类型都泛化为T，以提高通用性。比如用ceres::cos替代cos，因为后者仅处理double类型；常量也要强制转换为T，等等。

实际上Ceres有不同的求解方式，分别具有不同的函数形式。常见的类型以及它们的核心函数的参数列表如下:

AutoDiffCostFunction: (const T* const x, T* residual) 
NumericDiffCostFunction : (const double* const x, double* residual)
QuadraticCostFunction: (double const* const* parameters, double* residuals, double** jacobians)

Auto与Numeric相比较来说，官方文档推荐使用Auto，称它有更好的表现。
Quadratic需要手动指定雅克比，即指定了求导规则(而Auto和Numeric可以自动求导)，官方称除非有特定理由，否则不推荐。此外它重载的是名为Evaluate的方法而非括号运算符，类的声明更复杂。本文档不做重点研究。

把残差项添加到问题中

代价函数类定义了代价函数的运算逻辑，而添加残差项的操作实现了相关的处理动作。

Problem problem;
CostFunction* cost_function =
      new AutoDiffCostFunction<MyFunc, output_size, x1_size, x2_size, ... >(new MyFunc);
problem.AddResidualBlock(cost_function, LossFunc, &x1, &x2, ...);

创建Problem对象，对Problem对象调用AddResidualBlock方法，添加每个残差项。每个残差项指定了自变量、代价函数、损失函数(以及可能的常数项)。从程序实现的角度，过程有：

1. 构造代价函数MyFunc实体，如果有需要可以传入常量。当然，注意C++中构造对象时的写法，在无参时是不需要括号的。
2. 根据MyFunc实例创建CostFunction对象指针，主体是AutoDiffCostFunction类型，尖括号内指定了代价函数的类名、各个自变量维度、输出向量维度。
3. 对Problem对象调用AddResidualBlock方法，参数包括上文创建的CostFunction指针和损失函数，并绑定各个自变量, 用来传入变量初值并传出结果。

配置选项并求解

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

ceres::Solver::Summary summary;                //Optimization
Solve(options, &problem, &summary);            //solve problem
std::cout << summary.BriefReport() << "\n";

此处仅仅是常见的选项，具体查阅官方文档。

范例解析

例：各个残差项配置不同的情况

求解：

$$\underset{x}{\mathrm{argmin}}\sum _i||F_i(x)|{|}^2$$

其中
$$\begin{gathered} F_1(x)=x_1+10x_2 \\ {F}_2(x)=\sqrt{10}(x_2-x_3) \end{gathered}$$

分析：此处自变量是集合$\{x_1,x_2,x_3\}$，残差项有两项。首先声明自变量并给定初值：

double x1 =  1.0, x2 = 2.0, x3 = 3.0;

对各个残差项分别处理。以下仅给出了核心函数参数列表和主要操作代码。

残差项1：自变量为$\{x_1,x_2\}$,常量10，代价函数为$F_1(\cdot )$

struct F1: (const T* const x1, const T* const x2, T* residual){
                residual[0] = x1[0] + T(10.0) * x2[0];
            }
problem.AddResidualBlock(
  new AutoDiffCostFunction<F1, 1, 1, 1>(new F1), NULL, &x1, &x2);

残差项2：自变量为$\{x_2,x_3\}$，代价函数为$F_2(\cdot )$

struct F2: (const T* const x2, const T* const x3, T* residual){
                residual[0] = T(sqrt(10.0)) * (x2[0] - x3[0]) ;
            }
problem.AddResidualBlock(
  new AutoDiffCostFunction<F2, 1, 1, 1>(new F2), NULL, &x2, &x3);

注意代码中对根号10进行转换为T。

后续配置参数并求解即可。

当然，这里说一句，不同残差块的代价函数不一样的情况不常见。常见的如下文中的模型拟合，它们的代价函数是一样的，仅常量不一样。

例：多维向量的情况

上述问题也可以写成:

$$\underset{x}{\mathrm{argmin}}(x_1+10x_2{)}^2+(\sqrt{10}(x2-x3){)}^2$$

等价

$$\underset{x}{\mathrm{argmin}}||f(x)|{|}^2$$

其中
$$f(x)=\left[\begin{array}{c}{x}_1+10x_2\\ \sqrt{10}(x2-x3)\end{array}\right],x=[x_1,x_2,x_3{]}^T$$

也就是说，自变量只有一个三维向量，输出是一个二维向量。

struct MyFunc{
	template <typename T> bool operator()(const T* const x, T* residual) const {
		residual[0] = x[0] + T(10.0) * x[1];
		residual[1] = T(sqrt(10.0)) * (x[1] - x[2]);
		return true; 
	}
};
double x[3] = {1.0, 2.0, 3.0};
Problem problem;
problem.AddResidualBlock(
  new AutoDiffCostFunction<MyFunc, 2, 3>(new MyFunc), NULL, x);

注意C++语法特性，例如下标从零开始，数组名表示首元素或者整个数组的地址，等等。

例：模型拟合

这个基本上是最常见的情况了：拟合模型
$$y=e^{mx+c}$$

已经有了多组数据$\{x_i,y_i\}$，尝试计算

$$\underset{m,c}{\mathrm{argmin}}\sum _i(y_i-e^{m{x}_i+c}{)}^2$$

可见：自变量是$\{m,c\}$，多个残差项，它们的自变量和代价函数、损失函数都一致，但常量不一样。因此只需定义一个包含常量的代价函数类，然后实例化时传入不同常量即可。

struct F {
  F(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] = T(y_) - exp(m[0] * T(x_) + c[0]);
    return true;
  }
 private:
  const double x_, y_;
};

这里的书写很严谨:常量都是私有const类型，实际猜测应该没那么严格。

疑问: 里面的exp没用ceres::，是因为没有么？还是说这里的实际场合就是double不是别的，不用转换？

double m = 0.0, c = 0.0;
Problem problem;
for (int i = 0; i < N; ++i) {
  CostFunction* cost_function =
       new AutoDiffCostFunction<F, 1, 1, 1>(new F(x[i], y[i]));
  problem.AddResidualBlock(cost_function, NULL, &m, &c);
}

例: 损失函数的使用

损失函数涉及到优化问题的内点和外点的概念。外点可以认为是训练数据集中错误的条目。错误的数据会对结果产生影响，而损失函数可以降低这种影响。具体代码中写为：

problem.AddResidualBlock(cost_function, new CauchyLoss(0.5) , &m, &c);

损失函数的数学含义此处不做讨论。

总结

Ceres作为求解优化问题的利器，此处仅做入门级讨论，不过应该够用了。