Merge PR #72: Add Tikhonov regularization for MVAR matrix inversion

Replaces direct matrix inverse with Tikhonov-regularized solve to prevent
LinAlgError exceptions when computing Granger causality with near-singular
matrices.

Changes:
- Replaced xp.linalg.inv() with regularized solve in MVAR computations
- Added TIKHONOV_REGULARIZATION_FACTOR constant (1e-12)
- Scale-aware regularization: λ = factor × mean(||H||²)
- Added stress test for near-singular matrices

Impact:
- Granger causality measures handle highly correlated signals gracefully
- No more LinAlgError crashes with near-singular transfer functions
- Maintains numerical accuracy for well-conditioned cases

Resolves: #72

# Conflicts:
#	CHANGELOG.md
#	tests/test_connectivity.py

Author: edeno

CHANGELOG.md

@@ -83,6 +83,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
     - New (correct): 513 bins, includes 500 Hz (Nyquist)
   - See PR #71 for detailed analysis and discussion
 
+- **Tikhonov regularization for MVAR matrix inversion stability**:
+  - Replaced direct matrix inverse (`xp.linalg.inv()`) with Tikhonov-regularized solve in MVAR computations
+  - Prevents `LinAlgError` exceptions when computing Granger causality with near-singular matrices
+  - Affected functions: `_MVAR_Fourier_coefficients` property and `_estimate_transfer_function` function
+  - Uses scale-aware regularization: λ = `TIKHONOV_REGULARIZATION_FACTOR` × mean(||H||²)
+  - Added module-level constant `TIKHONOV_REGULARIZATION_FACTOR = 1e-12` for consistency
+  - Solves `(H + λI)x = I` instead of computing `inv(H)` for better numerical stability
+  - Added stress test `test_mvar_regularized_inverse_near_singular()` validating near-singular cases
+  - All Granger causality measures now handle highly correlated signals gracefully
+  - See PR #72 for detailed numerical analysis
+
 - CHANGELOG.md to track version changes following Keep a Changelog format
 - Ruff linter configuration for faster, more comprehensive Python linting
 - Enhanced package metadata with additional project URLs (Changelog, Source Code, Issue Tracker)

spectral_connectivity/connectivity.py

@@ -71,6 +71,10 @@
     "time_trials_tapers": partial(xp.mean, axis=(0, 1, 2)),
 }
 
+# Tikhonov regularization factor for stabilizing matrix inversions
+# Used to prevent numerical instability with near-singular matrices
+TIKHONOV_REGULARIZATION_FACTOR = 1e-12
+
 
 def _asnumpy(connectivity_measure: Callable) -> Callable:
     """Transform cupy array to numpy array.
@@ -569,7 +573,13 @@ def _noise_covariance(self) -> NDArray[np.floating]:
 
     @property
     def _MVAR_Fourier_coefficients(self) -> NDArray[np.complexfloating]:
-        return xp.linalg.inv(self._transfer_function)
+        H = self._transfer_function
+        # Tikhonov regularization: solve(H + λI, I) instead of inv(H)
+        # Scale-aware regularization parameter
+        lam = TIKHONOV_REGULARIZATION_FACTOR * xp.mean(xp.real(xp.conj(H) * H))
+        identity = xp.eye(H.shape[-1], dtype=H.dtype)
+        regularized_H = H + lam * identity
+        return xp.linalg.solve(regularized_H, identity)
 
     @property
     def _expectation(self) -> Callable:
@@ -1711,9 +1721,13 @@ def _estimate_transfer_function(
 
     """
     inverse_fourier_coefficients = ifft(minimum_phase, axis=-3).real
-    return xp.matmul(
-        minimum_phase, xp.linalg.inv(inverse_fourier_coefficients[..., 0:1, :, :])
-    )
+    H_0 = inverse_fourier_coefficients[..., 0:1, :, :]
+    # Tikhonov regularization: solve(H_0 + λI, I) instead of inv(H_0)
+    lam = TIKHONOV_REGULARIZATION_FACTOR * xp.mean(H_0 * H_0)  # Scale-aware regularization for real matrix
+    identity = xp.eye(H_0.shape[-1], dtype=H_0.dtype)
+    regularized_H_0 = H_0 + lam * identity
+    H_0_inv = xp.linalg.solve(regularized_H_0, identity)
+    return xp.matmul(minimum_phase, H_0_inv)
 
 
 def _estimate_predictive_power(

tests/test_connectivity.py

@@ -652,6 +652,59 @@ def test_nyquist_bin_odd_n():
     ), f"Expected {expected_n_frequencies} frequencies, got {coherence.shape[-3]}"
 
 
+def test_mvar_regularized_inverse_near_singular():
+    """Test regularized inverse handles near-singular frequency bins."""
+    np.random.seed(42)
+    n_time_samples, n_trials, n_tapers, n_fft_samples, n_signals = (
+        1, 10, 1, 5, 3
+    )
+
+    # Create nearly singular Fourier coefficients by making signals
+    # highly correlated
+    fourier_coefficients = np.zeros(
+        (n_time_samples, n_trials, n_tapers, n_fft_samples, n_signals),
+        dtype=complex,
+    )
+
+    # Base signal
+    base_signal = np.random.randn(
+        n_time_samples, n_trials, n_tapers, n_fft_samples
+    ) + 1j * np.random.randn(
+        n_time_samples, n_trials, n_tapers, n_fft_samples
+    )
+
+    # Create near-singular cross-spectral matrix by making signals
+    # nearly dependent
+    fourier_coefficients[..., 0] = base_signal
+    fourier_coefficients[..., 1] = base_signal + 1e-10 * (
+        np.random.randn(n_time_samples, n_trials, n_tapers, n_fft_samples)
+        + 1j * np.random.randn(n_time_samples, n_trials, n_tapers, n_fft_samples)
+    )
+    fourier_coefficients[..., 2] = base_signal + 1e-10 * (
+        np.random.randn(n_time_samples, n_trials, n_tapers, n_fft_samples)
+        + 1j * np.random.randn(n_time_samples, n_trials, n_tapers, n_fft_samples)
+    )
+
+    # This should not raise LinAlgError with regularized inverse
+    conn = Connectivity(fourier_coefficients=fourier_coefficients)
+
+    # Test that MVAR coefficients are computed without error
+    mvar_coeffs = conn._MVAR_Fourier_coefficients
+    assert mvar_coeffs is not None
+    assert np.all(np.isfinite(mvar_coeffs))
+
+    # Test that transfer function is computed without error
+    transfer_func = conn._transfer_function
+    assert transfer_func is not None
+    assert np.all(np.isfinite(transfer_func))
+
+    # Test connectivity measures that depend on MVAR work
+    dtf = conn.directed_transfer_function()
+    assert np.all(np.isfinite(dtf))
+    assert np.all(dtf >= 0)  # DTF should be non-negative
+    assert np.all(dtf <= 1)  # DTF should be bounded by 1
+
+
 def test_connectivity_rejects_wrong_ndim():
     """Test that Connectivity rejects inputs with wrong number of dimensions."""
     import pytest