安阳市网站建设_网站建设公司_交互流畅度_seo优化

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 + cT 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.1for (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 + cT 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_time0  5.127046e+03    0.00e+00    4.99e+03   0.00e+00   0.00e+00  1.00e+04        0    7.05e-05    2.06e-041  8.359362e+34   -8.36e+34    4.99e+03   0.00e+00  -1.83e+31  5.00e+03        1    1.35e-04    6.28e-042  8.256445e+34   -8.26e+34    4.99e+03   0.00e+00  -1.81e+31  1.25e+03        1    1.89e-05    1.15e-033  7.666005e+34   -7.67e+34    4.99e+03   0.00e+00  -1.68e+31  1.56e+02        1    2.89e-05    1.56e-034  3.875752e+34   -3.88e+34    4.99e+03   0.00e+00  -8.50e+30  9.77e+00        1    1.43e-05    1.90e-035  2.984885e+30   -2.98e+30    4.99e+03   0.00e+00  -6.64e+26  3.05e-01        1    3.58e-05    2.18e-036  4.698247e+07   -4.70e+07    4.99e+03   0.00e+00  -2.43e+04  4.77e-03        1    6.20e-06    2.45e-037  5.075148e+03    5.19e+01    5.81e+03   0.00e+00   1.08e+00  1.43e-02        1    4.67e-05    2.84e-038  4.873439e+03    2.02e+02    9.07e+03   1.08e-01   1.26e+00  4.29e-02        1    2.08e-05    3.01e-039  3.776580e+03    1.10e+03    2.66e+04   2.86e-01   1.84e+00  1.29e-01        1    1.92e-05    3.29e-0310  9.518702e+02    2.82e+03    4.90e+04   3.75e-01   1.74e+00  3.86e-01        1    1.77e-05    3.48e-0311  1.461114e+02    8.06e+02    2.43e+04   1.29e-01   1.12e+00  1.16e+00        1    1.71e-05    3.69e-0312  3.388251e+01    1.12e+02    3.01e+03   5.23e-02   1.02e+00  3.48e+00        1    1.71e-05    3.88e-0313  2.579351e+01    8.09e+00    1.46e+03   2.50e-02   1.01e+00  1.04e+01        1    1.55e-05    4.20e-0314  1.662811e+01    9.17e+00    1.35e+03   3.54e-02   9.96e-01  3.13e+01        1    1.38e-05    4.56e-0315  6.626112e+00    1.00e+01    7.18e+02   4.17e-02   9.80e-01  9.39e+01        1    1.52e-05    4.70e-0316  1.633045e+00    4.99e+00    2.45e+02   2.66e-02   9.64e-01  2.82e+02        1    4.87e-05    5.15e-0317  3.715462e-01    1.26e+00    5.74e+01   2.98e-02   9.75e-01  8.45e+02        1    4.89e-05    5.78e-0318  2.207972e-01    1.51e-01    8.88e+00   1.71e-02   9.98e-01  2.53e+03        1    5.97e-05    6.43e-0319  2.156097e-01    5.19e-03    6.47e-01   3.87e-03   1.01e+00  7.60e+03        1    3.51e-05    7.07e-0320  2.155782e-01    3.15e-05    1.99e-02   3.20e-04   1.01e+00  2.28e+04        1    3.61e-05    7.69e-0321  2.155781e-01    3.38e-08    3.28e-04   1.05e-05   1.01e+00  6.84e+04        1    6.85e-05    8.13e-0322  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）	从开始到当前迭代的总时间

从迭代日志中可以看到，迭代过程大概分成三个阶段：

1. 初期爆炸（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 算法鲁棒性的体现。
2. 中期快速下降（iter 7–15）：cost 从 5e3 快速降至 6.6；tr_ratio ≈ 1，说明局部线性模型有效；tr_radius 开始扩大（1e-2 → 1e2），进入高效牛顿阶段。
3. 后期精细收敛（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::expresidual[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 做的工作包含且不局限于：

1. 线性代数求解器的深度优化：能够自动选择最优线性求解器，对于小规模稠密问题，可以配置 DENSE_NORMAL_CHOLESKY，使用普通 Cholesky 求解，对于大规模稀疏问题（如 SLAM、Bundle Adjustment），可以配置 SPARSE_NORMAL_CHOLESKY，实现稀疏 Cholesky 求解。
2. 支持多种高度优化的稀疏线性代数库：例如SuiteSparse、Eigen等，能处理数十万参数、上百万残差项的 Bundle Adjustment 问题。
3. 雅可比矩阵的高效构建与存储：对于大规模非线性最小二乘问题，雅可比矩阵非常占用内存，例如 1M 观测数据 × 1K 待定参数，存储 double 类型的雅可比矩阵需要 8GB 的内存空间。而 Ceres 可以实现对雅可比矩阵按需计算 + 稀疏表示，用户只需提供每个残差块对局部参数的导数（即 AutoDiffCostFunction 返回的小 Jacobian 块），自动拼接成全局雅可比矩阵，只存储非零块。
4. 信任域策略的工业级鲁棒性：成熟的信任域管理，智能的 \(\lambda\) 调整逻辑（基于实际/预测下降比 \(\rho\)）；在噪声大、初值差、目标函数非凸的情况下仍能稳定收敛。
5. 参数约束与边界处理：支持参数边界约束。
6. 损失函数（Loss Function）抗外点能力：内置多种鲁棒损失函数，例如Huber、Cauchy、SoftLOne、Tukey 等，自动降低大残差的权重，抑制 outlier 影响。
7. 并行化与性能工程：通过CPU多线程或者GPU并行计算残差和雅克比矩阵，提升问题优化效率。

5. 总结

不得不说，要理解一个具有专业背景的库不是那么容易的事情，即使你理解底层算法的原理，但是工程实践又是另外一回事。只有将算法原理与工程实践融汇贯通，才能真正解决问题并创造价值。