A central mystery in deep learning is how neural networks, despite being highly non-convex and heavily overparameterized, are able to achieve near-zero training error while still generalizing well to unseen data. This paradox has sparked a surge of research aimed at understanding the convergence and generalization behavior of neural networks [1, 2, 6, 7, 15, 38, 41, 49]. The Neural Tangent Kernel (NTK), introduced by [20], has become one of a foundational tool for understanding the behavior of training dynamics for neural networks, especially those trained using gradient-based methods such as gradient descent (GD) and stochastic gradient descent (SGD). The core idea here is to linearize the neural network around its random initialization, which enables the evolution of the network during training to be closely approximated by a kernel method associated with the corresponding NTK. This framework establishes a powerful connection between the evolution of a neural network during training process and the behavior of kernel methods in a reproducing kernel Hilbert space (RKHS) induced by the NTK, allowing insights from the kernel methods to inform our understanding of neural networks. Following this perspective, the influential work [34] showed that for regression problems, shallow neural networks trained by SGD can achieve generalization performance on par with their kernel counterparts.

The transformer's emergent ability to perform in-context learning (ICL) has sparked a wide range of studies designed to understand its underlying mechanisms. Existing works often study how training task diversity, defined either as the number of ICL training task vectors or as the number of function classes from which the task vectors are drawn, shapes both the learning dynamics and generalization capabilities of ICL. While both definitions have uncovered many interesting phenomena, many observations under the latter definition remain theoretically unexplained. This paper presents a minimal analytical model under which these phenomena provably emerge from the properties of the training data. By modeling the training task vectors as a mixture of low-rank Gaussians, we show how training task diversity, defined by the number of non-overlapping columns between subspaces that parameterize the covariance matrices, improves both the generalization and optimization trajectory of ICL with linear attention. In particular, we show that our model can explain (i) why training with task diversity shortens the ICL plateau and (ii) why ICL appears to achieve out-of-distribution generalization. We conclude by empirically demonstrating how our results extend to nonlinear transformers and nonlinear function classes. Overall, our work presents a tractable framework to unify existing observations.

While algorithmic stability is a central tool for understanding generalization of learning algorithms, existing high-probability guarantees typically rely on uniform boundedness or sub-Gaussian/sub-Weibull tail assumptions, which can be overly restrictive for modern settings with heavy-tailed or unbounded losses. We develop a stability-based framework that requires only a finite $L_p$ moment condition. Our first contribution is sharp concentration inequalities for functions of independent random variables under $L_p$ constraints, extending McDiarmid's bounded-differences techniques beyond the classical regime. Leveraging these results, we derive sharp high-probability generalization bounds across a range of learning paradigms, including empirical risk minimization, transductive regression, and meta-learning. These guarantees show that $L_p$ stability suffices for robust generalization even when boundedness fails, substantially weakening the standard assumptions in the stability literature.

We study a finite-sample change-detection problem for one-dimensional noisy dynamical systems using partition-based empirical approximations of stationary behaviour. Given observations from an interval-valued process, we partition the state space into finitely many intervals and estimate a transition matrix from observed transitions between partition elements. After a small Doeblin-type regularisation, the resulting matrix has a unique stationary distribution. This stationary distribution is used as a finite-dimensional approximation of the invariant density, or stationary law, of the observed regime. Using an initial reference segment, we compute a baseline empirical stationary distribution bπ0,ρ. For each subsequent sliding window, we compute a window-based empirical stationary distribution bπt,ρ and define the score St = bπt,ρ bπ0,ρ 1. Large values of St indicate that the stationary behaviour of the observed regime has changed relative to the baseline. The statistic is therefore a detector of changes in stationary behaviour. It is not, by itself, a detector of all possible changes in transition dynamics that preserve the invariant density.

Robust Subspace Recovery (RSR) aims to identify an underlying d-dimensional subspace from a dataset heavily corrupted by outliers. Complexity-theoretic results establish a threshold for the problem's computational hardness based on the dimensionscaled signal-to-noise ratio (DS-SNR): the problem is SSE-hard when the DS-SNR is strictly less than 1, and solvable via practical algorithms when it is greater than 1 under general position assumptions. However, the exact behavior of practical algorithms at the critical boundary DS-SNR = 1 has remained unknown. Specifically, we prove that TME converges exactly to the true subspace for DS-SNR 1 under a new stability condition, which is less restrictive than the general position assumptions used in prior literature. I. Introduction Robust Subspace Recovery (RSR) is a fundamental problem in robust statistics, machine learning, and computer vision. The primary goal of RSR is to identify an underlying low-dimensional linear subspace from a dataset that is heavily corrupted by outliers. The standard formulation of the noiseless RSR problem assumes a dataset X = {xi}Ni=1 RD consisting of n1 inliers lying exactly on a d-dimensional linear subspace L RD, and n0 outliers lying strictly off L . We refer to such a dataset as a noiseless inlier-outlier dataset, where the total number of points is N = n0 +n1. The central algorithmic question in noiseless RSR is under what conditions one can exactly and efficiently recover the underlying d-subspace L . A natural metric for characterizing the difficulty of this problem is the ratio of inliers to outliers, n1/n0, which can be viewed as a signal-to-noise ratio (SNR) [8], [11], [12]. This leads to the dimension-scaled SNR (DS-SNR), denoted by δS: δS:= n1/d n0/(D d) . Hardt and Moitra [5] established a fundamental lower bound, showing that when δS < 1, the noiseless RSR problem is Small Set Expansion (SSE)-hard, a property conjectured to be equivalent to NP-hardness [15]. In the special case of hyperplanes (d = D 1), they showed NP-hardness by invoking a result from [7]. The noiseless RSR problem is SSE-hard if δS < 1.

