Papers

Covariance matrices arise naturally in different scientific fields, including finance, genomics, and neuroscience, where they encode dependence structures and reveal essential features of complex multivariate systems. In this work, we introduce a comprehensive Bayesian framework for analyzing heterogeneous covariance data through both classical mixture models and a novel mixture-of-experts Wishart (MoE-Wishart) model. The proposed MoE-Wishart model extends standard Wishart mixtures by allowing mixture weights to depend on predictors through a multinomial logistic gating network. This formulation enables the model to capture complex, nonlinear heterogeneity in covariance structures and to adapt subpopulation membership probabilities to covariate-dependent patterns. To perform inference, we develop an efficient Gibbs-within-Metropolis-Hastings sampling algorithm tailored to the geometry of the Wishart likelihood and the gating network. We additionally derive an Expectation-Maximization algorithm for maximum likelihood estimation in the mixture-of-experts setting. Extensive simulation studies demonstrate that the proposed Bayesian and maximum likelihood estimators achieve accurate subpopulation recovery and estimation under a range of heterogeneous covariance scenarios. Finally, we present an innovative application of our methodology to a challenging dataset: cancer drug sensitivity profiles, illustrating the ability of the MoE-Wishart model to leverage covariance across drug dosages and replicate measurements. Our methods are implemented in the R package moewishart available at https://github.com/zhizuio/moewishart .

Statistical hypothesis tests typically use prespecified sample sizes, yet data often arrive sequentially. Interim analyses invalidate classical error guarantees, while existing sequential methods require rigid testing preschedules or incur substantial losses in statistical power. We introduce a simple procedure that transforms any fixed-sample hypothesis test into an anytime-valid test while ensuring Type-I error control and near-optimal power with substantial sample savings when the null hypothesis is false. At each step, the procedure predicts the probability that a classical test would reject the null hypothesis at its fixed-sample size, treating future observations as missing data under the null hypothesis. Thresholding this probability yields an anytime-valid stopping rule. In areas such as clinical trials, stopping early and safely can ensure that subjects receive the best treatments and accelerate the development of effective therapies.

We formulate Ensemble-Conditional Gaussian Processes (Ens-CGP), a finite-dimensional synthesis that centers ensemble-based inference on the conditional Gaussian law. Conditional Gaussian processes (CGP) arise directly from Gaussian processes under conditioning and, in linear-Gaussian settings, define the full posterior distribution for a Gaussian prior and linear observations. Classical Kalman filtering is a recursive algorithm that computes this same conditional law under dynamical assumptions; the conditional Gaussian law itself is therefore the underlying representational object, while the filter is one computational realization. In this sense, CGP provides the probabilistic foundation for Kalman-type methods as well as equivalent formulations as a strictly convex quadratic program (MAP estimation), RKHS-regularized regression, and classical regularization. Ens-CGP is the ensemble instantiation of this object, obtained by treating empirical ensemble moments as a (possibly low-rank) Gaussian prior and performing exact conditioning. By separating representation (GP -> CGP -> Ens-CGP) from computation (Kalman filters, EnKF variants, and iterative ensemble schemes), the framework links an earlier-established representational foundation for inference to ensemble-derived priors and clarifies the relationships among probabilistic, variational, and ensemble perspectives.

