1. 安装

略

2. 问题定义

以 $y = \exp \left( {a{x}^{2} + {bx} + c}\right)$ 为例

我们现在拥有一系列的观测值

x $x_1$ $x_2$ … $x_n$
y $y_1$ $y_2$ … $y_n$

目标是用函数 $y = \exp(ax^2 + bx + c)$ 拟合这些数据点，确定参数 $a$、$b$ 和 $c$ 的最佳值

在迭代的某一步，假设我们有参数的当前估计值 $a_k$、$b_k$ 和 $c_k$。

对于每个数据点 $x_i$，我们可以计算模型预测值 ：$\hat{y}_i$

$\hat{y}_i = \exp(a_k x_i^2 + b_k x_i + c_k)$

那么残差就是实际观测值与预测值的差

$r_i = y_i - \hat{y}_i = y_i - \exp(a_k x_i^2 + b_k x_i + c_k)$

我们不希望正负残差相互抵消，所以计算残差的平方，然后求和得到代价函数：

$J(a, b, c) = \sum_{i=1}^{n} r_i^2 = \sum_{i=1}^{n} [y_i - \exp(a x_i^2 + b x_i + c)]^2$

最小二乘法的目标是找到参数 $a$、$b$ 和 $c$ 的值，使得代价函数 $J(a, b, c)$ 最小。即：

$\min_{a, b, c} J(a, b, c) = \min_{a, b, c} \sum_{i=1}^{n} [y_i - \exp(a x_i^2 + b x_i + c)]^2$

3. 三种求解方式

3.1 方法一：解析求导

这种方法需要我们手动计算雅可比矩阵，并将其提供给 Ceres。

3.1.1 计算雅可比矩阵

我们有残差 $r(a, b, c) = y - \exp(ax^2 + bx + c)$
令 $E = \exp(ax^2 + bx + c)$
使用链式法则求偏导：

  • $\frac{\partial r}{\partial a} = - \frac{\partial E}{\partial a} = - E \cdot \frac{\partial (ax^2 + bx + c)}{\partial a} = - E \cdot x^2$
  • $\frac{\partial r}{\partial b} = - \frac{\partial E}{\partial b} = - E \cdot \frac{\partial (ax^2 + bx + c)}{\partial b} = - E \cdot x$
  • $\frac{\partial r}{\partial c} = - \frac{\partial E}{\partial c} = - E \cdot \frac{\partial (ax^2 + bx + c)}{\partial c} = - E \cdot 1$

3.1.2 实现 Ceres Cost Function

在Ceres Solver中实现解析求导，我们需要继承的是ceres::SizedCostFunction类。

SizedCostFunction是一个模板类，参数格式如下：

1
2
ceres::SizedCostFunction<残差数量, 参数块1大小, 参数块2大小, ...>

对于这个问题（拟合模型 y = exp(ax² + bx + c)）：

  • 残差数量 = 1（每个观测点产生1个残差值）
  • 参数块大小 = 3（我们优化的参数是a, b, c）

所以我们需要继承的类是：

1
2
ceres::SizedCostFunction<1, 3>

定义我们的CostFunction类

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// 解析求导的Cost Function
class ExponentialResidualAnalytic : public ceres::SizedCostFunction<1, 3> {
public:
    // 构造函数
    ExponentialResidualAnalytic(double x, double y)
        : x_(x), y_(y) {}

    // 必须重写的Evaluate函数
    virtual bool Evaluate(double const* const* parameters,
                          double* residuals,
                          double** jacobians) const override;

private:
    // 存储数据点的成员变量
    const double x_;  // x坐标
    const double y_;  // y坐标
};

实现Evaluate函数

这是最重要的一步，它包含两部分：

  1. 计算残差值（必须）
  2. 计算雅可比矩阵（可选，但解析求导的关键）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
