diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8651734f..ee8a8742 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,10 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+### Fixed
+- The ndarray implementation of `Array2::from_iterator` now preserves column-major input order, including for rectangular matrices that previously caused a panic. Column-major input produces a column-major array (#482).
+
 ## [0.6.16]
 ### Fixed
 - `predict`, `predict_proba`, `predict_oob`, `predict_matrix` and `decision_function` now return `Err(Failed)` instead of panicking when `fit` has not run (#469, #470). Models return the same unfitted error even after deserialization, when the state fields are `None`. The guarded methods cover the decision trees, the random forests, the extra-trees regressor, the KNN classifier and regressor, `SVC`, `MultiClassSVC`, `SVR`, the linear models, the naive Bayes classifiers and `XGRegressor`.
diff --git a/src/linalg/ndarray/matrix.rs b/src/linalg/ndarray/matrix.rs
index 941d8062..ea540332 100644
--- a/src/linalg/ndarray/matrix.rs
+++ b/src/linalg/ndarray/matrix.rs
@@ -147,17 +147,14 @@ impl<T: Debug + Display + Copy + Sized> Array2<T> for ArrayBase<OwnedRepr<T>, Ix
     }
 
     fn from_iterator<I: Iterator<Item = T>>(iter: I, nrows: usize, ncols: usize, axis: u8) -> Self {
-        // `into_shape` was deprecated in ndarray 0.16; use `into_shape_with_order` instead.
-        let a = Array::from_iter(iter.take(nrows * ncols))
-            .into_shape_with_order(((nrows, ncols), Order::RowMajor))
-            .unwrap();
-        match axis {
-            0 => a,
-            _ => a
-                .reversed_axes()
-                .into_shape_with_order(((nrows, ncols), Order::RowMajor))
-                .unwrap(),
-        }
+        let order = if axis == 0 {
+            Order::RowMajor
+        } else {
+            Order::ColumnMajor
+        };
+        Array::from_iter(iter.take(nrows * ncols))
+            .into_shape_with_order(((nrows, ncols), order))
+            .unwrap()
     }
 
     fn transpose(&self) -> Self {
@@ -266,6 +263,38 @@ mod tests {
     use super::*;
     use ndarray::arr2;
 
+    #[test]
+    fn from_iterator_preserves_input_order() {
+        let row_major = <ndarray::Array2<i32> as Array2<i32>>::from_iterator(1..=6, 2, 3, 0);
+        assert_eq!(row_major, arr2(&[[1, 2, 3], [4, 5, 6]]));
+        assert_eq!(row_major.strides(), &[3, 1]);
+
+        let column_major = <ndarray::Array2<i32> as Array2<i32>>::from_iterator(1..=6, 2, 3, 1);
+        assert_eq!(column_major, arr2(&[[1, 3, 5], [2, 4, 6]]));
+        assert_eq!(column_major.strides(), &[1, 2]);
+
+        for (nrows, ncols) in [(0, 3), (3, 0)] {
+            for axis in [0, 1] {
+                let empty = <ndarray::Array2<i32> as Array2<i32>>::from_iterator(
+                    std::iter::empty(),
+                    nrows,
+                    ncols,
+                    axis,
+                );
+                assert_eq!(BaseArray::shape(&empty), (nrows, ncols));
+            }
+        }
+    }
+
+    #[test]
+    fn h_stack_preserves_column_order() {
+        let left = arr2(&[[1, 2], [4, 5]]);
+        let right = arr2(&[[3], [6]]);
+        let stacked = Array2::h_stack(&left, &right);
+
+        assert_eq!(stacked, arr2(&[[1, 2, 3], [4, 5, 6]]));
+    }
+
     #[test]
     fn test_dense_matrix_from_ndarray2() {
         let input = arr2(&[[1, 2, 3], [4, 5, 6]]);