SensorsCalibration（OpenCalib）：Ceres 解析导数（Analytical Derivatives）建模与性能优化实战指南

• 自动驾驶
• 计算机视觉

本篇基于 OpenCalib 多传感器标定工具包中随仓库内置的 Ceres Solver 官方文档章节 Analytic Derivatives，完整讲解如何以 Rat43 曲线拟合为例手工推导并实现CostFunction的解析雅可比，如何通过缓存中间子表达式将求值耗时降低约 2.8 倍，以及何时值得投入手工解析导数。读完后你可以：独立实现一个带解析雅可比的 Ceres 代价函数、理解Evaluate接口的参数布局与雅可比存储约定，并对照仓库中 lidar2car 模块 Fast-LOAM 风格的解析残差实现（EdgeAnalyticCostFunction、SurfNormAnalyticCostFunction）验证这些原则在激光雷达-车体标定中的真实落地。

Ceres 的三种导数策略：解析导数在其中的定位

Ceres 是一个基于梯度的非线性最小二乘求解器，定义目标函数及其雅可比矩阵是用户建模的核心任务。随仓库内置的 derivatives.rst 明确给出三条路线：

1. 解析导数（Analytical）：用户自己推导并实现导数，写入CostFunction；
2. 数值导数（Numerical）：Ceres 用有限差分计算；
3. 自动微分（Automatic）：Ceres 通过 C++ 模板与运算符重载自动计算解析级导数。

文档给出的高层建议很直接：优先用自动微分；某些情况下值得用解析导数；避免数值导数（仅在对接外部库等最后手段时使用）。解析导数章节正是对"什么情况值得手推"的完整论证，本文即围绕它展开。

问题建模：Rat43 曲线拟合

文档选取 NIST 非线性拟合数据集中的 Rat43 问题：给定数据 ${x_i, y_i},\ i=1,...,n$，寻找参数 $b_1, b_2, b_3, b_4$ 使曲线

$$y = \frac{b_1}{(1+e^{b_2-b_3x})^{1/b_4}}$$

最优地拟合数据。其目标函数为

$$ \begin{align} E(b_1, b_2, b_3, b_4) &= \sum_i f^2(b_1, b_2, b_3, b_4 ; x_i, y_i)\ &= \sum_i \left(\frac{b_1}{(1+e^{b_2-b_3x_i})^{1/b_4}} - y_i\right)^2 \end{align} $$

原文脚注特别强调：拟合"最优"与否取决于用于衡量拟合质量的目标函数选择，而后者取决于生成观测数据的噪声过程。仅当噪声为高斯分布时，最小化平方差之和才是正确的做法，此时参数最优解即为最大似然估计（MLE）。这一提醒在标定场景中同样关键——它解释了为什么 OpenCalib 的 lidar2car 优化在残差块上加了鲁棒核函数（见后文HuberLoss）。

用 Ceres 求解该问题，需要定义一个CostFunction，计算给定 $x, y$ 下的残差 $f$ 及其对 $b_1, b_2, b_3, b_4$ 的偏导数。

手工推导雅可比：四个偏导公式

对 $f(b_1,b_2,b_3,b_4;x,y)=\dfrac{b_1}{(1+e^{b_2-b_3x})^{1/b_4}}-y$ 做初等微分，可得：

$$ \begin{align} D_1 f &= \frac{1}{(1+e^{b_2-b_3x})^{1/b_4}}\ D_2 f &= \frac{-b_1e^{b_2-b_3x}}{b_4(1+e^{b_2-b_3x})^{1/b_4 + 1}} \ D_3 f &= \frac{b_1xe^{b_2-b_3x}}{b_4(1+e^{b_2-b_3x})^{1/b_4 + 1}} \ D_4 f &= \frac{b_1 \log(1+e^{b_2-b_3x})}{b_4^2(1+e^{b_2-b_3x})^{1/b_4}} \end{align} $$

第一版实现：直译公式的Rat43Analytic

拿到导数后，文档给出的直接实现如下（继承ceres::SizedCostFunction<1,4>，表示 1 个残差、1 个 4 维参数块）：

class Rat43Analytic : public SizedCostFunction<1,4> { public: Rat43Analytic(const double x, const double y) : x_(x), y_(y) {} virtual ~Rat43Analytic() {} virtual bool Evaluate(double const* const* parameters, double* residuals, double** jacobians) const { const double b1 = parameters[0][0]; const double b2 = parameters[0][1]; const double b3 = parameters[0][2]; const double b4 = parameters[0][3]; residuals[0] = b1 * pow(1 + exp(b2 - b3 * x_), -1.0 / b4) - y_; if (!jacobians) return true; double* jacobian = jacobians[0]; if (!jacobian) return true; jacobian[0] = pow(1 + exp(b2 - b3 * x_), -1.0 / b4); jacobian[1] = -b1 * exp(b2 - b3 * x_) * pow(1 + exp(b2 - b3 * x_), -1.0 / b4 - 1) / b4; jacobian[2] = x_ * b1 * exp(b2 - b3 * x_) * pow(1 + exp(b2 - b3 * x_), -1.0 / b4 - 1) / b4; jacobian[3] = b1 * log(1 + exp(b2 - b3 * x_)) * pow(1 + exp(b2 - b3 * x_), -1.0 / b4) / (b4 * b4); return true; } private: const double x_; const double y_; };

