Papers

Xu, Vardi and Safran (ICML 2026) prove that over-parameterized ridge regression over an unstructured random Gaussian feature map groks, with the delay between memorization and generalization growing as $1/\lambda$ in the weight decay. We show that on a structured feature map the same delay does not appear. For a band-limited Fourier feature map over $\mathbb{Z}_p^2$ carrying a single-character target that lies inside the expressible class, enlarging the band at fixed positive weight decay drives peak held-out accuracy monotonically from $1.00$ to $0.07$, with no memorize-then-generalize regime anywhere along the sweep. The degradation is not an interpolation effect. It sets in at capacity ratio $q/n = 0.638$, far below the interpolation threshold, on separate grounds from the exact null space that appears above it. What does have a sharp boundary is the active support. Holding the nominal dimension fixed and masking the band back to $1089$ active modes restores held-out accuracy of $1.000$ with zero variance across seeds, while the full $4225$-mode band collapses to $0.185$. The number of active modes acts through the teacher-weighted spectrum of the empirical Gram matrix and not through the capacity ratio, which makes this a statement about feature geometry and not a restatement of double descent.

Data attribution seeks to quantify how individual training examples shape a model's predictions and underpins problems including data valuation, machine unlearning, and model interpretability. Despite having a long line of work, computationally scalable methods often struggle to predict the effect of removing training data in neural networks due to their non-convex nature. To overcome this challenge, metagradient-based methods such as MAGIC (Ilyas and Engstrom, 2025) differentiate each prediction through the entire training run and compute its exact influence with respect to the training data, but require a separate run for every prediction. To reduce this cost, we cast budgeted attribution as estimating a large influence matrix from a small number of measurements. We show that the measurements most appropriate for recovering this matrix differ from those best suited for attribution itself. We then present two algorithms, MAGE and SPELL, suited for reconstruction and attribution respectively, that run on existing metagradient machinery at no extra cost. Empirical studies demonstrate strong performance over existing baselines across training scales and measurement budgets.

This paper aims to develop new and efficient distributed algorithms for solving a class of monotone inclusions, $0 \in \sum_{i=1}^n (G_ix + T_ix)$, over a connected network of $n$ agents, where the single-valued operator $G_i$ and the possibly multivalued operator $T_i$ remain private to agent $i$. Existing distributed algorithms for this problem class are primarily non-accelerated, and their exact convergence rates in the original primal space are largely unexplored. To bridge this gap, we propose two Decentralized Fast Fixed-Point-based algorithms, ND-DFFP and NI-DFFP, which integrate Nesterov-type acceleration with primal-dual techniques under two prominent settings: (i) \textit{Lipschitz continuity of $G_i$ and maximal monotonicity of $G_i+T_i$}; and (ii) \textit{co-coercivity of $G_i$ and maximal monotonicity of $T_i$}. While ND-DFFP utilizes a homogeneous network-dependent stepsize, NI-DFFP reformulates the problem into a three-operator inclusion to decouple the network topology, enabling heterogeneous network-independent stepsizes. Under appropriate assumptions, we establish an $\mathcal{O}(1/k)$ convergence rate for the consensus error and an $\mathcal{O}(1/k)$ rate for both the restricted gap function and the squared forward-backward splitting residual, with the latter two metrics evaluated at the network-average iterate or its projection onto the effective domain. Finally, numerical experiments on distributed bilinear matrix games and a virtual power plant problem demonstrate the competitive performance and computational efficiency of our methods over recent decentralized baselines in the literature.

For constant-stepsize stochastic approximation (SA), the iterates converge in distribution to a stationary law that depends on the stepsize $\alpha.$ Steady-state convergence (SSC) concerns the limit of the scaled stationary distribution as $\alpha \downarrow 0.$ Existing SSC theory requires i.i.d. or additive noise and global differentiability of the mean operator, and yields suboptimal rates. We develop a unified SSC theory for constant-stepsize contractive SA driven by Markovian, multiplicative noise, covering both locally differentiable and locally nondifferentiable mean operators. A key methodological contribution is a multi-step universality framework that progressively reduces the original stochastic recursion to tractable auxiliary dynamics while preserving its steady-state limit. Under local quadratic linearization at the fixed point, we obtain a Gaussian approximation of the scaled steady state at the optimal rate $O(\sqrt{\alpha})$ in Wasserstein-2 distance, which further gives finite-time Gaussian approximations for the raw iterates. In the locally nondifferentiable regime, we establish a general SSC result and show that the leading-order asymptotic bias can be of order $\sqrt{\alpha}$, in contrast to the $\alpha$-order bias in the smooth regime. We apply the theory to Markovian linear SA and asynchronous Q-learning, neither of which is covered by prior results. We further propose a bias-reduction scheme for Q-learning that requires no knowledge of the local smoothness regime, validated by numerical experiments.

