diff --git a/src/rbc/core/metrics/__init__.py b/src/rbc/core/metrics/__init__.py
index 81f32c19..5b301f53 100644
--- a/src/rbc/core/metrics/__init__.py
+++ b/src/rbc/core/metrics/__init__.py
@@ -4,6 +4,10 @@
 """
 
 from rbc.core.metrics.atlases import get_atlas
+from rbc.core.metrics.gradients import (
+    compute_gradients,
+    compute_gradients_from_files,
+)
 from rbc.core.metrics.timeseries import (
     compute_timeseries,
     correlation_matrix,
@@ -11,6 +15,8 @@
 )
 
 __all__ = [
+    "compute_gradients",
+    "compute_gradients_from_files",
     "compute_timeseries",
     "correlation_matrix",
     "extract_timeseries",
diff --git a/src/rbc/core/metrics/gradients.py b/src/rbc/core/metrics/gradients.py
new file mode 100644
index 00000000..ee091930
--- /dev/null
+++ b/src/rbc/core/metrics/gradients.py
@@ -0,0 +1,313 @@
+"""Connectivity gradient computation via diffusion map embedding.
+
+Computes continuous axes of brain organization from connectivity matrices
+using manifold learning.  The algorithm constructs an affinity matrix from
+pairwise correlations, then applies diffusion map embedding (Coifman &
+Lafon 2006) to recover low-dimensional gradients.
+
+Vendored from BrainSpace (Vos de Wael et al. 2020) to avoid pulling in
+VTK/matplotlib/nilearn.  Core math requires only numpy and scipy.
+"""
+
+from __future__ import annotations
+
+from pathlib import Path
+from typing import TYPE_CHECKING, NamedTuple
+
+import numpy as np
+from scipy.linalg import eigh
+
+if TYPE_CHECKING:
+    from typing import Literal
+
+    AffinityKernel = Literal["cosine", "normalized_angle", "gaussian"]
+
+
+class GradientResult(NamedTuple):
+    """Result of :func:`compute_gradients`."""
+
+    gradients: np.ndarray
+    """(n_rois, n_components) gradient scores."""
+    lambdas: np.ndarray
+    """(n_components,) eigenvalues (descending)."""
+    affinity: np.ndarray
+    """(n_rois, n_rois) affinity matrix used for embedding."""
+
+
+class GradientOutputs(NamedTuple):
+    """File outputs from :func:`compute_gradients_from_files`."""
+
+    gradients: Path
+    """Path to the gradients TSV file."""
+    lambdas: Path
+    """Path to the eigenvalues TSV file."""
+
+
+def _dominant_set(matrix: np.ndarray, k: int) -> np.ndarray:
+    """Keep the top *k* entries per row, zeroing the rest.
+
+    Matches BrainSpace's ``dominant_set`` (dense path).
+    """
+    n = matrix.shape[0]
+    idx = np.argpartition(matrix, n - k, axis=1)
+    row = np.arange(n)[:, None]
+    col = idx[:, -k:]
+    out = np.zeros_like(matrix)
+    out[row, col] = matrix[row, col]
+    return out
+
+
+def compute_affinity(
+    matrix: np.ndarray,
+    kernel: AffinityKernel = "normalized_angle",
+    sparsity: float = 0.9,
+) -> np.ndarray:
+    """Build a non-negative symmetric affinity matrix from a connectivity matrix.
+
+    Follows BrainSpace's default pipeline: sparsify the input matrix first,
+    then apply the kernel, then zero out negatives.
+
+    Args:
+        matrix: (n_rois, n_rois) correlation / connectivity matrix.
+        kernel: Similarity kernel -- ``"cosine"``, ``"normalized_angle"``
+            (default, matches BrainSpace), or ``"gaussian"``.
+        sparsity: Fraction of weakest connections to zero out per row.
+            Must be in [0, 1).  Default ``0.9`` retains the top 10%
+            of connections (BrainSpace default).
+
+    Returns:
+        (n_rois, n_rois) non-negative symmetric affinity matrix.
+
+    Raises:
+        ValueError: If *matrix* is not 2-D square, or *sparsity* is
+            outside [0, 1).
+    """
+    if matrix.ndim != 2 or matrix.shape[0] != matrix.shape[1]:
+        raise ValueError(f"Expected square 2D matrix, got shape {matrix.shape}")
+    if not (0.0 <= sparsity < 1.0):
+        raise ValueError(f"sparsity must be in [0, 1), got {sparsity}")
+
+    n = matrix.shape[0]
+    matrix = np.array(matrix, dtype=np.float64)
+    matrix = np.nan_to_num(matrix, nan=0.0)
+
+    # Sparsify input before kernel (BrainSpace default: pre_sparsify=True)
+    if sparsity > 0 and n > 1:
+        k = max(1, int(n * (1.0 - sparsity)))
+        matrix = _dominant_set(matrix, k)
+
+    # Compute pairwise cosine similarity between rows
+    norms = np.linalg.norm(matrix, axis=1, keepdims=True)
+    norms = np.where(norms == 0, 1.0, norms)
+    normed = matrix / norms
+    cosine_sim = normed @ normed.T
+    cosine_sim = np.clip(cosine_sim, -1.0, 1.0)
+
+    if kernel == "cosine":
+        affinity = cosine_sim
+    elif kernel == "normalized_angle":
+        affinity = 1.0 - np.arccos(cosine_sim) / np.pi
+    elif kernel == "gaussian":
+        # RBF kernel: exp(-gamma * ||x_i - x_j||^2), gamma = 1/n_features
+        norms_sq = np.sum(matrix**2, axis=1)
+        dist_sq = norms_sq[:, None] + norms_sq[None, :] - 2.0 * (matrix @ matrix.T)
+        dist_sq = np.maximum(dist_sq, 0.0)
+        gamma = 1.0 / matrix.shape[1]
+        affinity = np.exp(-gamma * dist_sq)
+    else:
+        raise ValueError(f"Unknown kernel {kernel!r}")
+
+    # Zero out negatives (BrainSpace: non_negative=True)
+    affinity[affinity < 0] = 0.0
+
+    return affinity
+
+
+def diffusion_map_embedding(
+    affinity: np.ndarray,
+    n_components: int = 10,
+    alpha: float = 0.5,
+    diffusion_time: int = 0,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Compute diffusion map embedding from an affinity matrix.
+
+    Follows BrainSpace's ``diffusion_mapping``.  Uses ``scipy.linalg.eigh``
+    on the symmetrized operator (rather than ``eigsh`` on the asymmetric
+    transition matrix) for dense numerical stability, then maps back to
+    the diffusion coordinate space with the same post-processing.
+
+    Args:
+        affinity: (n_rois, n_rois) non-negative symmetric affinity matrix.
+        n_components: Number of gradient components to return.
+        alpha: Anisotropic diffusion parameter (0 = classical, 1 = Laplace-
+            Beltrami normalization, 0.5 = default).
+        diffusion_time: Diffusion time.  0 uses multi-scale diffusion maps
+            (BrainSpace default).
+
+    Returns:
+        ``(gradients, lambdas)`` -- (n_rois, n_components) float64 gradient
+        scores and (n_components,) eigenvalues in descending order.
+
+    Raises:
+        ValueError: If *affinity* is not 2-D square or *n_components* exceeds
+            the number of ROIs minus one.
+    """
+    if affinity.ndim != 2 or affinity.shape[0] != affinity.shape[1]:
+        raise ValueError(f"Expected square 2D affinity, got shape {affinity.shape}")
+    n = affinity.shape[0]
+    if n_components >= n:
+        raise ValueError(f"n_components ({n_components}) must be < n_rois ({n})")
+
+    affinity = np.array(affinity, dtype=np.float64)
+
+    # Anisotropic diffusion: W_alpha = D^{-alpha} W D^{-alpha}
+    if alpha > 0:
+        d = affinity.sum(axis=1)
+        d = np.where(d == 0, 1.0, d)
+        d_alpha = d ** (-alpha)
+        affinity = affinity * d_alpha[:, None] * d_alpha[None, :]
+
+    # Form symmetric operator Ms = D^{-1/2} W_alpha D^{-1/2}
+    # (equivalent eigenvalues to transition matrix P = D^{-1} W_alpha,
+    # but eigh on the symmetric form is more stable for dense matrices)
+    d2 = affinity.sum(axis=1)
+    d2 = np.where(d2 == 0, 1.0, d2)
+    d_inv_sqrt = d2 ** (-0.5)
+    ms = affinity * d_inv_sqrt[:, None] * d_inv_sqrt[None, :]
+    ms = (ms + ms.T) / 2.0
+
+    eigenvalues, eigenvectors = eigh(ms)
+    eigenvalues = eigenvalues[::-1]
+    eigenvectors = eigenvectors[:, ::-1]
+
+    # Map symmetric eigenvectors back to transition-matrix space
+    v = eigenvectors * d_inv_sqrt[:, None]
+
+    # Normalize by first eigenvector (converts to diffusion coordinates)
+    v /= v[:, [0]]
+    eigenvalues /= eigenvalues[0]
+
+    # Drop first (trivial) component
+    eigenvalues = eigenvalues[1 : n_components + 1]
+    v = v[:, 1 : n_components + 1]
+
+    # Diffusion time scaling
+    if diffusion_time <= 0:
+        # Multi-scale diffusion map (BrainSpace default)
+        eigenvalues = eigenvalues / (1.0 - eigenvalues)
+    else:
+        eigenvalues = eigenvalues**diffusion_time
+
+    # Rescale eigenvectors by eigenvalues
+    v *= eigenvalues[None, :]
+
+    # Sign convention: largest-absolute-value element is positive
+    signs = np.sign(v[np.argmax(np.abs(v), axis=0), np.arange(v.shape[1])])
+    v *= signs[None, :]
+
+    return v, eigenvalues
+
+
+def compute_gradients(
+    matrix: np.ndarray,
+    n_components: int = 10,
+    kernel: AffinityKernel = "normalized_angle",
+    sparsity: float = 0.9,
+    alpha: float = 0.5,
+    diffusion_time: int = 0,
+) -> GradientResult:
+    """Compute connectivity gradients from a correlation matrix.
+
+    Combines :func:`compute_affinity` and :func:`diffusion_map_embedding`
+    into a single call.
+
+    Args:
+        matrix: (n_rois, n_rois) symmetric connectivity / correlation matrix.
+        n_components: Number of gradient components.
+        kernel: Affinity kernel (see :func:`compute_affinity`).
+        sparsity: Sparsity fraction (see :func:`compute_affinity`).
+        alpha: Anisotropic diffusion parameter.
+        diffusion_time: Diffusion time parameter.
+
+    Returns:
+        :class:`GradientResult` with ``gradients``, ``lambdas``, and
+        ``affinity``.
+
+    Raises:
+        ValueError: If *matrix* is not 2-D square, not symmetric, or has
+            fewer than 2 ROIs.
+    """
+    if matrix.ndim != 2 or matrix.shape[0] != matrix.shape[1]:
+        raise ValueError(f"Expected square 2D matrix, got shape {matrix.shape}")
+    n = matrix.shape[0]
+    if n < 2:
+        raise ValueError(f"Need at least 2 ROIs, got {n}")
+
+    mat = np.array(matrix, dtype=np.float64)
+    mat_no_nan = np.nan_to_num(mat, nan=0.0)
+    if not np.allclose(mat_no_nan, mat_no_nan.T, atol=1e-8):
+        raise ValueError("Matrix must be symmetric")
+
+    n_components = min(n_components, n - 1)
+
+    aff = compute_affinity(mat, kernel=kernel, sparsity=sparsity)
+    gradients, lambdas = diffusion_map_embedding(
+        aff, n_components=n_components, alpha=alpha, diffusion_time=diffusion_time
+    )
+
+    return GradientResult(gradients=gradients, lambdas=lambdas, affinity=aff)
+
+
+def compute_gradients_from_files(
+    in_file: str | Path,
+    out_dir: str | Path | None = None,
+    n_components: int = 10,
+    kernel: AffinityKernel = "normalized_angle",
+    sparsity: float = 0.9,
+    alpha: float = 0.5,
+    diffusion_time: int = 0,
+) -> GradientOutputs:
+    """Load a correlation matrix from TSV and write gradient results to disk.
+
+    Thin wrapper around :func:`compute_gradients` that handles file I/O.
+
+    Args:
+        in_file: Path to TSV file containing a correlation matrix (as
+            produced by :func:`~rbc.core.metrics.timeseries.compute_timeseries`).
+        out_dir: Output directory.  Defaults to the parent directory of
+            *in_file*.
+        n_components: Number of gradient components.
+        kernel: Affinity kernel.
+        sparsity: Sparsity fraction.
+        alpha: Anisotropic diffusion parameter.
+        diffusion_time: Diffusion time parameter.
+
+    Returns:
+        :class:`GradientOutputs` with paths to the gradient and eigenvalue
+        TSV files.
+    """
+    in_file = Path(in_file)
+    matrix = np.loadtxt(in_file, delimiter="\t")
+
+    result = compute_gradients(
+        matrix,
+        n_components=n_components,
+        kernel=kernel,
+        sparsity=sparsity,
+        alpha=alpha,
+        diffusion_time=diffusion_time,
+    )
+
+    if out_dir is None:
+        out_dir = in_file.parent
+    out_dir = Path(out_dir)
+    out_dir.mkdir(parents=True, exist_ok=True)
+
+    stem = in_file.stem
+    grad_path = out_dir / f"{stem}_gradients.tsv"
+    lambda_path = out_dir / f"{stem}_lambdas.tsv"
+
+    np.savetxt(grad_path, result.gradients, delimiter="\t")
+    np.savetxt(lambda_path, result.lambdas.reshape(1, -1), delimiter="\t")
+
+    return GradientOutputs(gradients=grad_path, lambdas=lambda_path)
diff --git a/tests/unit/core/test_gradients.py b/tests/unit/core/test_gradients.py
new file mode 100644
index 00000000..a4f22b1b
--- /dev/null
+++ b/tests/unit/core/test_gradients.py
@@ -0,0 +1,323 @@
+"""Unit tests for rbc.core.metrics.gradients."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pytest
+
+from rbc.core.metrics.gradients import (
+    compute_affinity,
+    compute_gradients,
+    compute_gradients_from_files,
+    diffusion_map_embedding,
+)
+
+if TYPE_CHECKING:
+    from pathlib import Path
+
+
+def _make_two_block_corr(n_per_block: int = 20, rng_seed: int = 42) -> np.ndarray:
+    """Create a correlation matrix with two clearly separated blocks.
+
+    Within-block correlations are high (~0.8), between-block correlations
+    are low (~0.1).
+    """
+    n = 2 * n_per_block
+    corr = np.eye(n)
+    rng = np.random.default_rng(rng_seed)
+
+    for i in range(n):
+        for j in range(i + 1, n):
+            same_block = (i < n_per_block) == (j < n_per_block)
+            val = rng.uniform(0.7, 0.9) if same_block else rng.uniform(0.0, 0.2)
+            corr[i, j] = val
+            corr[j, i] = val
+    return corr
+
+
+class TestComputeAffinity:
+    """Tests for compute_affinity."""
+
+    def test_output_shape_and_symmetry(self) -> None:
+        """Output should be (n, n) and symmetric."""
+        rng = np.random.default_rng(0)
+        mat = rng.standard_normal((10, 10))
+        mat = (mat + mat.T) / 2
+        aff = compute_affinity(mat, sparsity=0.0)
+        assert aff.shape == (10, 10)
+        np.testing.assert_allclose(aff, aff.T)
+
+    def test_non_negative(self) -> None:
+        """All affinity values should be >= 0."""
+        rng = np.random.default_rng(1)
+        mat = rng.standard_normal((8, 8))
+        mat = (mat + mat.T) / 2
+        aff = compute_affinity(mat)
+        assert np.all(aff >= -1e-15)
+
+    def test_sparsity_reduces_input(self) -> None:
+        """With sparsity > 0, some entries in the input are zeroed before kernel."""
+        n = 20
+        rng = np.random.default_rng(2)
+        mat = rng.standard_normal((n, n))
+        mat = (mat + mat.T) / 2
+
+        aff_full = compute_affinity(mat, sparsity=0.0)
+        aff_sparse = compute_affinity(mat, sparsity=0.9)
+
+        # Sparse version should differ from full
+        assert not np.allclose(aff_full, aff_sparse)
+
+    def test_sparsity_zero_keeps_all(self) -> None:
+        """sparsity=0 should not zero out any off-diagonal connections."""
+        n = 10
+        rng = np.random.default_rng(3)
+        mat = rng.standard_normal((n, n))
+        mat = (mat + mat.T) / 2
+        aff = compute_affinity(mat, sparsity=0.0)
+        off_diag = aff[~np.eye(n, dtype=bool)]
+        assert np.all(off_diag > 0)
+
+    def test_kernel_cosine(self) -> None:
+        """Cosine kernel should produce non-negative values (after clamping)."""
+        rng = np.random.default_rng(4)
+        mat = rng.standard_normal((8, 8))
+        mat = (mat + mat.T) / 2
+        aff = compute_affinity(mat, kernel="cosine", sparsity=0.0)
+        assert np.all(aff >= -1e-15)
+        assert np.all(aff <= 1.0 + 1e-15)
+
+    def test_kernel_normalized_angle(self) -> None:
+        """Normalized angle kernel should produce values in [0, 1]."""
+        rng = np.random.default_rng(5)
+        mat = rng.standard_normal((8, 8))
+        mat = (mat + mat.T) / 2
+        aff = compute_affinity(mat, kernel="normalized_angle", sparsity=0.0)
+        assert np.all(aff >= -1e-15)
+        assert np.all(aff <= 1.0 + 1e-15)
+
+    def test_kernel_gaussian(self) -> None:
+        """Gaussian kernel should produce values in (0, 1]."""
+        rng = np.random.default_rng(6)
+        mat = rng.standard_normal((8, 8))
+        mat = (mat + mat.T) / 2
+        aff = compute_affinity(mat, kernel="gaussian", sparsity=0.0)
+        assert np.all(aff >= -1e-15)
+        assert np.all(aff <= 1.0 + 1e-15)
+
+    def test_diagonal_is_one(self) -> None:
+        """Diagonal should be 1.0 (self-similarity) for normalized_angle."""
+        rng = np.random.default_rng(7)
+        mat = rng.standard_normal((6, 6))
+        mat = (mat + mat.T) / 2
+        aff = compute_affinity(mat, sparsity=0.0)
+        np.testing.assert_allclose(np.diag(aff), 1.0)
+
+    def test_identity_input(self) -> None:
+        """Identity matrix: one-hot rows have cosine sim = 0 off-diagonal."""
+        mat = np.eye(5)
+        aff = compute_affinity(mat, sparsity=0.0)
+        off_diag = aff[0, 1]
+        assert 0.4 < off_diag < 0.6
+
+    def test_rejects_non_square(self) -> None:
+        """Should reject non-square input."""
+        with pytest.raises(ValueError, match="square"):
+            compute_affinity(np.zeros((3, 4)))
+
+    def test_rejects_invalid_sparsity(self) -> None:
+        """Should reject sparsity outside [0, 1)."""
+        mat = np.eye(3)
+        with pytest.raises(ValueError, match="sparsity"):
+            compute_affinity(mat, sparsity=1.0)
+        with pytest.raises(ValueError, match="sparsity"):
+            compute_affinity(mat, sparsity=-0.1)
+
+    def test_nan_handling(self) -> None:
+        """NaN values in the matrix should be treated as 0."""
+        mat = np.ones((4, 4))
+        mat[0, 1] = np.nan
+        mat[1, 0] = np.nan
+        aff = compute_affinity(mat, sparsity=0.0)
+        assert not np.any(np.isnan(aff))
+
+
+class TestDiffusionMapEmbedding:
+    """Tests for diffusion_map_embedding."""
+
+    def _make_affinity(self, n: int = 20, seed: int = 0) -> np.ndarray:
+        """Build a simple positive affinity matrix."""
+        rng = np.random.default_rng(seed)
+        mat = rng.standard_normal((n, n))
+        mat = (mat + mat.T) / 2
+        return compute_affinity(mat, sparsity=0.5)
+
+    def test_output_shapes(self) -> None:
+        """Gradients should be (n, n_components), lambdas (n_components,)."""
+        aff = self._make_affinity(20)
+        grads, lams = diffusion_map_embedding(aff, n_components=5)
+        assert grads.shape == (20, 5)
+        assert lams.shape == (5,)
+
+    def test_eigenvalues_descending(self) -> None:
+        """Eigenvalues should be sorted in descending order."""
+        aff = self._make_affinity(20)
+        _, lams = diffusion_map_embedding(aff, n_components=5)
+        assert np.all(np.diff(lams) <= 1e-10)
+
+    def test_deterministic(self) -> None:
+        """Same input should produce identical output."""
+        aff = self._make_affinity(15)
+        g1, l1 = diffusion_map_embedding(aff, n_components=3)
+        g2, l2 = diffusion_map_embedding(aff, n_components=3)
+        np.testing.assert_array_equal(g1, g2)
+        np.testing.assert_array_equal(l1, l2)
+
+    def test_block_separable(self) -> None:
+        """Near-block-diagonal affinity with weak cross-links should separate."""
+        n = 10
+        aff = np.full((2 * n, 2 * n), 0.01)
+        aff[:n, :n] = 0.5
+        aff[n:, n:] = 0.5
+        np.fill_diagonal(aff, 1.0)
+
+        grads, _ = diffusion_map_embedding(aff, n_components=2)
+        g1 = grads[:, 0]
+        block_a = g1[:n]
+        block_b = g1[n:]
+        assert np.sign(block_a.mean()) != np.sign(block_b.mean())
+
+    def test_n_components_respected(self) -> None:
+        """Should return exactly n_components columns."""
+        aff = self._make_affinity(20)
+        for nc in [2, 5, 10]:
+            grads, lams = diffusion_map_embedding(aff, n_components=nc)
+            assert grads.shape[1] == nc
+            assert lams.shape[0] == nc
+
+    def test_alpha_zero(self) -> None:
+        """alpha=0 should skip the anisotropic normalization step."""
+        aff = self._make_affinity(15)
+        g0, _ = diffusion_map_embedding(aff, n_components=3, alpha=0.0)
+        g1, _ = diffusion_map_embedding(aff, n_components=3, alpha=1.0)
+        assert not np.allclose(g0, g1)
+
+    def test_diffusion_time_positive(self) -> None:
+        """Positive diffusion_time should scale eigenvalues as w^t."""
+        aff = self._make_affinity(15)
+        g0, _l0 = diffusion_map_embedding(aff, n_components=3, diffusion_time=0)
+        g1, _l1 = diffusion_map_embedding(aff, n_components=3, diffusion_time=1)
+        # Different scaling → different gradient magnitudes
+        assert not np.allclose(g0, g1)
+
+    def test_rejects_non_square(self) -> None:
+        """Should reject non-square affinity."""
+        with pytest.raises(ValueError, match="square"):
+            diffusion_map_embedding(np.zeros((3, 4)), n_components=1)
+
+    def test_rejects_too_many_components(self) -> None:
+        """Should reject n_components >= n_rois."""
+        aff = self._make_affinity(5)
+        with pytest.raises(ValueError, match="n_components"):
+            diffusion_map_embedding(aff, n_components=5)
+
+
+class TestComputeGradients:
+    """Tests for the high-level compute_gradients function."""
+
+    def test_output_types_and_shapes(self) -> None:
+        """Result should be a GradientResult with correct shapes."""
+        rng = np.random.default_rng(10)
+        mat = rng.standard_normal((15, 15))
+        mat = (mat + mat.T) / 2
+        result = compute_gradients(mat, n_components=3)
+
+        assert result.gradients.shape == (15, 3)
+        assert result.lambdas.shape == (3,)
+        assert result.affinity.shape == (15, 15)
+
+    def test_two_block_first_gradient_separates(self) -> None:
+        """Two-block correlation matrix: first gradient separates blocks."""
+        corr = _make_two_block_corr(n_per_block=20)
+        result = compute_gradients(corr, n_components=5, sparsity=0.8)
+
+        g1 = result.gradients[:, 0]
+        block_a = g1[:20]
+        block_b = g1[20:]
+
+        assert np.sign(block_a.mean()) != np.sign(block_b.mean())
+        gap = abs(block_a.mean() - block_b.mean())
+        within = max(block_a.std(), block_b.std())
+        assert gap > 2 * within
+
+    def test_rejects_non_symmetric(self) -> None:
+        """Should reject non-symmetric matrix."""
+        mat = np.array([[1.0, 0.5], [0.3, 1.0]])
+        with pytest.raises(ValueError, match="symmetric"):
+            compute_gradients(mat)
+
+    def test_rejects_too_few_rois(self) -> None:
+        """Should reject matrix with < 2 ROIs."""
+        with pytest.raises(ValueError, match="2 ROIs"):
+            compute_gradients(np.array([[1.0]]))
+
+    def test_clamps_n_components(self) -> None:
+        """n_components > n-1 should be silently clamped."""
+        mat = np.eye(5)
+        result = compute_gradients(mat, n_components=100)
+        assert result.gradients.shape[1] == 4  # n - 1
+
+    def test_correlation_matrix_input(self) -> None:
+        """Standard correlation matrix should work end-to-end."""
+        rng = np.random.default_rng(20)
+        ts = rng.standard_normal((10, 50))
+        corr = np.corrcoef(ts)
+        result = compute_gradients(corr, n_components=5)
+        assert result.gradients.shape == (10, 5)
+        assert not np.any(np.isnan(result.gradients))
+
+
+class TestComputeGradientsFromFiles:
+    """Tests for the file I/O wrapper."""
+
+    def test_round_trip(self, tmp_path: Path) -> None:
+        """Write TSV, compute, load results."""
+        corr = _make_two_block_corr(n_per_block=10)
+        in_file = tmp_path / "sub-01_correlation_matrix.tsv"
+        np.savetxt(in_file, corr, delimiter="\t")
+
+        result = compute_gradients_from_files(in_file, n_components=3)
+
+        assert result.gradients.exists()
+        assert result.lambdas.exists()
+
+        grads = np.loadtxt(result.gradients, delimiter="\t")
+        assert grads.shape == (20, 3)
+
+        lams = np.loadtxt(result.lambdas, delimiter="\t")
+        assert lams.size == 3
+
+    def test_output_naming(self, tmp_path: Path) -> None:
+        """Output files should follow <stem>_gradients.tsv convention."""
+        corr = _make_two_block_corr(n_per_block=5)
+        in_file = tmp_path / "sub-01_correlation_matrix.tsv"
+        np.savetxt(in_file, corr, delimiter="\t")
+
+        result = compute_gradients_from_files(in_file, n_components=2)
+
+        assert result.gradients.name == "sub-01_correlation_matrix_gradients.tsv"
+        assert result.lambdas.name == "sub-01_correlation_matrix_lambdas.tsv"
+
+    def test_custom_out_dir(self, tmp_path: Path) -> None:
+        """Should write to a custom output directory."""
+        corr = _make_two_block_corr(n_per_block=5)
+        in_file = tmp_path / "corr.tsv"
+        out_dir = tmp_path / "output"
+        np.savetxt(in_file, corr, delimiter="\t")
+
+        result = compute_gradients_from_files(in_file, out_dir=out_dir, n_components=2)
+
+        assert result.gradients.parent == out_dir
+        assert result.lambdas.parent == out_dir