KL Divergence (Kullback-Leibler Divergence)

1
#include <iostream>
2
#include <vector>
3
#include <cmath>
4
#include <numeric>
5
#include <iomanip>
6
#include <stdexcept>
7
8
// Normalize a nonnegative vector to sum to 1 (probability simplex)
9
std::vector<double> normalize(const std::vector<double>& v) {
10
    double sum = std::accumulate(v.begin(), v.end(), 0.0);
11
    if (sum <= 0.0) throw std::invalid_argument("Vector must have positive sum to normalize.");
12
    std::vector<double> out(v.size());
13
    for (size_t i = 0; i < v.size(); ++i) out[i] = v[i] / sum;
14
    return out;
15
}
16
17
// Additive smoothing: add eps to all entries, then renormalize
18
std::vector<double> smooth_and_normalize(const std::vector<double>& v, double eps) {
19
    if (eps < 0.0) throw std::invalid_argument("eps must be nonnegative");
20
    std::vector<double> w = v;
21
    for (double &x : w) x = std::max(0.0, x) + eps; // clamp negatives to 0, then add eps
22
    return normalize(w);
23
}
24
25
// Shannon entropy H(P) in nats
26
double entropy(const std::vector<double>& P) {
27
    double H = 0.0;
28
    for (double p : P) {
29
        if (p > 0.0) H -= p * std::log(p);
30
        // terms with p==0 contribute 0
31
    }
32
    return H;
33
}
34
35
// Cross-entropy H(P,Q) in nats; returns +inf if any Q_i==0 with P_i>0
36
double cross_entropy(const std::vector<double>& P, const std::vector<double>& Q) {
37
    if (P.size() != Q.size()) throw std::invalid_argument("Size mismatch");
38
    double H = 0.0;
39
    for (size_t i = 0; i < P.size(); ++i) {
40
        if (P[i] > 0.0) {
41
            if (Q[i] <= 0.0) return INFINITY; // undefined; infinite penalty
42
            H -= P[i] * std::log(Q[i]);
43
        }
44
    }
45
    return H;
46
}
47
48
// KL divergence D_KL(P||Q) in nats; returns +inf if Q has zero where P>0 (strict)
49
double kl_divergence(const std::vector<double>& P, const std::vector<double>& Q) {
50
    if (P.size() != Q.size()) throw std::invalid_argument("Size mismatch");
51
    double D = 0.0;
52
    for (size_t i = 0; i < P.size(); ++i) {
53
        double p = P[i], q = Q[i];
54
        if (p == 0.0) continue; // term is 0
55
        if (q <= 0.0) return INFINITY; // infinite KL
56
        D += p * (std::log(p) - std::log(q));
57
    }
58
    return D;
59
}
60
61
int main() {
62
    std::cout << std::fixed << std::setprecision(6);
63
64
    // Example distributions (possibly unnormalized counts)
65
    std::vector<double> P_raw{0.7, 0.2, 0.1};
66
    std::vector<double> Q1_raw{0.6, 0.3, 0.1};
67
    std::vector<double> Q2_raw{0.0, 0.5, 0.5}; // zero where P>0 -> infinite forward KL (strict)
68
69
    // Normalize (or smooth+normalize for safety)
70
    auto P = normalize(P_raw);
71
    auto Q1 = normalize(Q1_raw);
72
    auto Q2 = normalize(Q2_raw);
73
74
    // Entropy and cross-entropy
75
    double HP = entropy(P);
76
    double HPQ1 = cross_entropy(P, Q1);
77
    double HPQ2 = cross_entropy(P, Q2); // +inf expected
78
79
    // KL computations
80
    double D_P_Q1 = kl_divergence(P, Q1);
81
    double D_P_Q2 = kl_divergence(P, Q2); // +inf expected
82
    double D_Q1_P = kl_divergence(Q1, P); // reverse KL
83
84
    std::cout << "H(P): " << HP << " nats\n";
85
    std::cout << "H(P,Q1): " << HPQ1 << " nats, KL(P||Q1): " << D_P_Q1 << " nats\n";
86
    std::cout << "H(P,Q2): " << (std::isinf(HPQ2) ? std::string("inf") : std::to_string(HPQ2)) << ", KL(P||Q2): "
87
              << (std::isinf(D_P_Q2) ? std::string("inf") : std::to_string(D_P_Q2)) << "\n";
88
    std::cout << "KL(Q1||P): " << D_Q1_P << " nats (asymmetry check)\n";
89
90
    // Smoothed comparison to avoid infinities (for analysis/monitoring)
91
    double eps = 1e-12;
92
    auto P_s = smooth_and_normalize(P_raw, eps);
93
    auto Q2_s = smooth_and_normalize(Q2_raw, eps);
94
    double D_smoothed = kl_divergence(P_s, Q2_s);
95
    std::cout << "Smoothed KL(P||Q2) with eps=" << eps << ": " << D_smoothed << " nats\n";
96
97
    return 0;
98
}
99

