Papers

Let $\epsilon_{ij}, i,j\geq 1$ be independent Rademacher variables. We prove \begin{equation*} \mathbb{E} \max_{1\leq j\leq p}\left|\frac{1}{n}\sum_{i=1}^n\epsilon_{ij}\right| \geq \min\left\{\frac{255}{256},\frac{1}{\sqrt{2\log 2}}\sqrt{\frac{\log(2p)}{n}}\right\}. \end{equation*} The equality is attained, for instance, by $(n,p)=(2,1)$ and $(n,p)=(2,8).$ We also discuss the optimality of the numerical constants.

In many settings, observed multivariate time series are often nonstationary in nature, i.e., their second order properties vary over time. An additional feature is that their cross-channel dependencies are structured by an underlying network. Together, they give rise to complex interactions between temporal dynamics and network topology. We propose Locally Stationary Wavelet processes on Networks (Net-LSW), a new framework for modelling multiscale, time-varying dependencies that explicitly incorporates the network structure. Unlike traditional multivariate approaches, the Net-LSW process encodes the graph directly in the covariance structure of its random increments. We introduce the concept of local partial correlation graph, mathematically connecting absent edges to zero entries in the time-scale dependent inverse wavelet spectral structure. For inference on the local cross-nodal (partial) dependence, we develop a novel subprocess-based estimation scheme and establish its consistency properties. This new pipeline for network-based nonstationary process modelling, complete with estimation and simulation capabilities that extend outside time-varying vector autoregressive models, is shown to accurately recover evolving dependence structures whilst respecting the underlying graph topology. The analysis of daily stock price volatilities across a global bank network captures multiscale, highly nonstationary dependencies and identifies time-varying systemic shifts during major financial shocks, including Brexit and the COVID-19 pandemic.

Delayed generalization (\ie~grokking) refers to the phenomenon in which a neural network fits its training data early in training but only begins to generalize after a prolonged delay, often through an abrupt transition. Despite extensive empirical study, its underlying mechanism remains poorly understood. In this work, we first theoretically characterize a shell–core topological configuration of the reachable solution space induced by Adam's optimization dynamics with weight-shrinkage regularization, supported by empirical evidence. This optimization-induced topological configuration gives rise to grokking. In model's parameter space, random initialization solutions concentrate on a thin outer spherical shell, enclosing another spherical shell of memorization solutions, which in turn contains a core corresponding to the generalization solutions. Leveraging stopping-time theory, we then analyze the geometry of this topological configuration and the solution transition time at which optimization trajectories escape the memorization manifold and first reach the boundary of the generalization manifold. Our theoretical analysis derives grokking scaling laws for the learning rate, batch size, and $\ell_2$ regularization coefficient, which are further validated through experiments and shown to recover results from prior literature.

We study multiple testing under continuous global–local shrinkage priors, with a focus on the horseshoe prior in high-dimensional sparse settings. While such priors provide adaptive shrinkage and computational scalability, they do not induce exact zeros and hence do not directly yield posterior inclusion probabilities, making principled false discovery control nontrivial. We propose posterior–based decision rules for signal detection that are applicable across a broad class of continuous shrinkage priors and are calibrated to control the false discovery rate (FDR) while retaining high power. In the sparse normal means model, we show that the proposed procedures attain the optimal detection boundary and achieve frequentist asymptotic control of both FDR and false negative rate (FNR). The method is readily implementable via standard posterior sampling, and empirical studies indicate that the realised FDR and FNR closely track their theoretical targets. Applications to high-dimensional regression and Gaussian graphical models further illustrate the scope and practical effectiveness of the approach.