bool ExponentialResidualAnalytic::Evaluate(double const* const* parameters,
                                          double* residuals,
                                          double** jacobians) const {
    // 1. 获取参数值
    const double a = parameters[0][0];  // 第一个参数块的第一个元素
    const double b = parameters[0][1];  // 第一个参数块的第二个元素
    const double c = parameters[0][2];  // 第一个参数块的第三个元素

    // 2. 计算残差 r = y - exp(ax² + bx + c)
    const double exponent = a * x_ * x_ + b * x_ + c;
    const double prediction = std::exp(exponent);
    residuals[0] = y_ - prediction;

    // 3. 如果需要计算雅可比矩阵（这是一种优化策略，有时候迭代过程中不需要计算就会传入null）
    if (jacobians != nullptr && jacobians[0] != nullptr) {
        // 我们要计算残差对参数a, b, c的偏导数

        // 预先计算的exp项（提高效率）
        const double exp_term = std::exp(exponent);

        // dr/da = -exp(ax² + bx + c) * x²
        jacobians[0][0] = -exp_term * x_ * x_;

        // dr/db = -exp(ax² + bx + c) * x
        jacobians[0][1] = -exp_term * x_;

        // dr/dc = -exp(ax² + bx + c)
        jacobians[0][2] = -exp_term;
    }

    // 返回计算成功
    return true;
}

3.1.3 Evaluate函数参数详解

1. double const* const* parameters

  • 这是一个指向指针数组的指针
  • parameters[i]访问第i个参数块（我们只有一个参数块，所以使用parameters[0]）
  • parameters[0][j]访问第一个参数块的第j个元素

2. double* residuals

  • 用于存储计算的残差值
  • residuals[i]写入第i个残差值（我们只有一个残差，所以使用residuals[0]）

3. double** jacobians

  • 用于存储计算的雅可比矩阵
  • 如果为nullptr，表示优化器不需要计算雅可比矩阵
  • jacobians[i]指向与第i个参数块相关的雅可比矩阵
  • jacobians[0][j]存储残差对第一个参数块中第j个参数的偏导数

3.2 方法二：自动求导 (Automatic Differentiation)

这是 Ceres 推荐的方法之一，因为它结合了易用性和效率。我们只需要定义一个计算残差的模板化 “Functor”，Ceres 会自动计算导数。

3.2.1 Functor 结构基本形式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
struct ExponentialResidualAutoDiff {
    // 构造函数：保存每个残差项所需的输入数据
    ExponentialResidualAutoDiff(double x, double y) : x_(x), y_(y) {}

    // 模板化的 operator()：这是自动求导的关键
    template <typename T>
    bool operator()(const T* const parameters, T* residual) const {
        // 从参数数组中提取参数
        // 执行计算
        // 设置残差
        return true;
    }

private:
    // 存储每个数据点的信息
    const double x_;
    const double y_;
};

3.2.2 关键点详解

模板化 operator()

1
2
3
template <typename T>
bool operator()(const T* const parameters, T* residual) const

这是自动求导的核心。这里的类型 T 是关键，Ceres 会：

  • 使用 T = double 调用此函数来计算残差值
  • 使用 T = ceres::Jet<double, N> 调用此函数来计算导数

ceres::Jet 类型是 Ceres 的特殊类型，它存储了值和对应的导数信息。当函数使用 Jet 类型执行计算时，它自动跟踪导数。

参数和返回值

  • const T* const parameters：包含优化参数的数组
  • T* residual：输出残差的数组
  • 返回 bool：表示计算是否成功

3.2.3 详细示例：分步解析

对于模型 $y = \exp(ax^2 + bx + c)$，以下是详细的自动求导 Functor 实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
struct ExponentialResidualAutoDiff {
    // 构造函数：存储每个数据点的 x 和 y 值
    ExponentialResidualAutoDiff(double x, double y) : x_(x), y_(y) {}

    template <typename T>
    bool operator()(const T* const parameters, T* residual) const {
        // 步骤 1: 从参数数组中提取 a, b, c
        const T a = parameters[0];
        const T b = parameters[1];
        const T c = parameters[2];

        // 步骤 2: 将常量 x_ 和 y_ 转换为类型 T
        // 这是必要的，因为所有计算都需要在 T 类型的环境中进行
        const T x = T(x_);
        const T y = T(y_);

        // 步骤 3: 计算模型预测值
        // 注意：所有数学运算都使用 T 类型的变量
        T exponent = a * x * x + b * x + c;
        T prediction = exp(exponent);  // 使用与 T 兼容的 exp 函数

        // 步骤 4: 计算残差 (observed - predicted)
        residual[0] = y - prediction;

        // 计算成功
        return true;
    }

private:
    // 数据点的 x 和 y 值
    const double x_;
    const double y_;
};