1
#include <iostream>
2
#include <vector>
3
#include <cmath>
4
#include <stdexcept>
5
#include <iomanip>
6
7
// KL between N(mu0, diag(var0)) and N(mu1, diag(var1)) in d dimensions
8
// Formula: 0.5 * sum_i [ log(var1_i/var0_i) + (var0_i + (mu0_i - mu1_i)^2)/var1_i - 1 ]
9
double kl_diag_gaussians(const std::vector<double>& mu0,
10
                         const std::vector<double>& var0,
11
                         const std::vector<double>& mu1,
12
                         const std::vector<double>& var1,
13
                         double min_var = 1e-12) {
14
    size_t d = mu0.size();
15
    if (mu1.size() != d || var0.size() != d || var1.size() != d)
16
        throw std::invalid_argument("All vectors must have the same dimension");
17
    double D = 0.0;
18
    for (size_t i = 0; i < d; ++i) {
19
        double v0 = std::max(var0[i], min_var);
20
        double v1 = std::max(var1[i], min_var);
21
        double diff = mu0[i] - mu1[i];
22
        D += std::log(v1 / v0) + (v0 + diff * diff) / v1 - 1.0;
23
    }
24
    return 0.5 * D;
25
}
26
27
int main() {
28
    std::cout << std::fixed << std::setprecision(6);
29
30
    // 2D example
31
    std::vector<double> mu_p{0.0, 1.0};
32
    std::vector<double> var_p{1.0, 0.5};
33
    std::vector<double> mu_q{0.5, 0.5};
34
    std::vector<double> var_q{2.0, 0.25};
35
36
    double D_pq = kl_diag_gaussians(mu_p, var_p, mu_q, var_q);
37
    double D_qp = kl_diag_gaussians(mu_q, var_q, mu_p, var_p);
38
39
    std::cout << "KL(N_p || N_q) = " << D_pq << " nats\n";
40
    std::cout << "KL(N_q || N_p) = " << D_qp << " nats (asymmetry)\n";
41
42
    return 0;
43
}
44

