Gaussian Process Regression

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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310

"""
Complete Gaussian Process Regression Implementation
===================================================
A production-quality implementation covering all aspects of GP regression.
"""
 
import numpy as np
from scipy.linalg import cholesky, solve_triangular, cho_solve
from scipy.optimize import minimize
from typing import Tuple, Optional, Callable, Dict, Any
 
class GaussianProcessRegressor:
    """
    Gaussian Process Regressor with full functionality.
    
    Features:
    - RBF (Squared Exponential) kernel with customizable parameters
    - Automatic hyperparameter optimization via marginal likelihood
    - Numerically stable Cholesky-based inference
    - Both pointwise and full covariance predictions
    - Proper handling of noise and jitter
    """
    
    def __init__(
        self,
        kernel: str = 'rbf',
        length_scale: float = 1.0,
        signal_variance: float = 1.0,
        noise_variance: float = 1e-6,
        jitter: float = 1e-8,
        optimize_hyperparameters: bool = True,
        n_restarts: int = 5,
    ):
        """
        Initialize the GP regressor.
        
        Parameters:
        -----------
        kernel : str
            Kernel type. Currently supports 'rbf'.
        length_scale : float
            Initial length scale for the kernel.
        signal_variance : float
            Initial signal variance (kernel amplitude squared).
        noise_variance : float
            Observation noise variance. Set to small value for near-interpolation.
        jitter : float
            Small value added to diagonal for numerical stability.
        optimize_hyperparameters : bool
            Whether to optimize hyperparameters during fit.
        n_restarts : int
            Number of random restarts for optimization.
        """
        self.kernel_type = kernel
        self.length_scale = length_scale
        self.signal_variance = signal_variance
        self.noise_variance = noise_variance
        self.jitter = jitter
        self.optimize_hyperparameters = optimize_hyperparameters
        self.n_restarts = n_restarts
        
        # Fitted state
        self._is_fitted = False
        self._X_train = None
        self._y_train = None
        self._L = None  # Cholesky factor
        self._alpha = None  # Dual coefficients
        
    def _kernel(self, X1: np.ndarray, X2: np.ndarray) -> np.ndarray:
        """
        Compute the kernel matrix between X1 and X2.
        """
        if self.kernel_type == 'rbf':
            # Squared distances: ||x1 - x2||^2
            sq_dist = (
                np.sum(X1**2, axis=1, keepdims=True) 
                - 2 * X1 @ X2.T 
                + np.sum(X2**2, axis=1)
            )
            return self.signal_variance * np.exp(-0.5 * sq_dist / self.length_scale**2)
        else:
            raise ValueError(f"Unknown kernel: {self.kernel_type}")
    
    def _compute_kernel_with_params(
        self, X1: np.ndarray, X2: np.ndarray, 
        length_scale: float, signal_variance: float
    ) -> np.ndarray:
        """Compute kernel with specific hyperparameters."""
        sq_dist = (
            np.sum(X1**2, axis=1, keepdims=True) 
            - 2 * X1 @ X2.T 
            + np.sum(X2**2, axis=1)
        )
        return signal_variance * np.exp(-0.5 * sq_dist / length_scale**2)
    
    def _neg_log_marginal_likelihood(
        self, log_theta: np.ndarray, X: np.ndarray, y: np.ndarray
    ) -> float:
        """
        Compute negative log marginal likelihood for optimization.
        
        Parameters:
        -----------
        log_theta : array
            [log(length_scale), log(signal_variance), log(noise_variance)]
        """
        length_scale = np.exp(log_theta[0])
        signal_variance = np.exp(log_theta[1])
        noise_variance = np.exp(log_theta[2])
        
        n = len(y)
        
        # Compute kernel matrix
        K = self._compute_kernel_with_params(X, X, length_scale, signal_variance)
        K_y = K + (noise_variance + self.jitter) * np.eye(n)
        
        # Cholesky decomposition
        try:
            L = cholesky(K_y, lower=True)
        except np.linalg.LinAlgError:
            return 1e10  # Return large value if not positive definite
        
        # Solve for alpha
        alpha = cho_solve((L, True), y)
        
        # Log marginal likelihood components
        data_fit = -0.5 * y @ alpha
        complexity = -np.sum(np.log(np.diag(L)))  # -0.5 * log|K_y|
        constant = -0.5 * n * np.log(2 * np.pi)
        
        return -(data_fit + complexity + constant)
    
    def fit(self, X: np.ndarray, y: np.ndarray) -> 'GaussianProcessRegressor':
        """
        Fit the GP to training data.
        
        Parameters:
        -----------
        X : array of shape (n_samples, n_features)
            Training inputs.
        y : array of shape (n_samples,)
            Training targets.
        
        Returns:
        --------
        self : GaussianProcessRegressor
            Fitted estimator.
        """
        X = np.atleast_2d(X)
        y = np.asarray(y).ravel()
        
        if X.shape[0] != len(y):
            raise ValueError("X and y must have same number of samples")
        
        # Store training data
        self._X_train = X.copy()
        self._y_train = y.copy()
        n = len(y)
        
        # Optimize hyperparameters if requested
        if self.optimize_hyperparameters and n > 1:
            best_nlml = np.inf
            best_theta = None
            
            for _ in range(self.n_restarts):
                # Random initialization in log-space
                theta0 = np.array([
                    np.log(self.length_scale) + 0.5 * np.random.randn(),
                    np.log(self.signal_variance) + 0.5 * np.random.randn(),
                    np.log(self.noise_variance) + 0.5 * np.random.randn(),
                ])
                
                result = minimize(
                    self._neg_log_marginal_likelihood,
                    theta0,
                    args=(X, y),
                    method='L-BFGS-B',
                    bounds=[(-10, 10), (-10, 10), (-15, 5)],
                )
                
                if result.fun < best_nlml:
                    best_nlml = result.fun
                    best_theta = result.x
            
            if best_theta is not None:
                self.length_scale = np.exp(best_theta[0])
                self.signal_variance = np.exp(best_theta[1])
                self.noise_variance = np.exp(best_theta[2])
        
        # Compute kernel matrix and decomposition
        K = self._kernel(X, X)
        K_y = K + (self.noise_variance + self.jitter) * np.eye(n)
        
        # Cholesky decomposition
        self._L = cholesky(K_y, lower=True)
        
        # Compute dual coefficients
        self._alpha = cho_solve((self._L, True), y)
        
        self._is_fitted = True
        return self
    
    def predict(
        self, X: np.ndarray, return_std: bool = False, return_cov: bool = False
    ) -> Tuple[np.ndarray, ...]:
        """
        Make predictions at test points.
        
        Parameters:
        -----------
        X : array of shape (n_samples, n_features)
            Test inputs.
        return_std : bool
            Whether to return standard deviation.
        return_cov : bool
            Whether to return full covariance matrix.
        
        Returns:
        --------
        y_mean : array of shape (n_samples,)
            Predictive means.
        y_std : array of shape (n_samples,) (if return_std=True)
            Predictive standard deviations for latent function.
        y_cov : array of shape (n_samples, n_samples) (if return_cov=True)
            Predictive covariance for latent function.
        """
        if not self._is_fitted:
            raise RuntimeError("Must call fit() before predict()")
        
        X = np.atleast_2d(X)
        
        # Cross-covariance
        K_star = self._kernel(self._X_train, X)
        
        # Predictive mean
        y_mean = K_star.T @ self._alpha
        
        result = [y_mean]
        
        if return_std or return_cov:
            # Solve L @ V = K_star
            V = solve_triangular(self._L, K_star, lower=True)
            
            if return_cov:
                # Full covariance
                K_ss = self._kernel(X, X)
                y_cov = K_ss - V.T @ V
                # Ensure symmetry
                y_cov = 0.5 * (y_cov + y_cov.T)
                result.append(y_cov)
            
            if return_std:
                # Pointwise variance
                K_ss_diag = self.signal_variance * np.ones(X.shape[0])
                y_var = K_ss_diag - np.sum(V**2, axis=0)
                # Clamp to avoid numerical issues
                y_var = np.maximum(y_var, 0)
                result.append(np.sqrt(y_var))
        
        return tuple(result) if len(result) > 1 else result[0]
    
    def sample_y(self, X: np.ndarray, n_samples: int = 1) -> np.ndarray:
        """
        Sample functions from the posterior distribution.
        
        Parameters:
        -----------
        X : array of shape (n_points, n_features)
            Points at which to sample.
        n_samples : int
            Number of function samples.
        
        Returns:
        --------
        samples : array of shape (n_points, n_samples)
            Function samples.
        """
        y_mean, y_cov = self.predict(X, return_cov=True)
        
        # Cholesky of posterior covariance
        L_cov = cholesky(y_cov + self.jitter * np.eye(len(y_mean)), lower=True)
        
        # Sample
        z = np.random.randn(len(y_mean), n_samples)
        samples = y_mean[:, None] + L_cov @ z
        
        return samples
    
    def log_marginal_likelihood(self) -> float:
        """
        Compute log marginal likelihood of the fitted model.
        """
        if not self._is_fitted:
            raise RuntimeError("Must call fit() before log_marginal_likelihood()")
        
        n = len(self._y_train)
        
        data_fit = -0.5 * self._y_train @ self._alpha
        complexity = -np.sum(np.log(np.diag(self._L)))
        constant = -0.5 * n * np.log(2 * np.pi)
        
        return data_fit + complexity + constant
    
    def get_params(self) -> Dict[str, Any]:
        """Get current hyperparameters."""
        return {
            'length_scale': self.length_scale,
            'signal_variance': self.signal_variance,
            'noise_variance': self.noise_variance,
        }

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