这段代码能正确工作，但"冗长、难读、冗余很多"——1 + exp(b2 - b3 * x_)与pow(...)被反复计算。

Evaluate接口的关键约定

结合仓库内置的 cost_function.h 头文件注释，实现时需要注意四条约定：

• 参数布局：parameters[i]指向第 $i$ 个参数块，其长度为parameter_block_sizes_[i]，各参数块顺序与AddResidualBlock时登记的一致；
• 雅可比布局：jacobians[i]是num_residuals * parameter_block_sizes_[i]个元素的行主序数组，即jacobians[i][r * size + c] = d residual[r] / d parameters[i][c]；
• NULL 判定的双重含义：jacobians == NULL表示本次只算残差；jacobians[i] == NULL表示该参数块的雅可比不需要（例如参数被设为常量）。所以先判!jacobians再判!jacobian是标准写法；
• 返回值的双重用途：返回false可表示数值失败，也可用来施加参数约束——只要初始值可行，在约束被违反时返回false即可阻止求解器进入不可行区域。

而 sized_cost_function.h 显示，SizedCostFunction<kNumResiduals, Ns...>的价值在于构造函数中一次性把残差数与参数块尺寸写入基类（set_num_residuals/mutable_parameter_block_sizes），省去了动态尺寸基类的样板代码。

第二版实现：缓存子表达式，提速 2.8 倍

把重复出现的子表达式提为局部变量后，得到优化版：

class Rat43AnalyticOptimized : public SizedCostFunction<1,4> { public: Rat43AnalyticOptimized(const double x, const double y) : x_(x), y_(y) {} virtual ~Rat43AnalyticOptimized() {} virtual bool Evaluate(double const* const* parameters, double* residuals, double** jacobians) const { const double b1 = parameters[0][0]; const double b2 = parameters[0][1]; const double b3 = parameters[0][2]; const double b4 = parameters[0][3]; const double t1 = exp(b2 - b3 * x_); const double t2 = 1 + t1; const double t3 = pow(t2, -1.0 / b4); residuals[0] = b1 * t3 - y_; if (!jacobians) return true; double* jacobian = jacobians[0]; if (!jacobian) return true; const double t4 = pow(t2, -1.0 / b4 - 1); jacobian[0] = t3; jacobian[1] = -b1 * t1 * t4 / b4; jacobian[2] = -x_ * jacobian[1]; jacobian[3] = b1 * log(t2) * t3 / (b4 * b4); return true; } private: const double x_; const double y_; };

注意jacobian[2] = -x_ * jacobian[1]这种复用：$D_3$ 公式本身就是 $D_2$ 公式乘以 $-x$，缓存后连pow调用都省掉了。

性能对比

文档给出的实测耗时：

优化版比直译版快2.8 倍，文档明确指出这种量级的差距"并不罕见"，"要把解析导数的性能榨干，通常需要做公共子表达式优化"。

同目录下的 automatic_derivatives.rst 给出了同一 Rat43 问题的完整横向对比，可以佐证解析导数的收益边界：

自动微分（129 ns）只比手工优化后的解析版（92 ns）慢约 40%，却省去了全部手推工作——这正是前文"优先自动微分"建议的量化依据。

什么时候值得用解析导数

文档给出了五条适用判据，原文完整继承如下：

1. 表达式很简单（例如以线性项为主）；
2. 可以借助计算机代数系统（Maple、Mathematica 或 SymPy 等）对目标函数做符号微分，自动生成求值用的 C++ 代码；
3. 性能是最高优先级，且残差项存在可利用的代数结构，能取得优于自动微分的性能。但文档同时警告：榨干解析导数的性能需要不小的工作量，走这条路之前应先测量"雅可比求值占总求解时间的比例"，并牢记Amdahl 定律；
4. 没有其他办法计算导数，典型例子是求多项式根对变量的导数：

$$a_3(x,y)z^3 + a_2(x,y)z^2 + a_1(x,y)z + a_0(x,y) = 0$$

求 $z$ 关于 $x, y$ 的导数必须借助隐函数定理（Inverse Function Theorem）； 5. 你热爱链式法则，并且享受手推代数（原文如此）。

