[DiadkoShmek]

What does this implement/fix?

Preserve singleton batch dimensions in covariance stacks. The generic implementation squeezes every size-one axis of the result, including axes describing groups/windows, so a later batch reduction can silently average covariance rows instead.

For input (group=2, window=1, variable=2, observation=3), the documented result shape is (2, 1, 2, 2). Main returns (2, 2, 2). Consequently, mean(axis=1) produces [[1.5, 3], [6, 12]] for a simple two-group example instead of retaining the two covariance matrices [[1, 2], [2, 4]] and [[4, 8], [8, 16]].

Squeeze only the two trailing covariance axes. This retains the existing scalar-variance convention for one variable while preserving all leading batch axes. Parameterize the existing public batch test over singleton/multiple batch dimensions and one/three variables; clarify the one-variable return convention.

Additional information

Batch support was introduced in gh-507; this repairs its return-shape handling.

• The identical NumPy/read-only/strict/strictest covariance selection gives 24 failures / 111 passing controls / one existing skip on the exact base, and 135 passes / one skip after the fix.
• Separate per-population scalar checks cover real/complex observations, weights, corrections and downstream batch reductions. Independent review also used a pairwise-difference oracle, including zero-sized batches, zero/one observations and the one-variable scalar convention; no new blocker was found.
• Ruff check/format and git diff --check pass on all three changed files.

Optional/lazy/accelerator namespaces and the full project CI matrix were not executed locally.

AI Generation Disclosure

OpenAI Codex generated the implementation change, test parametrization, docstring clarification and this description, and performed source investigation, local execution and separate AI review.