The Variance Brain Foundation Models Forgot: Third-Order Statistics Predict Cognition Where Billion-Parameter Models Fail

TL;DR: Three state-of-the-art fMRI foundation models — billion-parameter-class Transformers pretrained on brain scans — predict individual cognition worse than plain linear regression on the ~80K numbers of a functional connectivity matrix, and the bigger BrainLM does worse than the smaller one. The authors diagnose why: masked-reconstruction pretraining is dominated by the high-variance, non-cognitive part of the fMRI signal, so it preserves second-order covariance passably but destroys the third-order co-skewness where cognitive signal lives. A GPU-free linear pipeline that computes connectivity inside the co-skewness tensor's Tucker subspace then beats everything, and a cumulant-aware finetuning loss restores the foundation model to the classical baseline — pinning the failure on the objective, not the architecture or scale.

The bet that brain foundation models are making

Resting-state fMRI gives you a matrix \mathbf{Z} \in \mathbb{R}^{P \times T}: activity for P \approx 400 brain regions over T timepoints. The classical way to predict an individual's traits from this is embarrassingly simple: compute functional connectivity (FC), the P \times P matrix of pairwise Pearson correlations between regions, and run kernel ridge regression (KRR) from its ~80K entries to a cognition score. This has been the workhorse of individual-differences neuroimaging for a decade.

Brain foundation models (BFMs) bet that this throws away too much. BrainLM is a ViT-MAE that reconstructs masked fMRI patches; Brain-JEPA predicts masked-patch representations in latent space; BrainMass does masked modeling on FC matrices directly. The reasoning is the standard foundation-model story: a deep Transformer with self-supervised pretraining on thousands of scans should learn richer structure than pairwise correlations — Transformers can represent arbitrary higher-order interactions, and FC is just the second moment.

Strikingly, BFM papers almost never benchmark against the FC baseline under a controlled protocol. This paper does, using the replicated nested cross-validation setup of Ooi et al. (2022): 20 repetitions × 10 folds, family-aware splits on HCP, decompositions fit on training folds only, 200 fold-level evaluations per method.

The result: it's not close, and scale makes it worse

Every readout the authors could extract from the pretrained models — CLS embeddings, patch tokens, embedding-similarity matrices, FC computed on the reconstructed timeseries — loses to FC computed on the model's own input:

On HCP, BrainLM-650M's embedding scores r = 0.003 and its reconstruction-FC scores r = -0.028 — noise — while FC on the same input reaches 0.393. And the 650M model underperforms the 111M model on both datasets and both readouts: inverse scaling, the same pathology McKenzie et al. documented in LLMs. (BrainMass, which ingests FC directly rather than raw timeseries, is a partial exception at r = 0.452 vs. its input-FC baseline of 0.504 — still below, but not catastrophic, which itself supports the diagnosis to come.)

An appendix result makes the failure vivid: raw FC identifies which of 955 subjects a scan segment came from with 87% accuracy; all three BFM embeddings sit at chance (~0.3%), even when FC is PCA-compressed to the same 768 dimensions (which retains 84%). The models haven't lost cognition specifically — they've lost individual identity wholesale.

The diagnosis: a variance allocation problem

Here is the "aha." Think of the joint distribution of regional activity through its cumulants \kappa_n: \kappa_2(z_i, z_j) = \mathbb{E}[z_i z_j] is the covariance (i.e., FC), and \kappa_3(z_i, z_j, z_k) = \mathbb{E}[z_i z_j z_k] is the three-way co-skewness — for a Gaussian it would be exactly zero, so any nonzero \kappa_3 is non-Gaussian structure invisible to FC. A growing neuroscience literature (Santoro et al., Luppi et al., Gatica et al.) says cognition-relevant signal concentrates precisely in these higher-order, synergistic interactions.

Now consider what an MSE reconstruction loss sees. fMRI variance is dominated by cardiac and respiratory cycles, vasomotor oscillations, head motion, scanner drift — big, structured, and cognitively irrelevant. Gradient descent on reconstruction error is an allocator of model capacity proportional to variance, so nearly all the gradient signal pushes the model to fit that dominant second-order surface. It's like training an audio autoencoder on recordings where 95% of the power is mains hum: the loss richly rewards modeling the hum and barely notices whether the speech survives.

The per-cumulant analysis (Section 4.2) is the paper's most direct evidence: BrainLM's reconstructions preserve \kappa_2 at 29–47% relative error, but the \kappa_3 NMSE is 2.6–6.2× the input's signal energy — the model doesn't just fail to capture co-skewness, it actively hallucinates wrong third-order structure. And the 650M model is uniformly worse than the 111M at both orders, which explains the inverse scaling: more capacity spent fitting the wrong objective better.

The fix without a model: FC in the co-skewness basis

If cognition lives in \kappa_3, why not just regress from the co-skewness tensor \mathcal{S} \in \mathbb{R}^{P \times P \times P}, S_{ijk} = \frac{1}{T}\sum_t z_i(t) z_j(t) z_k(t)? Because it has ~13 million unique entries for P = 424, estimated from ~200 timepoints — each entry is drowning in sampling noise. Appendix A confirms: raw co-skewness features under KRR score r \approx 0.23, below plain FC.

