Ceres Solvel非线性最小二乘法建模

简介

Ceres由两个部分组成,一个是建模API,它提供了非常丰富的工具,可以迅速构建一个优化问题模型。另一个是解算器API,用于管控最小化算法。这一章将围绕如何用Ceres进行优化问题建模展开。

自动微分AutoDiffCostFunction

定义一个CostFunction(例如使用数值微分法或者解析微分法)可能是一个繁琐且容易出错的过程,尤其是在计算导数的时候。为此,Ceres提供了AutoDiffCostFunction

template <typename CostFunctor,
        int kNumResiduals,  // Number of residuals, or ceres::DYNAMIC.
        int N0,       // Number of parameters in block 0.
        int N1 = 0,   // Number of parameters in block 1.
        int N2 = 0,   // Number of parameters in block 2.
        int N3 = 0,   // Number of parameters in block 3.
        int N4 = 0,   // Number of parameters in block 4.
        int N5 = 0,   // Number of parameters in block 5.
        int N6 = 0,   // Number of parameters in block 6.
        int N7 = 0,   // Number of parameters in block 7.
        int N8 = 0,   // Number of parameters in block 8.
        int N9 = 0>   // Number of parameters in block 9.
class AutoDiffCostFunction : public
SizedCostFunction<kNumResiduals, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9> {
    public:
        explicit AutoDiffCostFunction(CostFunctor* functor);
        // Ignore the template parameter kNumResiduals and use
        // num_residuals instead.
        AutoDiffCostFunction(CostFunctor* functor, int num_residuals);
};

为了获得一个可以自动微分的代价函数,必须定义一个类。这个类中带有模板函数:operator(),它使用模板T类型进行代价函数运算。自动微分将根据需要用Jet类型替代替代模板T。但这个是隐藏的,编程的时候要把这个T看作一个双精度浮点数。这个函数必须把计算结果以最后一个参数(唯一一个非常量参数)传递出来,并且返回True,告诉计算机运算成功完成。例如,现在有一个标量的偏差函数$e = k - x^Ty$。这里$x$和$y$都是二维向量参数,$k$是个常量 。这种类型的偏差,即一个常量和一个表达式的差值,在最小二乘法问题中很常见。例如,$x ^ T y$可能是一系列测量结果的期望值,那么每一次测量$K$都对应了一个代价函数类的实例。被加到Problem中的是$e ^ 2$或者$(k - x^Ty) ^ 2$。平方处理由Ceres优化框架完成。这个例子的具体代码如下:

这里有个疑问,之前自动微分算法的代码中没有使用类,而是使用的struct,不知道为什么这里是class,个人认为应该都可以,因为我们的代码中都使用的struct。

class MyScalarCostFunctor {
    MyScalarCostFunctor(double k): k_(k) {}

template <typename T>
bool operator()(const T* const x , const T* const y, T* e) const {
    e[0] = k_ - x[0] * y[0] - x[1] * y[1];
    return true;
}

private:
    double k_;
};

注意,在operator()的声明中,首先是输入参数,他们都是指向T类型数组的常指针。如果由更多的输入参数就跟在y后面。而输出值永远是最后一个参数,并且也是一个指向数组的指针。在上述例子中,e是标量,所以只赋值e[0]

然后给出这个类的定义,它的自动微分代价函数可以如下构造:

CostFunction* cost_function
    = new AutoDiffCostFunction<MyScalarCostFunctor, 1, 2, 2>(
        new MyScalarCostFunctor(1.0));              ^  ^  ^
                                                    |  |  |
                        Dimension of residual ------+  |  |
                        Dimension of x ----------------+  |
                        Dimension of y -------------------+

在这个例子中,对每次测量k都有一个实例。模板参数1,2,2将Functor描述为一个一维输出参数和两个二维输入参数。AutoDiffCostFunction也支持在运行时动态确定参数个数。例如下面的代码:

CostFunction* cost_function
    = new AutoDiffCostFunction<MyScalarCostFunctor, DYNAMIC,  2,  2>(
        new CostFunctorWithDynamicNumResiduals(1.0),    ^     ^  ^
        runtime_number_of_residuals);  <----+           |     |  |
                                            |           |     |  |
                                            |           |     |  |
           Actual number of residuals ------+           |     |  |
           Indicate dynamic number of residuals --------+     |  |
           Dimension of x ------------------------------------+  |
           Dimension of y ---------------------------------------+

Ceres目前支持代价函数最多有10个相互独立的变量,但是对每个变量有多少维度没有限制。

注意,新用户常常犯的一个错误就是把模板参数中的数字理解成参数的个数。但事实上,模板参数中数字的含义是每个参数的维度。这两个概念不能混淆。比如在这个例子中x,y都是二维变量,所以模板参数中有两个2。

Problem的构建

Problem保持了非线性最小二乘问题的强化的边界。要创建最小二乘问题,可以使用Problem::AddResidualBlock()Problem::AddParameterBlock()。例如,下面这个Problem包含了三个参数块,维度分别为3,4,5。同时有两个残差块,维度分别是2和6。

double x1[] = { 1.0, 2.0, 3.0 };
double x2[] = { 1.0, 2.0, 3.0, 5.0 };
double x3[] = { 1.0, 2.0, 3.0, 6.0, 7.0 };

Problem problem;
problem.AddResidualBlock(new MyUnaryCostFunction(...), x1);
problem.AddResidualBlock(new MyBinaryCostFunction(...), x2, x3);

Problem::AddResidualBlock(),顾名思义,就是把残差块加入到Problem中。它添加了一个CostFunction,一个LossFunction(非必要)并且把CostFunction链接到一系列的参数块上。代价函数包含了关于期望的参数块大小的信息。该函数检查他们是否与parameter_blocks一致。如果不匹配,程序将终止。LossFunction可以是Null,这种情况下这一项的代价就是残差的平方。

另外官网中给出了很多Problem的方法,可以查看一下。这里只对添加残差块的方法给出了实例。