From ecedab7ecbdc5fa1a60cc8bb80d600c954525259 Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Sat, 13 Jun 2026 17:29:53 +0200 Subject: [PATCH] fix: match MATLAB column-major average reference reref() computed the average reference as np.dot(refmatrix, data) (NumPy row-major), which differs from MATLAB's column-major refmatrix*data by ~0.06 uV in float32. ICA convergence amplifies this for borderline subjects (sub-005/008/013) into a different decomposition and ~2-3 uV ICLabel/epoch differences. Transposing the operands -- (data.T @ refmatrix).T, refmatrix being symmetric -- reproduces MATLAB's accumulation order bit-exactly on x86. reref now matches MATLAB to 0 uV across all 13 P300 parity subjects. --- src/eegprep/functions/sigprocfunc/reref.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/eegprep/functions/sigprocfunc/reref.py b/src/eegprep/functions/sigprocfunc/reref.py index bf0e601f..2edf4d44 100644 --- a/src/eegprep/functions/sigprocfunc/reref.py +++ b/src/eegprep/functions/sigprocfunc/reref.py @@ -80,15 +80,18 @@ def reref( else: # Average reference via matrix multiplication, matching MATLAB's # reref.m: refmatrix = eye(n) - ones(n)/n; data = refmatrix * data. - # This produces results numerically closer to MATLAB than the - # algebraically equivalent data - mean(data) because the float32 - # accumulation order in BLAS matrix multiplication matches MATLAB, while - # np.mean uses pairwise summation which can diverge enough to - # cascade through nonlinear downstream operations (e.g. ASR). + # MATLAB accumulates this product column-major; replicate that float32 + # accumulation order in NumPy (row-major) by transposing the operands + # (refmatrix is symmetric, so (data.T @ refmatrix).T == refmatrix @ data + # algebraically). The plain np.dot(refmatrix, data) form is ~0.06 uV off + # MATLAB in float32, which cascades through nonlinear downstream steps + # (ICA convergence) for borderline subjects; the transposed form is + # bit-exact to MATLAB on x86 BLAS. n = len(chansin) dt = original_dtype refmatrix = np.eye(n, dtype=dt) - np.ones((n, n), dtype=dt) / dt.type(n) - work[chansin_array, :] = np.dot(refmatrix, work[chansin_array, :].astype(dt)) + block = np.ascontiguousarray(work[chansin_array, :].astype(dt).T) + work[chansin_array, :] = (block @ refmatrix).T mean_data = None if locs is not None: