Kernel Design for Gaussian Processes

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

import numpy as np
from numpy.linalg import eigvalsh
 
def check_psd(K, tolerance=1e-8):
    """
    Check if a matrix is positive semi-definite.
    
    A matrix K is PSD if all eigenvalues are >= 0.
    In practice, we allow small negative eigenvalues due to numerical error.
    """
    eigenvalues = eigvalsh(K)  # Eigenvalues of symmetric matrix
    min_eig = np.min(eigenvalues)
    is_psd = min_eig >= -tolerance
    return is_psd, min_eig
 
def demonstrate_closure_properties():
    """
    Demonstrate that kernel operations preserve PSD property.
    """
    # Generate random test points
    np.random.seed(42)
    n = 50
    X = np.random.randn(n, 2)
    
    # Define two base kernels
    def rbf_kernel(X, length_scale=1.0):
        sq_dists = np.sum((X[:, np.newaxis, :] - X[np.newaxis, :, :])**2, axis=2)
        return np.exp(-sq_dists / (2 * length_scale**2))
    
    def linear_kernel(X):
        return X @ X.T
    
    K_rbf = rbf_kernel(X, length_scale=1.0)
    K_linear = linear_kernel(X)
    
    print("Base kernels:")
    is_psd, min_eig = check_psd(K_rbf)
    print(f"  RBF kernel: PSD={is_psd}, min eigenvalue={min_eig:.2e}")
    is_psd, min_eig = check_psd(K_linear)
    print(f"  Linear kernel: PSD={is_psd}, min eigenvalue={min_eig:.2e}")
    
    # Test addition
    K_sum = K_rbf + K_linear
    is_psd, min_eig = check_psd(K_sum)
    print(f"\nAddition (RBF + Linear): PSD={is_psd}, min eigenvalue={min_eig:.2e}")
    
    # Test multiplication (Hadamard product)
    K_prod = K_rbf * K_linear
    is_psd, min_eig = check_psd(K_prod)
    print(f"Multiplication (RBF × Linear): PSD={is_psd}, min eigenvalue={min_eig:.2e}")
    
    # Test scalar multiplication
    K_scaled = 5.0 * K_rbf
    is_psd, min_eig = check_psd(K_scaled)
    print(f"Scalar multiplication (5 × RBF): PSD={is_psd}, min eigenvalue={min_eig:.2e}")
    
    # Test exponentiation
    K_exp = np.exp(K_rbf)  # Note: exp of a valid kernel, not exp(-dist^2)
    is_psd, min_eig = check_psd(K_exp)
    print(f"Exponentiation (exp(RBF)): PSD={is_psd}, min eigenvalue={min_eig:.2e}")
    
    # WARNING: Test subtraction (NOT valid operation)
    K_diff = K_rbf - 0.5 * K_linear
    is_psd, min_eig = check_psd(K_diff)
    print(f"\n⚠️ Subtraction (RBF - 0.5×Linear): PSD={is_psd}, min eigenvalue={min_eig:.2e}")
    print("   Note: Subtraction does NOT preserve PSD property!")
    
demonstrate_closure_properties()

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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

import numpy as np
import matplotlib.pyplot as plt
 
def rbf_kernel(X1, X2, sigma_sq=1.0, length_scale=1.0):
    """Standard RBF (Squared Exponential) kernel."""
    X1, X2 = np.atleast_2d(X1).T, np.atleast_2d(X2).T
    sq_dists = np.sum((X1[:, np.newaxis] - X2.T[np.newaxis, :])**2, axis=2)
    return sigma_sq * np.exp(-sq_dists / (2 * length_scale**2))
 
def linear_kernel(X1, X2, sigma_sq=1.0, offset=0.0):
    """Linear (Bayesian linear regression) kernel."""
    X1 = np.atleast_1d(X1).flatten() - offset
    X2 = np.atleast_1d(X2).flatten() - offset
    return sigma_sq * np.outer(X1, X2)
 
def periodic_kernel(X1, X2, sigma_sq=1.0, period=1.0, length_scale=1.0):
    """Periodic kernel for cyclic patterns."""
    X1, X2 = np.atleast_1d(X1).flatten(), np.atleast_1d(X2).flatten()
    dists = np.abs(np.subtract.outer(X1, X2))
    return sigma_sq * np.exp(-2 * (np.sin(np.pi * dists / period) / length_scale)**2)
 