Rapid prototyping of algorithms is a critical step in modern machine learning. Most algorithms exploit linear algebra, creating a need for lightweight numerical routines which – while potentially sub-optimal for the task at hand – can be rapidly implemented. For the numerical solution of ill-conditioned linear systems of equations, the standard solution for prototyping is Tikhonov-regularised inversion using a nugget. However, selection of the size of nugget is often difficult, and the use of data-adaptive procedures precludes automatic differentiation, introducing instabilities into end-to-end training. Further, while data-adaptive procedures perform multiple linear solves to select the size of nugget, only the result of one such solve is returned, which we argue is wasteful. This paper aims to circumvent the above difficulties, presenting autonugget; a Python package for automatic and stable numerical solution of linear systems suitable for rapid prototyping, and fully compatible with automatic differentiation using JAX. autonugget combines multiple linear solves using Richardson extrapolation to determine the solution of the ill-conditioned system, improving in accuracy over approximations based on a single nugget.

In genome wide association studies (GWASs) based on a case-control design, single nucleotide polymorphisms (SNPs) are typically evaluated for an association test and a Hardy-Weinberg equilibrium (HWE) goodness-of-fit test. SNPs are then excluded from analysis based on a HWE cutoff to avoid false positives. In order to avoid cutoffs based on arbitrary threshold values, we propose a conditional genotype–based test that conditions the Pearson $\chi^2$-statistic in the 3x2 contingency table on the $\chi^2$-statistic for HWE in the control group, and develop the relevant asymptotic distribution theory. We show by simulations that our test in most scenarios is more powerful than two competing retrospective procedures. Another important advantage of the proposed method is a better ranking of SNPs in GWASs as HWE is accounted for in computing p-values of SNP association. We demonstrate this effect on a data set in an alopecia study. In conclusion, our test makes separate HWE testing superfluous by providing a unified framework and strictly improves on the standard procedure in terms of power and interpretability, thereby making replication more cost effective and improving subsequent fine mapping.\par

The Sliced Wasserstein (SW) distance has emerged as a computationally attractive alternative to the Wasserstein distance by leveraging one-dimensional optimal transport along random projections. Standard estimators of the SW distance rely on Monte Carlo averages of one-dimensional Wasserstein distances computed via quantile functions, which require sorting projected samples and access to full datasets. In this work, we introduce a new class of estimators for the Sliced Wasserstein distance based on cumulative distribution functions (CDFs) of projected measures, that avoid sorting and scale via massive dataset parallelism. This class includes several estimators, some of them being indexed by hyperparameters controlling their variance or smoothness. We show that they are especially well suited to scenarios in which CDFs are more tractable than quantile functions, such as mixtures of Gaussians, and moreover that they are also naturally compatible with federated learning, since CDFs of projected data can be computed and aggregated locally without requiring the exchange of raw samples.

Differential equation models are widely used to describe, interpret, and predict dynamical phenomena across science and engineering. In practice, however, the governing dynamics are rarely fully known and must be inferred from observational data. Traditionally, inverse problems in differential equation modelling have focused on estimating unknown parameter values. In this setting, structural identifiability determines whether parameter values can, in principle, be uniquely recovered from ideal observations and is, therefore, a prerequisite for meaningful inference. More recently, the integration of machine learning with mechanistic modelling has enabled the discovery of unknown equations, functions, and constitutive relationships, substantially expanding the space of admissible models. This raises a fundamental question: under what conditions can unknown functional components be uniquely recovered from data? In this paper, we generalise the classical notion of structural parameter identifiability to functional identifiability. We first identify broad classes of models for which unique functional recovery is impossible. We then show how functional identifiability can be assessed for differential equation models using differential algebra-based techniques which are well-established as a means of assessing structural identifiability for ordinary differential equation-based models. Our framework reveals new phenomena that arise in the transition from parametric to functional inference and have no analogue in the classical setting. Finally, we characterise functional identifiability in several common model classes. Taken together, our results demonstrate that functional identifiability provides a theoretical foundation for modern inverse problems in differential equation modelling, particularly those that use machine learning representations of unknown system components.

