make logNormalizationConstant return -log(k)

release/4.3a0
Varun Agrawal 2024-09-21 05:15:35 -04:00
parent 821b22f6f8
commit ecbf3d872e
11 changed files with 49 additions and 43 deletions

View File

@ -89,7 +89,7 @@ bool Conditional<FACTOR, DERIVEDCONDITIONAL>::CheckInvariants(
// normalization constant // normalization constant
const double error = conditional.error(values); const double error = conditional.error(values);
if (error < 0.0) return false; // prob_or_density is negative. if (error < 0.0) return false; // prob_or_density is negative.
const double expected = conditional.logNormalizationConstant() - error; const double expected = -(conditional.logNormalizationConstant() + error);
if (std::abs(logProb - expected) > 1e-9) if (std::abs(logProb - expected) > 1e-9)
return false; // logProb is not consistent with error return false; // logProb is not consistent with error
return true; return true;

View File

@ -34,9 +34,10 @@ namespace gtsam {
* `logProbability` is the main methods that need to be implemented in derived * `logProbability` is the main methods that need to be implemented in derived
* classes. These two methods relate to the `error` method in the factor by: * classes. These two methods relate to the `error` method in the factor by:
* probability(x) = k exp(-error(x)) * probability(x) = k exp(-error(x))
* where k is a normalization constant making \int probability(x) == 1.0, and * where k is a normalization constant making
* logProbability(x) = K - error(x) * \int probability(x) = \int k exp(-error(x)) == 1.0, and
* i.e., K = log(K). The normalization constant K is assumed to *not* depend * logProbability(x) = -(K + error(x))
* i.e., K = -log(k). The normalization constant k is assumed to *not* depend
* on any argument, only (possibly) on the conditional parameters. * on any argument, only (possibly) on the conditional parameters.
* This class provides a default logNormalizationConstant() == 0.0. * This class provides a default logNormalizationConstant() == 0.0.
* *
@ -163,8 +164,9 @@ namespace gtsam {
} }
/** /**
* All conditional types need to implement a log normalization constant to * All conditional types need to implement a
* make it such that error>=0. * (negative) log normalization constant
* to make it such that error>=0.
*/ */
virtual double logNormalizationConstant() const; virtual double logNormalizationConstant() const;

View File

@ -246,14 +246,15 @@ namespace gtsam {
double GaussianBayesNet::logNormalizationConstant() const { double GaussianBayesNet::logNormalizationConstant() const {
/* /*
normalization constant = 1.0 / sqrt((2*pi)^n*det(Sigma)) normalization constant = 1.0 / sqrt((2*pi)^n*det(Sigma))
logConstant = -0.5 * n*log(2*pi) - 0.5 * log det(Sigma) logConstant = -log(normalizationConstant)
= 0.5 * n*log(2*pi) + 0.5 * log(det(Sigma))
log det(Sigma)) = -2.0 * logDeterminant() log(det(Sigma)) = -2.0 * logDeterminant()
thus, logConstant = -0.5*n*log(2*pi) + logDeterminant() thus, logConstant = 0.5*n*log(2*pi) - logDeterminant()
BayesNet logConstant = sum(-0.5*n_i*log(2*pi) + logDeterminant_i()) BayesNet logConstant = sum(0.5*n_i*log(2*pi) - logDeterminant_i())
= sum(-0.5*n_i*log(2*pi)) + sum(logDeterminant_i()) = sum(0.5*n_i*log(2*pi)) - sum(logDeterminant_i())
= sum(-0.5*n_i*log(2*pi)) + bn->logDeterminant() = sum(0.5*n_i*log(2*pi)) - bn->logDeterminant()
*/ */
double logNormConst = 0.0; double logNormConst = 0.0;
for (const sharedConditional& cg : *this) { for (const sharedConditional& cg : *this) {

View File

@ -181,24 +181,24 @@ namespace gtsam {
/* ************************************************************************* */ /* ************************************************************************* */
// normalization constant = 1.0 / sqrt((2*pi)^n*det(Sigma)) // normalization constant = 1.0 / sqrt((2*pi)^n*det(Sigma))
// log = - 0.5 * n*log(2*pi) - 0.5 * log det(Sigma) // neg-log = 0.5 * n*log(2*pi) + 0.5 * log det(Sigma)
double GaussianConditional::logNormalizationConstant() const { double GaussianConditional::logNormalizationConstant() const {
constexpr double log2pi = 1.8378770664093454835606594728112; constexpr double log2pi = 1.8378770664093454835606594728112;
size_t n = d().size(); size_t n = d().size();
// Sigma = (R'R)^{-1}, det(Sigma) = det((R'R)^{-1}) = det(R'R)^{-1} // Sigma = (R'R)^{-1}, det(Sigma) = det((R'R)^{-1}) = det(R'R)^{-1}
// log det(Sigma) = -log(det(R'R)) = -2*log(det(R)) // log det(Sigma) = -log(det(R'R)) = -2*log(det(R))
// Hence, log det(Sigma)) = -2.0 * logDeterminant() // Hence, log det(Sigma)) = -2.0 * logDeterminant()
// which gives log = -0.5*n*log(2*pi) - 0.5*(-2.0 * logDeterminant()) // which gives neg-log = 0.5*n*log(2*pi) + 0.5*(-2.0 * logDeterminant())
// = -0.5*n*log(2*pi) + (0.5*2.0 * logDeterminant()) // = 0.5*n*log(2*pi) - (0.5*2.0 * logDeterminant())
// = -0.5*n*log(2*pi) + logDeterminant() // = 0.5*n*log(2*pi) - logDeterminant()
return -0.5 * n * log2pi + logDeterminant(); return 0.5 * n * log2pi - logDeterminant();
} }
/* ************************************************************************* */ /* ************************************************************************* */
// density = k exp(-error(x)) // density = 1/k exp(-error(x))
// log = log(k) - error(x) // log = -log(k) - error(x)
double GaussianConditional::logProbability(const VectorValues& x) const { double GaussianConditional::logProbability(const VectorValues& x) const {
return logNormalizationConstant() - error(x); return -logNormalizationConstant() - error(x);
} }
double GaussianConditional::logProbability(const HybridValues& x) const { double GaussianConditional::logProbability(const HybridValues& x) const {

View File

@ -133,8 +133,10 @@ namespace gtsam {
/// @{ /// @{
/** /**
* normalization constant = 1.0 / sqrt((2*pi)^n*det(Sigma)) * Return normalization constant in negative log space.
* log = - 0.5 * n*log(2*pi) - 0.5 * log det(Sigma) *
* normalization constant k = 1.0 / sqrt((2*pi)^n*det(Sigma))
* -log(k) = 0.5 * n*log(2*pi) + 0.5 * log det(Sigma)
*/ */
double logNormalizationConstant() const override; double logNormalizationConstant() const override;

View File

@ -257,13 +257,13 @@ double Gaussian::logDeterminant() const {
/* *******************************************************************************/ /* *******************************************************************************/
double Gaussian::logNormalizationConstant() const { double Gaussian::logNormalizationConstant() const {
// log(det(Sigma)) = -2.0 * logDetR // log(det(Sigma)) = -2.0 * logDetR
// which gives log = -0.5*n*log(2*pi) - 0.5*(-2.0 * logDetR()) // which gives neg-log = 0.5*n*log(2*pi) + 0.5*(-2.0 * logDetR())
// = -0.5*n*log(2*pi) + (0.5*2.0 * logDetR()) // = 0.5*n*log(2*pi) - (0.5*2.0 * logDetR())
// = -0.5*n*log(2*pi) + logDetR() // = 0.5*n*log(2*pi) - logDetR()
size_t n = dim(); size_t n = dim();
constexpr double log2pi = 1.8378770664093454835606594728112; constexpr double log2pi = 1.8378770664093454835606594728112;
// Get 1/log(\sqrt(|2pi Sigma|)) = -0.5*log(|2pi Sigma|) // Get -log(1/\sqrt(|2pi Sigma|)) = 0.5*log(|2pi Sigma|)
return -0.5 * n * log2pi + logDetR(); return 0.5 * n * log2pi - logDetR();
} }

View File

@ -274,8 +274,8 @@ namespace gtsam {
/** /**
* @brief Method to compute the normalization constant * @brief Method to compute the normalization constant
* for a Gaussian noise model k = \sqrt(1/|2πΣ|). * for a Gaussian noise model k = \sqrt(1/|2πΣ|).
* We compute this in the log-space for numerical accuracy, * We compute this in the negative log-space for numerical accuracy,
* thus returning log(k). * thus returning -log(k).
* *
* @return double * @return double
*/ */

View File

@ -76,11 +76,11 @@ TEST(GaussianBayesNet, Evaluate1) {
// the normalization constant 1.0/sqrt((2*pi*Sigma).det()). // the normalization constant 1.0/sqrt((2*pi*Sigma).det()).
// The covariance matrix inv(Sigma) = R'*R, so the determinant is // The covariance matrix inv(Sigma) = R'*R, so the determinant is
const double constant = sqrt((invSigma / (2 * M_PI)).determinant()); const double constant = sqrt((invSigma / (2 * M_PI)).determinant());
EXPECT_DOUBLES_EQUAL(log(constant), EXPECT_DOUBLES_EQUAL(-log(constant),
smallBayesNet.at(0)->logNormalizationConstant() + smallBayesNet.at(0)->logNormalizationConstant() +
smallBayesNet.at(1)->logNormalizationConstant(), smallBayesNet.at(1)->logNormalizationConstant(),
1e-9); 1e-9);
EXPECT_DOUBLES_EQUAL(log(constant), smallBayesNet.logNormalizationConstant(), EXPECT_DOUBLES_EQUAL(-log(constant), smallBayesNet.logNormalizationConstant(),
1e-9); 1e-9);
const double actual = smallBayesNet.evaluate(mean); const double actual = smallBayesNet.evaluate(mean);
EXPECT_DOUBLES_EQUAL(constant, actual, 1e-9); EXPECT_DOUBLES_EQUAL(constant, actual, 1e-9);

View File

@ -492,9 +492,10 @@ TEST(GaussianConditional, LogNormalizationConstant) {
VectorValues x; VectorValues x;
x.insert(X(0), Vector3::Zero()); x.insert(X(0), Vector3::Zero());
Matrix3 Sigma = I_3x3 * sigma * sigma; Matrix3 Sigma = I_3x3 * sigma * sigma;
double expectedLogNormalizingConstant = log(1 / sqrt((2 * M_PI * Sigma).determinant())); double expectedLogNormalizationConstant =
-log(1 / sqrt((2 * M_PI * Sigma).determinant()));
EXPECT_DOUBLES_EQUAL(expectedLogNormalizingConstant, EXPECT_DOUBLES_EQUAL(expectedLogNormalizationConstant,
conditional.logNormalizationConstant(), 1e-9); conditional.logNormalizationConstant(), 1e-9);
} }
@ -516,7 +517,7 @@ TEST(GaussianConditional, Print) {
" d = [ 20 40 ]\n" " d = [ 20 40 ]\n"
" mean: 1 elements\n" " mean: 1 elements\n"
" x0: 20 40\n" " x0: 20 40\n"
" logNormalizationConstant: -4.0351\n" " logNormalizationConstant: 4.0351\n"
"isotropic dim=2 sigma=3\n"; "isotropic dim=2 sigma=3\n";
EXPECT(assert_print_equal(expected, conditional, "GaussianConditional")); EXPECT(assert_print_equal(expected, conditional, "GaussianConditional"));
@ -531,7 +532,7 @@ TEST(GaussianConditional, Print) {
" S[x1] = [ -1 -2 ]\n" " S[x1] = [ -1 -2 ]\n"
" [ -3 -4 ]\n" " [ -3 -4 ]\n"
" d = [ 20 40 ]\n" " d = [ 20 40 ]\n"
" logNormalizationConstant: -4.0351\n" " logNormalizationConstant: 4.0351\n"
"isotropic dim=2 sigma=3\n"; "isotropic dim=2 sigma=3\n";
EXPECT(assert_print_equal(expected1, conditional1, "GaussianConditional")); EXPECT(assert_print_equal(expected1, conditional1, "GaussianConditional"));
@ -547,7 +548,7 @@ TEST(GaussianConditional, Print) {
" S[y1] = [ -5 -6 ]\n" " S[y1] = [ -5 -6 ]\n"
" [ -7 -8 ]\n" " [ -7 -8 ]\n"
" d = [ 20 40 ]\n" " d = [ 20 40 ]\n"
" logNormalizationConstant: -4.0351\n" " logNormalizationConstant: 4.0351\n"
"isotropic dim=2 sigma=3\n"; "isotropic dim=2 sigma=3\n";
EXPECT(assert_print_equal(expected2, conditional2, "GaussianConditional")); EXPECT(assert_print_equal(expected2, conditional2, "GaussianConditional"));
} }

View File

@ -55,7 +55,7 @@ TEST(GaussianDensity, FromMeanAndStddev) {
double expected1 = 0.5 * e.dot(e); double expected1 = 0.5 * e.dot(e);
EXPECT_DOUBLES_EQUAL(expected1, density.error(values), 1e-9); EXPECT_DOUBLES_EQUAL(expected1, density.error(values), 1e-9);
double expected2 = density.logNormalizationConstant()- 0.5 * e.dot(e); double expected2 = -(density.logNormalizationConstant() + 0.5 * e.dot(e));
EXPECT_DOUBLES_EQUAL(expected2, density.logProbability(values), 1e-9); EXPECT_DOUBLES_EQUAL(expected2, density.logProbability(values), 1e-9);
} }

View File

@ -810,9 +810,9 @@ TEST(NoiseModel, NonDiagonalGaussian)
TEST(NoiseModel, LogNormalizationConstant1D) { TEST(NoiseModel, LogNormalizationConstant1D) {
// Very simple 1D noise model, which we can compute by hand. // Very simple 1D noise model, which we can compute by hand.
double sigma = 0.1; double sigma = 0.1;
// For expected values, we compute 1/log(sqrt(|2πΣ|)) by hand. // For expected values, we compute -log(1/sqrt(|2πΣ|)) by hand.
// = -0.5*(log(2π) + log(Σ)) (since it is 1D) // = 0.5*(log(2π) + log(Σ)) (since it is 1D)
double expected_value = -0.5 * log(2 * M_PI * sigma * sigma); double expected_value = 0.5 * log(2 * M_PI * sigma * sigma);
// Gaussian // Gaussian
{ {
@ -839,7 +839,7 @@ TEST(NoiseModel, LogNormalizationConstant1D) {
auto noise_model = Unit::Create(1); auto noise_model = Unit::Create(1);
double actual_value = noise_model->logNormalizationConstant(); double actual_value = noise_model->logNormalizationConstant();
double sigma = 1.0; double sigma = 1.0;
expected_value = -0.5 * log(2 * M_PI * sigma * sigma); expected_value = 0.5 * log(2 * M_PI * sigma * sigma);
EXPECT_DOUBLES_EQUAL(expected_value, actual_value, 1e-9); EXPECT_DOUBLES_EQUAL(expected_value, actual_value, 1e-9);
} }
} }
@ -850,7 +850,7 @@ TEST(NoiseModel, LogNormalizationConstant3D) {
size_t n = 3; size_t n = 3;
// We compute the expected values just like in the LogNormalizationConstant1D // We compute the expected values just like in the LogNormalizationConstant1D
// test, but we multiply by 3 due to the determinant. // test, but we multiply by 3 due to the determinant.
double expected_value = -0.5 * n * log(2 * M_PI * sigma * sigma); double expected_value = 0.5 * n * log(2 * M_PI * sigma * sigma);
// Gaussian // Gaussian
{ {
@ -879,7 +879,7 @@ TEST(NoiseModel, LogNormalizationConstant3D) {
auto noise_model = Unit::Create(3); auto noise_model = Unit::Create(3);
double actual_value = noise_model->logNormalizationConstant(); double actual_value = noise_model->logNormalizationConstant();
double sigma = 1.0; double sigma = 1.0;
expected_value = -0.5 * n * log(2 * M_PI * sigma * sigma); expected_value = 0.5 * n * log(2 * M_PI * sigma * sigma);
EXPECT_DOUBLES_EQUAL(expected_value, actual_value, 1e-9); EXPECT_DOUBLES_EQUAL(expected_value, actual_value, 1e-9);
} }
} }