Papers

Interval localization in multiple change-point detection aims to identify intervals containing true change-points while controlling false discoveries. We formulate this task as a multiple testing problem over detector-reported intervals and develop an algorithm-agnostic and distribution-free framework that controls the interval-level FDR, the expected proportion of retained intervals containing no true change-points. The proposed permutation-based order-preserving interval selection (POIS) uses order-preserving sample splitting to separate detection from permutation testing. To improve power, POIS+ reuses the detection-side evidence and introduces detection-conditional $e$-values as unnormalized weights for permutation $p$-values. Both procedures achieve finite-sample FDR control. Under the stated conditions, POIS attains asymptotically full power while POIS+ has asymptotically no lower power than POIS. Simulation results show that POIS maintains high power and POIS+ further improves power with empirical FDR below the nominal level for both procedures. Two real data applications demonstrate the practical effectiveness of the proposed methods.

This article generalizes the de-initialization framework, proposed by (Roberts and Rosenthal, 2001) for total variation distance, to Wasserstein distances. In essence, de-initialization captures the phenomenon that in the presence of certain structural features, the convergence behaviour of a Markov chain can be reduced to that of a "simpler" stochastic process. Our results allow to measure this in terms of a general Wasserstein distance. We apply them to analyse the Wasserstein convergence behaviour of Gibbs sampling, linchpin variable sampling and slice sampling.

Generative artificial intelligence (GenAI) is reshaping scholarly research faster than journals have developed stable norms for its use. This article presents an edited thematic account of a 2026 Joint Statistical Meetings panel that brought together editorial perspectives from mathematical statistics, data science, biomedical statistics, and general statistical scholarship. The discussion examines journal policies, disclosure, authorship and research integrity, peer-review confidentiality, researcher training, editorial workload, access, and possible future models of scholarly publishing. Panelists shared commitments to human accountability, the protection of confidential submissions, and disclosure of consequential assistance, while offering different recommendations on assistance with research ideas and proofs, disclosure requirements, automated review, and policy enforcement. By distinguishing shared principles from unresolved implementation questions, the article clarifies the choices facing statistical publishing and outlines an agenda for evaluating policies and practices as GenAI evolves. The account seeks to foster continued discussion of GenAI in scientific communication and encourage statistical organizations to develop more robust operational standards.

Evaluating "neutral launches" (e.g., infrastructure upgrades) using traditional confidence interval overlap is flawed: it is dangerously permissive with scarce data and excessively restrictive with abundant data. To resolve this, this paper introduces Expected Bayesian Loss (EBL), a continuous metric that quantifies both the probability and expected severity of metric degradation. Computable directly from standard frequentist estimates, EBL explicitly penalizes empirical noise and high-variance experiments. Validated against expert decisions, EBL provides experimentation platforms with a rigorous, tunable guardrail that aligns statistical safety with institutional risk appetite.

Wasserstein gradient flow extends gradient descent to probability measures. Its Hessian-guided perturbed variant (PWGF) adds Gaussian perturbations to escape saddle points in nonconvex problems. We investigate when its approximation by finitely many interacting particles remains accurate over growing time horizons. Our analysis retains the curvature accumulated along the population-driven reference path: negative curvature can amplify approximation errors, while subsequent positive curvature can damp their influence. This captures favorable scenarios in which temporary instability is compatible with accurate tracking over growing horizons. Under regularity assumptions and a prescribed common perturbation schedule, we prove particle and objective-value tracking bounds on a high-probability event for reference paths satisfying explicit conditions on accumulated curvature. To handle state-dependent Gaussian jumps, we construct a population-first coupling that preserves the reference particles' conditional independence and reduces jump errors to covariance comparison. We verify the conditions in a variance-plus-cosine model, where curvature recovery yields a growing-horizon tracking guarantee. We also establish local attraction, transverse descent, and positive second variation in two regions of a regularized matrix-factorization model, motivating a positive-negative-positive curvature pattern.

