Papers

Adversarial training is one of the most reliable defenses against adversarial attacks, but its high computational cost must generally be paid anew for each task. Robust foundation models offer a promising alternative: adversarially pretrain a model once and then transfer its robustness to downstream tasks through lightweight adaptation. However, a fundamental question remains open: can robustness acquired during pretraining transfer to unseen tasks without further adversarial training? In this study, we answer this question affirmatively. A single model adversarially pretrained at scale can achieve optimal robustness on new tasks without additional task-specific training. Specifically, we show that, for a family of Gaussian-mixture classification tasks, a sufficiently deep linear transformer adversarially trained across tasks can asymptotically attain the robust Bayes error on previously unseen tasks through in-context learning from clean demonstrations. By contrast, a standardly trained model cannot. We further analyze convergence under gradient flow, an accuracy–robustness trade-off, and demonstration complexity.

We study the online calibration of multidimensional forecasts over an arbitrary convex set $Y\subseteq\mathbb{R}^d$ relative to an arbitrary error norm $\|\cdot\|_{L}$. For forecasting $d$ binary outcomes simultaneously ($Y=[0,1]^d$), we give the first algorithm that achieves $\varepsilon$-calibration in a number of rounds that is polynomial in $d$ for every fixed accuracy. It requires $d^{O(1/\varepsilon)}$ rounds, exponentially improving the dimension dependence of previous bounds. For multi-class forecasting ($Y=\Delta_d$), we obtain the same $d^{O(1/\varepsilon)}$ rate, improving the $d^{\widetilde{O}(1/\varepsilon^2)}$ bounds of Peng and Fishelson et al. Our algorithm is simple: on each round, it outputs a harmonically weighted distribution over harmonically smoothed past outcomes. The same algorithm works for every forecast set and norm. More generally, it achieves $\varepsilon$-calibration after $\exp(O(\gamma(Y,L)/\varepsilon))$ rounds, where $\gamma(Y,L)$ is a geometric parameter defined by a matrix discrepancy problem. The harmonic weights are motivated by the fact that the discrete Hilbert transform matrix achieves the optimal discrepancy up to a universal constant, simultaneously for every $L$. This optimality result may be of independent interest.

We study fair multi-armed bandits under the Nash Social Welfare (NSW) objective, which measures performance via the geometric mean of accumulated rewards. Existing work defines Nash regret as $\mathrm{NR}_T = \mu^\star - (\prod_{t=1}^T \mathbb{E}\mu_{I_t})^{1/T}$, where $\mu_{I_t}$ is the mean reward of the recommended arm $I_t$ and $T$ is the horizon. Since it applies the geometric mean to per-round marginal expectations, it ignores the joint distribution of rewards across rounds, leaving the NSW fairness motivation unaddressed at the trajectory level. We propose trajectory-wise Nash regret $\widetilde{\mathrm{NR}}_T = \mu^\star - \mathbb{E}[(\prod_{t=1}^T \mu_{I_t})^{1/T}]$, which computes the geometric mean over complete sample paths before taking expectations, capturing NSW fairness more faithfully. By Jensen's inequality, $\widetilde{\mathrm{NR}}_T \geq \mathrm{NR}_T$, making it a strictly stronger metric. We also introduce high probability Nash regret $\widehat{\mathrm{NR}}_T = \mu^\star - (\prod_t \mu_{I_t})^{1/T}$, giving the first high probability regret bounds in fair bandits. Our two-phase algorithm, Round Robin Nash Confidence Bound (RR-NCB), combines round robin exploration with a Nash confidence bound index policy. We show $\widetilde{\mathrm{NR}}_T \leq \widetilde{\mathcal{O}}(\sqrt{k\log T/T})$ and, with probability $1-\delta$, $\widehat{\mathrm{NR}}_T \leq \widetilde{\mathcal{O}}(\sqrt{k\log(kT/\delta)/T})$, matching the optimal $\widetilde{\mathcal{O}}(\sqrt{k/T})$ rate despite the stronger metrics. Optimality follows from a lower bound via AM-GM and standard $k$-armed bandit minimax arguments. Simulations validate our theory.