In the sciences, regression tasks often require predicting high-dimensional outputs from few training examples. Multi-output Gaussian processes excel in low-data regimes but typically struggle with high-dimensional outputs. Compress-then-predict pipelines such as PCA-GP (principal component analysis plus Gaussian process regression) handle high dimensionality, but rely on bases optimized for reconstruction rather than prediction. To address this gap, we propose a model that represents each output as a linear-Gaussian decoding of a low-dimensional latent state drawn from a Gaussian process prior. By analytically marginalizing the decoder weights, we couple compression and prediction in a single objective that scales to high-dimensional outputs. We refer to this model as Gaussian process latent factor regression (GPLFR). We demonstrate GPLFR by building the first spatially resolved emulator of global climate models for rocky exoplanets.

Predicting outputs that are located in non-Euclidean spaces, such as probability distributions, networks, and symmetric positive-definite matrices, is becoming increasingly important in modern data analysis, particularly when inputs are high-dimensional. We propose DeSI (Deep Single-Index Fréchet Regression), a semiparametric framework for regression with metric space-valued outputs and multivariate inputs that assumes a single-index structure for the conditional Fréchet mean. DeSI estimates an interpretable index direction, which quantifies the relative importance of inputs, using a deep neural network, and performs Fréchet regression along the resulting one-dimensional index in the target metric space. This structure mitigates the curse of dimensionality while retaining interpretability, which stands in contrast to standard deep neural networks. We establish theoretical guarantees for DeSI, including uniform approximation and convergence rates, and demonstrate its strong predictive performance through simulations on distributions, networks, and symmetric positive-definite matrices, as well as an application to compositional mood data from New Jersey.

We study the minimax estimation error for distributed covariance matrix estimation in the vertical-split (feature-split) setting, where two agents each observe different coordinates of $m$ i.i.d. sub-Gaussian samples and communicate a limited number of bits to a central server. While Rahmani et al. [2025] established nearly tight bounds for dense (unstructured) cross-covariance matrices, we investigate whether imposing elementwise $s$-sparsity on the cross-covariance $C_{21}$ can reduce the required communication and sample complexity. In contrast to the horizontal-split setting, where Braverman et al. [2016] showed that sparsity does not reduce communication cost for mean estimation, we prove that sparsity does help for cross-covariance estimation in the vertical split. Specifically, we establish minimax lower bounds showing that the communication budget per agent scales as $B_k = Ω(σ^4 d_k\, s' \log(d_1 d_2/s')/\varepsilon^2)$ and the sample complexity for cross-covariance estimation as $m = Ω(σ^4\, s' \log(d_1 d_2/s')/\varepsilon^2)$, where $s' = s \wedge d_{\min}$. For the $1$-sparse case, this yields an exponential improvement from $d_1 d_2$ to $\log(d_1 d_2)$ compared to the dense rate. Our lower bounds are established via Fano's method with an explicit sparse packing using a Varshamov--Gilbert-type argument for signed partial permutation matrices combined with the Conditional Strong Data Processing Inequality of Rahmani et al. [2025]. We show the bounds are tight with a matching achievable scheme, based on covering-net quantization and entry-wise hard thresholding, that attains the $s$-sparse lower bound up to polylogarithmic factors.

Recent progress has been made in understanding the statistical generalization performance of gradient descent methods for overparameterized neural networks within the neural tangent kernel (NTK) regime. However, most of the existing work on regression problems is limited to shallow network architectures, leaving a notable gap in the theory of deep neural networks. This paper addresses this gap by presenting a comprehensive generalization analysis for deep ReLU networks trained using gradient descent (GD) and stochastic gradient descent (SGD). Specifically, we establish the first known minimax-optimal rates of excess population risk for both GD and SGD with deep ReLU networks, under the assumption that the network width scales polynomially with respect to the network depth and training sample size. Our results demonstrate that with sufficient width, gradient descent methods for deep ReLU networks can achieve optimal generalization rates on par with kernel methods.

Decision-making in complex systems often requires understanding counterfactuals of general, potentially highdimensional, interventions with limited data. Collecting sufficient data for every counterfactual in complex systems may be near impossible due to cost or ethical reasons. With the recent growth in expressivity and power in generative modeling, generative models that can synthesize counterfactual outcomes under generalized interventions stand as a viable solution for supporting robust decision-making in real-world systems. In an ideal world, we may simply train a generative model with the data we have, and sample from the generator under the intervention of interest. Counterfactual generative modeling may fail with such an approach due to confounding bias. Correlations observed in the sampled data may be mistaken for true causal effects, yielding incorrect downstream decisions. For example, generating medical images under changes in intervention dose can help track disease progression and identify optimal dosing strategies. However, if the training data primarily consisted of those who were responsive to intervention (e.g., younger populations), then the generator would identify the ranges in the data as effective even if this does not hold for different populations (e.g.