Shapley value (SV)-based methods are the prevailing framework for feature attribution in machine learning, yet existing population-level Shapley estimators generally assume that observations used to evaluate the coalitional game are fully observed under a common feature space. This assumption is routinely violated in multi-site studies across biomedicine, social science, and environmental monitoring, where institutions record different features under different protocols, producing systematic blockwise missingness across sources. We first show that the standard remedy of imputing missing features before computing Shapley values introduces systematic, coalition-dependent bias into the resulting attributions. We then propose FUSHAP (Fusion Shapley Attribution from Partially-observed data), a method that leverages partially-observed auxiliary sites to reduce the variance of a preliminary single-site Shapley estimate without imputation. A permutation-based screening step detects and excludes sites whose data distributions are incompatible with the target population. In synthetic experiments, FUSHAP achieves $3$–$8\times$ lower MSE than the single-site estimator and $2$–$3\times$ lower MSE than imputation baselines without incurring imputation-induced bias, and the screening procedure identifies misaligned sites with $82\%$ power at moderate misalignment and $100\%$ for strong misalignment. On multi-site air quality and multi-center clinical data, FUSHAP reduces MSE by approximately $3$–$7\times$ relative to the single-site estimator; in the clinical application, standard imputation can increase MSE above the single-site baseline.

This paper introduces a novel framework for Regularized Multivariate Functional Principal Component Analysis (ReMFPCA) via Functional Singular Value Decomposition (SVD). The proposed method extends existing MFPCA approaches by incorporating a generalized functional SVD within a Hilbert space framework, enabling simultaneous regularization of both functional principal components (PCs) and their associated PC scores. A key innovation of this framework is the inclusion of a sparsity penalty on the PC scores, which enhances interpretability by filtering out irrelevant subject-specific variations. This dual-penalization strategy represents a significant advancement beyond existing covariance-based eigen decomposition methods, which penalize only the functional PCs. Two power algorithm implementations, sequential and joint, are proposed, together with a cross-validation approach based on iterative regression for optimal smoothing parameter selection. Comprehensive simulation studies and real data applications demonstrate that the proposed framework substantially improves the extraction of informative and interpretable components, offering methodological and practical benefits for analyzing multivariate functional data across diverse domains.

Importance weights are essential in domain adaptation under label shift, yet their utility is often undermined by the finite sample uncertainty associated with their estimation. Existing methods typically analyze this uncertainty through Gaussian elimination on interval-valued linear systems, which leads to overly conservative confidence regions and inefficient downstream applications. We propose a paradigm shift from inversion-based inference to a direct matrix constraint framework. We use this framework to define a joint confidence region and extract marginal intervals via linear programming, deriving provably tighter bounds for importance weights while maintaining exact finite-sample validity. Furthermore, we analyze the confidence region's geometry and provide the theoretical results for its diameter bounds. Evaluated across text, image, multimodal benchmarks, including AGNews, MNIST, CIFAR-10, N24News, and a real-world autonomous driving dataset, nuImages, our approach consistently yields shorter confidence intervals and smaller prediction sets than inversion-based methods.

This study addresses the PU classification problem under violations of the SCAR assumption. We investigate logistic regression-based approaches, namely the cluster method and its extensions with strict and non-strict Lasso regularization. The primary contribution of this work is the integration of the SMOTE technique to alleviate class imbalance and systematically assess its impact on the performance of the considered algorithms. SMOTE is first applied to rebalance the training dataset. Next, cleaning labels are derived via 2-means clustering. Logistic regression is then trained on the cleaned data, where identified positive instances are augmented with additional true positives and the remaining observations are treated as negative. The experimental evaluation is conducted on 13 real benchmark datasets and one synthetic dataset. For comparison, we include the naive approach and the Spy-EM method. The results demonstrate that incorporating SMOTE improves classification performance when the SCAR condition is violated and indicate moderate robustness of the LassoJoint method in this setting.