仓库实证：lidar2car 中真正的解析残差

OpenCalib 仓库里恰好有一个"判据 4 变体"的真实案例：SensorX2car/lidar2car的 Fast-LOAM 风格里程计/标定模块需要把 SE(3) 位姿上的点变换求导，这类几何导数无法用模板自动微分"自动"得到，必须手写。其残差实现位于 lidarOptimization.cpp，声明位于 lidarOptimization.h：

class EdgeAnalyticCostFunction : public ceres::SizedCostFunction<1, 7> { public: EdgeAnalyticCostFunction(Eigen::Vector3d curr_point_, Eigen::Vector3d last_point_a_, Eigen::Vector3d last_point_b_); virtual bool Evaluate(double const *const *parameters, double *residuals, double **jacobians) const; // ... };

其Evaluate的写法与 Rat43 文档的约定完全一致（参见 lidarOptimization.cpp）：

• 用Eigen::Map把 7 维参数块解释为四元数 + 平移：q_last_curr取parameters[0]，t_last_curr取parameters[0] + 4；
• 残差是"投影点到直线段的距离"：residuals[0] = nu.norm()/de_norm；
• 雅可比部分先if(jacobians != NULL)再if(jacobians[0] != NULL)双重判空；
• 通过预构建 3×6 矩阵dp_by_se3（左上角为-skew(lp)，右上角为单位阵）把链式法则"物化"成一次块乘法，再经Eigen::Map<Matrix<double,1,7,RowMajor>>以行主序写入jacobians[0]——这正是cost_function.h所述行主序约定的实际运用。

SurfNormAnalyticCostFunction（点到拟合平面的距离残差，lidarOptimization.cpp）同理：残差为plane_unit_norm.dot(point_w) + negative_OA_dot_norm，雅可比为plane_unit_norm^T * dp_by_se3。

配套的PoseSE3Parameterization（继承ceres::LocalParameterization）则示范了另一个细节：7 维四元数全局参数块需要自定义LocalParameterization把 6 维李代数增量映射回 7 维空间，并同时手写其 Jacobian（ComputeJacobian）。这是"局部参数化 + 解析导数"组合在姿态标定中的完整拼图。

残差块如何进入 Problem：调用链与求解器配置

在 odomEstimationClass.cpp 中可以追踪完整调用链：

• addEdgeCostFactor/addSurfCostFactor逐点构造EdgeAnalyticCostFunction/SurfNormAnalyticCostFunction，并通过problem.AddResidualBlock(cost_function, loss_function, parameters)注册（odomEstimationClass.cpp、L200-L201）；
• 参数块注册时挂上局部参数化：problem.AddParameterBlock(parameters, 7, new PoseSE3Parameterization())（L57）；
• 鲁棒核函数：ceres::HuberLoss(0.1)（L53）——呼应前文脚注中"噪声非理想高斯时应换目标函数/鲁棒核"的提醒；
• 求解器选项：linear_solver_type = ceres::DENSE_QR、max_num_iterations = 4、check_gradients = false（L62-L67）。

从源码结构看，check_gradients = false的选择与解析导数的场景吻合：手写雅可比在开发阶段应打开梯度检查验证正确性，确认无误后在实时性敏感的激光雷达里程计内环中关闭，以换取文档反复强调的"Amdahl 定律式"时间预算意识。

小结：解析导数的工程决策清单

结合本文与仓库实例，可提炼出一条可操作的决策清单：

1. 先写残差，导数用自动微分起步：Ceres 的cost_function.h头注释本身就建议"优先定义模板残差 functor 配合AutoDiffCostFunction"，通常代码更短、执行更快；
2. 确认瓶颈再优化：按 Amdahl 定律测量雅可比求值占比，占比小就不必手工优化；
3. 必须手写时，按 Rat43 范式走：先直译公式保证正确（开check_gradients验证），再引入t1/t2/...式子表达式缓存复用（参考Rat43AnalyticOptimized的 2.8 倍提速）；
4. 严格遵守Evaluate约定：行主序雅可比、两级 NULL 判空、false返回值可用于约束；
5. 几何类导数（SE(3)/SO(3)）无自动微分可依赖，参考 lidar2car 中EdgeAnalyticCostFunction+PoseSE3Parameterization的组合方式，把 skew 矩阵与链式法则块矩阵化，兼顾正确性与性能。

相关源码入口索引：

• 文档本体：analytical_derivatives.rst
• 三策略总览：derivatives.rst；自动微分对照表：automatic_derivatives.rst
• 接口定义：cost_function.h、sized_cost_function.h
• 仓库内真实解析残差：lidarOptimization.h、lidarOptimization.cpp、odomEstimationClass.cpp

• 自动驾驶
• 计算机视觉