Rényi Entropy & Divergence

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Compute log(sum_i exp(x_i)) stably.
5
static double logSumExp(const vector<double>& xs) {
6
    if (xs.empty()) return -numeric_limits<double>::infinity();
7
    double m = -numeric_limits<double>::infinity();
8
    for (double x : xs) m = max(m, x);
9
    if (!isfinite(m)) return m; // all -inf
10
    long double s = 0.0L;
11
    for (double x : xs) s += expl((long double)x - (long double)m);
12
    return m + log((double)s);
13
}
14
15
// Shannon entropy in arbitrary log base (default base 2 -> bits)
16
static double shannonEntropy(const vector<double>& p, double logBase = 2.0) {
17
    long double acc = 0.0L;
18
    for (double pi : p) {
19
        if (pi > 0.0) acc -= (long double)pi * log((long double)pi);
20
    }
21
    return (double)(acc / log(logBase));
22
}
23
24
// Rényi entropy H_alpha with robust handling of special cases.
25
// Assumes p is a valid PMF (nonnegative and sums to ~1).
26
// alpha > 0. alpha != 1. Use alpha = INFINITY for min-entropy.
27
static double renyiEntropy(const vector<double>& p, double alpha, double logBase = 2.0) {
28
    // Handle special cases explicitly
29
    if (alpha == 1.0) {
30
        return shannonEntropy(p, logBase);
31
    }
32
    if (alpha == 0.0) {
33
        // Hartley entropy: log of support size (count of pi > 0)
34
        size_t k = 0; for (double pi : p) if (pi > 0.0) ++k;
35
        return log((double)k) / log(logBase);
36
    }
37
    if (isinf(alpha)) {
38
        // Min-entropy: -log(max pi)
39
        double mp = 0.0; for (double pi : p) mp = max(mp, pi);
40
        if (mp <= 0.0) return numeric_limits<double>::infinity();
41
        return -log(mp) / log(logBase);
42
    }
43
    if (alpha <= 0.0) {
44
        throw invalid_argument("alpha must be > 0 for Rényi entropy");
45
    }
46
47
    // Compute log(sum_i p_i^alpha) via log-sum-exp on alpha * log p_i
48
    vector<double> terms; terms.reserve(p.size());
49
    for (double pi : p) {
50
        if (pi > 0.0) {
51
            terms.push_back(alpha * log(pi));
52
        }
53
        // If pi == 0, term is 0 in the sum; skip to avoid log(0)
54
    }
55
    double lsum = logSumExp(terms); // log of \sum p_i^alpha (in natural log)
56
57
    double H = (1.0 / (1.0 - alpha)) * lsum; // still in nats
58
    return H / log(logBase); // convert to chosen base
59
}
60
61
int main() {
62
    // Example distributions
63
    vector<double> uniform = {0.25, 0.25, 0.25, 0.25};
64
    vector<double> skewed  = {0.7, 0.1, 0.1, 0.1};
65
66
    // Demonstrate special cases and general α
67
    vector<double> alphas = {0.0, 0.5, 1.0, 2.0, numeric_limits<double>::infinity()};
68
69
    cout.setf(ios::fixed); cout << setprecision(6);
70
    cout << "Uniform distribution H_α (bits):\n";
71
    for (double a : alphas) {
72
        double H = renyiEntropy(uniform, a, 2.0);
73
        cout << "  alpha=" << (isinf(a)? string("inf"): to_string(a)) << ": " << H << "\n";
74
    }
75
76
    cout << "\nSkewed distribution H_α (bits):\n";
77
    for (double a : alphas) {
78
        double H = renyiEntropy(skewed, a, 2.0);
79
        cout << "  alpha=" << (isinf(a)? string("inf"): to_string(a)) << ": " << H << "\n";
80
    }
81
82
    // Sanity checks: monotonicity H_alpha decreasing in alpha
83
    // and uniform has maximum entropy for given support size.
84
    return 0;
85
}
86

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
static double logSumExpStreamed(const vector<double>& a) {
5
    // Streaming two-pass LSE without storing exponents: first pass finds max, second accumulates
6
    double m = -numeric_limits<double>::infinity();
7
    for (double x : a) m = max(m, x);
8
    if (!isfinite(m)) return m;
9
    long double s = 0.0L;
10
    for (double x : a) s += expl((long double)x - (long double)m);
11
    return m + log((double)s);
12
}
13
14
// Compute D_alpha(P||Q) for alpha>0, alpha!=1, with log base control (default bits).
15
// P and Q must be same length and sum to 1 (approximately).
16
static double renyiDivergence(const vector<double>& P, const vector<double>& Q, double alpha, double logBase = 2.0) {
17
    if (P.size() != Q.size()) throw invalid_argument("P and Q must have same size");
18
    if (alpha <= 0.0 || alpha == 1.0) throw invalid_argument("alpha must be >0 and !=1");
19
20
    // Domain check: for alpha>1, require absolute continuity of P w.r.t Q
21
    if (alpha > 1.0) {
22
        for (size_t i = 0; i < P.size(); ++i) {
23
            if (Q[i] == 0.0 && P[i] > 0.0) {
24
                return numeric_limits<double>::infinity();
25
            }
26
        }
27
    }
28
29
    vector<double> logTerms; logTerms.reserve(P.size());
30
    for (size_t i = 0; i < P.size(); ++i) {
31
        double p = P[i], q = Q[i];
32
        if (p == 0.0) {
33
            // term is zero regardless of q; skip to avoid -inf logs
34
            continue;
35
        }
36
        if (q == 0.0) {
37
            // If we reach here, p>0. For alpha<1, q^{1-alpha}=0 => term 0 (skip). For alpha>1 handled above as inf.
38
            if (alpha < 1.0) continue;
39
        } else {
40
            // x_i = alpha*log p + (1-alpha)*log q, contributes exp(x_i) to the inner sum
41
            logTerms.push_back(alpha * log(p) + (1.0 - alpha) * log(q));
42
        }
43
    }
44
45
    if (logTerms.empty()) {
46
        // Sum inside log is zero => divergence is -inf scaled? In valid probability settings this occurs
47
        // only if P puts all mass where Q=0 and alpha<1; D_alpha tends to +infinity as alpha->1-, but for alpha<1 the sum can be 0.
48
        // Mathematically, if the inner sum is 0, log(0) = -inf and after scaling 1/(alpha-1)<0, result is +inf.
49
        return numeric_limits<double>::infinity();
50
    }
51
52
    double lsum = logSumExpStreamed(logTerms); // natural log
53
    double D = (1.0 / (alpha - 1.0)) * lsum; // in nats
54
    return D / log(logBase); // convert to chosen base
55
}
56
57
int main() {
58
    vector<double> P = {0.7, 0.2, 0.1};
59
    vector<double> Q = {0.5, 0.3, 0.2};
60
61
    cout.setf(ios::fixed); cout << setprecision(6);
62
63
    // Compare divergence at different alpha
64
    for (double a : {0.5, 0.9, 2.0}) {
65
        double D = renyiDivergence(P, Q, a, 2.0);
66
        cout << "D_" << a << "(P||Q) = " << D << " bits\n";
67
    }
68
69
    // Edge case: Q has zero where P has mass and alpha>1 => infinite divergence
70
    vector<double> P2 = {0.8, 0.2};
71
    vector<double> Q2 = {0.0, 1.0};
72
    cout << "D_2(P2||Q2) = " << renyiDivergence(P2, Q2, 2.0, 2.0) << " (inf expected)\n";
73
    return 0;
74
}
75

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Build empirical PMF from integer-labeled samples and compute H_alpha.
5
// This is a plug-in estimator: p_i = count_i / N.
6
static double renyiEntropyFromSamples(const vector<int>& samples, double alpha, double logBase = 2.0) {
7
    if (samples.empty()) return 0.0; // undefined in theory; return 0 pragmatically
8
9
    unordered_map<int, long long> freq;
10
    freq.reserve(samples.size()*2);
11
    for (int x : samples) ++freq[x];
12
13
    vector<double> p; p.reserve(freq.size());
14
    const long double N = (long double)samples.size();
15
    for (auto& kv : freq) p.push_back((double)(kv.second / N));
16
17
    // Reuse the robust implementation from Example 1 (inline here for brevity)
18
    auto shannon = [&](const vector<double>& probs) {
19
        long double acc = 0.0L;
20
        for (double pi : probs) if (pi > 0.0) acc -= (long double)pi * log((long double)pi);
21
        return (double)(acc / log(logBase));
22
    };
23
24
    if (alpha == 1.0) return shannon(p);
25
26
    if (alpha == 0.0) {
27
        size_t k = 0; for (double pi : p) if (pi > 0.0) ++k;
28
        return log((double)k) / log(logBase);
29
    }
30
31
    if (isinf(alpha)) {
32
        double mp = 0.0; for (double pi : p) mp = max(mp, pi);
33
        return -log(mp) / log(logBase);
34
    }
35
36
    if (alpha <= 0.0) throw invalid_argument("alpha must be > 0");
37
38
    // log-sum-exp over alpha*log p_i
39
    double m = -numeric_limits<double>::infinity();
40
    for (double pi : p) if (pi > 0.0) m = max(m, alpha * log(pi));
41
    long double s = 0.0L;
42
    for (double pi : p) if (pi > 0.0) s += expl((long double)(alpha * log(pi) - m));
43
    double lsum = m + log((double)s);
44
45
    double H = (1.0 / (1.0 - alpha)) * lsum; // nats
46
    return H / log(logBase);
47
}
48
49
int main() {
50
    // Simulate samples from a skewed 3-category distribution
51
    vector<int> samples;
52
    mt19937 rng(123);
53
    discrete_distribution<int> dist({0.7, 0.2, 0.1});
54
    const int N = 50000;
55
    samples.reserve(N);
56
    for (int i = 0; i < N; ++i) samples.push_back(dist(rng));
57
58
    cout.setf(ios::fixed); cout << setprecision(6);
59
    for (double a : {0.5, 1.0, 2.0}) {
60
        double Hhat = renyiEntropyFromSamples(samples, a, 2.0);
61
        cout << "Estimated H_" << a << " = " << Hhat << " bits\n";
62
    }
63
64
    // Note: This plug-in estimator is biased for finite N; bias can be noticeable for α>1.
65
    return 0;
66
}
67