We study ridge regression from exactly $s$ distinct rows of a fixed design. Responses are fixed, and only the subset is random. The determinant law and selected ridge fit share one positive definite penalty. Established mean identities and exponential-family duality give the unique penalty that matches a prescribed full-data ridge fit in expectation. It exists exactly when $s$ exceeds the target's effective dimension. Our main result concerns centered covariance risk normalized by full-data penalized loss. For balanced signed coordinate replicas, a strict sector inequality gives the sharp risk and all maximizing responses at every budget from the dimension to one below the row count. This holds for any nonzero positive semidefinite query. With the target and query fixed, the maximizing response space is unchanged across these budgets. For general designs, we characterize attainment of a leave-one-out envelope. For existing real equiangular tight frames, flat row query energy characterizes when every nonzero residual response maximizes at two deletions. At three deletions, we give the sharp risk and complete maximizing space for isotropic queries, using unequal triangle weights. The balanced geometry yields a same-sample unbiased ridge–Horvitz–Thompson mixture with lower sharp risk and an exact mean-share improvement boundary. Under full recalibration after feature changes, we prove quadratic regret from searching the complete old maximizing space and a query-uniform bound on the mixture's risk gain. The strongest sector inequalities have exact computer-assisted proofs.

Transformers have exhibited impressive empirical success across various domains, but their theoretical foundations remain less developed. This work constitutes a mathematical study of the measure-to-measure operators defined by transformers. We show that transformers map sub-Gaussian inputs to sub-Gaussian outputs; this ensures that taking arbitrary-length compositions of the softmax operator is well-defined. We then show that transformers are Hölder continuous with respect to the 1-Wasserstein distance on appropriate spaces of sub-Gaussian inputs. This allows us to establish estimates on the error propagation along a transformer between a sub-Gaussian input and its empirical approximation. We also study a mean-field analog of the cross-attention mechanism, which is an operator from a pair of probability measures to a single probability measure. We show that cross-attention exhibits different Hölder regularity and sample-complexity in its two input arguments. Last, we apply our results to deduce approximation guarantees for measure-to-measure transformers. Together, these results provide a firm stability and finite-sample theory for transformers on sub-Gaussian data.

Discrete diffusion models have emerged as a practically successful framework for generative modeling on discrete product spaces, yet their statistical generalization properties remain poorly understood. Discrete real-world data such as text or biological sequences often concentrate on a small fraction of the astronomically large ambient space because of semantic or physical constraints, but existing bounds fail to capture this distributional structure and instead scale with the size of the ambient space, giving rise to almost vacuous error bounds. We address this gap for uniform discrete diffusion, one of the two dominant discrete diffusion paradigms alongside masking diffusion, by deriving statistical guarantees governed by the effective support size $s_n(P_0)$, a sample-size-dependent measure of distributional complexity. Given $n$ independent and identically distributed (i.i.d.) samples from an unknown data distribution $P_0$ on $[K]^d$, we show that, with appropriate choices of network size and hyperparameters, the expected total variation (TV) loss scales as $O(\sqrt{s_n(P_0)/n})$, while the expected Kullback–Leibler (KL) divergence is bounded by $O(\frac{1}{n}s_n(P_0)\log(eK^d/s_n(P_0))\log n)$. Furthermore, we show that the TV rate is minimax optimal and that the KL rate is minimax optimal up to a factor of $\log n$. Together, these upper and lower bounds show that uniform discrete diffusion successfully avoids the curse of dimensionality for distributions with small effective support size: the TV error rate depends on the ambient state-space size only through $s_n(P_0)$, while the corresponding KL rate incurs only an additional logarithmic dependence on the ambient state-space size.

Many real-world datasets exhibit unusually large values far more frequently than predicted by Gaussian models. Heavy-tailed distributions capture this behavior, yet evaluating learning performance under them remains challenging because rare, large feature entries retain non-vanishing effects even in high dimensions. Even in the canonical setting of empirical risk minimization for linear regression with entry-wise i.i.d. symmetric $\alpha$-stable data, a precise asymptotic characterization of prediction has been lacking. In this work, we introduce a functional order parameter that describes the random effective problem associated with each coefficient. Using the replica method, we fully characterize the generalization error in the proportional high-dimensional limit where the sample size and feature dimension diverge at a fixed ratio. Additionally, this analysis establishes a heavy-tail universality law, scaling laws relating typical errors to prediction reliability, and the Bayes-optimal prediction error. In addition to characterizing the effects of extreme entries on the learning process, our method applies broadly to other systems with persistent local heterogeneity.