Modern online experimentation faces two bottlenecks: scarce traffic forces tough choices on which variants to test, and post-hoc insight extraction is manual, inconsistent, and often content-agnostic. Meanwhile, organizations underuse historical A/B results and rich content embeddings that could guide prioritization and creative iteration. We present a unified framework to (i) prioritize which variants to test, (ii) explain why winners win, and (iii) surface targeted opportunities for new, higher-potential variants. Leveraging treatment embeddings and historical outcomes, we train a CTR ranking model with fixed effects for contextual shifts that scores candidates while balancing value and content diversity. For better interpretability and understanding, we project treatments onto curated semantic marketing attributes and re-express the ranker in this space via a sign-consistent, sparse constrained Lasso, yielding per-attribute coefficients and signed contributions for visual explanations, top-k drivers, and natural-language insights. We then compute an opportunity index combining attribute importance (from the ranker) with under-expression in the current experiment to flag missing, high-impact attributes. Finally, LLMs translate ranked opportunities into concrete creative suggestions and estimate both learning and conversion potential, enabling faster, more informative, and more efficient test cycles. These components have been built into a real Adobe product, called Experimentation Accelerator, to provide AI-based insights and opportunities to scale experimentation for customers. We provide an evaluation of the performance of the proposed framework on some real-world experiments by Adobe business customers that validate the high quality of the generation pipeline.

We propose a sequential test for detecting arbitrary distribution shifts that allows conformal test martingales (CTMs) to work under a fixed, reference-conditional setting. Existing CTM detectors construct test martingales by continually growing a reference set with each incoming sample, using it to assess how atypical the new sample is relative to past observations. While this design yields anytime-valid type-I error control, it suffers from test-time contamination: after a change, post-shift observations enter the reference set and dilute the evidence for distribution shift, increasing detection delay and reducing power. In contrast, our method avoids contamination by design by comparing each new sample to a fixed null reference dataset. Our main technical contribution is a robust martingale construction that remains valid conditional on the null reference data, achieved by explicitly accounting for the estimation error in the reference distribution induced by the finite reference set. This yields anytime-valid type-I error control together with guarantees of asymptotic power one and bounded expected detection delay. Empirically, our method detects shifts faster than standard CTMs, providing a powerful and reliable distribution-shift detector.

Large language models spend most of their inference cost on attention over long contexts, yet empirical behavior suggests that only a small subset of tokens meaningfully contributes to each query. We formalize this phenomenon by modeling attention as a projection onto the convex hull of key vectors and analyzing its entropic (softmax-like) relaxation. Our main theoretical contribution is a face-stability theorem showing that, under a strict complementarity margin (a support gap (\Delta) certified by KKT multipliers), entropic attention concentrates on a constant-size active face: the total mass assigned to inactive tokens decays exponentially as (\exp(-\Omega(\Delta/\varepsilon))), while the error on the active face scales linearly in the temperature/regularization parameter (\varepsilon). This yields a practical criterion for when sparse long-context decoding is safe and provides a principled knob to trade accuracy for compute. Building on these guarantees, we introduce Vashista Sparse Attention, a drop-in mechanism that maintains a small candidate set per query through a paging-style context selection strategy compatible with modern inference stacks. Across long-context evaluations, we observe stable constant-size effective support, strong wall-clock speedups, and minimal quality degradation in the regimes predicted by the support-gap diagnostics. Finally, we discuss deployment implications for privacy-sensitive and air-gapped settings, where interchangeable attention modules enable predictable latency and cost without external retrieval dependencies.

Missing values in datasets are common in applied statistics. For regression problems, theoretical work thus far has largely considered the issue of missing covariates as distinct from missing responses. However, in practice, many datasets have both forms of missingness. Motivated by this gap, we study linear regression with a labelled dataset containing missing covariates, potentially alongside an unlabelled dataset. We consider both structured (blockwise-missing) and unstructured missingness patterns, along with sparse and non-sparse regression parameters. For the non-sparse case, we provide an estimator based on imputing the missing data combined with a reweighting step. For the high-dimensional sparse case, we use a modified version of the Dantzig selector. We provide non-asymptotic upper bounds on the risk of both procedures. These are matched by several new minimax lower bounds, demonstrating the rate optimality of our estimators. Notably, even when the linear model is well-specified, our results characterise substantial differences in the minimax rates when unlabelled data is present relative to the fully supervised setting. Particular consequences of our sparse and non-sparse results include the first matching upper and lower bounds on the minimax rate for the supervised setting when either unstructured or structured missingness is present. Our theory is coupled with extensive simulations and a semi-synthetic application to the California housing dataset.