1
#include <iostream>
2
#include <random>
3
#include <cmath>
4
#include <functional>
5
#include <vector>
6
#include <limits>
7
#include <iomanip>
8
9
const double PI_CONST = 3.14159265358979323846;
10
11
// Log-density of univariate Normal N(mu, sigma^2)
12
double log_norm_pdf(double x, double mu, double sigma) {
13
    double z = (x - mu) / sigma;
14
    return -0.5 * std::log(2.0 * PI_CONST) - std::log(sigma) - 0.5 * z * z;
15
}
16
17
// Stable log-sum-exp for two values
18
inline double logsumexp2(double a, double b) {
19
    double m = std::max(a, b);
20
    return m + std::log(std::exp(a - m) + std::exp(b - m));
21
}
22
23
// Define P as a 2-component mixture: 0.5*N(-2, 0.5^2) + 0.5*N(2, 1^2)
24
double log_p(double x) {
25
    double l1 = std::log(0.5) + log_norm_pdf(x, -2.0, 0.5);
26
    double l2 = std::log(0.5) + log_norm_pdf(x,  2.0, 1.0);
27
    return logsumexp2(l1, l2);
28
}
29
30
// Define Q as a Normal N(m, s^2)
31
struct QNormal { double m; double s; };
32
33
double log_q(double x, const QNormal& q) {
34
    return log_norm_pdf(x, q.m, q.s);
35
}
36
37
// Samplers
38
struct RNG { std::mt19937_64 gen; RNG(uint64_t seed=42) : gen(seed) {} };
39
40
// Sample from P mixture
41
double sample_P(RNG& rng) {
42
    std::uniform_real_distribution<double> U(0.0, 1.0);
43
    std::bernoulli_distribution B(0.5);
44
    bool pick_first = B(rng.gen);
45
    if (pick_first) {
46
        std::normal_distribution<double> N(-2.0, 0.5);
47
        return N(rng.gen);
48
    } else {
49
        std::normal_distribution<double> N(2.0, 1.0);
50
        return N(rng.gen);
51
    }
52
}
53
54
// Sample from Q
55
double sample_Q(RNG& rng, const QNormal& q) {
56
    std::normal_distribution<double> N(q.m, q.s);
57
    return N(rng.gen);
58
}
59
60
// Monte Carlo estimate of forward KL D_KL(P||Q): E_P[log p(x) - log q(x)]
61
double mc_forward_KL(size_t N, const QNormal& q, RNG& rng) {
62
    double sum = 0.0;
63
    for (size_t i = 0; i < N; ++i) {
64
        double x = sample_P(rng);
65
        sum += (log_p(x) - log_q(x, q));
66
    }
67
    return sum / static_cast<double>(N);
68
}
69
70
// Monte Carlo estimate of reverse KL D_KL(Q||P): E_Q[log q(x) - log p(x)]
71
double mc_reverse_KL(size_t N, const QNormal& q, RNG& rng) {
72
    double sum = 0.0;
73
    for (size_t i = 0; i < N; ++i) {
74
        double x = sample_Q(rng, q);
75
        sum += (log_q(x, q) - log_p(x));
76
    }
77
    return sum / static_cast<double>(N);
78
}
79
80
int main() {
81
    std::cout << std::fixed << std::setprecision(6);
82
    RNG rng(123);
83
84
    // Baseline Q
85
    QNormal q{0.0, 1.0};
86
    size_t N = 200000; // samples for MC
87
88
    double D_fwd = mc_forward_KL(N, q, rng);
89
    double D_rev = mc_reverse_KL(N, q, rng);
90
91
    std::cout << "Forward KL D_KL(P||Q) with Q=N(0,1): " << D_fwd << " nats\n";
92
    std::cout << "Reverse KL D_KL(Q||P) with Q=N(0,1): " << D_rev << " nats\n";
93
94
    // Simple 1D grid search over mean m to illustrate mode-covering vs mode-seeking
95
    double best_m_fwd = 0.0, best_val_fwd = std::numeric_limits<double>::infinity();
96
    double best_m_rev = 0.0, best_val_rev = std::numeric_limits<double>::infinity();
97
98
    for (int i = -30; i <= 30; ++i) {
99
        double m = i * 0.2; // grid from -6 to 6
100
        QNormal qtest{m, 1.0}; // fix s=1 for simplicity
101
        RNG rng_local(1000 + i); // different seeds to reduce correlation
102
        double val_fwd = mc_forward_KL(40000, qtest, rng_local);
103
        RNG rng_local2(2000 + i);
104
        double val_rev = mc_reverse_KL(40000, qtest, rng_local2);
105
        if (val_fwd < best_val_fwd) { best_val_fwd = val_fwd; best_m_fwd = m; }
106
        if (val_rev < best_val_rev) { best_val_rev = val_rev; best_m_rev = m; }
107
    }
108
109
    std::cout << "Minimizer of forward KL over m (s=1): m* = " << best_m_fwd
110
              << ", KL = " << best_val_fwd << "\n";
111
    std::cout << "Minimizer of reverse KL over m (s=1): m* = " << best_m_rev
112
              << ", KL = " << best_val_rev << "\n";
113
114
    std::cout << "(Expected: reverse KL tends to sit on one mode near -2 or 2; forward KL tends to sit between modes.)\n";
115
    return 0;
116
}
117