Consider a collaborative decision-making setting in which team members with different reward structures need to make a joint decision. Repeated collaborative decisions can reflect individual differences in team members' negotiation ability. However, extracting such information is challenging because each observed response arises from a collaborative process involving multiple participants and their reward structures. In this paper, we propose a two-stage probabilistic measurement model for structured collaborative choice tasks. In the first stage, known member-specific rewards and latent participant negotiation traits jointly determine team-option utilities, from which a latent team choice is generated. In the second stage, individual responses are modeled conditional on the latent team choice, with possible deviations that depend on participants' own rewards. Participant-level covariates can also be incorporated through a structural model. A simulation study shows satisfactory recovery of the participant traits and model parameters across different sample sizes and item sizes. An application to data from simulated collaborative negotiation tasks shows choice prediction accuracy well above the random-choice baseline and meaningful associations between the estimated traits and external criterion variables. The proposed framework provides a psychometric approach to measuring participant-level negotiation traits with collaborative choice data.

Probability forecasts are calibrated when predicted probabilities match empirical outcome frequencies: among events assigned a probability $p$, we'd hope that the fraction of positive outcomes is close to $p$. We study the problem of sequential forecasting of binary outcomes. The classical $O(T^{2/3})$ bound on expected cumulative $\ell_1$-calibration error established by Foster and Vohra stood for over two decades until Dagan et al. reduced the exponent $2/3$ by an unspecified constant. We establish a new two-phase recursive labeling strategy for the sign-preservation-with-reuse game that yields the bound $O(n^{\alpha}t^\beta)$ for all choices of space and time. We then sharpen the reduction from upper bounds on sign preservation to calibration by modifying the equivalence of Dagan et al. to use only $O(\log T)$ instances of the sign-preservation-with-reuse game. This lets us establish an explicit bound of $O(T^{0.662942288})$, the first explicit exponent below $2/3$ for sequential calibration, by combining both improvements and choosing explicit feasible parameters.

Inverse inequalities are essential for controlling derivative growth and ensuring stable sampling in kernel approximation. Existing whole-space and single-scale estimates provide an established analytical basis, but their extension to multiscale analytic spaces on finite intervals requires uniform control of boundary effects and cross-scale cancellation. This paper establishes unweighted interval \(L^2\) inverse inequalities for spaces generated by primitives of kernels with positive exponential spectra. By retaining exterior diagonal energy and controlling off-diagonal interactions through one-sided estimates with endpoint-distance weights, we derive uniform coefficient bounds for normalized high derivatives. A synthesis-operator perturbation argument extends these bounds to positive spectral mixtures, while adjacent-derivative comparison and interval interpolation yield an inverse bound of order \(h_*^{-1}\), where \(h_*\) is the smallest scale. Under geometric scale separation and within-scale centre separation, the constant is independent of the number of scales and centres and their distances from the endpoints, allowing endpoint centres and coincident centres across scales. The coefficient lower bound approaches the optimal limit \(1/2\). Two-centre constructions establish the sharpness of the scale exponent for every fixed kernel in the class and, for a single exponential spectrum, a sharp deficit of order \(q^{-1/2}\) from \(1/2\) for the best uniform coefficient lower bound, where \(q\) is the derivative order. Numerical illustrations for Cauchy and logistic kernels demonstrate coefficient stability in several multiscale configurations. These results extend interval inverse estimates to multiscale analytic kernel spaces with endpoint centres and provide quantitative conditions for deterministic sampling stability.

Joint-event forecasts combine an estimated dependence model with current marginal forecasts. Retaining historical transforms at the marginal fits available when each outcome arrived creates overlapping estimation errors. We derive their joint influence with the terminal estimates in a stable Vine Copula VAR with normal innovation margins and a fixed, correctly specified Gaussian or positive Clayton vine. An intercept identity and the stable VAR filter reduce the historical correction to harmonically weighted innovation moments. Terminal slope uncertainty remains in the forecast gradient. A feasible covariance estimator yields asymptotically valid repeated-sample intervals for fixed one-sided event probabilities at the realized forecast state. A same-family terminal refit provides the efficiency comparison. Retaining dated transforms adds a positive semidefinite limiting covariance term, whose event-level cost depends on dependence sensitivity and forecast-state leverage. Monte Carlo comparisons show a small accuracy cost of retaining the transforms and a larger inferential role for current marginal estimation. Allowing Clayton dependence improves downside forecasts relative to Gaussian dependence under each transform rule. A fourteen-variable U.S. real-time forecasting application and simulations with twenty and fifty variables illustrate feasible joint-event estimation and analytic inference for fixed-depth vines.

