Add RNG to random projection structs

RNG defaults to Xoshiro256Plus if not provided by user.
Also added tests for minimum dimension using values from scikit-learn.

algorithms/linfa-reduction/Cargo.toml

@@ -45,6 +45,7 @@ rand = { version = "0.8", features = ["small_rng"] }
 linfa = { version = "0.7.0", path = "../.." }
 linfa-kernel = { version = "0.7.0", path = "../linfa-kernel" }
 sprs = "0.11.1"
+rand_xoshiro = "0.6.0"
 
 [dev-dependencies]
 ndarray-npy = { version = "0.8", default-features = false }

algorithms/linfa-reduction/examples/gaussian_projection.rs

@@ -6,14 +6,16 @@ use linfa_trees::{DecisionTree, SplitQuality};
 
 use mnist::{MnistBuilder, NormalizedMnist};
 use ndarray::{Array1, Array2};
-use rand::thread_rng;
+use rand::SeedableRng;
+use rand_xoshiro::Xoshiro256Plus;
 
 /// Train a Decision tree on the MNIST data set, with and without dimensionality reduction.
 fn main() -> Result<(), Box<dyn Error>> {
     // Parameters
     let train_sz = 10_000usize;
     let test_sz = 1_000usize;
     let reduced_dim = 100;
+    let rng = Xoshiro256Plus::seed_from_u64(42);
 
     let NormalizedMnist {
         trn_img,
@@ -54,10 +56,8 @@ fn main() -> Result<(), Box<dyn Error>> {
     println!("Training reduced model...");
     let start = Instant::now();
     // Compute the random projection and train the model on the reduced dataset.
-    let rng = thread_rng();
-    let proj = GaussianRandomProjection::<f32>::params()
+    let proj = GaussianRandomProjection::<f32>::params_with_rng(rng)
         .target_dim(reduced_dim)
-        .with_rng(rng)
         .fit(&train_dataset)?;
     let reduced_train_ds = proj.transform(&train_dataset);
     let reduced_test_data = proj.transform(&test_data);

algorithms/linfa-reduction/examples/sparse_projection.rs

@@ -6,13 +6,16 @@ use linfa_trees::{DecisionTree, SplitQuality};
 
 use mnist::{MnistBuilder, NormalizedMnist};
 use ndarray::{Array1, Array2};
+use rand::SeedableRng;
+use rand_xoshiro::Xoshiro256Plus;
 
 /// Train a Decision tree on the MNIST data set, with and without dimensionality reduction.
 fn main() -> Result<(), Box<dyn Error>> {
     // Parameters
     let train_sz = 10_000usize;
     let test_sz = 1_000usize;
     let reduced_dim = 100;
+    let rng = Xoshiro256Plus::seed_from_u64(42);
 
     let NormalizedMnist {
         trn_img,
@@ -53,7 +56,7 @@ fn main() -> Result<(), Box<dyn Error>> {
     println!("Training reduced model...");
     let start = Instant::now();
     // Compute the random projection and train the model on the reduced dataset.
-    let proj = SparseRandomProjection::<f32>::params()
+    let proj = SparseRandomProjection::<f32>::params_with_rng(rng)
         .target_dim(reduced_dim)
         .fit(&train_dataset)?;
     let reduced_train_ds = proj.transform(&train_dataset);

algorithms/linfa-reduction/src/random_projection/common.rs

@@ -5,7 +5,34 @@
 /// - [D. Achlioptas, JCSS](https://www.sciencedirect.com/science/article/pii/S0022000003000254)
 /// - [Li et al., SIGKDD'06](https://hastie.su.domains/Papers/Ping/KDD06_rp.pdf)
 pub(crate) fn johnson_lindenstrauss_min_dim(n_samples: usize, eps: f64) -> usize {
-    let log_samples = (n_samples as f64).log2();
-    let value = 4. * log_samples * (eps.powi(2) / 2. - eps.powi(3) / 3.);
-    value.ceil() as usize
+    let log_samples = (n_samples as f64).ln();
+    let value = 4. * log_samples / (eps.powi(2) / 2. - eps.powi(3) / 3.);
+    value as usize
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    /// Test against values computed by the scikit-learn implementation
+    /// of `johnson_lindenstrauss_min_dim`.
+    fn test_johnson_lindenstrauss() {
+        assert_eq!(johnson_lindenstrauss_min_dim(100, 0.05), 15244);
+        assert_eq!(johnson_lindenstrauss_min_dim(100, 0.1), 3947);
+        assert_eq!(johnson_lindenstrauss_min_dim(100, 0.2), 1062);
+        assert_eq!(johnson_lindenstrauss_min_dim(100, 0.5), 221);
+        assert_eq!(johnson_lindenstrauss_min_dim(1000, 0.05), 22867);
+        assert_eq!(johnson_lindenstrauss_min_dim(1000, 0.1), 5920);
+        assert_eq!(johnson_lindenstrauss_min_dim(1000, 0.2), 1594);
+        assert_eq!(johnson_lindenstrauss_min_dim(1000, 0.5), 331);
+        assert_eq!(johnson_lindenstrauss_min_dim(5000, 0.05), 28194);
+        assert_eq!(johnson_lindenstrauss_min_dim(5000, 0.1), 7300);
+        assert_eq!(johnson_lindenstrauss_min_dim(5000, 0.2), 1965);
+        assert_eq!(johnson_lindenstrauss_min_dim(5000, 0.5), 408);
+        assert_eq!(johnson_lindenstrauss_min_dim(10000, 0.05), 30489);
+        assert_eq!(johnson_lindenstrauss_min_dim(10000, 0.1), 7894);
+        assert_eq!(johnson_lindenstrauss_min_dim(10000, 0.2), 2125);
+        assert_eq!(johnson_lindenstrauss_min_dim(10000, 0.5), 442);
+    }
 }

algorithms/linfa-reduction/src/random_projection/gaussian/algorithms.rs

@@ -4,7 +4,8 @@ use ndarray_rand::{
     rand_distr::{Normal, StandardNormal},
     RandomExt,
 };
-use rand::{prelude::Distribution, rngs::SmallRng, Rng};
+use rand::{prelude::Distribution, Rng, SeedableRng};
+use rand_xoshiro::Xoshiro256Plus;
 
 use super::super::common::johnson_lindenstrauss_min_dim;
 use super::hyperparams::GaussianRandomProjectionParamsInner;
@@ -28,6 +29,7 @@ where
     fn fit(&self, dataset: &linfa::DatasetBase<Rec, T>) -> Result<Self::Object, ReductionError> {
         let n_samples = dataset.nsamples();
         let n_features = dataset.nfeatures();
+        let mut rng = self.rng.clone();
 
         let n_dims = match &self.params {
             GaussianRandomProjectionParamsInner::Dimension { target_dim } => *target_dim,
@@ -39,22 +41,31 @@ where
         let std_dev = F::cast(n_features).sqrt().recip();
         let gaussian = Normal::new(F::zero(), std_dev)?;
 
-        let proj = match self.rng.clone() {
-            Some(mut rng) => Array::random_using((n_features, n_dims), gaussian, &mut rng),
-            None => Array::random((n_features, n_dims), gaussian),
-        };
+        let proj = Array::random_using((n_features, n_dims), gaussian, &mut rng);
 
         Ok(GaussianRandomProjection { projection: proj })
     }
 }
 
 impl<F: Float> GaussianRandomProjection<F> {
+    /// Create new parameters for a [`GaussianRandomProjection`] with default value
+    /// `precision = 0.1` and a [`Xoshiro256Plus`] RNG.
+    pub fn params() -> GaussianRandomProjectionParams<Xoshiro256Plus> {
+        GaussianRandomProjectionParams(GaussianRandomProjectionValidParams {
+            params: GaussianRandomProjectionParamsInner::Precision { precision: 0.1 },
+            rng: Xoshiro256Plus::seed_from_u64(42),
+        })
+    }
+
     /// Create new parameters for a [`GaussianRandomProjection`] with default values
-    /// `precision = 0.1` and no custom [`Rng`] provided.
-    pub fn params() -> GaussianRandomProjectionParams<SmallRng> {
+    /// `precision = 0.1` and the provided [`Rng`].
+    pub fn params_with_rng<R>(rng: R) -> GaussianRandomProjectionParams<R>
+    where
+        R: Rng + Clone,
+    {
         GaussianRandomProjectionParams(GaussianRandomProjectionValidParams {
             params: GaussianRandomProjectionParamsInner::Precision { precision: 0.1 },
-            rng: None,
+            rng,
         })
     }
 }

algorithms/linfa-reduction/src/random_projection/gaussian/hyperparams.rs

@@ -43,12 +43,10 @@ impl<R: Rng + Clone> GaussianRandomProjectionParams<R> {
     }
 
     /// Specify the random number generator to use to generate the projection matrix.
-    ///
-    /// Optional: if no RNG is specified, uses the default RNG in [ndarray_rand::RandomExt].  
     pub fn with_rng<R2: Rng + Clone>(self, rng: R2) -> GaussianRandomProjectionParams<R2> {
         GaussianRandomProjectionParams(GaussianRandomProjectionValidParams {
             params: self.0.params,
-            rng: Some(rng),
+            rng,
         })
     }
 }
@@ -68,7 +66,7 @@ impl<R: Rng + Clone> GaussianRandomProjectionParams<R> {
 #[derive(Debug, Clone, PartialEq)]
 pub struct GaussianRandomProjectionValidParams<R: Rng + Clone> {
     pub(super) params: GaussianRandomProjectionParamsInner,
-    pub(super) rng: Option<R>,
+    pub(super) rng: R,
 }
 
 /// Internal data structure that either holds the dimension or the embedding,
@@ -107,8 +105,8 @@ impl<R: Rng + Clone> GaussianRandomProjectionValidParams<R> {
         self.params.eps()
     }
 
-    pub fn rng(&self) -> Option<&R> {
-        self.rng.as_ref()
+    pub fn rng(&self) -> &R {
+        &self.rng
     }
 }
 

algorithms/linfa-reduction/src/random_projection/mod.rs

@@ -39,23 +39,24 @@ pub use sparse::{
 #[cfg(test)]
 mod tests {
     use super::*;
-    use rand::rngs::SmallRng;
+
+    use rand_xoshiro::Xoshiro256Plus;
 
     #[test]
     fn autotraits_gaussian() {
         fn has_autotraits<T: Send + Sync + Sized + Unpin>() {}
         has_autotraits::<GaussianRandomProjection<f64>>();
         has_autotraits::<GaussianRandomProjection<f32>>();
-        has_autotraits::<GaussianRandomProjectionValidParams<SmallRng>>();
-        has_autotraits::<GaussianRandomProjectionParams<SmallRng>>();
+        has_autotraits::<GaussianRandomProjectionValidParams<Xoshiro256Plus>>();
+        has_autotraits::<GaussianRandomProjectionParams<Xoshiro256Plus>>();
     }
 
     #[test]
     fn autotraits_sparse() {
         fn has_autotraits<T: Send + Sync + Sized + Unpin>() {}
         has_autotraits::<SparseRandomProjection<f64>>();
         has_autotraits::<SparseRandomProjection<f32>>();
-        has_autotraits::<SparseRandomProjectionValidParams>();
-        has_autotraits::<SparseRandomProjectionParams>();
+        has_autotraits::<SparseRandomProjectionValidParams<Xoshiro256Plus>>();
+        has_autotraits::<SparseRandomProjectionParams<Xoshiro256Plus>>();
     }
 }

algorithms/linfa-reduction/src/random_projection/projection.rs

@@ -1,5 +1,5 @@
 /// Macro that implements [`linfa::traits::Transformer`]
-/// for [`GaussianRandomProjection`] and [`SparseRandomProjection`],
+/// for [`super::GaussianRandomProjection`] and [`super::SparseRandomProjection`],
 /// to avoid some code duplication.
 #[macro_export]
 macro_rules! impl_proj {

algorithms/linfa-reduction/src/random_projection/sparse/algorithms.rs

@@ -1,7 +1,9 @@
 use linfa::{prelude::Records, traits::Fit, Float};
 use ndarray::Ix2;
 use ndarray_rand::rand_distr::StandardNormal;
-use rand::{distributions::Bernoulli, prelude::Distribution, thread_rng, Rng};
+use rand::SeedableRng;
+use rand::{distributions::Bernoulli, prelude::Distribution, Rng};
+use rand_xoshiro::Xoshiro256Plus;
 use sprs::{CsMat, TriMat};
 
 use super::super::common::johnson_lindenstrauss_min_dim;
@@ -14,17 +16,19 @@ pub struct SparseRandomProjection<F: Float> {
     projection: CsMat<F>,
 }
 
-impl<F, Rec, T> Fit<Rec, T, ReductionError> for SparseRandomProjectionValidParams
+impl<F, Rec, T, R> Fit<Rec, T, ReductionError> for SparseRandomProjectionValidParams<R>
 where
     F: Float,
     Rec: Records<Elem = F>,
     StandardNormal: Distribution<F>,
+    R: Rng + Clone,
 {
     type Object = SparseRandomProjection<F>;
 
     fn fit(&self, dataset: &linfa::DatasetBase<Rec, T>) -> Result<Self::Object, ReductionError> {
         let n_samples = dataset.nsamples();
         let n_features = dataset.nfeatures();
+        let mut rng = self.rng.clone();
 
         let n_dims = match &self.params {
             SparseRandomProjectionParamsInner::Dimension { target_dim } => *target_dim,
@@ -36,7 +40,6 @@ where
         let scale = (n_features as f64).sqrt();
         let p = 1f64 / scale;
         let dist = SparseDistribution::new(F::cast(scale), p);
-        let mut rng = thread_rng();
 
         let (mut row_inds, mut col_inds, mut values) = (Vec::new(), Vec::new(), Vec::new());
         for row in 0..n_features {
@@ -90,10 +93,23 @@ impl<F: Float> Distribution<Option<F>> for SparseDistribution<F> {
 
 impl<F: Float> SparseRandomProjection<F> {
     /// Create new parameters for a [`SparseRandomProjection`] with default value
-    /// `precision = 0.1`.
-    pub fn params() -> SparseRandomProjectionParams {
+    /// `precision = 0.1` and a [`Xoshiro256Plus`] RNG.
+    pub fn params() -> SparseRandomProjectionParams<Xoshiro256Plus> {
         SparseRandomProjectionParams(SparseRandomProjectionValidParams {
             params: SparseRandomProjectionParamsInner::Precision { precision: 0.1 },
+            rng: Xoshiro256Plus::seed_from_u64(42),
+        })
+    }
+
+    /// Create new parameters for a [`SparseRandomProjection`] with default values
+    /// `precision = 0.1` and the provided [`Rng`].
+    pub fn params_with_rng<R>(rng: R) -> SparseRandomProjectionParams<R>
+    where
+        R: Rng + Clone,
+    {
+        SparseRandomProjectionParams(SparseRandomProjectionValidParams {
+            params: SparseRandomProjectionParamsInner::Precision { precision: 0.1 },
+            rng,
         })
     }
 }

algorithms/linfa-reduction/src/random_projection/sparse/hyperparams.rs

@@ -1,6 +1,7 @@
 use std::fmt::Debug;
 
 use linfa::ParamGuard;
+use rand::Rng;
 
 use crate::ReductionError;
 
@@ -13,9 +14,11 @@ use crate::ReductionError;
 /// However, this lemma makes a very conservative estimate of the required dimension,
 /// and does not leverage the structure of the data, therefore it is also possible
 /// to manually specify the dimension of the embedding.
-pub struct SparseRandomProjectionParams(pub(crate) SparseRandomProjectionValidParams);
+pub struct SparseRandomProjectionParams<R: Rng + Clone>(
+    pub(crate) SparseRandomProjectionValidParams<R>,
+);
 
-impl SparseRandomProjectionParams {
+impl<R: Rng + Clone> SparseRandomProjectionParams<R> {
     /// Set the dimension of output of the embedding.
     ///
     /// Setting the target dimension with this function
@@ -35,6 +38,14 @@ impl SparseRandomProjectionParams {
 
         self
     }
+
+    /// Specify the random number generator to use to generate the projection matrix.
+    pub fn with_rng<R2: Rng + Clone>(self, rng: R2) -> SparseRandomProjectionParams<R2> {
+        SparseRandomProjectionParams(SparseRandomProjectionValidParams {
+            params: self.0.params,
+            rng,
+        })
+    }
 }
 
 /// Sparse random projection hyperparameters
@@ -47,8 +58,9 @@ impl SparseRandomProjectionParams {
 /// and does not leverage the structure of the data, therefore it is also possible
 /// to manually specify the dimension of the embedding.
 #[derive(Debug, Clone, PartialEq)]
-pub struct SparseRandomProjectionValidParams {
+pub struct SparseRandomProjectionValidParams<R> {
     pub(super) params: SparseRandomProjectionParamsInner,
+    pub(super) rng: R,
 }
 
 /// Internal data structure that either holds the dimension or the embedding,
@@ -78,18 +90,22 @@ impl SparseRandomProjectionParamsInner {
     }
 }
 
-impl SparseRandomProjectionValidParams {
+impl<R: Rng + Clone> SparseRandomProjectionValidParams<R> {
     pub fn target_dim(&self) -> Option<usize> {
         self.params.target_dim()
     }
 
     pub fn precision(&self) -> Option<f64> {
         self.params.eps()
     }
+
+    pub fn rng(&self) -> &R {
+        &self.rng
+    }
 }
 
-impl ParamGuard for SparseRandomProjectionParams {
-    type Checked = SparseRandomProjectionValidParams;
+impl<R: Rng + Clone> ParamGuard for SparseRandomProjectionParams<R> {
+    type Checked = SparseRandomProjectionValidParams<R>;
     type Error = ReductionError;
 
     fn check_ref(&self) -> Result<&Self::Checked, Self::Error> {