这段代码实现了一个数值微分工具,基于有限差分法计算函数(曲线或向量场)的导数。它是 Sophus 库的一部分,用于在欧几里得空间中对标量函数、向量值函数或矩阵值函数进行数值微分。
功能和用途
曲线的数值微分:
支持对标量到标量、标量到向量或标量到矩阵的函数进行微分。
使用二阶中心差分公式近似导数。
向量场的数值微分:
计算从向量空间映射到向量空间的函数的雅可比矩阵。
针对不同输入维度(单变量或多变量)优化了实现。
接口封装:
curveNumDiff
:计算曲线(标量到向量/矩阵)的导数。
vectorFieldNumDiff
:计算向量场的雅可比矩阵。
提供了两个高层函数:
核心组件
1. 曲线微分:details::Curve
核心函数
num_diff
:
2. 向量场微分:details::VectorField
3. 高层封装
优点
高精度:
使用二阶中心差分公式,精度优于前向差分。
步长默认选择为机器精度的平方根,平衡了舍入误差和截断误差。
类型安全:
通过静态断言(
static_assert
)确保输入类型和返回类型符合预期。
扩展性:
模板化设计支持任意浮点类型(如
float
、double
)。兼容标量、向量和矩阵。
通用性:
支持标量函数和向量场,适用于各种数值优化或微分问题。
应用场景
数值优化:
在无解析导数时,通过有限差分计算梯度和雅可比矩阵。
微分验证:
用于验证解析导数的正确性。
计算机视觉和机器人学:
对轨迹曲线、速度场或运动模型的导数进行数值求解。
示例用法
/// @file // 文件注释/// 使用有限差分进行数值微分
#pragma once // 防止头文件被多次包含
#include <functional> // 引入functional库以使用std::function#include <type_traits> // 引入type_traits库以使用类型特征#include <utility> // 引入utility库以使用std::move
#include "types.hpp" // 引入自定义的types.hpp头文件
namespace Sophus { // 开始命名空间Sophus
namespace details { // 开始命名空间details
template <class Scalar>class Curve { // 定义模板类Curve public: template <class Fn> static auto num_diff(Fn curve, Scalar t, Scalar h) -> decltype(curve(t)) { // 静态成员函数,计算数值微分,返回类型与函数curve相同 using ReturnType = decltype(curve(t)); // 获取返回类型 static_assert(std::is_floating_point<Scalar>::value, "Scalar must be a floating point type."); // 确保Scalar是浮点类型 static_assert(IsFloatingPoint<ReturnType>::value, "ReturnType must be either a floating point scalar, " "vector or matrix."); // 确保返回类型是浮点标量、向量或矩阵return (curve(t + h) - curve(t - h)) / (Scalar(2) * h); // 计算并返回数值微分 }};
template <class Scalar, int N, int M>class VectorField { // 定义模板类VectorField public: static Eigen::Matrix<Scalar, N, M> num_diff( std::function<Sophus::Vector<Scalar, N>(Sophus::Vector<Scalar, M>)> vector_field, Sophus::Vector<Scalar, M> const& a, Scalar eps) { // 静态成员函数,计算矢量场的数值微分,返回N行M列的矩阵 static_assert(std::is_floating_point<Scalar>::value, "Scalar must be a floating point type."); // 确保Scalar是浮点类型 Eigen::Matrix<Scalar, N, M> J; // 定义矩阵J Sophus::Vector<Scalar, M> h; // 定义向量h h.setZero(); // 将向量h初始化为零 for (int i = 0; i < M; ++i) { h[i] = eps; // 将h的第i个元素设为eps J.col(i) = (vector_field(a + h) - vector_field(a - h)) / (Scalar(2) * eps); // 计算第i列的数值微分 h[i] = Scalar(0); // 将h的第i个元素重置为0 }return J; // 返回矩阵J }};
template <class Scalar, int N>class VectorField<Scalar, N, 1> { // 特化版本,当M为1时 public: static Eigen::Matrix<Scalar, N, 1> num_diff( std::function<Sophus::Vector<Scalar, N>(Scalar)> vector_field, Scalar const& a, Scalar eps) { // 静态成员函数,计算矢量场的数值微分,返回N行1列的矩阵 return details::Curve<Scalar>::num_diff(std::move(vector_field), a, eps); // 调用Curve类的num_diff方法 }};} // namespace details
/// 计算曲线在点t处的导数。// 这里,曲线是一个从标量到欧几里得空间的函数。因此,它返回标量、向量或矩阵。///template <class Scalar, class Fn>auto curveNumDiff(Fn curve, Scalar t, Scalar h = Constants<Scalar>::epsilonSqrt()) -> decltype(details::Curve<Scalar>::num_diff(std::move(curve), t, h)) { // 模板函数,计算曲线在点t处的导数 return details::Curve<Scalar>::num_diff(std::move(curve), t, h); // 调用Curve类的num_diff方法}
/// 计算矢量场在点a处的导数。// 这里,矢量场是一个从一个向量空间到另一个向量空间的函数。///template <class Scalar, int N, int M, class ScalarOrVector, class Fn>Eigen::Matrix<Scalar, N, M> vectorFieldNumDiff( Fn vector_field, ScalarOrVector const& a, Scalar eps = Constants<Scalar>::epsilonSqrt()) { // 模板函数,计算矢量场在点a处的导数 return details::VectorField<Scalar, N, M>::num_diff(std::move(vector_field), a, eps); // 调用VectorField类的num_diff方法}
} // namespace Sophus