Learning an ensemble of GFlowNets to sample from a discrete target distribution has become a common approach for achieving better state space exploration and convergence than that of a monolithic sampler. However, these methods often add a substantial runtime overhead to the base model, and their conceptual connection remains elusive. To address this, we first propose a general-purpose theoretical framework for describing a mixture of GFlowNets, which we specialize into continuously (CI) and discretely indexed (DI) collections. On the one hand, we show CI GFlowNets can be interpreted through the lens of a random features expansion, provably boosting the sampler's expressivity in graph-structured tasks and reducing learning instability via spectral shifting. On the other hand, we demonstrate DI GFlowNets encompass prior approaches for GFlowNet training and provide the foundation for the newly proposed Stratum-Conditioned (SC) GFlowNets. This method, which is inspired by the Doob's h-transform of Markov chains, decomposes the state space according to a prescribed modular function and restricts each component to sample from a distinct subset of it. Importantly, SC GFlowNets support centralized and component-wise embarrassingly parallel training, and we show both of them significantly speed up learning convergence and mode coverage without introducing any non-negligible extra computation.

Learning Gaussian mixture models (GMMs) using the Expectation-Maximization (EM) algorithm and its gradient-based variants is a fundamental problem in machine learning. It is known that randomly initialized (gradient) EM fails to learn multi-component GMMs in the exact-parameterized setting, where the number of components matches that of the ground-truth GMM. Recently, global convergence of gradient EM has been established in the over-parameterized setting, where more components are used, provided that the ground-truth components are well separated. In particular, the minimum separation between ground-truth components is required to scale as $\Omega(\sqrt{d})$, where $d$ is the dimension. In this paper, we show that this dimensional dependence is unavoidable in high-dimensional settings. Specifically, we consider a hybrid EM algorithm that uses standard EM updates for the mixing weights and gradient EM updates for the component means. For any $\epsilon > 0$, we prove that when the dimension is sufficiently large, in the worst case a separation of order $\Omega(d^{0.5-\epsilon})$ is insufficient to guarantee global convergence of population gradient EM in sub-exponential time under random initialization, even in the over-parameterized regime. Our result establishes an almost optimal worst-case lower bound on the ground-truth separation required for learning Gaussian mixtures via gradient EM in high dimensions.

Two populations of graphs often have to be compared without any correspondence between their vertices, for instance when networks come from different communities, or when a graph generative model is evaluated against held-out graphs. We study how many graphs such an unaligned two-sample test needs, and which graph statistics can detect which differences. For an Erdős–Rényi null and a planted two-block difference that leaves every expected degree unchanged, we show that $m\asymp t^{-3}$ graphs per group are necessary and sufficient when the per-graph signal-to-noise ratio is $t<1$. Signed triangle counts attain this rate, and the lower bound holds for every graph size. With aligned vertices $m\asymp t^{-1}$ graphs suffice, so misalignment costs a factor of order $t^{-2}$. When the triangle signal cancels, the rate becomes $t^{-4}$ and $4$-cycles are needed. Statistics built from trees have exactly the same expectation under both hypotheses, and tests based on finitely many of them have asymptotically no power. In the graphon limit, this class includes degree distributions and message-passing graph neural network features. For a non-constant null, a generic difference is visible at first order, and a simple motif test attains the aligned order of sample size, suggesting that misalignment is costly mainly for differences that are invisible at low orders. We also give an exactly valid test for one or two graphs per group, at a cost in power. In our simulations, the fitted exponents are close to the predicted ones, and degree-based and random-GNN evaluation metrics stay at their level in a setting where signed triangles need about $65$ graphs.