Forecasting presents a complex estimation challenge, as it involves balancing multiple, often conflicting, priorities and objectives. Conventional forecast optimization methods typically emphasize a single metric–such as minimizing the mean squared error (MSE)–which may neglect other crucial aspects of predictive performance. To address this limitation, the recently developed Smooth Sign Accuracy (SSA) framework extends the traditional MSE approach by simultaneously accounting for sign accuracy, MSE, and the frequency of sign changes in the predictor. This addresses a fundamental trade-off–the so-called accuracy-smoothness (AS) dilemma–in prediction. We extend this approach to the multivariate M-SSA, leveraging the original criterion to incorporate cross-sectional information across multiple time series. As a result, the M-SSA criterion enables the integration of various design objectives related to AS forecasting performance, effectively generalizing conventional MSE-based metrics. To demonstrate its practical applicability and versatility, we explore the application of the M-SSA in three primary domains: forecasting, real-time signal extraction (nowcasting), and smoothing. These case studies illustrate the framework's capacity to adapt to different contexts while effectively managing inherent trade-offs in predictive modelling.

Particle smoothing enables state estimation in nonlinear and non-Gaussian state-space models, but its practical use is often limited by high computational cost. Backward smoothing methods such as the Forward Filter Backward Smoother (FFBS) and its marginal form (FFBSm) can achieve high accuracy, yet typically require quadratic computational complexity in the number of particles. This paper examines the accuracy–computational cost trade-offs of particle smoothing methods through a trend-estimation example. Fixed-lag smoothing, FFBS, and FFBSm are compared under Gaussian and heavy-tailed (Cauchy-type) system noise, with particular attention to O(m) approximations of FFBSm based on subsampling and local neighborhood restrictions. The results show that FFBS and FFBSm outperform fixed-lag smoothing at a fixed particle number, while fixed-lag smoothing often achieves higher accuracy under equal computational time. Moreover, efficient FFBSm approximations are effective for Gaussian transitions but become less advantageous for heavy-tailed dynamics.

We study parametric change-point detection under local differential privacy, focusing on the offline localization of a single change point in a sequence of independent observations.~In the non-private setting, we derive improved finite-sample accuracy guarantees for a generalized likelihood-ratio estimator using martingale methods.~We establish finite-sample minimax lower bounds in both the non-private and locally private settings, showing that localization accuracy is fundamentally limited by the Chernoff information between the pre- and post-change distributions.~Guided by this, we design locally private CPD algorithms based on aggregated randomized response and binary mechanisms, which maximize the Chernoff information between the induced output distributions in low- and high-privacy regimes, respectively, and derive their finite-sample guarantees. Together, our upper and lower bounds quantify the statistical cost of privacy in change-point localization.~As part of our analysis, we also show that strong data processing inequality coefficients of Rényi divergence and Jeffreys–Rényi divergence of order infinity can be computed by restricting to pairs of input distributions with a common support of at most two symbols.~We corroborate our results through empirical evaluation on synthetic data.

Stochastic differential equations (SDEs) offer powerful and accessible mathematical models for capturing both deterministic and probabilistic aspects of dynamic behavior across a wide range of physical, financial, and social systems. However, analytical solutions for many SDEs are often unavailable, necessitating the use of numerical approximation methods. The rate of convergence of such numerical methods is of great importance, as it directly influences both computational efficiency and accuracy. This paper presents a proposed theorem, along with its proof, that facilitates the numerical evaluation of the strong (and weak) order of convergence of a numerical scheme for an SDE when the analytical solution is unavailable. Additionally, we address the challenge of numerically computing the multiple stochastic integrals required by the Milstein method to achieve improved convergence rates for multidimensional SDEs. In this context, two newly proposed numerical techniques for computing these multiple stochastic integrals are introduced and compared with existing approaches in terms of efficiency and effectiveness. The methodologies are further illustrated through simulation studies and applications to widely used financial models.

