Scalar and Feature-Banded Regularization¶

Regularization controls the variance of the estimated kernel. Select it using training data only. A scalar penalty is the default; feature-banded penalties are useful when prespecified predictor groups need different shrinkage.

In [1]:

import matplotlib
matplotlib.use("module://matplotlib_inline.backend_inline")
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import fftconvolve

from fftrf import TRF, pearsonr, r2_score

import matplotlib matplotlib.use("module://matplotlib_inline.backend_inline") import matplotlib.pyplot as plt import numpy as np from scipy.signal import fftconvolve from fftrf import TRF, pearsonr, r2_score

Simulate two predictor groups¶

The first feature carries a strong early response and the second a weaker later response. We will allow one penalty per feature.

In [2]:

rng = np.random.default_rng(33)
fs = 128.0
tmin, tmax = 0.0, 0.25
times = np.arange(round(tmax * fs)) / fs
kernels = np.stack([
    0.9 * np.exp(-0.5 * ((times - 0.055) / 0.015) ** 2),
    -0.3 * np.exp(-0.5 * ((times - 0.15) / 0.03) ** 2),
])

stimulus, response = [], []
for _ in range(7):
    x = rng.standard_normal((2048, 2))
    y = sum(fftconvolve(x[:, i], kernels[i], mode="full")[:2048]
            for i in range(2))
    y += 0.25 * rng.standard_normal(2048)
    stimulus.append(x)
    response.append(y[:, None])

train_x, test_x = stimulus[:-1], stimulus[-1]
train_y, test_y = response[:-1], response[-1]

rng = np.random.default_rng(33) fs = 128.0 tmin, tmax = 0.0, 0.25 times = np.arange(round(tmax * fs)) / fs kernels = np.stack([ 0.9 * np.exp(-0.5 * ((times - 0.055) / 0.015) ** 2), -0.3 * np.exp(-0.5 * ((times - 0.15) / 0.03) ** 2), ]) stimulus, response = [], [] for _ in range(7): x = rng.standard_normal((2048, 2)) y = sum(fftconvolve(x[:, i], kernels[i], mode="full")[:2048] for i in range(2)) y += 0.25 * rng.standard_normal(2048) stimulus.append(x) response.append(y[:, None]) train_x, test_x = stimulus[:-1], stimulus[-1] train_y, test_y = response[:-1], response[-1]

Cross-validate a scalar ridge grid¶

Start here unless separate predictor penalties were part of the analysis plan.

In [3]:

ridge_grid = np.logspace(-4, 0, 5)
scalar_model = TRF(direction=1)
scalar_scores = scalar_model.train(
    train_x,
    train_y,
    fs=fs,
    tmin=tmin,
    tmax=tmax,
    regularization=ridge_grid,
    k=4,
    seed=7,
)
_, scalar_r = scalar_model.predict(test_x, test_y)
print("scalar lambda:", scalar_model.regularization)
print(f"scalar held-out r: {float(scalar_r):.3f}")

ridge_grid = np.logspace(-4, 0, 5) scalar_model = TRF(direction=1) scalar_scores = scalar_model.train( train_x, train_y, fs=fs, tmin=tmin, tmax=tmax, regularization=ridge_grid, k=4, seed=7, ) _, scalar_r = scalar_model.predict(test_x, test_y) print("scalar lambda:", scalar_model.regularization) print(f"scalar held-out r: {float(scalar_r):.3f}")

scalar lambda: 1.0
scalar held-out r: 0.989

Select one penalty per prespecified feature band¶

bands=[1, 1] assigns one feature to each band. With five scalar candidates this evaluates 25 coefficient pairs. Candidate growth is multiplicative, so keep the grouping scientifically motivated and the grid compact.

In [4]:

banded_model = TRF(direction=1)
banded_scores = banded_model.train(
    train_x,
    train_y,
    fs=fs,
    tmin=tmin,
    tmax=tmax,
    regularization=ridge_grid,
    bands=[1, 1],
    k=4,
    seed=7,
)
_, banded_r = banded_model.predict(test_x, test_y)
score_grid = np.asarray(banded_scores).reshape(len(ridge_grid), len(ridge_grid))

print("selected band coefficients:", banded_model.regularization)
print("expanded feature penalties:", banded_model.feature_regularization)
print(f"banded held-out r: {float(banded_r):.3f}")

banded_model = TRF(direction=1) banded_scores = banded_model.train( train_x, train_y, fs=fs, tmin=tmin, tmax=tmax, regularization=ridge_grid, bands=[1, 1], k=4, seed=7, ) _, banded_r = banded_model.predict(test_x, test_y) score_grid = np.asarray(banded_scores).reshape(len(ridge_grid), len(ridge_grid)) print("selected band coefficients:", banded_model.regularization) print("expanded feature penalties:", banded_model.feature_regularization) print(f"banded held-out r: {float(banded_r):.3f}")

selected band coefficients: (1.0, 0.0001)
expanded feature penalties: [1.e+00 1.e-04]
banded held-out r: 0.989

Inspect the selection surface¶

A broad plateau implies the exact winning pair is not very stable. Report the grid and selection procedure, not only the chosen coefficients.

In [5]:

fig, ax = plt.subplots(figsize=(6, 5))
image = ax.imshow(score_grid, origin="lower", aspect="auto", cmap="viridis")
ax.set_xticks(range(len(ridge_grid)), [f"{v:.0e}" for v in ridge_grid])
ax.set_yticks(range(len(ridge_grid)), [f"{v:.0e}" for v in ridge_grid])
ax.set(xlabel="Feature 2 penalty", ylabel="Feature 1 penalty",
       title="Banded-ridge CV score")
fig.colorbar(image, ax=ax, label="Mean CV score")
plt.show()

fig, ax = plt.subplots(figsize=(6, 5)) image = ax.imshow(score_grid, origin="lower", aspect="auto", cmap="viridis") ax.set_xticks(range(len(ridge_grid)), [f"{v:.0e}" for v in ridge_grid]) ax.set_yticks(range(len(ridge_grid)), [f"{v:.0e}" for v in ridge_grid]) ax.set(xlabel="Feature 2 penalty", ylabel="Feature 1 penalty", title="Banded-ridge CV score") fig.colorbar(image, ax=ax, label="Mean CV score") plt.show()

See the Regularization and CV guide for leakage, fold, metric, and trial-weighting considerations.