"""
Gaussian Process Regression - Usage Examples
=============================================
"""
 
import numpy as np
 
# Assuming GaussianProcessRegressor is imported from the previous code
 
np.random.seed(42)
 
# =============================================================================
# Example 1: Basic Regression
# =============================================================================
print("Example 1: Basic Regression")
print("=" * 50)
 
# Generate synthetic data
def true_function(x):
    return np.sin(2 * x) + 0.5 * np.cos(4 * x)
 
X_train = np.sort(np.random.uniform(-3, 3, 15)).reshape(-1, 1)
y_train = true_function(X_train.flatten()) + 0.2 * np.random.randn(15)
 
# Fit GP
gp = GaussianProcessRegressor(optimize_hyperparameters=True, n_restarts=3)
gp.fit(X_train, y_train)
 
print(f"Optimized hyperparameters:")
params = gp.get_params()
for k, v in params.items():
    print(f"  {k}: {v:.4f}")
 
# Predict
X_test = np.linspace(-4, 4, 100).reshape(-1, 1)
y_mean, y_std = gp.predict(X_test, return_std=True)
 
print(f"\nPrediction range: [{y_mean.min():.3f}, {y_mean.max():.3f}]")
print(f"Uncertainty range: [{y_std.min():.4f}, {y_std.max():.3f}]")
print(f"Log marginal likelihood: {gp.log_marginal_likelihood():.2f}")
 
