FloatingPointError: overflow encountered in matmul with NumPy 2.x and SciPy 1.13+

[clivehywong]

python-picard fails with NumPy 2.x due to numerical overflow errors in scipy.linalg.expm during the optimization process.

Environment

• python-picard: 0.8.1 (also tested 0.8)
• NumPy: 2.2.6 (also tested 2.2.4)
• SciPy: 1.13.1 - 1.16.2
• Python: 3.10.13
• OS: macOS (also affects Linux)

Minimal Reproducible Example

import numpy as np
from picard import picard

# Create test data
np.random.seed(42)
n_channels, n_samples = 20, 5000
X = np.random.randn(n_channels, n_samples) * 0.01

# Run PICARD - FAILS with NumPy 2.x
K, W, Y = picard(X, max_iter=100)

Error Output

FloatingPointError: overflow encountered in matmul

Traceback:
  File "picard/solver.py", line 206, in picard
    Y, W, infos = core_picard(X1, **kwargs)
  File "picard/_core_picard.py", line 78, in core_picard
    current_loss = _loss(Y, W, density, signs, ortho, extended)
  File "picard/_core_picard.py", line 180, in _loss
    loss = - np.linalg.slogdet(W)[1]
  File "numpy/linalg/_linalg.py", line 2325, in slogdet
    sign, logdet = _umath_linalg.slogdet(a, signature=signature)
FloatingPointError: overflow encountered in slogdet

Alternative error location (in line search):

File "picard/_core_picard.py", line XX, in _line_search
    transform = expm(alpha * direction)
File "scipy/linalg/_matfuncs_expm.pyx", line 131, in scipy.linalg._matfuncs_expm.pick_pade_structure
FloatingPointError: overflow encountered in matmul

Root Cause Analysis

The PICARD algorithm creates poorly conditioned matrices during optimization, which cause numerical overflow when computing:

1. Matrix exponentials via scipy.linalg.expm
2. Log-determinants via np.linalg.slogdet

NumPy 1.x behavior: Silently returned inf/nan and continued (hiding the bug)

NumPy 2.x behavior: Raises FloatingPointError on overflow in compiled code (exposing the bug)

The error occurs in Cython-compiled code (scipy.linalg._matfuncs_expm), so it cannot be suppressed using np.seterr() or Python warnings.

Testing Details

Works:

• NumPy 1.26.x with SciPy 1.11.x ✅
• MNE-Python's FastICA (alternative algorithm) ✅

Fails:

• NumPy 2.x (all versions: 2.0.0 - 2.2.6) ❌
• All SciPy versions compatible with NumPy 2.x (1.13.0+) ❌
• Both picard() function and when called via MNE or eegprep ❌

Impact

This breaks compatibility with:

• Modern scientific Python stacks (NumPy 2.x is now standard)
• MNE-Python ICA with method='picard'
• eegprep (EEGLAB Python port) ICA
• Any package using picard with NumPy 2.x