Motivated by applications in tissue-wide association studies (TWAS), we develop a flexible and theoretically grounded empirical Bayes approach for integrating %vector-valued outcomes data obtained from different sources. We propose a linear shrinkage estimator that effectively shrinks singular values of a data matrix. This problem is closely connected to estimating covariance matrices under a specific loss, for which we develop asymptotically optimal estimators. The basic linear shrinkage estimator is then extended to a local linear shrinkage estimator, offering greater flexibility. Crucially, the proposed method works under sparse/dense or low-rank/non low-rank parameter settings unlike well-known sparse or reduced rank estimators in the literature. Furthermore, the empirical Bayes approach offers greater scalability in computation compared to intensive full Bayes procedures. The method is evaluated through an extensive set of numerical experiments, and applied to a real TWAS data obtained from the Genotype-Tissue Expression (GTEx) project.

The win ratio (WR) is a widely used metric to compare treatments in randomized clinical trials with hierarchically ordered endpoints. Counting-based approaches, such as Pocock's algorithm, are the standard for WR estimation. However, this algorithm treats participants with censored or missing data inadequately, which may lead to biased and inefficient estimates, particularly in the presence of heterogeneous censoring or missing data between treatment groups. Although recent extensions have addressed some of these limitations for hierarchical time-to-event endpoints, no existing methods – aside from the computationally intensive multiple imputation approach – can accommodate settings that include non-survival endpoints that are subject to missing data. In this paper, we propose a simple nonparametric maximum likelihood estimator (NPMLE) of WR for two hierarchical endpoints that are subject to censoring and missing data. Our method uses all observed data, avoids strong parametric assumptions, and comes with a closed-form asymptotic variance estimator. We demonstrate its performance using simulation studies and two data examples, based on the HEART-FID and ISCHEMIA trials. The proposed method provides a consistent estimator, improves estimation efficiency, and is robust under non-informative censoring and missing at random (MAR) assumptions, offering a flexible alternative to existing WR estimation methods. A user-friendly R package, WinRS, is available to facilitate implementation.

Standard direct forecasting models typically rely on point-wise objectives such as Mean Squared Error, which fail to capture the complex spatio-temporal dependencies inherent in graph-structured signals. While recent frequency-domain approaches such as FreDF mitigate temporal autocorrelation, they often overlook spatial and cross spatio-temporal interactions. To address this limitation, we propose FreST Loss, a frequency-enhanced spatio-temporal training objective that extends supervision to the joint spatio-temporal spectrum. By leveraging the Joint Fourier Transform (JFT), FreST Loss aligns model predictions with ground truth in a unified spectral domain, effectively decorrelating complex dependencies across both space and time. Theoretical analysis shows that this formulation reduces estimation bias associated with time-domain training objectives. Extensive experiments on six real-world datasets demonstrate that FreST Loss is model-agnostic and consistently improves state-of-the-art baselines by better capturing holistic spatio-temporal dynamics.

There is an intense and partly recent literature focussing on the problem of selecting the bandwidth parameter for kernel density estimators. Available methods are largely `very nonparametric', in the sense of not requiring any knowledge about the underlying density, or `very parametric', like the normality-based reference rule. This report aims at widening the scope towards the inclusion of many semiparametric bandwidth selectors, via Hermite type expansions aroundthe normal distribution. The resulting bandwidths may be seen as carrying out suitable corrections on the normal reference rule, requiring a low number of extra coefficients to be estimated from data. The present report introduces and discusses some basic ideas and develops the necessary initial theory, but modestly chooses to stop short of giving precise recommendations for specific procedures among the many possible constructions. This will require some further analysis, numerical work, and some simulation-based exploration.

