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

import numpy as np
from scipy.special import gamma
import matplotlib.pyplot as plt
 
def ball_volume(r: float, d: int) -> float:
    """Volume of a d-dimensional ball of radius r."""
    return (np.pi ** (d/2) / gamma(d/2 + 1)) * r**d
 
def cube_volume(d: int) -> float:
    """Volume of d-dimensional unit cube [0,1]^d."""
    return 1.0  # Always 1
 
def fraction_in_ball(r: float, d: int) -> float:
    """Fraction of unit cube volume covered by ball of radius r at center."""
    # Ball might extend outside cube, so we approximate with direct calculation
    # For small r relative to cube, this is approximately ball_volume / cube_volume
    return min(1.0, ball_volume(r, d))
 
# How the "reach" of a local neighborhood grows
print("The Emptiness of High-Dimensional Space")
print("=" * 60)
print("\nTo capture 1% of volume (and ~1% of uniformly distributed data):")
print("-" * 60)
print(f"{'Dimension d':<15} | {'Required radius r':<20} | {'r as % of [0,1]'}")
print("-" * 60)
 
for d in [1, 2, 3, 5, 10, 20, 50, 100]:
    # Find r such that ball volume ≈ 0.01 * cube volume
    # V_d(r) = 0.01 → r = (0.01 * Γ(d/2+1) / π^(d/2))^(1/d)
    target_fraction = 0.01
    r = (target_fraction * gamma(d/2 + 1) / np.pi**(d/2)) ** (1/d)
    
    print(f"{d:<15} | {r:<20.4f} | {100*r:<.1f}%")
 
print("\n⚠ In 100D, a 'local' ball capturing 1% of data extends ~90% ")
print("  of the way along each axis! 'Local' becomes 'global'.")

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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99

import numpy as np
from typing import Callable, List, Tuple
 