Training-free diffusion inverse solvers typically choose between local measurement guidance and costly clean-space posterior updates. Independent posterior refresh can improve global correction by sampling a clean conditional and re-noising it, but its practical realization requires probability-flow ODE integration and clean-space Markov chain Monte Carlo (MCMC). We propose Direct Conditional Transition Sampling (DCTS), a direct stochastic-flow approximation to the same ideal refresh target. Rather than explicitly drawing a clean sample, DCTS estimates the measurement-conditioned clean mean along a short inner path and transports Gaussian source noise directly to the next noisy state. A denoiser-compatible sufficient statistic and a covariance-scaled operator update enable this conditional-mean estimation. Experiments on four inverse problems demonstrate that DCTS achieves competitive reconstruction quality with up to $16.8\times$ speedups over competing methods.

Constraining reaction rate coefficients is a central challenge in the development of explicit atmospheric chemical mechanisms, particularly for autoxidation systems where many reaction pathways are only indirectly observed through high-resolution mass spectrometry. In this study, we evaluate rate-coefficient optimisation methods for a toy-case autoxidation mechanism using synthetic data with known ground truth. Two complementary approaches are compared: ODE-constrained neural-network optimisation, which provides efficient point estimates of uncertain rate coefficients, and the Markov Chain Monte Carlo (MCMC) approach, which samples the posterior distribution of rate coefficients and quantifies parameter uncertainty. The methods are tested using direct concentration observations and mass-spectral observations under different noise levels. For unperturbed and low-noise synthetic observations, both methods converged towards the known rate coefficients, with the neural-network optimiser providing faster point estimates. Under high-noise conditions (with the signal-to-noise ratio approximately S / N = 1), however, MCMC was substantially more robust in recovering the rate coefficients. The posterior analysis shows that mass-spectral aggregation broadens credible intervals even at low noise, and that high-noise mass spectra can leave many individual reaction rates weakly identifiable. Posterior predictive validation nevertheless shows how broad parameter uncertainty constrained by MCMC remains consistent with accurate reproduction of the observable mass spectrum. These results demonstrate that point-estimation and Bayesian sampling methods provide complementary information: neural-network optimisation is effective for informative data, whereas MCMC is essential for diagnosing uncertainty, non-uniqueness, and identifiability in noisy or aggregated inverse problems.

Many modern prediction tasks involve data from multiple heterogeneous sources, while the test distribution may differ substantially from any individual source. Although heterogeneity poses challenges, it also offers an opportunity: different sources may provide complementary information, with some regions of the feature space better represented in one source than another. We propose Multi-Source Randomly Localized Conformal Prediction (MS-RLCP), which builds on the local coverage properties of randomly localized conformal prediction (RLCP) (Hore and Barber, 2025) and extends it to multiple sources through data-adaptive source selection. Under the widely adopted assumption of a shared response distribution conditional on the features across sources and the test population, we establish finite-sample coverage bounds using an interpretable notion of envelope distribution that captures their aggregate feature-space representation. Our analysis allows the test feature distribution to be absolutely continuous with respect to the envelope, extending beyond mixtures of source distributions. Under additional regularity conditions, we also establish asymptotic test-conditional coverage. Simulations and real-world experiments demonstrate the effectiveness of MS-RLCP across varying levels of data heterogeneity.

Variance penalization is a principled approach to risk-sensitive reinforcement learning (RL) that explicitly trades expected return for policy stability. Existing methods require a dedicated second critic to estimate return variance online, adding architectural complexity and compounding estimation error during learning. We propose a nonparametric variance-penalized actor-critic (VPAC) framework that replaces the variance critic with statistically grounded online estimators based on bootstrapping and random scaling, techniques drawn from the statistical inference literature for stochastic approximation. These estimators require no auxiliary network, maintain a single-critic architecture, and produce variance penalties that are bounded by construction, enabling clean convergence analysis. We establish almost-sure convergence for both a variance-penalized Q-learning algorithm and a two-timescale actor-critic variant via the ordinary differential equation (ODE) method, requiring only that variance estimates remain bounded rather than consistent. Empirically, we evaluate across discrete and continuous stochastic environments, demonstrating that the proposed methods match or exceed the variance reduction achieved by the existing dual-critic VPAC baseline while eliminating the overhead of a second critic. We further validate on a high-temperature superconductor (HTS) manufacturing case study, where VPAC-RS (Random Scaling) achieves a 74% reduction in steady-state critical current variability and a 63% reduction in episode return standard deviation, translating directly to improved yield consistency. Our results establish nonparametric statistical inference as a practical and theoretically sound alternative to auxiliary critics for risk-sensitive RL.