Bayesian additive regression trees (BART) can require many splits to approximate boundaries misaligned with the predictor axes. RoBART assigns each tree a rotation shared by all internal nodes, retaining axis-aligned splits in rotated coordinates and constant leaves. We jointly propose a Givens rotation sequence and cutpoints on the resulting grid by Metropolis-Hastings and establish reversibility with respect to the conditional posterior with leaf means integrated out. For additive functions with component-specific rotations and anisotropic Hölder smoothness, we prove posterior contraction in empirical $L_2$ distance and for the noise standard deviation. Under the stated prior, design, and grid conditions, with fixed numbers of predictors, trees, and components and no more components than trees, the rate is a sum of componentwise rates determined by smoothness and the number of rotated coordinates used. We also establish a posterior contraction lower bound showing that there exist functions for which RoBART adapts to the intrinsic dimension but axis-aligned BART does not.

Links are disjoint unions of circles smoothly embedded in $S^3$. We use reinforcement learning and Bayesian optimisation to obtain new upper bounds on several link invariants that are not known to be algorithmically computable: the slice genus and the unknotting number for links, and the strong slice genus for algebraically split links. We also compute lower bounds using known invariants. Combining the upper and lower bounds, we obtain new exact values in many cases. Our unknotting agents can reproduce the non-additivity of the unknotting number for several counterexamples due to Brittenham and Hermiller, in some cases finding new unknotting trajectories.

Switching stochastic differential equations (SSDEs) describe continuous-time dynamics whose parameters switch according to a latent regime process that follows a continuous-time Markov chain (CTMC). By allowing dynamics to change between regimes, SSDEs represent heterogeneous system behavior and have been applied across diverse fields. However, Bayesian inference for SSDEs remains difficult, and existing SSDE inference methods have limited applicability, with restrictions such as noise-free observations, univariate states, linear drift, or state-independent diffusion. In this study, we propose an approximate Markov chain Monte Carlo sampler for SSDEs using uniformization and factorized neural likelihood estimation (FNLE), a simulation-based inference method. Uniformization provides an exact representation of the CTMC but requires SDE transition densities over arbitrary time intervals. We approximate these densities by training a time-conditioned FNLE model. The resulting sampler is broadly applicable to SSDEs without requiring analytically tractable transition densities. In synthetic-data experiments, our method recovered regime paths and parameters for three SSDE models for which previous methods have limited applicability. We also applied our method to a real dataset and detected a regime transition.

We analyze the first-order deterministic diffusion sampler of Karras et al. (2022), termed EDM, in 2-Wasserstein distance by separating two sources of error: local discretization error and its amplification by subsequent learned steps. We prove that local error admits a universal bound: for any data distribution with finite second moment, the one-step discretization error is quadratic in the step size, with an explicit constant that does not depend on the data distribution. Error propagation, in contrast, depends on the learned network. At high noise levels, we exploit the network parametrization of EDM to derive an explicit contraction criterion. At low noise levels, we measure propagation through the amplification realized on the distributions transported by the sampler; this realized amplification can be arbitrarily smaller than the worst-case Lipschitz constant. This analysis yields an $O(e^{\Lambda_K}/K)$ global discretization error for $K$ sampling steps, where $\Lambda_K$ is the low-noise log-amplification. Experiments on a one-dimensional Gaussian mixture show how measured amplification accounts for slower error decay on finite sampling grids. Diagnostics on a pretrained CIFAR-10 model illustrate related stability mechanisms without certifying the global assumptions.