We develop a novel debiased machine learning (DML) estimator to analyze count data, which are common in clinical or biomedical studies. Specifically, we estimate the effect of a binary treatment on count outcomes using a partially linear Poisson model. A key advantage of this model is that it represents the complex covariate structure through a nuisance function, which is estimated flexibly using neural networks. We use a Neyman-orthogonal score function to construct the estimator, reducing its sensitivity to errors in estimating the nuisance function. Cross-fitting is used to mitigate overfitting bias in neural networks. Under mild regularity conditions, this DML estimator is asymptotically normal with root-$n$ convergence. We derive a closed-form variance estimator and construct a Wald confidence interval for the treatment effect. Extensive simulations demonstrate that the proposed estimator reduces bias and root mean square error compared with the generalized linear model and the augmented inverse probability weighting (AIPW) estimator in most settings, while achieving a coverage probability near the nominal level. The proposed procedure is implemented in the R package $\texttt{PoissonDML}$. We apply our approach to a synthetic HIV cohort data to investigate the effect of treatment on AIDS-defining event counts of 7,141 analyzed patients. Our approach provides practical and valid estimation and inference for treatment effects under the partially linear Poisson model while combining modern machine learning methods to estimate complex covariate structures.

Muon combines current and past gradients into matrix momentum. For $M=U\Sigma V^\top$, the idealized polar update $Q=UV^\top$ gives every singular direction the same weight. We refer to this as the flat profile. Several recent optimizers replace this flat profile with fine-grained spectral maps that give each direction its own gain. We ask how much of this spectral detail a Muon update needs. Our spectral diagnostics show that approximately $94$–$97\%$ of measured singular modes lie below an estimated noise edge, yet collectively align positively with a reference gradient. We introduce BulkBoost, a two-band spectral reweighting framework with fixed-rank and noise-calibrated variants. The latter uses split-minibatch gradient differences to calibrate a Marchenko–Pastur reference edge for Muon's Nesterov input, separating the bulk below the edge from the spikes above it. Both variants increase the bulk's relative weight through one shared gain while preserving the Frobenius norm of each matrix's unreweighted direction. For a fixed partition, our theory gives the first-order condition under which moving weight toward the bulk lowers the loss. It also quantifies the fraction of the maximal first-order improvement rate, over all per-mode reallocations, that two bands can capture. Across 30 continued-pretraining settings spanning Pythia-14M to 410M and six corpora, two-band reweighting is competitive with the fine-grained power-law profile of Freon and outperforms Spectra. Measured against Muon's flat profile, Freon reduces final loss by $0.022\%$ of the pre-adaptation loss on average, whereas the two-band variants achieve reductions of $0.073$–$0.147\%$. These observations suggest that useful departures from the flat profile are surprisingly low-dimensional: a single bulk-to-spike gain captures at least as much benefit as the fine-grained spectral profiles.

Data assimilation (DA) is the process of combining forecasts from a model with observations in order to optimally estimate the state of a system. This is critical for chaotic systems, such as the atmosphere, since if observations are not continually assimilated the model will quickly lose skill. DA is routinely performed (usually every 6 hours) at operational forecasting centres around the world. In this article we discuss the interface of machine learning (ML) and DA. This is still an emerging and quickly developing field, and this article tries to give an overview of some of the main topics and methods.

Combinatorial Thompson sampling (CTS) draws independent posterior samples for every arm, so its exploration dynamics ignore any relation among arms. We study a minimal change to those dynamics: an LLM is queried once for a partition of the arms, the partition becomes a positive-definite correlation matrix $\Sigma$ through an RBF kernel on cluster ranks, and the per-round posterior sample is drawn with covariance $\Sigma$ while the Beta posteriors are updated from real rewards only, so the LLM shapes how the sampler moves, not what it believes. We give a self-contained Bayesian regret bound for the idealized Gaussian sampler whose information gain splits into a $K\log T$ term from the $K$-cluster structure and a ridge term that grows to $d\log T$: the $\sqrt{d/K}$ improvement over independent sampling is a finite-horizon transient, exact only as the within-cluster correlation tends to one. The correlated sampler reduces regret by 19% over CTS on 16 synthetic Bernoulli families at $T=2{,}500$ (6-7% at $T=25{,}000$ with data-adaptive kernels) and by 41% on the Microsoft MIND-small news benchmark ($d=200$ real articles), while pseudo-observation warm starts give nothing. An LLM-free ablation with a simulated oracle of controlled quality shows that on unstructured instances the gain is a property of the kernel shape (a random partition, or a plain tempering of the sampling noise, reproduces it), while belief injection at matched oracle quality never helps.