def white_noise_kernel(n, sigma_sq=0.1):
    """White noise (diagonal) kernel for observation noise."""
    return sigma_sq * np.eye(n)
 
def demonstrate_additive_decomposition():
    """
    Demonstrate how additive kernels decompose into interpretable components.
    """
    X = np.linspace(0, 10, 200)
    n = len(X)
    
    # Define component kernels
    K_trend = linear_kernel(X, X, sigma_sq=0.5, offset=5.0)
    K_smooth = rbf_kernel(X, X, sigma_sq=1.0, length_scale=2.0).reshape(n, n)
    K_periodic = periodic_kernel(X, X, sigma_sq=0.5, period=2.0, length_scale=0.5)
    K_noise = white_noise_kernel(n, sigma_sq=0.1)
    
    # Combined kernel
    K_total = K_trend + K_smooth + K_periodic + K_noise
    K_total += 1e-6 * np.eye(n)  # Numerical stability
    
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    np.random.seed(42)
    
    # Sample each component independently
    components = {
        'Trend (Linear)': K_trend,
        'Smooth (RBF)': K_smooth,
        'Periodic': K_periodic,
    }
    
    samples = {}
    col = 0
    for name, K in components.items():
        K_stable = K + 1e-6 * np.eye(n)
        L = np.linalg.cholesky(K_stable)
        z = np.random.randn(n)  # Save for reuse
        samples[name] = L @ z
        
        axes[0, col].plot(X, samples[name], 'b-', linewidth=2)
        axes[0, col].set_title(f'Component: {name}')
        axes[0, col].set_xlabel('x')
        axes[0, col].set_ylabel('f(x)')
        axes[0, col].grid(True, alpha=0.3)
        col += 1
    
    # Combined sample (sum of components + noise)
    noise = np.sqrt(0.1) * np.random.randn(n)
    f_total = samples['Trend (Linear)'] + samples['Smooth (RBF)'] + samples['Periodic'] + noise
    
    axes[1, 0].plot(X, f_total, 'k-', linewidth=2)
    axes[1, 0].set_title('Combined Function: Trend + Smooth + Periodic + Noise')
    axes[1, 0].set_xlabel('x')
    axes[1, 0].set_ylabel('f(x)')
    axes[1, 0].grid(True, alpha=0.3)
    
    # Show the stacked contributions
    axes[1, 1].fill_between(X, 0, samples['Trend (Linear)'], alpha=0.7, label='Trend')
    axes[1, 1].fill_between(X, samples['Trend (Linear)'], 
                            samples['Trend (Linear)'] + samples['Smooth (RBF)'], 
                            alpha=0.7, label='Smooth')
    axes[1, 1].fill_between(X, samples['Trend (Linear)'] + samples['Smooth (RBF)'],
                            samples['Trend (Linear)'] + samples['Smooth (RBF)'] + samples['Periodic'],
                            alpha=0.7, label='Periodic')
    axes[1, 1].set_title('Stacked Component Contributions')
    axes[1, 1].set_xlabel('x')
    axes[1, 1].set_ylabel('f(x)')
    axes[1, 1].legend()
    axes[1, 1].grid(True, alpha=0.3)
    
    # Covariance matrix of combined kernel
    im = axes[1, 2].imshow(K_total, extent=[0, 10, 10, 0], cmap='viridis')
    axes[1, 2].set_title('Combined Covariance Matrix')
    axes[1, 2].set_xlabel("x'")
    axes[1, 2].set_ylabel('x')
    plt.colorbar(im, ax=axes[1, 2])
    
    plt.suptitle('Additive Kernel Composition: Decomposing Complex Functions', fontsize=14)
    plt.tight_layout()
    return fig
 