Split Gibbs sampling (SGS) is a popular framework for posterior sampling in Bayesian imaging inverse problems. It decouples a Gaussian data-fidelity term from a complex prior through an auxiliary variable, so the data variable is updated exactly and only the prior-side conditional is hard to sample. Existing samplers treat this conditional in one of two ways. Plug-and-play SGS runs a multi-step diffusion denoiser at every iteration, which is expensive and lacks non-asymptotic guarantees. Langevin-within-SGS takes cheap overdamped Langevin steps but needs many iterations. We propose RED-KLwSGS, which keeps the exact Gaussian update for the data variable and updates the auxiliary variable with underdamped (kinetic) Langevin diffusions driven by a one-shot denoising score, at the same per-iteration cost as Langevin-within-SGS. We prove non-asymptotic Wasserstein-2 convergence in continuous and discrete time for strongly log-concave priors. We also introduce Joint-RED-KLwSGS, which applies kinetic Langevin diffusions to both variables. Experiments with Denoising diffusion probabilistic models as diffusion priors on FFHQ and ImageNet datasets show faster convergence and high-quality image reconstruction.

Bayesian optimization (BO) is widely used as a standard approach for expensive black-box optimization. However, BO algorithms often involve parameters that must be specified in advance, and their performance can strongly depend on these choices. We propose a framework for optimizing such parameters using sample paths drawn from a Gaussian process (GP) inferred from the information available at the start of BO. We use cumulative regret as the performance metric for a BO algorithm. By running the BO algorithm on the generated sample paths, we obtain an empirical estimate of its expected cumulative regret for a given parameter configuration. Optimizing this estimate allows us to identify parameter configurations that, given the currently available information, are expected to achieve low cumulative regret. Since this parameter optimization is itself a black-box optimization problem, we employ another BO procedure to solve it, which we refer to as outer BO. Through experiments, we demonstrate that the proposed framework can effectively select parameter configurations that achieve strong performance among a range of candidate configurations.

Split conformal prediction uses prediction errors on held-out (calibration) data to determine how wide the prediction intervals should be. It guarantees distribution-free finite-sample coverage when these errors and the error at the target site are exchangeable. This assumption may fail under spatial dependence and nonrandom sampling geometry. Existing spatial methods use fitting residuals to remove the predictable part of spatial variation from calibration and target errors. However, the spatial variation that only the calibration residuals can predict remains in both the target and calibration errors, reducing the efficiency and stability of the interval. We address this by additionally conditioning on the calibration residuals sequentially, which scales to large networks through nearest-neighbour approximations. Under a correct working covariance and an elliptical residual law, the resulting interval has exact finite-sample coverage under any spatial design, and under further conditions it is asymptotically oracle efficient. We also bound coverage loss under covariance misspecification and develop a diagnostic that identifies regions at risk of undercoverage. In simulated data, our method produces narrower and more stable intervals than global and localized state-of-the-art alternatives. In a national PM2.5 application, it produces narrower intervals within the network and identifies regions at risk of coverage failure.

Policy learning in digital experimentation faces three challenges: weak signal-to-noise ratios, rich covariate spaces, and massive data volumes. We formalize this regime by modeling treatment-effect estimates from increasingly fine covariate partitions as Gaussian observations with bounded signal-to-noise ratios. We establish that, in general, the optimal treatment policy is not learnable in this setting. Even learning the optimal policy value suffers from impractically slow rates. However, when treatment effects vary smoothly, we derive minimax-adaptive policies based on linear smoothers that achieve vanishing welfare regret. We demonstrate the practical value of our framework by applying it to large-scale real-world experiments at Netflix, showing that personalized linear-smoothing policies can dominate unpersonalized policies even in this challenging empirical setting.

Offline evals increasingly guide changes to LLM-powered products, but showing that an eval agrees with human judgment does not tell a team which changes it can ship without an A/B test. That decision depends on whether eval treatment effects predict online treatment effects across a class of interventions. We call this requirement trial-level surrogacy and show that it is logically independent of unit-level (Prentice) surrogacy. In historical launch logs, sampling noise in both effect estimates attenuates the observed relationship. Under a bivariate normal working model, subtracting the reported within-intervention variances recovers the between-intervention covariance, and the corrected moments give a closed-form ship/test/kill rule that bounds the probability of shipping a change with a non-positive online effect. LLM judges add a separate distortion. When the judge is equally accurate on both arms, eval effects are only rescaled, but arm-dependent accuracy can reverse their sign. Once deployed, the eval-to-outcome relationship can no longer be observed in the ship and kill regions, so drift there becomes undetectable without randomized audit experiments. On 183 interventions from the Upworthy Research Archive, the pooled eval-to-outcome relationship is weak despite precise measurement on both sides, and the illustrative gate sends every intervention to an A/B test.