Accurate forecasting of pathological outcomes is a central problem in psychology. To do so, psychologists often collect intensive longitudinal data. However, in such studies, the desire to acquire a large number of variables for the sake of accurate prediction is often counteracted by the need to minimize participant burden. Acquiring more variables per occasion can yield better predictions, but having too many acquisitions increase the risk of non-response and attrition. Longitudinal Active Feature Acquisition (LAFA) is a principled approach to resolve this conundrum. Instead of requiring responses to every item at every acquisition occasion, LAFA produces a policy that seeks to optimally select dynamic subsets of items to be acquired at each timepoint while preserving our ability to forecast a specific outcome. However, existing LAFA methods are mostly based on Neural Networks (NN) that are difficult to interpret in practice. In this work, we introduce a tree distillation method for learning an interpretable policy from NN-based LAFA networks. We validated our method through both a simulation and an empirical EMA dataset on forecasting daily alcohol consumption. In both cases, we find that we can meaningfully reduce the number of items acquired at each occasion with minimal loss in accuracy. Networks (NN) that are difficult to interpret in practice. In this work, we introduce a tree distillation method for learning an interpretable policy from NN-based LAFA networks. We validated our method through both a simulation and an empirical EMA dataset on forecasting daily alcohol consumption. In both cases, we find that we can meaningfully reduce the number of items acquired at each occasion with minimal loss in accuracy.

Wildfire regimes in Northern California exhibit strong spatio-temporal variability, heavy-tailed fire size distributions, and sensitivity to climatic conditions. We develop a hierarchical statistical framework to model wildfire occurrence, geometry, and burned area in Northern California from 1984-2023. Fire centroids are modeled as a spatio-temporal point process with covariate-driven intensity, extended by a multiplicative gamma-based stochastic shock to better capture the distribution of monthly fire counts. Conditional on fire location, burned area polygons are represented using a parametric ellipse model that separates scale, shape, and orientation, with ellipse parameters estimated from standardized fire perimeters. Total burned area is modeled using a heteroskedastic lognormal regression incorporating climatic and spatial covariates, including vapor pressure deficit anomalies. The framework enables Monte Carlo simulation of wildfire processes and provides a probabilistic tool for assessing wildfire risk and cumulative burned area under observed climate variability.

With the growing prevalence of streaming data, renewable estimation has become an important tool for updating statistical analyses without retaining historical raw observations. We study estimation and linear hypothesis testing in high-dimensional online M-estimation. Existing renewable procedures based on quadratic approximations can incur first-order approximation errors because the gradient of the unpenalized loss generally does not vanish at a preceding penalized estimator, posing a particular challenge for valid inference. To address this issue, we propose a corrected renewable loss that preserves historical first-order information and combine it with partial folded concave regularization and an iterative local linear approximation. We establish finite sample $\ell_1$- and $\ell_2$-error bounds, together with contraction bounds and strong oracle properties for both the constrained and unconstrained iterative estimators. Building on these results, we derive oracle Bahadur representations and establish asymptotically valid Wald and score tests for general linear hypotheses, allowing both the testing-oracle dimension and the number of restrictions to diverge. We further characterize their asymptotic power under local alternatives. The proposed procedure requires only recursively updated summaries while recovering the first-order inferential behavior of its pooled-data oracle counterpart. Simulation studies and real-data applications demonstrate the finite-sample performance of the proposed estimation and testing procedures.

Bayesian Optimization (BO) has become an established methodology for minimizing black-box functions of a vector input. Often, however, this parameter vector arises from the discretization of an inherently functional relationship. Several recent articles have considered the Functional Bayesian Optimization (FBO) setting, in which the variable to be optimized is not a member of a finite dimensional vector space, but rather an infinite dimensional function space. In this work, we propose $L^0$ Manifold Optimization (L0MO), a simple approach to FBO which searches the subset of a Reproducing Kernel Hilbert Space (RKHS) consisting of functions with a sparse representation in the kernel functions, optimizing both the kernel locations and their coefficients. We discuss in detail the relationship between our method and existing ones, providing a unifying lens through which to view prior works. To assess our method against the state of the art, we conduct an extensive computational study, and along the way develop a novel set of benchmark test functions which port standard finite-dimensional ones to the infinite dimensional domain. Our experiments demonstrate that, on balance, the proposed method achieves superior performance across a wide range of test benchmarks.

