sklears_decomposition/
pca.rs

1//! Principal Component Analysis and dimensionality reduction utilities
2//!
3//! This module provides comprehensive PCA implementations using truncated SVD
4//! for exact principal components computation.
5
6use scirs2_core::ndarray::{Array1, Array2, Axis};
7use scirs2_linalg::svd;
8#[cfg(feature = "serde")]
9use serde::{Deserialize, Serialize};
10use sklears_core::{
11    error::{Result, SklearsError},
12    traits::{Estimator, Fit, Transform, Untrained},
13    types::Float,
14};
15
16/// Basic PCA configuration
17#[derive(Debug, Clone)]
18#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
19pub struct PcaConfig {
20    pub n_components: Option<usize>,
21    pub whiten: bool,
22    pub svd_solver: String,
23    pub tol: Float,
24    pub iterated_power: usize,
25    pub random_state: Option<u64>,
26}
27
28impl Default for PcaConfig {
29    fn default() -> Self {
30        Self {
31            n_components: None,
32            whiten: false,
33            svd_solver: "auto".to_string(),
34            tol: 0.0,
35            iterated_power: 5,
36            random_state: None,
37        }
38    }
39}
40
41/// Basic PCA implementation
42#[derive(Debug, Clone)]
43#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
44pub struct PCA<State = Untrained> {
45    config: PcaConfig,
46    state: std::marker::PhantomData<State>,
47}
48
49/// Trained PCA model
50#[derive(Debug, Clone)]
51#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
52pub struct PcaTrained {
53    config: PcaConfig,
54    /// Principal components (loadings)
55    pub components: Array2<Float>,
56    /// Explained variance
57    pub explained_variance: Array1<Float>,
58    /// Explained variance ratio
59    pub explained_variance_ratio: Array1<Float>,
60    /// Mean of training data
61    pub mean: Array1<Float>,
62    /// Number of components
63    pub n_components: usize,
64}
65
66impl PCA<Untrained> {
67    pub fn new(config: PcaConfig) -> Self {
68        Self {
69            config,
70            state: std::marker::PhantomData,
71        }
72    }
73
74    pub fn builder() -> PcaBuilder {
75        PcaBuilder::default()
76    }
77}
78
79impl Estimator for PCA<Untrained> {
80    type Config = PcaConfig;
81    type Error = SklearsError;
82    type Float = Float;
83
84    fn config(&self) -> &Self::Config {
85        &self.config
86    }
87}
88
89impl Fit<Array2<Float>, ()> for PCA<Untrained> {
90    type Fitted = PcaTrained;
91
92    fn fit(self, x: &Array2<Float>, _y: &()) -> Result<Self::Fitted> {
93        let (n_samples, n_features) = x.dim();
94
95        if n_samples == 0 || n_features == 0 {
96            return Err(SklearsError::InvalidInput(
97                "Input array must be non-empty".to_string(),
98            ));
99        }
100
101        // Determine number of components
102        let n_components = self
103            .config
104            .n_components
105            .unwrap_or(n_features.min(n_samples));
106
107        if n_components > n_features.min(n_samples) {
108            return Err(SklearsError::InvalidInput(
109                "n_components cannot be larger than min(n_samples, n_features)".to_string(),
110            ));
111        }
112
113        // Compute mean
114        let mean = x
115            .mean_axis(Axis(0))
116            .expect("array should have elements for mean computation");
117
118        // Center data: X_c = X - mean
119        let x_centered = x - &mean.clone().insert_axis(Axis(0));
120
121        // Full SVD of the centered data matrix: X_c = U S Vt
122        // Principal components are rows of Vt (right singular vectors).
123        // Singular values s relate to explained variance: var_i = s_i^2 / (n_samples - 1).
124        let (_u, singular_values, vt) = svd(&x_centered.view(), false, None)
125            .map_err(|e| SklearsError::NumericalError(format!("SVD failed in PCA.fit: {}", e)))?;
126
127        // Take the top n_components principal components (rows of Vt)
128        let components = vt
129            .slice(scirs2_core::ndarray::s![..n_components, ..])
130            .to_owned();
131
132        // Explained variance: s_i^2 / (n - 1), where s_i are singular values
133        let denom = if n_samples > 1 {
134            (n_samples - 1) as Float
135        } else {
136            1.0
137        };
138        let explained_variance: Array1<Float> = singular_values
139            .iter()
140            .take(n_components)
141            .map(|&s| s * s / denom)
142            .collect();
143
144        // Total variance = sum of all squared singular values / (n - 1)
145        let total_variance: Float = singular_values.iter().map(|&s| s * s / denom).sum();
146
147        let explained_variance_ratio = if total_variance > Float::EPSILON {
148            explained_variance.mapv(|v| v / total_variance)
149        } else {
150            Array1::zeros(n_components)
151        };
152
153        Ok(PcaTrained {
154            config: self.config,
155            components,
156            explained_variance,
157            explained_variance_ratio,
158            mean,
159            n_components,
160        })
161    }
162}
163
164impl Transform<Array2<Float>, Array2<Float>> for PcaTrained {
165    fn transform(&self, x: &Array2<Float>) -> Result<Array2<Float>> {
166        // Center data
167        let x_centered = x - &self.mean.clone().insert_axis(Axis(0));
168
169        // Project onto components
170        let x_transformed = x_centered.dot(&self.components.t());
171
172        Ok(x_transformed)
173    }
174}
175
176/// PCA builder pattern
177#[derive(Debug, Clone, Default)]
178pub struct PcaBuilder {
179    config: PcaConfig,
180}
181
182impl PcaBuilder {
183    pub fn n_components(mut self, n_components: Option<usize>) -> Self {
184        self.config.n_components = n_components;
185        self
186    }
187
188    pub fn whiten(mut self, whiten: bool) -> Self {
189        self.config.whiten = whiten;
190        self
191    }
192
193    pub fn random_state(mut self, random_state: Option<u64>) -> Self {
194        self.config.random_state = random_state;
195        self
196    }
197
198    pub fn build(self) -> PCA<Untrained> {
199        PCA::new(self.config)
200    }
201}
202
203// Placeholder types for compatibility with existing code
204pub type PcaProcessor = PCA<Untrained>;
205pub type PcaValidator = PCA<Untrained>;
206pub type PcaEstimator = PCA<Untrained>;
207pub type PcaTransformer = PcaTrained;
208pub type PcaAnalyzer = PcaTrained;
209pub type DimensionalityReducer = PcaTrained;
210pub type DecompositionEngine = PcaTrained;
211pub type StandardPcaAnalyzer = PcaTrained;
212pub type PrincipalComponentAnalysis = PcaTrained;
213pub type VarianceExplained = Array1<Float>;
214pub type ComponentExtractor = PcaTrained;
215pub type DimensionalityReduction = PcaTrained;
216pub type PcaProjection = Array2<Float>;
217pub type PcaInverseTransform = Array2<Float>;
218
219// Additional placeholder types that might be referenced elsewhere
220#[derive(Debug, Clone)]
221#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
222pub struct RobustPCA;
223
224#[derive(Debug, Clone)]
225#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
226pub struct SparsePCA;
227
228#[derive(Debug, Clone)]
229#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
230pub struct ProbabilisticPCA;
231
232// SVD solver related types
233#[derive(Debug, Clone)]
234#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
235pub struct SvdSolver;
236
237#[derive(Debug, Clone)]
238#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
239pub struct SvdSolverConfig;
240
241// More placeholder types for comprehensive compatibility
242pub type RobustPcaConfig = PcaConfig;
243pub type SparsePcaConfig = PcaConfig;
244pub type ProbabilisticPcaConfig = PcaConfig;
245pub type IncrementalPcaConfig = PcaConfig;
246
247#[cfg(test)]
248mod tests {
249    use super::*;
250    use scirs2_core::ndarray::{array, Array2};
251
252    /// Build a simple 2D dataset with clear principal directions.
253    fn two_feature_data() -> Array2<Float> {
254        array![
255            [1.0, 0.5],
256            [2.0, 1.0],
257            [3.0, 1.5],
258            [4.0, 2.0],
259            [5.0, 2.5],
260            [6.0, 3.0],
261        ]
262    }
263
264    #[test]
265    fn test_pca_output_shape() {
266        let x = two_feature_data();
267        let config = PcaConfig {
268            n_components: Some(1),
269            ..Default::default()
270        };
271        let trained = PCA::new(config).fit(&x, &()).expect("fit should succeed");
272        let transformed = trained.transform(&x).expect("transform should succeed");
273        assert_eq!(transformed.shape(), &[6, 1]);
274    }
275
276    #[test]
277    fn test_pca_components_shape() {
278        let x = two_feature_data();
279        let config = PcaConfig {
280            n_components: Some(2),
281            ..Default::default()
282        };
283        let trained = PCA::new(config).fit(&x, &()).expect("fit should succeed");
284        assert_eq!(trained.components.shape(), &[2, 2]);
285    }
286
287    #[test]
288    fn test_pca_explained_variance_ratio_sums_to_one() {
289        let x = two_feature_data();
290        let config = PcaConfig {
291            n_components: Some(2),
292            ..Default::default()
293        };
294        let trained = PCA::new(config).fit(&x, &()).expect("fit should succeed");
295        let ratio_sum: Float = trained.explained_variance_ratio.sum();
296        assert!((ratio_sum - 1.0).abs() < 1e-10, "ratio sum = {}", ratio_sum);
297    }
298
299    #[test]
300    fn test_pca_explained_variance_ratio_non_negative() {
301        let x = two_feature_data();
302        let config = PcaConfig {
303            n_components: Some(2),
304            ..Default::default()
305        };
306        let trained = PCA::new(config).fit(&x, &()).expect("fit should succeed");
307        for &r in trained.explained_variance_ratio.iter() {
308            assert!(r >= 0.0, "negative ratio {}", r);
309        }
310    }
311
312    #[test]
313    fn test_pca_explained_variance_descending() {
314        // For a dataset with clear variance, PC1 should explain more than PC2.
315        let x = two_feature_data();
316        let config = PcaConfig {
317            n_components: Some(2),
318            ..Default::default()
319        };
320        let trained = PCA::new(config).fit(&x, &()).expect("fit should succeed");
321        assert!(
322            trained.explained_variance[0] >= trained.explained_variance[1],
323            "PC1 variance {} < PC2 variance {}",
324            trained.explained_variance[0],
325            trained.explained_variance[1]
326        );
327    }
328
329    #[test]
330    fn test_pca_transform_finiteness() {
331        let x = two_feature_data();
332        let config = PcaConfig {
333            n_components: Some(1),
334            ..Default::default()
335        };
336        let trained = PCA::new(config).fit(&x, &()).expect("fit should succeed");
337        let t = trained.transform(&x).expect("transform should succeed");
338        for &v in t.iter() {
339            assert!(v.is_finite(), "non-finite value {}", v);
340        }
341    }
342
343    #[test]
344    fn test_pca_n_components_validation() {
345        let x = two_feature_data();
346        // n_components > min(n_samples, n_features)
347        let config = PcaConfig {
348            n_components: Some(100),
349            ..Default::default()
350        };
351        assert!(PCA::new(config).fit(&x, &()).is_err());
352    }
353
354    #[test]
355    fn test_pca_empty_input() {
356        let x: Array2<Float> = Array2::zeros((0, 3));
357        let config = PcaConfig::default();
358        assert!(PCA::new(config).fit(&x, &()).is_err());
359    }
360
361    #[test]
362    fn test_pca_first_component_captures_most_variance() {
363        // The dataset has a strong first principal component (y ≈ 0.5x)
364        let x = two_feature_data();
365        let config = PcaConfig {
366            n_components: Some(1),
367            ..Default::default()
368        };
369        let trained = PCA::new(config).fit(&x, &()).expect("fit should succeed");
370        // With 1 component, variance ratio must be close to 1 (data is nearly 1-D)
371        let ratio = trained.explained_variance_ratio[0];
372        assert!(ratio > 0.99, "1-PC ratio = {}", ratio);
373    }
374}
sklears_decomposition/pca.rs

sklears_decomposition/
pca.rs
