Information Theory for Machine Learning

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

# Mutual Information: Equivalent Definitions
# ===========================================
 
# Definition 1: Reduction in uncertainty
# "How much does knowing Y reduce uncertainty about X?"
I(X; Y) = H(X) - H(X|Y)
 
# Definition 2: Symmetric form
# "How much does knowing either reduce uncertainty about the other?"
I(X; Y) = H(X) + H(Y) - H(X, Y)
 
# Definition 3: KL divergence from independence
# "How far is the joint from the product of marginals?"
I(X; Y) = D_KL(P(X,Y) || P(X)P(Y))
        = ∑_x ∑_y P(x,y) log[P(x,y) / (P(x)P(y))]
 
# Definition 4: Expected pointwise mutual information
# "Average surprise about joint vs independent co-occurrence"
I(X; Y) = E_{P(X,Y)}[log(P(X,Y) / (P(X)P(Y)))]
 
# All definitions are mathematically equivalent!
 
# Python Implementation:
import numpy as np
 
def entropy(probs):
    """Compute entropy of a distribution."""
    probs = np.array(probs).flatten()
    probs = probs[probs > 0]
    return -np.sum(probs * np.log2(probs))
 
def mutual_information(joint_probs):
    """
    Compute mutual information from joint probability matrix.
    
    Args:
        joint_probs: 2D array where joint_probs[i,j] = P(X=i, Y=j)
    
    Returns:
        Mutual information I(X; Y) in bits
    """
    joint = np.array(joint_probs)
    assert np.abs(joint.sum() - 1.0) < 1e-6, "Joint probs must sum to 1"
    
    # Marginal distributions
    p_x = joint.sum(axis=1)  # Sum over Y
    p_y = joint.sum(axis=0)  # Sum over X
    
    # Method 1: I(X;Y) = H(X) + H(Y) - H(X,Y)
    h_x = entropy(p_x)
    h_y = entropy(p_y)
    h_xy = entropy(joint)
    mi_method1 = h_x + h_y - h_xy
    
    # Method 2: I(X;Y) = Σ P(x,y) log[P(x,y)/(P(x)P(y))]
    mi_method2 = 0
    for i in range(joint.shape[0]):
        for j in range(joint.shape[1]):
            if joint[i,j] > 0:
                mi_method2 += joint[i,j] * np.log2(joint[i,j] / (p_x[i] * p_y[j]))
    
    return mi_method1, mi_method2
 
# Example: Weather and umbrella carrying
joint = np.array([
    [0.63, 0.07],  # Sunny: [no umbrella, umbrella]
    [0.03, 0.27]   # Rainy: [no umbrella, umbrella]
])
 
mi1, mi2 = mutual_information(joint)
print(f"Mutual Information (method 1): {mi1:.4f} bits")
print(f"Mutual Information (method 2): {mi2:.4f} bits")
 