def posterior_decomposition_example():
    """
    Show how posterior inference decomposes additive components.
    After observing data, we can attribute variance to each source.
    """
    # Generate synthetic data from known components
    np.random.seed(42)
    X_train = np.sort(np.random.uniform(0, 10, 30))
    
    # True components
    true_trend = 0.5 * (X_train - 5)
    true_periodic = 1.5 * np.sin(2 * np.pi * X_train / 3)
    true_smooth = np.zeros_like(X_train)  # RBF contribution
    noise = 0.5 * np.random.randn(len(X_train))
    
    y_train = true_trend + true_periodic + noise
    
    # Fit GP with additive kernel and decompose posterior
    # (Full implementation would compute posterior for each component)
    
    print("Posterior decomposition allows:")
    print("1. Estimating the trend component separately from oscillations")
    print("2. Quantifying uncertainty in each component")
    print("3. Forecasting each component independently")
    print("4. Identifying which component explains most variance")
    
    return X_train, y_train
 
demonstrate_additive_decomposition()

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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148

import numpy as np
import matplotlib.pyplot as plt
 
def locally_periodic_kernel(X1, X2, sigma_sq=1.0, period=1.0, 
                            length_scale_periodic=0.5, length_scale_decay=5.0):
    """
    Locally Periodic = Periodic × RBF
    
    Creates patterns that repeat but fade with distance.
    The RBF acts as an envelope that modulates the periodic structure.
    """
    X1 = np.atleast_1d(X1).flatten()
    X2 = np.atleast_1d(X2).flatten()
    
    # Periodic component
    dists = np.abs(np.subtract.outer(X1, X2))
    K_periodic = np.exp(-2 * (np.sin(np.pi * dists / period) / length_scale_periodic)**2)
    
    # RBF decay envelope
    sq_dists = np.subtract.outer(X1, X2)**2
    K_rbf = np.exp(-sq_dists / (2 * length_scale_decay**2))
    
    # Product kernel
    return sigma_sq * K_periodic * K_rbf
 
