Eliminate BLAS dependency for ALL crates (#207)

* Integrate ndarray-linalg-rs into linfa-clustering

* Restructure linfa-datasets to fix warnings

* Add support for linfa-ica

* Integrate into PLS (bugs pending)

* Integration with linfa-reduction complete

* Fix PLS CCA test

* Use invh instead of invc in elasticnet since the input might not be positive definite

* Integrate into linfa-preprocessing

* Integrate into linfa-linear

* Integrate into linfa-elasticnet

* Optimize gaussian mixture solve_triangular call

* Fix example Cargo invocations to not use blas features

* Remove BLAS from CI steps and add BLAS test job

* Enable most CI steps even in draft PRs

* Fix test-blas workflow name

* Eliminate BLAS from linfa-logistic

* Fix more examples invocations in READMEs

* Add non_exhaustive to modified error types so the blas feature is additive

* Change linfa-reduction to use truncated SVD and LOBPCG

* Switch to linfa-linalg

* Fix system BLAS to static BLAS in CI

* Centralized README docs on BLAS feature

Author: YuhanLiin

.github/workflows/benching.yml

@@ -6,10 +6,6 @@ jobs:
   testing:
     name: benching
     runs-on: ubuntu-18.04
-    if: github.event.pull_request.draft == false
-    container:
-      image: rustmath/mkl-rust:1.43.0
-      options: --security-opt seccomp=unconfined
 
     steps:
       - name: Checkout sources
@@ -26,4 +22,4 @@ jobs:
         uses: actions-rs/cargo@v1
         with:
           command: bench
-          args: iai --all --features intel-mkl-system
+          args: iai --all

.github/workflows/checking.yml

@@ -6,7 +6,6 @@ jobs:
   check:
     name: check-${{ matrix.toolchain }}-${{ matrix.os }}
     runs-on: ${{ matrix.os }}
-    if: github.event.pull_request.draft == false
     strategy:
       fail-fast: false
       matrix:

.github/workflows/codequality.yml

@@ -7,7 +7,6 @@ jobs:
   codequality:
     name: codequality
     runs-on: ubuntu-latest
-    if: github.event.pull_request.draft == false
     strategy:
       matrix:
         toolchain:
@@ -42,9 +41,6 @@ jobs:
     name: coverage
     runs-on: ubuntu-18.04
     if: github.event.pull_request.draft == false
-    container:
-      image: rustmath/mkl-rust:1.43.0
-      options: --security-opt seccomp=unconfined
 
     steps:
       - name: Checkout sources
@@ -74,7 +70,7 @@ jobs:
 
       - name: Generate code coverage
         run: |
-          cargo tarpaulin --verbose --features intel-mkl-system --timeout 120 --out Xml --all --release
+          cargo tarpaulin --verbose --timeout 120 --out Xml --all --release
       - name: Upload to codecov.io
         uses: codecov/codecov-action@v1
         with:

.github/workflows/testing.yml

@@ -6,7 +6,6 @@ jobs:
   testing:
     name: testing-${{ matrix.toolchain }}-${{ matrix.os }}
     runs-on: ${{ matrix.os }}
-    if: github.event.pull_request.draft == false
     strategy:
       fail-fast: false
       matrix:
@@ -28,11 +27,38 @@ jobs:
           toolchain: ${{ matrix.toolchain }}
           override: true
 
-      - name: Log active toolchain
-        run: rustup show
+      - name: Run cargo test
+        uses: actions-rs/cargo@v1
+        with:
+          command: test
+          args: --release --workspace
+
+  testing-blas:
+    name: testing-with-BLAS-${{ matrix.toolchain }}-${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        toolchain:
+          - 1.54.0
+          - stable
+        os:
+          - ubuntu-18.04
+          - windows-2019
+
+    steps:
+      - name: Checkout sources
+        uses: actions/checkout@v2
+
+      - name: Install toolchain
+        uses: actions-rs/toolchain@v1
+        with:
+          profile: minimal
+          toolchain: ${{ matrix.toolchain }}
+          override: true
 
-      - name: Run cargo test in release mode
+      - name: Run cargo test with BLAS enabled
         uses: actions-rs/cargo@v1
         with:
           command: test
-          args: --all --release --features intel-mkl-static
+          args: --release --workspace --features intel-mkl-static,linfa-clustering/blas,linfa-ica/blas,linfa-reduction/blas,linfa-linear/blas,linfa-preprocessing/blas,linfa-pls/blas,linfa-elasticnet/blas

CHANGELOG.md

@@ -2,11 +2,12 @@ Unreleased
 ========================
 
 Changes
-----------------------
-* remove `SeedableRng` trait bound from `KMeans` and `GaussianMixture`
+-----------
+ * remove `SeedableRng` trait bound from `KMeans` and `GaussianMixture`
+ * BLAS backend no longer required to build Linfa
 
 Breaking Changes
-----------------------
+-----------
  * parametrize `AsTargets` by the dimensionality of the targets and introduce `AsSingleTargets` and `AsMultiTargets`
  * 1D target arrays are no longer converted to 2D when constructing `Dataset`s
  * `Dataset` and `DatasetView` can now be parametrized by target dimensionality, with 2D being the default

README.md

@@ -53,19 +53,17 @@ If this strikes a chord with you, please take a look at the [roadmap](https://gi
 
 ## BLAS/Lapack backend
 
-At the moment you can choose between the following BLAS/LAPACK backends: `openblas`, `netblas` or `intel-mkl`
+Some algorithm crates need to use an external library for linear algebra routines. By default, we use a pure-Rust implementation. However, you can also choose an external BLAS/LAPACK backend library instead, by enabling the `blas` feature and a feature corresponding to your BLAS backend. Currently you can choose between the following BLAS/LAPACK backends: `openblas`, `netblas` or `intel-mkl`.
 
 |Backend  | Linux | Windows | macOS |
 |:--------|:-----:|:-------:|:-----:|
 |OpenBLAS |✔️      |-        |-      |
 |Netlib   |✔️      |-        |-      |
 |Intel MKL|✔️      |✔️        |✔️      |
 
-For example if you want to use the system IntelMKL library for the PCA example, then pass the corresponding feature:
-```
-cd linfa-reduction && cargo run --release --example pca --features linfa/intel-mkl-system
-```
-This selects the `intel-mkl` system library as BLAS/LAPACK backend. On the other hand if you want to compile the library and link it with the generated artifacts, pass `intel-mkl-static`.
+Each BLAS backend has two features available. The feature allows you to choose between linking the BLAS library in your system or statically building the library. For example, the features for the `intel-mkl` backend are `intel-mkl-static` and `intel-mkl-system`.
+
+An example set of Cargo flags for enabling the Intel MKL backend on an algorithm crate is `--features blas,linfa/intel-mkl-system`. Note that the BLAS backend features are defined on the `linfa` crate, and should only be specified for the final executable.
 
 # License
 Dual-licensed to be compatible with the Rust project.

algorithms/linfa-bayes/Cargo.toml

@@ -11,7 +11,7 @@ keywords = ["factorization", "machine-learning", "linfa", "unsupervised"]
 categories = ["algorithms", "mathematics", "science"]
 
 [dependencies]
-ndarray = { version = "0.15" , features = ["blas", "approx"]}
+ndarray = { version = "0.15" , features = ["approx"]}
 ndarray-stats = "0.5"
 thiserror = "1.0"
 

algorithms/linfa-clustering/Cargo.toml

@@ -18,6 +18,7 @@ categories = ["algorithms", "mathematics", "science"]
 
 [features]
 default = []
+blas = ["ndarray-linalg", "linfa/ndarray-linalg"]
 serde = ["serde_crate", "ndarray/serde", "linfa-nn/serde"]
 
 [dependencies.serde_crate]
@@ -29,15 +30,16 @@ features = ["std", "derive"]
 
 [dependencies]
 ndarray = { version = "0.15", features = ["rayon", "approx"]}
-ndarray-linalg = "0.14"
+linfa-linalg = { version = "0.1", default-features = false }
+ndarray-linalg = { version = "0.14", optional = true }
 ndarray-rand = "0.14"
 ndarray-stats = "0.5"
 num-traits = "0.2"
 rand_xoshiro = "0.6"
 space = "0.12"
 thiserror = "1.0"
 partitions = "0.2.4"
-linfa = { version = "0.5.0", path = "../..", features = ["ndarray-linalg"] }
+linfa = { version = "0.5.0", path = "../.." }
 linfa-nn = { version = "0.5.0", path = "../linfa-nn" }
 noisy_float = "0.2.0"
 

algorithms/linfa-clustering/README.md

@@ -23,6 +23,10 @@ Implementation choices, algorithmic details and a tutorial can be found
 
 **WARNING:** Currently the Approximated DBSCAN implementation is slower than the normal DBSCAN implementation. Therefore DBSCAN should always be used over Approximated DBSCAN.
 
+## BLAS/Lapack backend
+
+See [this section](../../README.md#blaslapack-backend) to enable an external BLAS/LAPACK backend.
+
 ## License
 Dual-licensed to be compatible with the Rust project.
 

algorithms/linfa-clustering/src/gaussian_mixture/algorithm.rs

@@ -3,13 +3,14 @@ use crate::gaussian_mixture::hyperparams::{
     GmmCovarType, GmmInitMethod, GmmParams, GmmValidParams,
 };
 use crate::k_means::KMeans;
-use linfa::{
-    dataset::{WithLapack, WithoutLapack},
-    prelude::*,
-    DatasetBase, Float,
-};
+#[cfg(feature = "blas")]
+use linfa::dataset::{WithLapack, WithoutLapack};
+use linfa::{prelude::*, DatasetBase, Float};
+#[cfg(not(feature = "blas"))]
+use linfa_linalg::{cholesky::*, triangular::*};
 use ndarray::{s, Array, Array1, Array2, Array3, ArrayBase, Axis, Data, Ix2, Ix3, Zip};
-use ndarray_linalg::{cholesky::*, triangular::*, Lapack, Scalar};
+#[cfg(feature = "blas")]
+use ndarray_linalg::{cholesky::*, triangular::*};
 use ndarray_rand::rand::Rng;
 use ndarray_rand::rand_distr::Uniform;
 use ndarray_rand::RandomExt;
@@ -264,10 +265,18 @@ impl<F: Float> GaussianMixtureModel<F> {
         let n_features = covariances.shape()[1];
         let mut precisions_chol = Array::zeros((n_clusters, n_features, n_features));
         for (k, covariance) in covariances.outer_iter().enumerate() {
-            let decomp = covariance.with_lapack().cholesky(UPLO::Lower)?;
-            let sol = decomp
-                .solve_triangular(UPLO::Lower, Diag::NonUnit, &Array::eye(n_features))?
-                .without_lapack();
+            #[cfg(feature = "blas")]
+            let sol = {
+                let decomp = covariance.with_lapack().cholesky(UPLO::Lower)?;
+                decomp
+                    .solve_triangular_into(UPLO::Lower, Diag::NonUnit, Array::eye(n_features))?
+                    .without_lapack()
+            };
+            #[cfg(not(feature = "blas"))]
+            let sol = {
+                let decomp = covariance.cholesky()?;
+                decomp.solve_triangular_into(Array::eye(n_features), UPLO::Lower)?
+            };
 
             precisions_chol.slice_mut(s![k, .., ..]).assign(&sol.t());
         }
@@ -461,7 +470,7 @@ impl<F: Float, R: Rng + Clone, D: Data<Elem = F>, T> Fit<ArrayBase<D, Ix2>, T, G
     }
 }
 
-impl<F: Float + Lapack + Scalar, D: Data<Elem = F>> PredictInplace<ArrayBase<D, Ix2>, Array1<usize>>
+impl<F: Float, D: Data<Elem = F>> PredictInplace<ArrayBase<D, Ix2>, Array1<usize>>
     for GaussianMixtureModel<F>
 {
     fn predict_inplace(&self, observations: &ArrayBase<D, Ix2>, targets: &mut Array1<usize>) {
@@ -473,7 +482,7 @@ impl<F: Float + Lapack + Scalar, D: Data<Elem = F>> PredictInplace<ArrayBase<D,
 
         let (_, log_resp) = self.estimate_log_prob_resp(observations);
         *targets = log_resp
-            .mapv(Scalar::exp)
+            .mapv(F::exp)
             .map_axis(Axis(1), |row| row.argmax().unwrap());
     }
 
@@ -486,10 +495,18 @@ impl<F: Float + Lapack + Scalar, D: Data<Elem = F>> PredictInplace<ArrayBase<D,
 mod tests {
     use super::*;
     use approx::{abs_diff_eq, assert_abs_diff_eq};
+    #[cfg(feature = "blas")]
     use lax::error::Error;
     use linfa_datasets::generate;
     use ndarray::{array, concatenate, ArrayView1, ArrayView2, Axis};
+
+    #[cfg(not(feature = "blas"))]
+    use linfa_linalg::LinalgError;
+    #[cfg(not(feature = "blas"))]
+    use linfa_linalg::Result as LAResult;
+    #[cfg(feature = "blas")]
     use ndarray_linalg::error::LinalgError;
+    #[cfg(feature = "blas")]
     use ndarray_linalg::error::Result as LAResult;
     use ndarray_rand::rand::prelude::ThreadRng;
     use ndarray_rand::rand::SeedableRng;
@@ -514,7 +531,10 @@ mod tests {
     }
     impl MultivariateNormal {
         pub fn new(mean: &ArrayView1<f64>, covariance: &ArrayView2<f64>) -> LAResult<Self> {
+            #[cfg(feature = "blas")]
             let lower = covariance.cholesky(UPLO::Lower)?;
+            #[cfg(not(feature = "blas"))]
+            let lower = covariance.cholesky()?;
             Ok(MultivariateNormal {
                 mean: mean.to_owned(),
                 covariance: covariance.to_owned(),
@@ -603,16 +623,18 @@ mod tests {
             .with_rng(rng.clone())
             .fit(&dataset);
 
-        assert!(
-            match gmm.expect_err("should generate an error with reg_covar being nul") {
-                GmmError::LinalgError(e) => match e {
-                    LinalgError::Lapack(Error::LapackComputationalFailure { return_code: 2 }) =>
-                        true,
-                    _ => panic!("should be a lapack error 2"),
-                },
-                _ => panic!("should be a linear algebra error"),
+        match gmm.expect_err("should generate an error with reg_covar being nul") {
+            GmmError::LinalgError(e) => {
+                #[cfg(feature = "blas")]
+                assert!(matches!(
+                    e,
+                    LinalgError::Lapack(Error::LapackComputationalFailure { return_code: 2 })
+                ));
+                #[cfg(not(feature = "blas"))]
+                assert!(matches!(e, LinalgError::NotPositiveDefinite));
             }
-        );
+            e => panic!("should be a linear algebra error: {:?}", e),
+        }
         // Test it passes when default value is used
         assert!(GaussianMixtureModel::params(3)
             .with_rng(rng)
@@ -632,16 +654,18 @@ mod tests {
             .reg_covariance(0.)
             .fit(&dataset);
 
-        assert!(
-            match gmm.expect_err("should generate an error with reg_covar being nul") {
-                GmmError::LinalgError(e) => match e {
-                    LinalgError::Lapack(Error::LapackComputationalFailure { return_code: 1 }) =>
-                        true,
-                    _ => panic!("should be a lapack error 1"),
-                },
-                _ => panic!("should be a linear algebra error"),
+        #[cfg(feature = "blas")]
+        match gmm.expect_err("should generate an error with reg_covar being nul") {
+            GmmError::LinalgError(e) => {
+                assert!(matches!(
+                    e,
+                    LinalgError::Lapack(Error::LapackComputationalFailure { return_code: 1 })
+                ));
             }
-        );
+            e => panic!("should be a linear algebra error: {:?}", e),
+        }
+        #[cfg(not(feature = "blas"))]
+        gmm.expect_err("should generate an error with reg_covar being nul");
 
         // Test it passes when default value is used
         assert!(GaussianMixtureModel::params(1).fit(&dataset).is_ok());