Hierarchical Archimedean copulas (HACs) admit redundant tree representations: when adjacent internal nodes have equal parameters, contracting the connecting edge leaves the copula unchanged. This creates a structural overparametrization. For homogeneous one-parameter families, we remove this redundancy and identify the resulting reduced model class with the ultrametric fan. We treat the Clayton, Frank, Ali-Mikhail-Haq, Gumbel and Joe families. Using the negative logit of normalized Kendall's tau as a common node-height coordinate, reduced hierarchies form polyhedral strata, binary hierarchies index maximal cones, and multifurcating hierarchies lie on common contraction faces. The induced cophenetic metric generates the same topology as uniform convergence of copulas. We also show that informative pairwise dependence summaries act componentwise on cophenetic coordinates and hence provide global coordinates for the reduced model. This representation yields a global estimator obtained by inverting estimated pairwise summaries and projecting onto the fan, without preselecting a hierarchy. We prove consistency and convergence-rate transfer. At binary truths, the unthresholded estimator recovers the hierarchy with probability tending to one; at multifurcating truths, it asymptotically selects a binary refinement, while vanishing thresholding consistently recovers the reduced hierarchy. Under a joint central limit theorem for the pairwise summaries, the projected estimator has a root-$n$ limit given by Euclidean projection of a Gaussian vector onto the tangent fan. The limit is Gaussian at binary hierarchies and generally non-Gaussian at multifurcating ones. For empirical Kendall's tau, we verify the required regularity conditions and obtain the corresponding singular limit.

We show upper and lower bounds on the regret of $m$-set adversarial bandits for different utilities (winner reward or sum of rewards) and feedback models (winner index, winner reward, sum of rewards, and their combinations). By comparing to standard bounds for combinatorial and MNL bandits, our results reveal how subtle changes in the setting can have a dramatic impact on the learning rates. Our main technical contributions are the information-theoretic lower bounds on the regret. Experiments on synthetic data confirm our theoretical analyses.

Gaussian processes (GPs) are widely used across machine learning, spatial statistics, time-series analysis, optimization, Bayesian statistics, and scientific applications. A central component of a GP model is its kernel, which is typically specified through a parametric family. Among the most widely used choices is the radial basis function (RBF), also known as the squared exponential or Gaussian kernel, owing to its simple form, smoothness, and flexibility. In practice, the kernel parameters are routinely estimated by the maximum likelihood estimators (MLEs), as implemented by standard GP software. Despite this widespread use, the asymptotic behavior of the MLEs remains poorly understood under fixed-domain asymptotics, even for the RBF kernel. The main difficulty arises from the increasingly strong dependence among densely sampled observations and the nonlinear dependence of the covariance matrix on the kernel parameters. In this paper, we address this gap by providing, to the best of our knowledge, the first complete asymptotic characterization of the joint MLE of the spatial variance, lengthscale, and nugget variance under fixed-domain asymptotics. We establish consistency, derive convergence rates for all three parameters, prove joint asymptotic normality, and show that these rates are minimax optimal.