# Verify: marginals
p_weather = joint.sum(axis=1)  # [0.7, 0.3]
p_umbrella = joint.sum(axis=0)  # [0.66, 0.34]
print(f"\nP(Weather): {p_weather}")
print(f"P(Umbrella): {p_umbrella}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

import numpy as np
 
def mi_from_joint(joint):
    """Compute MI from joint probability matrix."""
    joint = np.array(joint)
    p_x = joint.sum(axis=1)
    p_y = joint.sum(axis=0)
    
    mi = 0
    for i in range(joint.shape[0]):
        for j in range(joint.shape[1]):
            if joint[i,j] > 0:
                mi += joint[i,j] * np.log2(joint[i,j] / (p_x[i] * p_y[j] + 1e-15))
    return mi
 
def entropy(probs):
    probs = np.array(probs)
    probs = probs[probs > 0]
    return -np.sum(probs * np.log2(probs))
 
# Property 1: Symmetry
print("Property: Symmetry")
joint = np.array([[0.4, 0.1], [0.1, 0.4]])
joint_transposed = joint.T
print(f"I(X; Y) = {mi_from_joint(joint):.4f}")
print(f"I(Y; X) = {mi_from_joint(joint_transposed):.4f}")
print()
 
# Property 2: Non-negativity
print("Property: Non-negativity (random tests)")
for _ in range(3):
    # Random joint distribution
    joint = np.random.dirichlet(np.ones(9)).reshape(3, 3)
    mi = mi_from_joint(joint)
    print(f"  I(X; Y) = {mi:.4f} >= 0: {mi >= -1e-10}")
print()
 
# Property 3: Zero for independence
print("Property: Zero for independent variables")
p_x = np.array([0.3, 0.4, 0.3])
p_y = np.array([0.5, 0.5])
joint_independent = np.outer(p_x, p_y)  # Product of marginals
print(f"Independent joint: I(X; Y) = {mi_from_joint(joint_independent):.6f} ≈ 0")
print()
 
# Property 4: Self-information
print("Property: I(X; X) = H(X)")
p = np.array([0.25, 0.25, 0.25, 0.25])
# Diagonal joint: P(X=x, X=x) = P(X=x), off-diagonal = 0
joint_self = np.diag(p)
i_xx = mi_from_joint(joint_self)
h_x = entropy(p)
print(f"I(X; X) = {i_xx:.4f}")
print(f"H(X) = {h_x:.4f}")
print()
 
# Property 5: Upper bound
print("Property: Upper bound")
joint = np.array([[0.25, 0.25], [0.25, 0.25]])  # Very dependent
p_x = joint.sum(axis=1)
p_y = joint.sum(axis=0)
mi = mi_from_joint(joint)
print(f"I(X; Y) = {mi:.4f}")
print(f"min(H(X), H(Y)) = min({entropy(p_x):.4f}, {entropy(p_y):.4f}) = {min(entropy(p_x), entropy(p_y)):.4f}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

import numpy as np
from sklearn.feature_selection import mutual_info_regression
 
np.random.seed(42)
n = 5000
 
print("Comparing Correlation vs Mutual Information")
print("=" * 60)
print()
 
# Case 1: Linear relationship (both work)
x1 = np.random.normal(0, 1, n)
y1 = 2 * x1 + np.random.normal(0, 0.5, n)
 
corr1 = np.corrcoef(x1, y1)[0, 1]
mi1 = mutual_info_regression(x1.reshape(-1, 1), y1)[0]
 
print("Case 1: Linear relationship (Y = 2X + noise)")
print(f"  Correlation: {corr1:.4f}")
print(f"  Mutual Information: {mi1:.4f} nats")
print()
 
# Case 2: Quadratic relationship (correlation fails!)
x2 = np.random.uniform(-3, 3, n)
y2 = x2**2 + np.random.normal(0, 0.5, n)
 
corr2 = np.corrcoef(x2, y2)[0, 1]
mi2 = mutual_info_regression(x2.reshape(-1, 1), y2)[0]
 
print("Case 2: Quadratic relationship (Y = X² + noise)")
print(f"  Correlation: {corr2:.4f} (near zero!)")
print(f"  Mutual Information: {mi2:.4f} nats (high!)")
print()
 
# Case 3: Sinusoidal relationship
x3 = np.random.uniform(0, 4*np.pi, n)
y3 = np.sin(x3) + np.random.normal(0, 0.2, n)
 
corr3 = np.corrcoef(x3, y3)[0, 1]
mi3 = mutual_info_regression(x3.reshape(-1, 1), y3)[0]
 
print("Case 3: Sinusoidal relationship (Y = sin(X) + noise)")
print(f"  Correlation: {corr3:.4f} (near zero!)")
print(f"  Mutual Information: {mi3:.4f} nats (high!)")
print()
 
# Case 4: XOR-like relationship (binary)
x4 = np.random.randint(0, 2, n)
y4 = np.random.randint(0, 2, n)
z4 = (x4 + y4) % 2  # XOR
 
corr4 = np.corrcoef(x4, z4)[0, 1]
# For discrete, use discrete MI
from sklearn.metrics import mutual_info_score
mi4 = mutual_info_score(x4, z4)
 
print("Case 4: XOR relationship (Z = X XOR Y)")
print(f"  Correlation(X, Z): {corr4:.4f}")
print(f"  Mutual Information: {mi4:.4f} nats")
print()
 
# Case 5: Independence (both should be zero)
x5 = np.random.normal(0, 1, n)
y5 = np.random.normal(0, 1, n)
 
corr5 = np.corrcoef(x5, y5)[0, 1]
mi5 = mutual_info_regression(x5.reshape(-1, 1), y5)[0]
 
print("Case 5: Independence (X and Y unrelated)")
print(f"  Correlation: {corr5:.4f}")
print(f"  Mutual Information: {mi5:.4f} nats (≈0)")
 
print()
print("Key insight: Correlation can be zero even with strong dependence!")
print("MI captures ALL dependencies, linear and nonlinear.")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

import numpy as np
from sklearn.datasets import make_classification
from sklearn.feature_selection import mutual_info_classif, SelectKBest
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
 
# Create dataset with known structure
np.random.seed(42)
n_samples = 1000
n_informative = 5
n_redundant = 3
n_useless = 12
 
# Generate data
X, y = make_classification(
    n_samples=n_samples,
    n_features=n_informative + n_redundant + n_useless,
    n_informative=n_informative,
    n_redundant=n_redundant,
    n_clusters_per_class=2,
    flip_y=0.05,
    random_state=42
)
 
feature_names = [f"F{i}" for i in range(X.shape[1])]
 
# Compute MI for each feature
mi_scores = mutual_info_classif(X, y, random_state=42)
 
# Rank features
ranking = np.argsort(mi_scores)[::-1]
 
print("Feature Ranking by Mutual Information")
print("=" * 50)
print(f"{'Rank':<6} {'Feature':<10} {'MI Score':<12} {'Type'}")
print("-" * 50)
 
for i, idx in enumerate(ranking):
    if idx < n_informative:
        ftype = "Informative"
    elif idx < n_informative + n_redundant:
        ftype = "Redundant"
    else:
        ftype = "Useless"
    print(f"{i+1:<6} {feature_names[idx]:<10} {mi_scores[idx]:<12.4f} {ftype}")
 
print()
 
# Select top features and compare performance
print("Model Performance Comparison")
print("-" * 50)
 
for k in [5, 10, 20]:
    selector = SelectKBest(mutual_info_classif, k=k)
    X_selected = selector.fit_transform(X, y)
    
    clf = RandomForestClassifier(n_estimators=50, random_state=42)
    scores = cross_val_score(clf, X_selected, y, cv=5)
    
    print(f"Top {k:2d} features: Accuracy = {scores.mean():.4f} ± {scores.std():.4f}")
 
# Full features
clf = RandomForestClassifier(n_estimators=50, random_state=42)
scores_full = cross_val_score(clf, X, y, cv=5)
print(f"All {X.shape[1]:2d} features: Accuracy = {scores_full.mean():.4f} ± {scores_full.std():.4f}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

import numpy as np
from sklearn.feature_selection import mutual_info_classif
from sklearn.metrics import mutual_info_score
 
def mrmr_feature_selection(X, y, n_features=5):
    """
    Minimum Redundancy Maximum Relevance feature selection.
    
    Score = I(f; y) - (1/|S|) * Σ I(f; s) for s in S
    
    Where S is the set of already selected features.
    """
    n_total = X.shape[1]
    
    # Compute relevance: I(feature; target)
    relevance = mutual_info_classif(X, y, random_state=42)
    
    # For redundancy, we'll discretize continuous features
    def discretize(x, n_bins=10):
        return np.digitize(x, np.percentile(x, np.linspace(0, 100, n_bins)))
    
    X_discrete = np.apply_along_axis(discretize, 0, X)
    
    selected = []
    remaining = list(range(n_total))
    
    for i in range(n_features):
        if i == 0:
            # First feature: pure relevance
            scores = relevance[remaining]
        else:
            # Subsequent: relevance - avg(redundancy with selected)
            scores = []
            for f in remaining:
                rel = relevance[f]
                # Average MI with already selected features
                redundancy = np.mean([
                    mutual_info_score(X_discrete[:, f], X_discrete[:, s])
                    for s in selected
                ])
                scores.append(rel - redundancy)
            scores = np.array(scores)
        
        # Select feature with best score
        best_idx = np.argmax(scores)
        best_feature = remaining[best_idx]
        
        selected.append(best_feature)
        remaining.remove(best_feature)
        
        print(f"Step {i+1}: Selected F{best_feature} "
              f"(relevance={relevance[best_feature]:.4f})")
    
    return selected
 
# Apply mRMR
print("mRMR Feature Selection")
print("=" * 50)
selected_features = mrmr_feature_selection(X, y, n_features=5)
print(f"\nSelected features: {['F'+str(i) for i in selected_features]}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

import numpy as np
from sklearn.feature_selection import mutual_info_regression
from scipy.stats import entropy
 
def mi_histogram(x, y, bins=10):
    """
    Histogram-based MI estimation.
    Simple but sensitive to binning.
    """
    # Create 2D histogram
    hist_2d, x_edges, y_edges = np.histogram2d(x, y, bins=bins)
    
    # Normalize to probabilities
    pxy = hist_2d / hist_2d.sum()
    px = pxy.sum(axis=1)
    py = pxy.sum(axis=0)
    
    # Compute MI
    mi = 0
    for i in range(len(px)):
        for j in range(len(py)):
            if pxy[i, j] > 0 and px[i] > 0 and py[j] > 0:
                mi += pxy[i, j] * np.log(pxy[i, j] / (px[i] * py[j]))
    
    return mi  # in nats
 
def mi_ksg(x, y, k=3):
    """
    KSG estimator (using sklearn implementation).
    More robust than histogram methods.
    """
    return mutual_info_regression(x.reshape(-1, 1), y, n_neighbors=k)[0]
 
# Generate test data with known relationship
np.random.seed(42)
n = 1000
 
# Linear relationship
x_linear = np.random.normal(0, 1, n)
y_linear = x_linear + np.random.normal(0, 0.5, n)
 
# Theoretical MI for jointly Gaussian: I(X;Y) = -0.5 * log(1 - ρ²)
rho = np.corrcoef(x_linear, y_linear)[0, 1]
mi_theoretical = -0.5 * np.log(1 - rho**2)
 
print("MI Estimation for Linear Gaussian (Y = X + noise)")
print("=" * 60)
print(f"Theoretical MI (Gaussian formula): {mi_theoretical:.4f} nats")
print(f"Correlation: {rho:.4f}")
print()
 
# Compare estimators
print(f"{'Method':<30} {'Estimate':<12} {'Error':<12}")
print("-" * 60)
 
# Histogram with different bin sizes
for bins in [5, 10, 20, 50]:
    est = mi_histogram(x_linear, y_linear, bins=bins)
    err = est - mi_theoretical
    print(f"Histogram (bins={bins}){'':<15} {est:<12.4f} {err:+.4f}")
 
# KSG with different k
for k in [3, 5, 10]:
    est = mi_ksg(x_linear, y_linear, k=k)
    err = est - mi_theoretical
    print(f"KSG (k={k}){'':<20} {est:<12.4f} {err:+.4f}")
 
print()
print("Note: KSG is generally more accurate and robust than histogram methods.")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

import torch
import torch.nn.functional as F
 
def info_nce_loss(query, positive_key, negative_keys, temperature=0.07):
    """
    InfoNCE loss for contrastive learning.
    
    This is a lower bound on I(query; key):
    I(Q; K) >= log(N) - L_NCE
    
    Args:
        query: Query representations (batch_size, dim)
        positive_key: Positive key for each query (batch_size, dim)
        negative_keys: Negative keys (num_negatives, dim) or (batch_size, num_neg, dim)
        temperature: Softmax temperature (lower = sharper)
    
    Returns:
        InfoNCE loss (lower is better for training)
    """
    batch_size = query.size(0)
    
    # Normalize representations
    query = F.normalize(query, dim=1)
    positive_key = F.normalize(positive_key, dim=1)
    
    # Positive logits: q · k+ 
    positive_logits = torch.sum(query * positive_key, dim=1, keepdim=True)
    positive_logits = positive_logits / temperature
    
    # Negative logits
    if negative_keys.dim() == 2:
        # Shared negatives across batch
        negative_keys = F.normalize(negative_keys, dim=1)
        negative_logits = query @ negative_keys.T / temperature
    else:
        # Per-sample negatives
        negative_keys = F.normalize(negative_keys, dim=2)
        negative_logits = torch.bmm(query.unsqueeze(1), 
                                     negative_keys.transpose(1, 2)).squeeze(1)
        negative_logits = negative_logits / temperature
    
    # Concatenate: [positive, negatives]
    logits = torch.cat([positive_logits, negative_logits], dim=1)
    
    # Labels: positive is always at index 0
    labels = torch.zeros(batch_size, dtype=torch.long, device=query.device)
    
    # Cross-entropy loss
    loss = F.cross_entropy(logits, labels)
    
    return loss
 
# Example usage
batch_size = 32
dim = 128
num_negatives = 256
 
# Simulated representations
query = torch.randn(batch_size, dim)
positive_key = query + torch.randn(batch_size, dim) * 0.1  # Similar
negative_keys = torch.randn(num_negatives, dim)  # Random
 
loss = info_nce_loss(query, positive_key, negative_keys)
print(f"InfoNCE Loss: {loss.item():.4f}")
 
# Lower bound on MI
mi_lower_bound = np.log(num_negatives + 1) - loss.item()
print(f"MI Lower Bound: {mi_lower_bound:.4f} nats")
 
# With more negatives
for n_neg in [16, 64, 256, 1024]:
    neg_keys = torch.randn(n_neg, dim)
    loss = info_nce_loss(query, positive_key, neg_keys)
    mi_bound = np.log(n_neg + 1) - loss.item()
    print(f"N={n_neg:4d}: Loss={loss.item():.4f}, MI bound={mi_bound:.4f}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

# Conditional Mutual Information
# ================================
 
# Definition:
I(X; Y | Z) = H(X | Z) - H(X | Y, Z)
            = H(X, Z) + H(Y, Z) - H(Z) - H(X, Y, Z)
            = E_Z[I(X; Y | Z=z)]
 
# Properties:
# 1. Chain rule: I(X; Y, Z) = I(X; Z) + I(X; Y | Z)
# 2. Non-negative: I(X; Y | Z) >= 0
# 3. Reduces to MI when Z is constant: I(X; Y | ∅) = I(X; Y)
 
# Key insight: I(X; Y | Z) can be LESS than I(X; Y)
# This happens when Z "explains" some of the dependence between X and Y
 
# Example: Confounding
# X = Ice cream sales
# Y = Drowning deaths
# Z = Temperature
# I(X; Y) > 0 (correlated!)
# I(X; Y | Z) ≈ 0 (independent given temperature)
 
import numpy as np
 
def conditional_mi_discrete(x, y, z):
    """
    Compute I(X; Y | Z) from discrete samples.
    """
    from collections import defaultdict
    
    n = len(x)
    
    # Count joint occurrences
    xyz_counts = defaultdict(int)
    xz_counts = defaultdict(int)
    yz_counts = defaultdict(int)
    z_counts = defaultdict(int)
    
    for i in range(n):
        xyz_counts[(x[i], y[i], z[i])] += 1
        xz_counts[(x[i], z[i])] += 1
        yz_counts[(y[i], z[i])] += 1
        z_counts[z[i]] += 1
    
    # Convert to probabilities and compute CMI
    cmi = 0
    for (xi, yi, zi), count in xyz_counts.items():
        p_xyz = count / n
        p_xz = xz_counts[(xi, zi)] / n
        p_yz = yz_counts[(yi, zi)] / n
        p_z = z_counts[zi] / n
        
        if p_xyz > 0 and p_xz > 0 and p_yz > 0 and p_z > 0:
            # I(X;Y|Z) = Σ p(x,y,z) log[p(x,y,z)p(z) / (p(x,z)p(y,z))]
            cmi += p_xyz * np.log2((p_xyz * p_z) / (p_xz * p_yz))
    
    return cmi
 
# Example: Confounding
np.random.seed(42)
n = 5000
 
# Z causes both X and Y
z = np.random.randint(0, 3, n)  # Low, Medium, High temperature
 
# X depends on Z (ice cream sales)
x = z + np.random.binomial(1, 0.3, n)
 
# Y depends on Z (drowning - more swimming in hot weather)  
y = z + np.random.binomial(1, 0.2, n)
 
# Compute MIs
from sklearn.metrics import mutual_info_score
 
mi_xy = mutual_info_score(x, y)
cmi_xy_z = conditional_mi_discrete(x, y, z)
 
print("Confounding Example: Ice Cream Sales (X) vs Drownings (Y)")
print("=" * 60)
print(f"I(X; Y) = {mi_xy:.4f} nats  <- Appears correlated!")
print(f"I(X; Y | Z=Temperature) = {cmi_xy_z:.4f} nats  <- Much less after conditioning!")
print()
print("Conditioning on temperature 'explains away' the spurious correlation.")