Networks of modern industrial systems are increasingly monitored by distributed sensors, where each system comprises multiple subsystems that generate high-dimensional time-series. These subsystems are often interdependent, making it important to understand how temporal patterns observed at one subsystem relate to those observed across others. This problem is particularly challenging in decentralized settings where raw measurements cannot be shared, and client observations are heterogeneous. Moreover, in practice, each subsystem (client) operates a fixed, proprietary model that cannot be modified or retrained, limiting the applicability of existing approaches. Separately, nonlinear dynamics make cross-client temporal interdependencies difficult to characterize, as they are embedded implicitly within nonlinear state-transition functions. This paper presents a vertical federated learning framework for learning one-step-ahead cross-client temporal interdependencies under these constraints. Each client maps high-dimensional local observations to low-dimensional latent states using a nonlinear state-space model, while a Graph Attention Network (GAT)-based central server learns a state-transition model over the communicated states. Given a candidate graph skeleton, the server learns state-dependent attention coefficients that capture cross-client interactions. We analytically relate the GAT attention coefficients to the Jacobian of the learned state-transition model, yielding a first-order characterization of cross-client interdependencies. We further establish bounded closeness of the learned dynamics and Jacobians to a corresponding centralized oracle. Synthetic experiments validate the theoretical results on controlled systems with known graph skeletons. Experiments on real-world datasets with unknown graph skeletons demonstrate superior federated root cause analysis over baselines.

The hazard ratio from the Cox proportional hazards model is a ubiquitous summary of treatment effect. However, when hazards are non-proportional, the hazard ratio can lose a stable causal interpretation and become study-dependent because it effectively averages time-varying effects with weights determined by follow-up and censoring. We consider the average hazard (AH) as an alternative causal estimand: a population-level person-time event rate that remains well-defined and interpretable without assuming proportional hazards. Although AH can be estimated nonparametrically and regression-style adjustments have been proposed, existing approaches do not provide a general framework for flexible, high-dimensional nuisance estimation with valid sqrt{n} inference. We address this gap by developing a semiparametric, doubly robust framework for covariate-adjusted AH. We establish pathwise differentiability of AH in the nonparametric model, derive its efficient influence function, and construct cross-fitted, debiased estimators that leverage machine learning for nuisance estimation while retaining asymptotically normal, sqrt{n}-consistent inference under mild product-rate conditions. Simulations demonstrate that the proposed estimator achieves small bias and near-nominal confidence-interval coverage across proportional and non-proportional hazards settings, including crossing-hazards regimes where Cox-based summaries can be unstable. We illustrate practical utility in comparative effectiveness research by comparing immunotherapy regimens for advanced melanoma using SEER-Medicare linked data.

We derive novel concentration inequalities for the operator norm of the sum of self-adjoint operators that do not explicitly depend on the underlying dimension of the operator, but rather an intrinsic notion of it. Our analysis leads to tighter results (in terms of constants) and simplified proofs. Our results unify the current intrinsic-dimension and ambient-dimension inequalities under independence, strictly improving both categories of bounds (such as by Tropp and Minsker). We present a general master theorem that we instantiate to obtain specific sub-Gaussian, Hoeffding, Bernstein, Bennett, and sub-exponential type inequalities. We also establish widely applicable concentration bounds under martingale dependence that provide tighter control than existing results.

Random-restart heuristics are widely used in nonconvex optimization and equilibrium computation: practitioners run a local algorithm from many initial conditions and interpret repeated convergence to the same output as evidence that the result is robust, dominant, or even unique. Despite its widespread use, this reasoning is usually informal. We provide a probabilistic framework for interpreting restart evidence. We give broad, easy-to-verify sufficient conditions under which repeated runs of a solver can be treated as independent draws from a categorical distribution induced by random initial conditions. Within this framework, we develop Bayesian inference from repeated identical outputs. We derive posterior concentration rates for basin size and uniqueness. These rates demonstrate that uniqueness is inherently harder than learning basin size: posterior concentration for uniqueness is polynomial, whereas basin size concentrates exponentially fast. We also provide a verification protocol for checking whether a given problem fits our framework. We demonstrate the protocol on a widely used equilibrium solver for mixed-logit demand with multi-product firms, and complement the verification exercise with posterior tables that apply to any restart experiment satisfying the protocol. We conclude by delineating limits of restart-based inference, including failures induced by solver–problem mismatch and limited visibility of alternative outcomes.