Calibration is an essential requirement for probabilistic predictions to be useful for decision making. While state-of-the-art prediction methods often yield miscalibrated predictive distributions, several post-hoc recalibration schemes have been proposed to generate calibrated predictions. However, popular recalibration schemes can conceal miscalibration in specific regions of the outcome space. Since particular outcomes, such as extreme events, often matter most for decision making, probabilistic predictions should be calibrated when evaluation is restricted to these outcomes. Hence, in this paper, we introduce outcome-conditional recalibration, a post-hoc method to recalibrate probabilistic predictions on user-defined regions of the outcome space. The method is simple, easy to implement, and can be applied to arbitrary predictive distributions. It works by applying the quantile recalibration approach of Kuleshov et al. (2018) to forecast conditional distributions, before rescaling these conditional distributions so that forecast event probabilities match empirical occurrence frequencies. This produces valid and continuous predictive distributions that are calibrated within each region of interest. Across regression benchmarks, we demonstrate that existing recalibration schemes do not necessarily yield calibrated predictions when interest is on particular outcomes, and that our approach improves outcome-conditional calibration relative to existing conditional and unconditional recalibration methods, while retaining competitive calibration overall. In an application to day-ahead electricity price forecasting, the approach substantially improves calibration when predicting negative prices, at negligible cost to forecast accuracy.

Sampling transitions between metastable states is a central problem in dynamical systems theory and molecular dynamics in particular. A key challenge is the existence of high free-energy barriers that separate the states, making transitions extremely rare. Recent machine learning-based methods cast transition path sampling (TPS) as an optimal stochastic control (OSC) problem over a fixed time horizon, and parameterize the drift bias via a neural network trained by simulation-in-the-loop, requiring repeated biased rollouts. To address computational and performance guarantee issues of these models, we propose a new approach for the problem based on Koopman operators. Because Koopman operators are linear, their leading eigenfunctions reveal the metastable sets and provide an estimate of the committor function with no transition path information required. Furthermore, we formulate TPS as an OSC problem up to an exit time. Our time horizon is the first hitting time of the target set, and our running cost penalizes time spent in nonreactive regions by encoding the estimated committor function. We derive the optimal controller in closed form and approximate it in a reproducing kernel Hilbert space (RKHS). This reduces the problem of constructing the optimal controller to solving a single equality-constrained quadratic program, whose solution can be characterized by a linear Karush-Kuhn-Tucker (KKT) system. On the two-channel double well and alanine dipeptide, our controller increases the fraction of trajectories reaching the target from 0% to 99.8% within 1000 steps, and from 0% to 93% within 1ps, respectively.

A theoretical understanding of multi-head self-attention is fundamental to the study of modern neural networks. Using random matrix theory, we establish Gaussian equivalence for multi-head self-attention: replacing softmax attention with rescaled scores plus Gaussian noise preserves the limiting spectral law of the centered output. This equivalence also covers value and output projections that depend on the keys. The resulting laws separate the effects of head allocation and projection widths, and distinguish spectrum-preserving across-head sharing from within-head key–value dependence.

In this paper, we establish two nonasymptotic Berry–Esseen bounds over convex sets for the Gaussian approximation of multivariate nonlinear statistics. The statistics of interest can be written as a sum of independent centered random vectors plus a remainder that may depend on all observations. The first bound retains the classical factor $d^{1/4}$ in the contribution of the independent sum, where $d$ is the dimension, while controlling the remainder through its size and its sensitivity to replacing a single observation. The second bound expresses the contribution of the independent sum in terms of fourth moments and can allow the dimension to grow faster with the sample size. For sums of independent random vectors, it removes the logarithmic factor from an existing fourth moment bound without imposing additional moment assumptions. As applications, we apply these results to Polyak–Ruppert averaging for nonsmooth stochastic approximation, temporal difference learning with linear function approximation, and multivariate $U$-statistics. The resulting bounds provide explicit Gaussian approximation errors and sufficient conditions under which these errors converge to zero as the dimension grows with the sample size.

