diff --git a/Cargo.toml b/Cargo.toml
index 9a92dedd..a74d39cf 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -73,6 +73,7 @@ matrixcompare = "0.1.3"
 # RNG for certain tests. However, different versions of proptest may give different sequences of numbers,
 # which may cause more brittle tests (although ideally they should take enough samples for it not to matter).
 proptest = { version = "=0.10.1" }
+itertools = "0.9"
 
 [workspace]
 members = [ "nalgebra-lapack", "nalgebra-glm" ]
diff --git a/tests/proptest/mod.rs b/tests/proptest/mod.rs
index a5068344..40d61864 100644
--- a/tests/proptest/mod.rs
+++ b/tests/proptest/mod.rs
@@ -1,10 +1,13 @@
 //! Tests for proptest-related functionality.
 use nalgebra::base::dimension::*;
 use nalgebra::proptest::{matrix, DimRange, MatrixStrategy};
-use nalgebra::{DMatrix, DVector, Dim, Matrix3, Matrix4, MatrixMN, Vector3};
+use nalgebra::{DMatrix, DVector, Dim, Matrix3, MatrixMN, Vector3};
 use proptest::prelude::*;
 use proptest::strategy::ValueTree;
 use proptest::test_runner::TestRunner;
+use itertools::Itertools;
+use std::iter::repeat;
+use std::collections::HashSet;
 
 /// Generate a proptest that tests that all matrices generated with the
 /// provided rows and columns conform to the constraints defined by the
@@ -97,35 +100,51 @@ fn matrix_samples_all_possible_outputs() {
     // We use a deterministic test runner to make the test "stable".
     let mut runner = TestRunner::deterministic();
 
-    let strategy = matrix(0..=2usize, 0..=3, 0..=3);
+    // This number needs to be high enough so that we with high probability sample
+    // all possible cases
+    let num_generated_matrices = 200000;
 
-    // We use flags to record whether values and combinations of dimensions were encountered.
-    // For example, if we encounter value 1, we set the value flag of 1 to true,
-    // and if we encounted matrix dimensions 4x3, we set the flag of [4, 3] to true.
-    let mut value_encountered = Vector3::new(false, false, false);
-    let mut dimensions_encountered = Matrix4::repeat(false);
+    let values = -1..=1;
+    let rows = 0..=2;
+    let cols = 0..=3;
+    let strategy = matrix(values.clone(), rows.clone(), cols.clone());
 
-    for _ in 0..1000 {
+    // Enumerate all possible combinations
+    let mut all_combinations = HashSet::new();
+    for nrows in rows {
+        for ncols in cols.clone() {
+            // For the given number of rows and columns
+            let n_values = nrows * ncols;
+
+            if n_values == 0 {
+                // If we have zero rows or columns, the set of matrices with the given
+                // rows and columns is a single element: an empty matrix
+                all_combinations.insert(DMatrix::from_row_slice(nrows, ncols, &[]));
+            } else {
+                // Otherwise, we need to sample all possible matrices.
+                // To do this, we generate the values as the (multi) Cartesian product
+                // of the value sets. For example, for a 2x2 matrices, we consider
+                // all possible 4-element arrays that the matrices can take by
+                // considering all elements in the cartesian product
+                //  V x V x V x V
+                // where V is the set of eligible values, e.g. V := -1 ..= 1
+                for matrix_values in repeat(values.clone()).take(n_values).multi_cartesian_product() {
+                    all_combinations.insert(DMatrix::from_row_slice(nrows, ncols, &matrix_values));
+                }
+            }
+        }
+    }
+
+    let mut visited_combinations = HashSet::new();
+    for _ in 0..num_generated_matrices {
         let tree = strategy
             .new_tree(&mut runner)
             .expect("Tree generation should not fail");
         let matrix = tree.current();
-
-        dimensions_encountered[(matrix.nrows(), matrix.ncols())] = true;
-
-        for &value in matrix.iter() {
-            value_encountered[value] = true;
-        }
+        visited_combinations.insert(matrix.clone());
     }
 
-    assert!(
-        value_encountered.iter().all(|v| *v),
-        "Did not sample all possible values."
-    );
-    assert!(
-        dimensions_encountered.iter().all(|v| *v),
-        "Did not sample all possible matrix dimensions."
-    );
+    assert_eq!(visited_combinations, all_combinations, "Did not sample all possible values.");
 }
 
 #[test]