Quantifying the complexity of feed-forward neural networks (FFNNs) remains challenging due to their nonlinear, hierarchical structure and numerous parameters. We apply generalized degrees of freedom (GDF) to measure model complexity in FFNNs with binary outcomes, adapting the algorithm for discrete responses. We compare GDF with both the effective number of parameters derived via log-likelihood cross-validation and the null degrees of freedom of Landsittel et al. Through simulation studies and a real data analysis, we demonstrate that GDF provides a robust assessment of model complexity for neural network models, as it depends only on the sensitivity of fitted values to perturbations in the observed responses rather than on assumptions about the likelihood. In contrast, cross-validation-based estimates of model complexity and the null degrees of freedom rely on the correctness of the assumed likelihood and may exhibit substantial variability. We find that GDF, cross-validation-based measures, and null degrees of freedom yield similar assessments of model complexity only when the fitted model adequately represents the data-generating mechanism. These findings highlight GDF as a stable and broadly applicable measure of model complexity for neural networks in statistical modeling.

Computation in biological systems is fundamentally energy-constrained, yet standard theories of computation treat energy as freely available. Here, we argue that variational free energy minimization under a Poisson assumption offers a principled path toward an energy-aware theory of computation. Our key observation is that the Kullback-Leibler (KL) divergence term in the Poisson free energy objective becomes proportional to the prior firing rates of model neurons, yielding an emergent metabolic cost term that penalizes high baseline activity. This structure couples an abstract information-theoretic quantity – the *coding rate* – to a concrete biophysical variable – the *firing rate* – which enables a trade-off between coding fidelity and energy expenditure. Such a coupling arises naturally in the Poisson variational autoencoder (P-VAE) – a brain-inspired generative model that encodes inputs as discrete spike counts and recovers a spiking form of *sparse coding* as a special case – but is absent from standard Gaussian VAEs. To demonstrate that this metabolic cost structure is unique to the Poisson formulation, we compare the P-VAE against Grelu-VAE, a Gaussian VAE with ReLU rectification applied to latent samples, which controls for the non-negativity constraint. Across a systematic sweep of the KL term weighting coefficient $\beta$ and latent dimensionality, we find that increasing $\beta$ monotonically increases sparsity and reduces average spiking activity in the P-VAE. In contrast, Grelu-VAE representations remain unchanged, confirming that the effect is specific to Poisson statistics rather than a byproduct of non-negative representations. These results establish Poisson variational inference as a promising foundation for a resource-constrained theory of computation.

We develop a worst-case complexity theory for stochastically preconditioned stochastic gradient descent (SPSGD) and its accelerated variants under heavy-tailed noise, a setting that encompasses widely used adaptive methods such as Adam, RMSProp, and Shampoo. We assume the stochastic gradient noise has a finite $p$-th moment for some $p \in (1,2]$, and measure convergence after $T$ iterations. While clipping and normalization are parallel tools for stabilizing training of SGD under heavy-tailed noise, there is a fundamental separation in their worst-case properties in stochastically preconditioned settings. We demonstrate that normalization guarantees convergence to a first-order stationary point at rate $\mathcal{O}(T^{-\frac{p-1}{3p-2}})$ when problem parameters are known, and $\mathcal{O}(T^{-\frac{p-1}{2p}})$ when problem parameters are unknown, matching the optimal rates for normalized SGD, respectively. In contrast, we prove that clipping may fail to converge in the worst case due to the statistical dependence between the stochastic preconditioner and the gradient estimates. To enable the analysis, we develop a novel vector-valued Burkholder-type inequality that may be of independent interest. These results provide a theoretical explanation for the empirical preference for normalization over clipping in large-scale model training.