In actuarial science and quantitative finance, convex order provides a natural way to compare risks with the same mean. In dimension one, convex order is well understood through several characterisations. In higher dimensions, a natural approach is to compare all one-dimensional projections, but, although necessary, the resulting condition is in general not sufficient for multivariate convex order. A simple counterexample due to Pinelis exploits Popoviciu's inequality. We revisit this example using gauge functions, and show that gauge functions themselves do not suffice to characterise multivariate convex order for symmetric distributions, even in dimension two. We then show that the radial case makes the source of the failure transparent and construct radial counterexamples in every dimension at least two. Nevertheless, projection tests are sufficient for some important families of distributions. In particular, this is the case for linear transformations of a common radially distributed random vector. This framework includes elliptical distributions generated from a common radial distribution and, in particular, centred Gaussian distributions. Beyond comparisons between Gaussian distributions, Jourdain and Pagès recently studied the comparison of a Gaussian distribution with a Gaussian mixture, obtaining necessary conditions and sufficient conditions whose equivalence was left open in general. In this paper, we settle the comparison of a centred Gaussian distribution with a mixture of two centred Gaussian distributions: in every finite dimension, the Gaussian is dominated by the mixture in multivariate convex order if and only if all one-dimensional projections are ordered in convex order.

This paper studies simultaneous change point detection in the first- and second-order structures of high-dimensional matrix-valued time series. Specifically, we extend the framework of main effect factor models to accommodate multiple change points so that shifts in the mean and covariance structures are characterized by changes in the trend-stationary main effects main effects and in the factor-driven common components, respectively. Under this framework, the two types of changes are not necessarily aligned and can be arbitrarily large without masking the effects of each other, and each of the detected change points is identified either with the row or column categorizations which improves their interpretability. We propose a unified procedure based on moving-sum statistics to detect and localize both types of changes. Under general regularity conditions that allow for temporal and cross-sectional dependence, we derive detection and localization guarantees for the proposed method. The finite-sample performance of the method is demonstrated through extensive experiments, and its practical usefulness is illustrated by a real data application to New York City Yellow Taxi trip records.

We present a new perspective on the problem of likelihood approximation within the framework of simulation-based inference that promotes scalable and controllable simulation routines for large-scale data analysis, allows efficient parameter space exploration or smooth interpolation in high-dimensions and, thus, supports valid statistical treatments of hypothesis testings as well as uncertainty quantification. In particular, we consider a chain of $n$-aggregated normalizing flows for likelihood approximation scheme, where a set of upfront replicated observation datasets from the forward complex simulation model pass through the first set of bijective transformations, and then subsequently pass to the other sets of bijective transformations. Here, we assume that, for any $k \in \{1,\,2, \ldots, n\}$, the parameters corresponding to the first $k$ sets of bijective transformations are estimated sequentially, in some sense of optimality, for constructing flexible probability distributions, regardless of the remaining $(n-k)$ sets of bijective transformations. Moreover, our objects of interest are to highlight two complementary mathematical arguments that leverage an informatics-theoretic formalization, based-on empirical likelihood estimators under moment restrictions, and a sequential decision-making paradigm, with mixing distributions, for updating and aggregating the estimated parameters of the overall normalizing flows. As a by-product, the framework provides a reliable surrogate model, conditioned on the model parameters defining the forward computational simulation, that allows samples generation, with statistical powers, and facilitates computationally tractable scheme in the Bayesian paradigm for inference, hypothesis testings and uncertainty quantification.