Distributed multiple testing asks $N$ sites, each holding p-values for its own hypotheses, to control the false discovery rate (FDR) of the discoveries made across the whole network while communicating only a small number of bits. Greedy interval aggregation (Pournaderi and Xiang, IEEE TSIPN, 2023) meets the communication budget but controls FDR only asymptotically, and its FDR can exceed the target at finite sample sizes; the same interval counts both select the rejection regions and calibrate the stopping rule, which biases the selected interval densities upward. We propose budgeted BONuS-GA: each site mixes synthetic uniform p-values into its data, the center ranks candidate p-value intervals across sites using the pooled counts, and each site's false discoveries are estimated from its own synthetic counts under a per-site share of the level. The procedure needs no knowledge of the null proportions, uses every p-value for both selection and inference, and controls $FDR\le\alpha$ at every sample size for any fixed assignment of hypotheses to sites, assuming only that the null p-values are independent uniforms, independent of the non-null p-values. It keeps the original $O(\sqrt m\log m)$ communication, $m$ being the total number of p-values, when $N=O(\sqrt m)$. We also give a sample-splitting alternative, cross-fit greedy aggregation, with finite-sample control under a random-site model, and we quantify the selection bias behind the original procedure's failure. In a monitoring-network simulation the budgeted procedure retains most of the power of its heuristic counterpart at moderate-to-large site sizes.

We study a class of distributionally robust optimization (DRO) problems for the statistical risk problem, formulated as minimax problems over the product of a Euclidean space and a Riemannian manifold. Because the resulting minimax landscape is nonconvex nonconcave in general, no globally convergent first order method is known to be available. We instead introduce the notion of a basin saddle point, a Nash equilibrium defined locally on the Cartesian product of a $\delta$ basin around a connected component of the local minima critical set and a geodesic ball on the measure manifold. We develop an abstract convergence framework for a Riemannian gradient ascent multistep descent iteration to a basin saddle point under a local Łojasiewicz type growth condition, with exponent $\beta \in (1,2]$, in the $\delta$ basin around connected components of the local minima critical sets. Under Lipschitz regularity of critical sets we establish linear convergence for $\beta = 2$ and polynomial convergence for $\beta \in (1,2)$ to a basin saddle point, with explicit dependence on the sectional curvature of the manifold. We then instantiate this framework for the statistical risk DRO problem over Gaussian measures, where the ambiguity set is naturally modeled as the product of Euclidean space and the Bures Wasserstein manifold of covariance matrices, which we relax to a penalized DRO formulation. We derive nonasymptotic Hessian estimates for the resulting Lagrangian, establish existence and local uniqueness of its maximizer, and prove that an alternating Riemannian gradient scheme converges to a basin saddle point of the penalized DRO problem, recovering the linear and polynomial rates of the abstract theory with all constants explicit in terms of data dimension, loss moments, and the reference covariance.

