Ceres Solver 官方教程學習筆記（Ⅹ）——自動微分法介面Interfacing with Automatic Differentiation
阿新 • • 發佈：2019-01-22
在成本函式的有一個顯式表示式的情況下，自動微分演算法很容易使用。但有時候這不太現實。通常程式都需要與外部的程式或資料進行互動。在這一章中，我們將考慮幾種不同的方法來在這些特殊情況下使用自動微分法。
現在我們考慮一個優化問題。尋找引數和使
在這裡，是一個二維旋轉矩陣，依賴於旋轉角度。是一個二維向量，表示位移。是一個外部畸變函式。
首先考慮這種情況，我們有一個模板函式TemplatedComputeDistortion 可以計算函式 。然後，對應的殘差的程式碼實現如下：
template <typename T> T TemplatedComputeDistortion(const T r2) {
const double k1 = 0.0082;
const double k2 = 0.000023;
return 1.0 + k1 * y2 + k2 * r2 * r2;
}
struct Affine2DWithDistortion {
Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
x[0] = x_in[0 但現在讓我們考慮三種特殊情況。如果函式不能直接用於自動區分，常見的比如：
1. 是一個非模板求值函式。
2. 是一個可以計算值和微分的非模板函式。
3. 是一個待插值的值表函式。
下面我們依次探討這些情況。
返回值的非模板函式
假設我們有一個函式，其宣告如下：
double ComputeDistortionValue(double r2);函式的具體內部實現不重要。將這個函式對接到Affine2DWithDistortion中需要三步：
- 把
ComputeDistortionValue封裝成函式ComputeDistortionValueFunctor。 - 對
ComputeDistortionValueFunctor使用NumericDiffCostFunction進行數值微分，從而建立CostFunction. - 使用
CostFunctionToFunctor封裝CostFunction。封裝後得到一個帶有模板化操作符operator()的函式。這個操作符operator()方法可以將NumericDiffCostFunction計算出的雅可比矩陣變成Jet物件。
以上步驟的具體程式碼如下：
struct ComputeDistortionValueFunctor { // 第一步
bool operator()(const double* r2, double* value) const {
*value = ComputeDistortionValue(r2[0]);
return true;
}
};
struct Affine2DWithDistortion {
Affine2DWithDistortion(const double x_in[2], const double y_in[2]) { // 建構函式，在初始化過程中完成轉化
x[0] = x_in[0];
x[1] = x_in[1];
y[0] = y_in[0];
y[1] = y_in[1];
compute_distortion.reset(new ceres::CostFunctionToFunctor<1, 1>( // 第三步（外層函式）
new ceres::NumericDiffCostFunction<ComputeDistortionValueFunctor, // 第二步（內層函式）
ceres::CENTRAL,
1,
1>(
new ComputeDistortionValueFunctor)));
}
template <typename T>
bool operator()(const T* theta, const T* t, T* residuals) const {
const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
const T r2 = q_0 * q_0 + q_1 * q_1;
T f;
(*compute_distortion)(&r2, &f); // 變成一個模板類compute_distortion
residuals[0] = y[0] - f * q_0;
residuals[1] = y[1] - f * q_1;
return true;
}
double x[2];
double y[2];
std::unique_ptr<ceres::CostFunctionToFunctor<1, 1> > compute_distortion;//先定義
};返回值和微分的非模板函式
現在假設我們有一個函式ComputeDistortionValue，可以得到它的值，並且可以根據需要獲取其雅可比矩陣。其函式宣告如下：
void ComputeDistortionValueAndJacobian(double r2,
double* value,
double* jacobian);同樣，函式的實際實現並不重要。處理這個函式需要兩步：
與第一種情況相比，這裡直接可以求出雅可比矩陣，所以可以直接構建CostFunction。而不需要先準備Functor
- 把
ComputeDistortionValueAndJacobian封裝到一個CostFunction物件內。這個CostFunction物件我們稱為ComputeDistortionFunction。 - 用
CostFunctionToFunctor封裝剛剛得到的ComputeDistortionFunction物件。得到一個帶有模板操作符operator()方法的Functor，它將由NumericDiffCostFunction計算出的雅可比矩陣變成適合Jet物件。
程式碼如下：
class ComputeDistortionFunction : public ceres::SizedCostFunction<1, 1> { // 第一步
public:
virtual bool Evaluate(double const* const* parameters,
double* residuals,
double** jacobians) const {
if (!jacobians) { // 如果不需要雅可比矩陣
ComputeDistortionValueAndJacobian(parameters[0][0], residuals, NULL);
} else { // 如果需要雅可比矩陣
ComputeDistortionValueAndJacobian(parameters[0][0], residuals, jacobians[0]);
}
return true;
}
};
struct Affine2DWithDistortion {
Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {// 建構函式，在初始化過程中完成轉化
x[0] = x_in[0];
x[1] = x_in[1];
y[0] = y_in[0];
y[1] = y_in[1];
compute_distortion.reset( // 第二步
new ceres::CostFunctionToFunctor<1, 1>(new ComputeDistortionFunction));
}
template <typename T>
bool operator()(const T* theta,
const T* t,
T* residuals) const {
const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
const T r2 = q_0 * q_0 + q_1 * q_1;
T f;
(*compute_distortion)(&r2, &f); // 變成一個模板類compute_distortion
residuals[0] = y[0] - f * q_0;
residuals[1] = y[1] - f * q_1;
return true;
}
double x[2];
double y[2];
std::unique_ptr<ceres::CostFunctionToFunctor<1, 1> > compute_distortion; //先定義
};定義為值表的函式
最後一個例子是，函式是一個被定義在區間的值表，每個整數都有一個對應的輸出值。其本質就是一個向量。
vector<double> distortion_values;有很多方法可以插入一個值表。也許最簡單和最常用的方法是線性插值。但在這裡線性插值不是個好辦法，因為插值函式在抽樣點處是不可微的。
另一個簡單但是效能優異的可微插值方法是 Cubic Hermite Spline(中文埃爾米特插值) Ceres提供Cubic和Bi-Cubic插值的整個流程，並且可以很方便的應用自動微分演算法。
使用 Cubic插值，首先需要構造一個Grid1D物件來包裝值表，然後構造一個CubicInterpolator·物件來使用它。程式碼如下：
struct Affine2DWithDistortion {
Affine2DWithDistortion(const double x_in[2],
const double y_in[2],
const std::vector<double>& distortion_values) {
x[0] = x_in[0];
x[1] = x_in[1];
y[0] = y_in[0];
y[1] = y_in[1];
grid.reset(new ceres::Grid1D<double, 1>(
&distortion_values[0], 0, distortion_values.size()));
compute_distortion.reset(
new ceres::CubicInterpolator<ceres::Grid1D<double, 1> >(*grid));
}
template <typename T>
bool operator()(const T* theta,
const T* t,
T* residuals) const {
const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
const T r2 = q_0 * q_0 + q_1 * q_1;
T f;
compute_distortion->Evaluate(r2, &f);
residuals[0] = y[0] - f * q_0;
residuals[1] = y[1] - f * q_1;
return true;
}
double x[2];
double y[2];
std::unique_ptr<ceres::Grid1D<double, 1> > grid;//先定義
std::unique_ptr<ceres::CubicInterpolator<ceres::Grid1D<double, 1> > > compute_distortion;//先定義
};在上面的例子中，我們使用了Grid1D和CubicInterpolator來插入一個一維的值表。Grid2D``與CubicInterpolator相結合可以用於插入二維值表。注意，無論是Grid1D還是Grid2D“`都不侷限於標量值函式，它們也與向量值函式一起工作。