The following is only used for illustration purposes. Ceres Solver ships with an optimized, production grade QuaternionManifold implementation.
作为一个具体的例子,我们来看看 Quaternions 要如何实现。四元数是嵌入到 R4 的三维流形,即它们的局部维数为 4,切空间维数为 3。下面的函数定义了四元数流形上的 Plus 和 Minus 运算。它假定四元数在内存中的布局为 [w,x,y,z],即实数或标量部分是第一坐标。
struct QuaternionFunctor {
template <typename T>
bool Plus(const T* x, const T* delta, T* x_plus_delta) const {
const T squared_norm_delta =
delta[0] * delta[0] + delta[1] * delta[1] + delta[2] * delta[2];
T q_delta[4];
if (squared_norm_delta > T(0.0)) {
T norm_delta = sqrt(squared_norm_delta);
const T sin_delta_by_delta = sin(norm_delta) / norm_delta;
q_delta[0] = cos(norm_delta);
q_delta[1] = sin_delta_by_delta * delta[0];
q_delta[2] = sin_delta_by_delta * delta[1];
q_delta[3] = sin_delta_by_delta * delta[2];
} else {
// We do not just use q_delta = [1,0,0,0] here because that is a
// constant and when used for automatic differentiation will
// lead to a zero derivative. Instead we take a first order
// approximation and evaluate it at zero.
q_delta[0] = T(1.0);
q_delta[1] = delta[0];
q_delta[2] = delta[1];
q_delta[3] = delta[2];
}
QuaternionProduct(q_delta, x, x_plus_delta);
return true;
}
template <typename T>
bool Minus(const T* y, const T* x, T* y_minus_x) const {
T minus_x[4] = {x[0], -x[1], -x[2], -x[3]};
T ambient_y_minus_x[4];
QuaternionProduct(y, minus_x, ambient_y_minus_x);
T u_norm = sqrt(ambient_y_minus_x[1] * ambient_y_minus_x[1] +
ambient_y_minus_x[2] * ambient_y_minus_x[2] +
ambient_y_minus_x[3] * ambient_y_minus_x[3]);
if (u_norm > 0.0) {
T theta = atan2(u_norm, ambient_y_minus_x[0]);
y_minus_x[0] = theta * ambient_y_minus_x[1] / u_norm;
y_minus_x[1] = theta * ambient_y_minus_x[2] / u_norm;
y_minus_x[2] = theta * ambient_y_minus_x[3] / u_norm;
} else {
We do not use [0,0,0] here because even though the value part is
a constant, the derivative part is not.
y_minus_x[0] = ambient_y_minus_x[1];
y_minus_x[1] = ambient_y_minus_x[2];
y_minus_x[2] = ambient_y_minus_x[3];
}
return true;
}
};
根据这一结构,使用自动微分的四元数流形可以构造为
Manifold* manifold = new AutoDiffManifold<QuaternionFunctor, 4, 3>;