We derive exact local responses for attention interventions, allowing candidate edits to be scored from a cached baseline and one backward pass. The starting point is the RoPE derivative $\partial_p z(p) = A z(p)$: its integral gives the finite positional displacement, which we carry through the softmax without linearising either rotation or normalisation. The resulting predictions achieve 95.36–96.52% sign accuracy across 92,160 executed positional edits on 768 held-out prompt sets, reducing answer-margin MAE by 73.6–82.5% against the positional Jacobian and by 36.2–50.9% against zero. For simultaneous key and value edits, the same divided-difference calculus isolates the interaction term $C_{KV} = \sum_j (p'_j - p_j)\,\varepsilon_j$, which is omitted by adding separate attributions. Retaining it reduces downstream margin MAE by more than a factor of nine in every setting of a 5,120-intervention sweep across two Qwen sizes, two tasks, and multiple layers; reductions against a quadratic interaction correction are 75.9–98.5%. Exactness concerns the edited attention write; downstream predictions contract that response with a baseline gradient and are evaluated by native execution. The calculus also yields a KL certificate for local approximation error, an exact query-conditioned gradient-step representation whose curvature identifies attention-preserving query directions, and minimum-norm query control. Sparse evaluation supports candidate ranking and cache decisions under explicit local distortion criteria.

Variable importance scores from tree-based methods such as random forests favor continuous predictors over categorical ones. We present a theoretical analysis of this bias and propose a simple remedy: add a small amount of noise to each categorical predictor. The correction is demonstrated on a variety of simulated and real-world datasets and combined with integrated path stability selection to perform variable selection with false discovery control for mixed data.

Ensemble sampling offers a practical approach to randomized exploration by maintaining a collection of models, but how small an ensemble can be while retaining strong regret guarantees remains unresolved. In particular, the existing guarantees use an ensemble size of $\Theta(d\log T)$, leaving a logarithmic gap in the horizon $T$ relative to the intrinsic $\Omega(d)$ ensemble-size barrier. We aim to narrow this gap by proposing an ensemble sampling algorithm that refreshes the ensemble only when the regularized Gram matrix changes substantially. This mechanism localizes the perturbation analysis to epochs with controlled Gram-matrix drift and reduces the sufficient ensemble size to $\Theta(d\log d+d\log\log T)$, while preserving the state-of-the-art $\tilde O(d^{3/2}\sqrt T)$ regret for ensemble sampling with arbitrary bounded arm sets. We further show that, when the arm set is finite of cardinality $K$, the proposed algorithm achieves the sharper regret bound $\tilde O(d\sqrt{T\log K})$. To the best of our knowledge, this is the first ensemble-sampling guarantee that simultaneously recovers both canonical regret scalings known for randomized linear bandit algorithms: the $\tilde O(d^{3/2}\sqrt{T})$ rate for arbitrary bounded arm sets and the $\tilde O (d\sqrt{T\log K})$ rate for finite arm sets. The algorithm also admits an anytime implementation without resetting past data, and experiments show that it remains competitive with baselines while using substantially smaller ensembles.

This work studies the stability and convergence of augmented primal-dual dynamics when constraint values are estimated from samples. Unbiased constraint observations can produce a biased augmented multiplier signal, shifting the equilibria of the mean dynamics. For componentwise inequalities, we give a necessary and sufficient condition for preserving the Karush-Kuhn-Tucker (KKT) equilibria and construct a convex example with a locally exponentially stable equilibrium that violates complementarity. To address this bias, constraint values are estimated recursively before forming the augmented multiplier signal. For smooth convex conic problems, a joint energy analysis establishes boundedness of the primal, dual, and estimation states, vanishing estimation error, and almost sure convergence of the primal-dual iterates to a single KKT point under global regularity and bounded conditional second moments. The result allows nonunique solutions and multipliers while keeping the number of samples per iteration fixed. Numerical studies illustrate the predicted equilibrium bias and examine convergence with nonunique KKT points and nonlinear constraints.

Operator learning is formulated on function spaces, but training data are typically available only through finite-dimensional representations. In encoder–decoder architectures, a matrix-valued kernel on the encoded space induces an operator-valued kernel on the original function spaces, and the corresponding reproducing kernel Hilbert spaces are isometrically isomorphic. As the input and output resolutions increase, the induced kernels converge to a limiting kernel, in the sense of operator-norm convergence of their associated integral operators, allowing regularity assumptions to be stated independently of the encoding resolution. For regularized stochastic gradient descent, we establish upper bounds for decreasing and fixed step sizes, separating the encoding and regularization terms from optimization terms of order \(t^{-\theta}\) and \(T^{-\theta'}\), respectively, for any \(\theta,\theta'\in(0,1)\). We further prove lower bounds showing that these encoding-induced terms are generally unavoidable. The analysis is further extended to encoder–decoder neural networks through the limiting neural tangent kernel (NTK), yielding error bounds with an additional finite-width term and polynomial parameter and sample complexity guarantees when the encoding errors decay algebraically. The framework covers matrix-valued kernels constructed from radial and dot product kernels, NTKs arising from wide encoder–decoder neural networks, and encoder–decoder pairs based on Fourier, Legendre polynomial, wavelet, PCA, or pointwise sampling representations.

We establish the equivalence between the stochastic optimal control and path space formulations of the Schrödinger bridge problem (SBP) for the kinematic equation on a compact connected Lie group. Using the geometric concepts of horizontal lift and stochastic anti-development, we derive a Girsanov-type change-of-measure result, and show that the expected control energy equals the relative entropy of the controlled path law with respect to the reference Wiener measure. Thus, the SBP is equivalently a path space relative entropy minimization problem subject to prescribed endpoint marginals.Our result has three useful implications. From an analytic viewpoint, the shown equivalence helps prove the existence and uniqueness of the SB. From a probabilistic viewpoint, it helps interpret the SB as the most probable deviation of the uncontrolled stochastic dynamics consistent with the endpoint constraints. From a computational viewpoint, it allows using static Sinkhorn recursions to directly solve the relative entropy minimization problem and compute the optimal path measure. We illustrate the equivalence numerically on the torus $\mathbb{T}^2$. The code is publicly available at: https://github.com/gradslab/LargeDeviationSBP

On evolving graphs, node classifiers must satisfy two key requirements: inductive generalization to newly arriving nodes under distribution shift and calibrated uncertainty for safety-sensitive applications. Standard graph neural networks (GNNs) are typically trained once and address neither requirement. We adapt the Bayesian last-layer (BLL) model by placing random last-layer parameters on top of a deterministic GNN encoder for uncertainty quantification. The categorical softmax likelihood required for classification breaks Gaussian conjugacy, so neither the training posterior nor the test-time streaming update has a closed-form solution. To address both challenges, we introduce a variational Bayesian last-layer (VBLL) objective that jointly trains the encoder and an approximate last-layer posterior by maximizing an evidence lower bound with a Monte Carlo expected log-likelihood. At test time, we freeze the encoder and apply an online Laplace update to the last-layer posterior. This update corresponds to a power-prior Bayesian model with exponential forgetting and a Kullback-Leibler anchor to the training posterior. Across five node-classification benchmarks under distribution shift, online GVBLL is the only method to achieve the best accuracy and negative log-likelihood on every dataset. It improves accuracy by up to 17 percentage points on Cora and 14 percentage points on ogbn-arxiv over the strongest non-GVBLL baseline, while remaining competitive in calibration with MC Dropout, Deep Ensembles, Temperature Scaling, and Gaussian-process classifiers.

This work addresses decentralized online Riemannian optimization on Hadamard manifolds. Prior work under geodesic convexity (g-convexity) may require curvature information in the optimization analysis, typically through a finite lower bound on the sectional curvature. Curvature may also enter the step size or contraction factor of tangent-space Riemannian consensus schemes. In this work, we relax the curvature dependence for a narrower class of horospherical convex (h-convex) functions. We study Distributed Riemannian Online Gradient Descent (D-ROGD), which combines local Riemannian h-subgradient updates with an implicit Fréchet-mean consensus. For h-convex and strongly h-convex local objectives, we establish $O(\sqrt{T})$ and $O(\log T)$ static regret, respectively, matching the corresponding Euclidean rates with respect to $T$, with network dependence governed solely by the spectral gap. To our knowledge, these are the first curvature-independent regret guarantees for decentralized online optimization on Hadamard manifolds. Experiments on hyperbolic embeddings corroborate the predicted rates, with no observable degradation due to curvature.

Network analysis for multivariate time series is popular in many fields, from neuroscience to seismology. The inverse spectral density is a common choice for time series network analysis due to its representation of the frequency domain correlation between two variables after removing the best linear predictor of all other variables. In many applications, the goal is to study how these networks change across different conditions. For example, in neuroscience, one might be interested in how the brain connectivity network changes before and after stimulation. Towards this goal, we develop an inference framework based on a direct estimate of the difference in two high-dimensional inverse spectral densities. We develop a new Gaussian approximation error bound for any de-biased D-trace estimation procedure which is then leveraged to both inform optimal window sizes of Welch's estimators of the spectral density and establish asymptotic normality of our de-biased D-trace estimator. Moreover, we develop an efficient algorithm based on a generalized D-trace estimation procedure to overcome the computational complexity of high-dimensional inference. The method is illustrated on synthetic data experiments and on experiments with electroencephalography data.

We develop a causal framework for matrix completion under missing not at random (MNAR) data. Drawing on synthetic controls from the econometric panel data literature, our approach relaxes two assumptions common in MNAR matrix completion: positivity and independence of observation indicators. Unlike traditional panel data models, which often require prescribed block-sparse geometries, our framework accommodates flexible, heterogeneous observation patterns through target-specific local information structures. We propose synthetic nearest neighbors (SNN), a local synthetic-controls-inspired estimator, and establish finite-sample entrywise error bounds and consistency for mean recovery under suitable conditions. We further derive asymptotic normality under heteroskedastic noise and develop feasible entrywise inference. To estimate entry-specific noise variances, we apply the same local principle to squared outcomes, obtaining consistency under bounded noise and asymptotic unbiasedness under general subgaussian noise. Simulation studies corroborate the theoretical findings across a range of missingness designs and observation patterns.

Objectives in scientific machine learning are often prescribed as a sum of several terms, such as the residual, boundary, initial, and data losses of a physics-informed neural network. In the pullback-corrected scalar auxiliary variable (PB–SAV) method, one scalar tracks the shifted objective while the component gradients build a positive semidefinite curvature correction of rank at most the number of components. We carry that correction into an optimizer with momentum and an adaptive mobility, applying it to the gradient and the stored momentum in a single implicit solve. A mobility that is nonincreasing in the Loewner order yields an exact modified energy law, covering Euclidean and AMSGrad-type choices; the corresponding identity for momentum appended after the solve carries a cross term of indefinite sign. For a fixed mobility we give a necessary and sufficient condition for local stability at a stationary point, depending on the Hessian minus twice the correction, and show that it also gives local geometric convergence for every scalar relaxation sequence. The implicit solve reduces to a dense system whose order is the number of components. In the forward Burgers comparison, four components reduce the mean tail objective by 64.7% and the final solution error by 50.2% relative to one component at the same learning rate and momentum settings.

We study switching regret in adversarial multi-armed bandits, where the learner competes with an arm sequence that changes at most $S$ times. When $S$ is known, an optimal expected regret of $\widetilde{\mathcal{O}}(\sqrt{(S+1)KT})$ is obtainable [Auer et al., 2002]. However, when $S$ is unknown, Marinov and Zimmert [2021] show that this guarantee is impossible under an adaptive adversary. In this paper, we show that a single algorithm achieves $\widetilde{\mathcal{O}}(\sqrt{(S+1)KT})$ expected regret for every $S$ against an oblivious adversary, resolving an open problem of Auer et al. [2019b]. Our algorithm combines a fixed-share learner initialized with a small learning rate and dyadic-interval subroutines that search for local improvements using randomized learning rates and implicit exploration. Importantly, a non-uniform prior favors following the main learner, keeping the cost of maintaining many subroutines small. When the subroutines accumulate sufficient improvement over the main learner, its learning rate doubles, allowing adaptation to the unknown number of comparator switches $S$.

As compute scales, models evolve, and training algorithms advance, our ability to explain the increasingly powerful AI systems they enable is eroding. To help safeguard interpretability, we introduce a specialized training algorithm (MACCHIATO) that jointly constructs (i) an explicitly structured $\operatorname{ReLU}$-MLP from partial truth-table observations and (ii) an explicit Boolean circuit over signed literals with $\{\operatorname{AND},\operatorname{OR},\operatorname{XOR}\}$ gates certifying what its subnetworks compute and how they compose. Intuitively, we iteratively project the residuals of a Boolean function onto low-dimensional $\{\operatorname{AND},\operatorname{OR},\operatorname{XOR}\}$-circuit classes and exactly compile the resulting circuit into a $\operatorname{ReLU}$-MLP; we combine $\operatorname{ReLU}$-MLP circuit compilation, ESPRESSO logic minimization, and influence-based variable selection. Roughly speaking, our interpretability certificate is complemented by a statistical guarantee: under the theorem's influence-recovery conditions, if each of the $m$ stage-wise residuals depends on at most $\log_2(B)$ bits, a sample-splitting variant of our algorithm trained on $T$ observations returns a six-layer $\operatorname{ReLU}$-MLP (counting the input layer) of width $\mathcal{O}(mB)$ with truth-table error $\mathcal{O}\bigl(\sqrt{m(B+\log(m/\delta))/T}\bigr)$. On synthetic random-junta tasks, our networks outperform depth- and hidden-width-matched Adam-trained MLPs in several data-sparse or projection-aligned regimes, while the trained ReLU-MLPs are stronger in others. Moreover, in our explicit PyEDA truth-table implementation, the iterative procedure completes in regimes where flat ambient-dimensional ESPRESSO exceeds the three-hour computational budget.

Multi-objective Bayesian optimisation (MOBO) is a sample-efficient approach for optimising expensive black-box functions with multiple objectives. In MOBO, the goal is to adequately approximate the Pareto front; that is, to obtain a high-quality solution set with 1) good convergence (closeness to the Pareto front) and 2) good diversity (spread across the Pareto front). Existing MOBO methods typically aim to accomplish these two tasks simultaneously, i.e., driving the search towards the Pareto front while maintaining a diverse set of nondominated solutions, such that the solutions, ideally, can gradually approach the entire front. When sufficient search budgets are available, this approach is effective. However, considering both convergence and diversity throughout the search is not easy and requires careful design. Under very tight budgets, there may not be enough solutions generated to be able to simultaneously approach the entire Pareto front. To address this issue, this paper proposes a converge-then-diversify (CTD) approach that decouples convergence and diversity into two stages. In the first stage, CTD focuses on convergence, aiming to quickly drive the search toward a single point on the Pareto front. In the second stage, CTD focuses on diversity, aiming to spread solutions across the front. We present two simple instantiations of CTD by using widely adopted acquisition functions in the area. Experimental results show that, across all 446 pairwise comparisons, CTD statistically outperforms state-of-the-art methods in 72.9% of the cases, performs equivalently in 21.1%, and is statistically worse in only 6.1%, with the advantage being particularly evident in settings with very tight evaluation budgets or in high-dimensional problems.

