Fix covariance weighting across ensembles and Covobs sources

[s-kuberski]

This pull request fixes a bug in the current implementation of the computation of a covariance matrix. It keeps intact the current behavior of the construction regarding the treatment of autocorrelation.

We had written

WARNING: This function should be used with care, especially for observables with support on multiple
ensembles with differing autocorrelations.

in the docstring of the function, but I think that this was not sufficient, given that the covariance was basically wrong when having multiple ensembles.

Problem

covariance() currently combines pairwise zero-lag correlations across ensembles before rescaling them by the total errors. This loses the relative size of the ensemble contributions. It also mixes Covobs covariances, which have units of variance, into that correlation calculation.

There is a related issue when observables have support on different configurations of the same ensemble. The current pairwise normalization uses only configurations shared by each pair. Those pairwise correlations need not form a positive-semidefinite matrix and do not describe the covariance of averages computed from the full available samples. I always had in mind that the old choice was better, because the correlation was estimated only on the common set and thus nicely resolved. However, the violation of positive semidefiniteness is something that we don't want in this function. This is an edge case anyways...

Both effects can be seen with small examples:

import numpy as np
import pyerrors as pe

# Two ensembles with contributions of different sizes.
data = np.arange(100, dtype=float)
x = pe.Obs([data], ["e1"])
z = pe.Obs([data], ["e2"])
y1, y2 = x + 0.1 * z, x - 0.1 * z
for obs in (y1, y2):
    obs.gamma_method(S=0)

print(round(pe.covariance([y1, y2], correlation=True)[0, 1], 3))
# Current: 0.000; with this fix: 0.980

# A and C are measured on disjoint halves; B is measured on both.
data = np.array([1.0, -1.0] * 5)
a = pe.Obs([data], ["e"], idl=[range(1, 11)])
b = pe.Obs([np.tile(data, 2)], ["e"], idl=[range(1, 21)])
c = pe.Obs([data], ["e"], idl=[range(11, 21)])
for obs in (a, b, c):
    obs.gamma_method(S=0)

corr = pe.covariance([a, b, c], correlation=True)
print(np.round(np.linalg.eigvalsh(corr), 3))
# Current: [-0.414, 1.000, 2.414]; with this fix: [0.000, 1.000, 2.000]

Fix

Construct the zero-lag covariance matrix separately for each ensemble. Its $(i,j)$ element sums products of fluctuations on configurations shared by observables $i$ and $j$, while its normalization uses each observable’s full sample count on that ensemble. This gives a single zero-lag matrix for the ensemble, rather than separately normalized correlations for each pair.

Normalize that matrix to obtain $R_e^{(0)}$, then rescale it with the ensemble-specific gamma-method errors:
$$C_{\mathrm{MC}}=\sum_e D_e R_e^{(0)}D_e.$$
Here $D_e$ is diagonal, with entries equal to the contribution of ensemble $e$ to each observable’s error. Add shared Covobs inputs in variance units:
$$C=C_{\mathrm{MC}}+\sum_c G_c\Sigma_cG_c^{\mathsf T},$$
where $\Sigma_c$ is the input covariance and $G_c$ contains the propagated gradients.

Each zero-lag ensemble matrix is a Gram matrix when fluctuations are placed on the common configuration space, with zero entries where an observable was not measured. Normalization and diagonal rescaling preserve positive semidefiniteness; so does adding the Covobs covariance matrices.

The routine retains its existing approximation for cross-autocorrelations: it uses zero-lag correlations and the gamma-method errors of the individual observables. It does not estimate correlations at nonzero lags. For one ensemble and one common configuration set, the result agrees with the previous construction. Results can change for unequal configuration support, differing replica contributions, or multiple error sources, as illustrated above.

The tests cover these cases, the covariance of identical observables, and linear propagation across independent ensembles and shared Covobs inputs.

Future

A more detailed estimator of cross-autocorrelations should be considered separately. I have seen cases where autocorrelation in the cross-covariance is not negligible and there might be better ways to include this effect without ending up with non-positive matrices all of the time. It could be an idea to have a switch in order to support several different constructions.

[Copilot]

Copilot review overview

🟡 Changes recommended

Same-name Covobs matrices need consistency validation, alongside the documented performance and coverage concerns.

Review effort: Balanced
Findings: 1 · 1 · 2

Open (4)

• Validate covariance definitions for same-name Covobs instances · New
• Restrict covariance pairing to observables supported by each ensemble · New
• Account for per-observable normalization in PSD formula · New
• Test PSD with three observables having different overlaps · New

What changed in this PR

Fixes covariance weighting across Monte Carlo ensembles, unequal configuration support, and shared Covobs sources.

Changes:

• Computes and rescales zero-lag covariance per ensemble.
• Adds Covobs covariance in variance units.
• Expands regression coverage for mixed sources and sample support.

File	Description
pyerrors/​obs.py	Reworks covariance construction and documentation.
tests/​obs_test.py	Adds and updates covariance regression tests.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[s-kuberski]

I've addressed all of copilot's suggestions.

[Copilot]

🟢 Approval recommended

The implementation matches the documented estimator and is covered by focused regression tests.

0 open findings

4 resolved since last review

• Validate covariance definitions for same-name Covobs instances
• Restrict covariance pairing to observables supported by each ensemble
• Test PSD with three observables having different overlaps
• Account for per-observable normalization in PSD formula

🧠 Review effort: Balanced