Mutual information (MI) provides an objective for suppressing or encouraging statistical dependence in implicit generative models. However, direct MI evaluation is challenging in implicit models due to typically intractable densities. A remedy is estimating the generator gradient from the difference between conditional and marginal scores. This score difference can, in turn, be estimated by differentiating a log density ratio learned through classification. This construction nevertheless faces two difficulties: (i) singular distributions need not admit the required score functions, and (ii) poor overlap can hinder density-ratio estimation. We therefore introduce Spread Mutual Information (SMI), a weighted integral of MI across noise levels obtained by applying a common spreading kernel to the generated variable. Gaussian spreading yields smooth, strictly positive conditional and marginal densities, extending the gradient construction to distributions that may originally be singular. Across a variaty of experiments, SMI consistently achieves effective dependence control among MI-based methods and remains competitive with established task-specific approaches.

The proposed approach relies on high-frequency multivariate time series collected at the Seine aval (SIAAP 1 ) treatment plant, including pH, temperature, conductivity, and total suspended solids, as well as exogenous variables related to precipitation and measurements from an upstream plant. The objective is to forecast wastewater quality over a 24-hour horizon. A multivariate Long Short-Term Memory (LSTM) recurrent neural network is implemented to capture complex temporal dependencies and nonlinear patterns in the data. The model is trained on one year of data and validated on four months of data (one per season), and is compared with a persistence model and SARIMA models. The evaluation shows an overall superiority of the LSTM model, particularly for variables exhibiting high levels of noise and nonlinearity.

While binary classification is one of the most extensively studied problems in machine learning, the regime in which the goal is to learn a classifier with an almost zero false negative rate remains largely unexplored. In this paper, we introduce the Extreme Binary Classification problem, where the objective is to learn a classifier whose false negative rate $\alpha$ is constrained by $\epsilon_{N_1}=o_{N_1\to\infty}(1/N_1)$, with $N_1$ denoting the number of positive examples in the training set. To address this problem, we propose a threshold adaptation method theoretically grounded in guarantees derived from Extreme Value Theory, together with a feature selection procedure based on a permutation test applied to sample maxima. Experimental results on four real-world datasets of varying sizes demonstrate that our approach compares favorably with state-of-the-art methods. In addition, we illustrate its interpretability through an application to a cancer screening dataset.

Probing the hypothesis space after seeing the data remains valid under possibilistic inferential models (IMs), provided the significance level stays fixed. The price is computation, as each plausibility is a supremum of the possibility contour over the hypothesis, and the contour itself is approximated at each queried parameter value. We propose a possibilistic radial transport, which hides the contour value of a parameter in the radius of its source point. When a transport that maximizes within-shell entropy is picked, sampling parameters covering a confidence cut becomes a matter of truncating the radius. We provide a deep learning algorithm that enforces the contour depth condition while maximizing the entropy within each shell. Our amortization makes coverage and power assessments of the learned approximation practical as well as predictive check of new datasets. We also use the sampler to construct a Bel-Pl spectrum for comparing and selecting interpretable hypotheses that satisfy a prescribed Bel-Pl decision criterion. In simulations the learned contours match or improve on ellipsoidal approximations to the cuts, while the coverage and power track the exact reference. Finally, we probe hypotheses about ovarian aging using synthetic AMH records, asking for each woman how many more years her median AMH level will remain above a specified reference value.

In statistical models for categorical responses, including binary, nominal and ordinal responses, a phenomenon called “separation” may occur which can render maximum likelihood estimate nonexistent and/or nonunique. Based on Sablica et al. (2026)'s general characterization of separation phenomena via structure vectors, in this article we discuss computational, numerical and practical aspects of separation phenomena in different software packages, including consequences, remedies and ways of diagnosing them. We further introduce the R package divoRce which provides the most comprehensive software suite to diagnose separation in existence. It allows to check for separation, to distinguish between types of separation, characterize the recession cone and identify variables and observations contributing to separation for almost any categorical response model in the literature, including models with baseline-category link, cumulative link, adjacent-category link and sequential link specifications. Diagnostics are based on methods for computational geometry and linear programming and can be done in floating-point or rational arithmetic with different solvers. The functionality is easily extendible by wrappers and S3 methods to accommodate implementations in third-party packages or of new categorical response models. We discuss our implementation in detail and show our software in action in a wide range of applications. We recommend to incorporate our checks for separation and existence of the MLE as a diagnostic routine in every categorical data analysis both in third-party software implementations and as part of a standard model diagnostic workflow.

