saberbud/ceres-solver
1
0
Fork 0
mirror of https://github.com/ceres-solver/ceres-solver.git synced 2026-08-30 09:00:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

support log1p and expm1 jet

Browse Source 
Currently, it is not possible to accurately evaluate the derivative of
d/dx log(1 + x) under all circumstances and significant deviations from
the actual derivative d/dx log1p(x) can occur. This changeset introduces
the necessary Jet overload and its inverse, expm1.

Change-Id: Ifcf88f6d684f61ba86bbe49f0d551b703f34ad0d
This commit is contained in:
Sergiu Deitsch
2021-11-14 17:28:30 +01:00
parent a668cabbcf
commit 2fba61434b
2 changed files with 71 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/ceres/jet.h
+15
View File
@@ -394,6 +394,7 @@ using std::erf;
using std::erfc;
using std::exp;
using std::exp2;
using std::expm1;
using std::floor;
using std::fmax;
using std::fmin;
@@ -403,6 +404,7 @@ using std::isinf;
using std::isnan;
using std::isnormal;
using std::log;
using std::log1p;
using std::log2;
using std::norm;
using std::pow;
@@ -471,6 +473,13 @@ inline Jet<T, N> log(const Jet<T, N>& f) {
return Jet<T, N>(log(f.a), f.v * a_inverse);
}
// log1p(a + h) ~= log1p(a) + h / (1 + a)
template <typename T, int N>
inline Jet<T, N> log1p(const Jet<T, N>& f) {
const T a_inverse = T(1.0) / (T(1.0) + f.a);
return Jet<T, N>(log1p(f.a), f.v * a_inverse);
}
// exp(a + h) ~= exp(a) + exp(a) h
template <typename T, int N>
inline Jet<T, N> exp(const Jet<T, N>& f) {
@@ -478,6 +487,12 @@ inline Jet<T, N> exp(const Jet<T, N>& f) {
return Jet<T, N>(tmp, tmp * f.v);
}
// expm1(a + h) ~= expm1(a) + exp(a) h
template <typename T, int N>
inline Jet<T, N> expm1(const Jet<T, N>& f) {
return Jet<T, N>(expm1(f.a), exp(f.a) * f.v);
}
// sqrt(a + h) ~= sqrt(a) + h / (2 sqrt(a))
template <typename T, int N>
inline Jet<T, N> sqrt(const Jet<T, N>& f) {
internal/ceres/jet_test.cc
+56
View File
@@ -128,6 +128,22 @@ TEST(Jet, Jet) {
ExpectJetsClose(w, x);
}
{ // Check that expm1(log1p(x)) == x.
J z = expm1(x);
J w = log1p(z);
VL << "z = " << z;
VL << "w = " << w;
ExpectJetsClose(w, x);
}
{ // Check that log1p(expm1(x)) == x.
J z = log1p(x);
J w = expm1(z);
VL << "z = " << z;
VL << "w = " << w;
ExpectJetsClose(w, x);
}
{ // Check that (x * y) / x == y.
J z = x * y;
J w = z / x;
@@ -631,6 +647,46 @@ TEST(Jet, Jet) {
NumericalTest("cbrt", cbrt<double, 2>, 1e-5);
NumericalTest("cbrt", cbrt<double, 2>, 1.0);
{ // Check that log1p(x) == log(1 + x)
J z = log1p(x);
J w = log(J{1} + x);
VL << "z = " << z;
VL << "w = " << w;
ExpectJetsClose(z, w);
}
{ // Check that log1p(x) does not loose precision for small x
J x = MakeJet(1e-16, 1e-8, 1e-4);
J z = log1p(x);
J w = MakeJet(9.9999999999999998e-17, 1e-8, 1e-4);
VL << "z = " << z;
VL << "w = " << w;
ExpectJetsClose(z, w);
// log(1 + x) collapes to 0
J v = log(J{1} + x);
EXPECT_TRUE(v.a == 0);
}
{ // Check that expm1(x) == exp(x) - 1
J z = expm1(x);
J w = exp(x) - J{1};
VL << "z = " << z;
VL << "w = " << w;
ExpectJetsClose(z, w);
}
{ // Check that expm1(x) does not loose precision for small x
J x = MakeJet(9.9999999999999998e-17, 1e-8, 1e-4);
J z = expm1(x);
J w = MakeJet(1e-16, 1e-8, 1e-4);
VL << "z = " << z;
VL << "w = " << w;
ExpectJetsClose(z, w);
// exp(x) - 1 collapes to 0
J v = exp(x) - J{1};
EXPECT_TRUE(v.a == 0);
}
{ // Check that exp2(x) == exp(x * log(2))
J z = exp2(x);
J w = exp(x * log(2.0));