Speculative decoding accelerates language model inference by using a fast drafter to propose candidate tokens that are then verified by a larger target model. Existing theory largely studies the stochastic, distribution-preserving setting, where the goal is to exactly sample from the target distribution. In contrast, many practical systems use greedy decoding, relaxed acceptance rules, or tree-based candidate sets, where success is governed by local ranking and threshold events rather than exact distributional equality. We develop a theory for these regimes. We identify that many common acceptance criteria have rejection regions that can be characterized as lower level sets of the target distribution. For these, we characterize the exact KL divergence required for rejection yielding exact certificates and sharp margin-based bounds for strict greedy decoding, additive and multiplicative relaxed acceptance, top-(m) relaxed criteria, and entropy-thresholded acceptance. We then extend the framework to greedy tree decoding, deriving exact and margin-only certificates for when the target greedy token remains covered by the drafter's top-(m) candidates. Finally, we evaluate the resulting certificates on Qwen3 models, showing that relaxed and tree-based criteria substantially enlarge the region of certified acceptance, especially on decoding steps with low target model distribution margin. These results complement existing distribution-preserving analyses of speculative decoding by characterizing the deterministic local acceptance events common in practical inference systems.

Valid possibilistic inferential models provide exact finite-sample calibration, but validity alone does not determine which valid procedure results in the most informative inferential summary. This paper proposes Choquet risk as a decision-theoretic criterion for comparing valid possibility measures in finite samples. Given a non-negative penalty functional, Choquet loss is defined as the Choquet integral of that penalty with respect to the data-dependent possibility measure, and Choquet risk as its sampling expectation. A key reduction expresses this risk through the nested $\alpha$-cuts of the contour, linking procedure-level efficiency to the expected performance of calibrated confidence sets. For concentration penalties, the criterion reduces to integrated expected set size, equivalently expected contour volume, so levelwise optimal confidence sets induce Choquet-risk optimal valid contours. The framework is developed along two classical routes to optimality. First, a possibilistic notion of unbiasedness is introduced and shown, under validity, to coincide with unbiasedness of the induced confidence sets and tests, allowing UMPU and most-accurate-unbiased results to be transferred to valid contours. Second, an equivariant minimax theory is developed, including a Gaussian-location result in which the Gaussian possibility contour is Choquet-risk minimax for radial distance-to-truth losses. The construction also extends confidence risk from additive confidence distributions to non-additive calibrated inferential-model output, with Choquet loss acting as a least-favourable confidence loss. Finally, the paper clarifies the penalty-dependence of efficiency comparisons and motivates invariant size criteria and divergence-based intrinsic losses connected locally to Fisher–Rao geometry.

These notes recapitulate the high level mathematical principles behind different techniques for generative modeling. I show the connections between optimal transport and standard techniques such as Schrödinger bridge and flow matching.

We determine the exact region of Chatterjee's rank correlation $\xi$ and Blomqvist's $\beta$, that is, the set of all pairs $(\xi(C),\beta(C))$ where $C$ ranges over all bivariate copulas. The region equals $\{(x,y)\in[0,1]\times[-1,1]:|y|^3\le2x\}$, so every copula satisfies the sharp inequality $|\beta(C)|^3\le2\xi(C)$. Its left boundary $\xi=|\beta|^3/2$ is attained uniquely by an explicit family of copulas that perturb independence by a signed tent function centered at the median, and its right boundary $\xi=1$ by deterministic interval-exchange copulas. The proof reduces the problem to a one-dimensional minimization of an $L^2$-norm over $1$-Lipschitz functions. We further show that the region does not shrink when restricted to radially symmetric copulas and that positively quadrant dependent copulas attain exactly its upper half. For copulas under which the second variable is stochastically increasing in the first, the exact region is $\{(x,y)\in[0,1]^2:y^3\le2x,\ 3(1-y)^2\le4(1-x)\}$, and its right boundary is attained uniquely by explicit singular two-branch copulas. In particular, under stochastic increasingness $\beta=0$ is compatible with every value $\xi\in[0,1/4]$, whereas $\xi>1/4$ forces $\beta>0$.