This paper proposes strategies for designing a system whose computational model is subject to aleatory and epistemic uncertainty. Aleatory variables, which are caused by randomness in physical parameters, are draws from a possibly unknown distribution; whereas epistemic variables, which are caused by ignorance in the value of fixed parameters, are free to take any value in a bounded set. Chance-constrained formulations enforcing the system requirements at a finite number of realizations of the uncertain parameters are proposed. These formulations trade off a lower objective value against a reduced robustness by eliminating an optimally chosen subset of such realizations. Risk-aware designs are obtained by accounting for the severity of the requirement violations resulting from this elimination process. Furthermore, we propose a computationally efficient design approach in which the training dataset is sequentially updated according to the results of high-fidelity reliability analyses of suboptimal designs. Robustness is evaluated by using Monte Carlo analysis and Robust Scenario Theory, with the latter approach accounting for the infinitely many values that the epistemic variables can take.

A key challenge in probabilistic regression is ensuring that predictive distributions accurately reflect true empirical uncertainty. Minimizing overall prediction error often encourages models to prioritize informativeness over calibration, producing narrow but overconfident predictions. However, in safety-critical settings, trustworthy uncertainty estimates are often more valuable than narrow intervals. Realizing the problem, several recent works have focused on post-hoc corrections; however, existing methods either rely on weak notions of calibration (such as PIT uniformity) or impose restrictive parametric assumptions on the nature of the error. To address these limitations, we propose a novel nonparametric re-calibration algorithm based on conditional kernel mean embeddings, capable of correcting calibration error without restrictive modeling assumptions. For efficient inference with real-valued targets, we introduce a novel characteristic kernel over distributions that can be evaluated in $\mathcal{O}(n \log n)$ time for empirical distributions of size $n$. We demonstrate that our method consistently outperforms prior re-calibration approaches across a diverse set of regression benchmarks and model classes.

Unbinned likelihood fits maximize the information extracted from experimental data, yet their application in realistic high-dimensional analyses has been fundamentally bottlenecked by the prohibitive computational cost of profiling systematic uncertainties. Furthermore, current machine learning-based inference methods typically estimate scalar parameters, discarding complex high-dimensional correlations. To address this, we propose a general Simulation-Based Inference (SBI) framework that elevates the fit target from scalar parameters to a multivariate Distribution of Interest (DoI), a learnable, invertible transformation of the feature space. We employ Factorizable Normalizing Flows to model systematic variations as parametric deformations, preserving tractability without combinatorial explosion. Crucially, we develop an amortized training strategy that learns the conditional dependence of the DoI on nuisance parameters in a single optimization process, bypassing repetitive training during likelihood scans. To capture the finite-sample statistical variance of the neural network DoI, we introduce a Poisson-bootstrap ensemble, which we marginalize through an averaged likelihood to deliver a complete statistical-plus-systematic uncertainty budget within a single unbinned likelihood. Validated on a synthetic dataset emulating a high-energy physics measurement, our method demonstrates that rigorous, fully profiled unbinned measurements can now be extended to complete differential distributions. By turning the fit into a functional measurement, this approach offers a powerful, unifying framework for a broad range of tasks conventionally treated as distinct problems, from detector calibration and differential cross-sections to unfolding and continuous parameter estimation.

Finite-state mean-field games (MFGs) arise as limits of large interacting particle systems and are governed by an MFG system, a coupled forward-backward differential equation consisting of a forward Kolmogorov-Fokker-Planck (KFP) equation describing the population distribution and a backward Hamilton-Jacobi-Bellman (HJB) equation defining the value function. Solving MFG systems efficiently is challenging, with the structure of each system depending on an initial distribution of players and the terminal cost of the game. We propose an operator learning framework that solves parametric families of MFGs, enabling generalization without retraining for new initial distributions and terminal costs. We provide theoretical guarantees on the approximation error, parametric complexity, and generalization performance of our method, based on a novel regularity result for an appropriately defined flow map corresponding to an MFG system. We demonstrate empirically that our framework achieves accurate approximation for two representative instances of MFGs: a cybersecurity example and a high-dimensional quadratic model commonly used as a benchmark for numerical methods for MFGs.