Joint Models simultaneously model longitudinal and survival outcomes, leveraging patterns in patients' longitudinal trajectory to improve the prediction of survival outcomes. The classical parametric joint models, however, rely on fixed parametric assumptions, making them susceptible to bias under model misspecification and smaller sample sizes. We propose a deep joint model, DeepAJM, that does not require any parametric assumptions, while retaining a partially interpretable, per-longitudinal-outcome association structure. The joint model uses an encoder-decoder (sequence-to-sequence) architecture to learn the latent structure in patients' time-varying covariate trajectories. The model links the longitudinal processes to the survival processes through a learned interpretable association structure, in which each longitudinal output from the decoder gets remodulated by baseline covariates before it contributes to the risk scores from the survival head of the architecture. The model was evaluated on three datasets ( a cardiovascular-disease EHR cohort, a primary biliary cirrhosis (PBC2) dataset, and a simulated dataset) against a classical parametric joint model, TransformerJM, DA-LSTM and a Cox-based survival-only model. All models were assessed using C-index, integrated brier score (IBS), time-dependent AUROC, and time-dependent AUPRC. Our model achieved the best discrimination in terms of the C-index, time-dependent AUROC, and AUPRC across all datasets.

Synthetic patient data generation is a promising solution to the dual challenge of data scarcity and privacy constraints in healthcare machine learning. Realistic synthesis of patient-level clinical data requires jointly modeling heterogeneous static covariates, irregularly sampled longitudinal trajectories, and informative observation times - three tightly coupled components in practice yet rarely addressed together. We propose HyperNSDE, a continuous-time generative model that conditions a latent Neural SDE on static patient representations through a hypernetwork, allowing baseline characteristics to shape trajectory evolution beyond the initial condition without requiring a trajectory encoder, while stochastic latent dynamics capture realistic variability in generated paths. Observation times are modeled jointly through a latent-state-dependent intensity process, and training on irregular stochastic paths is stabilized via a deterministic-stochastic path decomposition with a non-adversarial signature-kernel objective. Experiments on simulated and real clinical datasets show improved observation-time fidelity and competitive performance, while matched-grid analyses reveal that forecasting and correlation metrics are affected by observation-grid regularity and trajectory smoothness.

We present a methodology for analyzing node-level outcomes while experimenting with edge-level treatments in a population connected by an undirected graph. Under our design, nodes are randomly assigned to test or control, and each edge inherits the treatment of its endpoints, with conflicts resolved by randomization. We use each node's assigned status as an instrument for its treatment exposure. We show that the Wald estimator is consistent for the global average treatment effect (GATE), even when the edge weights used to construct the exposure are misspecified. We formalize the assumptions needed both in terms of the linearity of potential outcomes and the sparsity of the graph, and prove the asymptotic normality of the Wald estimator. The estimator is straightforward to implement, requiring no assignment simulations over the graph. Our Monte Carlo study demonstrates the strong performance of our approach in the context of a social platform.

The goal of the partial information decomposition (PID) is to quantify the redundant and synergistic information that multiple sources provide about a target. PID has many applications in machine learning, neuroscience, and other fields, but defining and computing it for high-dimensional continuous systems remains challenging. Here, we define a PID for multivariate Gaussian systems based on the Blackwell order, which formalizes when one channel is more informative than another. We prove that Gaussian channels are optimal for extracting both redundant and union information, yielding an intuitive geometric interpretation and an efficient numerical algorithm for the PID. Our union information and synergy coincide with the well-known BROJA measures, and we derive closed-form expressions for both in the case of two sources. We also argue that Blackwell redundancy (which differs from BROJA) is the only existing redundancy measure that satisfies a set of natural desiderata. We demonstrate the scalability of our method on systems with up to a thousand dimensions or sources. Our approach is illustrated on an optimal control problem, where it identifies redundant and synergistic interactions between sensor and memory.

Estimating extreme conditional quantiles faces two major challenges: understanding complex nonlinear relationships between variables and accurately extrapolating into the tails of the distribution, where data are sparse. To address both issues, we introduce a novel approach that combines the theoretical framework of block maxima with the predictive power of generalized random forests. The conditional distribution of maxima is modeled by a generalized extreme value (GEV) distribution whose parameters depend on the covariates and are estimated via a weighted maximum likelihood procedure, with weights derived from generalized random forests to capture complex high-dimensional structures. The conditional quantile estimator follows from the inversion of the estimated GEV distribution function. We establish theoretical results ensuring the existence and consistency of the proposed weighted maximum likelihood estimator. The method provides a robust and flexible framework for extreme quantile regression. The performance of the proposed approach is illustrated using simulated data, as well as an application to financial portfolio losses from stocks listed on the NYSE, AMEX, and NASDAQ.