We investigate adaptive increasingly rare Markov chain Monte Carlo algorithms and the associated time-average estimator for approximating expectations. Under a simultaneous Wasserstein contraction assumption on the underlying family of Markov kernels we derive explicit bounds for the mean squared error. We illustrate the applicability of our estimate through adaptive stereographic algorithms and Metropolis-Hastings schemes that employ normalizing flows for adaptation. We also consider a generic adaptive algorithm for doubly intractable problems and provide a corresponding cost analysis to achieve a desired precision.

We consider nonparametric estimation of the regression function in a model where individuals share a common noise component and repeated measurements are available for each individual. We propose a projection estimator which minimizes a least-squares contrast that accounts for the covariance structure resulting from the common noise. We analyze its risk measured either as the expectation of the empirical norm or as the expectation of the theoretical norm associated with the contrast. We discuss how the number of repeated measurements affects the estimation rates in the common noise model, and precisely characterize the dependence on the number of repetitions. In addition, we propose a data-driven projection estimator and establish risk bounds in terms of the expected empirical norm. The results are illustrated with some simulation experiments.

We are interested in conducting political polls sequentially, so that one can stop acquiring data as soon as possible while safely yielding statistically significant results. Building off e-values, which have recently become a useful tool to create sequential testing methods, we develop a theory of posterior optimal e-values. We use voting as a convenient example on which to illustrate our method. First, we design statistical tests for Condorcet and Borda voting system, and also for Schulze voting system which we are the first to tackle statistically. Then, we study the construction of optimal sequential e-values in the deceptively simple setting of multivariate Bernoulli data, with general composite null and alternative hypothesis sets $\mathcal{H}_0$ and $\mathcal{H}_1$. We give a way to compute these e-values using an efficient Frank-Wolfe algorithm, giving a pretty general way to compute Reverse Information Projections, even when $\mathcal{H}_0$ corresponds to a non-convex parameter set. Finally, we illustrate the efficiency, both in terms of power and sample size of our method. We compare with state of the art in both simulated and real data experiments, with application to French 2022 presidential election data.

We propose an extension of the ordered stereotype model (OSM) for ordinal time series data, referred to as the Autoregressive OSM (AR-OSM). The model captures serial dependence by incorporating lagged values of the response as covariates in the systematic component. In contrast to existing regression models for ordinal time series, the AR-OSM does not assume equidistant categories, but instead allows the data to determine their relative spacing. This property makes the model particularly suitable for applications where the equidistance assumption is unrealistic. Such a case is illustrated through the analysis of infant sleep state data. Additionally, a comprehensive simulation study is conducted to assess the performance of the model under varying sample sizes and to investigate how parameter values influence the induced serial dependence structure.

This research reviews several crucial aspects of the universal model for the Lorenz curve proposed by Sitthiyot and Holasut (2023) (hereafter, SH (2023)). A first issue concerns the mathematical definition of the proposed curves. The four functional forms introduced by SH (2023) do not satisfy the necessary and sufficient conditions for a valid Lorenz curve. We propose corrected versions of the previous curves and derive analytical expressions for some measures of inequality.

Decision theory provides a formal framework for how agents should make choices under uncertainty, drawing on ideas from philosophy, probability, and causality. Despite significant progress, the field still lacks a unified modeling language, and key concepts - such as the distinction between subjective and objective elements, or what it means for a decision theory to perform well - are often left implicit. This can make it difficult to evaluate and compare competing theories, particularly in controversial cases. In this paper, we address these issues by introducing a formal framework for decision theory based on nonparametric structural equation models (NPSEMs), a well-established tool in causal inference. NPSEMs provide a unified foundation for representing agents, counterfactuals, and causal relationships, allowing for unambiguous definitions of EDT and CDT. Building on this foundation, we propose a novel decision theory - personal decision theory - which instructs agents to maximize a subjective model of their own counterfactual utility. We introduce a formal performance metric based on hypothetical interventions that enforce a given decision theory across a population - such as might be achieved through education or policy – and show that, under certain assumptions, personal decision theory is optimal with respect to this metric. Throughout, we use the smoking lesion problem as a running example and conclude with a formal analysis of Newcomb's problem. Our aim is to provide decision theory with a clearer modeling language and firmer evaluative ground, thereby enabling more rigorous comparisons and facilitating conceptual progress in the field.

