1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
|
// Copyright (c) 2007, 2008, 2009 libmv authors.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#ifndef LIBMV_NUMERIC_DERIVATIVE_H
#define LIBMV_NUMERIC_DERIVATIVE_H
#include <cmath>
#include "libmv/numeric/numeric.h"
#include "libmv/logging/logging.h"
namespace libmv {
// Numeric derivative of a function.
// TODO(keir): Consider adding a quadratic approximation.
enum NumericJacobianMode {
CENTRAL,
FORWARD,
};
template<typename Function, NumericJacobianMode mode = CENTRAL>
class NumericJacobian {
public:
typedef typename Function::XMatrixType Parameters;
typedef typename Function::XMatrixType::RealScalar XScalar;
typedef typename Function::FMatrixType FMatrixType;
typedef Matrix<typename Function::FMatrixType::RealScalar,
Function::FMatrixType::RowsAtCompileTime,
Function::XMatrixType::RowsAtCompileTime>
JMatrixType;
NumericJacobian(const Function &f) : f_(f) {}
// TODO(keir): Perhaps passing the jacobian back by value is not a good idea.
JMatrixType operator()(const Parameters &x) {
// Empirically determined constant.
Parameters eps = x.array().abs() * XScalar(1e-5);
// To handle cases where a paremeter is exactly zero, instead use the mean
// eps for the other dimensions.
XScalar mean_eps = eps.sum() / eps.rows();
if (mean_eps == XScalar(0)) {
// TODO(keir): Do something better here.
mean_eps = 1e-8; // ~sqrt(machine precision).
}
// TODO(keir): Elimininate this needless function evaluation for the
// central difference case.
FMatrixType fx = f_(x);
const int rows = fx.rows();
const int cols = x.rows();
JMatrixType jacobian(rows, cols);
Parameters x_plus_delta = x;
for (int c = 0; c < cols; ++c) {
if (eps(c) == XScalar(0)) {
eps(c) = mean_eps;
}
x_plus_delta(c) = x(c) + eps(c);
jacobian.col(c) = f_(x_plus_delta);
XScalar one_over_h = 1 / eps(c);
if (mode == CENTRAL) {
x_plus_delta(c) = x(c) - eps(c);
jacobian.col(c) -= f_(x_plus_delta);
one_over_h /= 2;
} else {
jacobian.col(c) -= fx;
}
x_plus_delta(c) = x(c);
jacobian.col(c) = jacobian.col(c) * one_over_h;
}
return jacobian;
}
private:
const Function &f_;
};
template<typename Function, typename Jacobian>
bool CheckJacobian(const Function &f, const typename Function::XMatrixType &x) {
Jacobian j_analytic(f);
NumericJacobian<Function> j_numeric(f);
typename NumericJacobian<Function>::JMatrixType J_numeric = j_numeric(x);
typename NumericJacobian<Function>::JMatrixType J_analytic = j_analytic(x);
LG << J_numeric - J_analytic;
return true;
}
} // namespace libmv
#endif // LIBMV_NUMERIC_DERIVATIVE_H
|