3.3 方法三：数值求导 (Numeric Differentiation)

这种方法最容易实现，因为我们只需要提供计算残差的函数。Ceres 会通过微小地改变参数值并观察残差的变化来近似计算雅可比矩阵。但它通常比自动或解析求导慢，并且可能不太精确。

3.3.1 实现 Functor (与自动求导类似，但非模板化)

创建一个结构体（或类），重载 operator()，但这次只使用 double 类型。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// 数值求导的 Functor
struct ExponentialResidualNumeric {
    ExponentialResidualNumeric(double x, double y) : x_(x), y_(y) {}

    // 注意：这里 operator() 不是模板化的
    bool operator()(const double* const parameters, double* residual) const {
        const double a = parameters[0];
        const double b = parameters[1];
        const double c = parameters[2];

        const double x = x_;
        const double y = y_;

        // 计算预测值 prediction = exp(ax^2 + bx + c)
        double exponent_term = a * x * x + b * x + c;
        double prediction = std::exp(exponent_term);

        // 计算残差 residual = y - prediction
        residual[0] = y - prediction;

        return true;
    }

private:
    const double x_; // 观测数据 x
    const double y_; // 观测数据 y
};

4. 后续步骤

流程：

  1. 定义残差函数：使用解析求导、自动求导或数值求导

  2. 创建和配置 Problem：

    在 Ceres 中，Problem 是整个优化问题的容器，它包含了所有残差块和参数块。

    1
    2
    ceres::Problem problem;

    上面这行代码创建了一个空的 Problem 对象。接下来，我们需要向它添加残差块（即我们前面定义的成本函数）。

    • 添加残差块

    对于每个数据点 $(x_i, y_i)$，我们需要添加一个残差块。这是通过 AddResidualBlock 方法完成的：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    for (size_t i = 0; i < x_data.size(); ++i) {
        ceres::CostFunction* cost_function =
            new ceres::AutoDiffCostFunction<ExponentialResidualAutoDiff, 1, 3>(
                new ExponentialResidualAutoDiff(x_data[i], y_data[i])
            );
        problem.AddResidualBlock(cost_function,
                                 nullptr, // 使用默认的平方损失函数
                                 parameters);
    }

    这段代码的详细解释：

    1. 创建成本函数：
      • 对于自动求导，我们使用 AutoDiffCostFunction<> 类
      • 模板参数 ExponentialResidualAutoDiff 是我们的函数对象类型
      • 1 表示残差的维度（每个观测产生一个标量残差）
      • 3 表示参数块的大小（我们的参数是 a, b, c 三个值）
      • 构造函数参数是我们的函数对象实例，它保存了特定的 $x_i$ 和 $y_i$ 值
    2. 添加残差块：

      • AddResidualBlock 将成本函数添加到问题中

      • 第一个参数是成本函数指针

      • 第二个参数是损失函数指针（nullptr 表示使用默认的平方损失函数）

        损失函数用于减轻离群点（outliers）的影响。当为 nullptr 时，Ceres 使用默认的平方损失函数，即直接对残差平方求和。

        如果需要使用鲁棒损失函数，比如 Huber 损失，可以这样写：

        1
        2
        3
        ceres::LossFunction* loss_function = new ceres::HuberLoss(1.0);
        problem.AddResidualBlock(cost_function, loss_function, parameters);

      • 第三个参数是指向我们要优化的参数数组的指针

    • 设置参数约束（如有需要）

      除了添加残差块外，还可以对参数设置一些约束：

      1. 固定某些参数（使它们在优化过程中不变）：

        1
        2
        problem.SetParameterBlockConstant(parameters + 2); // 固定第三个参数 c

      2. 设置参数的上下界：

        1
        2
        3
        problem.SetParameterLowerBound(parameters, 0, 0.0); // 设置 a >= 0
        problem.SetParameterUpperBound(parameters, 1, 0.0); // 设置 b <= 0

  3. 配置求解器选项

    求解器负责实际运行优化算法。我们通过 Solver::Options 对象来配置求解器：

    1
    2
    3
    4
    5
    ceres::Solver::Options options;
    options.linear_solver_type = ceres::DENSE_QR;
    options.minimizer_progress_to_stdout = true;
    options.max_num_iterations = 50;

    详细解释：

    1. linear_solver_type：
      • 在每次迭代中，非线性最小二乘问题会被线性化，产生一个线性子问题
      • 这个选项决定如何求解这个线性子问题
      • DENSE_QR 适合小规模问题（参数和残差数量都不太多）
      • 对于大规模问题，可以考虑 SPARSE_NORMAL_CHOLESKY 或 ITERATIVE_SCHUR 等
    2. minimizer_progress_to_stdout：
      • 设为 true 时会在优化过程中打印进度信息
      • 非常有用，可以看到成本函数的下降情况和收敛过程
    3. max_num_iterations：
      • 最大迭代次数
      • 防止算法在难以收敛的情况下无限循环

    其他常用选项包括：

    • function_tolerance：当成本函数相对变化小于此值时停止
    • gradient_tolerance：当梯度范数小于此值时停止
    • parameter_tolerance：当参数相对变化小于此值时停止
    • num_threads：使用的线程数量（并行化）

  4. 运行求解器

    设置好 Problem 和 Options 后，调用 Solve 函数执行优化：

    1
    2
    ceres::Solver::Summary summary;
    ceres::Solve(options, &problem, &summary);

    这会执行优化算法并将结果存储在 summary 对象中。优化过程完成后，最优参数值会被直接写入到我们提供给 Problem 的 parameters 数组中。

  5. 分析结果

    优化完成后，通常需要：

    1. 检查优化状态：

      1
      std::cout << summary.BriefReport() << "\n";

      这会打印一个简短报告，包括终止状态、迭代次数、求解时间等

    2. 查看最终参数值：

      1
      2
      3
      4
      std::cout << "a = " << parameters[0] << " (True: " << a_true << ")\n";
      std::cout << "b = " << parameters[1] << " (True: " << b_true << ")\n";
      std::cout << "c = " << parameters[2] << " (True: " << c_true << ")\n";

    3. 计算最终代价：

      1
      2
      3
      4
      5
      6
      7
      8
      double final_cost = 0;
      for (size_t i = 0; i < x_data.size(); ++i) {
          double residual = y_data[i] - std::exp(parameters[0] * x_data[i] * x_data[i] +
                                              parameters[1] * x_data[i] + parameters[2]);
          final_cost += residual * residual;
      }
      std::cout << "Final Total Cost: " << final_cost << std::endl;