We exhibit a simple deterministic one-dimensional convex composite optimization problem for which AdaGrad scheme does not achieve the classical convergence rate $\mathcal{O}(n^{-(1+\nu)/2})$ associated with Hölder-smooth objectives. The example highlights a basic mismatch between classical AdaGrad accumulation and composite optimality. A main insight is that the gradient of the smooth term may not vanish at the optimum, causing AdaGrad to keep reducing its stepsize excessively and converge more slowly. We also discuss why alternative accumulation mechanisms based on gradient mappings or on successive gradient differences, avoid this pathology.

Predictive models are increasingly embedded in operational decision-making, yet standard explanation methods typically explain forecasts rather than the decisions those forecasts induce. This distinction is important in predict-then-optimize systems: large forecast changes may leave the optimizer's action unchanged, while small changes can alter the selected decision and its realized value. We propose Decision Value Attribution (DVA), a Shapley-based framework for attributing the value of a fixed prediction–optimization pipeline. The framework defines cooperative games whose payoff is the downstream decision value, allowing the players to be information sources, optimization or design parameters, or both. We present three variants: InfoDVA attributes value to features, DesignDVA attributes value to operational configurations, and Decision-Value Interactions (DVI) quantifies how information and design jointly create value. We further distinguish post-DVA, which evaluates decisions using realized outcomes, from pre-DVA, which evaluates decisions under the model's full prediction. This separation turns attribution into a decision-level diagnostic of whether the model's operational beliefs align with realized performance. The resulting attributions are expressed in the units of the operational objective and decompose the gain or loss relative to a baseline. Case studies in electricity storage arbitrage and emergency medical service coverage show that predictive explanations can be poor proxies for operational value, that DVA can guide targeted information-control interventions, and that optimization configurations determine when predictive information is decision-relevant.

We introduce a method for visual and auditory feedback when exploring the fit of a model to data. Starting with a best-fit curve fit to data, the user can drag the curve to a new position and the computer will emit a squeal, becoming louder and more unpleasant as the discrepancy between curve and data increases. We demonstrate with four examples: a two-parameter curve fit to golf putting data, a four-parameter curve fit to dilution assays, a fit to cosmological data sensitive to the parameters of the Big Bang model, and a nonparametric Gaussian process fit to temperature readings.

In 2020, the U.S. Census Bureau adopted differential privacy for the Decennial Census by injecting integer-valued Gaussian noise into published census tabulations. Exactly evaluating the privacy guarantees of these data releases would enable the Bureau to determine the absolute minimum noise required to satisfy a given privacy budget, preventing the injection of unnecessary excess noise and thereby substantially enhancing the statistical utility of the data for downstream applications such as federal funding allocation and political redistricting. In this paper, we introduce a computationally efficient and mathematically rigorous quadrature method to evaluate the exact privacy profile of practical, large-scale census releases under the composition of heterogeneous discrete Gaussian mechanisms. Mathematically, this problem reduces to evaluating the tail probabilities of high-dimensional convolutions of integer-valued random variables sampled from heterogeneous discrete Gaussian distributions under exceptionally stringent numerical error tolerances (e.g., $10^{-35}$). By recasting the exact privacy accounting as a numerical integration problem via the discrete Fourier transform, we explicitly exploit the exponential convergence of the trapezoidal rule for complex analytic, periodic characteristic functions. Furthermore, to overcome the computational bottleneck of evaluating highly oscillatory integrands in high dimensions, we develop a sieve algorithm that identifies and prunes negligible quadrature nodes, accelerating the computation by three orders of magnitude. Taken together, these numerical innovations enable the first exact, assumption-free privacy accounting for the 2020 Census Demographic and Housing Characteristics File, achieving a 1,824-fold speedup over prior methods while maintaining census-mandated error tolerances.