We propose a new model and estimation framework for spatiotemporal streamflow exceedances above a threshold that flexibly captures asymptotic dependence and independence in the tail of the distribution. We model streamflow using a mixture of processes with spatial, temporal and spatiotemporal asymptotic dependence regimes. A censoring mechanism allows us to use only observations above a threshold to estimate marginal and joint probabilities of extreme events. As the likelihood is intractable, we use simulation-based inference powered by random forests to estimate model parameters from summary statistics of the data. Simulations and modeling of streamflow data from the U.S. Geological Survey illustrate the feasibility and practicality of our approach.

Neural networks, particularly message-passing neural networks (MPNNs), are increasingly used as heuristics for hard combinatorial optimization problems. Yet many learning-based methods rely on supervision, reinforcement learning, or gradient estimators, causing high computational cost, unstable training, or limited guarantees. Classical approximation algorithms provide worst-case guarantees but are non-differentiable and cannot adapt to structure in natural input distributions. We study this tradeoff through Uniform Facility Location (UniFL), a problem with applications in clustering, summarization, logistics, and supply chains. We propose a fully differentiable MPNN that incorporates approximation-algorithmic principles without solver supervision or discrete relaxations. The model has provable approximation guarantees and empirically improves on standard approximation algorithms, narrowing the gap to integer linear programming.

We develop methodology to detect structural breaks in the slope function of a concurrent functional linear regression model for functional time series in $C[0,1]$. Our test is based on a CUSUM process of regressor-weighted OLS residual functions. To accommodate both global and local changes, we propose $L^2$- and sup-norm versions, with the sup-norm particularly sensitive to spike-like changes. Under Hölder regularity and weak dependence conditions, we establish a functional strong invariance principle, derive the asymptotic null distribution, and show that the resulting tests are consistent against a broad class of alternatives with breaks in the slope function. Simulation studies illustrate finite-sample size and power. We apply the method to sports data obtained via body-worn sensors from running athletes, focusing on hip and knee joint-angle trajectories recorded during a fatiguing run. As fatigue accumulates, runners adapt their movement patterns, and sufficiently pronounced adjustments are expected to appear as a change point in the regression relationship. In this manner, we illustrate how the proposed tests support interpretable inference for biomechanical functional time series.

Vanilla gradient methods are often highly sensitive to the choice of stepsize, which typically requires manual tuning. Adaptive methods alleviate this issue, with AdaGrad being a prominent example. In this work, we first show that, in composite convex optimization, AdaGrad may fail to exploit the Hölder-smoothness of the differentiable term, as its stepsizes may decrease too rapidly, leading to convergence rates arbitrarily close to the nonsmooth $\mathcal O(n^{-1/2})$ rate. Motivated by this limitation, we study AdaGrad-Diff, a block-wise extension of the UniSgd method considered by \cite{rodomanov2024universality}, which accumulates successive gradient differences rather than gradients themselves. In the deterministic Lipschitz-smooth setting, we prove summability of squared gradient differences without assuming bounded iterates, yielding weak convergence to a minimizer in Hilbert spaces. Moreover, we show that, under $\nu$-Hölder-smoothness, stepsizes decay no faster than $n^{-(1-\nu)/2}$, thereby avoiding AdaGrad's failure mechanism and recovering the standard $\mathcal O(n^{-(1+\nu)/2})$ rate. The same approach recovers known stochastic convergence rates under a unified condition covering bounded variance and interpolation. Numerical experiments demonstrate that the proposed method is more robust than AdaGrad in several practically relevant settings.