Probabilistic forecasting is central to decision-making under uncertainty, yet its methodological landscape has become increasingly fragmented across temporal and spatiotemporal forecasting, statistical modeling, machine learning, and deep generative modeling. This survey develops a unified perspective by organizing probabilistic forecasting methods according to where and how uncertainty is introduced into the forecasting pipeline. Our taxonomy connects model-agnostic approaches including ensembles and distribution-free calibration, with model-intrinsic approaches spanning Bayesian modeling, parametric predictive distributions, distributional regression, and modern generative models, and further examines the emerging role of time series foundation models. Beyond methodological synthesis, we identify the assumptions, computational demands, and forms of uncertainty represented by different paradigms, and translate these distinctions into data-driven and domain-specific guidance for method selection. We complement the survey with a cross-paradigm empirical study on univariate, multivariate, and spatiotemporal forecasting tasks. The results reveal that no single uncertainty-quantification paradigm dominates across settings. Calibration, sharpness, predictive accuracy, and computational efficiency can lead to substantially different model preferences, while expressive generative models and zero-shot foundation models exhibit markedly different accuracy-efficiency trade-offs. Lastly, we identify unresolved challenges surrounding uncertainty in evolving dependency structures, physics-informed predictive distributions, forecasting extreme events, handling count-valued, directional, and continuous-time series, and the development of unified software resources. Our survey provides both a conceptual framework and a practical roadmap for probabilistic forecasting research.

Stochastic gradient descent (SGD) admits diffusion approximations that replace the complicated randomness of stochastic gradients by Gaussian noise, providing a powerful tool for understanding its dynamics and long-time behavior. We investigate whether an analogous approximation principle holds for optimization over probability measures, where the objective is a functional defined on the Wasserstein space P2. The nonlinear geometry and infinite-dimensional nature of P2 prevent a direct extension of the classical Euclidean theory. Using Lions differentiability, we lift the problem to a linear Hilbert space, where higher-order differential calculus becomes available. We then construct a Gaussian random-field approximation whose velocity field matches the mean and covariance of the original stochastic gradient. By exploiting this moment matching through higher-order Taylor expansions, we show that the Gaussian approximation captures the SGD dynamics with second-order weak accuracy. Our result provides a rigorous foundation for replacing sample-driven randomness by analytically tractable Gaussian fluctuations in stochastic optimization over probability measures.