This study proposes coarse-to-fine downscaling (CF-DS), a scalable spatial downscaling method extending coarse-to-fine spatial modeling. Unlike conventional spatial-statistical downscaling methods such as area-to-point kriging, CF-DS does not require covariance matrix inversion or likelihood evaluation. Instead, it represents latent spatial processes through the synthesis of multi-scale local models, substantially reducing computational cost while approximately satisfying the aggregation constraint. Monte Carlo experiments show that CF-DS achieves predictive accuracy comparable to area-to-point kriging with dramatically shorter computation times, particularly for large datasets. An application to downscaling electricity consumption in the Tokyo metropolitan area further demonstrates its practical usefulness. The results suggest that CF-DS provides an efficient alternative for large-scale spatial downscaling problems. CF-DS is implemented in an R package spCF (https://cran.r-project.org/web/packages/spCF/index.html).

Outlier detection (OD) aims to identify anomalous instances by learning the underlying structure of normal data (inliers), and is particularly challenging in fully unsupervised settings where no information about anomalies is available during training. Recent advances have leveraged the inlier-memorization (IM) effect, a phenomenon in which deep models memorize inlier patterns earlier than those of outliers, as a powerful signal for distinguishing outliers. However, despite its empirical success, the theoretical understanding of the IM effect remains limited. In this work, we present a theoretical study of the IM effect. Focusing on a simple autoencoder, we show that, under mild assumptions, the model can successfully memorize inliers while failing to memorize outliers during certain stages of early training. In particular, we characterize not only the emergence of the IM effect, but also its strength and persistence, and analyze how these properties depend on the data distribution and parameter initialization. In addition, building on these insights, we derive simple yet practical guidelines for enhancing the IM effect, including data preprocessing and parameter initialization schemes, achieving state-of-the-art performance on the ADBench datasets. Our findings provide a theoretical foundation for the IM effect and offer actionable directions for improving IM-based outlier detection methods.

Reliable generative AI models critically rely on expert human annotations to evaluate output quality, yet these "gold" labels are expensive to collect and limited in quantity. Organizations thus often turn to collecting vast but noisy "silver" labels from crowdsourced workers or vendor annotators as proxies for gold labels. Because gold remains the evaluation target, naively aggregating noisy silver labels may introduce bias, and estimators built on sparsely observed gold labels may have high variance to resolve the model performance gaps that guide practical decisions. Model evaluation has become an ongoing operational practice rather than a one-time exercise, with evaluation rounds repeating across model versions, releases, and content domains. A natural question is whether the previous historical evaluation data can be used to improve each new round of evaluation. We introduce HERO (History Enhanced RObust model evaluation), a novel framework that uses historical data to suppress bias (improve reliability) and reduce variance (improve sensitivity) in model performance evaluation. HERO calibrates silver labelers' performance learned from historical gold annotations, and stabilizes the resulting estimator by anchoring it to covariate information measured with high precision in the historical data. HERO can be broadly applied across multiple common evaluation tasks, and remains valid when only a subset of historical labelers appears in the current round. We establish conditions under which the bias and variance reductions hold, showcase HERO's performance in simulation studies, and demonstrate its effectiveness on real-world model evaluation benchmarking datasets.

Large language models (LLMs) are interactive stochastic systems whose most consequential behaviors are still only partially understood. This discussion argues that statistics curricula should treat LLMs not only as tools, but as objects of inquiry: students can probe variability, bias, and prompt sensitivity by designing small experiments and analyzing distributions of outputs. Building on the Veridical Data Science framework and Predictability-Computability-Stability (PCS) principles, this discussion outlines how to organize critical LLM engagement across educational levels and propose four curricular examples, from introductory “ask it twice” activities to graduate PCS stability audits of LLM-based analysis workflows.

Classical hypothesis testing frameworks break down in contemporary settings in which null hypotheses are increasingly abstract, the same data are used to both generate and test hypotheses, and minimal assumptions about the underlying data are made. In this work, we propose a new framework for conducting valid hypothesis tests in broad contexts. We propose to add and subtract external noise generated from a symmetric shift-family to our data, $X$, to partition it into two pieces, $X^{(1)}$ and $X^{(2)}$. We provide a generic strategy for orthogonalizing $X^{(2)}$ against $X^{(1)}$ under the null hypothesis $H_0$, then show that testing whether the orthogonalization was successful provides a valid test of $H_0$ under mild assumptions. Remarkably, this framework extends naturally to the post-selection inference setting: we simply select a hypothesis on $X^{(1)}$, then perform orthogonalization under the selected null. As our approach neither requires pre-specification of the selection mechanism, nor is restricted to a small class of data-generating distributions, it dramatically expands the settings for which valid post-selection inference can be conducted. We showcase the flexibility of our proposal in several case studies involving challenging pre-specified null hypotheses and post-selection inference scenarios.

Bogomolny, Bohigas and Schmit (BBS) found that the spectrum of the pairwise distance matrix on N points sampled from a smooth d-dimensional manifold encodes a signature of the underlying geometry. We develop I-BBS (Inference-BBS), a coordinate-free method that identifies a low-dimensional latent sub-manifold embedded in a high-dimensional ambient distance matrix alone, without accessing an ambient high-dimensional vector space. It therefore applies even when that space is only partly observable or undefined. We model the ambient embedding by two classes of generative noise, model-based and model-free. The noise mixes the latent signal with off-manifold components, so the eigenvalues reorganise collectively and the latent geometry cannot be read off eigenvalue by eigenvalue. We recover it instead from two integer-stable signatures that survive the noise: the multiplicity of the top non-Perron multiplet, which fixes $d$, and a parameter-free law for how the multiplet positions shrink as the noise grows. On synthetic spheres $S^1$, $S^2$ and $S^3$ these integer signatures are far more stable under noise than the continuous spectral slope, and a blind test recovers both the manifold and the noise model from a single distance matrix. Applications to neural-network representations and to the dynamic training regime are developed in two companion papers.

This paper examines how metric adjustments to Multidimensional Scaling (MDS) can enhance its effectiveness as a visual tool for pattern recognition. The distance under consideration, referred to as Max-D-SW, is an adjustment of the Max-Sliced Wasserstein distance. In contrast to the original formulation, which optimizes over single unit directions, Max-D-SW aggregates contributions over orthonormal bases. This modification provides a clear numerical advantage in MDS outcomes, particularly when applied to heavy-tailed distributions. We also establish sample-complexity bounds showing that Max-D-SW remains statistically tractable, with rates comparable to those of its max-sliced counterpart. Moreover, we show that a better sample complexity for a metric does not necessarily translate into better performance when the metric is used as an input for MDS.

Modern optimizers combine gradients from the current mini-batch with historical optimization state, such as momentum or adaptive moments. While effective, this standard practice can produce parameter updates that actively increase the loss of individual samples. We term this phenomenon per-sample interference and propose redefining the parameter update as an optimization problem that explicitly minimizes it. Because the exact formulation of the problem is computationally prohibitive, we introduce a highly efficient surrogate. By reducing the problem's dimensionality to the batch size and restricting the optimization to the last linear layer, we overcome memory and speed bottlenecks. This strategy hinges on our unexpected finding that this layer alone can reliably capture core second-order statistics of the full network. The resulting surrogate problem integrates readily into standard optimizers like SGD and AdamW, and can be solved using a small number of GPU-friendly iterations. Crucially, the method exhibits favorable scaling properties, as the relative computational overhead shrinks as the model size or input grows. Experiments on image classification benchmarks confirm reduced per-sample interference and improved generalization.