5. 完整代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
// 用于生成模拟数据的真实参数
const double a_true = 1.0;
const double b_true = 2.0;
const double c_true = 1.0;

// 生成带噪声的数据
void generate_data(std::vector<double>& x_data, std::vector<double>& y_data, int num_points = 100, double noise_stddev = 0.1) {
    x_data.resize(num_points);
    y_data.resize(num_points);
    std::default_random_engine generator(123); // Use a fixed seed for reproducibility
    std::normal_distribution<double> distribution(0.0, noise_stddev); // 高斯噪声

    std::cout << "Generated Data:" << std::endl;
    for (int i = 0; i < num_points; ++i) {
        double x = static_cast<double>(i) / (num_points / 10.0); // x 从 0 到 9.9 or similar range
        x_data[i] = x;
        double y_clean = std::exp(a_true * x * x + b_true * x + c_true);
        double noise = distribution(generator);
        y_data[i] = y_clean + noise;
        // Optional: Print generated data to verify
        // std::cout << "x: " << x << ", y_clean: " << y_clean << ", y_observed: " << y_data[i] << std::endl;
    }
     std::cout << "Data generation complete. " << num_points << " points generated." << std::endl;
}

int main(int argc, char** argv) {
    google::InitGoogleLogging(argv[0]); // 初始化 glog

    // 1. 生成数据
    std::vector<double> x_data;
    std::vector<double> y_data;
    generate_data(x_data, y_data, 100, 0.2); // 使用 100 个点，噪声标准差 0.2

    // 2. 设置初始参数估计值
    double initial_a = 0.9;
    double initial_b = 1.1;
    double initial_c = 0.9;
    double parameters[3] = {initial_a, initial_b, initial_c};

    // 3. 构建 Ceres 问题 (Problem)
    ceres::Problem problem;

    // --- 选择一种求导方法并将残差块添加到 Problem ---
    // --- 一次只能取消注释一个方法块来测试 ---

    // 方法一：使用解析求导
    /*
    std::cout << "\nUsing Analytic Differentiation:\n";
    for (size_t i = 0; i < x_data.size(); ++i) {
        ceres::CostFunction* cost_function =
            new ExponentialResidualAnalytic(x_data[i], y_data[i]);
        problem.AddResidualBlock(cost_function,
                                 nullptr, // 使用默认的平方损失函数
                                 parameters);
    }
    */

    // 方法二：使用自动求导
    std::cout << "\nUsing Automatic Differentiation:\n";
    for (size_t i = 0; i < x_data.size(); ++i) {
        ceres::CostFunction* cost_function =
            new ceres::AutoDiffCostFunction<ExponentialResidualAutoDiff, 1, 3>(
                new ExponentialResidualAutoDiff(x_data[i], y_data[i])
            );
        problem.AddResidualBlock(cost_function,
                                 nullptr, // 使用默认的平方损失函数
                                 parameters);
    }

    // 方法三：使用数值求导
    /*
    std::cout << "\nUsing Numeric Differentiation:\n";
    for (size_t i = 0; i < x_data.size(); ++i) {
        ceres::CostFunction* cost_function =
            new ceres::NumericDiffCostFunction<ExponentialResidualNumeric, ceres::CENTRAL, 1, 3>(
                new ExponentialResidualNumeric(x_data[i], y_data[i])
            );
        problem.AddResidualBlock(cost_function,
                                 nullptr, // 使用默认的平方损失函数
                                 parameters);
    }
    */

    // 4. 配置求解器 (Solver)
    ceres::Solver::Options options;
    options.linear_solver_type = ceres::DENSE_QR; // 对于小型问题，DENSE_QR 通常足够
    options.minimizer_progress_to_stdout = true;  // 将优化过程输出到控制台
    options.max_num_iterations = 50;             // 最大迭代次数

    // 5. 运行求解器
    ceres::Solver::Summary summary;
    ceres::Solve(options, &problem, &summary);

    // 6. 输出结果
    std::cout << "\nSolver Summary:\n";
    std::cout << summary.BriefReport() << "\n";
    // std::cout << summary.FullReport() << "\n"; // 更详细的报告

    std::cout << "\nEstimated parameters:\n";
    std::cout << "a = " << parameters[0] << " (True: " << a_true << ")\n";
    std::cout << "b = " << parameters[1] << " (True: " << b_true << ")\n";
    std::cout << "c = " << parameters[2] << " (True: " << c_true << ")\n";

    // 计算最终的总代价 (Sum of Squared Errors)
    double final_cost = 0;
     for (size_t i = 0; i < x_data.size(); ++i) {
         double residual = y_data[i] - std::exp(parameters[0] * x_data[i] * x_data[i] + parameters[1] * x_data[i] + parameters[2]);
         final_cost += residual * residual;
     }
     std::cout << "Final Total Cost (Sum of Squared Errors): " << final_cost << std::endl;

    return 0;
}

6. 如何选择求导方法

  • 自动求导 (AutoDiff): 通常是首选。它易于实现（只需写 Functor），准确（不像数值求导有近似误差），并且通常比数值求导快。Ceres 对其进行了高度优化。
  • 解析求导 (Analytic): 如果你能正确无误地推导出雅可比矩阵并实现它，这通常是最快的方法。但是，手动求导和实现很容易出错，尤其对于复杂的模型。调试起来也比较困难。
  • 数值求导 (Numeric): 最容易实现，因为你只需要提供计算残差的函数。但它通常是最慢的，并且其精度受限于步长选择（Ceres 会尝试选择合适的步长），可能在某些情况下导致收敛性问题。适合快速原型设计或当解析/自动求导难以实现时。