# =============================================================================
# Example 2: Posterior Sampling
# =============================================================================
print("\nExample 2: Posterior Sampling")
print("=" * 50)
 
# Sample functions from posterior
samples = gp.sample_y(X_test, n_samples=5)
 
print(f"Generated {samples.shape[1]} posterior function samples")
print(f"Sample 1 range: [{samples[:, 0].min():.3f}, {samples[:, 0].max():.3f}]")
 
# =============================================================================
# Example 3: Interpolation (Low Noise)
# =============================================================================
print("\nExample 3: Interpolation (Low Noise)")
print("=" * 50)
 
# Create GP with very low noise for interpolation
gp_interp = GaussianProcessRegressor(
    noise_variance=1e-10,
    optimize_hyperparameters=False,
    length_scale=1.0,
    signal_variance=1.0,
)
gp_interp.fit(X_train, y_train)
 
y_at_train = gp_interp.predict(X_train)
residuals = np.abs(y_at_train - y_train)
 
print(f"Max residual at training points: {residuals.max():.2e}")
print("(Should be near machine precision for interpolation)")
 
# =============================================================================
# Example 4: Multi-Dimensional Input
# =============================================================================
print("\nExample 4: Multi-Dimensional Input")
print("=" * 50)
 
# 2D Branin function (classic optimization benchmark)
def branin(X):
    x1, x2 = X[:, 0], X[:, 1]
    a, b, c = 1, 5.1/(4*np.pi**2), 5/np.pi
    r, s, t = 6, 10, 1/(8*np.pi)
    return a*(x2 - b*x1**2 + c*x1 - r)**2 + s*(1-t)*np.cos(x1) + s
 
# Generate training data
X_train_2d = np.random.uniform(low=[-5, 0], high=[10, 15], size=(30, 2))
y_train_2d = branin(X_train_2d) + 10 * np.random.randn(30)
 
# Fit GP
gp_2d = GaussianProcessRegressor(optimize_hyperparameters=True, n_restarts=3)
gp_2d.fit(X_train_2d, y_train_2d)
 
print(f"2D GP hyperparameters:")
for k, v in gp_2d.get_params().items():
    print(f"  {k}: {v:.4f}")
 
# Predict on grid
xx, yy = np.meshgrid(np.linspace(-5, 10, 20), np.linspace(0, 15, 20))
X_grid = np.column_stack([xx.ravel(), yy.ravel()])
y_pred, y_std_2d = gp_2d.predict(X_grid, return_std=True)
 
print(f"2D prediction range: [{y_pred.min():.1f}, {y_pred.max():.1f}]")
print(f"2D uncertainty range: [{y_std_2d.min():.2f}, {y_std_2d.max():.2f}]")
 
# =============================================================================
# Example 5: Extrapolation Behavior  
# =============================================================================
print("\nExample 5: Extrapolation Behavior")
print("=" * 50)
 
# Fit GP on limited range
X_limited = np.linspace(-2, 2, 10).reshape(-1, 1)
y_limited = np.sin(X_limited.flatten())
 
gp_extrap = GaussianProcessRegressor(
    optimize_hyperparameters=False,
    length_scale=1.0,
    signal_variance=1.0,
    noise_variance=0.01,
)
gp_extrap.fit(X_limited, y_limited)
 
# Predict in and out of training range
X_in = np.array([[0.0]])  # In range
X_out = np.array([[5.0]])  # Out of range
 
y_in_mean, y_in_std = gp_extrap.predict(X_in, return_std=True)
y_out_mean, y_out_std = gp_extrap.predict(X_out, return_std=True)
 
print(f"In-range (x=0): mean={y_in_mean[0]:.3f}, std={y_in_std[0]:.4f}")
print(f"Out-of-range (x=5): mean={y_out_mean[0]:.3f}, std={y_out_std[0]:.3f}")
print("Note: Extrapolation reverts to prior (mean→0, std→signal_variance)")