KL Divergence

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Normalize a vector to sum to 1 (if the sum is positive). Negative entries are clamped to 0.
5
static void normalize(vector<double>& v) {
6
    double sum = 0.0;
7
    for (double& x : v) {
8
        if (x < 0.0) x = 0.0; // clamp negatives
9
        sum += x;
10
    }
11
    if (sum > 0.0) {
12
        for (double& x : v) x /= sum;
13
    }
14
}
15
16
// Compute KL(P || Q) for discrete distributions represented as probabilities.
17
// - P, Q: vectors of probabilities (need not be normalized; we normalize internally)
18
// - log_base: base of logarithm (e.g., exp(1) for nats, 2.0 for bits)
19
// - smoothing: additive smoothing applied to each component before normalization (>= 0)
20
// Returns +infinity if support mismatch remains (p_i > 0 and q_i == 0) after smoothing.
21
double kl_divergence(vector<double> P, vector<double> Q, double log_base = std::exp(1.0), double smoothing = 0.0) {
22
    if (P.size() != Q.size() || P.empty()) {
23
        throw invalid_argument("P and Q must have the same non-zero length");
24
    }
25
    // Add optional smoothing to avoid zeros (Laplace/additive smoothing)
26
    if (smoothing < 0.0) smoothing = 0.0;
27
    for (size_t i = 0; i < P.size(); ++i) {
28
        P[i] = max(0.0, P[i]) + smoothing;
29
        Q[i] = max(0.0, Q[i]) + smoothing;
30
    }
31
    // Normalize to make them valid distributions
32
    normalize(P);
33
    normalize(Q);
34
35
    // Convert natural log to chosen base by dividing by ln(base)
36
    const double inv_log_base = 1.0 / log(log_base);
37
38
    // Compute KL sum, skipping terms where p_i == 0
39
    long double kl = 0.0L;
40
    for (size_t i = 0; i < P.size(); ++i) {
41
        const double p = P[i];
42
        const double q = Q[i];
43
        if (p == 0.0) continue; // term is zero
44
        if (q == 0.0) {
45
            return numeric_limits<double>::infinity(); // support mismatch after smoothing
46
        }
47
        kl += (long double)p * (log(p) - log(q)) * inv_log_base;
48
    }
49
    return (double)kl;
50
}
51
52
int main() {
53
    // Example 1: Simple 2-category distributions in nats
54
    vector<double> P = {0.4, 0.6};
55
    vector<double> Q = {0.5, 0.5};
56
    cout << fixed << setprecision(6);
57
    cout << "KL(P || Q) (nats) = " << kl_divergence(P, Q) << "\n";
58
59
    // Example 2: Handling zeros via smoothing; compute in bits
60
    vector<double> P2 = {0.0, 1.0, 0.0, 0.0};
61
    vector<double> Q2 = {0.25, 0.25, 0.25, 0.25};
62
    double kl_bits = kl_divergence(P2, Q2, /*base=*/2.0, /*smoothing=*/1e-12);
63
    cout << "KL(P2 || Q2) (bits) = " << kl_bits << "\n";
64
65
    // Example 3: If Q assigns zero probability where P has mass (after smoothing), KL is infinite
66
    vector<double> P3 = {0.5, 0.5};
67
    vector<double> Q3 = {1.0, 0.0}; // Q3[1] = 0 while P3[1] > 0
68
    double kl_inf = kl_divergence(P3, Q3, /*base=*/exp(1.0), /*smoothing=*/0.0);
69
    cout << "KL(P3 || Q3) (nats) = " << (isinf(kl_inf) ? string("infinity") : to_string(kl_inf)) << "\n";
70
71
    return 0;
72
}
73

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Build frequency map from integer samples
5
unordered_map<int, long long> counts(const vector<int>& samples) {
6
    unordered_map<int, long long> c;
7
    for (int x : samples) ++c[x];
8
    return c;
9
}
10
11
// Compute KL from empirical samples using additive (Laplace) smoothing with parameter alpha >= 0
12
// P_samples ~ P, Q_samples ~ Q. The union of observed categories defines bins.
13
double kl_from_samples(const vector<int>& P_samples, const vector<int>& Q_samples, double alpha = 1.0, double log_base = std::exp(1.0)) {
14
    if (P_samples.empty() || Q_samples.empty()) {
15
        throw invalid_argument("Sample sets must be non-empty");
16
    }
17
    auto cP = counts(P_samples);
18
    auto cQ = counts(Q_samples);
19
20
    // Union of keys determines the support for smoothing
21
    vector<int> keys;
22
    keys.reserve(cP.size() + cQ.size());
23
    for (auto& kv : cP) keys.push_back(kv.first);
24
    for (auto& kv : cQ) keys.push_back(kv.first);
25
    sort(keys.begin(), keys.end());
26
    keys.erase(unique(keys.begin(), keys.end()), keys.end());
27
    const size_t K = keys.size();
28
29
    const long long nP = (long long)P_samples.size();
30
    const long long nQ = (long long)Q_samples.size();
31
32
    const double denomP = nP + alpha * (double)K;
33
    const double denomQ = nQ + alpha * (double)K;
34
35
    const double inv_log_base = 1.0 / log(log_base);
36
37
    long double kl = 0.0L;
38
    for (int key : keys) {
39
        const double p = ((double)cP[key] + alpha) / denomP;
40
        const double q = ((double)cQ[key] + alpha) / denomQ;
41
        // With alpha > 0, q and p are strictly positive, avoiding infinite KL
42
        kl += (long double)p * (log(p) - log(q)) * inv_log_base;
43
    }
44
    return (double)kl;
45
}
46
47
int main() {
48
    // Simulate two categorical sources over integers {0,1,2}
49
    vector<int> P_samples = {0,0,0,1,1,2,0,1,0,2,0,1,1,0,2}; // P biased toward 0 and 1
50
    vector<int> Q_samples = {0,1,2,2,2,2,1,2,1,2,2,1,2,2,2}; // Q biased toward 2
51
52
    cout << fixed << setprecision(6);
53
    double kl_nats = kl_from_samples(P_samples, Q_samples, /*alpha=*/1.0, /*log_base=*/exp(1.0));
54
    cout << "KL(P || Q) from samples (nats, alpha=1) = " << kl_nats << "\n";
55
56
    double kl_bits = kl_from_samples(P_samples, Q_samples, /*alpha=*/1.0, /*log_base=*/2.0);
57
    cout << "KL(P || Q) from samples (bits, alpha=1) = " << kl_bits << "\n";
58
59
    return 0;
60
}
61

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Log-PDF of univariate normal N(mu, sigma^2)
5
double log_normal_pdf(double x, double mu, double sigma) {
6
    const double var = sigma * sigma;
7
    const double log_norm = -0.5 * log(2.0 * M_PI * var);
8
    const double quad = -0.5 * (x - mu) * (x - mu) / var;
9
    return log_norm + quad; // natural log (nats)
10
}
11
12
// Closed-form KL for univariate Gaussians in nats
13
// KL( N(mu0, sigma0^2) || N(mu1, sigma1^2) )
14
double kl_gaussian_closed_form(double mu0, double sigma0, double mu1, double sigma1) {
15
    if (sigma0 <= 0.0 || sigma1 <= 0.0) {
16
        throw invalid_argument("Standard deviations must be positive");
17
    }
18
    double term1 = log(sigma1 / sigma0);
19
    double term2 = (sigma0 * sigma0 + (mu0 - mu1) * (mu0 - mu1)) / (2.0 * sigma1 * sigma1);
20
    return term1 + term2 - 0.5;
21
}
22
23
// Monte Carlo estimator of KL by sampling x ~ N(mu0, sigma0^2)
24
// Returns estimate in nats using natural logs
25
double kl_gaussian_monte_carlo(double mu0, double sigma0, double mu1, double sigma1, size_t N, unsigned seed = 42) {
26
    mt19937 rng(seed);
27
    normal_distribution<double> distP(mu0, sigma0);
28
    long double acc = 0.0L;
29
    for (size_t i = 0; i < N; ++i) {
30
        double x = distP(rng);
31
        acc += (long double)(log_normal_pdf(x, mu0, sigma0) - log_normal_pdf(x, mu1, sigma1));
32
    }
33
    return (double)(acc / (long double)N);
34
}
35
36
int main() {
37
    cout << fixed << setprecision(6);
38
39
    double mu0 = 0.0, sigma0 = 1.0;
40
    double mu1 = 1.0, sigma1 = 2.0;
41
42
    double kl_cf = kl_gaussian_closed_form(mu0, sigma0, mu1, sigma1);
43
    cout << "Closed-form KL (nats) = " << kl_cf << "\n";
44
    cout << "Closed-form KL (bits) = " << kl_cf / log(2.0) << "\n";
45
46
    // Monte Carlo check (increase N for higher accuracy)
47
    size_t N = 100000; // 1e5 samples
48
    double kl_mc = kl_gaussian_monte_carlo(mu0, sigma0, mu1, sigma1, N, 123);
49
    cout << "Monte Carlo KL (nats, N=100000) = " << kl_mc << "\n";
50
51
    return 0;
52
}
53