def demonstrate_product_vs_sum():
    """
    Compare additive vs multiplicative composition for Periodic + RBF.
    """
    X = np.linspace(0, 20, 300)
    n = len(X)
    
    # Define base kernels
    def periodic(X1, X2, period=3.0, ell=0.5):
        dists = np.abs(np.subtract.outer(X1, X2))
        return np.exp(-2 * (np.sin(np.pi * dists / period) / ell)**2)
    
    def rbf(X1, X2, ell=5.0):
        sq_dists = np.subtract.outer(X1, X2)**2
        return np.exp(-sq_dists / (2 * ell**2))
    
    K_per = periodic(X, X)
    K_rbf = rbf(X, X)
    
    # Additive: Periodic + RBF
    K_sum = K_per + K_rbf
    
    # Multiplicative: Periodic × RBF
    K_prod = K_per * K_rbf
    
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    np.random.seed(42)
    
    # Sample from periodic kernel
    L_per = np.linalg.cholesky(K_per + 1e-6*np.eye(n))
    f_per = L_per @ np.random.randn(n)
    axes[0, 0].plot(X, f_per)
    axes[0, 0].set_title('Periodic Kernel Samples')
    axes[0, 0].grid(True, alpha=0.3)
    
    # Sample from RBF kernel
    L_rbf = np.linalg.cholesky(K_rbf + 1e-6*np.eye(n))
    f_rbf = L_rbf @ np.random.randn(n)
    axes[0, 1].plot(X, f_rbf)
    axes[0, 1].set_title('RBF Kernel Samples')
    axes[0, 1].grid(True, alpha=0.3)
    
    # Covariance comparison
    axes[0, 2].plot(X, K_per[n//2, :], label='Periodic', linewidth=2)
    axes[0, 2].plot(X, K_rbf[n//2, :], label='RBF', linewidth=2)
    axes[0, 2].axvline(X[n//2], color='gray', linestyle='--', alpha=0.5)
    axes[0, 2].set_title('Covariance from Midpoint')
    axes[0, 2].legend()
    axes[0, 2].grid(True, alpha=0.3)
    
    # Sample from SUM kernel (Periodic + RBF)
    L_sum = np.linalg.cholesky(K_sum + 1e-6*np.eye(n))
    for i in range(5):
        f_sum = L_sum @ np.random.randn(n)
        axes[1, 0].plot(X, f_sum, alpha=0.7)
    axes[1, 0].set_title('SUM: Periodic + RBF\n(Independent components)')
    axes[1, 0].grid(True, alpha=0.3)
    
    # Sample from PRODUCT kernel (Periodic × RBF)
    L_prod = np.linalg.cholesky(K_prod + 1e-6*np.eye(n))
    for i in range(5):
        f_prod = L_prod @ np.random.randn(n)
        axes[1, 1].plot(X, f_prod, alpha=0.7)
    axes[1, 1].set_title('PRODUCT: Periodic × RBF\n(Locally periodic = decaying oscillations)')
    axes[1, 1].grid(True, alpha=0.3)
    
    # Covariance comparison for composed kernels
    axes[1, 2].plot(X, K_sum[n//2, :], label='Sum', linewidth=2)
    axes[1, 2].plot(X, K_prod[n//2, :], label='Product', linewidth=2)
    axes[1, 2].axvline(X[n//2], color='gray', linestyle='--', alpha=0.5)
    axes[1, 2].set_title('Composed Kernels: Covariance from Midpoint')
    axes[1, 2].legend()
    axes[1, 2].grid(True, alpha=0.3)
    
    plt.suptitle('Addition vs Multiplication: Fundamentally Different Behaviors', fontsize=14)
    plt.tight_layout()
    return fig
 
def ard_kernel_example():
    """
    Automatic Relevance Determination: Product of per-dimension kernels.
    
    This is crucial for feature selection in high-dimensional problems.
    """
    print("ARD Kernel (Automatic Relevance Determination):")
    print("="*50)
    print()
    print("k_ARD(x, x') = σ² · exp(-Σᵢ (xᵢ - x'ᵢ)² / (2ℓᵢ²))")
    print()
    print("Each dimension i has its own length-scale ℓᵢ")
    print()
    print("Feature Selection Interpretation:")
    print("- Large ℓᵢ → dimension i is irrelevant (flat kernel)")
    print("- Small ℓᵢ → dimension i is highly relevant (sensitive)")
    print()
    print("In practice: optimize all ℓᵢ via marginal likelihood")
    print("Irrelevant features get ℓᵢ → ∞ automatically")
    
    # Demonstrate in 2D with different relevances
    np.random.seed(42)
    n = 100
    X = np.random.randn(n, 2)
    
    # Dimension 1 relevant (small length-scale), Dimension 2 irrelevant (large)
    length_scales = np.array([0.5, 10.0])
    
    def ard_kernel(X1, X2, length_scales, sigma_sq=1.0):
        """ARD-RBF kernel."""
        X1_scaled = X1 / length_scales
        X2_scaled = X2 / length_scales
        sq_dists = np.sum((X1_scaled[:, np.newaxis, :] - X2_scaled[np.newaxis, :, :])**2, axis=2)
        return sigma_sq * np.exp(-sq_dists / 2)
    
    K = ard_kernel(X, X, length_scales)
    
    print(f"\nExample: length_scales = {length_scales}")
    print(f"Dimension 1 (ℓ={length_scales[0]}): RELEVANT")
    print(f"Dimension 2 (ℓ={length_scales[1]}): IRRELEVANT")
    
    return K
 
demonstrate_product_vs_sum()

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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import expit  # Sigmoid function
 
def changepoint_kernel(X1, X2, k1_fn, k2_fn, changepoint=5.0, steepness=2.0):
    """
    Changepoint kernel: different behavior before/after a transition.
    
    k_CP(x,x') = k₁(x,x')·σ(x)σ(x') + k₂(x,x')·(1-σ(x))(1-σ(x'))
    
    Parameters:
    -----------
    k1_fn, k2_fn : callables
        Kernel functions for before/after changepoint
    changepoint : float
        Location of the transition
    steepness : float
        How sharp the transition is (higher = sharper)
    """
    X1 = np.atleast_1d(X1).flatten()
    X2 = np.atleast_1d(X2).flatten()
    
    # Transition function (sigmoid)
    sigma1 = expit(steepness * (X1 - changepoint))  # n1
    sigma2 = expit(steepness * (X2 - changepoint))  # n2
    
    # Create weight matrices
    w1 = np.outer(sigma1, sigma2)          # After changepoint
    w2 = np.outer(1 - sigma1, 1 - sigma2)  # Before changepoint
    
    # Combine kernels
    K1 = k1_fn(X1, X2)
    K2 = k2_fn(X1, X2)
    
    return K1 * w2 + K2 * w1
 
def demonstrate_changepoint():
    """
    Show how a changepoint kernel transitions between regimes.
    """
    X = np.linspace(0, 10, 200)
    n = len(X)
    changepoint = 5.0
    
    # Define two different kernels for before/after
    def k_smooth(X1, X2):
        """Smooth, slowly-varying behavior."""
        sq_dists = np.subtract.outer(X1, X2)**2
        return np.exp(-sq_dists / (2 * 2.0**2))  # Large length-scale
    
    def k_wiggly(X1, X2):
        """Rapidly varying behavior."""
        sq_dists = np.subtract.outer(X1, X2)**2
        return np.exp(-sq_dists / (2 * 0.3**2))  # Small length-scale
    
    K_cp = changepoint_kernel(X, X, k_smooth, k_wiggly, changepoint=5.0, steepness=2.0)
    K_cp += 1e-6 * np.eye(n)
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 4))
    
    # Samples from changepoint kernel
    np.random.seed(42)
    L = np.linalg.cholesky(K_cp)
    for i in range(5):
        sample = L @ np.random.randn(n)
        axes[0].plot(X, sample, alpha=0.7)
    axes[0].axvline(changepoint, color='red', linestyle='--', label='Changepoint')
    axes[0].set_title('Changepoint Kernel: Smooth → Wiggly')
    axes[0].set_xlabel('x')
    axes[0].set_ylabel('f(x)')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # Covariance matrix
    im = axes[1].imshow(K_cp, extent=[0, 10, 10, 0], cmap='viridis')
    axes[1].axhline(changepoint, color='red', linestyle='--')
    axes[1].axvline(changepoint, color='red', linestyle='--')
    axes[1].set_title('Covariance Matrix')
    axes[1].set_xlabel("x'")
    axes[1].set_ylabel('x')
    plt.colorbar(im, ax=axes[1])
    
    # Transition weights
    sigma = expit(2.0 * (X - changepoint))
    axes[2].plot(X, 1 - sigma, 'b-', linewidth=2, label='Weight for k₁ (smooth)')
    axes[2].plot(X, sigma, 'r-', linewidth=2, label='Weight for k₂ (wiggly)')
    axes[2].axvline(changepoint, color='gray', linestyle='--', alpha=0.5)
    axes[2].set_title('Transition Weights')
    axes[2].set_xlabel('x')
    axes[2].set_ylabel('Weight')
    axes[2].legend()
    axes[2].grid(True, alpha=0.3)
    
    plt.suptitle('Changepoint Kernel: Regime-Switching Behavior', fontsize=14)
    plt.tight_layout()
    return fig
 
def multi_output_kernel():
    """
    Multi-output GPs using kernel composition.
    
    For modeling multiple related outputs (y₁, y₂, ...) jointly:
    K(x,x') = B ⊗ k(x,x')
    
    Where B is a positive semi-definite matrix of output correlations.
    """
    print("Multi-Output (Vector-Valued) GPs via Kernel Composition")
    print("="*55)
    print()
    print("For D outputs, construct a DxD output correlation matrix B:")
    print("  B[i,j] = correlation between output i and output j")
    print()
    print("The full kernel is a Kronecker product:")
    print("  K_full((x,i), (x',j)) = B[i,j] × k_base(x, x')")
    print()
    print("This models:")
    print("  - Each output as a GP with kernel k_base")
    print("  - Cross-output correlations via B")
    print()
    print("Example: Modeling temperature and humidity jointly")
    print("  - k_base = Matérn captures spatial correlation")
    print("  - B captures that temperature and humidity are related")
    
    # Simple example: 2 outputs
    n_locations = 5
    n_outputs = 2
    
    # Correlation between outputs
    B = np.array([[1.0, 0.7],
                  [0.7, 1.0]])
    
    # Base kernel for locations
    X = np.array([0, 1, 2, 3, 4])
    K_base = np.exp(-np.subtract.outer(X, X)**2)
    
    # Full kernel via Kronecker product
    K_full = np.kron(B, K_base)
    
    print(f"\nBase kernel shape: {K_base.shape}")
    print(f"Output correlation matrix B shape: {B.shape}")
    print(f"Full multi-output kernel shape: {K_full.shape}")
    
    return K_full
 
demonstrate_changepoint()

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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

import numpy as np
from itertools import product
 
class KernelExpression:
    """
    Abstract syntax tree for kernel expressions.
    Enables composition and natural language interpretation.
    """
    pass
 
class BaseKernel(KernelExpression):
    def __init__(self, name, description):
        self.name = name
        self.description = description
    
    def __repr__(self):
        return self.name
    
    def describe(self):
        return self.description
 
class SumKernel(KernelExpression):
    def __init__(self, left, right):
        self.left = left
        self.right = right
    
    def __repr__(self):
        return f"({self.left} + {self.right})"
    
    def describe(self):
        return f"{self.left.describe()}, plus independently {self.right.describe()}"
 
class ProductKernel(KernelExpression):
    def __init__(self, left, right):
        self.left = left
        self.right = right
    
    def __repr__(self):
        return f"({self.left} × {self.right})"
    
    def describe(self):
        return f"{self.right.describe()}, but only where {self.left.describe()} applies"
 
# Define base kernels with interpretations
RBF = BaseKernel("RBF", "smooth, locally correlated variations")
LIN = BaseKernel("Linear", "a linear trend")
PER = BaseKernel("Periodic", "repeating periodic patterns")
WN = BaseKernel("WhiteNoise", "independent noise at each point")
 
def generate_depth_1_kernels(bases):
    """Generate all kernels of depth 1 (single base)."""
    return bases.copy()
 
def generate_depth_2_kernels(bases):
    """Generate all kernels of depth 2 (base op base)."""
    results = []
    for b1, b2 in product(bases, repeat=2):
        results.append(SumKernel(b1, b2))
        results.append(ProductKernel(b1, b2))
    return results
 
def demonstrate_grammar():
    """Show how the grammar generates interpretable expressions."""
    print("Kernel Grammar: Compositional Structure Discovery")
    print("="*55)
    print()
    
    bases = [RBF, LIN, PER]
    
    print("Base Kernels:")
    for b in bases:
        print(f"  {b.name}: {b.description}")
    
    print()
    print("Example Composed Kernels and Their Interpretations:")
    print("-"*55)
    
    # Manual examples of useful compositions
    examples = [
        SumKernel(LIN, RBF),
        ProductKernel(PER, RBF),
        SumKernel(LIN, ProductKernel(PER, RBF)),
        ProductKernel(LIN, PER),
        SumKernel(SumKernel(LIN, PER), RBF),
    ]
    
    for expr in examples:
        print(f"\n  Expression: {expr}")
        print(f"  Meaning: The function exhibits {expr.describe()}")
    
    print()
    print("-"*55)
    print("KEY INSIGHT: The grammar enables automatic model search!")
    print("Algorithms can enumerate expressions and select the best via")
    print("marginal likelihood, then explain findings in natural language.")
 
def greedy_kernel_search(X, y, max_depth=3, bases=None):
    """
    Simplified greedy kernel structure search.
    
    Real implementations (like GPy's kernel_search or AutoGP) use:
    - Efficient marginal likelihood computation
    - Beam search or MCTS for exploration
    - Regularization to prevent overfitting
    """
    if bases is None:
        bases = [RBF, LIN, PER]
    
    print("Greedy Kernel Search Algorithm:")
    print("="*40)
    print("1. Start with each base kernel")
    print("2. Evaluate marginal likelihood for current best")
    print("3. Try all single expansions: k+base, k×base")
    print("4. Keep best expansion if it improves")
    print("5. Repeat until max depth or no improvement")
    print()
    print("This finds interpretable structure automatically!")
    
    # Pseudocode for actual implementation:
    """
    current_best = argmax(bases, key=compute_log_ml)
    for depth in range(max_depth):
        candidates = []
        for base in bases:
            candidates.append(SumKernel(current_best, base))
            candidates.append(ProductKernel(current_best, base))
        new_best = argmax(candidates, key=compute_log_ml)
        if log_ml(new_best) <= log_ml(current_best) + threshold:
            break  # No improvement
        current_best = new_best
    return current_best
    """
 
demonstrate_grammar()