The insight is that the tensor's geometry is stable even when its entries are not. Tucker decomposition (higher-order SVD) factors \mathcal{S} \approx \mathcal{G} \times_1 \mathbf{U} \times_2 \mathbf{U} \times_3 \mathbf{U}, where \mathbf{U} \in \mathbb{R}^{P \times R} spans the spatial subspace that best captures third-order structure — the exact analogue of what PCA does for variance. The pipeline: fit \mathbf{U} on training subjects, project the timeseries to \tilde{\mathbf{Z}} = \mathbf{U}^\top \mathbf{Z}, compute FC inside that subspace, run KRR. Second-order features in a third-order-informed basis.

This beats FC-full on all four dataset × parcellation cells: +0.065 on AOMIC AAL (d = 0.90), +0.060 on AOMIC Schaefer (d = 0.82), +0.038 and +0.028 on HCP. On HCP Schaefer-400 it reaches r = 0.571 \pm 0.089, exceeding the best published rs-FC numbers on the same cognition factor (Kong et al. 2023: 0.520) — with no pretraining and no GPU. The crucial control is PCA at matched dimensionality: same feature count, same pipeline, only the basis criterion differs (\kappa_2 vs. \kappa_3). Tucker wins with d = 1.57 and 1.21 on the AAL cells; PCA at K = 117 actually underperforms FC-full. Variance magnitude is not behavioral relevance. Ablations show the gain is spatial not temporal (matched temporal reduction is a wash), isn't a resolution artifact (Tucker at R = 80, a 25× compression, beats Schaefer-400's full 79,800 features), isn't post-hoc rank selection (fully nested CV preserves the win everywhere), and isn't motion loading (Tucker components are no more motion-correlated than PCA components).

The fix with the model: change the loss, keep the architecture

The cleanest way to prove the objective is the bottleneck is to change only the objective. The authors finetune BrainLM end-to-end with a dual-moment loss: a Log-Cholesky distance (the natural Riemannian metric on SPD matrices) between input and reconstructed FC, computed inside the frozen rank-80 Tucker subspace, lightly mixed (\lambda = 10^{-3}) with ambient FC distance. No cognition labels are ever seen.

The dual-moment finetune recovers the raw-FC ceiling essentially exactly, at both scales, in both directions of cross-dataset transfer — and the inverse-scaling gap between 111M and 650M disappears (both hit r = 0.392 on HCP to three decimals). A direct \kappa_3-MSE surrogate confirms the gain comes from third-order structure specifically, though its tiny-magnitude gradients train poorly on their own.

One honest wrinkle the authors foreground: only the reconstruction readout recovers; the CLS embedding stays near zero under every FC-preservation loss, and even direct supervised regression to cognition barely moves it (r = 0.137). The latent geometry set during pretraining does not reorganize when the loss reaches it only through the decoder. So the finetuned BFM matches, but never exceeds, raw FC — the result is diagnostic, not a better model. Coupling the \kappa_3 target directly to the latent is the open problem they name.

How much to believe

The effect sizes are the strength; the p-values are the honest weakness. Under the conservative Nadeau–Bengio corrected test with FDR, Tucker-over-FC-full is significant only in the two AOMIC cells (q = 0.034); the HCP cells (d = 0.49–0.62) survive on effect size and replicated direction, and the authors compute they'd need ~2,000 subjects for HCP-Schaefer significance. The matched-dimensionality Tucker-vs-PCA claim — the one that isolates the third-order mechanism — is significant in 3 of 4 cells. Appendix B lays all this out with unusual candor.

Other caveats: only two cognitive targets (an HCP cognition factor, an AOMIC IQ composite), both composites of healthy, Western, well-educated cohorts; clinical and developmental phenotypes may have different variance structures. Per-variable analysis shows Tucker helps on 17 of 57 HCP behavioral measures and FC-full on 20 — the composite-level win is not universal across phenotypes. The BFM-side evaluation is thorough (four readouts per model, plus supervised finetuning probes), but one can't fully exclude that some other extraction or heavier task-specific adaptation would rescue the pretrained models. And the paper stops at \kappa_3; whether \kappa_4+ adds anything is open (estimation cost grows fast).

Why this matters beyond fMRI

The general lesson is sharper than "brain foundation models are overhyped." Reconstruction and joint-embedding objectives allocate capacity by variance — Balestriero & LeCun made this argument in vision — and in any domain where dominant variance is not task-relevant (physiological signals, climate, single-cell dynamics, finance), scaling a self-supervised model can be actively counterproductive, and per-cumulant preservation analysis is a concrete diagnostic for localizing what a pretrained model discarded. The prescriptions are actionable: always benchmark foundation models against strong classical baselines under a controlled protocol; when you know the task-relevant statistical order, either extract it explicitly (FC-Tucker is a drop-in, CPU-only default for brain-behavior prediction) or bake it into an auxiliary loss.

If you read two sections, make it Section 4.3 (the Tucker-vs-PCA controlled comparison that isolates the third-order mechanism) and Appendix A (the elegant negative result showing the co-skewness tensor's entries are useless while its basis is gold — the single observation the whole method hinges on). Appendix B is worth five minutes for a model of how to report CV statistics without overclaiming.