diff --git a/benchmarks/queries/clickbench/README.md b/benchmarks/queries/clickbench/README.md
index 2fd31c0b4c630..c0ee7bc327f29 100644
--- a/benchmarks/queries/clickbench/README.md
+++ b/benchmarks/queries/clickbench/README.md
@@ -260,7 +260,6 @@ FROM hits
 WHERE "URL" < 'zzzz';
 ```
 
-
 ### Q14: Grouped `COUNT(DISTINCT <string>)` beside a non-distinct `COUNT(*)`
 
 **Question**: "For the top 10 search phrases by hit count, how many distinct
@@ -301,6 +300,34 @@ LIMIT 10;
 
 
 
+### Q15: Correlation of ad slot dimensions per user
+
+**Question**: "How correlated are the width and height of the screens each user browses from?"
+
+**Important Query Properties**: `covar_samp` over a very high cardinality `GROUP BY` (about 17.6M distinct `UserID`s). Q16 is the low cardinality version of this query.
+
+```sql
+SELECT MAX(c) FROM (
+    SELECT covar_samp("ResolutionWidth", "ResolutionHeight") as c
+    FROM hits
+    GROUP BY "UserID"
+);
+```
+
+### Q16: Correlation of ad slot dimensions per region
+
+**Question**: "How correlated are the width and height of the screens browsed from in each region?"
+
+**Important Query Properties**: the same `covar_samp` aggregate as Q15, but over about 9,000 distinct `RegionID`s instead of 17.6M. It is the low cardinality version of Q15.
+
+```sql
+SELECT MAX(c) FROM (
+    SELECT covar_samp("ResolutionWidth", "ResolutionHeight") as c
+    FROM hits
+    GROUP BY "RegionID"
+);
+```
+
 ## Data Notes
 
 Here are some interesting statistics about the data used in the queries
diff --git a/benchmarks/queries/clickbench/extended/q15.sql b/benchmarks/queries/clickbench/extended/q15.sql
new file mode 100644
index 0000000000000..e648b8dd97992
--- /dev/null
+++ b/benchmarks/queries/clickbench/extended/q15.sql
@@ -0,0 +1,8 @@
+-- Must set for ClickBench hits_partitioned dataset. See https://github.com/apache/datafusion/issues/16591
+-- set datafusion.execution.parquet.binary_as_string = true
+
+SELECT MAX(c) FROM (
+    SELECT covar_samp("ResolutionWidth", "ResolutionHeight") as c
+    FROM hits
+    GROUP BY "UserID"
+);
diff --git a/benchmarks/queries/clickbench/extended/q16.sql b/benchmarks/queries/clickbench/extended/q16.sql
new file mode 100644
index 0000000000000..f3db0ae2df8ec
--- /dev/null
+++ b/benchmarks/queries/clickbench/extended/q16.sql
@@ -0,0 +1,8 @@
+-- Must set for ClickBench hits_partitioned dataset. See https://github.com/apache/datafusion/issues/16591
+-- set datafusion.execution.parquet.binary_as_string = true
+
+SELECT MAX(c) FROM (
+    SELECT covar_samp("ResolutionWidth", "ResolutionHeight") as c
+    FROM hits
+    GROUP BY "RegionID"
+);
diff --git a/datafusion/functions-aggregate-common/src/aggregate/groups_accumulator.rs b/datafusion/functions-aggregate-common/src/aggregate/groups_accumulator.rs
index 41917329b2df8..d49cae6a161c3 100644
--- a/datafusion/functions-aggregate-common/src/aggregate/groups_accumulator.rs
+++ b/datafusion/functions-aggregate-common/src/aggregate/groups_accumulator.rs
@@ -235,14 +235,23 @@ impl GroupsAccumulatorAdapter {
 
         assert_eq!(values[0].len(), group_indices.len());
 
+        // groups_with_rows holds a list of group indexes that have any rows
+        // that need to be accumulated, stored in order of first appearance in
+        // this batch
+        let mut groups_with_rows = vec![];
+
         // figure out which input rows correspond to which groups.
-        // Note that self.state.indices starts empty for all groups
-        // (it is cleared out below)
+        // Note that self.state.indices starts empty for all groups (it is
+        // cleared out below), so an empty one is exactly a group that this
+        // batch has not reached yet.
         // Charge retained scratch capacity only when a push grows a vector,
         // avoiding another pass over every group after indexing.
         let mut indices_allocation_delta = 0;
         for (idx, group_index) in group_indices.iter().enumerate() {
             let indices = &mut self.states[*group_index].indices;
+            if indices.is_empty() {
+                groups_with_rows.push(*group_index);
+            }
             if indices.len() < indices.capacity() {
                 indices.push(idx as u32);
             } else {
@@ -253,27 +262,17 @@ impl GroupsAccumulatorAdapter {
         }
         self.add_allocation(indices_allocation_delta);
 
-        // groups_with_rows holds a list of group indexes that have
-        // any rows that need to be accumulated, stored in order of
-        // group_index
-
-        let mut groups_with_rows = vec![];
-
         // batch_indices holds indices into values, each group is contiguous
-        let mut batch_indices = vec![];
+        let mut batch_indices = Vec::with_capacity(group_indices.len());
 
         // offsets[i] is index into batch_indices where the rows for
-        // group_index i starts
-        let mut offsets = vec![0];
+        // groups_with_rows[i] start
+        let mut offsets = Vec::with_capacity(groups_with_rows.len() + 1);
+        offsets.push(0);
 
         let mut offset_so_far = 0;
-        for (group_index, state) in self.states.iter_mut().enumerate() {
-            let indices = &state.indices;
-            if indices.is_empty() {
-                continue;
-            }
-
-            groups_with_rows.push(group_index);
+        for &group_index in groups_with_rows.iter() {
+            let indices = &self.states[group_index].indices;
             batch_indices.extend_from_slice(indices);
             offset_so_far += indices.len();
             offsets.push(offset_so_far);