def simulate_kernel_regression_error(
    d: int,
    n_values: List[int],
    n_trials: int = 20,
    noise_std: float = 0.3
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Simulate kernel regression in d dimensions and measure MSE.
    
    True function: f(x) = sum of sines (smooth, d-dimensional)
    """
    mse_means = []
    mse_stds = []
    
    def f_true(X):
        """Smooth additive function for testing."""
        return np.sum(np.sin(2 * X), axis=1)
    
    def kernel_regression_predict(X_train, y_train, X_test, h):
        """Simple Nadaraya-Watson in d dimensions."""
        n_test = len(X_test)
        y_pred = np.zeros(n_test)
        
        for i in range(n_test):
            # Euclidean distances
            diff = X_train - X_test[i]
            distances = np.sqrt(np.sum(diff**2, axis=1))
            
            # Gaussian kernel
            weights = np.exp(-distances**2 / (2 * h**2))
            weight_sum = np.sum(weights)
            
            if weight_sum > 1e-10:
                y_pred[i] = np.sum(weights * y_train) / weight_sum
            else:
                y_pred[i] = np.mean(y_train)
        
        return y_pred
    
    for n in n_values:
        trial_mses = []
        
        for trial in range(n_trials):
            # Generate data
            X_train = np.random.uniform(0, 2*np.pi, (n, d))
            y_train = f_true(X_train) + np.random.normal(0, noise_std, n)
            
            # Test points
            n_test = 100
            X_test = np.random.uniform(0, 2*np.pi, (n_test, d))
            y_test_true = f_true(X_test)
            
            # Use rule-of-thumb bandwidth scaled for dimension
            sigma_x = 2*np.pi / np.sqrt(12)  # std of uniform on [0, 2π]
            h = sigma_x * n**(-1/(4 + d))  # Optimal rate scaling
            
            # Predict
            y_pred = kernel_regression_predict(X_train, y_train, X_test, h)
            
            mse = np.mean((y_pred - y_test_true)**2)
            trial_mses.append(mse)
        
        mse_means.append(np.mean(trial_mses))
        mse_stds.append(np.std(trial_mses))
    
    return np.array(mse_means), np.array(mse_stds)
 
# =============================================================================
# Run simulation for different dimensions
# =============================================================================
np.random.seed(42)
 
dimensions = [1, 2, 5, 10]
n_values = [50, 100, 200, 500, 1000, 2000]
 
print("Kernel Regression MSE vs Sample Size and Dimension")
print("=" * 70)
print(f"{'n':<8}", end="")
for d in dimensions:
    print(f"{'d=' + str(d):<12}", end="")
print()
print("-" * 70)
 
results = {}
for d in dimensions:
    mse_means, mse_stds = simulate_kernel_regression_error(d, n_values, n_trials=10)
    results[d] = mse_means
 
for i, n in enumerate(n_values):
    print(f"{n:<8}", end="")
    for d in dimensions:
        print(f"{results[d][i]:<12.4f}", end="")
    print()
 
print("\n⚠ Note how MSE barely improves with n in high dimensions!")
print("  This is the curse of dimensionality in action.")

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

import numpy as np
 
def demonstrate_distance_concentration(dimensions: list, n_points: int = 1000):
    """
    Show how distances concentrate in high dimensions.
    """
    print("Distance Concentration in High Dimensions")
    print("=" * 70)
    print("Random points in unit hypercube [0,1]^d")
    print("-" * 70)
    print(f"{'Dim d':<8} | {'Mean NN dist':<15} | {'Mean all-pairs dist':<20} | {'Ratio':<10}")
    print("-" * 70)
    
    for d in dimensions:
        # Generate random points
        X = np.random.uniform(0, 1, (n_points, d))
        
        # Compute all pairwise distances
        all_distances = []
        nn_distances = []
        
        for i in range(min(n_points, 500)):  # Limit for speed
            dists = np.sqrt(np.sum((X - X[i])**2, axis=1))
            dists[i] = np.inf  # Exclude self
            
            nn_distances.append(np.min(dists))
            all_distances.extend(dists[dists < np.inf])
        
        mean_nn = np.mean(nn_distances)
        mean_all = np.mean(all_distances)
        ratio = mean_nn / mean_all
        
        print(f"{d:<8} | {mean_nn:<15.4f} | {mean_all:<20.4f} | {ratio:<10.4f}")
    
    print("\n⚠ As d increases, nearest neighbors become almost as far as average!")
    print("  'Nearest' loses meaning—all points are effectively far from each other.")
 
def demonstrate_volume_in_shell(dimensions: list):
    """
    Show how volume concentrates in thin shells in high dimensions.
    """
    print("\n\nVolume Concentration in Thin Shells")
    print("=" * 70)
    print("Fraction of unit ball volume in outer 10% shell (r ∈ [0.9, 1.0])")
    print("-" * 70)
    
    for d in dimensions:
        inner_frac = 0.9 ** d  # Volume within r=0.9 as fraction of r=1.0
        shell_frac = 1 - inner_frac  # Volume in shell
        
        print(f"d = {d:>3}: {100*shell_frac:>6.2f}% of volume in outer 10% shell")
    
    print("\n⚠ In high dimensions, almost ALL volume is near the surface!")
 
# Run demonstrations
np.random.seed(42)
demonstrate_distance_concentration([1, 2, 5, 10, 20, 50, 100])
demonstrate_volume_in_shell([1, 2, 5, 10, 20, 50, 100, 500])

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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110

import numpy as np
from typing import Tuple
 
def additive_kernel_regression(
    X: np.ndarray,
    y: np.ndarray,
    X_test: np.ndarray,
    h: float = None
) -> np.ndarray:
    """
    Additive model via backfitting with univariate kernel smoothers.
    
    f(x_1, ..., x_d) = μ + f_1(x_1) + ... + f_d(x_d)
    
    This ESCAPES the curse of dimensionality by solving d 
    one-dimensional problems instead of one d-dimensional problem.
    """
    n, d = X.shape
    n_test = X_test.shape[0]
    
    if h is None:
        h = 1.06 * np.std(X) * n**(-0.2)  # Rule of thumb for 1D
    
    def smooth_1d(x, y, x_new, h):
        """Nadaraya-Watson in 1D."""
        y_hat = np.zeros(len(x_new))
        for i, x0 in enumerate(x_new):
            u = (x - x0) / h
            w = np.exp(-u**2 / 2)  # Gaussian kernel
            if np.sum(w) > 1e-10:
                y_hat[i] = np.sum(w * y) / np.sum(w)
            else:
                y_hat[i] = np.mean(y)
        return y_hat
    
    # Initialize: overall mean
    mu = np.mean(y)
    f_components = [np.zeros(n) for _ in range(d)]
    
    # Backfitting iterations
    max_iter = 10
    for iteration in range(max_iter):
        for j in range(d):
            # Partial residual: y minus all other components
            partial_residual = y - mu - sum(f for k, f in enumerate(f_components) if k != j)
            
            # Smooth partial residual against x_j
            f_components[j] = smooth_1d(X[:, j], partial_residual, X[:, j], h)
            
            # Center the component
            f_components[j] -= np.mean(f_components[j])
    
    # Predict on test set
    y_test = np.full(n_test, mu)
    for j in range(d):
        # Train smoothed function maps X[:,j] -> f_j
        # We need to interpolate for test points
        partial_residual = y - mu - sum(f for k, f in enumerate(f_components) if k != j)
        y_test += smooth_1d(X[:, j], partial_residual, X_test[:, j], h)
    
    return y_test
 
# =============================================================================
# Demonstration: Additive vs Full Nonparametric in High Dimensions
# =============================================================================
np.random.seed(42)
 
d = 10  # 10 dimensions - challenging for full nonparametric
n_train = 500
n_test = 100
 
# True additive function (sum of sines)
f_true = lambda X: np.sum(np.sin(X), axis=1)
 
# Generate data
X_train = np.random.uniform(0, 2*np.pi, (n_train, d))
y_train = f_true(X_train) + np.random.normal(0, 0.3, n_train)
X_test = np.random.uniform(0, 2*np.pi, (n_test, d))
y_test_true = f_true(X_test)
 
# Additive model (exploits structure)
y_pred_additive = additive_kernel_regression(X_train, y_train, X_test, h=0.5)
 
# Full kernel regression (doesn't exploit structure)
def full_kernel_regression(X_train, y_train, X_test, h):
    y_pred = np.zeros(len(X_test))
    for i in range(len(X_test)):
        diff = X_train - X_test[i]
        dist = np.sqrt(np.sum(diff**2, axis=1))
        w = np.exp(-dist**2 / (2*h**2))
        if np.sum(w) > 1e-10:
            y_pred[i] = np.sum(w * y_train) / np.sum(w)
        else:
            y_pred[i] = np.mean(y_train)
    return y_pred
 
# Use a reasonable bandwidth for d=10
h_full = 2.0 * n_train**(-1/14)  # Optimal rate scaling for d=10
y_pred_full = full_kernel_regression(X_train, y_train, X_test, h_full)
 
mse_additive = np.mean((y_pred_additive - y_test_true)**2)
mse_full = np.mean((y_pred_full - y_test_true)**2)
 
print("Escaping the Curse via Additive Structure (d=10)")
print("=" * 60)
print(f"Additive Model MSE: {mse_additive:.4f}")
print(f"Full Kernel MSE:    {mse_full:.4f}")
print(f"Improvement ratio:  {mse_full / mse_additive:.1f}x")
print("\n⚠ When the truth IS additive, the additive model wins dramatically!")
print("  It solves 10 easy 1D problems instead of one impossible 10D problem.")