diff --git a/Cargo.toml b/Cargo.toml
index d284d10da..61b8475ec 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -38,11 +38,11 @@ serde = ["serde_crate", "ndarray/serde"]
 [dependencies]
 num-traits = "0.2"
 rand = { version = "0.8", features = ["small_rng"] }
-approx = "0.4"
+approx = "0.5"
 
-ndarray = { version = "0.15", features = ["approx"] }
-ndarray-linalg = { version = "0.16", optional = true }
-sprs = { version = "=0.11.1", default-features = false }
+ndarray = { version = "0.16", features = ["approx"] }
+ndarray-linalg = { version = "0.17", optional = true }
+sprs = { version = "0.11", default-features = false }
 
 thiserror = "1.0"
 
@@ -56,7 +56,7 @@ default-features = false
 features = ["std", "derive"]
 
 [dev-dependencies]
-ndarray-rand = "0.14"
+ndarray-rand = "0.15"
 linfa-datasets = { path = "datasets", features = [
     "winequality",
     "iris",
diff --git a/algorithms/linfa-bayes/Cargo.toml b/algorithms/linfa-bayes/Cargo.toml
index 5030a7e45..ae7e48654 100644
--- a/algorithms/linfa-bayes/Cargo.toml
+++ b/algorithms/linfa-bayes/Cargo.toml
@@ -21,12 +21,14 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = { version = "0.15" , features = ["approx"]}
-ndarray-stats = "0.5"
+ndarray = { version = "0.16", features = ["approx"] }
+ndarray-stats = "0.6"
 thiserror = "1.0"
 
 linfa = { version = "0.7.1", path = "../.." }
 
 [dev-dependencies]
-approx = "0.4"
-linfa-datasets = { version = "0.7.1", path = "../../datasets", features = ["winequality"] }
+approx = "0.5"
+linfa-datasets = { version = "0.7.1", path = "../../datasets", features = [
+    "winequality",
+] }
diff --git a/algorithms/linfa-clustering/Cargo.toml b/algorithms/linfa-clustering/Cargo.toml
index 11ba07da9..ebb211ce9 100644
--- a/algorithms/linfa-clustering/Cargo.toml
+++ b/algorithms/linfa-clustering/Cargo.toml
@@ -35,11 +35,11 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = { version = "0.15", features = ["rayon", "approx"] }
-linfa-linalg = { version = "0.1", default-features = false }
-ndarray-linalg = { version = "0.16", optional = true }
-ndarray-rand = "0.14"
-ndarray-stats = "0.5"
+ndarray = { version = "0.16", features = ["rayon", "approx"] }
+linfa-linalg = { version = "0.2", default-features = false }
+ndarray-linalg = { version = "0.17", optional = true }
+ndarray-rand = "0.15"
+ndarray-stats = "0.6"
 num-traits = "0.2"
 rand_xoshiro = "0.6"
 space = "0.12"
@@ -50,14 +50,14 @@ linfa-nn = { version = "0.7.1", path = "../linfa-nn" }
 noisy_float = "0.2.0"
 
 [dev-dependencies]
-ndarray-npy = { version = "0.8", default-features = false }
+ndarray-npy = { version = "0.9", default-features = false }
 linfa-datasets = { version = "0.7.1", path = "../../datasets", features = [
     "generate",
 ] }
 criterion = "0.4.0"
 serde_json = "1"
-approx = "0.4"
-lax = "0.15.0"
+approx = "0.5"
+lax = "0.17.0"
 linfa = { version = "0.7.1", path = "../..", features = ["benchmarks"] }
 
 [[bench]]
diff --git a/algorithms/linfa-clustering/src/gaussian_mixture/algorithm.rs b/algorithms/linfa-clustering/src/gaussian_mixture/algorithm.rs
index 5d3da7524..54ddf222e 100644
--- a/algorithms/linfa-clustering/src/gaussian_mixture/algorithm.rs
+++ b/algorithms/linfa-clustering/src/gaussian_mixture/algorithm.rs
@@ -390,7 +390,7 @@ impl<F: Float> GaussianMixtureModel<F> {
         let n_clusters = matrix_chol.shape()[0];
         let log_diags = &matrix_chol
             .to_owned()
-            .into_shape((n_clusters, n_features * n_features))
+            .into_shape_with_order((n_clusters, n_features * n_features))
             .unwrap()
             .slice(s![.., ..; n_features+1])
             .to_owned()
diff --git a/algorithms/linfa-elasticnet/Cargo.toml b/algorithms/linfa-elasticnet/Cargo.toml
index 2148c0d11..aedf9fa7b 100644
--- a/algorithms/linfa-elasticnet/Cargo.toml
+++ b/algorithms/linfa-elasticnet/Cargo.toml
@@ -29,12 +29,12 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = { version = "0.15", features = ["approx"] }
-linfa-linalg = { version = "0.1", default-features = false }
-ndarray-linalg = { version = "0.16", optional = true }
+ndarray = { version = "0.16", features = ["approx"] }
+linfa-linalg = { version = "0.2", default-features = false }
+ndarray-linalg = { version = "0.17", optional = true }
 
 num-traits = "0.2"
-approx = "0.4"
+approx = "0.5"
 thiserror = "1.0"
 
 linfa = { version = "0.7.1", path = "../.." }
@@ -44,5 +44,5 @@ linfa-datasets = { version = "0.7.1", path = "../../datasets", features = [
     "diabetes",
     "linnerud",
 ] }
-ndarray-rand = "0.14"
+ndarray-rand = "0.15"
 rand_xoshiro = "0.6"
diff --git a/algorithms/linfa-ftrl/Cargo.toml b/algorithms/linfa-ftrl/Cargo.toml
index 716fb6476..daf95440c 100644
--- a/algorithms/linfa-ftrl/Cargo.toml
+++ b/algorithms/linfa-ftrl/Cargo.toml
@@ -24,7 +24,7 @@ version = "1.0"
 features = ["derive"]
 
 [dependencies]
-ndarray = { version = "0.15.4", features = ["serde"] }
+ndarray = { version = "0.16", features = ["serde"] }
 ndarray-rand = "0.14.0"
 argmin = { version = "0.9.0", default-features = false }
 argmin-math = { version = "0.3", features = ["ndarray_v0_15-nolinalg"] }
@@ -36,7 +36,7 @@ linfa = { version = "0.7.1", path = "../.." }
 
 [dev-dependencies]
 criterion = "0.4.0"
-approx = "0.4"
+approx = "0.5"
 linfa-datasets = { version = "0.7.1", path = "../../datasets", features = [
     "winequality",
 ] }
diff --git a/algorithms/linfa-hierarchical/Cargo.toml b/algorithms/linfa-hierarchical/Cargo.toml
index f8aafb7cd..bf658cd20 100644
--- a/algorithms/linfa-hierarchical/Cargo.toml
+++ b/algorithms/linfa-hierarchical/Cargo.toml
@@ -10,11 +10,17 @@ license = "MIT OR Apache-2.0"
 repository = "https://github.com/rust-ml/linfa"
 readme = "README.md"
 
-keywords = ["hierachical", "agglomerative", "clustering", "machine-learning", "linfa"]
+keywords = [
+    "hierachical",
+    "agglomerative",
+    "clustering",
+    "machine-learning",
+    "linfa",
+]
 categories = ["algorithms", "mathematics", "science"]
 
 [dependencies]
-ndarray = { version = "0.15" }
+ndarray = { version = "0.16" }
 kodama = "0.2"
 thiserror = "1.0.25"
 
@@ -23,5 +29,7 @@ linfa-kernel = { version = "0.7.1", path = "../linfa-kernel" }
 
 [dev-dependencies]
 rand = "0.8"
-ndarray-rand = "0.14"
-linfa-datasets = { version = "0.7.1", path = "../../datasets", features = ["iris"] }
+ndarray-rand = "0.15"
+linfa-datasets = { version = "0.7.1", path = "../../datasets", features = [
+    "iris",
+] }
diff --git a/algorithms/linfa-ica/Cargo.toml b/algorithms/linfa-ica/Cargo.toml
index 7dcdb7561..d33da9be6 100644
--- a/algorithms/linfa-ica/Cargo.toml
+++ b/algorithms/linfa-ica/Cargo.toml
@@ -25,11 +25,11 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = { version = "0.15" }
-linfa-linalg = { version = "0.1", default-features = false }
-ndarray-linalg = { version = "0.16", optional = true }
-ndarray-rand = "0.14"
-ndarray-stats = "0.5"
+ndarray = { version = "0.16" }
+linfa-linalg = { version = "0.2", default-features = false }
+ndarray-linalg = { version = "0.17", optional = true }
+ndarray-rand = "0.15"
+ndarray-stats = "0.6"
 num-traits = "0.2"
 rand_xoshiro = "0.6"
 thiserror = "1.0"
@@ -37,7 +37,7 @@ thiserror = "1.0"
 linfa = { version = "0.7.1", path = "../.." }
 
 [dev-dependencies]
-ndarray-npy = { version = "0.8", default-features = false }
+ndarray-npy = { version = "0.9", default-features = false }
 paste = "1.0"
 criterion = "0.4.0"
 linfa = { version = "0.7.1", path = "../..", features = ["benchmarks"] }
diff --git a/algorithms/linfa-kernel/Cargo.toml b/algorithms/linfa-kernel/Cargo.toml
index e6792a91c..87da7486c 100644
--- a/algorithms/linfa-kernel/Cargo.toml
+++ b/algorithms/linfa-kernel/Cargo.toml
@@ -24,9 +24,9 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = "0.15"
+ndarray = "0.16"
 num-traits = "0.2"
-sprs = { version = "=0.11.1", default-features = false }
+sprs = { version = "0.11", default-features = false }
 
 linfa = { version = "0.7.1", path = "../.." }
 linfa-nn = { version = "0.7.1", path = "../linfa-nn" }
diff --git a/algorithms/linfa-linear/Cargo.toml b/algorithms/linfa-linear/Cargo.toml
index 0c9f7ec6d..0421b49d9 100644
--- a/algorithms/linfa-linear/Cargo.toml
+++ b/algorithms/linfa-linear/Cargo.toml
@@ -29,9 +29,9 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = { version = "0.15", features = ["approx"] }
-linfa-linalg = { version = "0.1", default-features = false }
-ndarray-linalg = { version = "0.16", optional = true }
+ndarray = { version = "0.16", features = ["approx"] }
+linfa-linalg = { version = "0.2", default-features = false }
+ndarray-linalg = { version = "0.17", optional = true }
 num-traits = "0.2"
 argmin = { version = "0.9.0", default-features = false }
 argmin-math = { version = "0.3", features = ["ndarray_v0_15-nolinalg"] }
@@ -43,7 +43,7 @@ linfa = { version = "0.7.1", path = "../.." }
 linfa-datasets = { version = "0.7.1", path = "../../datasets", features = [
     "diabetes",
 ] }
-approx = "0.4"
+approx = "0.5"
 criterion = "0.4.0"
 statrs = "0.16.0"
 linfa = { version = "0.7.1", path = "../..", features = ["benchmarks"] }
diff --git a/algorithms/linfa-logistic/Cargo.toml b/algorithms/linfa-logistic/Cargo.toml
index a0d25223b..0b5130225 100644
--- a/algorithms/linfa-logistic/Cargo.toml
+++ b/algorithms/linfa-logistic/Cargo.toml
@@ -23,7 +23,7 @@ optional = true
 version = "1.0"
 
 [dependencies]
-ndarray = { version = "0.15", features = ["approx"] }
+ndarray = { version = "0.16", features = ["approx"] }
 ndarray-stats = "0.5.0"
 num-traits = "0.2"
 argmin = { version = "0.9.0", default-features = false }
@@ -34,7 +34,7 @@ thiserror = "1.0"
 linfa = { version = "0.7.1", path = "../.." }
 
 [dev-dependencies]
-approx = "0.4"
+approx = "0.5"
 linfa-datasets = { version = "0.7.1", path = "../../datasets", features = [
     "winequality",
 ] }
diff --git a/algorithms/linfa-logistic/src/lib.rs b/algorithms/linfa-logistic/src/lib.rs
index 833dbab0b..49f068aca 100644
--- a/algorithms/linfa-logistic/src/lib.rs
+++ b/algorithms/linfa-logistic/src/lib.rs
@@ -476,9 +476,9 @@ fn logistic_loss<F: Float, A: Data<Elem = F>>(
 ) -> F {
     let n_features = x.shape()[1];
     let (params, intercept) = convert_params(n_features, w);
-    let yz = x.dot(&params.into_shape((params.len(), 1)).unwrap()) + intercept;
+    let yz = x.dot(&params.into_shape_with_order((params.len(), 1)).unwrap()) + intercept;
     let len = yz.len();
-    let mut yz = yz.into_shape(len).unwrap() * y;
+    let mut yz = yz.into_shape_with_order(len).unwrap() * y;
     yz.mapv_inplace(log_logistic);
     -yz.sum() + F::cast(0.5) * alpha * params.dot(&params)
 }
@@ -492,9 +492,9 @@ fn logistic_grad<F: Float, A: Data<Elem = F>>(
 ) -> Array1<F> {
     let n_features = x.shape()[1];
     let (params, intercept) = convert_params(n_features, w);
-    let yz = x.dot(&params.into_shape((params.len(), 1)).unwrap()) + intercept;
+    let yz = x.dot(&params.into_shape_with_order((params.len(), 1)).unwrap()) + intercept;
     let len = yz.len();
-    let mut yz = yz.into_shape(len).unwrap() * y;
+    let mut yz = yz.into_shape_with_order(len).unwrap() * y;
     yz.mapv_inplace(logistic);
     yz -= F::one();
     yz *= y;
@@ -523,7 +523,10 @@ fn multi_logistic_prob_params<'a, F: Float, A: Data<Elem = F>>(
     let h = x.dot(&params) + intercept;
     // This computes `H - log(sum(exp(H)))`, which is equal to
     // `log(softmax(H)) = log(exp(H) / sum(exp(H)))`
-    let log_prob = &h - log_sum_exp(&h, Axis(1)).into_shape((h.nrows(), 1)).unwrap();
+    let log_prob = &h
+        - log_sum_exp(&h, Axis(1))
+            .into_shape_with_order((h.nrows(), 1))
+            .unwrap();
     (log_prob, params)
 }
 
diff --git a/algorithms/linfa-nn/Cargo.toml b/algorithms/linfa-nn/Cargo.toml
index f599668b3..9a8c3e699 100644
--- a/algorithms/linfa-nn/Cargo.toml
+++ b/algorithms/linfa-nn/Cargo.toml
@@ -24,8 +24,8 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = { version = "0.15", features = ["approx"]}
-ndarray-stats = "0.5"
+ndarray = { version = "0.16", features = ["approx"] }
+ndarray-stats = "0.6"
 num-traits = "0.2.0"
 noisy_float = "0.2.0"
 order-stat = "0.1.3"
@@ -36,10 +36,10 @@ kdtree = "0.6.0"
 linfa = { version = "0.7.1", path = "../.." }
 
 [dev-dependencies]
-approx = "0.4"
+approx = "0.5"
 criterion = "0.4.0"
 rand_xoshiro = "0.6"
-ndarray-rand = "0.14"
+ndarray-rand = "0.15"
 linfa = { version = "0.7.1", path = "../..", features = ["benchmarks"] }
 
 [[bench]]
diff --git a/algorithms/linfa-pls/Cargo.toml b/algorithms/linfa-pls/Cargo.toml
index 13002b8ad..368a825a3 100644
--- a/algorithms/linfa-pls/Cargo.toml
+++ b/algorithms/linfa-pls/Cargo.toml
@@ -25,11 +25,11 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = { version = "0.15" }
-linfa-linalg = { version = "0.1", default-features = false }
-ndarray-linalg = { version = "0.16", optional = true }
-ndarray-stats = "0.5"
-ndarray-rand = "0.14"
+ndarray = { version = "0.16" }
+linfa-linalg = { version = "0.2", default-features = false }
+ndarray-linalg = { version = "0.17", optional = true }
+ndarray-stats = "0.6"
+ndarray-rand = "0.15"
 num-traits = "0.2"
 paste = "1.0"
 thiserror = "1.0"
@@ -40,7 +40,7 @@ linfa = { version = "0.7.1", path = "../..", features = ["benchmarks"] }
 linfa-datasets = { version = "0.7.1", path = "../../datasets", features = [
     "linnerud",
 ] }
-approx = "0.4"
+approx = "0.5"
 rand_xoshiro = "0.6"
 criterion = "0.4.0"
 statrs = "0.16.0"
diff --git a/algorithms/linfa-preprocessing/Cargo.toml b/algorithms/linfa-preprocessing/Cargo.toml
index f6878e081..638bd29eb 100644
--- a/algorithms/linfa-preprocessing/Cargo.toml
+++ b/algorithms/linfa-preprocessing/Cargo.toml
@@ -19,17 +19,17 @@ serde = ["serde_crate", "ndarray/serde", "serde_regex"]
 
 [dependencies]
 linfa = { version = "0.7.1", path = "../.." }
-ndarray = { version = "0.15", features = ["approx"] }
-ndarray-linalg = { version = "0.16", optional = true }
-linfa-linalg = { version = "0.1", default-features = false }
-ndarray-stats = "0.5"
+ndarray = { version = "0.16", features = ["approx"] }
+ndarray-linalg = { version = "0.17", optional = true }
+linfa-linalg = { version = "0.2", default-features = false }
+ndarray-stats = "0.6"
 thiserror = "1.0"
-approx = { version = "0.4" }
-ndarray-rand = { version = "0.14" }
+approx = "0.5"
+ndarray-rand = "0.15"
 unicode-normalization = "0.1.8"
 regex = "1.4.5"
 encoding = "0.2"
-sprs = { version = "=0.11.1", default-features = false }
+sprs = { version = "0.11", default-features = false }
 
 serde_regex = { version = "1.1", optional = true }
 
diff --git a/algorithms/linfa-reduction/Cargo.toml b/algorithms/linfa-reduction/Cargo.toml
index c1de0c2f9..c658a957c 100644
--- a/algorithms/linfa-reduction/Cargo.toml
+++ b/algorithms/linfa-reduction/Cargo.toml
@@ -34,25 +34,25 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = { version = "0.15", features = ["approx"] }
-linfa-linalg = { version = "0.1" }
-ndarray-linalg = { version = "0.16", optional = true }
-ndarray-rand = "0.14"
+ndarray = { version = "0.16", features = ["approx"] }
+linfa-linalg = { version = "0.2" }
+ndarray-linalg = { version = "0.17", optional = true }
+ndarray-rand = "0.15"
 num-traits = "0.2"
 thiserror = "1.0"
 rand = { version = "0.8", features = ["small_rng"] }
 
 linfa = { version = "0.7.1", path = "../.." }
 linfa-kernel = { version = "0.7.1", path = "../linfa-kernel" }
-sprs = "=0.11.1"
+sprs = { version = "0.11" }
 rand_xoshiro = "0.6.0"
 
 [dev-dependencies]
-ndarray-npy = { version = "0.8", default-features = false }
+ndarray-npy = { version = "0.9", default-features = false }
 linfa-datasets = { version = "0.7.1", path = "../../datasets", features = [
     "iris",
     "generate",
 ] }
-approx = { version = "0.4" }
+approx = "0.5"
 mnist = { version = "0.6.0", features = ["download"] }
 linfa-trees = { version = "0.7.1", path = "../linfa-trees" }
diff --git a/algorithms/linfa-svm/Cargo.toml b/algorithms/linfa-svm/Cargo.toml
index c99784495..d220cf9f0 100644
--- a/algorithms/linfa-svm/Cargo.toml
+++ b/algorithms/linfa-svm/Cargo.toml
@@ -24,8 +24,8 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = { version = "0.15" }
-ndarray-rand = "0.14"
+ndarray = { version = "0.16" }
+ndarray-rand = "0.15"
 num-traits = "0.2"
 thiserror = "1.0"
 
@@ -38,4 +38,4 @@ linfa-datasets = { version = "0.7.1", path = "../../datasets", features = [
     "diabetes",
 ] }
 rand_xoshiro = "0.6"
-approx = "0.4"
+approx = "0.5"
diff --git a/algorithms/linfa-svm/examples/noisy_sin_svr.rs b/algorithms/linfa-svm/examples/noisy_sin_svr.rs
index ca8cc64fb..595d12aeb 100644
--- a/algorithms/linfa-svm/examples/noisy_sin_svr.rs
+++ b/algorithms/linfa-svm/examples/noisy_sin_svr.rs
@@ -23,7 +23,7 @@ fn main() -> Result<()> {
         .filter(|(i, _)| i % 5 == 0)
         .for_each(|(_, y)| *y = 3. * (0.5 - rng.gen::<f64>()));
 
-    let x = x.into_shape((40, 1)).unwrap();
+    let x = x.into_shape_with_order((40, 1)).unwrap();
     let dataset = DatasetBase::new(x, y);
     let model = Svm::params()
         .c_svr(100., Some(0.1))
diff --git a/algorithms/linfa-svm/src/regression.rs b/algorithms/linfa-svm/src/regression.rs
index 77aa50a6f..7bc8b2b1d 100644
--- a/algorithms/linfa-svm/src/regression.rs
+++ b/algorithms/linfa-svm/src/regression.rs
@@ -221,7 +221,7 @@ pub mod tests {
     fn test_epsilon_regression_linear() -> Result<()> {
         // simple 2d straight line
         let targets = Array::linspace(0f64, 10., 100);
-        let records = targets.clone().into_shape((100, 1)).unwrap();
+        let records = targets.clone().into_shape_with_order((100, 1)).unwrap();
         let dataset = Dataset::new(records, targets);
 
         let model = Svm::params()
@@ -245,7 +245,7 @@ pub mod tests {
     fn test_nu_regression_linear() -> Result<()> {
         // simple 2d straight line
         let targets = Array::linspace(0f64, 10., 100);
-        let records = targets.clone().into_shape((100, 1)).unwrap();
+        let records = targets.clone().into_shape_with_order((100, 1)).unwrap();
         let dataset = Dataset::new(records, targets);
 
         // Test the precomputed dot product in the linear kernel case
@@ -268,9 +268,12 @@ pub mod tests {
     #[test]
     fn test_epsilon_regression_gaussian() -> Result<()> {
         let records = Array::linspace(0f64, 10., 100)
-            .into_shape((100, 1))
+            .into_shape_with_order((100, 1))
+            .unwrap();
+        let sin_curve = records
+            .mapv(|v| v.sin())
+            .into_shape_with_order((100,))
             .unwrap();
-        let sin_curve = records.mapv(|v| v.sin()).into_shape((100,)).unwrap();
         let dataset = Dataset::new(records, sin_curve);
 
         let model = Svm::params()
@@ -285,8 +288,13 @@ pub mod tests {
     #[test]
     fn test_nu_regression_polynomial() -> Result<()> {
         let n = 100;
-        let records = Array::linspace(0f64, 5., n).into_shape((n, 1)).unwrap();
-        let sin_curve = records.mapv(|v| v.sin()).into_shape((n,)).unwrap();
+        let records = Array::linspace(0f64, 5., n)
+            .into_shape_with_order((n, 1))
+            .unwrap();
+        let sin_curve = records
+            .mapv(|v| v.sin())
+            .into_shape_with_order((n,))
+            .unwrap();
         let dataset = Dataset::new(records, sin_curve);
 
         let model = Svm::params()
diff --git a/algorithms/linfa-trees/Cargo.toml b/algorithms/linfa-trees/Cargo.toml
index a270ec6af..b0907c339 100644
--- a/algorithms/linfa-trees/Cargo.toml
+++ b/algorithms/linfa-trees/Cargo.toml
@@ -24,16 +24,18 @@ default-features = false
 features = ["std", "derive"]
 
 [dependencies]
-ndarray = { version = "0.15" , features = ["rayon", "approx"]}
-ndarray-rand = "0.14"
+ndarray = { version = "0.16", features = ["rayon", "approx"] }
+ndarray-rand = "0.15"
 
 linfa = { version = "0.7.1", path = "../.." }
 
 [dev-dependencies]
 rand = { version = "0.8", features = ["small_rng"] }
 criterion = "0.4.0"
-approx = "0.4"
-linfa-datasets = { version = "0.7.1", path = "../../datasets/", features = ["iris"] }
+approx = "0.5"
+linfa-datasets = { version = "0.7.1", path = "../../datasets/", features = [
+    "iris",
+] }
 linfa = { version = "0.7.1", path = "../..", features = ["benchmarks"] }
 
 [[bench]]
diff --git a/algorithms/linfa-tsne/Cargo.toml b/algorithms/linfa-tsne/Cargo.toml
index 2cb7ebe7d..0d860d3c2 100644
--- a/algorithms/linfa-tsne/Cargo.toml
+++ b/algorithms/linfa-tsne/Cargo.toml
@@ -15,8 +15,8 @@ categories = ["algorithms", "mathematics", "science"]
 
 [dependencies]
 thiserror = "1.0"
-ndarray = { version = "0.15" }
-ndarray-rand = "0.14"
+ndarray = { version = "0.16" }
+ndarray-rand = "0.15"
 bhtsne = "0.4.0"
 pdqselect = "=0.1.0"
 
@@ -24,9 +24,11 @@ linfa = { version = "0.7.1", path = "../.." }
 
 [dev-dependencies]
 rand = "0.8"
-approx = "0.4"
+approx = "0.5"
 
-linfa-datasets = { version = "0.7.1", path = "../../datasets", features = ["iris"] }
+linfa-datasets = { version = "0.7.1", path = "../../datasets", features = [
+    "iris",
+] }
 linfa-reduction = { version = "0.7.1", path = "../linfa-reduction" }
 
 [target.'cfg(not(target_family = "windows"))'.dev-dependencies]
diff --git a/datasets/Cargo.toml b/datasets/Cargo.toml
index 77dbdd587..5ded9426f 100644
--- a/datasets/Cargo.toml
+++ b/datasets/Cargo.toml
@@ -9,14 +9,14 @@ repository = "https://github.com/rust-ml/linfa"
 
 [dependencies]
 linfa = { version = "0.7.1", path = ".." }
-ndarray = { version = "0.15" }
-ndarray-csv = "=0.5.1"
+ndarray = { version = "0.16" }
+ndarray-csv = "0.5"
 csv = "1.1"
 flate2 = "1.0"
-ndarray-rand = { version = "0.14", optional = true }
+ndarray-rand = { version = "0.15", optional = true }
 
 [dev-dependencies]
-approx = "0.4"
+approx = "0.5"
 statrs = "0.16.0"
 
 [features]
diff --git a/src/composing/multi_target_model.rs b/src/composing/multi_target_model.rs
index 31341d184..86c73d3e1 100644
--- a/src/composing/multi_target_model.rs
+++ b/src/composing/multi_target_model.rs
@@ -40,16 +40,21 @@ impl<L: Default, F, D: Data<Elem = F>> PredictInplace<ArrayBase<D, Ix2>, Array2<
             &[arr.nrows(), self.models.len()],
             "The number of data points must match the number of output targets."
         );
+        assert!(
+            targets.is_standard_layout(),
+            "targets not in row-major layout"
+        );
         *targets = self
             .models
             .iter()
             .flat_map(|model| {
                 let mut targets = Array1::default(arr.nrows());
                 model.predict_inplace(arr, &mut targets);
-                targets.into_raw_vec()
+                let (v, _) = targets.into_raw_vec_and_offset();
+                v
             })
             .collect::<Array1<L>>()
-            .into_shape((self.models.len(), arr.len_of(Axis(0))))
+            .into_shape_with_order((self.models.len(), arr.len_of(Axis(0))))
             .unwrap()
             .reversed_axes();
     }
diff --git a/src/dataset/impl_dataset.rs b/src/dataset/impl_dataset.rs
index 239a144e8..f36668c2f 100644
--- a/src/dataset/impl_dataset.rs
+++ b/src/dataset/impl_dataset.rs
@@ -124,7 +124,7 @@ impl<X, Y> Dataset<X, Y> {
     // Convert 2D targets to 1D. Only works for targets with shape of form [X, 1], panics otherwise.
     pub fn into_single_target(self) -> Dataset<X, Y, Ix1> {
         let nsamples = self.records.nsamples();
-        let targets = self.targets.into_shape(nsamples).unwrap();
+        let targets = self.targets.into_shape_with_order(nsamples).unwrap();
         let features = self.records;
         Dataset::new(features, targets)
     }
@@ -960,7 +960,7 @@ impl<F, E, I: TargetDim> Dataset<F, E, I> {
         let feature_names = self.feature_names();
 
         // split records into two disjoint arrays
-        let mut array_buf = self.records.into_raw_vec();
+        let (mut array_buf, _) = self.records.into_raw_vec_and_offset();
         let second_array_buf = array_buf.split_off(n1 * nfeatures);
 
         let first = Array2::from_shape_vec((n1, nfeatures), array_buf).unwrap();
@@ -969,7 +969,7 @@ impl<F, E, I: TargetDim> Dataset<F, E, I> {
         // split targets into two disjoint Vec
         let dim1 = self.targets.raw_dim().nsamples(n1);
         let dim2 = self.targets.raw_dim().nsamples(n2);
-        let mut array_buf = self.targets.into_raw_vec();
+        let (mut array_buf, _) = self.targets.into_raw_vec_and_offset();
         let second_array_buf = array_buf.split_off(dim1.size());
 
         let first_targets = Array::from_shape_vec(dim1, array_buf).unwrap();
@@ -977,7 +977,7 @@ impl<F, E, I: TargetDim> Dataset<F, E, I> {
 
         // split weights into two disjoint Vec
         let second_weights = if self.weights.len() == n1 + n2 {
-            let mut weights = self.weights.into_raw_vec();
+            let (mut weights, _) = self.weights.into_raw_vec_and_offset();
 
             let weights2 = weights.split_off(n1);
             self.weights = Array1::from(weights);