[Bug] Numerical instability in `PearsonCorrelation` due to naive variance formula

[Prathamesh8989]

The current implementation of PearsonCorrelation in ignite.metrics.regression uses the naive "sum of squares" algorithm to compute variance and covariance:

[
Var(X) = E[X^2] - (E[X])^2
]

Although mathematically valid, this approach is numerically unstable when the input values have large magnitudes relative to their variance.

This leads to catastrophic cancellation, where two very large numbers (E[X^2] and (E[X])^2) are subtracted, causing loss of precision in float32. As a result, the metric can produce incorrect results, such as returning 0.0 when the true correlation is ≈ 0.707.

---

Steps to Reproduce

import torch
from ignite.metrics.regression import PearsonCorrelation
import numpy as np

# Magnitude that triggers precision loss in float32
offset = 1e8 

y_true = torch.tensor([1.0, 2.0, 3.0, 4.0, 5.0]) + offset
y_pred = torch.tensor([1.1, 2.1, 3.1, 4.1, 5.1]) + offset

metric = PearsonCorrelation()
metric.update((y_pred, y_true))
ignite_res = metric.compute()

# Ground truth using PyTorch's stable corrcoef
combined = torch.stack([y_pred, y_true])
torch_res = torch.corrcoef(combined)[0, 1].item()

print(f"Offset used: {offset}")
print(f"Ignite Result: {ignite_res}")
print(f"Torch Ground Truth: {torch_res}")

---

actual output

Offset used: 100000000.0
Ignite Result: 0.0
Torch Ground Truth: 0.7071067690849304

---

Expected Behavior

The Pearson correlation metric should be invariant to constant offsets in the input data.

Therefore, the result should be consistent with numerically stable implementations such as:

• torch.corrcoef
• scipy.stats.pearsonr

Expected Result

≈ 0.7071067

---

Proposed Fix

To improve numerical stability, the implementation can be updated as follows:

Use a Numerically Stable Online Algorithm

Implement Welford’s Online Algorithm (or a similar one-pass algorithm) to compute:

• Mean
• Variance
• Covariance

incrementally.

---

Use Higher Precision Accumulators

Maintain internal accumulators using torch.float64 to prevent precision loss during batch updates.

---

Ensure Compatibility with Ignite Metrics

The updated implementation should still support batch-wise streaming updates used by Ignite metrics.

---

Benefit

These changes will prevent catastrophic cancellation and ensure correct results even when the data has:

• Large offsets
• High magnitude values

---

Hi @vfdev-5, I would like to work on a PR to resolve this issue if this looks good to you!

[TahaZahid05]

@Prathamesh8989 Thanks for the issue! The issue of instability persists across other metrics as well such as in metrics/regression/r2_score.py, and metrics/gan/fid.py. I'll say maybe a shared utility fix would be better.

@vfdev-5 thoughts?

[Prathamesh8989]

Hey @TahaZahid05,
Thanks for pointing that out! You're right that the same numerical instability could affect other metrics such as r2_score.py and fid.py.

A shared utility for numerically stable variance/covariance computation would probably make the implementation cleaner and avoid duplicating fixes across modules. I can take a look at how these metrics currently compute their statistics and explore introducing a shared helper (possibly using a stable online algorithm like Welford’s).

If that sounds reasonable, I can start by prototyping the change for PearsonCorrelation and then see how it can be generalized for the other affected metrics.

[Prathamesh8989]

Hi @TahaZahid05, @vfdev-5,

Following up on the discussion regarding a shared utility, I’ve spent some time researching a robust implementation that supports Ignite's distributed architecture.

To resolve the numerical instability across PearsonCorrelation, R2Score, and FID, I propose implementing a Parallelized Welford-based Accumulator.

---

The Technical Challenge

The current naive variance formulation

$$ E[X^2] - (E[X])^2 $$

is prone to catastrophic cancellation because the magnitude of the squared sum grows at

$$ O(N \cdot \mu^2) $$

In float32, this leads to a total loss of precision when values have large offsets (as seen in the $10^8$ case).

While a standard online Welford algorithm works for serial streams, Ignite requires distributed synchronization via idist.all_reduce.

---

Proposed Implementation

I suggest a utility (potentially in ignite.metrics.utils) that implements Chan et al.'s formula for merging statistics from distributed nodes.

Mean Merge

$$ \mu_{A \cup B} = \mu_A + \frac{n_B}{n_A + n_B}(\mu_B - \mu_A) $$

Variance Merge

$$ M2_{A \cup B} = M2_A + M2_B + \frac{n_A n_B}{n_A + n_B}(\mu_B - \mu_A)^2 $$

---

Architectural Plan

1. High-Precision Accumulation

Maintain internal state

$$ (\mu, M2, C_{xy}) $$

in torch.float64 to prevent overflow/underflow during large batch updates.
Cast back to the user's preferred dtype only during .compute().

---

2. Shared Logic Across Metrics

This utility can be reused across multiple metrics:

• PearsonCorrelation
  Use for stable

  • $Cov(X, Y)$
  • $Var(X)$
  • $Var(Y)$

• R2Score
  Use for the total sum of squares

$$ SS_{tot} $$

• FID
  Use for stable covariance matrix estimation across high-dimensional feature vectors.

---

Validation

I have verified this merging logic locally. It successfully handles the $10^8$ offset case provided in the reproduction script.

Now that I have submitted the Spearman PR (#3663 ), I am ready to prototype this shared utility if this approach aligns with the maintainers' vision.

[joemunene-by]

Hi! I'd like to work on this. I'll implement the fix using Welford's online algorithm to replace the naive variance formula, which should resolve the catastrophic cancellation issue in float32. I'll submit a PR shortly.

[TahaZahid05]

First, let us close #3672. Then, if @Prathamesh8989 wants, he can create the shared utility else someone else can take up on it.