Reversible Jump Markov Chain Monte Carlo (RJMCMC) is a principled framework for Bayesian model comparison, but its practical use is often limited due to the difficulty of designing between-model proposals that reach regions of high posterior probability. In addition, its inherently sequential nature limits efficient use of modern GPU architectures. We address both challenges by introducing Anchored Reversible Jump Sequential Monte Carlo. Our approach combines RJMCMC with Sequential Monte Carlo (SMC) to explore the model space in parallel. Importantly, the particle population enables effective between-model proposals without problem-specific knowledge: particles are used to construct kernel density approximations of the target distributions within each model, which are then used to generate trans-dimensional proposals. Numerical experiments show that the proposed method is computationally efficient and competitive with state-of-the-art RJMCMC approaches, demonstrating its potential for scalable Bayesian model comparison.

Sample complexity is a widely used metric in sequential decision-making problems, defined as the number of suboptimal decisions during the interaction between the agent and an environment. We study the sample complexity of stochastic multi-armed bandit problems and introduce the expected sample complexity performance measure, analyzing it in a novel framework called approximately correct in expectation (ACE). We show that ACE guarantees imply almost sure convergence to the optimal expected reward, in contrast to high-probability guarantees found in other frameworks, and also show how to convert ACE guarantees into explicit expected regret bounds. We further show that, in contrast to existing measures, deterministic algorithms cannot obtain favorable ACE bounds, and analyze stochastic algorithms in two settings: when the allowed suboptimality level $\epsilon$ is known to the algorithm and when it is unknown. In the former, we devise an explore-then-$\epsilon$-greedy algorithm, and in the latter, we analyze the expected sample complexity of Thompson sampling. Finally, we establish nearly matching lower bounds for both settings, showing that the algorithms are tight in $\epsilon$ and proving a performance separation between the two regimes.

Hessian spectra at trained models in deep learning exhibit a persistent pattern: eigenvalues organize into distinct clusters, including a large bulk near zero and a few isolated outliers. This paper shows that a natural account of these spectral phenomena emerges when the original setting is understood as a departure from a nearby, otherwise hidden, highly symmetric reference. Modifications, including changes to the architecture, data distribution, or parameter metric, expose a nearby reference configuration whose Hessian exhibits rich invariances-ones not accounted for by weight symmetries. There, symmetry enables a precise description of the spectra, forcing high-dimensional kernels and eigenvalues of large multiplicity. Returning to the original configuration breaks the Hessian symmetry and thereby produces the observed hierarchy of clusters and outliers. The framework is developed in some generality, with a detailed analysis of three-layer ReLU networks and applications to convolutional, graph, and transformer models, as well as to the NTK. The same mechanism is further shown to yield analogous spectral structures in layerwise Hessians and the Gauss-Newton matrix.

Assessing compound weather risks requires forecasts representing dependence between variables. CLARA (Calibrated Advection-Routing Attention) learns joint Gaussian predictive distributions of five surface variables from station observations alone, without numerical weather prediction or reanalysis; the approximately 28,000-parameter model supports CPU training and prediction. Across six multi-year folds on 96 stations, its lead-mean energy score is 4.9% lower than that of a learned comparator with matched temporal inputs (4.7% with a similar parameter count) and 11-65% lower than those of statistical baselines. Holding marginal variances fixed, removing learned correlations worsens joint negative log-likelihood by 1.0-2.8 nats per station. A covariance-scale estimator, proved consistent under stated assumptions, improves short-lead calibration but over-corrects at long leads. Synthetic interventions show an attention-bias coefficient alone does not measure forecast influence. Retrained in ten regions on six continents, CLARA outperforms persistence in all 60 multi-year region-lead comparisons and a similarly sized learned model in 57 of 60.