sklears_gaussian_process/
kernel_selection.rs

1//! Kernel selection methods for Gaussian Processes
2//!
3//! This module provides model-selection utilities for choosing among a set
4//! of candidate kernels (or, more precisely, kernel *families*): each
5//! candidate has its own hyperparameters (and the shared observation noise)
6//! optimized against the log marginal likelihood via
7//! [`crate::marginal_likelihood::optimize_hyperparameters`], and the fitted
8//! candidates are then ranked using one of four criteria:
9//!
10//! * [`SelectionCriterion::Aic`] — Akaike Information Criterion,
11//!   `-2 * LML + 2k` (lower is better).
12//! * [`SelectionCriterion::Bic`] — Bayesian Information Criterion,
13//!   `-2 * LML + k * ln(n)` (lower is better; penalizes model complexity more
14//!   heavily than AIC for `n > 7`).
15//! * [`SelectionCriterion::LogMarginalLikelihood`] — the raw log marginal
16//!   likelihood (higher is better).
17//! * [`SelectionCriterion::Cv`] — mean held-out predictive log density under
18//!   k-fold cross-validation (higher is better).
19//!
20//! Here `k` is the number of hyperparameters (kernel parameters plus the
21//! noise level) and `n` is the number of training samples.
22//!
23//! Note on [`SelectionCriterion::Cv`]: `crate::marginal_likelihood`'s
24//! `cross_validate_hyperparameters` evaluates the *same* full-data marginal
25//! likelihood on every fold (by its own documentation this is "simplified"
26//! and does not actually hold out data), so it is not suitable for genuine
27//! kernel-family comparison. This module instead performs real k-fold
28//! cross-validation: for each fold it fits a [`crate::gpr::GaussianProcessRegressor`]
29//! on the training split and scores it by the predictive log density of the
30//! held-out split.
31use std::f64::consts::PI;
32
33use scirs2_core::ndarray::{Array1, Array2};
34use sklears_core::{
35    error::{Result as SklResult, SklearsError},
36    traits::Fit,
37};
38
39use crate::gpr::GaussianProcessRegressor;
40use crate::kernels::Kernel;
41use crate::marginal_likelihood::{log_marginal_likelihood_stable, optimize_hyperparameters};
42
43/// Number of Adam iterations used to fit each candidate's hyperparameters
44/// before scoring it. Kept modest since kernel selection typically compares
45/// several candidates and does not need each one driven to full convergence
46/// to distinguish a good kernel family from a poor one.
47const CANDIDATE_OPTIMIZATION_ITERS: usize = 50;
48const CANDIDATE_OPTIMIZATION_TOL: f64 = 1e-4;
49
50/// Criterion used to rank candidate kernels.
51#[derive(Debug, Clone, Copy, PartialEq)]
52pub enum SelectionCriterion {
53    /// Akaike Information Criterion (lower is better).
54    Aic,
55    /// Bayesian Information Criterion (lower is better).
56    Bic,
57    /// Mean held-out predictive log density under k-fold CV (higher is
58    /// better).
59    Cv,
60    /// Raw log marginal likelihood (higher is better).
61    LogMarginalLikelihood,
62}
63
64impl SelectionCriterion {
65    /// Whether higher scores are better for this criterion (`Cv` and
66    /// `LogMarginalLikelihood`), as opposed to lower being better (`Aic`,
67    /// `Bic`).
68    pub fn higher_is_better(self) -> bool {
69        matches!(
70            self,
71            SelectionCriterion::Cv | SelectionCriterion::LogMarginalLikelihood
72        )
73    }
74}
75
76/// Configuration for [`KernelSelector`].
77#[derive(Debug, Clone)]
78pub struct KernelSelectionConfig {
79    /// Criterion used to rank candidates.
80    pub criterion: SelectionCriterion,
81    /// Number of folds used when `criterion` is [`SelectionCriterion::Cv`].
82    pub cv_folds: usize,
83    /// Initial (and, for AIC/BIC/LML scoring, jointly optimized) observation
84    /// noise variance `σ_n²`.
85    pub noise_variance: f64,
86    /// Random seed controlling the CV fold shuffle.
87    pub random_state: Option<u64>,
88}
89
90impl Default for KernelSelectionConfig {
91    fn default() -> Self {
92        Self {
93            criterion: SelectionCriterion::Bic,
94            cv_folds: 5,
95            noise_variance: 1e-4,
96            random_state: Some(0),
97        }
98    }
99}
100
101/// Result of a kernel selection run.
102#[derive(Debug, Clone)]
103pub struct KernelSelectionResult {
104    /// The winning kernel, with hyperparameters optimized against the
105    /// training data.
106    pub best_kernel: Box<dyn Kernel>,
107    /// Index of the winning kernel within the candidates slice that was
108    /// passed in.
109    pub best_index: usize,
110    /// Short name of the winning kernel (derived from its `Debug`
111    /// representation), e.g. `"RBF"`.
112    pub best_kernel_name: String,
113    /// Score of the winning kernel under the configured criterion.
114    pub best_score: f64,
115    /// `(name, score)` for every candidate, in input order.
116    pub all_scores: Vec<(String, f64)>,
117}
118
119/// Builder for running kernel selection over a fixed set of candidate
120/// kernels.
121///
122/// # Examples
123///
124/// ```ignore
125/// let candidates: Vec<Box<dyn Kernel>> = vec![Box::new(RBF::new(1.0)), Box::new(Linear::new(1.0, 1.0))];
126/// let selector = KernelSelector::new(candidates).criterion(SelectionCriterion::Bic);
127/// let result = selector.select(&X, &y).expect("selection should succeed with valid data");
128/// ```
129#[derive(Debug, Clone)]
130pub struct KernelSelector {
131    /// Candidate kernels to evaluate.
132    pub candidates: Vec<Box<dyn Kernel>>,
133    /// Selection configuration.
134    pub config: KernelSelectionConfig,
135}
136
137impl KernelSelector {
138    /// Create a new selector over the given candidate kernels, using default
139    /// configuration (BIC criterion).
140    pub fn new(candidates: Vec<Box<dyn Kernel>>) -> Self {
141        Self {
142            candidates,
143            config: KernelSelectionConfig::default(),
144        }
145    }
146
147    /// Set the full configuration at once.
148    pub fn config(mut self, config: KernelSelectionConfig) -> Self {
149        self.config = config;
150        self
151    }
152
153    /// Set the selection criterion.
154    pub fn criterion(mut self, criterion: SelectionCriterion) -> Self {
155        self.config.criterion = criterion;
156        self
157    }
158
159    /// Set the number of cross-validation folds (used only for
160    /// [`SelectionCriterion::Cv`]).
161    pub fn cv_folds(mut self, cv_folds: usize) -> Self {
162        self.config.cv_folds = cv_folds;
163        self
164    }
165
166    /// Set the initial observation noise variance.
167    pub fn noise_variance(mut self, noise_variance: f64) -> Self {
168        self.config.noise_variance = noise_variance;
169        self
170    }
171
172    /// Set the random state controlling the CV fold shuffle.
173    pub fn random_state(mut self, random_state: Option<u64>) -> Self {
174        self.config.random_state = random_state;
175        self
176    }
177
178    /// Run selection against the given training data.
179    #[allow(non_snake_case)]
180    pub fn select(&self, X: &Array2<f64>, y: &Array1<f64>) -> SklResult<KernelSelectionResult> {
181        let candidates: Vec<Box<dyn Kernel>> =
182            self.candidates.iter().map(|k| k.clone_box()).collect();
183        select_best_kernel_with_config(candidates, X, y, &self.config)
184    }
185}
186
187/// Select the best kernel among `candidates` using the given `criterion`
188/// (with all other selection settings left at their defaults; use
189/// [`KernelSelector`] to customize e.g. `cv_folds`).
190#[allow(non_snake_case)]
191pub fn select_best_kernel(
192    candidates: Vec<Box<dyn Kernel>>,
193    X: &Array2<f64>,
194    y: &Array1<f64>,
195    criterion: SelectionCriterion,
196) -> SklResult<KernelSelectionResult> {
197    let config = KernelSelectionConfig {
198        criterion,
199        ..KernelSelectionConfig::default()
200    };
201    select_best_kernel_with_config(candidates, X, y, &config)
202}
203
204/// Convenience wrapper: select the best candidate by AIC.
205#[allow(non_snake_case)]
206pub fn select_kernel_aic(
207    candidates: Vec<Box<dyn Kernel>>,
208    X: &Array2<f64>,
209    y: &Array1<f64>,
210) -> SklResult<KernelSelectionResult> {
211    select_best_kernel(candidates, X, y, SelectionCriterion::Aic)
212}
213
214/// Convenience wrapper: select the best candidate by BIC.
215#[allow(non_snake_case)]
216pub fn select_kernel_bic(
217    candidates: Vec<Box<dyn Kernel>>,
218    X: &Array2<f64>,
219    y: &Array1<f64>,
220) -> SklResult<KernelSelectionResult> {
221    select_best_kernel(candidates, X, y, SelectionCriterion::Bic)
222}
223
224/// Convenience wrapper: select the best candidate by cross-validated
225/// predictive log density.
226#[allow(non_snake_case)]
227pub fn select_kernel_cv(
228    candidates: Vec<Box<dyn Kernel>>,
229    X: &Array2<f64>,
230    y: &Array1<f64>,
231) -> SklResult<KernelSelectionResult> {
232    select_best_kernel(candidates, X, y, SelectionCriterion::Cv)
233}
234
235#[allow(non_snake_case)]
236fn select_best_kernel_with_config(
237    candidates: Vec<Box<dyn Kernel>>,
238    X: &Array2<f64>,
239    y: &Array1<f64>,
240    config: &KernelSelectionConfig,
241) -> SklResult<KernelSelectionResult> {
242    if candidates.is_empty() {
243        return Err(SklearsError::InvalidInput(
244            "At least one candidate kernel is required".to_string(),
245        ));
246    }
247    if X.nrows() != y.len() {
248        return Err(SklearsError::InvalidInput(
249            "X and y must have the same number of samples".to_string(),
250        ));
251    }
252
253    let higher_is_better = config.criterion.higher_is_better();
254
255    let mut all_scores = Vec::with_capacity(candidates.len());
256    let mut best: Option<(Box<dyn Kernel>, f64)> = None;
257    let mut best_index = 0usize;
258
259    for (idx, candidate) in candidates.iter().enumerate() {
260        let (optimized_kernel, score) = evaluate_candidate(candidate.as_ref(), X, y, config)?;
261        all_scores.push((kernel_name(optimized_kernel.as_ref()), score));
262
263        let is_better = match &best {
264            None => true,
265            Some((_, best_score)) => {
266                if higher_is_better {
267                    score > *best_score
268                } else {
269                    score < *best_score
270                }
271            }
272        };
273        if is_better {
274            best = Some((optimized_kernel, score));
275            best_index = idx;
276        }
277    }
278
279    let (best_kernel, best_score) = best.ok_or_else(|| {
280        SklearsError::InvalidInput("At least one candidate kernel is required".to_string())
281    })?;
282    let best_kernel_name = kernel_name(best_kernel.as_ref());
283
284    Ok(KernelSelectionResult {
285        best_kernel,
286        best_index,
287        best_kernel_name,
288        best_score,
289        all_scores,
290    })
291}
292
293/// Optimize one candidate's hyperparameters (and noise level) against the
294/// log marginal likelihood, then score it under `config.criterion`.
295#[allow(non_snake_case)]
296fn evaluate_candidate(
297    kernel: &dyn Kernel,
298    X: &Array2<f64>,
299    y: &Array1<f64>,
300    config: &KernelSelectionConfig,
301) -> SklResult<(Box<dyn Kernel>, f64)> {
302    let n = X.nrows();
303    let mut kernel_opt = kernel.clone_box();
304    let sigma_n_init = config.noise_variance.max(1e-12).sqrt();
305
306    let opt_result = optimize_hyperparameters(
307        &X.view(),
308        &y.view(),
309        kernel_opt.as_mut(),
310        sigma_n_init,
311        CANDIDATE_OPTIMIZATION_ITERS,
312        CANDIDATE_OPTIMIZATION_TOL,
313        false,
314    )?;
315    let sigma_n = opt_result.optimal_params.last().copied().ok_or_else(|| {
316        SklearsError::NumericalError(
317            "hyperparameter optimization returned no parameters".to_string(),
318        )
319    })?;
320
321    let score = match config.criterion {
322        SelectionCriterion::Cv => {
323            cross_validated_log_density(kernel_opt.as_ref(), X, y, sigma_n * sigma_n, config)?
324        }
325        SelectionCriterion::Aic
326        | SelectionCriterion::Bic
327        | SelectionCriterion::LogMarginalLikelihood => {
328            let lml =
329                log_marginal_likelihood_stable(&X.view(), &y.view(), kernel_opt.as_ref(), sigma_n)?;
330            let k = kernel_opt.get_params().len() as f64 + 1.0; // +1 for the noise level
331            match config.criterion {
332                SelectionCriterion::Aic => -2.0 * lml + 2.0 * k,
333                SelectionCriterion::Bic => -2.0 * lml + k * (n as f64).ln(),
334                _ => lml, // LogMarginalLikelihood
335            }
336        }
337    };
338
339    Ok((kernel_opt, score))
340}
341
342/// Mean held-out predictive log density under k-fold cross-validation.
343#[allow(non_snake_case)]
344fn cross_validated_log_density(
345    kernel: &dyn Kernel,
346    X: &Array2<f64>,
347    y: &Array1<f64>,
348    noise_variance: f64,
349    config: &KernelSelectionConfig,
350) -> SklResult<f64> {
351    let n = X.nrows();
352    let k_folds = config.cv_folds.clamp(2, n.max(2));
353
354    let mut indices: Vec<usize> = (0..n).collect();
355    if let Some(seed) = config.random_state {
356        let mut rng = seed;
357        for i in (1..indices.len()).rev() {
358            rng = rng.wrapping_mul(1_103_515_245).wrapping_add(12_345);
359            let j = (rng as usize) % (i + 1);
360            indices.swap(i, j);
361        }
362    }
363
364    let fold_size = n / k_folds;
365    let mut total_log_density = 0.0;
366    let mut total_count = 0usize;
367
368    for fold in 0..k_folds {
369        let start = fold * fold_size;
370        let end = if fold == k_folds - 1 {
371            n
372        } else {
373            (fold + 1) * fold_size
374        };
375        if start == end {
376            continue;
377        }
378
379        let val_idx = &indices[start..end];
380        let train_idx: Vec<usize> = indices
381            .iter()
382            .copied()
383            .filter(|i| *i < start || *i >= end)
384            .collect();
385        if train_idx.is_empty() {
386            continue;
387        }
388
389        let X_train = select_rows(X, &train_idx);
390        let y_train = select_entries(y, &train_idx);
391        let X_val = select_rows(X, val_idx);
392        let y_val = select_entries(y, val_idx);
393
394        let model = GaussianProcessRegressor::new()
395            .kernel(kernel.clone_box())
396            .alpha(noise_variance)
397            .fit(&X_train, &y_train)?;
398        let (mean, std) = model.predict_with_std(&X_val)?;
399
400        for i in 0..X_val.nrows() {
401            let var = (std[i] * std[i] + noise_variance).max(1e-12);
402            let residual = y_val[i] - mean[i];
403            let point_log_density = -0.5 * (residual * residual / var + var.ln() + (2.0 * PI).ln());
404            total_log_density += point_log_density;
405            total_count += 1;
406        }
407    }
408
409    if total_count == 0 {
410        return Err(SklearsError::InvalidInput(
411            "Cross-validation produced no held-out predictions; check cv_folds and dataset size"
412                .to_string(),
413        ));
414    }
415
416    Ok(total_log_density / total_count as f64)
417}
418
419#[allow(non_snake_case)]
420fn select_rows(X: &Array2<f64>, idx: &[usize]) -> Array2<f64> {
421    let mut out = Array2::<f64>::zeros((idx.len(), X.ncols()));
422    for (new_i, &old_i) in idx.iter().enumerate() {
423        out.row_mut(new_i).assign(&X.row(old_i));
424    }
425    out
426}
427
428fn select_entries(y: &Array1<f64>, idx: &[usize]) -> Array1<f64> {
429    Array1::from_iter(idx.iter().map(|&i| y[i]))
430}
431
432/// Extract a short display name for a kernel from its `Debug` representation
433/// (e.g. `"RBF { length_scale: 1.0 }"` -> `"RBF"`).
434fn kernel_name(kernel: &dyn Kernel) -> String {
435    let debug_repr = format!("{kernel:?}");
436    match debug_repr.find(['{', '(']) {
437        Some(idx) => debug_repr[..idx].trim().to_string(),
438        None => debug_repr,
439    }
440}
441
442#[cfg(test)]
443#[allow(non_snake_case)]
444mod tests {
445    use super::*;
446    use crate::kernels::{Linear, RBF};
447    use crate::utils::robust_cholesky;
448
449    /// Deterministic standard-normal samples via a seeded 64-bit LCG +
450    /// Box-Muller transform, mirroring the RNG-free sampling style already
451    /// used elsewhere in this crate (see
452    /// `regression.rs::MultiOutputGaussianProcessRegressor::uniform`).
453    fn standard_normal_vec(n: usize, seed: u64) -> Array1<f64> {
454        let mut state = seed.max(1);
455        let mut next_uniform = || {
456            state = state
457                .wrapping_mul(6_364_136_223_846_793_005)
458                .wrapping_add(1_442_695_040_888_963_407);
459            ((state >> 11) as f64 / (1u64 << 53) as f64).clamp(1e-12, 1.0 - 1e-12)
460        };
461        Array1::from_iter((0..n).map(|_| {
462            let u1 = next_uniform();
463            let u2 = next_uniform();
464            (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
465        }))
466    }
467
468    /// Sample `y = L @ z` where `L` is the Cholesky factor of a known RBF
469    /// Gram matrix (plus a small noise variance on the diagonal), so the
470    /// data is drawn from a GP with a known, correct kernel family.
471    fn rbf_generated_data(
472        n: usize,
473        length_scale: f64,
474        noise_std: f64,
475        seed: u64,
476    ) -> (Array2<f64>, Array1<f64>) {
477        let X = Array2::from_shape_fn((n, 1), |(i, _)| i as f64 * 10.0 / n as f64);
478        let kernel = RBF::new(length_scale);
479        let mut K = kernel
480            .compute_kernel_matrix(&X, None)
481            .expect("kernel matrix computation should succeed");
482        for i in 0..n {
483            K[[i, i]] += noise_std * noise_std + 1e-8;
484        }
485        let L = robust_cholesky(&K).expect("Cholesky of RBF Gram matrix should succeed");
486        let z = standard_normal_vec(n, seed);
487        let y = L.dot(&z);
488        (X, y)
489    }
490
491    #[test]
492    fn test_select_kernel_aic_prefers_rbf_on_rbf_data() {
493        let (X, y) = rbf_generated_data(24, 1.0, 0.05, 7);
494        let candidates: Vec<Box<dyn Kernel>> =
495            vec![Box::new(RBF::new(1.0)), Box::new(Linear::new(1.0, 1.0))];
496        let result = select_best_kernel(candidates, &X, &y, SelectionCriterion::Aic)
497            .expect("kernel selection should succeed");
498        assert_eq!(result.best_kernel_name, "RBF");
499        assert_eq!(result.best_index, 0);
500    }
501
502    #[test]
503    fn test_select_kernel_bic_prefers_rbf_on_rbf_data() {
504        let (X, y) = rbf_generated_data(24, 1.0, 0.05, 11);
505        let candidates: Vec<Box<dyn Kernel>> =
506            vec![Box::new(RBF::new(1.0)), Box::new(Linear::new(1.0, 1.0))];
507        let result =
508            select_kernel_bic(candidates, &X, &y).expect("kernel selection should succeed");
509        assert_eq!(result.best_kernel_name, "RBF");
510        assert_eq!(result.all_scores.len(), 2);
511    }
512
513    #[test]
514    fn test_select_kernel_cv_runs_and_scores_all_candidates() {
515        let (X, y) = rbf_generated_data(24, 1.0, 0.05, 13);
516        let candidates: Vec<Box<dyn Kernel>> =
517            vec![Box::new(RBF::new(1.0)), Box::new(Linear::new(1.0, 1.0))];
518        let result =
519            select_kernel_cv(candidates, &X, &y).expect("CV kernel selection should succeed");
520        assert_eq!(result.all_scores.len(), 2);
521        assert!(result.best_score.is_finite());
522    }
523
524    #[test]
525    fn test_kernel_selector_builder() {
526        let candidates: Vec<Box<dyn Kernel>> = vec![Box::new(RBF::new(1.0))];
527        let selector = KernelSelector::new(candidates).criterion(SelectionCriterion::Bic);
528        assert_eq!(selector.config.criterion, SelectionCriterion::Bic);
529    }
530
531    #[test]
532    fn test_select_best_kernel_rejects_empty_candidates() {
533        let (X, y) = rbf_generated_data(10, 1.0, 0.05, 3);
534        let result = select_best_kernel(Vec::new(), &X, &y, SelectionCriterion::Aic);
535        assert!(result.is_err());
536    }
537
538    #[test]
539    fn test_selection_criterion_direction() {
540        assert!(!SelectionCriterion::Aic.higher_is_better());
541        assert!(!SelectionCriterion::Bic.higher_is_better());
542        assert!(SelectionCriterion::Cv.higher_is_better());
543        assert!(SelectionCriterion::LogMarginalLikelihood.higher_is_better());
544    }
545}
sklears_gaussian_process/kernel_selection.rs

sklears_gaussian_process/
kernel_selection.rs
