最小二乘问题详解9：使用Ceres求解非线性最小二乘
1 引言
在上一篇文章《最小二乘问题详解8：Levenberg-Marquardt方法》中，笔者使用 Eigen 实现了求解非线性最小二乘问题的 Levenberg-Marquardt 方法。不过在实际的工程实践中，更多的是使用像 Ceres Solver 这样成熟的、专门用于求解大规模非线性最小二乘问题的库。不过，因为有着极强的专业性，像 Ceres Solver 这样的库使用起来并不容易。如果是初次接触这方面知识的读者，非常建议先读一读本系列的前置文章。
2 实例
如果要构建安装 Ceres Solver，可以参考文章《CMake构建学习笔记30-Ceres Solver库的构建》。Ceres Solver 的构建过程还是挺麻烦的，推荐直接找预编译版本，比如 GISBasic3rdParty。
还是求解与《最小二乘问题详解8：Levenberg-Marquardt方法》一样的最小二乘问题，模型函数为：\(f(x; \boldsymbol{\theta}) = \exp(a x^2 + b x + c)\)，具体代码如下：
#include <ceres/autodiff_cost_function.h>
#include <ceres/ceres.h>
#include <ceres/covariance.h>
#include <iostream>
#include <random>
#include <vector>
using namespace std;
// 残差计算结构体（用于自动微分）
struct ExpModelResidual {
ExpModelResidual(double x, double y) : x_(x), y_(y) {}
// 模板 operator() 支持自动微分
template <typename T>
bool operator()(const T* const a, const T* const b, const T* const c,
T* residual) const {
// 计算指数部分: a*x^2 + b*x + c
T exponent = (*a) * T(x_) * T(x_) + (*b) * T(x_) + (*c);
// 防止 exp 溢出（Ceres 内部对梯度也有保护，但这里加一层更稳）
const T kMaxExp = T(300.0);
if (exponent > kMaxExp) exponent = kMaxExp;
if (exponent < -kMaxExp) exponent = -kMaxExp;
T y_pred = ceres::exp(exponent);
residual[0] = T(y_) - y_pred;
return true;
}
private:
const double x_;
const double y_;
};
int main() {
// ========================
// 1. 真实参数
// ========================
double a_true = 0.05, b_true = -0.4, c_true = 1.0;
cout << "真实参数: a=" << a_true << ", b=" << b_true << ", c=" << c_true
<< endl;
// ========================
// 2. 生成带噪声的数据
// ========================
int N = 50;
vector<double> x_data(N), y_data(N);
random_device rd;
mt19937 gen(rd());
uniform_real_distribution<double> x_dis(-5.0, 5.0);
normal_distribution<double> noise_gen(0.0, 0.1); // 噪声标准差 0.1
for (int i = 0; i < N; ++i) {
x_data[i] = x_dis(gen);
double y_true =
exp(a_true * x_data[i] * x_data[i] + b_true * x_data[i] + c_true);
y_data[i] = y_true + noise_gen(gen);
}
// ========================
// 3. 初始化参数（猜测）
// ========================
double a = 0.0, b = 0.0, c = 0.0;
cout << "初始猜测: a=" << a << ", b=" << b << ", c=" << c << endl;
// ========================
// 4. 构建 Ceres 问题
// ========================
ceres::Problem problem;
for (int i = 0; i < N; ++i) {
ceres::CostFunction* cost_function =
new ceres::AutoDiffCostFunction<ExpModelResidual, 1, 1, 1, 1>(
new ExpModelResidual(x_data[i], y_data[i]));
problem.AddResidualBlock(cost_function, nullptr, &a, &b, &c);
}
// ========================
// 5. 配置并运行求解器
// ========================
ceres::Solver::Options options;
options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;
options.minimizer_progress_to_stdout = true; // 打印迭代过程
options.max_num_iterations = 50;
options.function_tolerance = 1e-12;
options.gradient_tolerance = 1e-10;
options.parameter_tolerance = 1e-8;
ceres::Solver::Summary summary;
ceres::Solve(options, &problem, &summary);
cout << summary.BriefReport() << "\n";
// ========================
// 6. 输出结果
// ========================
cout << "--- 拟合完成 ---" << endl;
cout << "估计参数: a=" << a << ", b=" << b << ", c=" << c << endl;
cout << "真实参数: a=" << a_true << ", b=" << b_true << ", c=" << c_true
<< endl;
// 计算最终残差平方和
double final_cost = 0.0;
for (int i = 0; i < N; ++i) {
double pred = exp(a * x_data[i] * x_data[i] + b * x_data[i] + c);
double r = y_data[i] - pred;
final_cost += r * r;
}
cout << "最终残差平方和: " << final_cost << endl;
// ========================
// 7. （可选）计算参数协方差与标准差
// ========================
ceres::Covariance::Options cov_options;
ceres::Covariance covariance(cov_options);
vector<pair<const double*, const double*>> covariance_blocks;
covariance_blocks.emplace_back(&a, &a);
covariance_blocks.emplace_back(&b, &b);
covariance_blocks.emplace_back(&c, &c);
if (!covariance.Compute(covariance_blocks, &problem)) {
cerr << "协方差计算失败！" << endl;
return 1;
}
double cov_a, cov_b, cov_c;
covariance.GetCovarianceBlock(&a, &a, &cov_a);
covariance.GetCovarianceBlock(&b, &b, &cov_b);
covariance.GetCovarianceBlock(&c, &c, &cov_c);
// 估计噪声方差 σ² = RSS / (N - p)
double sigma2 = final_cost / (N - 3);
double std_a = sqrt(cov_a * sigma2);
double std_b = sqrt(cov_b * sigma2);
double std_c = sqrt(cov_c * sigma2);
cout << "\n参数标准差 (近似):" << endl;
cout << "a: ±" << std_a << endl;
cout << "b: ±" << std_b << endl;
cout << "c: ±" << std_c << endl;
return 0;
}
可以看到，Ceres Solver 的实现思路与原生的基于 Eigen 的实现思路完全不同。原生的基于 Eigen 的实现只要按照 Levenberg-Marquardt 方法的原理来实现即可，反而更容易理解；但是 Ceres Solver 的实现则需要考虑通用性，提供的接口更加抽象。所以这也是笔者写那么多前置文章的原因，对于一个需要专业知识的程序库，如果你不理解其中的原理，很可能也没法正确地使用它。
2.1 构建问题
从前置文章中可以知道，无论是 Gauss-Newton 法、梯度下降法还是 Levenberg-Marquardt 方法，关键的问题在于求解问题模型的雅可比矩阵。但是，对于用户来说更加熟悉的是最小二乘问题的残差，也就是\(y - f(x; \boldsymbol{\theta})\)。因此，Ceres Solver 的核心思想就是以残差为中心组织代码，关于雅可比矩阵的数值计算则作为 Ceres Solver 优化器的内部机制。参看如下关键代码：
// ========================
// 4. 构建 Ceres 问题
// ========================
ceres::Problem problem;
for (int i = 0; i < N; ++i) {
ceres::CostFunction* cost_function =
new ceres::AutoDiffCostFunction<ExpModelResidual, 1, 1, 1, 1>(
new ExpModelResidual(x_data[i], y_data[i]));
problem.AddResidualBlock(cost_function, nullptr, &a, &b, &c);
}
这里的 CostFunction（代价函数）是什么意思呢？简而言之，可以将其理解成一个残差项\(r_i(\theta) = y_i - f(x_i; \theta)\)；但是更准确地说，它是一个“残差块生成器”——封装了“如何计算一个残差向量及其对参数的导数”的完整逻辑。在 Ceres Solver 中：ceres::CostFunction 是一个接口类，它的核心方法是：
bool Evaluate(double const* const* parameters,
double* residuals,
double** jacobians);
从这个接口参数可以看到，ceres::CostFunction不仅可以计算 residuals（即 \(r_i\)），还可以计算 jacobians（即 \(\frac{\partial r_i}{\partial \theta}\)）。所以，它不仅仅是残差函数，而是残差 + 导数的联合计算单元。
Ceres Solver 支持三种雅可比矩阵的计算方式：
- 自动微分（
AutoDiffCostFunction）：精度高、无需手写Jacobian。例如这里的ExpModelResidual。 - 数值微分（
NumericDiffCostFunction）：用有限差分近似导数，用于黑盒函数。 - 解析微分（继承
CostFunction）：用户自定义Jacobian，适用于需要极致性能的情况。
但无论哪种方式，残差的定义不变，一个 ceres::CostFunction代表一个残差项，最后都通过 AddResidualBlock接口将残差项对象添加到定义好的问题模型 ceres::Problem中。大概来说，可以这样理解：
| 数学概念 | Ceres 实现 |
|---|---|
| 残差函数\(r_i(\theta)\) | ExpModelResidual::operator()（只算值） |
| 残差 + Jacobian 计算 | ceres::AutoDiffCostFunction<...>（自动微分包装后） |
| 可被优化器调用的完整残差项 | ceres::CostFunction* |
展开自定义的残差计算单元 ExpModelResidual：
// 残差计算结构体（用于自动微分）
struct ExpModelResidual {
ExpModelResidual(double x, double y) : x_(x), y_(y) {}
// 模板 operator() 支持自动微分
template <typename T>
bool operator()(const T* const a, const T* const b, const T* const c,
T* residual) const {
// 计算指数部分: a*x^2 + b*x + c
T exponent = (*a) * T(x_) * T(x_) + (*b) * T(x_) + (*c);
// 防止 exp 溢出（Ceres 内部对梯度也有保护，但这里加一层更稳）
const T kMaxExp = T(300.0);
if (exponent > kMaxExp) exponent = kMaxExp;
if (exponent < -kMaxExp) exponent = -kMaxExp;
T y_pred = ceres::exp(exponent);
residual[0] = T(y_) - y_pred;
return true;
}
private:
const double x_;
const double y_;
};
可以看到，ExpModelResidual::operator() 的实现本质上就是残差函数 \(r_i(\theta)\) 的计算过程。当它被 ceres::AutoDiffCostFunction 包装后，Ceres 能够自动完成求导，从而生成对应的雅可比矩阵项。其关键在于这一行代码：
T y_pred = ceres::exp(exponent);
这里使用的是 ceres::exp，而非标准库中的 std::exp。这是因为在自动微分（AutoDiffCostFunction）过程中，Ceres 会将模板参数 T 实例化为一种特殊的数值类型——ceres::Jet<T, N>。该类型不仅存储函数值，还同时携带其对若干参数的偏导数。为了支持这种类型的运算，Ceres 为常见的数学函数提供了重载版本，这些版本能够正确传播导数信息（即应用链式法则）。
虽然底层涉及 C++ 模板和编译期“计算图录制”的机制，理解起来有一定门槛，但用户只需遵循一个简单规则：在残差函数的实现中，所有数学运算必须使用 ceres:: 命名空间下的函数，而不能使用 std:: 或全局命名空间中的版本。具体的需要替换的常见函数对照表如下：
| 标准写法（❌ 错误） | Ceres 正确写法（✅ 必须使用） | 说明 |
|---|---|---|
std::exp(x) |
ceres::exp(x) |
指数函数 |
std::log(x) |
ceres::log(x) |
自然对数 |
std::log10(x) |
ceres::log10(x) |
常用对数 |
std::sqrt(x) |
ceres::sqrt(x) |
平方根 |
std::pow(x, y) |
ceres::pow(x, y) |
幂函数（注意：y 也应为模板类型 T） |
std::sin(x) |
ceres::sin(x) |
正弦函数 |
std::cos(x) |
ceres::cos(x) |
余弦函数 |
std::tan(x) |
ceres::tan(x) |
正切函数 |
std::asin(x) |
ceres::asin(x) |
反正弦函数 |
std::acos(x) |
ceres::acos(x) |
反余弦函数 |
std::atan(x) |
ceres::atan(x) |
反正切函数 |
std::sinh(x) |
ceres::sinh(x) |
双曲正弦 |
std::cosh(x) |
ceres::cosh(x) |
双曲余弦 |
std::tanh(x) |
ceres::tanh(x) |
双曲正切 |
std::abs(x) / std::fabs(x) |
ceres::abs(x) |
绝对值（⚠️ 在 0 处不可导，慎用） |
std::max(a, b) |
ceres::max(a, b) |
最大值（⚠️ 在相等点不可导，需谨慎） |
std::min(a, b) |
ceres::min(a, b) |
最小值（⚠️ 在相等点不可导，需谨慎） |
只要遵守上述规范，Ceres 就能在运行时自动、高效且精确地计算出残差及其雅可比矩阵，无需用户手动推导或实现导数逻辑。
2.2 配置参数
接下来看一下核心配置部分的代码：
// ========================
// 5. 配置并运行求解器
// ========================
ceres::Solver::Options options;
options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;
options.minimizer_progress_to_stdout = true; // 打印迭代过程
options.max_num_iterations = 50;
options.function_tolerance = 1e-12;
options.gradient_tolerance = 1e-10;
options.parameter_tolerance = 1e-8;
逐项详解求解器的配置参数：
options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;：指定每次迭代中求解线性子问题（即 \((J^T J + \lambda D)\Delta\theta = J^T r\)）所用的线性代数求解器类型。常见的选项如下所示：
| 类型 | 适用场景 | 说明 |
|---|---|---|
DENSE_NORMAL_CHOLESKY |
小规模稠密问题（\(p < 100\)） | 构造正规方程\(A = J^T J\)，用 Cholesky 分解求解。 |
DENSE_QR |
小规模，但\(J\) 可能病态 | 直接对\(J\) 做 QR 分解，更稳定但稍慢 |
SPARSE_NORMAL_CHOLESKY |
大规模稀疏问题 （如 SLAM、Bundle Adjustment） |
利用\(J^T J\) 的稀疏性 |
ITERATIVE_SCHUR / CGNR |
超大规模问题 | 迭代法，内存友好 |
options.minimizer_progress_to_stdout = true;：是否将优化过程的迭代信息打印到标准输出。通过打印内容可以看迭代过程是否收敛、是否震荡、步长是否合理。options.max_num_iterations = 50;：最大迭代次数。即使未收敛，达到此次数后强制停止。防止算法陷入死循环，控制计算时间上限。对于简单问题，20~50 次就足够了； 复杂 BA 这样的复杂问题，可能需要 100~200 次。options.function_tolerance = 1e-12;：目标函数（cost）的相对变化容差，即\(\frac{S(\theta_k) - S(\theta_{k+1})}{S(\theta_k)} < \texttt{function\_tolerance}\)。当连续两次迭代的 cost 相对变化小于此值时，认为收敛。高精度拟合可设置为1e-12~1e-15；一般应用可设置为1e-6~1e-9。options.gradient_tolerance = 1e-10;：梯度范数的绝对容差。当 \(\| \nabla S(\theta) \| = \| J^T r \|\) 小于此值时，认为达到极小点。通常设为1e-8~1e-12，不受 cost 绝对大小影响，比function_tolerance更可靠，是最常用的收敛判据之一。options.parameter_tolerance = 1e-8;：参数更新量的绝对容差。当 \(\| \Delta \theta \|\) 小于此值时，认为参数不再显著变化。注意如果参数尺度差异大（如 a≈0.01, b≈1000），此判据可能不合理，根据参数实际量级调整。
2.3 优化报告
案例运行的结果如下所示：
真实参数: a=0.05, b=-0.4, c=1
初始猜测: a=0, b=0, c=0
iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time
0 5.127046e+03 0.00e+00 4.99e+03 0.00e+00 0.00e+00 1.00e+04 0 7.05e-05 2.06e-04
1 8.359362e+34 -8.36e+34 4.99e+03 0.00e+00 -1.83e+31 5.00e+03 1 1.35e-04 6.28e-04
2 8.256445e+34 -8.26e+34 4.99e+03 0.00e+00 -1.81e+31 1.25e+03 1 1.89e-05 1.15e-03
3 7.666005e+34 -7.67e+34 4.99e+03 0.00e+00 -1.68e+31 1.56e+02 1 2.89e-05 1.56e-03
4 3.875752e+34 -3.88e+34 4.99e+03 0.00e+00 -8.50e+30 9.77e+00 1 1.43e-05 1.90e-03
5 2.984885e+30 -2.98e+30 4.99e+03 0.00e+00 -6.64e+26 3.05e-01 1 3.58e-05 2.18e-03
6 4.698247e+07 -4.70e+07 4.99e+03 0.00e+00 -2.43e+04 4.77e-03 1 6.20e-06 2.45e-03
7 5.075148e+03 5.19e+01 5.81e+03 0.00e+00 1.08e+00 1.43e-02 1 4.67e-05 2.84e-03
8 4.873439e+03 2.02e+02 9.07e+03 1.08e-01 1.26e+00 4.29e-02 1 2.08e-05 3.01e-03
9 3.776580e+03 1.10e+03 2.66e+04 2.86e-01 1.84e+00 1.29e-01 1 1.92e-05 3.29e-03
10 9.518702e+02 2.82e+03 4.90e+04 3.75e-01 1.74e+00 3.86e-01 1 1.77e-05 3.48e-03
11 1.461114e+02 8.06e+02 2.43e+04 1.29e-01 1.12e+00 1.16e+00 1 1.71e-05 3.69e-03
12 3.388251e+01 1.12e+02 3.01e+03 5.23e-02 1.02e+00 3.48e+00 1 1.71e-05 3.88e-03
13 2.579351e+01 8.09e+00 1.46e+03 2.50e-02 1.01e+00 1.04e+01 1 1.55e-05 4.20e-03
14 1.662811e+01 9.17e+00 1.35e+03 3.54e-02 9.96e-01 3.13e+01 1 1.38e-05 4.56e-03
15 6.626112e+00 1.00e+01 7.18e+02 4.17e-02 9.80e-01 9.39e+01 1 1.52e-05 4.70e-03
16 1.633045e+00 4.99e+00 2.45e+02 2.66e-02 9.64e-01 2.82e+02 1 4.87e-05 5.15e-03
17 3.715462e-01 1.26e+00 5.74e+01 2.98e-02 9.75e-01 8.45e+02 1 4.89e-05 5.78e-03
18 2.207972e-01 1.51e-01 8.88e+00 1.71e-02 9.98e-01 2.53e+03 1 5.97e-05 6.43e-03
19 2.156097e-01 5.19e-03 6.47e-01 3.87e-03 1.01e+00 7.60e+03 1 3.51e-05 7.07e-03
20 2.155782e-01 3.15e-05 1.99e-02 3.20e-04 1.01e+00 2.28e+04 1 3.61e-05 7.69e-03
21 2.155781e-01 3.38e-08 3.28e-04 1.05e-05 1.01e+00 6.84e+04 1 6.85e-05 8.13e-03
22 2.155781e-01 1.07e-11 3.90e-06 1.87e-07 1.01e+00 2.05e+05 1 2.97e-05 8.33e-03
Ceres Solver Report: Iterations: 23, Initial cost: 5.127046e+03, Final cost: 2.155781e-01, Termination: CONVERGENCE
--- 拟合完成 ---
估计参数: a=0.0500655, b=-0.400392, c=0.996087
真实参数: a=0.05, b=-0.4, c=1
最终残差平方和: 0.431156
参数标准差 (近似):
a: ±0.000357974
b: ±0.00223108
c: ±0.00464517
大部分内容是 Ceres 输出的优化报告：
ceres::Solver::Summary summary;
ceres::Solve(options, &problem, &summary);
cout << summary.BriefReport() << "\n";
先看看最终汇总行的内容：
Ceres Solver Report: Iterations: 23, Initial cost: 5.127046e+03, Final cost: 2.155781e-01, Termination: CONVERGENCE
汇总行中字段的含义是：
| 字段 | 含义 | 结果解读 |
|---|---|---|
| Iterations: 23 | 总迭代次数（从第 0 次到第 22 次） | 共执行 23 轮优化 |
| Initial cost | 初始目标函数值（即\(\frac{1}{2} \sum r_i^2\)） | 初值(0,0,0) 下误差很大（≈5127） |
| Final cost | 最终目标函数值 | 收敛到 ≈0.2156，非常小 |
| Termination: CONVERGENCE | 终止原因 | ✅ 正常收敛（满足梯度或步长停止条件） |
注意 Ceres 中\(cost = \frac{1}{2} \sum r_i^2\)，这是为了在计算 Jacobian 和梯度的的时候无需额外除以 2 或乘以 2，使用 LM 方法求解就能使用更为简洁的形式。
而在具体的迭代日志中，字段的含义如下：
| 列名 | 全称 / 数学含义 | 单位 / 类型 | 说明 |
|---|---|---|---|
iter |
迭代序号 | 整数 | 从 0 开始计数 |
cost |
当前总代价\(C(\theta) = \frac{1}{2} \sum_{i=1}^N r_i(\theta)^2\) | 标量（浮点） | 越小越好；理想情况趋近于 0 |
cost_change |
上次 cost − 本次 cost | 浮点 | • 正值：代价下降（好）• 负值：代价上升（异常，通常因步长过大）• 接近 0：趋于收敛 |
gradient |
梯度范数\(|\nabla_\theta C(\theta)|_2\) | 浮点 | 衡量当前点是否接近极小值：• 大 → 远离最优• 小（如 < 1e-6）→ 接近收敛 |
step |
参数更新步长\(|\Delta \theta|_2\) | 浮点 | • 大：参数剧烈变化• 小（如 < 1e-6）→ 参数几乎不变，可能收敛 |
tr_ratio |
信任域比率（Trust Region Ratio）$$\rho = \dfrac{\text{实际代价下降}}{\text{局部模型预测下降}}$$ | 无量纲 | •\(\rho \approx 1\)：模型准确，可增大步长• \(\rho \ll 0\)：预测失败，需缩小步长• \(\rho < 0\)：代价反而上升（步长过大） |
tr_radius |
信任域半径（Trust Region Radius） | 参数空间尺度 | 控制最大允许步长：• 小 → 保守更新• 大 → 激进更新（接近牛顿法） |
ls_iter |
线搜索迭代次数 | 整数 | 在 Levenberg-Marquardt（LM）模式下恒为 1，因为 LM 是信任域方法，不使用线搜索 |
iter_time |
本次迭代耗时 | 秒（s） | 通常为微秒级（1e-5 ~ 1e-4 s） |
total_time |
累计总耗时 | 秒（s） | 从开始到当前迭代的总时间 |
从迭代日志中可以看到，迭代过程大概分成三个阶段：
- 初期爆炸（iter 1–6）：
cost从 5e3 飙升到 8e34。原因是初始猜测(a,b,c)=(0,0,0)导致模型 \(y = e^0 = 1\)，但真实数据可能远大于 1。第一步尝试大步长，但exp(a x² + ...)对参数极度敏感，导致数值溢出。Ceres 自动缩小tr_radius（1e4 → 1e-3），稳住优化——这是 LM 算法鲁棒性的体现。 - 中期快速下降（iter 7–15）：
cost从 5e3 快速降至 6.6；tr_ratio ≈ 1，说明局部线性模型有效；tr_radius开始扩大（1e-2 → 1e2），进入高效牛顿阶段。 - 后期精细收敛（iter 16–22）：
gradient从 245 降至3.9e-6；step降至1.87e-7；tr_ratio ≈ 1.01；直到满足默认收敛条件（比如梯度足够小），正常终止。
2.4 输出精度
精度是用户最为关心的问题，ceres::Covariance是 Ceres 提供的一个后处理工具类，用于在优化完成后估计参数的协方差矩阵，从而得到每个参数的不确定性（标准差）。
首先是使用默认选项创建协方差计算器，其中 cov_options可设置算法（如是否使用 sparse）、内存策略等：
// 1. 创建 Covariance 对象
ceres::Covariance::Options cov_options;
ceres::Covariance covariance(cov_options);
为了提升效率，Covariance::Compute 不会自动计算所有参数的完整协方差矩阵，因此需要显式声明只关心 a、b、c 各自的方差，即对角线元素：
// 2. 指定要计算哪些参数块之间的协方差
vector<pair<const double*, const double*>> covariance_blocks;
covariance_blocks.emplace_back(&a, &a); // a 与 a 的协方差 → 方差
covariance_blocks.emplace_back(&b, &b); // b 的方差
covariance_blocks.emplace_back(&c, &c); // c 的方差
Ceres 会在当前参数值（即优化后的 a, b, c）处重新计算 Jacobian，然后构建并求解 \((J^\top J)^{-1}\) 的指定子块：
// 3. 执行协方差计算（在 problem 的当前解处线性化）
if (!covariance.Compute(covariance_blocks, &problem)) {
cerr << "协方差计算失败！" << endl;
return 1;
}
提取各参数的归一化方差：
// 4. 提取各参数的“归一化”方差（即 (JᵀJ)⁻¹ 的对角元）
double cov_a, cov_b, cov_c;
covariance.GetCovarianceBlock(&a, &a, &cov_a); // cov_a = [(JᵀJ)⁻¹]_{aa}
covariance.GetCovarianceBlock(&b, &b, &cov_b);
covariance.GetCovarianceBlock(&c, &c, &cov_c);
计算噪声尺度 \(\sigma^2\)：
// 5. 估计噪声方差 σ² = RSS / (N - p)
double sigma2 = final_cost / (N - 3);
计算最终标准差，也就是最终的精度：
// 6. 计算最终标准差：std = sqrt(σ² × (JᵀJ)⁻¹_ii)
double std_a = sqrt(cov_a * sigma2);
double std_b = sqrt(cov_b * sigma2);
double std_c = sqrt(cov_c * sigma2);
3 补充
3.1 资源管理
回到构建 Ceres 问题的关键代码：
// ========================
// 4. 构建 Ceres 问题
// ========================
ceres::Problem problem;
for (int i = 0; i < N; ++i) {
ceres::CostFunction* cost_function =
new ceres::AutoDiffCostFunction<ExpModelResidual, 1, 1, 1, 1>(
new ExpModelResidual(x_data[i], y_data[i]));
problem.AddResidualBlock(cost_function, nullptr, &a, &b, &c);
}
可以看到这里使用了 new但是没有对应的 delete，不是代码中存在遗漏，而是 Ceres 的 API 设计是基于裸指针 + 显式所有权转移来实现的。这种设计在 C++ 中非常常见（比如Qt），简而言之，就是一个对象会托管住其数据成员的所有权。在这里就是 Problem会在自身析构时自动 delete所有它拥有的 CostFunction,AutoDiffCostFunction又会托管住一个 new出来的 ExpModelResidual*。像这样一层一层嵌套，只用管理 Problem对象就可以控制整个链条的资源。
3.2 自动求导
在 ExpModelResidual中，虽然定义了残差计算:
\[r_i(a,b,c) = y_i - exp(a*x_i² + b*x_i + c)
\]
但 Ceres 需要：
- 残差值：用于计算总代价 \(\sum r_i^2\)
- 雅可比矩阵：用于构建 \(J^T J\) 和梯度 \(J^T r\)
Ceres 的解决方案是通过模板参数类型 T 的不同，让同一个 operator() 既能算值又能算导数:
template <typename T>
bool operator()(const T* a, const T* b, const T* c, T* residual) const {
T exponent = (*a) * T(x_) * T(x_) + (*b) * T(x_) + (*c);
T y_pred = ceres::exp(exponent); // ← 重要！用 ceres::exp
residual[0] = T(y_) - y_pred;
return true;
}
在这里 T可以是 double（纯数值），也可以是 ceres::Jet<double, 3>（数值+导数）。自动微分包装器 AutoDiffCostFunction使用了这个 ExpModelResidual作为模板参数：
new ceres::AutoDiffCostFunction<ExpModelResidual, 1, 1, 1, 1>(...)
这些模板参数含义是：
| 参数 | 含义 |
|---|---|
ExpModelResidual |
你的残差计算类 |
1 |
残差维度（每个观测产生 1 个残差） |
1, 1, 1 |
三个参数块的维度（a、b、c 各 1 维） |
当 Ceres 执行 LM 算法时，会分两步调用 AutoDiffCostFunction。当需要计算残差值时，Ceres 内部调用：
// T = double
ExpModelResidual r(x_i, y_i);
double residual_val;
r(&a_val, &b_val, &c_val, &residual_val); // 正常数值计算
当需要计算雅可比矩阵时。Ceres 构造特殊的 Jet 数值：
// T = ceres::Jet<double, 3>
ceres::Jet<double, 3> a_jet(0.05, {1, 0, 0}); // 当前 a=0.05，∂a/∂a=1
ceres::Jet<double, 3> b_jet(-0.4, {0, 1, 0}); // ∂b/∂b=1
ceres::Jet<double, 3> c_jet(1.0, {0, 0, 1}); // ∂c/∂c=1
然后调用：
ExpModelResidual r(x_i, y_i);
ceres::Jet<double, 3> residual_jet;
r(&a_jet, &b_jet, &c_jet, &residual_jet); // 关键！
注意这里 ceres::Jet<T, N> 是一个模板结构体，大概定义如下：
template<typename T, int N>
struct Jet {
T a; // 函数值（value）
Eigen::Matrix<T, N, 1> v; // 对 N 个参数的偏导数（gradient）
Jet() : a(0.0), v(Eigen::Matrix<T, N, 1>::Zero()) {}
Jet(const T& value) : a(value), v(Eigen::Matrix<T, N, 1>::Zero()) {}
Jet(const T& value, const Eigen::Matrix<T, N, 1>& derivatives)
: a(value), v(derivatives) {}
};
如果使用 Ceres 的内置函数计算 Jet 类型的数值，就可以自动传播导数。例如这里使用的 ceres::exp：
template<typename T, int N>
inline Jet<T, N> exp(const Jet<T, N>& x) {
T exp_a = std::exp(x.a); // 标量指数值
return Jet<T, N>(exp_a, exp_a * x.v); // 导数 = exp(x) * dx/dθ
}
代入到 ExpModelResidual的 operator()，其运算的过程就是：
// exponent = a*x² + b*x + c
// Jet 运算规则：值相加，导数也相加
exponent = a_jet * 4.0 + b_jet * 2.0 + c_jet; // x_i=2.0
// → exponent.a = 0.05*4 + (-0.4)*2 + 1.0 = 0.9
// → exponent.v = [4.0, 2.0, 1.0] ← 这是 ∂exponent/∂[a,b,c]
// y_pred = ceres::exp(exponent)
y_pred = ceres::exp(exponent);
// → y_pred.a = exp(0.9)
// → y_pred.v = exp(0.9) * [4.0, 2.0, 1.0] ← 链式法则！
// residual = y_true - y_pred
residual_jet = Jet(y_i, [0,0,0]) - y_pred;
// → residual.a = y_i - exp(0.9)
// → residual.v = [0,0,0] - exp(0.9)*[4,2,1] = -exp(0.9)*[4,2,1]
最终 residual_jet.v 就是雅可比行向量：
\[\left[ \frac{\partial r_i}{\partial a}, \frac{\partial r_i}{\partial b}, \frac{\partial r_i}{\partial c} \right] = -e^{ax^2+bx+c} \cdot [x^2, x, 1]
\]
4. 实践
尽管笔者在上一篇文章《最小二乘问题详解8：Levenberg-Marquardt方法》中手写实现了 Levenberg-Marquardt（LM）算法，但是求解非线性最小二乘问题是一个很复杂的工程，实践中更倾向于使用 Ceres 这样稳健、高效、通用的成熟库。Ceres 做的工作包含且不局限于：
- 线性代数求解器的深度优化：能够自动选择最优线性求解器，对于小规模稠密问题，可以配置
DENSE_NORMAL_CHOLESKY，使用普通 Cholesky 求解，对于大规模稀疏问题（如 SLAM、Bundle Adjustment），可以配置SPARSE_NORMAL_CHOLESKY，实现稀疏 Cholesky 求解。 - 支持多种高度优化的稀疏线性代数库：例如SuiteSparse、Eigen等，能处理数十万参数、上百万残差项的 Bundle Adjustment 问题。
- 雅可比矩阵的高效构建与存储：对于大规模非线性最小二乘问题，雅可比矩阵非常占用内存，例如 1M 观测数据 × 1K 待定参数，存储 double 类型的雅可比矩阵需要 8GB 的内存空间。而 Ceres 可以实现对雅可比矩阵按需计算 + 稀疏表示，用户只需提供每个残差块对局部参数的导数（即
AutoDiffCostFunction返回的小 Jacobian 块），自动拼接成全局雅可比矩阵，只存储非零块。 - 信任域策略的工业级鲁棒性：成熟的信任域管理，智能的 \(\lambda\) 调整逻辑（基于实际/预测下降比 \(\rho\)）；在噪声大、初值差、目标函数非凸的情况下仍能稳定收敛。
- 参数约束与边界处理：支持参数边界约束。
- 损失函数（Loss Function）抗外点能力：内置多种鲁棒损失函数，例如Huber、Cauchy、SoftLOne、Tukey 等，自动降低大残差的权重，抑制 outlier 影响。
- 并行化与性能工程：通过CPU多线程或者GPU并行计算残差和雅克比矩阵，提升问题优化效率。
5. 总结
不得不说，要理解一个具有专业背景的库不是那么容易的事情，即使你理解底层算法的原理，但是工程实践又是另外一回事。只有将算法原理与工程实践融汇贯通，才能真正解决问题并创造价值。
浙公网安备 33010602011771号